© The Authors 2023

1 Introduction

During the early prestellar core stage of star formation, the temperature is low (10 K) and the density is high, that is, higher than 10⁵ cm⁻³ in the central region (Caselli & Ceccarelli 2012). This cold, dense environment encourages the formation of the deuterated species H₂D⁺ (first detected by Stark et al. 1999), which is an exothermic process (Gerlich et al. 2002): $${\text{H}}_3^{+}+\text{HD}\rightleftharpoons {\text{H}}_2{\text{D}}^{+}+{\text{H}}_2+232\text{\hspace{0.17em}K},$$(1)

while also giving rise to the depletion of CO due to carbon-, oxygen- and nitrogen-bearing species being frozen-out onto dust grains (Caselli et al. 1999, 2002). The rate of formation of H₂D⁺ is slowed down through the preferential reaction of the ${\text{H}}_3^{+}$ ion with CO to form HCO⁺ instead of with HD to form H₂D⁺: $${\text{H}}_3^{+}+\text{CO}\to {\text{HCO}}^{+}+{\text{H}}_2,$$(2)

and therefore the absence of CO in these environments also favours the production of the deuterated species. Furthermore, we can go beyond this with the production of D₂H⁺ (first detected by Vastel et al. 2004) and ${\text{D}}_3^{+}$.

The protostellar phase is initiated because of cloud collapse as a result of gravity overcoming pressure. This generally occurs when the cloud mass reaches a limit known as the Jeans mass, although this value can change during the collapse process, leading to the formation of multiple cores. During the evolutionary process, the temperature increases up to approximately 100 K; CO is released back into the gas phase from the dust grains when the temperature surpasses the critical value of approximately 25 K. Therefore, as the embedded protostar evolves and more species are injected into the gas phase via sublimation from the grain surface, the destruction of H₂D⁺, D₂H⁺ and ${\text{D}}_3^{+}$ occurs through gas phase reactions with neutrals. In particular, reactions with CO and N₂ lead to an increase in DCO⁺ and N₂D⁺, respectively; all of these deuterated ions eventually recombine in the gas phase to form, alongside other species, D atoms, which lead to the formation of neutral deuterated species and an increase in the D/H ratio in the gas phase. Recombination also takes place on the icy surfaces of dust grains, allowing deuteration to also occur in this phase. Formaldehyde (H₂CO) is already present in the core by this stage and therefore, in theory, deuterium in its increased abundance will form deuterated formaldehyde.

Multiple observations have confirmed the presence of both HDCO and D₂CO in protostellar cores (e.g. Loinard et al. 2002; Parise et al. 2006; Persson et al. 2018; Sahu et al. 2018). The D/H ratio of protostellar cores varies depending on the species, but molecular D/H enhancements of up to 13 orders of magnitude compared to the cosmic value of 10⁻⁵ have been observed (Loinard et al. 2002; Parise et al. 2004; Ceccarelli et al. 2007); including D/H ratios of up to 30% in D₂CO (Loinard et al. 2002). More recently, Persson et al. (2018) found the formaldehyde D/H ratio in the low-mass protostar IRAS 16293–2422 B to be 0.03 for HDCO and 0.08 for D₂CO.

The exact deuteration pathway for formaldehyde is not conclusively known. Two scenarios have been presented because undeuterated H₂CO can be formed both in the gas phase (Parise et al. 2009; Yamamoto 2017) and on grain surfaces (Watanabe 2005; Roberts & Millar 2007; Bergman et al. 2011). Therefore, deuteration pathways in both phases have been explored. The grain surface pathway proposed by Rodgers & Charnley (2002) suggests that (deuterated) formaldehyde is formed in intermediate steps of the (deuterated) CH₃OH formation pathway involving (D-addition or) hydrogenation. Meanwhile, the gas phase pathway has been proposed to involve formaldehyde deuteration either directly from reaction with H₂D⁺ or via the deuteration of ${\text{CH}}_3^{+}$ (Roberts & Millar 2007; Roueff et al. 2007), which showed consistency with subsequent modelling predictions (Bergman et al. 2011). However, observational analyses (Taquet et al. 2012; Ceccarelli et al. 2014; Fontani et al. 2014; Persson et al. 2018; Manigand et al. 2020; Zahorecz et al. 2021), as well as chemical modelling (Rodgers & Charnley 2002), favour the grain surface deuteration pathway as the most likely for formaldehyde.

Despite all of this previous research, our knowledge of the fractionation processes in H₂CO is not yet clear, especially at small linear scales, because previous results were mostly obtained through single-dish observations and therefore represent the average values of regions known to have density and temperature variations at small scales. Interferometric observations are therefore required in order to probe deeper into protostellar objects and reveal the chemistry within. In particular, the target studied here, [BHB2007] 11, is highly understudied compared to other similar sources, such as IRAS 16293–2422 B (Ceccarelli et al. 2001; Bottinelli et al. 2004; Crimier et al. 2010; Persson et al. 2018). With this in mind, [BHB2007] 11 was chosen as one of 13 Class 0 and Class I sources to be observed with the Atacama Millimeter/Submillimeter Array (ALMA) within the context of the Large Program Fifty AU STudy of the chemistry in the disk/envelope system of Solar-like protostars (FAUST¹; PI: Satoshi Yamamoto; see also Codella et al. 2021), with particular setups chosen, in part, in order to observe multiple deuterated species. Additionally, observing the core with multiple configurations enables both large- and small-scale analyses of our source. These observations will provide a new perspective to the chemical analysis of a protostellar core. In addition, the present work represents the first deuteration study towards [BHB2007] 11.

