High Precision Dynamical Masses of Pre-Main Sequence Stars with ALMA and Gaia

TL;DR

This paper demonstrates how combining ALMA and Gaia data allows for highly precise dynamical mass measurements of pre-main sequence stars, providing critical tests for stellar evolutionary models.

Contribution

It introduces a method using ALMA and Gaia data to measure stellar masses with ~2% precision, and compares these with model predictions for young stars.

Findings

Masses of three young stars measured with ~2% accuracy.

Broad consistency found between dynamical and evolutionary model masses.

Potential discrepancies suggest areas for refinement in models.

Abstract

The Keplerian rotation in protoplanetary disks can be used to robustly measure stellar masses at very high precision if the source distance is known. We present Atacama Large Millimeter/submillimeter Array (ALMA) observations of spatially and spectrally resolved $^{12}$CO (2--1) emission towards the disks around 2MASS J16262774--2527247 (the tertiary companion to ROXs 12 at 5100 au), CT Cha, and DH Tau. We employ detailed modeling of the Keplerian rotation profile, coupled with accurate distances from $Gaia$, to directly measure the stellar masses with $\sim$2\% precision. We also compare these direct mass measurements with the masses inferred from evolutionary models, determined in a statistically rigorous way. We find that 2MASS J16262774--2527247 has a mass of $0.535^{+0.006}_{-0.007}~M_{\odot}$ and CT Cha has a mass of $0.796^{+0.015}_{-0.014}~M_{\odot}$, broadly consistent with evolutionary models, although potentially significant differences remain. DH Tau has a mass of $0.101^{+0.004}_{-0.003}~M_{\odot}$, but it suffers from strong foreground absorption that may affect our mass estimate. The combination of ALMA, $Gaia$, and codes like \texttt{pdspy}, presented here, can be used to infer the dynamical masses for large samples of young stars and substellar objects, and place constraints on evolutionary models.

Figures

Tables

Table 1. Table 1: Log of ALMA Observations

Source	Observation Date	Baselines	Integration Time	Beam	RMS	Calibrators
	(UT)	(m)	(min)		(mJy)	(Flux/Bandpass, Gain)
DH Tau	Sep. 14 2016	15 - 3247	13	0.30” $\times$ 0.17”	4	J0510+1800, J0433+2905
2M1626-2527	Sep. 16 2016	15 - 3143	12	0.21” $\times$ 0.18”	3.1	J1517-2422, J1634-2058
CT Cha	Sep. 27 2016	15 - 3247	13.5	0.27” $\times$ 0.14”	3.8	J1107-4449, J1058-8003

Table 2. Table 2: Best-fit Model Parameters

Source	$M_{\ast }$	$\log M_{disk}$	$R_{in}$	$R_{disk}$	$\gamma$	$T_0$	$q$	$a_{turb}$	$v_{sys}$	$i$	$P.A.$	$d_{pc}$
	$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$	[AU]	[AU]		[K]		[km s-1]	[km s-1]	${[}^{\circ }]$	${[}^{\circ }]$	[pc]
CT Cha	${0.796}_{-0.014}^{+0.015}$	$>-2.89$		${9.3}_{-5.1}^{+12.0}$	${0.9}_{-0.4}^{+0.2}$	${106}_{-14}^{+22}$	${0.44}_{-0.03}^{+0.04}$		${4.608}_{-0.014}^{+0.014}$	${54.0}_{-0.9}^{+1.0}$	${49.5}_{-0.6}^{+0.6}$	${191.8}_{-0.8}^{+0.7}$
DH Tau	${0.101}_{-0.003}^{+0.004}$	$-{4.54}_{-0.36}^{+0.32}$			${1.4}_{-0.2}^{+0.1}$	${237}_{-24}^{+33}$	${0.53}_{-0.04}^{+0.04}$		${5.652}_{-0.010}^{+0.009}$	${48.4}_{-1.5}^{+1.4}$	${2.5}_{-0.9}^{+0.7}$	${135.4}_{-1.2}^{+1.2}$
2M1626-2527	${0.535}_{-0.007}^{+0.006}$	$-{1.07}_{-0.87}^{+0.49}$			${1.3}_{-0.1}^{+0.1}$	${57}_{-4}^{+5}$	${0.20}_{-0.02}^{+0.02}$	${0.137}_{-0.010}^{+0.010}$	${4.150}_{-0.006}^{+0.005}$	${54.0}_{-0.4}^{+0.5}$	${95.5}_{-0.2}^{+0.2}$	${138.2}_{-1.3}^{+1.4}$

Table 3. Table 3: Stellar Photosphere Model Parameters & Stellar Parameters Inferred from Evolutionary Tracks

