A Case of Simultaneous Star and Planet Formation

BHB1 was observed with ALMA at 226 GHz as part of project code 2013.1.00291.S (PI: F. O. Alves). The total continuum-dedicated bandwidth is 2.4 GHz. CASA 4.5 was used for calibration and imaging, where details of the calibration can be found in Alves et al. (2017). These observations used 44 antennae with baseline ranging from 15 to 1460 m. In this configuration, the maximum recoverable angular scale is ∼21''. We performed self-calibration of the uv visibilities by interpolating over decreasing solution intervals the phase gains (scan and integration time intervals) and amplitude gains (scan intervals). The reprocessed image has an rms noise of 0.095 mJy beam⁻¹, intensity peak of 19.5 mJy beam⁻¹ and flux density of 0.16 ± 0.02 Jy. The continuum map was produced using Briggs robust parameter 0.5 and has a synthesized beam of 024 × 020 and position angle (PA) of 77° (east of north).

¹²CO and C¹⁸O J = 2 → 1 were simultaneously observed covering the source and its surroundings. The spectral line data have velocity resolution of 0.35 km s⁻¹ (∼270 kHz) and peaks at channel −0.95 km s⁻¹, where intensity reaches 200 mJy beam⁻¹. Images were produced using a 0.5 Briggs robust parameter. The final map has an rms noise of 3.5 mJy beam⁻¹ and a synthesized beam of 025 × 020 (PA  ∼ 72°). Both continuum and spectral line maps were primary beam corrected, important because the source is offset by  ∼ 46 from the phase center.

We have used the VLA in its most extended configuration to observe the 22.2 GHz continuum (K-band, λ = 1.35 cm). A detailed description of the correlator setup is found in Alves et al. (2019), who reported on a distinct source with the same set of VLA observations used for BHB1. The VLA image was obtained using CASA clean and a robust weighting of 1, yielding a synthesized beam of 0197 × 0087 (≃21 au) with a PA of 103 and an rms noise of 7.5 μJy beam⁻¹.

The disk appears clearly in the millimeter data, with a quasi-symmetric morphology and three distinctive peaks, interpreted as an edge-on view of a gapped disk (Figure 1(a)). The disk has a radius of 107 au and is oriented  ∼15° east of north. The inner and outer disk are separated by a gap with ∼70 au width. The mm dust brightness temperature peaks at ∼10 K.

Zoom In Zoom Out Reset image size

The disk's molecular gas traced by ¹²CO is more extended than the mm emission (Figure 1(b)). The spatial distribution of the ¹²CO and C¹⁸O velocity components show a clear Keplerian rotation pattern and a flared morphology. The ¹²CO has cloud contamination in the central channels; however, the C¹⁸O (Figure 1(c)) has less contamination and we find that the source velocity is consistent with previous estimates of 3.6 km s⁻¹ (Onishi et al. 1999). Using the uncontaminated channels and following the prescription of Seifried et al. (2016), we fit a Keplerian model to the position–velocity diagram of the CO emission and find a mass of 2.23 ± 0.04 M_⊙ for the embedded protostar.

Molecular emission reveals large-scale (∼4000 au) and narrow (∼80–300 au) bipolar filaments connecting the ambient gas and the source (Figure 1(d)). These structures are distributed in a north–south orientation, similar to the major axis of the disk. The northern (redshifted) and southern (blueshifted) filaments exhibit a velocity shift of ∼±2–3 km s⁻¹ with respect to the source ambient velocity.

No clear velocity gradients are seen along the filaments, whose entire structures are visible over a ∼1 km s⁻¹ velocity range in each lobe. This low-velocity width implies we are likely seeing the filament along the plane of the sky. Interestingly, the pair of filamentary structures have global redshifts and blueshifts. The northern filament is redshifted as the disk is with a similar velocity. The same is true for the southern filament, which is largely blueshifted, consistent with the disk rotation on the southern side of the disk.