The object known as [BHB2007] 11 is a Class 0/I protostellar core (Brooke et al. 2007; Forbrich et al. 2009; Sandell et al. 2021) located 163 ± 5 pc away (Dzib et al. 2018) in the Barnard 59 (B59) core. The source is among at least 20 low-mass young stellar objects (YSOs) forming a protocluster within B59. While the members of this protocluster have been shown to exhibit diversity in evolutionary stage (Onishi et al. 1999), the YSO [BHB2007] 11 is believed to be the youngest among the cluster at approximately 0.1–0.2 Myr old (Brooke et al. 2007; Riaz et al. 2009). This object has a bolometric temperature of 60–70 K, a luminosity of 2.2–4.5 L_⊙ and systemic velocity of 3.6 km s⁻¹ (Onishi et al. 1999; Brooke et al. 2007; Forbrich et al. 2009; Sandell et al. 2021). Previous observations revealed [BHB2007] 11 to be a protobinary system with two central objects (hereafter denoted A (α(2000)= 17h11m23.1058s, δ(2000) = −27°24′32.828″) and B (α(2000)= 17h11m23.1015s, δ(2000) = −27°24′33.987″), respectively, located approximately 28 au apart (in projection) with A to the north of B (Alves et al. 2019). Source B has been shown to exhibit preferential accretion from the circumbinary disk, which was traced by infalling gas within filaments (Alves et al. 2019), while Source A is the most massive of the two protostars. Preferential accretion onto the least massive companion of a protobinary system has been predicted by simulations (Bate & Bonnell 1997; Bate et al. 2002; Matsumoto et al. 2019).

In addition, recent observations revealed the presence of narrow structures from beyond the core down to disk scales in free-fall, which are named streamers (see Pineda et al. 2023 fora review), while observations of CO, C¹⁸O and H₂CO showed that this core features a symmetrical bipolar outflow. It has also been determined that this outflow is launched from outside the disk edge (at about 100 au; Alves et al. 2017), which is unusual for this type of object. Using H₂CO observations; Alves et al. (2017) also detected the existence of a centrifugal barrier marking a sharp change in kinematics within the spiral structure, which itself is revealed in continuum (Alves et al. 2017). Position–velocity (PV) diagrams towards [BHB2007] 11 have varying profiles depending on the species; while higher energy transitions (E_u ~ 70 K) of H₂CO showed a smooth Keplerian PV plot profile (Alves et al. 2019), CH₃OH and lower energy transitions of H₂CO (E_u ~ 20 K) showed a more fragmented profile consistent with multiple velocity components (Alves et al. 2017; Vastel et al. 2022). This CH₃OH emission has been detected in three identified positions within 30 au of A and B in the circumstellar region (in contrast to the more ubiquitously detected H₂CO) and is believed to originate from the sublimation of dust grains due to non-thermal processes, such as the existence of a shocked region. These shocks would most likely be caused by the interaction between quiescent gas and previously detected streamers connecting the circumbinary disk to the circumstellar region (Alves et al. 2019; Vastel et al. 2022).

The structure of this paper is set out as follows: Sect. 2 describes the observations, Sect. 3 describes the emission morphology, Sect. 4 discusses the radiative transfer analysis and Sect. 5 contains the conclusions.

2 Observations

Within the context of the Large Program FAUST, the Class 0/I protostellar core [BHB2007] 11 was observed in the ALMA large project 2018.1.01205.L between 2018 and 2020 with a spectral resolution of 0.14 MHz (equivalent to a velocity resolution of 0.2 km s⁻¹) in setup 1 and 0.12 MHz (equivalent to a velocity resolution of 0.2 km s⁻¹) in setup 2, along with an angular resolution of 0.25″ (equivalent to a spatial resolution of ~40 au) across the setups. This source was observed with the antennas configured in three different ways: two configurations with the main array (12 m) and one configuration with the Atacama Compact Array (ACA, 7 m antennas). In this paper, we present observations derived from the combination of these 12 m and 7 m configurations using the ‘mosaic’ gridder in tclean within CASA. Our combination process involves identifying line-free continuum channels common to all datasets and configurations by hand, ad subsequently using self-calibration of the resulting continuum to align both phases and amplitudes to take out position and amplitude calibration offsets between execution blocks. During this process, we take care to make sure that the resulting continuum model is as complete as possible in order to avoid introducing any scaling offsets that can suppress the overall amplitude scale. The final robust weighting used for the continuum is one chosen to deliver 50 au resolution (as required for the FAUST project) but natural weighting is used during the self-calibration to ensure all the emission is recovered by the sky model during the imaging, resulting in a beamsize of 0.43″ × 0.41″. The FAUST program provides observations in Bands 3 and 6 divided into three spectral setups, for a total of 13 spectral windows (spws). In the present paper, we analyse lines from Band 6 only, because the frequency range covers the formaldehyde isotopologues that we are targeting. Therefore, in Table 1 we list only the frequency ranges of the windows included in these setups. Table 1 also contains useful observational information, such as the angular resolution and the maximum recoverable scale achieved in each configuration and setup, as well as water vapour conditions. The system temperature ranged from 50-160 K. For 12 m observations in setups 1 and 2, the quasars J1427-4206 and J2056-4714 were used as bandpass calibrators, while J1700-2610 was used to calibrate phase and flux. The 1.3 mm continuum was observed with spectral resolutions of 0.98 MHz (setup 1) and 1.13 MHz (setup 2).

Table 2 shows the spectroscopic parameters of all of the formaldehyde transitions, including singly (HDCO) and doubly (D₂CO) deuterated transitions, which were identified within the FAUST dataset using CASSIS² software.

The data were reduced, calibrated and cleaned using the Common Astronomy Software Application (CASA³, bf CJC: version 5.6.1-8) using a pipeline developed by the FAUST data reduction team (see Vastel et al. 2022 for more details), with a Briggs robust uv-weighting of 0.5.

Fig. 1 Comparing continuum and line emission. Left panel: H2CO 3(0,3)−2(0,2) moment 0 map integrated between −5 and 15 km s−1. This velocity range was chosen as this represents the full range of emission detected at or above 4σ. Right panel: 216 GHz continuum emission (colour scale) taken from the FAUST dataset with H2CO 3(0,3)−2(0,2) overlaid in pink contours. In both images, contours start at 4σ and end at 44σ in 8σ increments (σ = 17.6 mJy beam−1 km s−1), the positions of protostars A and B are represented by the black stars and the ALMA synthesised beam is shown in the lower left corner.

3 Emission morphology

The clean maps were analysed with the GILDAS⁴ software. The line channels were identified and used to create moment 0 maps for each of the transitions. In order to investigate both large-and small-scale structure, we combined observations performed using both 12 m and 7 m ALMA configurations in all moment 0 maps shown in this section. For the analysis of all lines in our sample, the 4σ contour (the minimum to define a detection) of the map of the HDCO transition at 259 GHz was used to define a polygon for extraction. This transition was chosen to estimate the D/H ratio because it is detected at a signal-to-noise ratio (S/N) of 8σ and because the HDCO transition at 246 GHz shows possible signs of blending with the CH₃OCHO (19-18) transition.

The spectra of the four transitions were extracted in flux density units (S_v) within the aforementioned polygon; these were converted to synthesised beam temperature units (K) in the GILDAS package CLASS using conversion factors calculated using Eq. (3), which was taken from the ALMA handbook (Warmels et al. 2018): $$\left(\frac{T}{1\text{\hspace{0.17em}K}}\right)=\left(\frac{S_v}{1\text{\hspace{0.17em}Jy}}\right)\left[13.6{\left(\frac{300\text{\hspace{0.17em}GHz}}{v}\right)}^2\left(\frac{1″}{{\theta }_{\text{max}}}\right)\left(\frac{1″}{{\theta }_{\text{min}}}\right)\right],$$(3)

where θ_min and θ_max represent the half power beam widths along the minor and major axes, respectively.

3.1 H₂ CO 3(0,3)−2(0,2)

The emission morphology of H₂CO 3(0,3)−2(0,2) can be seen in the moment 0 map shown in the left panel of Fig. 1. Kinematic analyses performed in previous studies (Alves et al. 2017) and also later in this paper show that this is a disk+outflow system.

The axis of the extension seems to be aligned with that seen in the continuum, which has a position angle of 138° (Alves et al. 2017); this is illustrated in the right panel of Fig. 1, which overlays the H₂CO 3(0,3)−2(0,2) emission morphology onto the 1.3 mm dust emission, which are clearly coincident.

Figure 2 shows the centroid velocity maps – obtained from a pixel-by-pixel Gaussian analysis using the pyspeckit package (Ginsburg et al. 2022) – overlaid with contours from 1.3 mm emission from the FAUST data. For this Gaussian fitting, we include pixels with velocities of between 3.5 and 6 km s⁻¹ (red-shifted lobe; left panel of Fig. 2) and between 0 and 3.5 km s⁻¹ (blueshifted lobe; right panel of Fig. 2) and a minimum intensity of 5σ. The redshifted and blueshifted high-velocity extensions of the disk emission are clearly visible to the NW and SE, along with lower-velocity extensions seen closer to the E–W axis. The high-velocity emission coincides with the morphology of the continuum emission and therefore appears to be originating from the disk, while the lower-velocity emission is believed to be tracing dynamic material extending beyond the circumbinary disk, potentially originating from incoming streamer material or an outgoing disk wind. The kinematics of the formaldehyde emission suggest the presence of a second outflow (see Appendix A for details). A threshold of 5σ was chosen for the centroid velocity maps because, as explored in the Appendix, this is the optimal S/N to model dynamic emission without dilution from surrounding emission. To obtain the maps in Fig. 2, the velocity limits for integration were chosen as 0 and 6 km s⁻¹ because the extended structure is not seen outside these velocities and therefore this velocity range is optimal for detecting and discerning both compact (disk) and extended (incoming streamer or outgoing disk wind) emission simultaneously. It should be noted that, despite observing the same transition, the higher-spatial-resolution H₂CO observations published in Alves et al. (2017) only appear to trace the disk emission aligned with the continuum emission and not the circumbinary large-scale flow. The accretion filaments surrounding the binary system itself are a distinct feature at a much smaller scale. The reason for this is that our shortest baseline and overall uv coverage are more sensitive to extended emission than those in Alves et al. (2017). Notably, as illustrated in Fig. 2, the velocity structure does not show alignment with the outflow seen in CO 2-1 (Alves et al. 2017) either in terms of morphology or velocity. The interferometric nature of these observations does lead to the possibility that some low-velocity CO emission may be filtered out; however, the region affected by this is spatially much more extended than the region where deuterated species are detected (as described in subsequent sections) and therefore the main conclusions of this work are unaffected. The features appear almost perpendicular to one another and at differing velocities, with the feature seen in H₂CO appearing closer to the systemic velocity of 3.6 km s⁻¹ while the feature traced by CO peaks at more extreme redshifted and blueshifted velocities. However, it should be noted that the beamsize of the observations presented in Alves et al. (2017) is 0.20″ × 0.25″, a mismatch with the present FAUST data, which has a beamsize of 0.43″ × 0.41″.

As can be seen in Fig. 3, the H₂CO spectrum shows three clear emission components at ~−5 to 0 km s⁻¹, ~1 to 8 km s⁻¹ and ~ 8 to 15 km s⁻¹. In order to support this finding and to show the origin of each component, Fig. 3 also shows an individual moment 0 map for each of the three emission components (integrated between −5 and 0 km s⁻¹, 1 and 8 km s⁻¹ and between 8 and 15 km s⁻¹). The three emission components are broadly consistent with those seen in CH₃OH towards this object, which were seen at approximately −2, 2.8 and 9.9 km s⁻¹ (Vastel et al. 2022). Evidently, as seen in CH₃OH, the two H₂CO emission components at more extreme velocities are originating from the rotating disk as emission from material moving towards Source B (previously observed at high velocities by Alves et al. 2019), while the central component originates from Source A and the surrounding envelope. In Sect. 3.5, we explore and compare the origins of CH₃OH and H₂CO emission in this source in more detail. Additionally, H₂CO exhibits an absorption feature (~2 to 5 km s⁻¹) due to absorption in the cold envelope of emission originating from warm inner regions, which, according to our pixel-by-pixel analysis, is seen in a circular region that extends to approximately 500 au from the centre and is centred on the dust peak. We note that in order to model this absorption feature, a continuum value of 10.3 K has been added to the spectrum shown in the upper left panel of Fig. 3. This value was obtained by measuring the flux of the 216 GHz continuum emission.

Fig. 2 Colour scale: H2CO 3(0,3)−2(0,2) centroid velocity map from pixel-by-pixel Gaussian analysis showing pixels with intensity above 5σ in the redshifted (left panel; integrated between 3.5 and 6 km s−1) and blueshifted (right panel; integrated between 0 and 3.5 km s−1) regimes; the systemic velocity is ~3.6 km s−1. Blue contours: 1.3 mm continuum contours starting at 4σ and ending at 300σ in intervals of 15σ (σ = 0.19 mJy beam−1) indicating disk emission. Red (left panel) and cyan (right panel) contours: Redshifted and blueshifted (respectively) CO 2−1 emission from Alves et al. (2017) detailing outflow emission seen at 8.67 km s−1 and −0.78 km s−1, respectively, representing the channels with the most extended outflow CO 2−1 emission. Contours are at 4, 6, 8, 10, 15, 20, 30, 40, 50, 60, 70, 80 and 90 times σ – which is 2.8 mJy beam−1 – in order to match the CO 2−1 observations, which have a beamsize of 0.2″ × 0.25″. The ALMA synthesised beam for the H2CO observations, with a beamsize of 0.43″ × 0.41″ is shown in the lower left corner.

3.2 HDCO 4(2,2)−3(2,1)

The emission morphology for HDCO 4(2,2)−3(2,1) at 259 GHz is shown in Fig. 4. The compact disk emission is clearly seen. However, in contrast to the H₂CO emission in Fig. 1, the full extended structure is not visible in this map. Figure 4 shows a moment 1 map integrated between the same velocity range as Fig. 2, confirming that, of the two velocity structures seen in H₂CO (see Fig. 2), only the more compact velocity structure representing disk emission is fully present in HDCO. However, there is some asymmetry present in the morphology of this transition that would not be expected from the disk. Comparison with Fig. 2 reveals the possibility that this is arising at the base of the dynamic feature in the redshifted regime.

As can be seen in Fig. 5, we extracted the HDCO 4(2,2)−3(2,1) transition, although this spectral window (spw) also contains a much stronger transition of H¹³CN (3−2). However, this transition (rest frequency 259.0118 GHz) has a central peak that is well separated from the HDCO transition (rest frequency 259.0349 GHz) and therefore does not interfere at all with our analysis. The HDCO transition shows three emission components between the velocities of ~−5 and 12 km s⁻¹, which likely originate from the disk. Indeed, as shown by the contours of Fig. 4, parts of the HDCO emission, notably those contributing to the aforementioned asymmetric morphology, fall outside the 4σ detection threshold denoted by the cyan contour and this supports the notion that only disk emission is included in the extracted spectrum seen in Fig. 5. More specifically, comparing Fig. 3 – which shows the morphology of each velocity component relative to Sources A and B – to Fig. 5, the HDCO component between ~−5 and 0 km s⁻¹ and the HDCO component with a lower limit of 6 km s⁻¹ are likely originating from Source B, while the HDCO component between ~1 and 6 km s⁻¹ likely originates from Source A, similarly to the H₂CO emission. However, the elevated baseline visible in this spw complicates the process of defining the upper limit on the velocity range of the emission. The spectrum is flat between ~12 and 20 km s⁻¹ but is consistently at an intensity of approximately 0.1 K. The flatness of the spectrum seems to be too consistent for any lines to be present and therefore we postulated that this may be an effect of the wing emission of the nearby strong H¹³CN transition. In order to check this, the spw was recleaned using the same method (the FAUST pipeline) but with a reduced velocity window excluding the H¹³CN transition. In this way, we can see whether this elevated baseline is still present when the H¹³ CN transition is disregarded, by only considering velocities ~15 km s⁻¹ either side of the centre. When comparing the narrower spectrum to the spectrum covering the full ~50 km s⁻¹ velocity range (shown in Fig. 5), the baseline around 12 to 20 km s⁻¹ is still at a consistent intensity of 0.1 K. Additionally, not only does the frequency range of continuum subtraction cover this spw but also there was no evidence of issues of bandpass calibration, which could have caused such an artefact to appear. We therefore conclude that the cause of the elevated baseline is most likely line emission of HDCO. However, Fig. 4 is integrated between a velocity range that excludes the affected component so that it is certain that only confirmed HDCO 4(2,2)−3(2,1) emission is included, while the spectrum shown in Fig. 5 has been corrected for this raised baseline.

Fig. 3 Comparing the multiple velocity components of the H2CO 3(0,3)−2(0,2) emission. Upper left: H2CO 3(0,3)−2(0,2) spectrum (black) with velocity resolution 0.17 km s−1 covering the full velocity range of −5 to 15 km s−1, extracted from the region where HDCO 4(2,2)−3(2,1) emission is above 4σ, overlaid with the best fit (pink) obtained using four Gaussian fits with FWHM, υLSR and Tex constrained using values obtained from fitting the HDCO transitions. A continuum value of 10.3 K has been added to this spectrum in order to model the absorption feature located at the systemic velocity of 3.6 km s−1. Upper right, lower left, lower right: H2CO 3(0,3)−2(0,2) moment 0 map integrated between −5 and 0 km s−1, between 0 and 8 km s−1 and between 8 and 15 km s−1, respectively, in order to individually show the velocity components of emission shown in the upper left panel. Contours begin at 4σ and are in increments of 8σ (σ = 17.6 mJy beam−1 km s−1). The pink contour illustrates the area within which the spectrum in the upper left panel was extracted, defined as the 4σ contour of the HDCO 4(2,2)−3(2,1) moment 0 map (see Fig. 4). The ALMA synthesised beam is in the lower left corner of all three maps.

3.3 HDCO 4(1,4)−3(1,3)

As illustrated in Fig. 6, the HDCO 4(1,4)−3(1,3) spectrum is also consistent with three emission components at ~−5 to 0 km s⁻¹, ~1 to 6 km s⁻¹ and a third component with v_LSR ≥ 6 km s⁻¹ that is affected by blending. The line-identification software in CASSIS was used to check the transitions that could potentially cause this blending; based on the Einstein coefficients and E_u (upper energy level, using a maximum value of 150 K), the most likely culprit is CH₃OCHO 19(4,15)−18(4,14) A, as illustrated in Fig. 7. With the HDCO components showing similar velocity values to those seen in H₂CO (see Fig. 3), CO (Alves et al. 2017) and CH₃OH (Vastel et al. 2022), we conclude that, similar to these species, the high-velocity HDCO emission components are most likely originating from Source B while the central component likely arises from Source A.

Figure 8 shows the moment 0 map of the HDCO 4(1,4)−3(1,3) emission, along with the area used to extract the spectra seen in Figs. 3, 5 and 6. The emission is compact and – unlike the HDCO 4(2,2)−3(2,1) emission – symmetrical, suggesting that this emission is likely only originating from the disk.

Fig. 4 Emission maps for HDCO 4(2,2)−3(2,1). Left panel: HDCO 4(2,2)−3(2,1) moment 0 map integrated between −1 and 7 km s−1 that peaks at the position of the dust emission. Observations performed using both 12 m and 7 m ALMA configurations are combined in this map in order to investigate both large- and small-scale structure. Contours start at 4σ and end at 10σ in 3σ increments (σ = 4.1 mJy beam−1 km s−1). The 4σ contour is highlighted in cyan as this is the contour that served as the extraction polygon for all spectra shown in this work. Right panel: Moment 1 map of HDCO 4(2,2)−3(2,1) integrated between −1 and 7 km s−1, in order to directly compare with the moment 0 map and to show the velocity structure of the disk. In both maps, the locations of the protostars A and B are shown by black stars. The ALMA synthesised beam is in the lower left corner.

Fig. 5 HDCO spectrum (black) with a velocity resolution of 0.14 km s−1 for the 4(2,2)−3(2,1) transition overlaid with the best fit (pink) from our simultaneous radiative transfer analysis, extracted within the 4σ contour shown in Fig. 4. The analysis was run twice, assuming upper limits for the size of the emitting region of 1.0″ and 0.3″ for Sources A and B, respectively and lower limits of 0.5″ and 0.15″ for Sources A and B, respectively; the figure shown represents the latter. Dashed lines denote the previously identified velocity components at −2, 2.8 and 9.9 km s−1.

Fig. 6 HDCO spectrum (black) with a velocity resolution of 1.18 km s−1 for the 4(1,4)−3(1,3) transition overlaid with the best fit (pink) from our simultaneous radiative transfer analysis, extracted within the 4σ contour shown in Fig. 4. The analysis was run twice, assuming upper limits for the size of the emitting region of 1.0″ and 0.3″ for Sources A and B, respectively and lower limits of 0.5″ and 0.15″ for Sources A and B, respectively; the figure shown represents the latter. We note that the third component is not included in the fit because of the blending with a transition of CH3OCHO 19(4,15)−18(4,14) A, the emission components of which are clearly seen above 6 km s−1 and highlighted explicitly in Fig. 7. Dashed lines denote the previously identified velocity components at −2, 2.8 and 9.9 km s−1.

3.4 D₂CO 4(0,4)−3(0,3)

As can be seen in Fig. 9, a single peak is seen in D₂CO 4(0,4)−4(0,3), with emission seen between approximately 1 and 4 km s⁻¹. As this peak is visibly weak, we use Eq. (4) to calculate the integrated intensity (W; K km s⁻¹) value for the 3σ upper limit and compare this to the strength of our observations. $$W=\left(3\times \text{rms\hspace{0.17em}\hspace{0.17em}}\times FWHM\right)\sqrt{{\left(2\times cal\right)}^2+\left(2\times \frac{\text{Δv}}{FWHM}\right).}$$(4)

In this equation, cal represents the calibration error (taken in this case to be 10%), Δυ is the velocity resolution (which is 0.2 km s⁻¹) and the rms is 0.15 K. Therefore, using a measured FWHM line width obtained using Gaussian fitting within CASSIS of 1.4 km s⁻¹, we obtain a 3σ upper limit value of 0.36 K km s⁻¹. The integrated intensity of our D₂CO line is 0.826 K km s⁻¹, therefore yielding a ratio of 2.3 compared to the 3σ threshold. Thus, our D₂CO observations have a S/N ratio of approximately 7σ in integrated intensity.

Fig. 7 Extracted spectrum centred on the rest frequency of the CH3OCHO 19(4,15)−18(4,14) A transition at 246.9147 GHz with dashed lines showing the predicted locations of the velocity components assuming similar VLSR to methanol components (−2, 2.8 and 9.9 km s−1). The locations of the same velocity components (−2, 2.8 and 9.9 km s−1) in HDCO 4(1,4)−3(1,3) are shown by the pink dashed lines. The −2 km s−1 component of CH3OCHO 19(4,15)−18(4,14) A is blended with the 9.9 km s−1 component of HDCO 4(1,4)−3(1,3).

Fig. 8 Moment 0 map of HDCO 4(1,4)−3(1,3) integrated between -5 and 8 km s−1. This velocity range was chosen to include the two unblended HDCO emission components shown in Fig. 6. Contours start at 4σ and end at 16σ in 4σ increments (σ = 8.5 mJy beam−1 km s−1). The pink contour illustrates the area within which the spectrum in Fig. 6 was extracted, defined as the 4σ contour shown in Fig. 4. The locations of protostars A and B are shown by the black stars. The ALMA synthesised beam is shown in the lower left corner.

Fig. 9 D2CO 4(0,4)−3(0,3) spectrum (black) with a velocity resolution of 0.16 km s−1 extracted from the 4σ contour shown in Fig. 4 and smoothed three times with the Hanning algorithm overlaid with best fit (pink) obtained using single-component Gaussian fitting with FWHM, VLSR and Tex constrained using values obtained from fitting the HDCO transitions. The analysis was run twice, assuming upper limit values for the size of the emitting region of 0.3″ (covering source B emission) and 1.0″ (covering source A emission); the figure shown represents the latter. The rms level (0.15 K) is shown in green.

3.5 Emission origin

As introduced in Sect. 1 (see also Rodgers & Charnley 2002), the grain surface formation pathways for methanol and formaldehyde are intrinsically linked in that both form through hydrogenation reactions on the surface itself: $$\text{CO\hspace{0.17em}}+\text{\hspace{0.17em}H}\to \text{HCO},$$(5) $$\text{HCO\hspace{0.17em}}+\text{\hspace{0.17em}H}\to {\text{H}}_2\text{CO},$$(6) $${\text{H}}_2\text{CO\hspace{0.17em}}+\text{\hspace{0.17em}H}\to {\text{CH}}_3\text{O},$$(7) $${\text{CH}}_3+\text{\hspace{0.17em}H}\to {\text{CH}}_3\text{OH}.$$(8)

During any of these stages, it is possible for a D-addition reaction to occur rather than H-addition, creating a grain surface deuteration pathway. It is widely understood (Watanabe & Kouchi 2002; Rimola et al. 2014) that these grain surface formation processes occur at low temperatures during the cold prestellar phase during which the atoms and molecules are frozen in the icy mantles of dust grains. This is followed by subsequent injection of the particles into the gas phase due to the increase in temperature associated with the protostellar phase. More specifically, this temperature rise can be the result of thermal (e.g. the presence of hot corinos – Taquet et al. 2016) or non-thermal (e.g. shocks caused by interaction of material with quiescent gas – Jiménez-Serra et al. 2008) processes. Previous analyses of methanol (Vastel et al. 2022) and formaldehyde (Alves et al. 2019) emission towards [BHB2007] 11 showed that the dust is optically thick towards the protobinary system and that the emission of both formaldehyde and methanol shows an increasing velocity gradient towards protostar B. This evidence strongly favours the scenario in which the injection of formaldehyde and methanol into the gas phase occurred due to a shock, as the observations appear to rule out the suggestion that A and/or B are hot corinos; indeed; Vastel et al. (2022) found that the CH₃OH emission peak appears spatially shifted with respect to both sources. In addition, dynamic material, which may be involved in causing such a shock, has previously been observed in this source (Alves et al. 2019) and is seen in the present work. Figure 10, which overlaps our four spectra on a single plot, also shows dashed lines denoting the three velocity components first identified in CH₃OH emission towards [BHB2007] 11 as described in Vastel et al. (2022) at −2, 2.8 and 9.9 km s⁻¹. Notably, this previous research suggests that the −2 and 9.9 km s⁻¹ components are originating from higher-velocity warm gas streaming towards Source B (see also Alves et al. 2019, who resolved these components for the first time), while the 2.8 km s⁻¹ CH₃OH component is believed to originate from Source A. Figure 10 also shows multiple velocity components in H₂CO and HDCO, which show consistency with the three components first identified in CH₃OH emission by Vastel et al. (2022). Meanwhile, D₂CO shows a single peak with a velocity shift compared to the previously identified velocity components.

Comparing our H₂CO observations with the high-resolution CH₃OH observations presented in Vastel et al. (2022), clear differences are apparent in the morphology of the emission of these species. Notably, it is very clear that CH₃OH is very compact compared to the more extended H₂CO seen in the present observations. This is despite the fact that both species are injected into the gas phase from the frozen surfaces of dust grains due to the rising temperature of the system under similar conditions. The T_ex values obtained from non-LTE analysis of the CH₃ OH emission ranged from 110 to 130 K, suggesting that this is tracing hot gas in this object (Vastel et al. 2022). Infall motions suggest the presence of a shocked region, which likely enhances the production of CH₃OH on small scales, while the low population of the upper states of the observed CH₃OH lines makes this species harder to detect in the cold envelope region. All of these factors may explain why CH₃OH and H₂CO appear to be tracing different regions of [BHB2007] 11, although it is also possible to explain the observed morphology differences if H₂CO is also being produced in the gas phase.

Figure 11 shows the moment 1 map for H₂CO emission in colour scale overlaid with five pixel positions from which we extract individual spectra. As previously discussed, there appears to be separate velocity components along two different axes in this moment 1 map, with the more east–west aligned components showing lower velocities than the components aligned with the continuum emission extension. Four of the pixels shown represent the blueshifted and redshifted velocity components along both axes, while the fifth represents the central position and the peak of the emission. Each pixel outside the centre was chosen as the most intense within four velocity ranges: 0−2.5 km s⁻¹, 2.5−3.5 km s⁻¹, 3.5–4 km s⁻¹ and 4–6 km s⁻¹. The H₂CO, HDCO (both transitions) and D₂CO spectra extracted from these five pixels are shown in Fig. 12. It is clear from these spectra that, while the H₂CO emission is extended, the HDCO emission is compact and not seen in any of the spectra outside that extracted from the central pixel; indeed, HDCO is only seen within ∼300 au of the centre. From this, we postulate that HDCO likely originates from the shocked regions first identified in Vastel et al. (2022) located on either side of Source B and also near Source A. Moreover, the absorption feature is only seen in the H₂CO spectrum extracted from the central pixel, further illustrating that this feature is not seen outside ∼500 au. Meanwhile, the low S/N of the D₂CO spectra in these particular pixels prevents any conclusions from being drawn.

Fig. 10 Comparison of H2CO 3(0,3)–2(0,2), D2CO 4(0,4)–3(0,3), HDCO 4(1,4)–3(1,3) and HDCO 4(2,2)–3(2,1) lines with the baselines subtracted. The dashed lines represent velocities of –2, 2.8 and 9.9 km s−1. The CH3OCHO 19(4,15)–18(4,14) A transition (246.9147 GHz) that is near HDCO 4(1,4)–3(1,3) is shown in red, while the H13CN (3−2) transition (259.0118 GHz) in the same spw as HDCO 4(2,2)–3(2,1) is denoted in green.

Fig. 11 H2CO moment 1 map in colour scale with five pixels overlaid in black and white (for clarity in this case) triangles, chosen as the pixel at the location of the dust peak (A), along with the most intense pixels for each velocity range (B: 0−2 km s−1, C: 2.5−3.5 km s−1, D: 3.5−4 km s−1 and E: 5−6 km s−1). The ellipse in the bottom-left corner represents the ALMA synthesised beam.

4 Radiative transfer analysis

The molecular column densities were derived assuming local thermodynamic equilibrium (LTE) conditions. Therefore, we can use equations from the Formalism for CASSIS Software, specifically within the ‘Documentation and Data’ section. In particular, we use Eq. (9) for the calculation of the total column density (N_TOT): $$N_{\text{TOT}}=\frac{8\pi v^3\left(\text{Δ}V\right){\tau }_{0}Q\left(T_{\text{ex}}\right)\sqrt{\pi }}{2\sqrt{\ln 2}}\frac{e^{\frac{E_{\text{u}}}{T_{\text{ex}}}}}{c^3{A}_{\text{ul}}{g}_{\text{u}}\left(e^{\frac{hv}{k{T}_{\text{ex}}}}-1\right)},$$(9)

where c, h, and k represent the speed of light, the Boltzmann constant and the Planck constant, respectively, in cgs units, while τ₀ is the opacity at the centre of the line. We take the background temperature (T_bg) as that of the cosmic microwave background (CMB; i.e. 2.73 K). We obtain the FWHM (∆V in cm s⁻¹ in Eq. (9)) through Gaussian fitting of the spectra. Frequency (ν in Hz), partition function ((QT_ex)), Einstein coefficient (A_ul), energy of the upper level (E_u in K) and upper level population (ɡ_u) are taken from the Cologne Database for Molecular Spectroscopy (CDMS; Müller et al. 2001).

Fig. 12 Spectra extracted from the five pixels (labelled A, B, C, D and E) indicated in Fig. 11. Clockwise from upper left: H2CO (218 GHz), HDCO (259 GHz), D2CO (231 GHz) and HDCO (246 GHz).

4.1 HDCO

AS there are two HDCO transitions in this sample, simultaneous Levenberg-Marquardt fitting within CASSIS of both transitions was employed in order to obtain the initial Gaussian fit of the extracted spectra, yielding best-fit values for V_LSR, FWHM and intensity for each visible component. However, as previously mentioned, one of the components in HDCO 4(1,4)−3(1,3) shows blending with one of the components of CH₃OCHO 19(4,15)−18(4,14) A, as predicted in Fig. 7. For this reason, the two unblended components were fitted simultaneously with the corresponding components in HDCO 4(2,2)−3(2,1) in order to constrain FWHM and V_LSR; the FWHM and V_LSR values for the blended component were thus fixed by fitting the 4(2,2)−3(2,1) transition alone; we apply the same FWHM and V_LSR values to the same component of both transitions in order to calculate N_TOT for this component. We also use the assumption that the HDCO emission is of a similar origin to that of CH₃OH to constrain our FWHM and V_LSR values. In order to check this assumption, we ran a 2D deconvolution of all three emission components in both HDCO transitions (with the exception of the blended component of the HDCO 4(1,4)−3(1,3) transition) applying the CASA imfit command to the peak emission channel in each case. The results of this deconvolution returned image-component size estimates (major axis FWHM x minor axis FWHM) of 0.34″ × 0.27″ for the component at approximately −1 km s⁻¹, 1.37″ × 0.68″ for the component at approximately 3 km s⁻¹ and 0.30″ × 0.20″ for the component at approximately 9 km s⁻¹. These findings agree with the 2D deconvolution presented in Vastel et al. (2022) in that the component associated with Source A is slightly more extended than those associated with Source B. Using these results, along with the consistency between emission components of CH₃OH and H₂CO, we therefore ran two LTE radiative transfer analyses, one assuming a lower limit emitting size of 0.15″ (for the components associated with Source B) and 0.50″ (for the Source A component) and the other assuming an upper limit emission area of 0.30″ (for the components associated with Source B) and 1.00″ (for the Source A component). Therefore, all N_TOT values for HDCO are corrected for beam dilution accordingly.

Figures 5 and 6 in Sect. 3 show the Gaussian fitting for the HDCO transitions. Within CASSIS, we used constraints provided by our initial Gaussian fitting in order to obtain best-fit values for FWHM, V_LSR, T_ex and N_TOT of our two unblended components at ∼−1 and ~ 3 km s⁻¹ using Monte Carlo Markovchain (MCMC) fitting, assuming a 10% calibration error. As we assume a common origin of the disk emission, V_LSR, FWHM and T_ex are also applicable to our H₂CO transition, enabling us to also obtain N_TOT values for this, which are reported in Table 3. As we are assuming LTE, all excitation temperatures for one molecule are the same; therefore, these best-fit parameters were then used to obtain an initial constraint for excitation temperature (T_ex) for the HDCO component at ∼−1 km s⁻¹, which was found to be approximately ${90}_{-8}^{+3}\text{\hspace{0.17em}K}$.

Following our analysis, we also obtained an upper-limit column density for absorption in HDCO of 6.6 × 10¹¹ cm⁻². We obtained this by noting the maximum T_ex for which absorption is observed after adding a continuum value of 10.3 K − obtained by measuring the peak intensity of the continuum observations of our dataset (see Fig. 1) – and running an LTE analysis at this temperature.

4.2 H₂CO

As we are assuming LTE conditions, we can also apply our best-fit T_ex value obtained via our HDCO analysis to fit our single H₂CO transition. Therefore, using a fixed T_ex of 90 K and our best-fit V_LSR and FWHM values from HDCO, we constrained Gaussian fits for each of the three emission components visible in the H₂CO spectrum. In addition to this, the absorption feature also required its own Gaussian fitting. To achieve this, we added a value for the continuum of 10.3 K from our observations and assumed an excitation temperature of 13 K to simulate the temperature of the cold envelope in which the absorption was occurring. This excitation temperature is the upper limit obtained for absorption in HDCO and, as such, the obtained column density for H₂CO in the absorption component is also an upper limit. We use this temperature in order to facilitate calculation of upper limits for the column density and D/H ratios in this component, which are shown in Table 3. We also note that the V_LSR value associated with the absorption component corresponds with the systemic velocity of the source (3.6 km s⁻¹). The resultant fitting of H₂CO can be seen in Fig. 3 in Sect. 3.1.

4.3 D₂CO

Levenberg-Marquardt fitting yields a FWHM of 2.8±0.6 km s⁻¹ and a V_LSR of 2.1±0.3 km s⁻¹ for the single detected D₂CO component. This V_LSR value is not consistent with those obtained for H₂CO and HDCO and therefore it is unclear which of the three identified emission components this D₂CO emission corresponds to. Furthermore, due to the weakness of the emission, the 2D deconvolution failed to converge, leading us to assume that the D₂CO emission has a similar emitting size to that of HDCO. Therefore, we used LTE analysis with the aforementioned V_LSR and FWHM values, an assumed excitation temperature of 90 K and an upper limit emitting size of both 0.3″ (covering Source B emission) and 1.0″ (covering Source A emission) to obtain upper limit values of NTOT for the single D₂CO component. Thus, once again, we take the effects of beam dilution into account.

4.4 D/H ratio derivation

Table 3 shows the best-fit parameters, which were constrained using Levenberg-Marquardt fitting and MCMC modelling, for FWHM, V_LSR and intensity, along with the resultant values for N_TOT, all with their relative uncertainties, for all transitions in our sample. Notably from our results, the FWHM of the components seen in the formaldehyde spectra show similarities with those of methanol (Vastel et al. 2022), again providing support for the theory that these have a common origin (refer to Sect. 3.5). However, it should be noted that the number of higher S/N CH₃OH transitions analysed means the FWHM values for these emission components are better constrained than for the formaldehyde emission presented in this work.

Due to the presence of two H atoms in H₂CO, either of which may be deuterated to form HDCO, we calculate the D/H ratios using the following equation: $$\frac{\text{D}}{\text{H}}=\frac{N\left(\text{HDCO}\right)}{N\left({\text{H}}_2\text{CO}\right)}\times 0.5,$$(10)

where N(HDCO) and N(H₂CO) represent the column densities (N_TOT) of HDCO and H₂CO, respectively. Using the N (HDCO) results shown in Table 3, we calculate D/H ratios of ${0.02}_{-0.01}^{+0.02}$ for the component at approximately −1 km s⁻¹, ${0.015}_{-0.004}^{+0.008}$ for the component at approximately 3 km s⁻¹ and ${0.01}_{-0.01}^{+0.03}$ for the component at approximately 9 km s⁻¹ (the latter computed using only the unblended transition). This gives an average D/H ratio for HDCO in this source of ${0.02}_{-0.01}^{+0.02}$ calculated using Eq. (10). However, for D₂CO, we cannot conclude with any certainty on where the emission is originating from or which of the multiple components detected in H₂CO and HDCO our single D₂CO emission component corresponds to. Therefore, while we can use our results to obtain upper limit values for N_TOT, we cannot constrain the D/H ratio for D₂CO with any kind of certainty.

For comparison, Persson et al. (2018) found a D/H ratio of 0.03±0.01 in HDCO towards the protobinary IRAS 16293−2422 B, within the uncertainty of the value found in the present work towards [BHB2007] 11. However, it should be noted that IRAS 16293−2422 B harbours a hot corino as interstellar complex organic molecules (iCOMs) have been observed towards the centre of the hot core (Bottinelli et al. 2004; Pineda et al. 2012; Jørgensen et al. 2016). It has been shown that emission of methanol, the simplest iCOM, originating from the protostars of [BHB2007] 11 is optically thick and is therefore not centred on the hot core as one would expect for a hot corino (Vastel et al. 2022). Meanwhile, Bianchi et al. (2017) found D/H for HDCO towards the Class I source SVS13-A to be 0.086, which is also higher than the average value obtained for [BHB2007] 11. Similarly, Zahorecz et al. (2021) obtained D/H ratios of 0.05 in HDCO towards a sample of high-mass star-forming regions, which is again higher than the average value obtained for [BHB2007] 11.

In summary, the D/H ratio in HDCO obtained in this work is consistently low compared with values obtained from other deuteration studies, although the value seen is within the error of that found towards the most similar source compared here: another protobinary, namely IRAS 16293 2422 B. Nevertheless, it is worth noting that our D/H calculation assumes that the H₂CO emission is uniform across our extraction region where HDCO is emitting. However, as previously discussed, the H₂CO emission shows similarity to that of CH₃OH, which would not show a uniform spectrum as compact emission regions have been detected (see Vastel et al. 2022). This means that although compact H₂CO emission regions are also likely to exist, our current resolution is insufficient to confirm this through direct observation. Therefore, to account for this, the D/H ratios for [BHB2007] 11 presented in this work are to be taken as upper limit values.

5 Conclusions

We present detections of formaldehyde and its deuterated isotopologues towards the protobinary system [BHB2007] 11, which we observed using the ALMA interferometer within the context of the FAUST Large Program. From this data, we present column density measurements of all species obtained using a radiative transfer analysis along with D/H ratios calculated using both HDCO and D₂CO for the first time. We also detect for the first time a large-scale velocity feature separate from that previously observed in CO 2−1 and, following the results of the application of a kinematic model, tentatively propose the detection of a second large-scale asymmetric molecular outflow launched by a wide-angle disk wind. The conclusions of this work can be summarised as follows:

• In this sample, we detect one H₂CO transition, 3(0,3)−2(0,2) at 218.222192 GHz; two HDCO transitions, 4(1,4)−3(1,3) at 246.9246 GHz and 4(2,2)−3(2,1) at 259.03491 GHz; and one D₂CO transition, 4(0,4)−3(0,3) at 231.410234 GHz. The undeuterated H₂CO isotopologue seems to be ubiquitous, with clear presence in both the disk and extended envelope of the source. There is also some absorption associated with this transition due to the cold envelope oriented along the line of sight. The emission also traces extended velocity components oriented along the NW–SE axis perpendicular to the previously detected outflow feature. The morphology of the deuterated isotopologue transitions are consistent with each other and show much more compact emission only originating from within the circumstellar disk.

• Radiative transfer analysis has shown that there are three velocity components associated with both deuterated and undeuterated formaldehyde in this source, which is consistent with other FAUST observations. The velocity values of these components are consistent with the conclusion that the two components showing strong redshift and blueshift are originating from the southern Source B and the low-velocity component is originating from the northern Source A, which is again consistent with other organic species. This is also consistent with the theory that Source B is associated with high-velocity small-scale streamer material and a preferential accretion given that it is believed to be the least massive of the two protostars.

• Radiative transfer analysis has also enabled the column densities of H₂CO, namely (3−8) × 10¹⁴ cm⁻², as well as both HDCO isotopologues, (0.8−2.9) × 10¹³ cm⁻², and the single D₂CO isotopologue, (2.6−4.3) × 10¹² cm⁻², to be obtained, enabling measurement of the D/H ratio towards this source for the first time; this was found to be ${0.02}_{-0.01}^{+0.02}$ averaged across the multiple HDCO velocity components. This D/H ratio in HDCO is consistent when compared to what was previously observed in a similar protobinary object.

• Following our analysis using a kinematic model (we refer readers to the Appendix for the full analysis), we propose that the extended features traced by the H₂CO emission are associated with a second asymmetric molecular outflow in addition to the one previously detected by Alves et al. (2017). Other molecular tracers of this feature need to be observed and analysed in order to distinguish between these two scenarios, with a study underway using FAUST observations of multiple dynamic tracers.

Acknowledgements

This project has received funding from (i) the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 811312 for the project “Astro-Chemical Origins” (ACO) and (ii) the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme, for the Project “The Dawn of Organic Chemistry” (DOC), grant agreement No 741002. This paper makes use of the following ALMA data: ADS/JAO.ALMA2018.1.01205.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST 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. I.J.-S. acknowledges support from grant No. PID2019-105552RB-C41 by the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe.” M.B. acknowledges the support from the European Research Council (ERC) Advanced Grant MOPPEX 833460.

Appendix A Kinematic streamer modelling

As mentioned in Sect. 3.1, our H₂CO emission clearly traces large-scale dynamic material extending along the NW–SE axis, which does not show alignment with any other feature previously detected in this source. It is already known that this source is dynamic as both small-scale streamers within the circumbinary disk (Alves et al. 2019) and a large-scale bipolar molecular outflow (Alves et al. 2017) have been seen; indeed, it is possible that interaction between the small-scale streamers and quiescent gas in the circumstellar region is producing a shocked region, causing the non-thermal sublimation of organic material (e.g. CH₃OH and H₂CO) from the dust grain surface (Vastel et al. 2022).

We implemented a streamline-generating model (Pineda et al. 2020) in order to probe this dynamic feature with the knowledge that, if the feature is indeed streamer-like in nature, kinematic information such as the free-fall timescale and infall rate could be obtained as well as the extent and velocity relation of the feature (as previously obtained for the small-scale streamers by Alves et al. (2019) using a similar method). However, due to a lack of a clear acceleration or deceleration, as well as a velocity relation that is flat or even outflowing, this investigation instead led to the conclusion that the feature is most likely not streamer-like in nature and is rather most likely a second molecular outflow in addition to the previously detected feature seen in CO 2-1 by Alves et al. (2017).

In order to investigate the kinematics of the dynamic features detected for the first time in this source, we used the same model that was previously used to investigate the streamer detected in the protobinary system Per-emb-2 in Pineda et al. (2020). In this source, the streamer is chemically fresh material composed of carbon-chain species originating from outside the dense core (more than 10,500 au) and infalling down to the disk scale. We focus on the region in [BHB2007] 11 with bright detections at large scales (approximately 2000 au) and therefore we require a minimum intensity threshold of 5σ for a pixel to be included in our modelling. To this end, prior to the modelling, we ran a pixel-by-pixel Gaussian fitting analysis using pyspeckit on the redshifted and blueshifted features (individually) in order to identify these pixels for the model; the resultant maps can be seen in Fig. A.1. As can be seen from these pixel maps, we are clearly tracing two distinct velocity features in both the redshifted and blueshifted regimes. The location and alignment of the high-velocity component suggest that this is originating from the disk. As confirmed in Sect. 3.1, this coincides with the high-velocity component of H₂CO emission (see also Fig. 2). Using this, we define polygons that exclude the high-velocity disk emission and focus on the region with a clear and cohesive dynamic component. These are the polygons that are used in the modelling.

The model generates streamlines assuming that material is free-falling (with conserved angular momentum) to a central object and uses analytic solutions from Mendoza et al. (2009) to generate Cartesian coordinates for the position and velocity of a parcel of mass within a streamline, where the xy plane is defined as the disk plane and the z axis is along the angular momentum vector. Along with the central mass of the system (taken from Vastel et al. 2022), the input parameters for the model are the spherical co-ordinates of the initial position and radial velocity of this parcel of mass within a cloud: r₀ (initial distance), v_r,0 (radial velocity), θ₀ (position angle with respect to the z axis) and ф₀ (initial angle within the disk plane), along with Ω₀ (initial angular velocity of the cloud itself). In order to obtain the position and velocity values with respect to the observer, two rotations are applied due to the inclination angle i and the disk position angle (both taken from Alves et al. 2017). r₀ is estimated by comparing the extent of the generated streamline to that seen in the integrated intensity maps shown in the upper panels of Figs. A.2 and A.3, ф₀ and θ₀ are estimated by comparing the trajectory of the generated streamlines to the data both in the plane of the sky (again, using the integrated intensity maps shown in the upper panels of Figs. A.2 and A.3) and through comparison between the velocity–distance relation of the generated streamline to the raw data (see the lower panels of Figs. A.2 and A.3, while v_r,0 is estimated by matching the propagation of velocity values along the length of the streamer, again using the velocity–distance maps shown in the lower panels of Figs. A.2 and A.3.

Fig. A.1 Centroid velocity maps from pixel-by-pixel Gaussian analysis showing pixels with intensity above 5σ for the redshifted (upper panel; integrated between 3.5 and 6 km s−1) and blueshifted (lower panel; integrated between 0 and 3.5 km s−1) features overlaid with continuum emission (blue contours) between 4σ and 300σ in 50σ intervals (σ=0.19 mJy beam−1). Polygons representing the pixels included in the kinematic model of the streamers are shown in red (upper panel) and cyan (lower panel). The ALMA synthesised beam is shown in the lower left corner.

Figures A.2 and A.3 show the results of the streamer modelling using these polygons with the initial parameters optimised to the data shown in Fig. A.1. The upper panel of each figure shows the modelled streamer superimposed onto the moment 0 map for H₂CO 3(0,3)−2(0,2) while the lower panel shows the streamline velocity as a function of distance from the source. We note that the modelled trajectory on the plane of the sky is reasonable, but that the velocity relation is not. The velocity relation is very hard to reproduce and it is possible that either a problem with the assumption of conservation of angular momentum or the feature being located too close to the plane of the sky could be the reason behind the limited ability to constrain the model in this case. From these lower panels, it should also be acknowledged that the velocity relation seen from these features does not exhibit the inward acceleration expected for a streamer feature. For the redshifted streamer, the velocity relation is almost completely flat, which could be another consequence of the proximity of the feature to the plane of the sky. For the blueshifted feature, we report a small outward velocity trend, which suggests an outflow rather than a streamer. However, as previously discussed, a bipolar outflow launched at symmetric positions with respect to the disk, approximately 110 au from the dust peak emission, has been traced by CO 2-1 in this source by Alves et al. (2017). This outflow is not spatially consistent with the features seen in H₂CO (refer to Fig. 2), while no evidence of the latter is seen in the outflow-tracing CO 2-1 observations. Assuming an infalling symmetric streamer, we would extract H₂CO spectra using the polygons shown in Figs. A.2 and A.3 in order to obtain values for the H₂CO column density within each modelled streamer and use this to obtain a value for the column density of H₂ (N(H₂) in cm⁻²) in each regime, from which we would use Eqs. A.1, A.2 and A.3 to obtain the free-fall timescale (t_ff) and infall rate (Ṁ_⋆), respectively: $$M_{TOT}=\mu m_H{d}^2\delta x\delta yN\left(H_2\right),$$(A.1) $$t_{ff}=\sqrt{\frac{R^3}{GM},}$$(A.2) $${\dot{M}}_{\star }=\sqrt{\frac{M_{TOT}}{t_{ff}},}$$(A.2)

where µ is the molecular weight (2.7, taking into account H₂, He and heavy element contributions), m_H is the atomic weight of hydrogen (1. 67 × 10⁻²⁴ g), d is the distance to the star (163 pc - Dzib et al. 2018), δ xδ y is the pixel size and M represents the mass enclosed by radius R (for which we use R = r₀ as found by the streamer model, G is the gravitational constant and L_bolrepresents the bolometric luminosity of the source (taken as 1 68 × 10²⁷ W - Sandell et al. 2021), all in cgs units. We could then compare our estimate for the accretion rate to the value of approximately 10⁻⁵ M_⊙ year⁻¹ previously obtained for this source by Alves et al. (2019) using a similar method. Therefore, material delivered via the potential infalling streamers might be comparable to the value delivered to the protostar. This is similar to what is seen in other objects (Pineda et al. 2020; Valdivia-Mena et al. 2022).

As a result, we propose that we are tracing a large-scale asymmetric molecular outflow formed as the molecular gas is displaced from the cavity and evacuated by a disk wind. It is difficult to pinpoint the driving points. We see an asymmetric outflow detected in CO launched at the disk edge (approximately 100 au), as well as asymmetric flowing material traced by H₂CO at a much larger (envelope) scale. The lack of blueshifted H₂CO emission in the northeastern region and the lack of red-shifted H₂CO emission in the southwestern region is puzzling but may be due to less gas being affected by the outflow due to lower densities and temperature. CO, which is tracing less dense material than H₂CO, may be hidden behind optically thick larger scale envelope emission which is being filtered out by the interferometer at velocities close to that of the systemic velocity. Recently, outflowing gas launched by an extended disk wind from a Keplerian disk with asymmetric, monopolar cavity walls has been detected (Bjerkeli et al. 2016), while another example of such a wide angle wind has been identified by Torrelles et al. (2011), although the latter towards a high-mass protostellar system. Therefore, wide-angle winds could be a natural outcome of the process of both low-mass and high-mass star formation.

Fig. A.2 Streamer modelling results for a blueshifted streamer. Upper panel: Streamline (blue) generated using best-fit initial parameters superimposed onto the H2CO 3(0,3)-2(0,2) moment 0 map first presented in the left panel of Fig. 1. Lower panel: Velocity of the streamer as a function of the distance from the source centre for raw data (black) and the model (blue), with the VLSR value shown in green.

Fig. A.3 Streamer modelling results for a redshifted streamer. Upper panel: Streamline (red) generated using best-fit initial parameters superimposed onto a moment 0 map. Lower panel: Velocity of the streamer as a function of the distance from the source centre for raw data (black) and the model (red), with the VLSR value shown in green.

Further, archival ALMA observations tracing CO 3-2 towards this source (project number 2019.1.01566) show morphological evidence of the dynamic feature along with the expected cavity that would be associated with such a wide-angle disk wind. Therefore, while it is not possible with the present data alone to conclude on the nature of the dynamic material traced for the first time here, other molecular tracers detected towards these features would be very useful in contributing further evidence; indeed, an investigation to this end will be published in an upcoming study (Martinez-Henares et al. in prep.). In particular, this study will combine two transitions of CO from the archive with FAUST observations of multiple dynamic tracers with the primary goal being to provide firm conclusions as to the nature of this newly detected dynamic feature of [BHB2007] 11.

References

All Tables

All Figures

Fig. 1 Comparing continuum and line emission. Left panel: H2CO 3(0,3)−2(0,2) moment 0 map integrated between −5 and 15 km s−1. This velocity range was chosen as this represents the full range of emission detected at or above 4σ. Right panel: 216 GHz continuum emission (colour scale) taken from the FAUST dataset with H2CO 3(0,3)−2(0,2) overlaid in pink contours. In both images, contours start at 4σ and end at 44σ in 8σ increments (σ = 17.6 mJy beam−1 km s−1), the positions of protostars A and B are represented by the black stars and the ALMA synthesised beam is shown in the lower left corner.

In the text

Fig. 2 Colour scale: H2CO 3(0,3)−2(0,2) centroid velocity map from pixel-by-pixel Gaussian analysis showing pixels with intensity above 5σ in the redshifted (left panel; integrated between 3.5 and 6 km s−1) and blueshifted (right panel; integrated between 0 and 3.5 km s−1) regimes; the systemic velocity is ~3.6 km s−1. Blue contours: 1.3 mm continuum contours starting at 4σ and ending at 300σ in intervals of 15σ (σ = 0.19 mJy beam−1) indicating disk emission. Red (left panel) and cyan (right panel) contours: Redshifted and blueshifted (respectively) CO 2−1 emission from Alves et al. (2017) detailing outflow emission seen at 8.67 km s−1 and −0.78 km s−1, respectively, representing the channels with the most extended outflow CO 2−1 emission. Contours are at 4, 6, 8, 10, 15, 20, 30, 40, 50, 60, 70, 80 and 90 times σ – which is 2.8 mJy beam−1 – in order to match the CO 2−1 observations, which have a beamsize of 0.2″ × 0.25″. The ALMA synthesised beam for the H2CO observations, with a beamsize of 0.43″ × 0.41″ is shown in the lower left corner.

In the text

Fig. 3 Comparing the multiple velocity components of the H2CO 3(0,3)−2(0,2) emission. Upper left: H2CO 3(0,3)−2(0,2) spectrum (black) with velocity resolution 0.17 km s−1 covering the full velocity range of −5 to 15 km s−1, extracted from the region where HDCO 4(2,2)−3(2,1) emission is above 4σ, overlaid with the best fit (pink) obtained using four Gaussian fits with FWHM, υLSR and Tex constrained using values obtained from fitting the HDCO transitions. A continuum value of 10.3 K has been added to this spectrum in order to model the absorption feature located at the systemic velocity of 3.6 km s−1. Upper right, lower left, lower right: H2CO 3(0,3)−2(0,2) moment 0 map integrated between −5 and 0 km s−1, between 0 and 8 km s−1 and between 8 and 15 km s−1, respectively, in order to individually show the velocity components of emission shown in the upper left panel. Contours begin at 4σ and are in increments of 8σ (σ = 17.6 mJy beam−1 km s−1). The pink contour illustrates the area within which the spectrum in the upper left panel was extracted, defined as the 4σ contour of the HDCO 4(2,2)−3(2,1) moment 0 map (see Fig. 4). The ALMA synthesised beam is in the lower left corner of all three maps.

In the text

Fig. 4 Emission maps for HDCO 4(2,2)−3(2,1). Left panel: HDCO 4(2,2)−3(2,1) moment 0 map integrated between −1 and 7 km s−1 that peaks at the position of the dust emission. Observations performed using both 12 m and 7 m ALMA configurations are combined in this map in order to investigate both large- and small-scale structure. Contours start at 4σ and end at 10σ in 3σ increments (σ = 4.1 mJy beam−1 km s−1). The 4σ contour is highlighted in cyan as this is the contour that served as the extraction polygon for all spectra shown in this work. Right panel: Moment 1 map of HDCO 4(2,2)−3(2,1) integrated between −1 and 7 km s−1, in order to directly compare with the moment 0 map and to show the velocity structure of the disk. In both maps, the locations of the protostars A and B are shown by black stars. The ALMA synthesised beam is in the lower left corner.

In the text

Fig. 5 HDCO spectrum (black) with a velocity resolution of 0.14 km s−1 for the 4(2,2)−3(2,1) transition overlaid with the best fit (pink) from our simultaneous radiative transfer analysis, extracted within the 4σ contour shown in Fig. 4. The analysis was run twice, assuming upper limits for the size of the emitting region of 1.0″ and 0.3″ for Sources A and B, respectively and lower limits of 0.5″ and 0.15″ for Sources A and B, respectively; the figure shown represents the latter. Dashed lines denote the previously identified velocity components at −2, 2.8 and 9.9 km s−1.

In the text

Fig. 6 HDCO spectrum (black) with a velocity resolution of 1.18 km s−1 for the 4(1,4)−3(1,3) transition overlaid with the best fit (pink) from our simultaneous radiative transfer analysis, extracted within the 4σ contour shown in Fig. 4. The analysis was run twice, assuming upper limits for the size of the emitting region of 1.0″ and 0.3″ for Sources A and B, respectively and lower limits of 0.5″ and 0.15″ for Sources A and B, respectively; the figure shown represents the latter. We note that the third component is not included in the fit because of the blending with a transition of CH3OCHO 19(4,15)−18(4,14) A, the emission components of which are clearly seen above 6 km s−1 and highlighted explicitly in Fig. 7. Dashed lines denote the previously identified velocity components at −2, 2.8 and 9.9 km s−1.

In the text

Fig. 7 Extracted spectrum centred on the rest frequency of the CH3OCHO 19(4,15)−18(4,14) A transition at 246.9147 GHz with dashed lines showing the predicted locations of the velocity components assuming similar VLSR to methanol components (−2, 2.8 and 9.9 km s−1). The locations of the same velocity components (−2, 2.8 and 9.9 km s−1) in HDCO 4(1,4)−3(1,3) are shown by the pink dashed lines. The −2 km s−1 component of CH3OCHO 19(4,15)−18(4,14) A is blended with the 9.9 km s−1 component of HDCO 4(1,4)−3(1,3).

In the text

Fig. 8 Moment 0 map of HDCO 4(1,4)−3(1,3) integrated between -5 and 8 km s−1. This velocity range was chosen to include the two unblended HDCO emission components shown in Fig. 6. Contours start at 4σ and end at 16σ in 4σ increments (σ = 8.5 mJy beam−1 km s−1). The pink contour illustrates the area within which the spectrum in Fig. 6 was extracted, defined as the 4σ contour shown in Fig. 4. The locations of protostars A and B are shown by the black stars. The ALMA synthesised beam is shown in the lower left corner.

In the text

Fig. 9 D2CO 4(0,4)−3(0,3) spectrum (black) with a velocity resolution of 0.16 km s−1 extracted from the 4σ contour shown in Fig. 4 and smoothed three times with the Hanning algorithm overlaid with best fit (pink) obtained using single-component Gaussian fitting with FWHM, VLSR and Tex constrained using values obtained from fitting the HDCO transitions. The analysis was run twice, assuming upper limit values for the size of the emitting region of 0.3″ (covering source B emission) and 1.0″ (covering source A emission); the figure shown represents the latter. The rms level (0.15 K) is shown in green.

In the text

Fig. 10 Comparison of H2CO 3(0,3)–2(0,2), D2CO 4(0,4)–3(0,3), HDCO 4(1,4)–3(1,3) and HDCO 4(2,2)–3(2,1) lines with the baselines subtracted. The dashed lines represent velocities of –2, 2.8 and 9.9 km s−1. The CH3OCHO 19(4,15)–18(4,14) A transition (246.9147 GHz) that is near HDCO 4(1,4)–3(1,3) is shown in red, while the H13CN (3−2) transition (259.0118 GHz) in the same spw as HDCO 4(2,2)–3(2,1) is denoted in green.

In the text

	Fig. 11 H2CO moment 1 map in colour scale with five pixels overlaid in black and white (for clarity in this case) triangles, chosen as the pixel at the location of the dust peak (A), along with the most intense pixels for each velocity range (B: 0−2 km s−1, C: 2.5−3.5 km s−1, D: 3.5−4 km s−1 and E: 5−6 km s−1). The ellipse in the bottom-left corner represents the ALMA synthesised beam.
In the text

	Fig. 12 Spectra extracted from the five pixels (labelled A, B, C, D and E) indicated in Fig. 11. Clockwise from upper left: H2CO (218 GHz), HDCO (259 GHz), D2CO (231 GHz) and HDCO (246 GHz).
In the text

Fig. A.1 Centroid velocity maps from pixel-by-pixel Gaussian analysis showing pixels with intensity above 5σ for the redshifted (upper panel; integrated between 3.5 and 6 km s−1) and blueshifted (lower panel; integrated between 0 and 3.5 km s−1) features overlaid with continuum emission (blue contours) between 4σ and 300σ in 50σ intervals (σ=0.19 mJy beam−1). Polygons representing the pixels included in the kinematic model of the streamers are shown in red (upper panel) and cyan (lower panel). The ALMA synthesised beam is shown in the lower left corner.

In the text

Fig. A.2 Streamer modelling results for a blueshifted streamer. Upper panel: Streamline (blue) generated using best-fit initial parameters superimposed onto the H2CO 3(0,3)-2(0,2) moment 0 map first presented in the left panel of Fig. 1. Lower panel: Velocity of the streamer as a function of the distance from the source centre for raw data (black) and the model (blue), with the VLSR value shown in green.

In the text

	Fig. A.3 Streamer modelling results for a redshifted streamer. Upper panel: Streamline (red) generated using best-fit initial parameters superimposed onto a moment 0 map. Lower panel: Velocity of the streamer as a function of the distance from the source centre for raw data (black) and the model (red), with the VLSR value shown in green.
In the text