Parameters	2M1626-2527	CT Cha	DH Tau
$T_{eff}$ [K]	${3746}_{-122}^{+139}$	${4402}_{-166}^{+151}$	${3628}_{-65}^{+81}$
$\log g$	${3.6}_{-0.7}^{+0.6}$	${4.5}_{-1.2}^{+0.7}$	${3.1}_{-0.4}^{+0.4}$
$R_{\ast }$ [$R_{\odot }$]	${1.17}_{-0.04}^{+0.06}$	${2.06}_{-0.05}^{+0.05}$	${1.66}_{-0.03}^{+0.04}$
$d_{pc}$ [pc]	${138.3}_{-1.4}^{+1.4}$	${191.8}_{-0.8}^{+0.8}$	${135.4}_{-1.3}^{+1.2}$
$A_V$	${2.5}_{-0.6}^{+0.6}$	${1.6}_{-0.3}^{+0.3}$	${0.5}_{-0.3}^{+0.3}$
$L_{\ast }$ [$L_{\odot }$]	${0.24}_{-0.04}^{+0.05}$	${1.41}_{-0.16}^{+0.17}$	${0.42}_{-0.04}^{+0.04}$
$M_{\ast }$ [M⊙]
Baraffe et al. 2015	${0.53}_{-0.10}^{+0.12}$	${1.06}_{-0.20}^{+0.20}$	${0.40}_{-0.03}^{+0.05}$
Bressan et al. 2012	${0.43}_{-0.08}^{+0.10}$	${0.87}_{-0.17}^{+0.18}$	${0.34}_{-0.03}^{+0.04}$
Chen et al. 2014	${0.71}_{-0.04}^{+0.01}$	${0.87}_{-0.13}^{+0.18}$	${0.51}_{-0.04}^{+0.05}$
Feiden 2016	${0.54}_{-0.09}^{+0.12}$	${0.90}_{-0.15}^{+0.17}$	${0.41}_{-0.03}^{+0.04}$
Feiden 2016 (magnetic)	${0.78}_{-0.12}^{+0.11}$	${1.41}_{-0.09}^{+0.06}$	${0.69}_{-0.06}^{+0.08}$
Age [Myr]
Baraffe et al. 2015	${4.41}_{-1.06}^{+1.43}$	${1.77}_{-0.49}^{+0.60}$	${1.21}_{-0.10}^{+0.14}$
Bressan et al. 2012	${3.07}_{-0.62}^{+0.92}$	${1.26}_{-0.32}^{+0.41}$	${1.04}_{-0.06}^{+0.07}$
Chen et al. 2014	${6.76}_{-1.14}^{+1.10}$	${1.26}_{-0.23}^{+0.41}$	${1.62}_{-0.17}^{+0.19}$
Feiden 2016	${4.56}_{-0.95}^{+1.41}$	${1.41}_{-0.30}^{+0.38}$	${1.43}_{-0.08}^{+0.09}$
Feiden 2016 (magnetic)	${10.56}_{-2.50}^{+3.15}$	${3.51}_{-0.39}^{+0.38}$	${3.23}_{-0.37}^{+0.47}$

Equations

$$\Sigma (r)={\Sigma }_0{\left(\frac{r}{r_c}\right)}^{-\gamma }\exp \left[-{\left(\frac{r}{r_c}\right)}^{-(2-\gamma )}\right],$$

$${\Sigma }_0=\frac{(2-\gamma )M_d}{2\pi r_c^2}.$$

$$\frac{\partial \log \rho }{\partial z}=-\left[\left(\frac{G{M}_{\ast }z}{(r^2+z^2{)}^{3/2}}\right)\left(\frac{\mu m_H}{k_b{T}_g}\right)+\frac{\partial \log T_g}{\partial z}\right],$$

$$T_g(r)=T_0{\left(\frac{r}{1\mathrm{au}}\right)}^{-q}$$

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

$$h(r)={\left[\frac{k_b{r}^3{T}_g(r)}{G{M}_{\ast }\mu m_H}\right]}^{1/2}.$$

$$n_{\mathrm{CO}}(r,z)=X_{\mathrm{CO}}\frac{\rho (r,z)}{\mu m_H},$$

$$v_{\varphi }(r)=\sqrt{\frac{G{M}_{\ast }}{r}}.$$

$$\ell (\hat{\theta })=\exp \left[-\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^M{\left|V(u_i,v_i,{\nu }_j)-M(u_i,v_i,{\nu }_j|\hat{\theta })\right|}^{2}W(u_i,v_i,{\nu }_j)\right],$$

Full text

protoplanetary disks — stars: fundamental parameters — stars: pre-main sequence — radio lines: stars — parallaxes

††software: pdspy (Sheehan, 2018), CASA (McMullin et al., 2007), RADMC-3D (Dullemond, 2012), emcee (Foreman-Mackey et al., 2013), matplotlib (Hunter, 2007), corner (Foreman-Mackey, 2016), GALARIO (Tazzari et al., 2018)

1 Introduction

Protoplanetary disks are found nearly ubiquitously around pre-main sequence stars (e.g., Hernández et al., 2008), and studies of these disks are crucial to understanding the environment in which stars and planets form. Now, thanks to the revolutionary power of the Atacama Large Millimeter/submillimeter Array (ALMA), these disks have been studied in great detail. One notable development has been the unprecedented sensitivity to molecular line emission, which was inaccessible by previous generations of millimeter interferometers for all but the brightest sources(e.g., Simon et al., 2000). Access to gas emission has opened up a number of new avenues of study of protoplanetary disks, including their chemistry (e.g., Öberg et al., 2015; Schwarz et al., 2016), direct measurements of gas disk sizes and masses (e.g., Williams & Best, 2014; Ansdell et al., 2016, 2018; Miotello et al., 2017), as well as constraints on the turbulence in disks (Flaherty et al., 2015, 2017, 2018).

