A DISK-BASED DYNAMICAL MASS ESTIMATE FOR THE YOUNG BINARY V4046 Sgr

Mass is the fundamental property that sets the evolutionary path of a star. The masses of young stars are of particular interest in many astrophysical problems: they provide unique information about the star formation process (e.g., accretion histories, the initial mass function; see Bastian et al. 2010) and are thought to have a substantial influence on the evolution of their circumstellar material, and therefore the efficiency of the planet formation process (e.g., see Alibert et al. 2011). Unfortunately, measurements of pre-main-sequence (pre-MS) star masses are difficult and accordingly rare. Unlike their more evolved counterparts, the location of a young star in the Hertzsprung–Russell (H-R) diagram does not provide a robust estimate of M_* (Hillenbrand & White 2004). The theoretical models for stellar evolution at these early stages are plagued with uncertainties related to rotation (Siess & Livio 1997; Mendes et al. 1999), accretion (Siess et al. 1997; Baraffe et al. 2009; Baraffe & Chabrier 2010), magnetic fields (D'Antona et al. 2000), atmosphere properties (e.g., convection, opacities; Baraffe et al. 2002), and unknown initial conditions. Ultimately, a more nuanced understanding of star and planet formation requires pre-MS evolution models that are empirically calibrated with direct, independent, and accurate M_* measurements.

The only direct methods available for measuring stellar masses are based on orbital dynamics. In sufficiently close pre-MS binary star systems, M_* can be estimated from the stellar orbits using multi-epoch radial velocity (RV) measurements (e.g., Mathieu et al. 1989, 1991, 1997) and/or long-term astrometric monitoring (Tamazian et al. 2002; Schaefer et al. 2003, 2006; Duchêne et al. 2006). For a double-lined spectroscopic binary (SB2), the RV method provides a robust estimate of M_*(sin i_*)³ for each component, where i_* is the inclination angle of the orbit projected on the sky. For "visual" binaries, the astrometric monitoring technique offers a constraint on the quantity M_totd³, where M_tot is the sum of the stellar masses and d is the distance to the binary. Generally, the stellar mass estimates from these methods are inherently uncertain due to their strong dependences on the unknown values of i_* or d. The {M_*, i_*} degeneracy is broken for the special case of an eclipsing SB2, but few pre-MS systems with such favorable orientations are known (e.g., Stassun et al. 2004; Morales-Calderón et al. 2012, and references therein). In a subset of ideal cases, the RV and astrometric techniques can be combined to alleviate the uncertainties related to {i_*, d} and extract accurate M_* values (e.g., Steffen et al. 2001; Schaefer et al. 2008; Boden et al. 2005, 2007, 2009, 2012).

These standard methods are only applicable for binary stars with a narrow range of orbital separations. Alternatively, for any isolated young star with a circumstellar disk, M_* can be determined from a single millimeter-wave interferometric observation of an optically thick emission line (with a linear dependence on d). This latter technique relies on modeling the spatially and spectrally resolved Keplerian rotation curve of the molecular gas disk that orbits the young star (Koerner et al. 1993; Dutrey et al. 1994, 1998; Simon et al. 2000). While this method has extraordinary value in its more general applicability, it has only been successfully employed for small samples. Attempts to expand its reach have been frustrated by molecular cloud contamination and observational limitations in resolution and sensitivity. Moreover, a reconstruction of the disk velocity field necessarily involves fitting such data with a relatively complicated model of the disk structure (Beckwith & Sargent 1993). Given that added complexity, there is naturally some concern about the absolute accuracy of the M_* estimates from this method (e.g., Gennaro et al. 2012), despite the impressive formal precision of the measurements (∼2%–3%; e.g., Piétu et al. 2007).

The young binary V4046 Sgr provides a rare opportunity to benchmark the disk kinematics method for estimating M_* against the more traditional RV technique for an SB2 system. V4046 Sgr is a nearly equal mass (q ≈ 0.94) pair of solar-type pre-MS stars in a circular (e ⩽ 0.01), non-eclipsing orbit with a 2.4 day period (asin i_* ≈ 5.1 R_☉; Byrne 1986; Quast et al. 2000; Stempels & Gahm 2004). The system is completely isolated from any known molecular clouds and has been kinematically associated with the ∼8–20 Myr old β Pic moving group; a moving cluster analysis suggests it is relatively nearby, d = 73 pc (Torres et al. 2006). Despite its advanced age, V4046 Sgr hosts a large and massive circumbinary disk that exhibits a rich molecular emission line spectrum (Kastner et al. 2008; Rodriguez et al. 2010; Öberg et al. 2011). Since the binary orbit is tight and circular, it has no dynamical impact on the disk structure outside a radius of ∼0.1 AU (e.g., see Artymowicz & Lubow 1994). With its central SB2 host, rare proximity to the Sun, lack of molecular cloud contamination, and intrinsically bright, spatially extended line emission, the V4046 Sgr disk is an ideal target to assess the accuracy of the disk kinematics technique for measuring M_*.

In this article, we build on some initial work by Rodriguez et al. (2010) and use high-quality spatially and spectrally resolved observations of the CO J = 2–1 emission line to measure the velocity field of the V4046 Sgr disk and extract the total mass of the close binary at its center. Our millimeter-wave observations with the Submillimeter Array (SMA) and data calibration procedures are described in Section 2. A detailed overview of the modeling analysis is provided in Section 3. The modeling results are presented and compared with the complementary RV analysis of the central SB2 by H. C. Stempels (2012, in preparation) in Section 4. These results are discussed in the context of the V4046 Sgr system in particular, pre-MS evolution models more generally, and future prospects for M_* estimates from the disk kinematics method in Section 5. Some key conclusions from this work are summarized in Section 6.

The V4046 Sgr circumbinary disk was observed at 225 GHz (1.3 mm) with SMA (Ho et al. 2004) on four occasions starting in 2009, using each of the available antenna configurations: sub-compact (baseline lengths of 9–25 m; 2011 March 18), compact (16–70 m; 2009 April 25), extended (28–226 m; 2009 February 23), and very extended (68–509 m; 2011 September 4). The SMA double sideband receivers were tuned to simultaneously observe the J = 2–1 transitions of ¹²CO, ¹³CO, and C¹⁸O at 230.538, 220.399, and 219.560 GHz, respectively, and the adjacent dust continuum. The correlator was configured to place those emission lines in separate 104 MHz spectral chunks and sample them finely with 512 channels per chunk, corresponding to a native velocity resolution of ∼0.25 km s⁻¹. The continuum was observed with a more coarse frequency sampling (in 3.25 MHz channels), with a total bandwidth of 1.6 and 3.6 GHz in 2009 and 2011, respectively. The observations cycled between V4046 Sgr and the quasars J1924-292 (15° away) and J1733-130 (22° from V4046 Sgr, 30° from J1924-292) with a total loop time of 10–20 minutes. The bright quasars 3C 84 and 3C 454.3 were also observed as bandpass calibrators, along with Uranus, Callisto, and Ceres for use in determining the absolute scaling of the amplitudes. All of the data were collected in outstanding weather conditions for this observing frequency, with an atmospheric zenith optical depth of only ∼0.05 (corresponding to a precipitable water vapor level of ∼1.0 mm).

Each individual data set was calibrated independently with the MIR software package. The bandpass response was corrected based on observations of bright quasars, and broadband continuum channels were generated from the central portions of all line-free spectral chunks. The visibility amplitude scale was derived by bootstrapping the gain calibrator (quasar) flux densities from the observations of Uranus, Callisto, or Ceres, with a systematic uncertainty estimated at 10%–15%. The antenna-based complex gain response of the system was determined with reference to J1924-292, and the quality of the phase transfer was assessed using the observations of J1733-130. That comparison suggests only a small amount of "seeing" (∼01) was introduced by atmospheric phase noise (or small baseline errors), consistent with the excellent observing conditions. After applying the appropriate (and small) phase shifts to account for the V4046 Sgr proper motion (μ_αcos δ = 0003 yr⁻¹, μ_δ = −0052 yr⁻¹; Zacharias et al. 2010) and confirming that the continuum amplitudes from different array configurations were consistent on overlapping baseline lengths, the visibility data sets from each observation were combined. The observations of the dust continuum and CO isotopologue emission will be presented in a separate article; the focus here will be solely on the ¹²CO J = 2–1 emission. Note that although the 2009 data were originally presented by Rodriguez et al. (2010), those data have been re-calibrated here for consistency (and modest improvements).

The CO visibilities were continuum subtracted and truncated outside a projected baseline length of 200 kλ to reduce the data volume and improve the signal-to-noise ratio. They were then Fourier inverted, deconvolved with the CLEAN algorithm, and restored with a synthesized beam using the MIRIAD package. The naturally weighted spectral images shown as channel maps in Figure 1 were synthesized on a 0.4 km s⁻¹ smoothed velocity scale with a 155 × 129 beam (at a position angle of 24°). The typical rms noise level in each channel is 70 mJy beam⁻¹. There is CO emission firmly detected (>3σ) out to ±4.8 km s⁻¹ from the systemic velocity, estimated to be v(LSR) = +2.87 ± 0.05 km s⁻¹ (corresponding to −6.26 ± 0.05 km s⁻¹ in the heliocentric frame), with an integrated intensity of 27.7 ± 2.8 Jy km s⁻¹ and a peak flux density of 1.48 ± 0.16 Jy beam⁻¹ (17 ± 2 K; a peak S/N = 21), including the calibration uncertainties. Figure 2 shows a map of the velocity-integrated CO intensities (zeroth moment; contours) overlaid on the intensity-weighted velocities (first moment; color scale), as well as a spatially integrated CO spectrum. These data exhibit a molecular gas disk with a clear rotation pattern, from east (blueshifted) to west (redshifted), and suggest a modest inclination angle to the line of sight. Near the systemic velocity, the CO emission subtends ∼5'' in radius (∼365 AU; see Rodriguez et al. 2010).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The fundamental goal here is to derive a dynamical estimate of the central stellar mass based on the kinematic properties of the V4046 Sgr gas disk. In order to extract M_* from a measurement of the disk velocity field, we need to construct a detailed physical model of the disk structure.³ We adopt a modeling formalism motivated by Beckwith & Sargent (1993), which makes three basic assumptions about disk properties. First, the disk material is assumed to be orbiting in Keplerian rotation around a point mass, meaning the central stars are treated as one object that dominates the disk velocity field. Previous work has found little evidence to support (simple parametric) deviations from Keplerian v-fields in disks (e.g., see Simon et al. 2000); the corollary criterion that the disk mass is only a small fraction of the stellar mass (M_d/M_* ≪ 1) can be confirmed a posteriori (see Section 4). Second, the disk is assumed to be geometrically thin at all radii. Although Piétu et al. (2007) argued that this may not be valid at very large radii (>800 AU), the higher M_* for the V4046 Sgr binary (which decreases the disk scale height; see below) and the limited extent of the CO emission from its disk (<400 AU) suggest it is a reasonable approximation in this case. And third, in the context of the disk structure, we assume that the gas is vertically isothermal and in hydrostatic pressure equilibrium. A more realistic model would include a temperature inversion (e.g., D'Alessio et al. 1998). However, excluding that kind of added complexity does not diminish the accuracy of an M_* estimate: the emission from a single CO line is generated in a narrow vertical layer of the disk atmosphere, which has a roughly constant temperature (Dartois et al. 2003).

The model for the gas density distribution is constructed in cylindrical coordinates (r, ϕ, z),

$$\begin{equation} \rho (r,z) = \frac{\Sigma }{\sqrt{2\pi } H_p} \exp \left[ -\frac{1}{2} \left(\frac{z}{H_p} \right)^2 \right], \end{equation} \tag{ 1 }$$

where Σ is the radial surface density profile and H_p is the pressure scale height at each radius. We assume a parametric version of the former that is appropriate for an accretion disk with a static, power-law viscosity profile (Lynden-Bell & Pringle 1974; Hartmann et al. 1998) and is currently the basis for most dust-based disk density measurements (e.g., Andrews et al. 2009, 2010),

$$\begin{equation} \Sigma = \Sigma _c \left(\frac{r}{r_c}\right)^{-\gamma } \exp \left[-\left(\frac{r}{r_c}\right)^{2-\gamma }\right], \end{equation} \tag{ 2 }$$

where γ is a density gradient, r_c is a characteristic radius, and Σ_c is a normalization equivalent to e · Σ(r_c). The scale heights are calculated with the explicit assumption of hydrostatic equilibrium and a constant vertical temperature profile, such that

$$\begin{equation} H_p = \frac{c_s}{\Omega } = \left(\frac{kT}{\mu m_{\rm H}} \cdot \frac{r^3}{G M_{\ast }} \right)^{1/2}, \end{equation} \tag{ 3 }$$

where c_s is the sound speed, Ω is the angular velocity, T is the radial temperature profile, k is the Boltzmann constant, G is the gravitational constant, m_H is the mass of a hydrogen atom, and μ = 2.37 is the mean molecular weight of the gas. We further assume a simple power-law behavior for the radial temperature profile,

$$\begin{equation} T = T_{10} \left(\frac{r}{{\rm 10\,AU}}\right)^{-q}, \end{equation} \tag{ 4 }$$

where q is a temperature gradient and T₁₀ is the gas temperature at a radius of 10 AU. The parametric descriptions in Equations (1)–(4) completely characterize the temperature and density structure of a model gas disk. The gas kinematics are described by a Keplerian velocity field,

$$\begin{equation} v_{\phi }(r) = v_k = \sqrt{\frac{G M_{\ast }}{r}} \,; \,\,\,\, v_r = v_z = 0, \end{equation} \tag{ 5 }$$

meaning there is only rotation, and no net motion in the radial or vertical dimensions.

Given a physical disk structure, we can construct a model of the CO J = 2–1 emission with the added assumption that the energy levels relevant for this transition are populated according to local thermodynamic equilibrium (LTE). Although the LTE approximation is not always valid in protoplanetary disks, Pavlyuchenkov et al. (2007) demonstrated that it is an appropriate simplification for the low energy and high optical depths associated with this particular species and transition. We model the emission line intensity distribution by integrating the radiative transfer equation along each sight line s, so that

$$\begin{equation} I_{\nu } = \int _0^{\infty } K_{\nu }(s) S_{\nu }(s) e^{-\tau _{\nu }(s)} ds, \end{equation} \tag{ 6 }$$

where K_ν(s) is the absorption coefficient, S_ν(s) = B_ν(T) is the source function (here equivalent to the Planck function), and τ_ν(s) = ∫^s₀K_ν(x)dx is the optical depth along the line of sight. The absorption coefficient is the sum of contributions from dust and gas. For the former, we assume K_ν(s)_dust = ρ(s) ζ κ_ν, where the dust-to-gas mass ratio is ζ = 0.01 and the grain opacity is κ_ν ≈ (ν/100 GHz) cm² g⁻¹ (Beckwith et al. 1990). In the particular case of the V4046 Sgr disk, the contribution of dust to the absorption coefficient is effectively negligible compared to the gas; the specific choices of ζ and κ_ν (within reasonable limitations) have no tangible effect on our results (moreover, continuum emission is removed from the data before our modeling analysis).

The absorption coefficient of the CO gas is calculated from the transition cross section weighted by the population of the lower energy level, ℓ, such that K_ν(s)_co = n_ℓ(s)σ_ν(s) (note that in this specific case, ℓ = 1). Since we assume that the line is thermally populated in LTE, the level populations are determined by the local disk temperature via the Boltzmann equation,

$$\begin{equation} n_{\ell }(s) = \frac{X_{\rm co} \, \rho (s)}{\mu {\rm m}_{\rm H}} \cdot \frac{g_{\ell }}{Z} \, \exp \left(- \frac{E_{\ell }}{kT(s)} \right), \end{equation} \tag{ 7 }$$

where E_ℓ is the transition energy, g_ℓ = 2ℓ + 1 is the statistical weight, Z is the partition function, and X_co is the CO fractional abundance (assumed to be constant everywhere in the disk). The details of the emission line model are encoded in the absorption cross section,

$$\begin{equation} \sigma _{\nu }(s) = \phi _{\nu }(s) \cdot \sigma _0 (1 - e^{-h \nu / k T(s)}), \end{equation} \tag{ 8 }$$

where ϕ_ν is the line profile function and the integrated cross section is

$$\begin{equation} \sigma _0 = \frac{h\nu }{4\pi } \cdot \frac{g_{\ell +1}}{g_{\ell }} B_{21} = \frac{c^2}{8 \pi \nu ^2} \cdot \frac{g_{\ell +1}}{g_{\ell }} A_{21}, \end{equation} \tag{ 9 }$$

where we have used the Einstein relation in the last equality to express the Einstein-B coefficient in terms of the Einstein-A coefficient provided by the LAMDA molecular database (Schöier et al. 2005). The line profile function naturally determines the shape of the emission line, where a given frequency ν corresponds to an intrinsic Doppler velocity v_d = (c/ν₀)(ν − ν₀) relative to the line center, ν₀. The gas in a given disk model has a projected, line-of-sight velocity field given by v_obs ≈ v_ϕ(r)cos ϕsin i_d, where i_d is the inclination of the disk relative to the observer (such that i_d = 90° is edge-on; for the geometry, see Isella et al. 2007). The line profile shape at each frequency is determined by the difference between the Doppler and line-of-sight velocities,

$$\begin{equation} \phi _{\nu }(s) = \frac{c}{\sqrt{\pi } \nu _0 \Delta v} \exp \left[ - \left(\frac{v_d-v_{\rm obs}}{\Delta v}\right)^2 \right], \end{equation} \tag{ 10 }$$

and the effective line width, Δv. The latter is comprised of the quadrature sum of contributions from thermal and non-thermal broadening terms,

$$\begin{equation} \Delta v = \left(\frac{2 k T(s)}{m_{\rm co}} + \xi ^2 \right)^{1/2}, \end{equation} \tag{ 11 }$$

where m_co is the mean molecular weight of CO and ξ is the contribution from microturbulence. We assume that the turbulent velocity width ξ is constant throughout the disk.

The Beckwith & Sargent (1993) model formalism highlighted above is admittedly complex. A single model is completely described by 14 parameters: four quantify the distribution of CO densities {X_co, Σ_c, R_c, γ}, two characterize the temperature structure (and therefore the vertical density structure) {T₁₀, q}, three play pivotal roles describing the disk kinematics {M_*, ξ, ν₀}, and the remaining five relate the model to the observations—i.e., convert from physical coordinates in the disk to the observed plane projected on the sky—including the disk inclination and major axis position angle {i_d, PA_d}, distance {d}, and disk center {α₀, δ₀}. Some of these parameters can be reliably determined in a simple way, and then fixed to facilitate the data modeling problem (with no tangible impact on the quality of the M_* determination). In this case, the systemic LSR velocity (effectively ν₀) was measured directly from the SMA data and set to v₀ = +2.87 km s⁻¹. The disk center coordinates were estimated from both the dust continuum and CO emission morphology at the systemic velocity (see Figure 1) and set to $\alpha _0 = 18^{\rm h}14^{\rm m}10\mbox{$.\!\!^{\mathrm s}$}48$, δ₀ = −32°47'3508 (J2000), coincident with the composite V4046 Sgr stellar position in the UCAC3 catalog (Zacharias et al. 2010). We adopted the distance inferred by Torres et al. (2006) from the kinematic parallax (moving cluster) technique, d = 73 pc (but see Section 5 regarding alternative values).

Some additional practical simplifications can be made, since the focus here is on the velocity field and not the disk structure details. Because the line opacity is so much greater than the dust opacity, the normalization parameters X_co and Σ_c are not independent: we can only constrain their product. In practice, we adopt a joint CO disk mass parameter, M_co = X_coΣ_c(2πr²_c)/(2 − γ), for that purpose. Moreover, after extensive experimentation, we concluded that two more parameters should be fixed. First, we found that the disk orientation was very well determined (within 1°) at PA_d = 76° (E of N), so that continued iteration on a more precise value was a waste of computational resources: fixing this parameter has negligible quantitative impact on the other model parameters or their uncertainties. And second, we fixed the surface density gradient, γ = 1. The key parameters that set the density profile, {γ, r_c}, are anti-correlated and strongly degenerate; statistically indistinguishable model fits can be found over a large range of these parameters. The degeneracy is remarkably narrow and could easily be missed with standard minimization algorithms (resulting in local χ² minima with misleadingly tight constraints on the gradient). Fortunately, this degeneracy has minimal impact on the parameters most relevant for characterizing the disk velocity field: the precision and accuracy of an M_* determination are not notably affected by the uncertainties in {γ, r_c}. With that in mind, we have made the modeling process more tractable with a fixed γ.

Making use of those simplifications, a synthetic CO spectral datacube can be calculated by specifying seven free parameters, {M_co, r_c, T₁₀, q, ξ, i_d, and M_*}. That model datacube is then resampled at the observed velocities and spatial frequencies and processed in the same way as the data (see Section 2) to produce a set of synthetic spectral visibilities. Those model visibilities are evaluated with respect to the data by computing a composite χ² value, summed over the real and imaginary components in 25 spectral channels. The best-fit set of model parameters and their associated uncertainties were determined with a Monte Carlo Markov Chain (MCMC) technique, using the affine invariant ensemble sampler developed by Goodman & Weare (2010; see also Foreman-Mackey et al. 2012), using a jump probability $\propto e^{-\Delta \chi ^2/2}$. To our knowledge, this is the first application of these more sophisticated MCMC methods for parameter estimation in this particular context: previous studies have relied on downhill simplex routines (e.g., Guilloteau & Dutrey 1998; Simon et al. 2000), sometimes with clever modifications to address asymmetric uncertainties (Piétu et al. 2007). Although comparatively the MCMC technique is computationally expensive, this Bayesian treatment has the distinct advantage of providing the posterior probability distributions for each parameter in the complex, multi-dimensional parameter-space of the underlying disk model.

The best-fit model derived from the SMA data is compiled in Table 1, which lists the mode (peak) of the posterior probability distribution for each parameter along with its uncertainty, based on the range of values that encompass 68% of the distribution area (i.e., 1σ uncertainties for a Gaussian distribution). Figure 3 makes a direct comparison between the data and the best-fit model in the spectral image plane, demonstrating the fit quality with only low-level (statistically insignificant) residuals present in the channel maps. The best-fit model has a reduced $\tilde{\chi }^2 = 967,427/977,000 = 0.9902$. A more detailed view of the multi-dimensional results of the modeling analysis is presented in Figure 4, where we have taken the 30,000 MCMC samples computed and marginalized over parameter subsets to display one- and two-parameter posterior probability distributions. The marginalized distributions for individual parameters are single peaked, while the paired distributions highlight internal degeneracies in the model. For example, there are tight correlations between normalizations and sizes (e.g., {M_co, r_c}, to conserve the integrated line intensity; this would be modified were γ a free parameter) or gradients (e.g., {T₁₀, q}, to maintain an appropriate temperature in the outer disk; Mundy et al. 1996; Andrews & Williams 2007), as well as more subtle associations among parameters related to line broadening (among r_c, T₁₀, and ξ; Guilloteau & Dutrey 1998; Hughes et al. 2011). The degeneracy of most interest here is the {M_*, i_d} anti-correlation, which is a natural consequence of reproducing the observed line-of-sight velocity pattern, $v_{\rm obs} \propto v_k \sin {i_d} \propto \sqrt{M_{\ast }} \sin {i_d}$ (i.e., M_*∝1/sin ²i_d).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The best-fit model parameters that describe the V4046 Sgr disk structure are typical for young protoplanetary disks. The adopted (fixed) surface density gradient (γ = 1) and characteristic radius (r_c = 45⁺⁵_{− 3} AU) lie near the median values for a survey of disk structures in the Ophiuchus region (Andrews et al. 2009, 2010); the temperature profile (T₁₀ = 115 ± 5 K, q = 0.63 ± 0.01; and therefore vertical density distribution) and turbulent linewidth (ξ = 0.14 ± 0.01 km s⁻¹) are comparable to those inferred from CO emission in other T Tauri disks (e.g., Guilloteau & Dutrey 1998; Piétu et al. 2007; Hughes et al. 2008, 2011). If we assume that the CO mass fraction found in dark clouds (X_co ≈ 10⁻⁴) is also applicable in the disk, the molecular density normalization parameter M_co implies a modest gas mass, 2.8^+3.7_{− 1.5} × 10⁻² M_☉, validating a posteriori our assumption that self-gravity is negligible (M_d/M_* ≈ 10⁻² ≪ 1). In any case, our simple treatment of the disk structure is really only used as a means to an end: the focus here is to place a firm constraint on the key parameter that determines the behavior of the disk velocity field—the central stellar mass, M_*.

We infer a total stellar mass of M_* = 1.75^+0.09_{− 0.06} M_☉ for the V4046 Sgr spectroscopic binary, a formally precise constraint with a relative uncertainty of 3%–5%. Figure 4 demonstrates that this estimate of M_* is essentially independent of the disk structure parameters in the model. The sole relevant degeneracy is with the disk inclination, which we infer to be i_d = 335^+0.7_{− 1.4} (a relative precision of ∼2%–5%). Within the quoted uncertainties, these {M_*, i_d} values are consistent with the previous determinations by Rodriguez et al. (2010), who used a simple χ² grid search and an altogether different underlying disk structure model. These {M_*, i_d} measurements based on the circumstellar gas disk kinematics can be directly compared with the joint (degenerate) constraints on {M_*, i_*} imposed by RV measurements of the spectroscopic binary itself. Although various groups have studied this SB2 system (Quast et al. 2000; Stempels & Gahm 2004; Donati et al. 2011), by far the most robust constraints on the stellar parameters come from the long-term, extensive, and high resolution optical spectroscopic monitoring campaign conducted by H. C. Stempels (2012, in preparation). In that work, the RV data are found to be best explained with a binary that has a total mass M_*(sin i_*)³ = 0.2923 ± 0.0007 M_☉. The inset in Figure 4 confirms that our disk kinematics constraints and the RV constraints by H. C. Stempels (2012, in preparation) coincide in the stellar-mass–inclination plane: the two completely independent methods find the same results well within their formal 1σ uncertainties. Adopting our best estimate of M_* and propagating the relevant uncertainties, the RV constraints suggest that i_* = 3342 ± 0.01. Therefore, we quantitatively confirm the finding of Rodriguez et al. (2010): the V4046 Sgr spectroscopic binary and its associated, large-scale circumbinary disk are co-planar within |i_d − i_*| ≈ 01 ± 1° across ∼4 orders of magnitude in radial scale.

We have presented spatially and spectrally resolved SMA observations of the ¹²CO J = 2–1 line emission from the isolated, nearby, and gas-rich circumstellar disk that orbits the close, pre-MS SB2 V4046 Sgr. Adopting a simple parametric model for the structure of a flared Keplerian disk, we employ an MCMC technique to infer the properties of the disk velocity field from these data and thereby provide a robust statistical estimate for the total mass of the central binary. We find that these CO line data are best reproduced for a binary mass M_* = 1.75^+0.09_{− 0.06} M_☉ and a disk viewing angle inclined by i_d = 335^+0.7_{− 1.4} from face-on. Those values are in excellent agreement with the completely independent inferences of H. C. Stempels (2012, in preparation), made from his extensive RV monitoring campaign of the spectroscopic binary itself. The orbital planes of the binary (a = 0.045 AU) and its associated circumbinary disk (on radial scales out to ∼400 AU) are aligned within ∼01 ± 1°. These results demonstrate that, despite its complexity, the disk-based kinematic method for estimating the masses of young stars is both precise (at the level of a few percent) and accurate in an absolute sense, as verified here by an entirely independent dynamical constraint.

The disk-based dynamical estimate of M_* can be combined with the RV constraints of H. C. Stempels (2012, in preparation) and other information from the literature to assess the predictions of pre-MS stellar evolution models. Coupling our M_* estimate with the stellar mass ratio (q) and mass function (M_*(sin i_*)³) determined by Stempels, we can derive the masses of the individual stellar components in the V4046 Sgr binary, 0.90 ± 0.05 and 0.85 ± 0.04 M_☉, and their orbital inclination, i_* = 3342 ± 0.01.⁴ Assuming that the stellar rotation axes are aligned with the binary orbital axis (a spin–orbit alignment presumably induced by tides in this close, circularized system; e.g., Zahn 1977; Hut 1981; Melo et al. 2001), the radial and rotational velocities measured for each stellar component by Stempels can be used to determine the stellar radii, 1.25 ± 0.04 and 1.21 ± 0.04 R_☉. Stempels & Gahm (2004) inferred spectral types of K5 and K7 in this system, corresponding to effective temperatures of 4350 ± 240 and 4060 ± 210 K for the standard conversion advocated by Schmidt-Kaler (1982); an ambiguity of one spectral subclass has been assumed to estimate the temperature uncertainties. Those measurements imply that the individual stellar luminosities are 0.50 ± 0.11 and 0.36 ± 0.08 L_☉ (a total luminosity of 0.86 ± 0.14 L_☉).

Figure 5 shows the composite optical/near-infrared spectrum expected from the V4046 Sgr binary, given the effective temperatures, stellar radii, and masses derived above from the combination of RV data and our disk-based measurements (black curve, with uncertainty in shaded gray). This spectral energy distribution prediction was generated from an interpolated grid of Lejeune et al. (1997) model spectra, assuming the appropriate stellar gravities (log g = 4.20 ± 0.02 in both cases), negligible interstellar reddening (A_V = 0), and our adopted d = 73 pc. Representative broadband photometric data are marked as the data points in Figure 5, constructed from the weighted mean magnitudes of Hutchinson et al. (1990), Strassmeier et al. (1993), Jensen & Mathieu (1997), and the Two Micron All Sky Survey (Skrutskie et al. 2006) and DENIS (Epchtein et al. 1997) surveys: uncertainties were taken as the quadrature sum of the standard deviation of the weighted mean and the maximal deviation from the weighted mean, in an effort to more faithfully represent potential uncertainties due to variability. A quick examination of Figure 5 demonstrates that the spectral morphology of the data and model prediction are a good match, although the normalizations appear discrepant. Scaling the (best-fit) predicted spectrum down by ∼30% ± 15% provides a good match to the observed spectrum; a composite luminosity of 0.60 ± 0.13 L_☉ is more appropriate (red curve). Therefore, the predicted and observed composite binary luminosities are marginally (∼1.4σ) different.

Zoom In Zoom Out Reset image size

Although formally the conflict in these luminosity values is not statistically significant, it is still interesting to consider some potential paths for a better reconciliation of all the data. Perhaps the most straightforward means of doing so is to modify the assumed distance to the system: taken at face value, a ∼30% increase to d ≈ 95 pc implies an inherently more luminous pair of stars with a spectrum that would be in good agreement with the observations. However, that same re-normalization impacts the total stellar mass inferred from the CO data, with a linear scaling to M_* ≈ 2.3 M_☉ that introduces a substantial (∼2σ) discrepancy in the {M_*, i}-plane between the disk-kinematics and RV techniques for estimating stellar parameters. Although it hardly seems worthwhile to trade one marginal discrepancy for another (which is technically less marginal), there is no a priori reason that the orbital planes of the disk and binary need to be aligned with the precision inferred in Section 4. Unfortunately, little guidance is provided in the form of an uncertainty on the kinematic parallax measurement of Torres et al. (2006). In any case, discrepancies in effective temperatures and luminosities are not necessarily uncommon for active young stars (Stassun et al. 2012), and particularly for those in close binary systems (e.g., Goméz Maqueo Chew et al. 2012).

Ultimately, dwelling on a marginal luminosity discrepancy is not well justified, especially given the independent distance estimate from the Torres et al. (2006) study. In the following, we adopt the stellar parameters inferred from the joint constraints of the CO disk spatio-kinematics and the optical RV studies of Stempels assuming d = 73 pc, but use individual stellar luminosities scaled down by 30% to best match the observed optical and near-infrared spectrum: 0.35 ± 0.10 L_☉ and 0.25 ± 0.08 L_☉ (note that each component is slightly less luminous than reported by Donati et al. 2011). With those results, we can make comparisons with the predictions of pre-MS evolution models to estimate the ages (t_*) and masses of the individual stars in the V4046 Sgr binary. To accomplish that we follow the Bayesian methodology of Jørgensen & Lindegren (2005), which employs a finely interpolated grid of pre-MS models in the H-R diagram to estimate the probability distributions of stellar mass and age given the measured values (and uncertainties) of effective temperatures and luminosities (see also Gennaro et al. 2012). If we define the observables as {x, y} = {log T, log L} with associated uncertainties {σ_x, σ_y} and the pre-MS model predictions {$\hat{x}(M_{\ast }$, t_*), $\hat{y}(M_{\ast }$, t_*)}, we can write the likelihood function as a multivariate Gaussian,

$$\begin{equation} \mathcal {L}(\hat{x},\hat{y}|x, y) = \frac{1}{2\pi \sigma _x \sigma _y } \exp \left\lbrace -\frac{1}{2} \left[ \frac{(x - \hat{x})^2}{\sigma _x^2} + \frac{(y- \hat{y})^2}{\sigma _y^2} \right]\right\rbrace. \vspace*{4pt} \end{equation} \tag{ 12 }$$

The best-fit {M_*, t_*} are then directly determined from the {$\hat{x}$, $\hat{y}$} that maximize the likelihood, with uncertainties that can be calculated from the shape of the likelihood distribution. This procedure was conducted for four different pre-MS evolution models, assuming solar composition, a fractional deuterium abundance of 2 × 10⁻⁵, and a convective mixing parameter α ≈ 1.7–2.0: D'Antona & Mazzitelli (1998), Baraffe et al. (1998), Siess et al. (2000), and Tognelli et al. (2011).

The results of these comparisons are presented in Table 2. Figure 6 shows the V4046 Sgr stellar properties that were inferred from each of the four reference sets of pre-MS evolutionary model tracks in the H-R diagram. All of the evolutionary models make predictions for the V4046 Sgr stellar masses that are remarkably consistent with each other and the dynamical masses inferred in Section 4. The stellar masses that formally maximize the likelihood in Equation (12) are systematically ∼2%–10% below the best-fit dynamical masses, an underprediction typical of pre-MS models regardless of the dynamical method used to estimate M_* (Hillenbrand & White 2004; Gennaro et al. 2012). However, those modest discrepancies are not statistically significant, given the uncertainties on {L, T} and the dynamical masses. An examination of Figure 6 demonstrates that the V4046 Sgr stars are nearly evolved off the Hayashi track, having presumably developed radiative cores as expected for their solar-like masses. The clear implication of their location in the H-R diagram is that the V4046 Sgr binary is comparatively old for a pre-MS system, as has been previously reported by Donati et al. (2011) and Kastner et al. (2011). The models used here suggest that a large range of ages are plausible, from ∼5 to 30 Myr, and confirm that the binary components are coeval. If we include a Gaussian prior representing the dynamical masses into the Bayesian analysis of the H-R diagrams described above (labeled as "+dyn" in Table 2), we can infer a smaller range of acceptable "dynamical" ages from each set of pre-MS models. The different pre-MS model predictions with these dynamical priors are shown together, and averaged (weighted by the posterior probability for each model; Hoeting et al. 1999), in Figure 7: the averaged results suggest ages of 12⁺¹⁷_{− 3} and 13⁺¹¹_{− 3} Myr for V4046 Sgr A and B, respectively. The corresponding coeval age is 13⁺⁸_{− 3} Myr, in good agreement with the age constraints from the putative far-flung companion(s) V4046 Sgr C[D] (at a separation of ∼12,350 AU or 28 on the sky; Kastner et al. 2011).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Torres et al. (2006) suggested that the observed space motion of V4046 Sgr is consistent with a high probability of membership in the β Pic moving group, a widespread and local (d ≈ 34 ± 21 pc) stellar association with an age range estimated to be ∼8–20 Myr (e.g., Zuckerman et al. 2001; Zuckerman & Song 2004; Torres et al. 2008). The kinematic parallax calculated by Torres et al. (d ≈ 73 pc) and component ages estimated here for the V4046 Sgr binary are certainly supportive of that conclusion. If that is the case, it is worth noting the uniqueness of the V4046 Sgr circumstellar environment: there are no other β Pic moving group members that are known to host such a massive, long-lived, gas-rich accretion disk. Regardless of its original association, the proximity and advanced age of V4046 Sgr make for a remarkable case study in both the long-term evolution of protoplanetary disk structure and the fundamental properties of pre-MS binary stars.

Although V4046 Sgr is a particularly striking example of the methodology behind disk-based dynamical estimates of stellar masses, we anticipate that these techniques will find considerably more use in the near future as the Atacama Large Millimeter Array (ALMA) project is completed. Even with vastly improved data quality, this simplified model to extract the gas disk rotation curve from such interferometric observations (see Section 3) remains complex. However, we have demonstrated that the method is robust, accurate, and precise, by using independent dynamical constraints from the V4046 Sgr spectroscopic binary RV monitoring results. In practice, these constraints and others like them (notably for UZ Tau E; see Simon et al. 2000; Prato et al. 2002) effectively serve as a check on the absolute calibration of the modeling technique described in Section 3. The consistency between the independent dynamical constraints on the V4046 Sgr stellar masses validates the application of this procedure to dynamically "weigh" single, isolated pre-MS stars based on their gas disk kinematics. Along with an investment in accurate stellar luminosity and temperature measurements, ALMA observations of the molecular gas kinematics in circumstellar disks will usher in a new era of precision in the fundamental astrophysical properties of young stars.

We have presented sensitive, high-resolution SMA observations of the ¹²CO J = 2–1 line emission from the massive, gas-rich disk orbiting the SB2 star V4046 Sgr. Using simple radiation transfer calculations for a disk structure model with a Keplerian velocity field, we fit the observed spectral visibilities using a stochastic model optimization technique that simultaneously infers model parameter values and their uncertainties. Our specific focus has been on the key model parameters that describe the disk velocity field, with a goal of placing a firm, dynamical constraint on the mass of the central binary. The key conclusions of our analysis are as follows.

• 1.  
  From modeling the CO line emitted by the circumbinary disk, we infer that the total stellar mass of the V4046 Sgr binary is M_* = 1.75^+0.09_{− 0.06} M_☉, assuming the kinematic parallax distance of 73 pc estimated by Torres et al. (2006). That measurement is in excellent agreement with the independent dynamical constraints imposed from the analysis of an optical spectroscopic monitoring campaign of the RV variations from the binary itself (H. C. Stempels 2012, in preparation).
• 2.  
  The mutual consistency of these distinct dynamical constraints on the stellar masses verifies that the mass determination based on the velocity field of the disk gas is accurate in an absolute sense, as well as remarkably precise (with a 3%–5% formal uncertainty on M_*).
• 3.  
  That same combination of constraints from millimeter and optical spectroscopic measurements confirms that the orbital planes of the stars and their accompanying circumbinary disk are coplanar, with inferred inclination angles (i_d = 335^+0.7_{− 1.4} and i_* = 334 ± 0.01) that differ by 01 ± 1° over roughly four decades in radius, from ∼0.04 to 400 AU.
• 4.  
  The inferred component masses of the V4046 Sgr binary are in uniformly good agreement with a variety of pre-MS evolution model predictions. When combined with our dynamical constraints, those same models confirm the coevality of the binary components and suggest an average age for the system of 13⁺⁸_{− 3} Myr. Therefore, V4046 Sgr hosts one of the oldest and nearest gas-rich primordial accretion disks currently known.

We are very grateful to Joel Kastner, Guillermo Torres, and Zachory Berta for some insightful discussions and suggestions, as well as to Germano Quast for being kind enough to provide additional information on their optical observations of the V4046 Sgr binary. The SMA is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica.