\section{Introduction}

A fundamental consequence of protostellar collapse is the formation of a viscous accretion disk that
transfers gas onto the stellar photosphere. Such disks are the progenitors of planetary systems and
their subsequent evolution and ultimate dispersal have extraordinary implications for the formation
of planets, orbital dynamics of planetary mass bodies, and the processing of dust grains and volatiles
within the disk. The timescale of disk dissipation within the terrestrial region has been reasonably
well-established by ground-based (Haisch et al. 2001; Mamajek et al. 2004) and {\it Spitzer Space
Telescope} (Uchida et al. 2004; Silverstone et al. 2006) observations to be $\le$10 Myr. By inference,
planetary systems must be in an advanced evolutionary stage within this period, before significant
quantities of gas and dust are depleted from the disk. Disk evolution is expected to proceed from the
interior outward (Dullemond \& Dominik 2005), however, sub-millimeter and mid-infrared observations by
Cieza et al. (2008) suggest that inner disks only begin to dissipate after the outer disk has been
significantly depleted of mass.

Pre-main sequence stars exhibiting spectral energy distributions (SED) that are indicative of an
optically thin disk interior surrounded by an optically thick outer disk are referred to as transition
disk objects (Strom et al. 1989; Najita et al. 2007; Muzerolle et al. 2010). Such disks are suggestive
of having experienced significant evolution from the continuous disk structures that are associated
with classical T Tauri stars (CTTS). The
duration of this transitional phase has been inferred from population statistics of star forming regions
to be only of order $\sim$10$^{5}$ yr (Hartmann 2009; Luhman et al. 2010). Transition-like SEDs,
however, can arise from various pathways and may not be representative of a direct evolutionary sequence
from primordial disk to debris disk (Najita et al. 2007). Dust grain growth and mid-plane settling
(Dullemond \& Dominik 2005), giant planet formation and dynamical clearing (Calvet et al. 2002),
disk photoevaporation (Alexander et al. 2006), and the presence of a stellar companion (Ireland \& Kraus 2009)
have been suggested as disk clearing mechanisms capable of producing transition-like SEDs.
An alternate evolutionary path from primordial disk to debris disk has been proposed by Currie et al. (2009),
who suggest that disks exhibiting reduced levels of near- and mid-infrared excess emission are
indicative of reduced masses of small dust grains at all disk radii. Such homologously depleted disks
counter the canonical inside-out disk evolutionary scenario.