While there has been a significant effort toward studies of disk structure, the detection of gas emission has also created an opportunity to study the central star. Detailed radiative transfer modeling of the Keplerian motion in spatially and spectrally resolved disks has long been used as a tool for directly measuring stellar masses (e.g., Dutrey et al. 1994; Simon et al. 2000; Dutrey et al. 2003; Czekala et al. 2015, 2016; MacGregor et al. 2017; Ricci et al. 2017; Czekala et al. 2017; Wu & Sheehan 2017; Simon et al. 2017). However, the precision of these measurements is limited almost entirely by the uncertainty in the distance to the star, typically $>$10%, thereby limiting their usefulness for constraining evolutionary models. Although VLBI has produced precise distance measurements for some sources, potentially enabling precise mass measurements (e.g., Simon et al., 2017), care must be taken when extrapolating those distances (and distance uncertainties) to the star-forming region as a whole (e.g., Ortiz-León et al., 2017). In the era of precision distance measurements with $Gaia$, however, dynamical masses should be measurable with significantly improved precision for large samples of pre-main sequence stars.

This development comes at a welcome time as, while theoretical evolutionary tracks have progressed significantly over the last few decades, major discrepancies remain between the models and the direct measurements of stellar radii and masses (Hillenbrand & White, 2004; Gennaro et al., 2012; Stassun et al., 2014; Bell, 2016). Up to this point, constraints on pre-main sequence tracks have primarily come from eclipsing binaries, however these systems are relatively rare, with only $\sim$20 having accurate mass measurements (e.g., Stassun et al., 2014). Astrometry of young binary systems provides another avenue for measuring masses (e.g., Rizzuto et al., 2016; Dupuy et al., 2016; Rodet et al., 2018), but requires precise measurements over a number of epochs. These measurements are also limited, by definition, to binary stars. On the contrary many (or most) studies of young stars focus on single stars, which are easier to study. While it is typically assumed that relations between mass and temperature, radius, or other stellar parameters are common for binary and single stars, we must actually measure masses for a sample of single stars to test this assumption.

Constraints on evolutionary tracks are important for understanding fundamental stellar parameters (e.g., Stassun et al., 2014; Bell, 2016), but also because such models are often used to draw conclusions about the star and planet formation processes. Cluster ages inferred from evolutionary tracks have been used to place constraints on protoplanetary disk lifetimes (e.g., Haisch et al., 2001; Hernández et al., 2008). More recently, a large amount of effort has gone into studying the disk-mass to stellar mass relation, how it evolves with time, and what that evolution means for disk evolution and planet formation (e.g., Andrews et al., 2013; Pascucci et al., 2016; Barenfeld et al., 2016; Ansdell et al., 2016; Eisner et al., 2018). However, if stellar properties are derived from pre-main sequence evolutionary tracks, they are always model-dependent (e.g., Andrews et al., 2013; Ward-Duong et al., 2018).

Here we demonstrate how ALMA 12CO(2–1) observations and precise distances with $Gaia$ can be used measure high precision ($\sim$2%) dynamical masses for young stars, using observations of the three pre-main sequence stars 2MASS J16262774–2527247, CT Cha, and DH Tau. These sources were originally observed as part of a program to look for disks around young, accreting planetary-mass companions in both continuum and spectral line emission. The results of this program as it pertains to the companions are presented in Wu et al. (2017b); Wu & Sheehan (2017); Wu et al. (2017a). The star 2MASS J16262774–2527247 (hereafter 2M1626–2527) was observed due to a mix-up propagated through the literature that identified its coordinates as those of the $\sim$17 $M_{\rm Jup}$ planetary-mass companion ROXs 12 B (Bowler et al., 2017). Interestingly, 2M1626–2527 is also shown to be a companion to ROXs 12 A (Bowler et al., 2017). Although these sources are not truly “single”, this technique can be equally well applied to any star, single or multiple, with a disk.

The paper is structured as follows. In Section 2 we describe the ALMA observations of our three targets. In Section 3.1 we present our modeling analysis of the spectral line observations, including their mass measurements, using our new code pdspy111The pdspy code can be found here: https://github.com/psheehan/pdspy (Sheehan, 2018). Then, in Section 3.2 we use theoretical evolutionary tracks to calculate stellar masses for our sources, and compare with our dynamically measured masses in Sections 4 and 5. Finally we provide extensive appendices to further document the modeling technique as well as the pdspy code. pdspy can handle both spectral line and continuum observations, and is now publicly available. The continuum portion of the code was used in Sheehan & Eisner (2017) but the code was not formally introduced at that time.

2 Observations & Data Reduction

Our targets were observed with ALMA Band 6 during Cycle 3, in September 2016. We include details of each individual observation in Table 1. For each observation, the Band 6 receiver was set up with three basebands configured for wideband continuum observations, each with 2 GHz of bandwidth and centered at 233 GHz, 246 GHz, and 248 GHz. The final baseband was configured for 12CO (2–1) observations, with 3840 0.122 MHz channels centered at 230.538 GHz (0.32 km s*-1* velocity resolution; Hanning smoothed) in order to spectrally resolve the molecular line emission. In addition to observing our science targets, each track included time spent on quasars for the purposes of calibrating our data.

The data were calibrated using the ALMA Pipeline in the CASA software package (McMullin et al., 2007). Following the bandpass, flux, and gain calibrations, we employed a single iteration of phase-only self-calibration, with one solution per scan, to the continuum basebands and applied these phase solutions to our spectral line window. After calibration, the data were imaged and deconvolved using the spectral cube mode of the CASA clean task. The imaging was done using natural weighting to increase the sensitivity of the channel maps, as the emission is often faint, and with 0.4 km s*-1* channels. We provide additional details of the images in Table 1, and for brevity we show moment 0 and 1 maps in Figure 1 and channel maps with a limited selection of channels in Figures 2, 3, and 4. We do, however, show the full channel maps in Appendix C.