If these features are indeed filaments, they appear to be orbiting the system like large propellers or gas streamers accreting into the disk (Section 5.1). Alternatively, we could be seeing a limb-brightened large-scale (>1000 au) outer disk or flattened envelope, where the central channels are too optically thin to observe in emission. Additional observations of higher critical density tracers with deeper observations would be needed to disentangle a filament or remnant-flattened envelope. What we can say with certainty is that we see large-scale emission associated with the disk based on its velocities that is moving too slowly to be an outflow, and is oriented parallel to the disk plane rather than along a conventional outflow axis.

The VLA data show a compact radio source at the protostellar position with a flux density of 103 ± 12 μJy (we refer to this source as VLA 1). In addition, the data reveal a marginal second peak in the southern gap at a separation of 03 (∼49 au) from the star. This tentative source has a flux density of 30 ± 8 μJy (for simplicity we refer to it as VLA 1b). The northern and southern gaps seen in the ALMA continuum map show intensity minima (5.21 and 5.28 mJy beam⁻¹) that within uncertainties (0.10 mJy beam⁻¹) are the same , i. e., there are no traces of VLA 1b at 1.3 mm. This puts an upper limit for the spectral index of VLA 1b of  ≲0.8 (assuming a 2σ upper limit at 1.3 mm), indicating that VLA 1b traces ionized gas.

Interestingly, VLA 1b coincides with a bright (∼100 K) and compact ¹²CO emission hot spot (Figures 2(a), (b)). This spot is probably localized in the southern gap and VLA 1b because (1) it is spatially compact and with a narrow velocity width—both features indicate that it is truly compact in three dimensions—and (2) the peak velocity of this feature is close to the highest velocity at the projected radius (see Figures 2(c), (d)), i.e., the spot arises at a radius similar to the projected radius on the plane of the sky.

Zoom In Zoom Out Reset image size

To estimate the size of the gap seen at mm-wavelengths, we explore parametric axisymmetric models using the TORUS radiative transfer code (Harries 2014). As the goal of this is to approximately constrain the size of the gap, we have adjusted the model parameters until we achieve a good fit by eye to the mm data and leave more complex modeling to a future analysis.

The functional form of the disk model follows the Lynden-Bell & Pringle (1974) formalism for a viscously evolving disk, where

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{g}(R)={{\rm{\Sigma }}}_{c}{\left(\displaystyle \frac{R}{{R}_{c}}\right)}^{-\gamma }\exp \left[-{\left(\displaystyle \frac{R}{{R}_{c}}\right)}^{2-\gamma }\right].\end{eqnarray} \tag{ 1 }$$

For the gas, Σ_c = 15 g cm⁻², R_c = 100 au, and γ = 1, and we do not include gas gaps. The subscript c indicates the critical radius at which point the surface density profile transitions from being power-law dominated to falling off with an exponential taper. We assume a vertically Gaussian density distribution with a scale-height h of 12 au at 100 au, and a flaring parameter of ψ = 1.2 where $h={h}_{100}{\tfrac{{R}_{\mathrm{au}}}{100}}^{\psi }$ and h is in au. For the dust, we assume a global gas to dust ratio of 100; however, the local gas/dust ratio in surface densities varies radially. The dust is split into two MRN population distributions with sizes up to 1 micron and 1 mm (see Cleeves et al. 2016), and both adopt pure astrosilicate compositions (Draine & Lee 1984). The small dust grains follow the gas and contain 1% of the dust mass. The mm-grains have the same functional form as Equation (1), but have a power law of γ = 0.2 and an R_c of 1000 au, well beyond the physical disk edge such that the mm dust surface density has a truncated power-law behavior. The large grains' scale-height is 10% of the gas/small dust to approximate settling. The inner disk edge is set to 0.1 au and extends to 20 au. We model an empty gap, and then an outer disk component starting at 93 au extending out to 160 au where it is truncated. We do not vary the depth of the gap as we only have one thermal continuum wavelength, and thus it is challenging to make accurate constraints on both the dust optical properties and the minimum dust mass in the gap, especially given its near edge-on nature. Instead we emphasize that this analysis is aimed to estimate the size of the gap, which we use to estimate possible companion masses for a range of uncertain physical properties (Section 4.2).