The Upper Scorpius OB association is a critically important region for studies of disk evolution. At
$\sim$145 pc distant, it is among the nearest OB associations to the Sun (Blaauw 1991; de Zeeuw et al. 1999)
and has a well-established age of $\sim$5 Myr (Preibisch \& Zinnecker 1999; Preibisch et al. 2002), when
most ($\sim$80\%) optically thick, protoplanetary disks have dissipated (Haisch et al. 2001; Hernandez
et al. 2007). Carpenter et al. (2006) conducted a {\it Spitzer} 4.5--16.0 $\mu$m photometric survey
of 218 confirmed association members. These sources were compiled from
{\it Hipparcos} astrometry (de Zeeuw et al. 1999), color-magnitude diagrams and \ion{Li}{1} $\lambda$6708
follow-up observations (Preibisch \& Zinnecker 1999; Preibisch et al. 2002), and X-ray detected late-type
stars with \ion{Li}{1} $\lambda$6708 follow-up (Walter et al. 1994; Kunkel 1999; K\"{o}hler et al. 2000).
Given that the membership selection criteria were based upon stellar properties unrelated to circumstellar
material, the sample is considered to be unbiassed toward the presence or absence of disks. Carpenter et al. (2006)
found that 24 of 127 (19\%) K- and M-type stars in their sample exhibit infrared excesses similar to
Class II sources in the Taurus-Auriga star forming region. For comparison, Luhman et al. (2010) carried out a comparable {\it Spitzer}
IRAC and MIPS survey of 348 Taurus-Auriga members, finding that the disk fraction steadily declines
from $\sim$75\% for solar mass stars to 45\% for low-mass stars and brown dwarfs (0.3--0.01 M$_{\odot}$).

To better characterize disk emission among the Upper Scorpius infrared excess stars identified by
Carpenter et al. (2006), Dahm \& Carpenter (2009) examined {\it Spitzer} mid-infrared spectra for a
substantial fraction (26 of 35) of these sources: 8 early-type (B+A) stars and 18 late-type (K+M) stars.
In general, excess emission among the late-type Upper Scorpius population becomes apparent between 2.2 and
4.5 $\mu$m, a region sampled only by Two-Micron All Sky Survey (2MASS) $K_{S}-$band photometry and the {\it Spitzer}
Infrared Array Camera (IRAC) [4.5] fluxes. Compared to Class II sources in Taurus-Auriga, Dahm \& Carpenter (2009)
found the disk population in Upper Scorpius to exhibit reduced levels of near- and mid-infrared excess emission and
an order of magnitude lower mass accretion rates. The apparent abundance of depleted inner disk systems in
Upper Scorpius relative to the number of such objects in Taurus-Auriga led Dahm \& Carpenter (2009) to suggest
that such disks represent a common evolutionary pathway. Near-infrared spectra spanning the 2.2--4.5 $\mu$m
region, however, were needed to complement the {\it Spitzer} IRAC and Infrared Spectrograph (IRS) observations
and to isolate the onset of continuum excess emission from the terrestrial disk regions.

To further examine emission arising from the depleted disk interiors of the Upper Scorpius sample, 0.8--5.2 $\mu$m
moderate-resolution spectra for 12 of 18 late-type disk-bearing stars from Dahm \& Carpenter (2009) were obtained
using SpeX on the NASA Infrared Telescope Facility (IRTF) on Mauna Kea. The SpeX spectroscopy permits the unambiguous
detection of continuum excess emission arising from hot dust in the disk interiors of the Upper Scorpius sample.
Merging these observations with {\it Spitzer} IRS spectra and Multi-band Imaging Photometer for {\it Spitzer} (MIPS)
24 and 70 $\mu$m photometry from Carpenter et al. (2009), the resulting 0.8--70 $\mu$m SEDs are compared with the
two-dimensional, radiative transfer accretion disk models of Robitaille et al. (2006). These models are used to
constrain inner disk radii for the observed Upper Scorpius excess sources and to examine their SEDs for evidence
of inner disk evolution.

In Section 2 the observed Upper Scorpius sample is described and placed into context with the greater late-type
stellar population of the association. Details of the SpeX observations and analysis are also provided. The
near-infrared continuum excess spectra and blackbody fits to these spectra are discussed in Section 3. Next
(Section 4) accretion luminosities and mass accretion rates are derived for suspected accretors using Pa$\beta$
and Br$\gamma$ emission line luminosities. The merged 0.8--70 $\mu$m SEDs are then compared in Section 5
with the accretion disk models of Robitaille et al. (2006). Constraints for inner disk radii are derived
and the effects of mid-plane settling are discussed using the best-fitting accretion disk models. Finally the
results of this paper are summarized in Section 6.

\section{Observations and Analysis}

\subsection{The Upper Scorpius Membership Sample}

The 12 Upper Scorpius infrared excess sources observed with SpeX were selected from the mid-infrared photometric
and spectroscopic surveys of Carpenter et al. (2006) and Dahm \& Carpenter (2009), respectively. These sources
represent half of the late-type, disk-bearing stars identified by Carpenter et al. (2006) and include all K-type
and all M0 through M2-type excess sources. To further place this sample into context with the low-mass stellar
population of Upper Scorpius, there are $\sim$250 known pre-main sequence stars in the mass range from $\sim$0.1
to $\sim$2.0 M$_{\odot}$ (Preibisch \& Mamajek 2008; Preibisch et al. 2002). If the infrared excess fraction
for late-type stars (19\%) derived by Carpenter et al. (2006) is assumed, $\sim$48 of these should host primordial
disks. The SpeX sample therefore represents $\sim$25\% of all disk-bearing stars expected among the known
late-type association members. Integrating the best-fitting mass function for Upper Scorpius, Preibisch \& Mamajek (2008)
estimate that the total low-mass (0.1--2.0 M$_{\odot}$) stellar population of the association exceeds 2400 stars,
corresponding to a total disk-bearing population of $\sim$450 late-type stars. While the SpeX sample considered
here does represent a substantial fraction of the late-type, disk-bearing stars identified by Carpenter et al. (2006),
it cannot be considered a statistically significant representation of the total low-mass stellar population of the
association.

Given that unresolved companions could provide an explanation for transition-like SEDs (Ireland \& Kraus 2008),
some knowledge of the binary frequency of the Upper Scorpius late-type stellar population is warranted.
K\"{o}hler et al. (2000) used speckle interferometry and direct imaging for 118 X-ray selected T Tauri stars
in the greater Scorpius-Centaurus OB association to identify companions with separations ranging from 0.13\arcsec\
to 6.0\arcsec. They found a multiplicity fraction of $\sim$32.6\%, which exceeds the binary fraction of
main sequence stars by a factor $\sim$1.6, but is slightly lower than that observed in Taurus-Auriga
(K\"{o}hler et al. 2000). Kraus et al. (2008) found a similar frequency of binary companions among 82
late-type (G0--M4) Upper Scorpius members,  $\sim$35$\pm$5\%. Two of the disk candidates observed with SpeX
are established binaries, [PZ99]J161411.1-230536 and ScoPMS 31, resolved by high angular resolution imaging.
With projected separations of $\sim$32 and 84 AU, respectively (Metchev \& Hillenbrand 2009; K\"{o}hler et al. 2000),
it is unlikely that these stellar companions have significantly influenced the inner disk evolutionary
timescales of their primaries. Kraus et al. (2008) observed J160823.2-193001 and J160900.7-190852
using an aperture mask interferometry technique, placing firm upper limits that preclude the possibility
of an undetected stellar companion within $\sim$5 AU of these stars. The multiplicity fraction among the
remaining Upper Scorpius disk-bearing sample is relatively unexplored, but is not expected to deviate
significantly from that derived by K\"{o}hler et al. (2000) and Kraus et al. (2008) for the pre-main sequence
population as a whole.

Shown in Figure 1 is the $J-H$, [8.0]$-$[4.5] color-color diagram for the 218 Upper Scorpius members included
in the Carpenter et al. (2006) {\it Spitzer} IRAC and IRS survey. The 12 sources observed with SpeX span the
full range of the 8.0 $\mu$m excess distribution, from moderate excess (e.g. J160643.8-190805, ScoPMS 31)
to significant (e.g. J160900.7-190852). Their masses, derived using the pre-main sequence models of Siess et al. (2000)
and assuming a distance of 145 pc, range from $\sim$0.25 to 2.0 M$_{\odot}$. General properties of the
stellar sample including spectral type, extinction ($A_{V}$), mass, radius, luminosity, and effective
temperature are presented in Table 1.

\subsection{SpeX Observations}

The SpeX (Rayner et al. 2003) observations were made on the nights of 2009 May 20--23 under variable cirrus
and seeing conditions (0.6--0.9''). The Upper Scorpius members were observed in the short cross-dispersed SXD
(0.84--2.4 $\mu$m) and long cross-dispersed LXD (2.2--5.2 $\mu$m) modes for complete near-infrared wavelength
coverage. Bright A0 V stars were observed at similar airmasses ($\Delta$$\chi$ $<$ 0.1) and within short
periods of time of the program stars for telluric correction. Arc-lamps and internal flat-field exposures
were obtained for each set of observations to account for instrument flexure. All observations were made
with the 0.5\arcsec\ slit yielding a nominal spectral resolution of $\lambda$/$\delta$$\lambda$$\sim$1500.

The spectra were reduced using {\it Spextool}, an IDL-based reduction package that provides for sky-subtraction,
flat-fielding, wavelength calibration, and optimal extraction (Cushing et al. 2004). For telluric corrections
{\it xtellcor}, an extension package of {\it Spextool}, was used to create a kernel for the telluric spectrum
using a model spectrum of $\alpha$ Lyra. {\it xtellcor} interpolates over broad hydrogen absorption lines in
the spectra of the A0 V stars using a technique developed by Vacca et al. (2003). The orders of the telluric
corrected spectra were then combined and the SXD and LXD spectra merged using routines available within
{\it Spextool}. The resulting continuum levels across all orders and between the SXD and LXD modes were found
to be remarkably consistent. The near infrared spectra were then merged with the {\it Spitzer} IRS spectra.
Details of the {\it Spitzer} IRS observations and their subsequent reduction and analysis can be found in
Dahm \& Carpenter (2009).

\section{Near-Infrared Continuum Excess Spectra}

The extinction-corrected, near-infrared (0.84--5.2 $\mu$m) spectra of the 12 Upper Scorpius sources are shown
in Figure 2. Regions of strong telluric absorption (atmospheric transmission $<$20\%) have been excised from
the figure. Also shown are the spectra of solar metallicity, main sequence stars of identical or closely matched
($\pm$1 sub-class) spectral type obtained from the IRTF spectral library (Rayner et al. 2009). These standards
have reliable spectral types and are directly linked to the MK classification system. The standard or template
spectra are scaled to the flux levels of the Upper Scorpius sources at 1.65 $\mu$m, near the peak of the
stellar SED and where extinction effects are minimized (D'Alessio et al. 1999; Furlan et al. 2006).

To characterize excess emission attributable to the inner disk rim, the scaled photospheric template spectra
are subtracted from the dereddened spectra of the Upper Scorpius sample. The resulting continuum excess
spectra are shown in Figure 3. Significant noise is present due to imperfect telluric correction, particularly
in the thermal region where the water column varies both temporally and with airmass. In general, the slopes
of the continuum excess spectra are increasing from $K-$band to $\sim$3.5 $\mu$m before turning over toward
redder wavelengths. The shapes of many of the continuum excess spectra are suggestive of having been produced by
single-temperature blackbodies. By fitting these excess spectra with Planck functions, the characteristic
temperatures of the blackbodies are found to range from the sublimation temperature for silicate dust, i.e.
$\sim$1400 K, to $\le$500 K. The slope of excess emission, particularly between $\sim$2.75 and 4.2 $\mu$m, was
found to be of critical importance when fitting the blackbody profiles to the continuum excess spectra.

To estimate the uncertainty associated with the derived blackbody temperatures, photospheric templates within
a range of $\pm$1 sub-class of the assigned spectral type were subtracted from the Upper Scorpius spectra.
The best-fitting blackbody curves for this range of photosphere-subtracted excess spectra suggest typical
uncertainties of $\sim$$\pm$200 K for most sources. The adopted blackbody curve as well as curves representing
the limits of uncertainty in the characteristic temperature are shown in Figure 3. The characteristic temperatures
of the best-fitting blackbodies are provided in Table 1 and are assumed to be representative of the dust
temperatures in the disk interiors.

The disk fractions of Upper Scorpius and Taurus-Auriga differ significantly for the range in stellar masses
considered here: 19\% for the former (Carpenter et al. 2006) and up to 75\% for the latter (Luhman et al. 2010).
Muzerolle et al. (2003) used SpeX in LXD-mode to observe 9 CTTSs in Taurus-Auriga, 8 of which span a range
of spectral types (G5--M1), luminosities (12.8--0.5 L$_{\odot}$), and stellar radii (3.6--1.8 R$_{\odot}$) that are
comparable with those of the Upper Scorpius disk sample. Restricting further comparison of these samples to
just those sources having similar physical properties (i.e. spectral type, luminosity), the spectral profiles
of continuum excess emission of the Upper Scorpius disks are found to differ from those of Class II sources
in Taurus-Auriga. For 6 of 9 members in the reduced Upper Scorpius sample, the characteristic temperatures
of excess emission fall below the dust sublimation temperature, a condition satisfied by
only one member (the binary DQ Tau) of the Taurus-Auriga sample. Furthermore, the excess profiles of
J160643.8-190805 and J160357.9-194210, which are not represented in Figure 3, exhibit no distinctive shape,
implying minimal continuum excess emission at near-infrared wavelengths. Shallow rises may be present near
$\sim$3.75 $\mu$m for these 2 sources, suggestive of warm ($\le$500 K) dust emission, but higher SN spectra
are needed for confirmation. The excess spectra of J161420.2-190648, [PZ99]J160357.6-203105, J160900.7-190852,
and J155829.8-231007 all have characteristic temperatures at or near the dust sublimation temperature,
implying that emission arises from dust in close proximity to the host star. Of these sources, all are
suspected accretors (Section 4).

\section{Accretion Luminosity and Mass Accretion Rates}

Dahm \& Carpenter (2009) identified 5 Upper Scorpius members in the present sample as probable accretors using
the H$\alpha$ velocity width criteria of White \& Basri (2003): [PZ99]J160357.6-203105 (K5), J160900.7-190852 (K7),
J161420.2-190648 (M0), ScoPMS 31 (M0.5), and J155829.8-231007 (M3). The near-infrared spectra of these sources
exhibit \ion{He}{1} $\lambda$10830, Pa$\gamma$, and Pa$\beta$ emission and for most, weak Br$\gamma$ emission.
Shown in Figure 4 are the normalized spectra of these sources centered upon Pa$\gamma$, Pa$\beta$, and
Br$\gamma$. Dahm \& Carpenter (2009) used veiling at $\lambda$6500 \AA\ as well as \ion{Ca}{2} $\lambda$8542
emission line luminosity to estimate mass accretion rates ($\dot{M}$) for these stars. In the near-infrared,
Muzerolle et al. (1998) found Pa$\beta$ and Br$\gamma$ emission line luminosity ($L_{Pa\beta}$ and $L_{Br\gamma}$)
to be well-correlated with accretion luminosity, $L_{acc}$. The resulting least-squares fits from Muzerolle et al. (1998)
are given by:

log$(\frac{L_{acc}}{L_{\odot}})=(1.14\pm0.16)$log$(\frac{L_{Pa\beta}}{L_{\odot}})+3.15\pm0.58$   (1)

\noindent and

log$(\frac{L_{acc}}{L_{\odot}})=(1.26\pm0.19)$log$(\frac{L_{Br\gamma}}{L_{\odot}})+4.43\pm0.79$   (2).

For the 5 suspected accretors in the Upper Scorpius sample, template spectra of identical or similar spectral
type from the IRTF spectral library of Rayner et al. (2009) were used to subtract Pa$\beta$ and Br$\gamma$
photospheric absorption. Pa$\beta$ and Br$\gamma$ line luminosities were then determined using their measured
equivalent widths and the extinction-corrected $J-$ and $K_{S}-$band magnitudes obtained from the 2MASS point
source catalog. The resulting emission line luminosities were then transformed into $L_{acc}$ using the above
linear relationships of Muzerolle et al. (1998). The derived $L_{acc}$ values are directly proportional to
$\dot{M}$ such that:

$L_{acc}\sim \frac{GM_{*}\dot{M}}{R_{*}}(1-\frac{R_{*}}{R_{in}})$   (3)

\noindent where the stellar mass and radius estimates used are those predicted by the pre-main sequence models
of Siess et al. (2000). The factor of ($1 - R_{*}/R_{in}$) is assigned a value of 0.8, which assumes an inner
disk radius ($R_{in}$) of 5 R$_{*}$ (Gullbring et al. 1998). Given the possibility that inner disk radii for
the Upper Scorpius sources may exceed those of typical Class II sources in Taurus-Auriga, this value could
be underestimated by a factor of $\sim$1.25. The derived $L_{acc}$ and $\dot{M}$ values with their associated
uncertainties are presented in Table 2. These uncertainties arise from multiple sources: error in the measured
equivalent width (assumed to be $\sim$20\%), spectral type uncertainty when correcting for photospheric
absorption ($\pm$1 sub-class), the uncertainty in each coefficient of the Muzerolle et al. (1998) relationships,
and the uncertainty in R$_{in}$ when determining $\dot{M}$ values.

In general the $\dot{M}$ values derived using $L_{Pa\beta}$ and $L_{Br\gamma}$ agree reasonably well with each
other and with the $\dot{M}$ values from the veiling and \ion{Ca}{2} $\lambda$8542 analysis of Dahm \& Carpenter (2009).
In summary, 5 of 12 of the Upper Scorpius sample are accreting, providing unambiguous evidence for the presence of
gas within the terrestrial regions of these disk-bearing systems.

\section{Accretion Disk Model Analysis}

To further examine inner disk structure and to constrain inner disk radii for the Upper Scorpius
disk-bearing sample, the observed SEDs are compared with those predicted by the accretion disk models
of Robitaille et al. (2006). To facilitate the comparison, the flux-calibrated SpeX observations
presented here are merged with the {\it Spitzer} IRS low-resolution spectra and MIPS 24 and 70 $\mu$m fluxes taken
from Dahm \& Carpenter (2009) and Carpenter et al. (2009), respectively. The resulting 0.8--70 $\mu$m
SEDs are well sampled between 2.2 and 24 $\mu$m, a spectral region dominated by disk emission originating
from the terrestrial region.

The Robitaille et al. (2006) grid of pre-computed, two-dimensional radiative transfer models consists of
20,000 young stellar objects in varying stages of evolution and viewed from 10 inclination angles. A total
of 14 parameters are randomly sampled that specify stellar (e.g. $M_{*}$, $R_{*}$, $T_{eff}$) as well as
disk (e.g. $M_{disk}$, $\dot{M}$, $R_{min}$) properties. The models have been successfully applied to the
SEDs of 30 spatially resolved Class I and II sources in Taurus-Auriga by Robitaille et al. (2007), including
the transition disk objects GM Aur and DM Tau.

The SEDs of 9 Upper Scorpius members were compared to the accretion disk models. The sources excluded from
the model fitting analysis were J160545.4-202308 and J155829.8-231007, which were not observed with IRS,
and [PZ99]J160421.7-213028, which exhibited significant (factor of $\sim$4) mid-infrared variability
(Dahm \& Carpenter 2009). To compare the observed SEDs with the accretion disk models, fluxes were measured
in $\sim$30 narrow passbands defined between 0.8 and 37.0 $\mu$m. The MIPS 24 and 70 $\mu$m fluxes were
also incorporated, with the 70 $\mu$m fluxes being used to constrain disk emission originating from beyond
the terrestrial region.

Two assumptions were made in the model fitting analysis: that all Upper Scorpius members lie between 125
and 185 pc distant, and that foreground extinction is no more than $A_V \sim 3$ mag, both reasonable
assumptions based upon substantial Upper Scorpius literature (e.g. Preibisch \& Zinnecker 1999;
Preibisch et al. 2002; Carpenter et al. 2006). A range of effective temperatures was defined for each
Upper Scorpius source based upon its adopted spectral type and assuming an uncertainty of $\pm$1
spectral sub-class. Only models having effective temperatures within this specified range were considered
when fitting the observed SEDs. The model fitting program was written in IDL and uses the $\chi^{2}$ statistic
to evaluate the goodness of fit of each model for each assumed distance, extinction value, and inclination
angle. Table 3 lists the range of $T_{eff}$ values considered for each Upper Scorpius source, the number of
models included within the specified range, the best-fitting model number, its corresponding $A_{V}$, distance,
and inclination angle, and the minimum reduced $\chi^{2}$ value achieved by the fit.

Shown in the upper panels of Figures 5a--i are the observed 0.8--70.0 $\mu$m SEDs for the 9 Upper Scorpius
sources included in the model fitting analysis. Superimposed in red are the best-fitting model profiles
placed at their respective distances and extinctions. The disk contributions to the model SEDs are represented
by the dotted black curves, and the model stellar atmospheres, taken from the Kurucz (1979) atlas, are
plotted in blue. In general, reasonable fits were achieved for many of the Upper Scorpius disk-bearing
stars despite the modest number of models (typically $\le$2000) examined for each source. The near-infrared
spectrum of the strongly accreting source J161420.2-190648 is suggestive of significant extinction,
$A_V\sim$6 mag, (Figure 5f). The optimal model fits for this source yielded the largest $\chi^{2}$ values
of the sample given the a priori assumption that $A_V \le 3$ mag.

Degeneracy among model SEDs having significantly different stellar and disk properties is an inherent
limitation of the SED fitting process. It is possible, however, to place constraints upon a specific
parameter by examining the range of values returned by the best-fitting models. To identify this subset
of models, $\chi^{2} - \chi^{2}_{min} < N$ is used as a statistical measure of goodness of fit, where
$\chi^{2}_{min}$ is the minimum $\chi^{2}$ value returned by the fitting process and $N$ is set to the
1-$\sigma$ confidence interval for $\nu$ degrees of freedom. In Table 4 the disk properties
of the Upper Scorpius sample obtained from these best-fitting subsets are summarized. Tabulated for each source are
the minimum, best-fitting, and maximum values of disk mass (M$_{disk}$), $\dot{M}$, and inner disk radius
(R$_{in}$). The ranges of M$_{disk}$ and $\dot{M}$ are significant, up to several orders of magnitude
separating the minimum and maximum values for a given source. The observed dispersion in M$_{disk}$ arises
in part from the lack of far-infrared or millimeter wavelength observations that could constrain outer
disk emission where substantial quantities of disk mass could remain unaccounted for. Mass accretion rates
calculated by the models were found by Robitaille et al. (2007) to be overestimated in their sample of
Taurus sources. In the Upper Scorpius sample, whether established accretors or not, $\dot{M}$ is
observed to range at least 3 orders of magnitude, and in some cases up to 5. For most suspected accretors
included in the model fitting analysis, the $\dot{M}$ values derived in Section 4 using the Pa$\beta$ and
Br$\gamma$ emission line luminosities are consistent with those predicted by the best-fitting models.

As might be expected, inner disk radius appears to be better constrained by the subsets of best-fitting
models than either M$_{disk}$ or $\dot{M}$. Several Upper Scorpius sources, e.g. [PZ99]J161411.0-230536,
[PZ99]J160357.6-203105, and ScoPMS 31, exhibit ranges in $R_{in}$ from minimum to maximum of an order of
magnitude or less. Five sources have best-fitting inner disk radii that exceed their respective dust
sublimation radii: [PZ99]J161411.0-230536 ($R_{in}$=0.87 AU), J160900.7-190852 (2.18 AU), ScoPMS 31 (9.62 AU),
J161115.3-175721 (0.29 AU), and J160357.9-194210 (0.17 AU). Given the poor quality of the $\chi^{2}$ fit
for the strongly accreting source J161420.2-190648, the best-fitting $R_{in}$ value of 0.8 AU is regarded
with skepticism. The remaining Upper Scorpius sources ([PZ99]J160357.6-203105, J160643.8-190805, and
J160823.2-193001) have best-fitting inner disk radii that are consistent with their sublimation radii.
The confidence levels for these inner disk radii, however, are best examined in terms of a probability
density function determined from the $\chi^{2}$ values of the accretion disk model fits.

\subsection{Constraints for Inner Disk Radii}

Flux from the terrestrial disk regions is well sampled by the SpeX and IRS observations, which
are dominated by disk emission extending from the sublimation radius to more than 20 AU for a typical
$\sim$0.7 M$_{\odot}$ pre-main sequence star. The Upper Scorpius disks are presumably at an advanced
evolutionary stage relative to those found around Class II sources in Taurus-Auriga and exhibit SEDs
that are consistent with reduced levels of near and mid-infrared disk emission (Dahm \& Carpenter 2009).
Robitaille et al. (2007) find that all Taurus sources, except for the known transition disk objects
(e.g GM Aur, DM Tau) can be fit by models having disks and envelopes with inner disk radii equal to
the dust sublimation limit. Approximately one-third of the Robitaille et al. (2006) models have inner
disk radii set to the dust destruction radius. The remaining models have increasing inner disk radii
that span from the dust destruction radius to 100 AU. The inner disk gaps are treated by the models as
being completely evacuated of dust (Robitaille et al. 2006).

To constrain R$_{in}$ for the Upper Scorpius disk sample a probability ($P$) is calculated for a given
inner disk radius using the returned $\chi^{2}$ values from the model fits of the observed SEDs:

$P=\frac{1}{\sqrt{2 \pi}} e^{\frac{-\chi^{2}}{2}}$   (4)

\noindent This probability density function assumes a normal distribution for the returned $\chi^{2}$
values of the individual fits. In the center panels of Figures 5a--i, $P$ is plotted as a function of
R$_{in}$ for each Upper Scorpius source. Superimposed in the figures as cross-hatched histograms are
the distributions of inner disk radii for the entire sample of models considered for each source.
Only 3 Upper Scorpius disks have predicted inner radii that are larger than their respective dust
sublimation radii at the 1$\sigma$ confidence level or greater: [PZ99]J161411.0-230536 (2-$\sigma$),
J160900.7-190852 (1-$\sigma$), and ScoPMS31 (3-$\sigma$). The remaining sources have inner disk radii
that are either most consistent with the sublimation radius or exhibit a broad dispersion among predicted
$R_{in}$ values (e.g. J161115.3-175721). The SED of J160643.8-190805,
which exhibits minimal near-infrared continuum excess emission (Section 3), is readily fit by a large number
of models having a range of inner disk radii. This source exhibits an SED that is reminiscent of those
defined as anemic (Lada et al. 2006), weak (Dahm \& Hillenbrand 2007), or homologously depleted
disk-bearing systems (Currie et al. (2009).

\subsection{Disk Mid-plane Settling}

The current disk evolution scenario suggests that dust grains coagulate and settle toward the mid-plane
prior to the formation of large planetesimals. The Robitaille et al. (2006) models use 2 disk structure
parameters to mimic the effects of dust settling: a disk flaring parameter ($\beta$) and a disk scale
height factor ($z$). If both of these parameters are low for a given model, this could be indicative of
dust mid-plane settling (Robitaille et al. 2006). In the lower panels of Figures 5a--i $\beta$ is plotted
as a function of $z$ for the entire sample of models considered for each source. Superimposed in red are
the best-fitting models as determined by the $\chi^{2} - \chi^{2}_{min} < N$ measure of goodness of fit.
In general the spread in both parameters for the best-fitting models is significant, suggesting that
neither is well-constrained. Some argument can be made that the SEDs of [PZ99]J161411.0-230536,
J160643.8-190805, and ScoPMS31 are better fit by models having lower than average values of $\beta$.
The disk scale height factors for these sources, however, are found to vary significantly. Evidence for
dust settling effects among the Upper Scorpius primordial disk sample remains inconclusive at best.

\section{Discussion and Conclusions}

The clearly different near- and mid-infrared color excess distributions of the Upper Scorpius and
Taurus-Auriga Class II populations found by Dahm \& Carpenter (2009) suggest that the Upper Scorpius
stars have experienced some degree of inner disk evolution relative to their presumably younger
counterparts in Taurus-Auriga. This is circumstantially supported by the characteristic temperatures
of the continuum excess emission for most Upper Scorpius sources, which are substantially lower than
the sublimation temperature for silicate dust. This condition is satisfied by only one member of the
Muzerolle et al. (2003) Taurus-Auriga sample of Class II sources having similar spectral types and
luminosities. The accretion disk model fitting results, however, are unable to effectively constrain
inner disk radii for the majority of the Upper Scorpius disk-bearing sample. This could in part result
from the limited number of models examined for each source, which do not sample the full range of
available parameter space of disk structure. Another limitation of the Robitaille et al. (2006) models
is the assumption that disk gaps are completely devoid of dust, a reduction in complexity that impacts
the predicted near- and mid-infrared excess distributions, critical for constraining disk emission
originating from the inner disk rim.

In summary 3 of 9 late-type, disk-bearing stars in the Upper Scorpius sample exhibit SEDs that are most
consistent with having inner disk radii that lie beyond the sublimation radius for silicate dust:
[PZ99]J161411.0-230536, J160900.7-190852, and ScoPMS31. The best-fitting models for 2 additional
sources: J161115.3-175721 and J160357.9-194210, have inner disk radii that substantially exceed their
respective dust sublimation radii, but the probability distributions of $R_{in}$ values for these
sources exhibit significant dispersion. Proposed indicators for the effects of mid-plane settling
in the best-fitting Robitaille et al. (2007) models of the Upper Scorpius sample (i.e. concurrent
decreased values of the disk flaring parameter and disk scale height factor) are inconclusive, leaving
open to question the nature of the remaining Upper Scorpius primordial disks and their evolutionary state.

Transition disks are believed to represent an early stage of disk clearing and may be in the process
of rapid grain growth, planetesimal and planet formation. Adopting the disk classification scheme of
Luhman et al. (2010), many Upper Scorpius sources would be classified as pre-transitional objects,
i.e. disks that exhibit reduced emission at wavelengths $\le$10 $\mu$m, but that still retain
significant disk emission at longer wavelengths. Two of the Upper Scorpius pre-transitional candidates
appear to exhibit gapped disk structure: J160900.7-190852 and ScoPMS31. Both sources are accreting
(Dahm \& Carpenter 2009 and Section 4) and both exhibit SEDs that are most consistent with having
substantial inner disk radii, $\sim$1--10 AU. At least one Upper Scorpius source, J160643.8-190805,
exhibits an infrared SED that is suggestive of an homologously depleted disk system as defined by
Currie et al. (2009).

Before attributing the reduced levels of infrared excess emission in the Upper Scorpius sample to
disk evolutionary processes, binarity must be considered as an explanation for the cleared out inner
cavities (e.g. CoKu Tau 4, Ireland \& Kraus 2008). Two of the disk candidates included in this analysis,
[PZ99]J161411.1-230536 and ScoPMS 31, are established binaries resolved by high angular resolution
imaging. It must also be acknowledged that the SpeX sample represents only $\sim$25\% of all primordial
disk-bearing stars expected among the $\sim$250 known late-type members of the Upper Scorpius OB association.
Observations of more late-type, disk-bearing systems are critically needed to confirm the results
presented here.

Improved modeling of the Upper Scorpius primordial disk-bearing sample is clearly needed to provide
better constraints for the inner disk structure of these presumably evolved primordial disk systems.
Monte Carlo three-dimensional, radiative transfer codes are now available that could be applied to
the SEDs of these sources. High angular resolution imaging is also needed to identify close ($\le$10 AU)
binary companions that may account for the reduced near- and mid-infrared excess emission of these
sources relative to Class II sources in Taurus-Auriga. Infrared interferometric observations or precision
radial velocity monitoring are also needed to identify tighter pairs or spectroscopic binaries that
would be capable of dynamically clearing the disk interiors of these systems.

\acknowledgments
This work is based on observations made with the {\it Spitzer} Space Telescope, which is operated by the
Jet Propulsion Laboratory (JPL), California Institute of Technology, under NASA contract 1407. The
Digitized Sky Surveys, which were produced at the Space Telescope Science Institute under
U.S. Government grant NAG W-2166, were used as were the the SIMBAD database operated at CDS,
Strasbourg, France, and the Two Micron All Sky Survey (2MASS), a joint project of the University of
Massachusetts and the Infrared Processing and Analysis Center (IPAC)/California Institute of Technology,
funded by NASA and the National Science Foundation. The author gratefully acknowledges John Carpenter and
an anonymous referee for many helpful suggestions that significantly improved this manuscript. The
author also thanks Elise Furlan, Thayne Currie, and Joan Najita for insightful discussions and
Thomas Robitaille for graciously providing direct access to the accretion disk models. Finally
the author thanks John Rayner for obtaining a second-epoch SpeX SXD spectrum of [PZ99]J160421.7-213028.


\section{Introduction}

Newton's root-finding method is as old as analysis, but still not
well understood, even in the fundamental case of finding all roots of
a polynomial in a single variable. Its local convergence properties
are well known; near simple roots convergence is quadratic and thus
extremely rapid. However, the global dynamical properties are
insufficiently understood so that numerical analysis algorithms often
use different global methods, and resort to Newton's method for a
final local ``polishing'' of the roots.

This article is a contribution towards a better understanding of the
global properties of Newton's method, applied to polynomials in a
single complex variable. Even for polynomials over the reals, and
even if all the roots are real, it is often preferable to use complex
methods; see Figure~\ref{Fig:NewtonDynamics}.

\begin{figure}
\includegraphics[width=0.49\textwidth]{NewtonDeg7Plane.pdf}
\includegraphics[width=0.49\textwidth]{NewtonDeg4Real.pdf}
\caption{Dynamical planes of Newton maps of two complex polynomials.
Different colors illustrate basins of attraction of different roots;
shades of color illustrate different speeds of convergence. It is
clearly visible that all immediate basins are unbounded and have one
or several channels to $\infty$ of different widths. Left: a
polynomial of degree $7$. Right: a polynomial of degree $11$ with all
roots real. Some of the roots are very close to each other; however, away
from the disk containing all the roots, the basins and their channels
all have almost uniform width, so that finding the real roots using
complex methods is much easier.}
\label{Fig:NewtonDynamics}
\end{figure}

\begin{figure}
\includegraphics[width=0.51\textwidth,trim=0 130 0 60,clip]{NewtonDeg11ThinRoots.pdf}
\caption{The dynamical plane of the polynomial $p(z)=z(z^{10}-1)$: the ten roots of unity each have one ``thick'' channel, while the root $z=0$ has 10 channels (red) which are all rather ``thin''. The deterministic method from \cite{HSS} would search for the individual thin channels and thus requires more points, while our method searches for the union of all thin channels, which together are much bigger.}
\label{Fig:NewtonDynamicsThinChannels}
\end{figure}

Among the difficulties with Newton's method are the following:
\begin{itemize}
\item
if an orbit under iteration comes close to a critical point of the
polynomial, the Newton map sends the orbit far away near $\infty$, so
that control of the dynamics is lost, and in any case a large number
of iterations are required until the orbit comes back to where the
roots are;
\item
there are polynomials with open sets of starting points that do not
converge to any root (Smale~\cite{Smale} asked, in 1984, for a
classification of such polynomials; an answer has recently been given by
Mikulich in current work \cite{ZhenyaThesis});
\item
the boundary of the basins of convergence for the roots may have
positive planar Lebesgue measure (this follows from recent work by
Buff and Ch\'eritat on the existence of Julia sets with positive
measure \cite{BuffCheritat}, combined with Douady and Hubbard's
renormalization theory \cite{PolyLike});
\item
even if almost every point in $\C$ converges to some root under the
Newton iteration, our goal is to find \emph{all} roots of the
polynomial, and with bounded complexity. Finding some roots and
deflating is usually not an option, because deflation is in general
numerically unstable (unless the roots are found in a specific
order), and because deflation might not be compatible with the way
the polynomial may be specified, or evaluated efficiently (for
instance, if the polynomial itself is given by an efficient iteration
procedure).
\end{itemize}
See \cite{Johannes} for a recent survey of known results on Newton's method.


\section{Introduction}
\label{sec:introduction}
Top-quark production is one of the most important processes at hadron
colliders. Within the Standard Model, top quarks can be produced in
pairs (via strong interaction) or individually (via electroweak
processes). Top-antitop pair production is known to have the largest
cross section, and has been studied extensively, both experimentally
and theoretically. In fact, the top quark was discovered in 1995 at
the Tevatron in $t\bar{t}$ production events.  The observation of
single-top electroweak production is instead more difficult and it was
announced by the Tevatron experiments more
recently~\cite{Aaltonen:2009jj,Abazov:2009ii}.  This was due not only
to the fact that single-top cross section is smaller than the
$t\bar{t}$ one, but also to the presence of large backgrounds, namely
$W+\mbox{jet}$ and $t\bar{t}$, that required highly non-trivial
analysis strategies.

Although experimentally challenging, single-top electroweak production
is particularly important since it provides a relatively clean place
to study the electroweak properties of the top quark. For instance, it
allows a direct measurement of the $V_{tb}$ CKM matrix
element~\cite{Alwall:2006bx}, which is important in testing the
unitarity of the CKM matrix. Moreover, electroweak-produced top-quarks
are highly polarized, hence angular correlations of top-quark decay
products are particularly sizeable~\cite{Mahlon:1996pn,Mahlon:1999gz},
providing a good probe of the spin structure of the W-t-b
vertex. Finally, since the top-quark mass is close to the scale of
electroweak symmetry breaking, new physics effects could be discovered
in the the top quark sector, and reactions where only a single-top is
produced are particularly sensitive to some BSM models (see, for
example, refs.~\cite{Tait:2000sh,Cao:2007ea}). Single-top production
(in the $Wt$-channel) is also an important background for some
Higgs-boson search channels, such as $H\to W^+
W^-$~\cite{Dittmar:1996ss}.

The hadroproduction of single-top quarks in proton-proton collisions
is traditionally classified according to the partonic processes
present at LO:
the s-channel processes ($q \bar{q} \to t b$) involve the exchange of
a time-like $W$ boson, the $t$-channel processes ($b q \to t
q^{\prime}$) involve the exchange of a space-like $W$ boson, while
associated $Wt$ production ($b g \to t W^-$) involves the production
of a top quark in association with a $W$ boson. The $t$-channel
process is the main source of single-top quarks, both at the Tevatron
and the LHC. At the Tevatron the $Wt$ contribution is negligible and
the $s$-channel cross section is roughly half of the $t$-channel's
one. At the LHC, instead, the $Wt$-production cross section is a
factor 3 less than the $t$-channel's one, while the $s$-channel is
negligible. Therefore, an accurate description of all the three
production channels is important: in particular, at the LHC, the
$Wt$-channel will play a significant role.

Given the above reasons, it is then desirable to reach a high
precision in the theoretical predictions, both for total rates and for
more exclusive distributions.  Higher order corrections to single-top
hadroproduction have been calculated in
refs.~\cite{Bordes:1994ki,Giele:1995kr,Stelzer:1997ns,
  Harris:2002md,Zhu:2002uj,Campbell:2004ch,Campbell:2005bb,
  Cao:2004ap,Cao:2005pq,Kidonakis:2006bu,Campbell:2009ss,
  Campbell:2009gj}.  Nowadays, one of the ways to go beyond this level
of accuracy
is to merge the fixed order accuracy of a NLO calculation with the
(Next-to)-Leading-Logarithm accuracy (and the flexibility) of a
(parton shower) Monte Carlo event generator. At present there are two
methods to interface NLO calculations with parton showers in a
consistent way: \MCatNLO{}~\cite{Frixione:2002ik} and
\POWHEG{}~\cite{Nason:2004rx,Frixione:2007vw}.\footnote{In literature,
  other proposals
  exist~\cite{Nagy:2005aa,Giele:2007di,Lavesson:2008ah}. The
  \MCatNLO{} and \POWHEG{} methods are the only ones where, currently,
  several full processes at lepton and hadron colliders have been
  implemented.} They have been applied to several processes at
lepton~\cite{LatundeDada:2006gx,LatundeDada:2007jg,LatundeDada:2008bv}
and hadron~\cite{Frixione:2003ei,Frixione:2005vw,Nason:2006hf,
  Frixione:2007nw,Alioli:2008gx,Hamilton:2008pd,Frixione:2008yi,
  Alioli:2008tz,Hamilton:2009za,Papaefstathiou:2009sr,
  LatundeDada:2009rr,Alioli:2009je,Nason:2009ai,Weydert:2009vr,Torrielli:2010aw}
colliders, in conjunction with the \HERWIG{}~\cite{Corcella:2000bw},
\PYTHIA{}~\cite{Sjostrand:2006za} and \HWpp{}~\cite{Bahr:2008pv}
parton shower algorithms.\footnote{In~\cite{Butterworth:2010ym} a
  compared study of the effects of different parton shower algorithms
  has also been performed, in case of Higgs-boson production via gluon
  fusion.}\,\footnote{An implementation of the \POWHEG{} method has
  also appeared recently within the \SHERPA{} event
  generator~\cite{Hoeche:2010pf}.} Despite the theoretical formulation
of the two methods being quite different, reasonable agreement has
been found as well as the reason for the (few) differences encountered
(see for example~\cite{Nason:2010ap} for a recent discussion).

In ref.~\cite{Alioli:2009je} the details of the implementation of
single-top $s$- and $t$-channel production processes in the \POWHEG{}
framework have already been described.  The aim of this paper is to
describe the implementation of the $Wt$-channel production process.
Therefore, results presented here can be seen as a completion of the
work described in the previous reference.

This paper is organized as follows. In sec.~\ref{sec:description} we
describe some of the technicalities of the implementation. We also
summarize how the problem of the interference between the NLO
corrections to this process and the $t\bar{t}$ calculation has been
dealt with. Here we anticipate that we have used the same strategy
adopted by the \MCatNLO{} authors, i.e.~two definitions have been
implemented, such that their difference can be considered a measure of
one of the theoretical uncertainties that affects this process.  In
sec.~\ref{sec:results} results are presented: in particular we compare
\POWHEG{} and \MCatNLO{} results for the two aforementioned
definitions of the NLO corrections. We also show results obtained with
typical cuts and results obtained using the \PYTHIA{} shower. Finally,
in sec.~\ref{sec:conclusions} we give our conclusions.

\section{The \POWHEG{} implementation}
\label{sec:description}

To match single-top $Wt$-channel production NLO corrections with a
parton shower using the \POWHEG{} method, the \BOX{} package has been
used. \BOX{} is a program that automates all the steps described in
ref.~\cite{Frixione:2007vw}, turning a NLO calculation into a
\POWHEG{} simulation. The details of how the program works have been
largely described in ref.~\cite{Alioli:2010xd}. In practice, starting
from some of the typical building blocks of the NLO calculation for
the process at hand, the program produces a set of partonic events
ready to be showered by a shower Monte Carlo program.

In this section we present some of the inputs that we calculated and
then coded in the format needed by the \BOX{} package to work. In
particular, we describe how the LO kinematics and the relevant
differential cross sections up to next-to-leading order in the strong
coupling $\as$ have been obtained. All quark masses have been set to
zero (except, of course, the top-quark mass) and the full
Cabibbo-Kobayashi-Maskawa~(CKM) matrix has been taken into account in
the calculation. However, for sake of illustration, through this
section we set the CKM matrix equal to the identity. Therefore, the
$d$-type quark connected to the top and the $W$ boson will be denoted
as $b$.  Furthermore, in this paper we always refer to top-quark
production: anti-top production is obtained simply by charge
conjugation.

With this convention, the LO partonic process for $Wt$ production is
\begin{equation}
\label{eq:wtch_lo}
b + g \to W^-\! + t\,.
\end{equation}
We denote with $\mathcal{B}_{b g}$ the (summed and averaged) squared
amplitude, divided by the flux factor. The process initiated by the
partons $g b$ is also present. For brevity, throughout the paper
processes that can be obtained from the written ones by simply
exchanging the order of the incoming partons will be omitted from the
formulae.

Before giving the relevant formulae, here we stress again that the NLO
corrections to the $Wt$ production channel are not well defined, due
to interference effects with the $t\bar{t}$ process. This problem is
well known, and several approaches have been
introduced~\cite{Belyaev:1998dn,Tait:1999cf,Campbell:2005bb}. To deal
with this theoretical issue, we used the same strategy first described
in the corresponding \MCatNLO{} publication~\cite{Frixione:2008yi} and
later extensively studied in ref.~\cite{White:2009yt}: two definitions
for the NLO corrections have been considered, and the relative
difference between the two can be interpreted as a measure of the
theoretical uncertainty in the definition of genuine NLO corrections
to the $Wt$-channel process. Since this issue concerns the radiative
part of the QCD NLO corrections, a more accurate discussion will be
given in sec.~\ref{sec:real}.

\subsection{Born kinematics}
\label{sec:kinematics}
Following the notation of ref.~\cite{Frixione:2007vw}, we denote with
$\bar{k}_{\splus}$ and $\bar{k}_{\sminus}$ the incoming parton
momenta, aligned along the plus and minus direction of the $z$ axis,
and by $\bar{k}_1$ and $\bar{k}_2$ the outgoing $W$-boson and
top-quark momenta, respectively. The top-quark and $W$-boson masses
are denoted by $m_t$ and $m_W$. During the \POWHEG{}-algorithm step
where the inclusive NLO cross section is evaluated (the calculation of
the $\bar{B}$ function), their virtualities are kept fixed and equal
to their masses.  If $K_\splus$ and $K_\sminus$ are the momenta of the
incoming hadrons, then we have
\begin{equation}
\bar{k}_\splusminus=\bar{x}_\splusminus K_\splusminus \,,
\end{equation}
where $\bar{x}_\splusminus$ are the momentum fractions, and momentum
conservation reads
\begin{equation}
  \bar{k}_\splus + \bar{k}_\sminus = \bar{k}_1 + \bar{k}_2\,.
\end{equation}
We now introduce the variables
\begin{equation}
  \bar{s} = (\bar{k}_\splus+\bar{k}_\sminus)^2 ,\qquad\quad \bar{Y} =
  \frac{1}{2}
  \log \frac{(\bar{k}_\splus + \bar{k}_\sminus)^0 + (\bar{k}_\splus +
  \bar{k}_\sminus)^3} {(\bar{k}_\splus + \bar{k}_\sminus)^0 -
  (\bar{k}_\splus + \bar{k}_\sminus)^3}\,,
\end{equation}
and $\bar{\theta}$, the angle between the outgoing top quark and the
$\bar{k}_{\splus}$ momentum, as seen in the partonic
center-of-mass~(CM) frame. We denote with $\bar{\phi}$ the azimuthal
angle of the outgoing top quark in the same reference frame.  Since
the differential cross sections do not depend on the overall azimuthal
orientation of the outgoing partons, we set this angle to zero.  At
the end of the generation of an event, we perform a uniform, random
azimuthal rotation of the whole event, in order to cover the whole
final-state phase space.  The set of variables
$\tmmathbf{\bar{\Phi}}_2 \equiv\lg
\bar{s},\bar{Y},\bar{\theta},\bar{\phi}\rg$ fully parameterizes the
Born kinematics. From them, we can reconstruct the momentum fractions
\begin{equation}
  \label{eq:xpxm_bar}
  \bar{x}_{\splus}  = \sqrt{\frac{\bar{s}}{S}}\, e^{\bar{Y}}, \qquad\quad
  \bar{x}_{\sminus} = \sqrt{\frac{\bar{s}}{S}}\, e^{- \bar{Y}}\,,
\end{equation}
where $S=(\Kplus + \Kminus)^2$ is the squared CM energy of
the hadronic collider.  The outgoing momenta are first reconstructed
in their longitudinal rest frame, where $\bar{Y}=0$. In this frame, their
energies are
\begin{equation}
  \bar{k}_1^0 |_{\sss \bar{Y}=0}=\frac{\bar{s}+m_W^2-m_t^2}{2\sqrt{\bar{s}}}
  \qquad\mbox{and}\qquad
  \bar{k}_2^0 |_{\sss \bar{Y}=0}=\frac{\bar{s}+m_t^2-m_W^2}{2\sqrt{\bar{s}}}\,.
\end{equation}
The two spatial momenta are opposite and their modulus
$|\vec{k}|_{\sss \bar{Y}=0}$ is found by using the on-shell condition
$m_W^2=(\bar{k}_1^0|_{\sss \bar{Y}=0})^{2} - (|\vec{k}|_{\sss
  \bar{Y}=0})^2$.  We fix the top-quark momentum to form an angle
$\bar{\theta}$ with the $\splus$ direction and to have zero azimuth
(i.e.~it lies in the $x z$ plane and has positive $x$ component). Both
$\bar{k}_1$ and $\bar{k}_2$ are then boosted back in the laboratory
frame, with boost rapidity $\bar{Y}$.  The Born phase space, in terms
of these variables, can be written as
\begin{eqnarray}
\label{eq:Bphsp}
d \tmmathbf{\bar{\Phi}}_2 &=& d\bar{x}_\splus \, d\bar{x}_\sminus
(2\pi)^4 \delta^4\!\(\bar{k}_\splus+ \bar{k}_\sminus -\bar{k}_1 -\bar{k}_2\)
\frac{d^3 \bar{k}_1}{(2\pi)^3 2 \bar{k}_1^0} \,
\frac{d^3 \bar{k}_2}{(2\pi)^3 2 \bar{k}_2^0} \nonumber\\
&=&\frac{1}{S}\frac{\beta}{16\pi}\ d\bar{s}\ d\bar{Y}\ d\cos{\bar{\theta}}\
\frac{d\bar{\phi}}{2\pi}\,,
\end{eqnarray}
where
\begin{equation}
  \beta=\sqrt{1-\rho}\,,\qquad\mbox{with}\qquad \rho=\frac{2\, (m_W^2+m_t^2)}{\bar{s}} - \frac{(m_W^2-m_t^2)^2}{\bar{s}^2}\,.
\end{equation}
The content of eq.~(\ref{eq:Bphsp}) and the procedure described above
to define the parton momenta from the variables set
$\tmmathbf{\bar{\Phi}}_2$ have been coded in the subroutine
\tmtexttt{born\_phsp} of the \BOX{} package~\cite{Alioli:2010xd}.

\subsection{Born and virtual contributions}
\label{sec:bornvirtual}
The squared matrix element for the lowest order subprocess $b g \to
W^- t$ has been obtained with MadGraph~\cite{Alwall:2007st}. It has
been checked with the expression reported in eq.~(3.5) of
ref.~\cite{Frixione:2008yi} and with an independent calculation
performed with FeynCalc~\cite{Mertig:1990an}, starting from the two
Feynman diagrams in fig.~\ref{fig:wtch_lo}. The top-quark width
$\Gamma_t$ has been set here to zero.
\begin{figure}[htb]
  \begin{center}
    \epsfig{file=JD_bg2tw,width=0.4\textwidth}
  \end{center}
  \caption{\label{fig:wtch_lo} Feynman diagrams for the LO process $b
    g \to W^- t$. The top quark is denoted with a double line.}
\end{figure}

To regularize soft and collinear divergences, the \BOX{} package uses
an automated implementation of the subtraction algorithm proposed by
\FKS{}~\cite{Frixione:1995ms,Frixione:1997np}.  The counterterms
needed for the numerical subtraction of the real squared amplitudes
soft divergences are obtained from the color-linked Born squared
amplitudes $\mathcal{B}_{ij}$ defined in eq.~(2.97) of
ref.~\cite{Frixione:2007vw}, using the eikonal approximation.  For the
process at hand, the $\mathcal{B}_{ij}$ matrix is proportional to the
full Born squared amplitude. In fact, since at the lowest order we
have only three color-connected partons, the equality
\begin{equation}
\label{eq:color_cons}
  {\bf T}_\splus + {\bf T}_\sminus +{\bf T}_2 = 0\,
\end{equation}
holds, where ${\bf T}_i$ is the color operator associated with the $i$
parton in the Born process~(\ref{eq:wtch_lo}).  Using the property
\begin{equation}
  {\bf T}_i\mydot{\bf T}_j =
  \frac{({\bf T}_i + {\bf T}_j)^2 - {\bf T}_i^2 -{\bf T}_j^2}{2}\,,
\end{equation}
and the fact that ${\bf T}_i^2=\CF$ if $i$ is a quark and ${\bf
  T}_i^2=\CA$ if $i$ is a gluon, we obtain ${\matB}_{ij} = (2\CF-\CA)
\matB_{bg}/2$ if neither $i$ or $j$ are gluons, and ${\matB}_{ij} =
\CA \matB_{bg}/2$ if $i$ or $j$ is a gluon.

Collinear counterterms are obtained using the collinear factorization.
However, since in the $Wt$-channel LO process~(\ref{eq:wtch_lo}) an
external gluon is present, the collinear limits associated with this
leg
do not factorize in terms of the Altarelli-Parisi unpolarized
splitting kernels times the Born contribution $\mathcal{B}_{bg}$. In
fact azimuthal correlations in the branching process are present, and
to build a local counterterm the \BOX{} package makes use of the spin
correlated Born cross sections $\mathcal{B}_{\mu\nu}$, defined in
eq.~(2.8) of ref.~\cite{Frixione:2007vw}. The FeynCalc program has
been used to calculate this matrix and to translate the result in a
Fortran routine.

One loop virtual contributions have been calculated and algebraically
reduced to scalar integrals using the Passarino-Veltman algorithm. We
used the same renormalization procedure described in
ref.~\cite{Frixione:2008yi}. We checked that the IR pole structure we
are left with corresponds to the singularities of the real
contributions. The finite part that enters in the soft-virtual
contribution $\mathcal{V}_{bg}$ (eq. (4.6) of
ref.\cite{Alioli:2010xd}) is computed numerically, using the package
QCDloop~\cite{Ellis:2007qk} to evaluate the finite part of the scalar
integrals.

\subsection{Real contributions}
\label{sec:real}
The real emission corrections can be classified as follows:
\begin{eqnarray}
  \label{eq:wtch_r_bg}
  b + g &\to& W^-\! + t + g\,,\\
  \label{eq:wtch_r_bq}
  b + q(\bar{q}) &\to& W^-\! + t + q(\bar{q})\,\qquad (q\ne b),\\
  \label{eq:wtch_r_gg}
  g + g &\to& W^-\! + t + \bar{b}\,,\\
  \label{eq:wtch_r_qq}
  q + \bar{q} &\to& W^-\! + t + \bar{b}\,.
\end{eqnarray}
We denote the corresponding contributions to the cross section as
$\mathcal{R}_{b g}$, $\mathcal{R}_{b q}$, $\mathcal{R}_{g g}$ and
$\mathcal{R}_{q \bar{q}}$ respectively, where we used again the
standard \POWHEG{} notation first introduced in
ref.~\cite{Frixione:2007vw}. The 3-body phase space is denoted as
$\tmmathbf{\Phi}_3$ and the corresponding momenta as $k_\splus$,
$k_\sminus$, $k_1$, $k_2$ and $k_3$, where $k_3$ is the momentum of
the outgoing light parton (the FKS parton) while the other momenta
correspond to those of the two incoming partons, the $W$-boson and the
top-quark. In fig.~\ref{fig:wtch_nlo} some representative diagrams are
shown.
\begin{figure}[htb]
  \begin{center}
    \epsfig{file=JD_tw_real_bg_bq,width=0.5\textwidth}\\
    \epsfig{file=JD_tw_real_gg,width=0.5\textwidth}\\
    \epsfig{file=JD_tw_real_qq,width=0.5\textwidth}
  \end{center}
  \caption{\label{fig:wtch_nlo} Representative Feynman diagrams for
    the processes~(\ref{eq:wtch_r_bg}),~(\ref{eq:wtch_r_bq})
    (up),~(\ref{eq:wtch_r_gg}) (center) and~(\ref{eq:wtch_r_qq})
    (bottom).}
\end{figure}

The processes~(\ref{eq:wtch_r_gg}) and~(\ref{eq:wtch_r_qq}) have a
final state that corresponds to $t\bar{t}$ production followed by the
decay $\bar{t}\to W^- \bar{b}$. A consequence of this fact is the well
known problem of interference between $Wt$ and $t\bar{t}$ production
and it is the reason why QCD NLO corrections to the $Wt$-channel are
not well defined. In the following we will state more precisely the
nature of the problem and explain how we dealt with it.

As one can see from the last four Feynman diagrams in
fig.~\ref{fig:wtch_nlo}, diagrams associated with the
subprocesses~(\ref{eq:wtch_r_gg}) and~(\ref{eq:wtch_r_qq}) can be
divided into two sets. Following the nomenclature of
ref.~\cite{Frixione:2008yi}, we call \emph{doubly-resonant} the
diagrams where a top-quark propagator goes on-shell when the invariant
mass of the system made by the $W$-boson and the outgoing
$\bar{b}$-quark ($m_{W\bar{b}}$) approaches $m_t$. The other diagrams
we call \emph{singly-resonant}. Since these two sets of diagrams have
to be summed at the amplitude level, interference effects are
present. This interference between NLO real corrections to $Wt$ and
lowest order $t\bar{t}$ production (followed by a decay) would not be
a problem if the contamination on $t\bar{t}$-like diagrams was
numerically negligible with respect to the size of singly-resonant
diagrams.  This is certainly not the case when one approaches the
region $(k_W+k_{\,\bar{b}})^2\equiv (k_1+k_3)^2\sim m_t^2$. In fact,
in this region the real corrections to doubly-resonant diagrams become
huge (an internal propagator goes on-shell), and the perturbative
expansion (in power of $\as$) for the NLO corrections to single-top
$Wt$-channel loses its meaning.

To deal with this problem, several approaches are possible:
\begin{itemize}
\item The more drastic approach is to consider that top quarks are not
  detectable particles. This approach would remove any interference
  issue, since processes would be unambiguously classified accordingly
  only to experimentally measurable (QCD) final states: one would have
  the processes $W^+W^-b\bar{b}$ and $W^-W^+b$. At present, the price
  to pay would be to neglect NLO corrections, since these are not
  available for processes where top-quarks are not
  on-shell.\footnote{We also recall that a calculation of NLO
    corrections to the production and the decay of top quarks
    performed in the narrow width approximation (as the one appeared
    in ref.~\cite{Campbell:2005bb}) would not avoid the interference
    problem. To our knowledge, the only calculations that fully
    include top offshellness effects are the ones described in
    ref.~\cite{Pittau:1996rp,Falgari:2010sf} for single-top $s$- and
    $t$-channel production respectively.}
\item An alternative approach is to give a prescription for removing
  the contamination from $t\bar{t}$-like contributions, while keeping
  the top-quark as a final state particle. This task can be
  accomplished using cuts to avoid the doubly-resonant region or
  including an \emph{ad~hoc} subtraction term to remove (or suppress)
  the $t\bar{t}$ contribution.

  The cut strategy was first adopted in ref.~\cite{Belyaev:1998dn},
  where an explicit cut on $m_{W\bar{b}}$ was used. In
  ref.~\cite{Campbell:2005bb} a $b$-jet veto was instead used,
  together with a careful choice of the factorization scale. Moreover,
  in ref.~\cite{Campbell:2005bb} for the first time differential
  results at the NLO were presented and QCD corrections to the
  top-quark decay have been also included, in the narrow-width
  approximation.

  In ref.~\cite{Tait:1999cf} a subtraction term was first introduced
  in the context of an inclusive LO(+LL) calculation. A similar method
  has been adopted in the fully inclusive NLO calculation presented in
  ref.~\cite{Zhu:2002uj}.

  Finally, in refs.~\cite{Frixione:2008yi,White:2009yt} the
  interference problem was reexamined at length, and it was shown that
  a separate treatment of $t\bar{t}$ and $Wt$ production is feasible,
  also when the NLO calculation is matched with a parton shower.
\end{itemize}

In this work we have used the same approach described by the
\MCatNLO{} authors. This relies on the observation that a meaningful
definition of the $Wt$-channel process (as a signal or a background)
is possible only through cuts on final state objects.  If interference
effects with $t\bar{t}$ are negligible \emph{after} these cuts are
applied, then it is possible to consider $Wt$-channel a well defined
process. Since cuts act differently in different phase-space regions,
one needs to quantify the interference between $Wt$ and $t\bar{t}$
locally in the phase-space. To this end, two definitions for the NLO
corrections were given. It was shown that by comparing the two results
it is possible to address the previous question, and it was found
that, for several sets of cuts, the theoretical uncertainty due to
interference effects is typically negligible with respect to other
theoretical errors.  In the following we will discuss briefly the two
subtraction strategies and how they have been implemented within the
\POWHEG{} framework. In sec.~\ref{sec:results} the corresponding
results will be shown.

The two definitions proposed in ref.~\cite{Frixione:2008yi} are known
as \emph{Diagram Removal} (DR) and \emph{Diagram Subtraction} (DS).
Their difference can be better understood by writing a generic
amplitude for the processes~(\ref{eq:wtch_r_gg})
and~(\ref{eq:wtch_r_qq}) as
\begin{equation}
  \label{eq:ampl_wt_tt}
  \mathcal{M}=\mathcal{M}^{Wt} + \mathcal{M}^{t\bar{t}}\,,
\end{equation}
where $\mathcal{M}^{Wt}$ and $\mathcal{M}^{t\bar{t}}$ denote
respectively the sum of all the singly- and doubly-resonant Feynman
diagrams for the partonic subprocess at hand.  In DR one defines the
real contribution $\mathcal{R}$ by eliminating the $t\bar{t}$
contribution $\mathcal{M}^{t\bar{t}}$ from $\mathcal{M}$ before
squaring the amplitude. Instead in DS one keeps the full squared
amplitude but subtracts from it a local counterterm
$\mathcal{C}^{\sss\rm SUB}$ in order to suppress the $t\bar{t}$
contribution at the cross section level.  In this respect, DS can be
seen as a refinement of the method proposed in ref.\cite{Tait:1999cf}.
Schematically, we have:
\begin{eqnarray}
  \mathcal{R}^{\sss\rm DR} &=& \frac{|\mathcal{M}^{Wt}|^2}{2s}\,,  \label{eq:dr}\\
  \mathcal{R}^{\sss\rm DS} &=& \frac{|\mathcal{M}^{Wt}+\mathcal{M}^{t\bar{t}}|^2
  - \mathcal{C}^{\sss\rm SUB}}{2s}\,,  \label{eq:ds}
\end{eqnarray}
where $s$ is the squared CM energy.  Some comments are due here:
\begin{itemize}
\item The DR method is not gauge invariant. This issue was studied in
  depth by the authors of ref.\cite{Frixione:2008yi}, and it was shown
  that the impact of gauge dependence in the DR calculation is
  numerically negligible.
\item In the DS approach one wants to build a gauge invariant
  subtraction term that exactly cancels the $t\bar{t}$ contribution
  when the doubly-resonant region is approached. Thus, the subtraction
  term $\mathcal{C}^{\sss\rm SUB}$ has to fulfill the following
  requirements:
  \begin{enumerate}
  \item gauge invariance.
  \item match exactly the $|\mathcal{M}^{t\bar{t}}|^2$ contribution
    when the doubly-resonant region is approached.
  \item possibly fall off (quickly) far from the doubly-resonant
    region.
  \end{enumerate}
  The third requirement is needed to keep the full NLO corrections
  unmodified away from the $t\bar{t}$ peak.  Apart from the three
  requirements above, there is some freedom in the definition of
  $\mathcal{C}^{\sss\rm SUB}$.
\item By taking the difference between eq.~(\ref{eq:ds})
  and~(\ref{eq:dr}), one finds
  \begin{equation}
    \mathcal{R}^{\sss\rm DS}-\mathcal{R}^{\sss\rm DR}=
    \frac{\mathcal{I} + |\mathcal{M}^{t\bar{t}}|^2 -
      \mathcal{C}^{\sss\rm SUB}}{2 s}\,,
  \end{equation}
  where
  $\mathcal{I}=2\,\Re{(\mathcal{M}^{Wt}\,{\mathcal{M}^{t\bar{t}}}^*)}$.
  Therefore, the difference between results obtained with DR and DS
  can be interpreted as a measure of the size of the interference
  $\mathcal{I}$, provided that the difference
  $|\mathcal{M}^{t\bar{t}}|^2 - \mathcal{C}^{\sss\rm SUB}$ is small.
\end{itemize}

We implemented in \POWHEG{} the two subtraction methods. The squared
amplitudes $|\mathcal{M}|^2$ and $|\mathcal{M}^{Wt}|^2$ have been
obtained using MadGraph. The subtraction term was chosen as in
ref.~\cite{Frixione:2008yi}:
\begin{equation}
  \mathcal{C}^{\sss\rm SUB}(\tmmathbf{\Phi}_3)=\frac{(m_t\Gamma_t)^2}
  {((k_1+k_2)^2-m_t^2)^2 + (m_t\Gamma_t)^2}
  \ {|\mathcal{M}^{t\bar{t}}(\tmmathbf{\Phi}_3^{\prime})|^2}\,,
\end{equation}
where $\tmmathbf{\Phi}_3^{\prime}$ is a point in the 3-body phase
space obtained by reshuffling the $\tmmathbf{\Phi}_3$ kinematics in
order to have $(k_1+k_3)^2=m_t^2$, i.e.~an exactly doubly-resonant
configuration. Hence, the choice of the subtraction term is the same
as the one used by the \MCatNLO{} authors. In fact, in spite of the
aforementioned freedom, as it was already pointed out in
ref.~\cite{Frixione:2008yi}, the choice of the amplitude
$\mathcal{M}^{t\bar{t}}$ evaluated at the point
$\tmmathbf{\Phi}_3^{\prime}$ is unavoidable if one wants to achieve
the exact cancellation of the doubly-resonant contribution while
retaining gauge invariance, since
$\mathcal{M}\sim\mathcal{M}^{t\bar{t}}$ when $(k_1+k_3)^2\to m_t^2$
and gauge invariance is preserved only if the internal $\bar{t}$
propagator is on shell.  The only real freedom is in the choice of the
prefactor, and the Breit-Wigner profile seems the more natural choice
if one wants the difference between $|\mathcal{M}^{t\bar{t}}|^2$ and
$\mathcal{C}^{\sss\rm SUB}$ to be close to zero as much as possible
away from the resonance.

The inclusion of DR in \POWHEG{} is straightforward. The procedure we
adopted for DS is more subtle, since $\mathcal{R}^{\sss\rm DS}$ is not
positive-definite. This affects two steps of the \POWHEG{} method. In
the following we describe how we proceeded in our implementation.  As
usual, we use the standard notation of ref.~\cite{Frixione:2007vw}.
\begin{itemize}
\item Having a real correction that is not always positive-definite
  increases the chances to have regions of $\tmmathbf{\bar{\Phi}}_2$
  where the $\tilde{B}$ function becomes negative. In these cases, the
  radiative event generated with \POWHEG{} starting from the
  underlying-Born configuration $\tmmathbf{\bar{\Phi}}_2$ is
  negative-weighted.  We have checked that this occurrence is rare, so
  that the benefits of having positive-weighted events are not
  spoilt. Moreover, in the \BOX{} package, a procedure (called
  ``folded'' integration) to reduce further the occurrence of these
  events is available, as explained in sec. 4.1 of
  ref.~\cite{Alioli:2010xd}. By using it, we have verified that the
  occurrence of negative-weighted events can be further reduced also
  for the $Wt$-channel case. This also implies that the function
  $\bar{B}$ is positive, as in the other \POWHEG{} implementations.

  Quantitatively, in the event sample we generated to produce the
  plots shown in sec.~\ref{sec:results}, the fraction of
  negative-weighted events was $0.05$. By using the aforementioned
  folded integration, this fraction can be reduced down to $0.005$,
  becoming therefore completely negligible.

\item For 3-body kinematic configurations where $\mathcal{R}^{\sss\rm
    DS}$ is not positive, another problem can occur during the
  generation of the hardest radiation. We recall that, to generate the
  hardest radiation, the \POWHEG{} algorithm works by finding an
  upper-bound for the ratio $\mathcal{R}/\mathcal{B}$, assuming that
  this ratio is always positive. Negative values can therefore spoil
  the accuracy of the method. We have checked that this happens with a
  certain frequency only close to the doubly-resonant region. Since
  this is the kinematic region where the separation of $t\bar{t}$ and
  $Wt$ is already particularly critical,

  we decided to explicitly avoid to generate radiative events when
  close to the doubly-resonant peak.  This has been obtained using a
  theta function that vanishes in this region, i.e.~in DS the
  following substitution is performed when $\mathcal{R}/\mathcal{B}$
  is evaluated:
  \begin{equation}
  \mathcal{R}^{\sss\rm DS} \to \mathcal{R}^{\sss\rm DS} \times \
  \theta\(|m_{W\bar{b}}-m_t| - \kappa \Gamma_t\) \,.
  \end{equation}
  We have tested the effects of this cutoff trying values of order 1
  for $\kappa$, in order to avoid introducing effects from this
  parameter in phase-space regions where it is not needed. The outcome
  of this check is that, at the end of the event generation
  (i.e.~after the shower and the hadronization stage), no problems
  caused by this cutoff are present, since results depend negligibly
  on the value of $\kappa$. The only observable where some dependence
  was observed is the differential distribution of $m_{W\bar{b}}$
  close to the doubly-resonant peak. When values for $\kappa$ in the
  region $3-5$ are used, the difference with the \MCatNLO{} result is
  minimized.  The results shown in sec.~\ref{sec:results} have been
  obtained with $\kappa=3$.

  Although we are aware that this solution may be considered
  unappealing, we found it a reasonable choice to handle with the
  problem of having negative values for $\mathcal{R}/\mathcal{B}$ in
  the \POWHEG{} framework. We also recall that the presence of this
  cutoff in our implementation does not affect the issue of negative
  weights, since $\kappa$ is used only in the generation of the
  hardest radiation.
\end{itemize}

We used MadGraph to calculate the term $\mathcal{C}^{\sss\rm SUB}$. We
also notice that to avoid the divergence of the internal top-quark
propagator when $(k_1+k_3)^2\to m_t^2$, in DS a non-vanishing value for
the top width is needed.\footnote{We recall that to obtain the proper
  cancellation of soft and collinear divergences, $\mathcal{M}^{Wt}$
  has to be calculated with $\Gamma_t=0$.}

\section{Results}
\label{sec:results}
In this section we present our results obtained after showering with
\HERWIG{}~6.510 and \PYTHIA{}~6.4.22 the partonic events generated
with \POWHEG{}. We considered top production at the LHC, with an
hadronic center-of-mass energy $\sqrt{S}=10 $ TeV.  All results have
been obtained assuming that the top-quark decays semileptonically
($t\to b\,\bar{\ell}\,\nu$) and that the $W$-boson involved in the
hard scattering decays leptonically ($W^-\to
\ell\,\bar{\nu}$). Branching ratios have been removed, so that plots
are normalized to the total cross section.

We have used the CTEQ6M~\cite{Pumplin:2002vw} set for the parton
distribution functions and the associated value of
\mbox{$\Lambda_{\scriptscriptstyle\overline{\rm
      MS}}^{(5)}=0.226$~GeV}.  Furthermore, as discussed in
refs.~\cite{Frixione:2007vw,Nason:2006hf}, we use a rescaled value
\mbox{$\Lambda_{\scriptscriptstyle\rm{MC}}=
  1.569\,\Lambda_{\scriptscriptstyle\overline{\rm MS}}^{(5)}$} in the
expression for $\as$ appearing in the Sudakov form factors, in order
to achieve next-to-leading logarithmic accuracy.

Although the matrix-element calculation has been performed in the
massless-quark limit,
the lower cutoff
in the generation of the radiation has been fixed according to the mass of
the emitting quark. The lower bound on the transverse momentum for the
emission off a massless emitter ($u$, $d$, $s$) has been set to the value
$\ptmin = \sqrt{5}\,\Lambda_{\scriptscriptstyle\rm{MC}}$. We instead choose
$\ptmin$ equal to $m_c$ or $m_b$ when the gluon is emitted by a charm or a
bottom quark, respectively.  We set $m_c=1.55$~GeV and $m_b= 4.95$~GeV.

The renormalization and factorization scales have been taken equal to
the transverse momentum of the radiated light parton during the
generation of radiation, in accordance with the \POWHEG{}
formalism. We have also taken into account properly the heavy-flavour
thresholds in the running of $\as$ and in the PDF's, by changing the
number of active flavours when the renormalization or factorization
scales cross a mass threshold.  In the $\bar{B}$ calculation, instead,
$\mu_R$ and $\mu_F$ have been chosen equal to the top-quark mass,
whose value has been fixed to $m_t=175$~GeV. In the DS approach, the
amplitudes where doubly-resonant graphs are present and the
subtraction term $\mathcal{C}^{\sss\rm SUB}$ have been calculated with
$\Gamma_t=1.7$ GeV.

To assess the validity of the approximations and the choices we made,
we compare our results (obtained both with DR and DS) with the
\MCatNLO{} outputs.

\begin{figure}[htb]
  \begin{center}
    \epsfig{file=cmp_ptt_etaw_wt_DR_LHC10-POW-MCNLO,width=\figwidth}\\
    \epsfig{file=cmp_pttw_phitw_wt_DR_LHC10-POW-MCNLO,width=\figwidth}
  \end{center}
  \captskip
  \caption{\label{fig:cmp1_DR-POW-MCNLO} Comparisons between \POWHEG{}
    (interfaced to \HERWIG{}) and \MCatNLO{} results at the LHC $pp$
    collider ($\sqrt{S}=10$ TeV), obtained with the DR prescription.
    NLO results are also shown in the lower panel.}
\end{figure}
In fig.~\ref{fig:cmp1_DR-POW-MCNLO} we show a comparison between
\POWHEG{} and \MCatNLO{}, obtained with the DR prescription, without
applying any cut on the final state particles. \POWHEG{} results have
been obtained using the \HERWIG{} parton shower, in order to minimize
differences arising from different shower algorithms and hadronization
models. In the upper panel of fig.~\ref{fig:cmp1_DR-POW-MCNLO} we show
the transverse momentum of the top quark ($\pt^{\,t}$) and the
pseudorapidity of the $W$-boson ($\eta^{\,W}$) produced in the hard
process (i.e.~not the $W$-boson present in the decay chain of the top
quark). As expected, we found very good agreement between the two
results, since the shape of these two distributions is due mainly to
the fixed order result. In the lower panel, we show instead $\pt^{(t
  W\!)}$, the transverse momentum of the system made by the top quark
and the $W$ boson, and $\Delta\phi_{t\mbox{-}W}$, the difference
between the azimuthal angles of the two particles. These two
quantities are significant because their shape is affected by Sudakov
suppression effects due to the resummation performed by parton
showers.  To stress the size of these effects, we also superimposed
the fixed order (NLO) prediction to the last two plots.  We observe
good agreement between showered results, and the expected difference
with the NLO curve, where the cancellation of soft and collinear
divergences takes place only at the edge of the distributions (at
$\pt^{(t W\!)}=0$ and $\Delta\phi_{t\mbox{-}W}=\pi$ respectively).
\begin{figure}[htb]
  \begin{center}
    \epsfig{file=cmp_ptt_etaw_wt_DS_LHC10-POW-MCNLO,width=\figwidth}\\
    \epsfig{file=cmp_pttw_phitw_wt_DS_LHC10-POW-MCNLO,width=\figwidth}
  \end{center}
  \captskip
  \caption{\label{fig:cmp1_DS-POW-MCNLO} Comparisons between \POWHEG{}
    (interfaced to \HERWIG{}) and \MCatNLO{} results at the LHC $pp$
    collider ($\sqrt{S}=10$ TeV), obtained with the DS prescription.
    NLO results are also shown in the lower panel.}
\end{figure}

In fig.~\ref{fig:cmp1_DS-POW-MCNLO} we show the same set of plots,
obtained with the DS prescription. The same considerations made above
are valid also for this case. We recall that the plots shown in
figs.~\ref{fig:cmp1_DR-POW-MCNLO} and~\ref{fig:cmp1_DS-POW-MCNLO} have
been obtained using the top-quark and the $W$-boson momenta extracted
from the parton shower history. Therefore, these quantities are not
measurable in a real detector. Nevertheless, since they are useful to
check the implementation, we have shown the corresponding results.

We have also included the generation of the top-quark and the
$W$-boson decay products, according to the method originally proposed
in ref.~\cite{Frixione:2007zp}. This enables the generation of events
in which spin correlation effects in the production-decay stage are
taken into account with leading-order accuracy.  In
figs.~\ref{fig:cmp3_DR-POW-MCNLO} and~\ref{fig:cmp3_DS-POW-MCNLO} we
show a comparison between \MCatNLO{} and \POWHEG{} for some relevant
leptonic distributions. We plot the transverse momenta of the hardest
($\pt^{\ell_1}$) and the second hardest ($\pt^{\ell_2}$) charged
lepton in the event. We also show $\pt^{(\ell_1 \ell_2\!)}$, the
transverse momentum of the system made by $\ell_1$ and $\ell_2$, and
$\Delta\phi_{\ell_1\mbox{-}\ell_2}$, the difference between the
azimuthal angles of the two leptons. Here, again, we found very good
agreement in both the DR and the DS case.  We recall that a quantity
like $\Delta\phi_{\ell_1\mbox{-}\ell_2}$ is sensitive to
spin-correlation effects, as we have verified by running the code with
the decay-generation procedure switched off and letting the shower
perform isotropic decays.
\begin{figure}[htb]
  \begin{center}
    \epsfig{file=cmp_ptl1_ptl2_wt_DR_LHC10-POW-MCNLO,width=\figwidth}\\
    \epsfig{file=cmp_ptl1l2_phil1l2_wt_DR_LHC10-POW-MCNLO,width=\figwidth}
  \end{center}
  \captskip
  \caption{\label{fig:cmp3_DR-POW-MCNLO} Comparisons between \POWHEG{}
    (interfaced to \HERWIG{}) and \MCatNLO{} results at the LHC $pp$
    collider ($\sqrt{S}=10$ TeV), obtained with the DR prescription,
    for leptonic quantities.}
\end{figure}
\begin{figure}[htb]
  \begin{center}
    \epsfig{file=cmp_ptl1_ptl2_wt_DS_LHC10-POW-MCNLO,width=\figwidth}\\
    \epsfig{file=cmp_ptl1l2_phil1l2_wt_DS_LHC10-POW-MCNLO,width=\figwidth}
  \end{center}
  \captskip
  \caption{\label{fig:cmp3_DS-POW-MCNLO} Comparisons between \POWHEG{}
    (interfaced to \HERWIG{}) and \MCatNLO{} results at the LHC $pp$
    collider ($\sqrt{S}=10$ TeV), obtained with the DS prescription,
    for leptonic quantities.}
\end{figure}

In ref.~\cite{Campbell:2005bb} the use of a $b$-veto as a method to
discriminate between the $Wt$ and the $t\bar{t}$ processes was
proposed. This idea was also reconsidered in
ref.~\cite{Frixione:2008yi}, where its applicability was studied in
the context of a NLO result merged with a parton shower.  We have
performed a similar exercise using the results obtained with \POWHEG{}
and showered with \HERWIG{}. The $b$-veto condition is defined as
follows: after having sorted in decreasing $\pt$ all the $b$-flavoured
hadrons present in the final state, an event is kept if the
second-hardest $b$-flavoured hadron among those with central
pseudorapidity ($\eta^B\le 2.5$) has $\pt^B<\pt^{(b\rm\mbox{-}veto)}$,
otherwise it is discarded.\footnote{If all the $b$-hadrons in the
  central pseudorapidity region have $\pt^B<\pt^{(b\rm\mbox{-}veto)}$,
  or if the second-hardest $b$-hadron has a large pseudorapidity, the
  event is kept.} In our analysis procedure we chose
$\pt^{(b\rm\mbox{-}veto)}=50$ GeV.

As expected, we observe that the difference between DR and DS results
is reduced when the $b$-veto is applied: in fact, although we have not
performed a detailed study including also uncertainties due to scale
variations, we found that the ratios between the total NLO cross
sections using DR and DS are as follows:
\begin{equation}
  \frac{\sigma_{(\rm\sss DS)}}{\sigma_{(\rm\sss DR)}}= 0.95\,, \ \
  \frac{\sigma_{(\rm\sss DS)}^{(b\rm\mbox{-}veto)}}{\sigma_{(\rm\sss DR)}^{(b\rm\mbox{-}veto)}}= 0.98\,.
\end{equation}

As the \MCatNLO{} authors already pointed out, despite the fact that
the $b$-veto reduces the difference between DR and DS total cross
sections, it is not guaranteed that all the differential distributions
are affected in the same way. To address this question, in
fig.~\ref{fig:cmpbveto_POWHER-DR-DS} we show the effect of the
$b$-veto cut for two transverse-momentum spectra. In the upper panel
we plot the results for $\pt^{(t W\!)}$, obtained with the DR and the
DS procedures, before and after imposing the $b$-veto. The curves in
the upper-left panel are obtained without cuts, those on the
upper-right have been obtained keeping only the events that fulfil the
veto condition. It can be seen that the mismatch between DR and DS in
the high-$\pt$ tail is less sizeable when the $b$-veto is applied,
which is the expected result. Since in ref.~\cite{Frixione:2008yi} it
was noticed that the transverse momentum of the system made by the two
hardest leptons turns out to be sensitive to the treatment of the
doubly-resonant region, in the lower panel of
fig.~\ref{fig:cmpbveto_POWHER-DR-DS} we show also the predictions for
$\pt^{(\ell_1 \ell_2\!)}$.  Also for this quantity the effect of the
$b$-veto is to reduce the differences between DR and DS, as the plot
in the lower-right panel shows, in accordance with what has been found
in ref.~\cite{Frixione:2008yi}.
\begin{figure}[htb]
  \begin{center}
    \epsfig{file=cmp_pttw_wt_LHC10_POWHER-DR-DS_bveto,width=\figwidth}\\
    \epsfig{file=cmp_ptl1l2_wt_LHC10_POWHER-DR-DS_bveto,width=\figwidth}
  \end{center}
  \captskip
  \caption{\label{fig:cmpbveto_POWHER-DR-DS} Comparisons between
    \POWHEG{} (interfaced to \HERWIG{}) results obtained with the DR
    and DS prescriptions.  Plots in the right panel have been obtained
    with the $b$-veto described in the text.}
\end{figure}

The conclusions of this analysis can be summarized as follows. We
checked that the $b$-veto reduces the differences between DR and
DS. This is the expected result, since the $b$-veto was originally
proposed by the authors of ref.~\cite{Campbell:2005bb} to reduce
interference effects from $t\bar{t}$ production. In fact, by requiring
a $b$-veto, one reduces the number of events with two hard and central
$b$-flavoured hadrons, which is indeed one of the typical signatures
of $t\bar{t}$ production.  The fact that DR and DS total cross
sections become closer when the $b$-veto is applied is also in
accordance with the interpretation of the difference between DR and DS
being a measure of the interference between $Wt$ and $t\bar{t}$
production. This interpretation is supported by the plots in
fig.~\ref{fig:cmpbveto_POWHER-DR-DS}, where it is shown that the
difference between DR and DS is reduced also for differential
distributions, when the $b$-veto is applied. However, the plots in
fig.~\ref{fig:cmpbveto_POWHER-DR-DS} show also that observables are
potentially affected by the $b$-veto non uniformly (i.e.~ratios
between DR and DS results can be bin-dependent). This suggests that
particular care should be taken when one performs a full analysis
where the contribution from the $Wt$-channel process is supposed to be
relevant, since the size of interference effects may depend on the
cuts applied. As it was already observed by the \MCatNLO{} authors, a
comparison between DR and DS predictions gives an estimate of the
theoretical uncertainty due to these effects.

We have also performed some comparisons between DR and DS results
obtained with the \PYTHIA{} shower.  In order to maximize consistency
with the \POWHEG{} prescriptions, we used the $\pt$-ordered shower. In
fig.~\ref{fig:cmp1cuts1_POWPYT-DR-DS} we show the pseudorapidity and
the transverse momentum of the system made by the two hardest leptons,
while in fig.~\ref{fig:cmp2cuts1_POWPYT-DR-DS} the transverse momentum
of the hardest non $b$-flavoured jet ($\pt^{\,j_1}$) and for the
hardest jet that contains a $b$-flavoured hadron ($\pt^{\,bj_1}$) are
shown.  The plots on the left have been obtained without imposing any
cut. Instead, the plots on the right have been obtained using cuts
similar to the ``$Wt$-signal cuts'' of ref.~\cite{White:2009yt}. Jets
have been defined according to the $\kt$
algorithm~\cite{Catani:1993hr}, as implemented in the {\tt FASTJET}
package~\cite{Cacciari:2005hq}, setting $R=0.7$ and imposing a lower
$10$~GeV cut on jet transverse momenta.  To accept an event, we
required the following properties:
\begin{itemize}
\item There is exactly one $b$-jet with $\pt^j>50$ GeV and
  $|\eta^j|<2.5$. A $b$-jet is defined as a jet that contains at least
  one $b$-flavoured hadron and has $\pt^j>25$ GeV and $|\eta^j|<2.5$.
\item There are at least two light-flavoured jets with $\pt^j>25$ GeV
  and $|\eta^j|<2.5$. The invariant mass of the system made by the two
  hardest jets among these light-flavoured jets has to lie within $55$
  and $85$ GeV.
\item There is one lepton with $\pt^{\ell}>25$ GeV and
  $|\eta^{\ell}|<2.5$.  This lepton has also to be isolated with
  respect to the $b$-jet and the two light-flavoured jets, i.e.~its
  distance from the jets in the $(\eta,\phi)$ plane has to be larger
  than $0.4$.
\item The missing transverse energy is larger than $25$ GeV.
\end{itemize}
Although we have not performed as detailed a study as the one of
ref.~\cite{White:2009yt}, from figs.~\ref{fig:cmp1cuts1_POWPYT-DR-DS}
and~\ref{fig:cmp2cuts1_POWPYT-DR-DS} we observe that the DR and the DS
predictions are consistent (within the statistical accuracy) also when
the above cuts are applied, as was observed in the aforementioned
work.

\begin{figure}[!htb]
  \begin{center}
    \epsfig{file=cmp_etal1l2_wt_LHC10_POWPYT-DR-DS_cuts1,width=\figwidth}\\
    \epsfig{file=cmp_ptl1l2_wt_LHC10_POWPYT-DR-DS_cuts1,width=\figwidth}
  \end{center}
  \captskip
  \caption{\label{fig:cmp1cuts1_POWPYT-DR-DS} Comparisons between
    \POWHEG{} (interfaced to \PYTHIA{}) results obtained with the DR
    and DS prescriptions, for leptonic quantities.  Plots in the right
    panel have been obtained with the ``$Wt$-cuts'' described in the
    text. Uncertainties indicated by the vertical bars are only
    statistical. }
\end{figure}
\begin{figure}[htb]
  \begin{center}
    \epsfig{file=cmp_ptj1_wt_LHC10_POWPYT-DR-DS_cuts1,width=\figwidth}\\
    \epsfig{file=cmp_ptbj1_wt_LHC10_POWPYT-DR-DS_cuts1,width=\figwidth}
  \end{center}
  \captskip
  \caption{\label{fig:cmp2cuts1_POWPYT-DR-DS} Comparisons between
    \POWHEG{} (interfaced to \PYTHIA{}) results obtained with the DR
    and DS prescriptions.  Plots in the right panel have been obtained
    with the ``$Wt$-cuts'' described in the text. Uncertainties
    indicated by the vertical bars are only statistical.}
\end{figure}
\cleardoublepage
\section{Conclusions}
\label{sec:conclusions}
In this paper we have described the implementation of $Wt$-channel
single-top production at next-to-leading order in QCD, in the
\POWHEG{} framework. We have used the \BOX{} package, which is a
program that automates the algorithm first proposed in
ref.~\cite{Nason:2004rx} and then carefully described in
ref.~\cite{Frixione:2007vw}.  Having used this package, we described
how we calculated the main needed ingredients.

Since NLO corrections to single-top production in the $Wt$-channel are
known not to be well-defined (real contributions interfere with
$t\bar{t}$ production), we decided to follow the same strategy
originally proposed by the \MCatNLO{} authors in
ref.~\cite{Frixione:2008yi}: we included two definitions of the NLO
corrections, known as DR (\emph{Diagram Removal}) and DS
(\emph{Diagram Subtraction}), both of which can be used when the
merging of the fixed order result with parton showers is performed.
Moreover, the difference of results obtained with these two
prescriptions gives an estimate of the size of interference
effects. We have described how we included the two prescriptions in
\POWHEG{}, and how we dealt with DS and the problem of its exact
implementation within the \POWHEG{} method.

To check the correctness of the whole implementation and to assess the
validity of the choices we made, results have been compared with the
\MCatNLO{} program for the LHC, where $Wt$ production is relevant: we
found very good agreement, both for DR and DS.
We have also compared DR and DS results when a $b$-veto is imposed. We
found that a $b$-veto reduces interference effects with $t\bar{t}$,
the difference between DR and DS results becoming smaller when the
veto is in place.  Moreover, we have also presented some results
obtained with the \PYTHIA{} shower and with typical defining cuts for
the $Wt$ signal. Although we have not performed as full an analysis as
the one reported in ref.~\cite{White:2009yt}, good agreement between
DR and DS has been found also when typical $Wt$-channel cuts are
applied.

The main purpose of this work was the completion of the work presented
in ref.~\cite{Alioli:2009je}, where the \POWHEG{} implementation of
single-top $s$- and $t$-channel was described. Therefore, at present
all single-top processes can be simulated in the context of a NLO+PS
approach with \POWHEG{} as well as with \MCatNLO{}. We also stress
that, in \POWHEG{}, $Wt$-channel was the missing process among the
ones relevant for $H\to WW$ searches: Higgs production via
gluon~\cite{Alioli:2008tz} and vector-boson fusion~\cite{Nason:2009ai}
are already available, as well as $t\bar{t}$~\cite{Frixione:2007nw}
and $VV$~\cite{Hamilton:VV} production, which (together with
single-top $Wt$) are the main backgrounds .

The computer code for this \POWHEG{} implementation will soon be
available within the public branch of the \BOX{} package, that
can be downloaded at the site
\begin{center}
  \url{http://virgilio.mib.infn.it/~nason/POWHEG}
\end{center}

\section*{Acknowledgements}
I would like to particularly thank C. White, for providing me with the
one-loop corrections for this process, for many useful discussions on
this project, and also for comments on the manuscript.  I would also
like to thank R. Frederix and G. Zanderighi for helpful hints on the
usage of MadGraph and QCDloop.  I am also grateful to S. Alioli and
K. Hamilton for comments on the manuscript.