We note that there is a distinct lack of redshifted emission towards DH Tau from $\sim$6–7.2 km s*-1* (see Figure 1, 4, and Appendix C), with emission reappearing at 7.4 km s*-1*. We suspect this lack of emission arises from extinction by molecular cloud in that velocity range, as in the case of GO Tau (Schaefer et al., 2009). We therefore exclude those channels from the channel map modeling described in the next section. The channels that were not considered are marked with an “X” in any relevant figures. A similar effect (and resolution) were adopted by Czekala et al. (2016).

3 Analysis

3.1 Dynamical Stellar Masses with the pdspy Code

Keplerian rotation imprints a unique pattern in the disk kinematics that can be seen in channel maps (the so-called “butterfly pattern”), and that pattern is strongly dependent on stellar mass. We generate synthetic Keplerian disk channel maps using radiative transfer modeling and fit these synthetic observations to the actual data using a Markov Chain Monte Carlo (MCMC) fitting routine. We describe our disk model, and our fitting procedure, in more detail below.

We assume a radial structure for our disk model that is motivated by models of viscous disk evolution (Lynden-Bell & Pringle, 1974), with the surface density given by

[TABLE]

where $r$ is the stellocentric radius in cylindrical coordinates and $\gamma$ is the surface density power law exponent. $r_{c}$ is the radius beyond which the disk surface density is exponentially tapered, and serves as a proxy for the disk radius ($R_{disk}=r_{c}$) as the surface density profile has no hard outer limit. We assume that the disk is truncated at an inner radius of $R_{in}$. $\Sigma_{0}$ is a proportionality constant related to the total disk gas mass,

[TABLE]

As total disk gas mass is easier to interpret, we leave it as a free parameter in place of $\Sigma_{0}$.

The vertical structure of the disk is governed by hydrostatic equilibrium,

[TABLE]

where $\rho$ is the gas density, $M_{*}$ is the stellar mass, and $T_{g}$ is the gas temperature. Moreover, $\mu=2.37$ is the mean molecular weight of the gas, appropriate for a molecular gas with Solar metallicity, and $z$ is the distance in the vertical direction in cylindrical coordinates. For simplicity, we assume that the disk is vertically isothermal, with a radial temperature dependence of

[TABLE]

where $T_{0}$ is the temperature at 1 au and $q$ is the temperature power law exponent. Under these assumptions, solving the equation of hydrostatic equilibrium finds that the density structure is given by,

[TABLE]

where $h(r)$ is the scale height of the gas. The scale height from hydrostatic equilibrium is found to be,

[TABLE]

Finally, the number density of the CO gas can be derived from,

[TABLE]

where $X_{\rm CO}$ is the abundance of CO relative to H2. We typically assume that $X_{\rm CO}=1.0\times 10^{-4}$.

We also assume that the disk velocity structure is determined by Keplerian rotation, with an azimuthal velocity of,

[TABLE]

The velocities in the vertical and radial directions are assumed to be zero. We do, however, include microturbulent line broadening ($\xi$, in units of km s*-1*) as a parameter in our fit.

This density, temperature, and velocity structure can be input into the RADMC-3D radiative transfer code (Dullemond, 2012) to generate synthetic channel maps for a given set of model parameters. We assume local thermodynamic equilibrium to generate images, which is appropriate for disks as their densities are much higher than the critical density for CO of $\sim 10^{3}$ cm*-3*. We also leave a number of viewing orientation parameters as free parameters in our fit: inclination ($i$), position angle ($P.A.$), centroid ($x_{0}$, $y_{0}$), systemic velocity ($v_{sys}$), and source distance ($d_{pc}$). The synthetic channel maps generated by raytracing in RADMC-3D are Fourier transformed to produce synthetic visibilities, sampled in the $uv$-plane at the same baselines as our datasets, using the GALARIO code (Tazzari et al., 2017, 2018).

In all, our model includes 14 free parameters: $M_{*}$, $M_{disk}$, $R_{in}$, $R_{disk}$, $\gamma$, $T_{0}$, $q$, $\xi$, $v_{sys}$, $i$, $d_{pc}$, $P.A.$, $x_{0}$, and $y_{0}$. In practice, however, our modeling code is flexible, and allows additional model features and parameters to be included, for example gaps in the gas density distribution, but here we restrict the model to these 14. We provide further details on additional parameters in Appendix B.

We fit the model channel map visibilities generated by RADMC-3D and GALARIO to our data in the $uv$-plane using the MCMC fitting code emcee (Foreman-Mackey et al., 2013). For brevity, and because emcee has been broadly used by the community, we opt to provide further details on the setup of our MCMC runs, and details of the convergence in Appendix A. The best-fit values are determined by discarding the burn-in MCMC steps and calculating the median of the sample of post-burn-in walker positions for each parameter. We also estimate the uncertainty on the measured parameters by considering the range around the median value that contains 68% of the post-burn-in walker positions. We provide further details on this process, as well as posterior probability distribution function triangle plots and plots of the walkers’ steps in each parameter, in Appendix A.

3.2 Masses and Stellar Parameters from Evolutionary Tracks