Based on this structure and the stellar parameters described in Sections 1 and 3.2, TORUS computes a synthetic image at 1.287 mm, which we convolve with the observed beam. We explored just over 100 models and arrived at a reasonable fit as shown in Figure 3. To arrive at this fit, we varied all parameters except the γ value of the gas, R_c for the gas, the inner edge of the inner disk from 0.1 au, the relative dust mass between small and large grains, and the relative scale height of the large grains compared to the small. We find that within our grid of models, the inner edge of the outer ring can reasonably reproduce the structure qualitatively within ±5 au. The gas+dust mass of this disk model is 0.1 M_⊙.

Zoom In Zoom Out Reset image size

Taking this structure, we have additionally run a few non-local thermodynamic equilibrium (non-LTE) LIME (Brinch & Hogerheijde 2010) models with our TORUS-derived structure to qualitatively investigate how much gas could be missing in the gap but not be clearly seen in the optically thick ¹²CO. Given the lower signal to noise of the C¹⁸O observations, we do not try to compare simulations with these data. We assume a simple uniform CO abundance of 10⁻⁴ relative to H above dust temperatures of 20 K and 10⁻¹⁰ per H below similar to the approach of Qi et al. (2006) to approximate CO freeze-out and photodissociation. We have artificially decreased the gas inside of the gap by varying factors. We find that a gap in the gas distribution is only visible when the gap depth is higher than a factor 100, otherwise it remains hidden due to the high optical depth of the ¹²CO emission. We emphasize that the true gap depth would require more detailed fitting, ideally with thermo-chemical models, as the physical properties of gas in the gaps can be different depending on the gas to dust ratio in the gap, the presence or absence of small grains in the gap, the stellar ultraviolet heating of gas, and many more parameters (e.g., Bruderer et al. 2012; Bruderer 2013; Bruderer et al. 2014; van der Marel et al. 2015, 2016, 2018; Facchini et al. 2018) that are not constrained here. Instead, these simple models tell us that there is some amount of gas still in the gap and it is not fully cleared. Future deeper observations combined with thermo-chemical simulations should be conducted to better constrain the degree of gas depletion in the gap, if it exists.

Such a wide gap in the dust may indicate the presence of a single or multiple planets. To estimate the mass of the potential planet, we adopt the approach used in the DSHARP survey (Zhang et al. 2018). First, based on the gap profile, we estimate the relative width of the gap (Δ parameter in Zhang et al. 2018, which is the gap width over the radius of the outer gap edge) as 0.78. If we adopt DSD1 dust size distribution in DSHARP (dust size follows a power-law distribution from 0.005 μm to 0.1 mm with a power-law index of −3.5) and the gas surface density derived above, the Stokes number of the biggest dust particle at 50 au is 8.6 × 10⁻⁴ and the derived K' parameter through the width-mass fitting (defined in Table 1 of Zhang et al. 2018) is 0.40. Thus, with ${\left(h/r\right)}_{50\mathrm{au}}$ = 0.1, we can derive that the planet-to-star mass ratio is 0.031 if the disk viscosity coefficient is α = 10⁻³ and 0.015 if α = 10⁻⁴. If we adopt DSD2 dust size distribution (dust size follows a power-law distribution from 0.005 μm to 1 cm with a power-law index of −2.5), the Stokes number at 50 au is 0.086 and the derived K' parameter is 0.05. The planet-to-star mass ratio is 0.0039 if α = 10⁻³, and 0.0019 if α = 10⁻⁴.

Taking into account all of these uncertain disk physical parameters, the possible planet-to-star mass ratio spans a wide range of [0.0019, 0.031]. Adopting the stellar mass computed in Section 3.2, the mass of the planet ranges between ∼4 and 70 M_Jup. We note that the dominant uncertainty in this determination is the maximum particle size and disk viscosity, rather than the details of the fit in Section 4.1. While the range is large, all masses point to a super-Jupiter sized companion, which would be consistent with the mass needed to launch jets capable of producing substantial free–free radio emission (Zhu et al. 2018); see also Section 3.4.

The large-scale filaments are aligned with the disk major axis and coherent with the protostellar disk's Keplerian velocity seen in CO. In addition, no evidence of acceleration at along the filament makes it unlikely that the observed gas traces outflowing material.