Although masses have been previously estimated for 2M1626–2527, CT Cha, and DH Tau using evolutionary models, they were done across separate studies that do not always use uniform methodologies. As such, we find it worthwhile to repeat those measurements in a uniform way, following the method outlined in Andrews et al. (2013) and Czekala et al. (2015).

To begin, we derived the effective temperature ($T_{eff}$) and bolometric luminosity ($L_{*}$) by fitting optical and near-infrared photometry collected from the literature (Vrba et al., 1989; Epchtein et al., 1999; Briceño et al., 2002; Ducati, 2002; Cutri et al., 2003; Monet et al., 2003; Zacharias et al., 2003; Torres et al., 2006; Robberto et al., 2012; Muiños & Evans, 2014; Henden et al., 2016; Bowler et al., 2017) with the BT-Settl atmospheric models (Allard et al., 2011) using the emcee package (Foreman-Mackey et al., 2013). For any set of input parameters, the model grid was interpolated to match $T_{eff}$ and surface gravity (log $g$). The overall brightness was scaled by the distance to the source ($d_{pc}$) and the radius of the star ($R_{*}$), and then the model was reddened by some amount ($A_{V}$) using the extinction curves from Weingartner & Draine (2001) with $R_{V}=5.5$. We assumed a Gaussian prior on $d_{pc}$ derived from $Gaia$ parallaxes (see Appendix A). Similarly we assumed a Gaussian prior on $T_{eff}$, using temperatures collected from the literature (Herczeg & Hillenbrand, 2014; Bowler et al., 2017; Manara et al., 2016) with a standard deviation of one spectral sub-class (150 K). Finally, the stellar luminosity ($L_{*}$) was calculated from $T_{eff}$ and $R_{*}$.

For each post-burn in MCMC sample, we mapped $T_{eff}$ and $L_{*}$ on to evolutionary tracks to build up a probabilistic sample of stellar masses ($M_{*}$) and ages ($\tau_{*}$). The stellar mass and age, and the associated uncertainties, were estimated from these samples using the median and range around the median containing 68% of the samples. This process was repeated for each of the sets of evolutionary tracks that we considered: the Baraffe et al. (2015) models, the PARSEC models (Bressan et al., 2012) and subsequent updates to low-mass ($<0.75~{}M_{\sun}$) stars (Chen et al., 2014), as well as both the non-magnetic and magnetic tracks from Feiden (2016).

4 Results

We show the best fit Keplerian disk models, as determined by our spectral line fitting procedure, to 2M1626–2527, CT Cha, and DH Tau in Figures 2, 3, and 4. In these figures we show only a subset of channels to demonstrate that our models provide excellent fits to the data, but the comparison of the models to the full channel map data can be found in Appendix C. For all three sources the models provide good fits to the data, with no (or few) significant, $>3\sigma$, residuals. The best fit parameters found from these fits are listed in Table 2. All three disks are small, with characteristic radii of $\lesssim 15$ au, although due to the exponential tapering our disk models do contain a small amount of mass at larger radii.

Table 3 lists the masses and other stellar parameters derived from the stellar photosphere models, and Figure 5 shows the best fit models compared with the data. Figure 6 demonstrates the process of mapping $T_{eff}$ and $L_{*}$ on to the evolutionary tracks from Baraffe et al. (2015) to estimate $M_{*}$ and $\tau_{*}$. The posterior distributions for stellar mass and age for each set of evolutionary tracks are shown in Figure 7, with a gray line representing the mass we derive from our spectral line fits.

We find that the dynamical masses we measure for 2M1626–2527 and CT Cha are in reasonable agreement with masses inferred from evolutionary tracks, although we discuss the agreement for each source in more detail below. The magnetic tracks from Feiden (2016), however, seem to universally over-estimate the masses for all three sources.

4.1 2M1626–2527

We find from our spectral line modeling that 2M1626–2527 has a dynamical mass of $0.535^{+0.006}_{-0.007}$ $M_{\sun}$, in very good agreement with the mass estimated by the Baraffe et al. (2015) and non-magnetic Feiden (2016) evolutionary tracks. Interestingly, though, the original PARSEC models by Bressan et al. (2012) under-predict the mass by $\sim 1\sigma$ while the revision to the PARSEC models by Chen et al. (2014), intended to improve the tracks for stellar masses $<0.75$ $M_{\sun}$, over-predict 2M1626–2527’s mass by $\gtrsim$ 3$\sigma$. Our mass measurement is also consistent with the value of 0.5 $\pm$ 0.1 $M_{\sun}$ derived by Bowler et al. (2017), who also used the Baraffe et al. (2015) tracks.

We also find that the spin axis of 2M1626–2527 is probably not orthogonal to the disk plane (disk inclination of 54° versus stellar inclination of 9°–27° found by Bowler et al. (2017)). A similar misalignment of disk and star was seen in our previous study of GQ Lup (Wu et al., 2017b).

4.2 CT Cha

CT Cha has a stellar mass of $0.796^{+0.015}_{-0.014}$ $M_{\odot}$, which, contrary to 2M1626–2527, is in very good agreement with both sets of PARSEC evolutionary tracks. The non-magnetic Feiden (2016) tracks are also in reasonable, $\lesssim 0.5\sigma$, agreement, while the Baraffe et al. (2015) evolutionary tracks predict a mass that is $\sim 1\sigma$ too high. Furthermore, the mass is in reasonable ($\lesssim 0.5\sigma$) agreement with the mass of $0.87^{+0.23}_{-0.18}~{}M_{\odot}$ estimated by Pascucci et al. (2016) using a combination of the Baraffe et al. (2015) and Feiden (2016) tracks.

4.3 DH Tau

We find that DH Tau has a mass of $0.101^{+0.004}_{-0.003}$ $M_{\odot}$, significantly lower than 0.3–0.5 $M_{\odot}$ we derive from evolutionary track models. It is also significantly lower than the 0.2–0.5 $M_{\odot}$ estimated by Hartigan et al. (1994), White & Ghez (2001), and Andrews et al. (2013). However, it is possible that our modeling is affected by the strong molecular cloud absorption obscuring much of the low velocity emission. In order to test whether foreground extinction could be affecting our mass measurement for DH Tau, we add a simple model for the extinction to our fit. We assume that each channel is extincted by an amount of $\exp{(-\tau_{\nu})}$. We assume that the optical depth, $\tau_{\nu}$, is Gaussian in shape with some central velocity, $v_{0,ext}$, width $\sigma_{ext}$, and peak value, $\tau_{0}$. We otherwise use the same model described above and fit the model to the data using the same methods, but with three additional parameters.

We find that this new model does increase the stellar mass to $0.141^{+0.016}_{-0.011}$ $M_{\sun}$. The difference in mass can be attributed to the lower inclination found by this model, of $38.0^{+2.4}_{-2.5}$ degrees. The best fit values of the new parameters are $v_{0,ext}=6.629^{+0.016}_{-0.018}$ km s*-1*, $\sigma_{ext}=0.389^{+0.027}_{-0.026}$ km s*-1*, and $\tau_{0}=6.15^{+1.71}_{-1.01}$, and all other parameters are consistent with the values we find for our base model. Although this simple model cannot fully reconcile the mass discrepancy between our measurement and the evolutionary tracks, it is possible that a more complicated extinction profile could increase the inferred stellar mass even more.

5 Discussion

The dynamical masses we measure from our Keplerian disk fits demonstrates the power of combining ALMA spectral line maps with precise distance measurements with $Gaia$. Previous Keplerian disk dynamical mass measurements exhibited extreme degeneracies between stellar mass and source distance (e.g., Czekala et al., 2015). As uncertainties in the source distance could be as large as 10–20% or more, they translated directly into 10–20% uncertainty on the measured stellar mass. Here, using precise estimates of distances from $Gaia$, our stellar mass estimates have uncertainties of only 2-4%.

As the uncertainty on $Gaia$ distances is $\lesssim 1\%$, it seems likely that these uncertainties are no longer dominated by the distance uncertainty. Indeed, the posterior probability density function plots shown in Appendix A show little degeneracy between stellar mass and source distance for DH Tau and CT Cha, suggesting that the measurement is limited primarily by the quality of the ALMA data. The ALMA data for 2M1626–2527 are good enough that the $Gaia$ distance uncertainty does affect the mass measurement, noticeable as slight a degeneracy between mass and distance, but the $Gaia$ uncertainty is only $\sim 1\%$ so it does not prevent us from making a high precision measurement of the stellar mass.

Interestingly, the comparison of our dynamical mass measurements with evolutionary track mass estimates shows a wide range in agreement. The PARSEC evolutionary tracks (Bressan et al., 2012; Chen et al., 2014) do the best job of reproducing the highest mass source in our sample, CT Cha, but fail to correctly reproduce the mass of 2M1626–2527. Conversely, the Baraffe et al. (2015) evolutionary tracks do an excellent job of reproducing the mass of 2M1626–2527, but overestimate the mass of CT Cha. The non-magnetic tracks from Feiden (2016) may provide the best balance of matching the masses of both CT Cha and 2M1626–2527, but do not match both perfectly. It is clear, though, that the magnetic tracks from Feiden (2016) do a poor job of reproducing the masses of both sources.

The case of DH Tau is particularly interesting because the mass we measure is discrepant from evolutionary track estimates, which are uniformly ($>4\sigma$) too high. As we demonstrated earlier, foreground extinction could affect our stellar mass measurement, although our simple model for extinction only increases the stellar mass to $\sim$0.14 $M_{\sun}$ and does not reconcile with evolutionary tracks. Of course, a more complicated extinction profile may increase the stellar mass by more. If the mass estimate is correct, however, then a relatively major shift in the evolutionary tracks for young, low-mass objects is needed. Higher signal-to-noise data may help to better understand the effects of cloud contamination on this source and better constrain its mass, as might more optically thin gas tracers that better see through the foreground cloud. For now, however, it remains unclear whether the discrepancy is a result of the limitations of the data, or whether there are real problems with the evolutionary tracks.

Moving forward, however, a much larger sample of pre-main sequence stars is needed to understand whether the trends suggested by our targets are more generally true, and thereby truly place constraints on evolutionary tracks. Our results here demonstrate that with relatively modest observations ($\sim$15 minutes on source; $0\farcs 2-0\farcs 3$ resolution; 0.4 km s*-1* channels) we can robustly measure stellar masses as low as $\sim$0.1 $M_{\sun}$. Moreover, with channel widths as small as 0.03–0.05 km s*-1* for CO (2–1) or CO (3–2), it should be possible to directly measure masses for young sub-stellar and planetary mass objects. Building a sample with several tens of pre-main sequence sources could therefore reasonably be done.

6 Conclusions