A more plausible explanation is that the filaments are streamers between the bulk of gas in the B59 core in the southeast and BHB1. The lack of a clear velocity gradient along the filaments suggests that the gas motion occurs mostly in the plane of the sky. Correspondingly, it is hard to assess which direction the gas flows; however, such low-velocity gas is typical of infalling streamers observed toward accreting protostars (Mottram et al. 2013). In this context, the filaments could be molecular gas nurturing the disk and the YSO, which is consistent with single-dish observations of B59 showing dust and molecular patches protruding from the core toward the west, where BHB1 is located (Duarte-Cabral et al. 2012).

Figure 4 displays the line-of-sight (LOS) velocity of CO gas as a function of distance to the protostar. At larger distances, the LOS velocity increases from ∼2000 to 1750 au at the redshifted lobe, and from ∼1650 to 1500 au at the blueshifted lobe. The LOS component then either decreases or becomes nearly constant at shorter distances. This velocity shift can be interpreted as material initially falling along the LOS, changing its direction and finally moving toward the source along the plane of the sky. The region of increasing velocity is seen in Figure 1(d) as bends in the endpoints of the filaments in both lobes. The bend is especially clear in the less-filtered redshifted lobe. Infalling gas parcels departing from rest with no energy dissipation have their LOS velocities determined as ${V}_{{LOS}}=\sqrt{2{{GM}}_{\star }/d}\times \cos \alpha$, where d is the distance to the protostar, α is the angle between the velocity vector and the LOS, and M_⋆ is the protostellar mass (Figure 4, inset). Large-scale accretion streamers have been recently reported by Pineda et al. (2020), who presented a scaled-up version of our filaments feeding a protostellar core. This indicates that the filamentary accretion streamers are the natural process in accreting objects regardless their age, as also indicated by time-dependent multi-scale simulations (Kuffmeier et al. 2017).

Zoom In Zoom Out Reset image size

While the feature still must be confirmed, we find three intriguing lines of evidence that point to the possibility of a single "dominant" planetary or potentially brown dwarf companion in the disk. The first is the very broad disk gap apparent in the mm-emission. With a width of over 70 au centered around approximately 50 au, a massive gas giant or low-mass brown dwarf would be necessary to clear such a gap. Based on the estimates presented in Section 4.2, its mass must be greater than 4 M_Jup up to 70 M_Jup. How such a massive object forms at this distance remains uncertain. Although there are cases of wide-orbit giant exoplanets (e.g., Marois et al. 2010), the relative occurrence is not significant, and they are as rare as wide-orbit brown Dwarfs (Nielsen et al. 2019).

In addition to the disk morphology, there are two lines of evidence indicating preferential structure inside the southern gap. The first is the presence of a CO "hot spot" at a position coinciding with the gap. Given that young planets are expected to be luminous and thereby heat their surroundings, this may be a sign of such localized heating (Wolf & D'Angelo 2005; Cleeves et al. 2015). The second line of evidence is the presence of marginal radio emission, VLA 1b, inside the southern gap. No such source is seen in the ALMA data, suggesting that the emission originates from ionized gas (see Section 3.4). It is common to see radio emission associated with strong winds and/or jets from protostellar objects (Hull et al. 2016; Anglada et al. 2018). There is a known correlation between radio luminosity and a YSO's bolometric luminosity. If VLA 1b is of the same kind, its radio flux would be consistent with a substellar object (Morata et al. 2015; Rodríguez et al. 2017; Ricci et al. 2017). If confirmed, VLA 1b would be a signpost that there is a substellar object significantly accreting gas from the disk. Alternatively, the radio emission could be produced by the a strong magnetosphere of a fast rotating (sub)stellar object. However, the radio luminosities observed in old brown dwarfs are 1−2 orders of magnitudes lower than the flux detected here (Berger et al. 2001; Kao et al. 2018). Yet, these nonthermal mechanism possibilities can be variable, and so repeat follow-up observations would likely be necessary in ascertaining the nature of VLA 1b. Near-infrared observations from scattered light produced by the stellar radiation are also being analyzed in order to constrain the nature of the companion (A. Zurlo et al., 2020, in preparation).