We have used our Keplerian disk radiative transfer code psdpy to fit ALMA 12CO (2–1) channel maps for three T Tauri stars 2M1626–2527, CT Cha, and DH Tau. With new, precise distance measurements for these sources from $Gaia$, we are able to make high precision ($\sim$2–4%) direct measurements of their masses. The mass uncertainty is no longer dominated by distance as in previous studies, but instead by the quality of ALMA data. We find that 2M1626–2527 has a mass of $0.535^{+0.006}_{-0.007}$ M*⊙, CT Cha has a mass of $0.796^{+0.015}_{-0.014}$ M⊙, and DH Tau has a mass of $0.101^{+0.004}_{-0.003}$ M⊙*. Comparing with the stellar masses estimated using multiple sets of evolutionary tracks, we find that both 2M1626–2527 and CT Cha are in reasonable agreement with most models ($<$ 2$\sigma$), although potentially significant differences remain for most sets of tracks. The large discrepancy for DH Tau may be caused by a significant local cloud obscuration, which completely absorbs the redshifted disk emission, although simple models of the foreground extinction only increase the mass to $0.141^{+0.016}_{-0.011}$ $M_{\sun}$.

Looking towards the future, the synergy of ALMA, $Gaia$, and detailed radiative transfer disk modeling will enable precise mass measurements of a large number of pre-main sequence stars, and even for substellar objects, and ultimately place stringent constraints on evolutionary models.

We would like the thank the referee for a number of comments that helped to significantly improve this work. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.00773.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Y.-L. W. is grateful to the support from the Heising-Simons Foundation. J.T. would like to acknowledge support from the Homer L. Dodge Endowed Chair. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. 2012115762. This work was supported by NSF AAG grant 1311910. The results reported herein benefitted from collaborations and/or information exchange within NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate. The computing for this project was performed at the OU Supercomputing Center for Education & Research (OSCER) at the University of Oklahoma (OU).

Appendix A The MCMC Fitting Procedure

We fit our radiative transfer model, as described in Section 3.1, to our data using the Markov Chain Monte Carlo (MCMC) fitting code emcee (Foreman-Mackey et al., 2013). emcee uses an implemenation of the Goodman & Weare (2010) Affine-Invariant ensemble sampler to sample parameter space with a large number of walkers. We use 200 separate MCMC “walkers”, which we spread out over a large range of parameter space. These walkers are allowed to move through parameter space, with their steps being towards another randomly selected walker and the size of the step dictated by a comparison of the goodness of fit of each set of walker parameters. The likelihood function for a given set of parameters $\hat{\theta}$ used for the fit is given by

[TABLE]

where $V(u_{i},v_{i},\nu_{j})$ are the complex visibility data at baseline $(u_{i},v_{i})$ and frequency $\nu_{j}$, $W(u_{i},v_{i},\nu_{j})$ are the visibility data weights, and $M(u_{i},v_{i},\nu_{j}\,|\,\hat{\theta})$ is the model with parameters $\hat{\theta}$ evaluated at baseline $(u_{i},v_{i})$ and frequency $\nu_{j}$. $N$ is the total number of baselines in the visibility data, and $M$ is the number of channels included in the fit. We consider a fit to be “converged” when the median walker position does not change significantly compared to the spread of the walkers over a significant number of steps. We show the steps taken by the walkers in each fit in Figures A1, A2, and A3.

We note that while we commonly refer to the 14 parameters of our model fit as $\hat{\theta}=\{M_{*},M_{disk},R_{in},R_{disk},$ $\gamma,T_{0},q,\xi,v_{sys},i,d_{pc},P.A.,x_{0},y_{0}\}$, reasonable ranges of values for a number of these parameters span multiple orders of magnitude. For such parameters, it makes practical sense to use the log of the parameter as the actual fit quantity. For this reason, the actual parameters being fit are $\hat{\theta}=\{\log_{10}{M_{*}},\log_{10}{M_{disk}},\log_{10}{R_{in}},\log_{10}{R_{disk}},\gamma,\log_{10}{T_{0}},q,\log_{10}{\xi},$ $v_{sys},i,d_{pc},P.A.,x_{0},y_{0}\}$. As the linear values are more intuitive to understand however, for most tables and figures and discussion in the text we convert the posterior Markov Chain Monte Carlo samples to linear values before evaluating best fit values. We assume a uniform prior on all parameters except for the source distance ($d_{pc}$) and stellar mass ($M_{*}$), with the following limits: $-2\leq\log_{10}{M_{star}}\leq 1$, $-10\leq\log_{10}{M_{disk}}\leq 0$, $-1\leq\log_{10}{R_{in}}\leq 1.5$, $0\leq\log_{10}{R_{disk}}\leq 2.5$, $-0.5\leq\gamma\leq 2$, $1\leq\log_{10}{T_{0}}\leq 3$, $0\leq q\leq 1$, $-1.5\leq\log_{10}{a_{turb}}\leq 1$, $0^{\circ}\leq i\leq 90^{\circ}$, $0^{\circ}\leq p.a.\leq 360^{\circ}$, $-0.1"\leq x_{0}\leq 0.1"$, $-0.1"\leq y_{0}\leq 0.1"$. In principle, sensitive ALMA spectral line data can distinguish between inclinations above and below $90^{\circ}$ (e.g., Czekala et al., 2015). However our initial tests found that the data does not have the required sensitivity, with the walkers split approximately evenly between the two solutions, so we limit $i<90^{\circ}$.

For a young star for which the stellar mass is not a priori known, the probability of that star having mass $M_{*}$ is given by the stellar initial mass function. As such we use the Chabrier IMF (Chabrier, 2003) as a prior on the stellar mass in our fit. Furthermore, while $d_{pc}$ is left as a free parameter in the fit, our ALMA data do not provide a direct constraint on its value. As stellar mass is highly degenerate with source distance (e.g., Czekala et al., 2015), placing some constraint on source distance is important. Moreover, strong constraints on source distance lead to higher precision measurements of stellar mass. Here, we use direct measurements of source distances from trigonometric parallaxes by $Gaia$ to place a prior constraint on the distance. We find from $Gaia$ that 2M1626–2527 has a parallax of $7.232\pm 0.074$ mas and a distance of $138.258^{+1.426}_{-1.396}$ pc, CT Cha has a parallax of $5.214\pm 0.021$ mas and a distance of $191.775^{+0.781}_{-0.774}$ pc, and DH Tau has a parallax of $7.388\pm 0.069$ mas and a distance of $135.355^{+1.282}_{-1.258}$ pc (Gaia Collaboration et al., 2018). We assume a Gaussian prior on the parallax, with a mean given by the measured $Gaia$ parallax and the standard deviation given by the $Gaia$ uncertainty. The model parallax is calculated as $1/d_{pc}$.

While in general the walkers tend to converge into one bunch, a handful of walkers appear to become lost in parameter space and have trouble finding their way to the main group (e.g., see Figure A1). We believe that this is the result of how the walkers move through parameter space. To determine the next step for a walker, a different walker is randomly drawn from the sample of walkers. The walker will (randomly) move towards or away from the randomly selected walker with a step size between 1/a and a times the distance between them (a = 2 by default in emcee). The walkers become “lost”, however, because the global minimum is quite narrow in some parameters (e.g., $v_{sys}$), so the likelihood that the step puts it into the main bunch (i.e. a = 1) gets very small. It therefore takes a large number of step proposals before the walker finally makes the jump. While we have confidence that, if given enough time, these walkers would converge as well, the computational demands of these fits are significant, so we opt not to continue running the fits. Instead, we discard any walkers with $(\chi^{2}-\chi^{2}_{min})>44$, appropriate for a $\chi^{2}$ distribution with 14 parameters and $p=0.00006$ ($\sim 4\sigma$). The autocorrelation time for these fits is typically about $\sim 1000$ steps, so with $5000$ post-burn-in steps in the chain, we have $\sim 5$ independent sets of walker positions. With 200 walkers, or a total of $\sim 1000$ independent samples, we would expect $\ll 1$ ($\sim 0.06$) walkers to have $(\chi^{2}-\chi^{2}_{min})>44$, so this cut will identify any “lost” walkers that are not likely representative of the underlying distribution. We show the steps taken by all walkers in Figures A1, A2, and A3. The discarded walkers shown in gray to demonstrate that the cuts are reasonable.

In Figures A4, A5, and A6 we show the two dimensional posterior probability density functions for each pair of parameters, after the burn-in steps and “lost” walkers have been discarded, for each source. We estimate the best fit value for a parameter as the median of the walker positions in that parameter after convergence is achieved, the “lost” walkers are discarded, and the walker steps are trimmed by the autocorrelation time. Similarly, we estimate the uncertainties on that value as the region around the median containing 68% of the walker positions after convergence and “lost” walkers are discarded. For parameters where the walkers move up against an edge of the allowed range ($a_{turb}$ for CT Cha and DH Tau, $R_{disk}$ for DH Tau and 2M1626–2527, $M_{disk}$ for CT Cha, $R_{in}$ for all), we instead use the value containing 99.7% of walkers above or below that value as a limit.

Appendix B Additional Model Parameters

pdspy is designed to be a flexible modeling tool that can handle a range of input parameters, priors, and models. For example, here we have left the surface density profile power law exponent, $\gamma=1$, as a free parameter. However, the walker plots show that $\gamma$ converges slowly and does not ultimately affect the measured value of $M_{*}$. As such, it may be wise to fix $\gamma=1$ (e.g., Czekala et al., 2015) to reduce the needed computational time, and this can easily be done in $pdspy$ by changing a single boolean value in the configuration file. This, of course, may not be wise if the goal is to measure disk structure, as parameters like $M_{disk}$ and $R_{disk}$ are degenerate with $\gamma$, but for studies measuring stellar mass this would be appropriate. In addition to adding or removing parameters from our model, pdspy also includes a variety of models that can be employed. Here we have employed an exponentially tapered surface density profile, however pdspy includes options for a truncated disk model as well. Moreover, up to three gaps and a cavity can be added as features of the model. Finally, simple Gaussian priors can be added to any parameter, and a few “special” priors are included in the code. As described above the distance can be constrained with a prior on the parallax, and the stellar mass can be constrained by using the IMF as a prior.

In all, pdspy currently has over 25 parameters that can be turned on or off at the users discretion, as well as options to use a truncated disk and to include continuum opacity/subtraction. We are also exploring adding additional functionality, including non-vertically isothermal disk temperatures and inclusion of envelope emission for embedded protostars. More information, including the current, up-to-date status of modeling options, can be found at https://github.com/psheehan/pdspy.

Appendix C Channel Map Images and Models

In Figures C1, C2, and C3 we show the best fit models and residuals for each of our sources, here including every channel in the maps.