© The Authors 2023

1 Introduction

Observing protostellar envelopes around very young protostars is fundamental to gaining a better understanding of the progression of the collapse of protostellar cores toward planetary systems. Class 0 objects are known to represent the youngest stage of protostellar evolution (André et al. 1993, 2000). Most of their mass is contained in a dense envelope that accretes to the central protostar during the main accretion phase (Maury et al. 2011; Evans et al. 2009). Protostars are deeply embedded in their parent cores, which may cause interactions between protostellar outflows and surrounding gas leading to complex morphologies. Detailing these structures will allow us to learn about the dynamics of protostellar evolution at an early stage. In light of these circumstances, it becomes essential to investigate in detail the earliest stages of star formation for specific sources.

The study of the extent and contribution of magnetic fields in star formation and the competition between magnetic and turbulent forces is still a highly debated topic in modern astronomy (e.g., Mac Low & Klessen 2004; McKee & Ostriker 2007; Crutcher 2012). However, in the star formation process, especially during the early stages, magnetic fields (B) are expected to play a crucial role, providing a source of non-thermal pressure against the gravitational pull (McKee & Ostriker 2007). In light of the fact that interstellar gases are often mildly ionized (Caselli et al. 1998), the matter is likely to be coupled with the magnetic field lines at envelope scales. Due to gravity, magnetic lines bend inward, thus producing an hourglass shape, and in low-mass star-forming regions this effect is not detected frequently (detected in 30 percent of young stellar objects in polarization; Hull & Zhang 2019, Pattle et al. 2022).

IRAS15398-3359 (hereafter IRAS15398) is a low-mass Class 0 protostar at a distance of 156 pc (Dzib et al. 2018), embedded in the Lupus I molecular cloud, α₂₀₀₀ = 15^h43^m02^s.2, δ₂₀₀₀ = –34°09'07.7". It has a bolometric temperature of 44 K (Jørgensen et al. 2013). The protostellar mass is ${0.007}_{-0.003}^{+0.004}{M}_{\odot }$ (Okoda et al. 2018). Lupus I is the least evolved component of the Lupus complex (Rygl et al. 2013), and optical polarization studies have demonstrated that Lupus I is threaded by a very ordered magnetic field that is perpendicular to its filamentary extension (Franco & Alves 2015). Therefore, it is an ideal place to study the kinematics of the early stages of low-mass star formation and the connection between the source kinematics and the strong large-scale magnetic field. By observing its CO emission line with the single-dish and interferometric observation, a molecular outflow was detected from this source (Tachihara et al. 1996; Bjerkeli et al. 2016; van Kempen et al. 2009). The core is embedded in a less dense (N(H₂) ~ 10²² cm⁻²) filamentary structure, which extends toward the northwest.

Based on magnetic field studies in protostellar core simulation analysis, more magnetized cores show strong alignments of the outflow axis with the magnetic field orientation, whereas less magnetized cores display more random alignment (Lee et al. 2017). Observational results present a mixture of cases. Galametz et al. (2018) used a sample of 12 low-mass Class 0 protostars to investigate the submillimeter polarized emission at scales of ~600–5000 au, and demonstrated a relation between the field morphology, the core rotational energy, and the multiplicity of the protostellar system. According to that paper’s analysis, the envelope scale magnetic field tends to be either aligned or perpendicular to the outflow direction, but for single sources the magnetic field is aligned along the outflow direction. Yen et al. (2021) studied 62 low-mass Class 0 and I protostars in nearby (<450 pc) star-forming regions with the orientations of the magnetic fields on 0.05–0.5 pc scales. They suggest that the outflows are likely to be misaligned with B-fields by 50 degrees in 3D space. While Hull & Zhang (2019) used Atacama Large Millimeter/submillimeter Array (ALMA) observations with spatial resolutions of up to ~100 au, and conclude that magnetic fields and outflows are randomly aligned in low-mass protostel-lar cores. The discrepancy between simulations and observations can be due to the limitations of the simulation setup. As an example, Lee et al. (2017) applied an ideal magnetohydrodynamic (MHD) simulation when in reality non-ideal MHD effect might be important (Wurster 2021).

Redaelli et al. (2019a) used polarimetric observations of the dust thermal emission at 1.4 THz obtained with the Stratospheric Observatory for Infrared Astronomy (SOFIA) telescope to investigate the magnetic field properties at the core scales toward IRAS15398. The authors found a uniform magnetic field consistent with the large-scale field derived from optical observations (Franco & Alves 2015). They suggested the core experienced a magnetically driven collapse and the core inherited the B-field morphology from the parental cloud during its evolution. The field lines pinch inward toward the central object, leading to the characteristic hourglass shape that is predicted by models of magnetically driven collapse. They showed that the mean direction of the magnetic field is aligned with the large-scale B-field and with the direction of outflow. Their prediction for a magnetic field strength of B = 78 μG is expected to be accurate within a factor of two. They calculated the mass-to-flux ratio, λ = 0.95, which means that the core is in a state of transition between supercritical and subcritical states.

In this paper we present new observational data that allow us to study the gas kinematics, which we compare to the magnetic field direction. The aim is to assess the importance of magnetic fields in the dynamical evolution of low-mass star-forming regions. The outline of the paper is as follows. The observations and results are described in detail in Sects. 2 and 3. In Sect. 4 we analyze the observed line profiles of C¹⁸O (2–1) and DCO⁺ (3–2) emission lines. The results are summarized in Sect. 5.

2 Observations

2.1 APEX

IRAS15398 was observed using the Atacama Pathfinder Experiment (APEX) single-dish antenna located at Llano de Chajnantor in the Atacama desert of Chile on 2019 September 14, 16, 17, 21, and 23. We used the PI230 receiver coupled with the FFTS4G backend in the on-the-fly mode. We used the highest spectral resolution that the FFTS4G could provide, 62.5 kHz (≈0.08 km s⁻¹ at the frequency of the DCO⁺ (3–2) line). The data were reduced to a pixel size of 9 arcsec. The broad bandwidth of the PI230 receiver can be set up to observe simultaneously C¹⁸O (2–1) and the DCO⁺ (3–2) transitions at 219.560 and 216.113 GHz, respectively, the target lines of this research. The angular resolution at these frequencies is ~28 arcsec, corresponding to 0.02 pc at the source distance 156 pc (Dzib et al. 2018). The data reduction was performed based on the standard procedure of the CLASS software, GILDAS¹. The antenna temperature, T_A, was converted to main-beam brightness temperature using the forward efficiency (η_fw) and main beam efficiency, $T_{\text{mb}}=T_{\text{A}}\frac{{\eta }_{\text{fw}}}{{\eta }_{\text{mb}}}$, given a main beam efficiency η_mb = 0.8².

Fig. 1 Integrated intensity of C18O (1–2) (top) and DCO+ (3–2) (bottom) toward IRAS 15398. The white contour levels are 10, 20, and 30 times of mean rms value for the DCO+ line and 40, 50, and 60 times of mean rms value for the C18O (2–1) line. The beam size is shown in the bottom left corner. The yellow contours in both images show H2 column density (levels:[1.0, 1.5, 2.0] 1022 cm−2). The black star gives the position of the protostar.

2.2 Herschel

We used archival data from the Gold Belt Survey, obtained with the Herschel space telescope to obtain the gas column density and the dust temperature (Bontemps et al. 2010; Rygl et al. 2013; Benedettini et al. 2018). The N(H₂) column density map has a resolution of ~38 arcsec.

3 Result

Figure 1 presents the moment 0 (integrated intensity) maps of C¹⁸O (2–1) and DCO⁺ (3–2) overlaid with the contours of H₂ column density. The mean rms in T_mb scale is 0.1 K and 0.08 K for C¹⁸O (2–1) and DCO⁺ (3–2), respectively. On the basis of emission-free channels only, we derived the mean rms per channel for each transition. This rms is used to associate uncertainties at each pixel to the integrated intensity, the mean value for C¹⁸O (2-1) and DCO+ (3-2) being 0.03 K km s⁻¹ and 0.02 K km s⁻¹, respectively. These lines are optically thin and do not present crowded hyperfine structures. In order to confirm this hypothesis, we evaluate the opacity as $$\tau =-\ln \left[1-\frac{T_{\text{MB}}}{J_v\left(T_{\text{ex}}\right)-J_v\left(T_{\text{bg}}\right)}\right],$$(1)

where T_ex = 12 K is the excitation temperature (obtained from the Herschel dust temperature map, see also Sect. 4.1), T_MB = 3.5 K is the peak main beam temperature, the function J_v is the equivalent Rayleigh-Jeans temperature, and T_bg = 2.73 K is the cosmic background temperature. Therefore, we obtain τ = 0.66. For the DCO⁺ (3–2) line we obtain τ = 0.44 using = 1.0 and T_ex = 7 K. Through a combination of excitation and abundance, distinctive species give complementary information on gas conditions. Due to its relatively high abundance and low critical density n_cr(C¹⁸O (2–1)) ~ 10⁴cm⁻³, computed with numbers in the LAMDA database³, C¹⁸O (2–1) is a sensitive tracer of relatively low-density material in the cloud, which traces the more extended gas in the filamentary structure. The DCO⁺ (3–2) molecule, on the other hand, presents a higher critical density n_cr(DCO⁺(3–2)) ~ 10⁶cm⁻³, which makes it more selective of dense gas closer to the central protostar. DCO⁺ (3–2) emission is also known as a remarkably sensitive tracer for gas properties during the early stages of protostellar evolution (e.g., Gerner et al. 2015). In the location of the protostar, we see a decrease in the C¹⁸O (2–1) integrated intensity, suggesting that the molecule is partially depleted onto the dust grains (Caselli et al. 1999; Bacmann et al. 2002). Figure 2 shows the channel maps of the C¹⁸O (2–1) line. The figure presents the signal emission at velocity intervals of ≈0.2 km s⁻¹.

Fig. 2 Channel maps of C18O (2–1) emission. The black cross indicates the position of the protostar core. The velocity of each channel is shown above each panel.

4 Analysis

4.1 Column density maps

To calculate the column density map of these two molecules we used the same procedure as used in Redaelli et al. (2019b) and Caselli et al. (2002) for an optically thin transition. The DCO⁺ (3–2) and C¹⁸O (2–1) lines are both optically thin, as we show in the previous section. The expression of the total column density derived by an optically thin transition is given by $$N_{\text{col}}=\frac{8\pi W{v}^3}{c^3{A}_{\text{ul}}}\frac{Q\left(T_{\text{ex}}\right)}{J_v\left(T_{\text{ex}}\right)-J_v\left(T_{\text{bg}}\right)}\frac{\frac{{\text{E}}_{\text{u}}}{e^{{\text{k}}_{\text{B}}\text{Tex}}}}{\left(\frac{\text{h}v}{e^{{\text{k}}_{\text{B}}\text{Tex}}-1}\right)},$$(2)

where T_ex is the excitation temperature, the function J_v is the equivalent Rayleigh-Jeans temperature, T_bg = 2.73 K is cosmic background temperature, E_u is the upper state energy, g_u is the degeneracy, A_ul is the Einstein coefficient, Q is the partition function, vis the line frequency, h is the Planck constant, k_B is the Boltzmann constant (see Table 1 for details)⁴, and W is the integrated intensity of the line. Since the DCO⁺ (3– 2) transition shows only one velocity component, we use the result of the Gaussian fit to compute the integrated intensity of this line (see Sect. 4.2 for more details) by calculating the area under the Gaussian profile. The C¹⁸O (2–1) emission, instead, shows signs of multiple velocity components along the line of sight. We therefore compute the integrated the intensity from the data cube, integrating emission over the velocity range [4–6.5] km s⁻¹, which contains the whole line profile.

We use the dust temperature map to approximate the excitation temperature for the C¹⁸O (2–1) line, obtained from Herschel data (Benedettini et al. 2018; Rygl et al. 2013) since we expect this line to be thermally excited. This assumption may induce some small errors as at the volume densities traced by the C¹⁸O (2–1) line the gas and dust are not thermally coupled (Goldsmith 2001). We therefore use the dust temperature as a proxy for the gas kinetic temperature. On the contrary, the DCO⁺ (3–2) line is the 3–>2 transition, which has a high critical density; therefore, the dust temperature is not a good approximation for excitation temperature because we expect the line to be sub-thermally excited, so we use excitation temperature equal to 7 K with a variation of 2 K for this line. The column density peak for the DCO⁺ and H₂ is found at the protostar position. For the protostar position, the C¹⁸O, DCO⁺, and H₂ column densities are (7.8 ± 0.1) × 10¹⁴cm⁻², ${1.4}_{-0.5}^{+2.24}\times {10}^{12}{\text{cm}}^{-2}$, and (4.2 ± 1.8) × 10²² cm⁻² (Roy et al. 2014), respectively. The column density for H₂ is obtained based on the Herschel map.

4.2 Spectral line analysis

In order to derive the kinematics parameter maps (e.g., V_lsr, σ_V), we perform a Gaussian fitting of the transitions using the PYSPECKIT package of python (Ginsburg & Mirocha 2011). For the DCO⁺ (3–2) data cube, we use a single-Gaussian component fit. The initial guesses are then 2.5 K, 5.2 km s⁻¹, and 0.2 km s⁻¹ for the amplitude, velocity dispersion, and width, respectively.

Instead, C¹⁸O (2–1) presents more comple× kinematics. It often shows two velocity components in its profiles. Since the line is optically thin, as shown in Sect. 3, we are confident that these are multiple velocity components and they are not due to self-absorption. In order to fit two Gaussian profiles on the C¹⁸O (2–1) data we perform a simple S/N cut, we mask those pixels where $\frac{T_{\text{MB}}^{\text{peak}}}{\text{rms}}<20$, then fit one Gaussian profile to all the unmasked pixels (65% of the pixels). Then we fit two Gaussians for those pixels that had a residual larger than 2 × rms and those Gaussian fits with a width broader than 0.25 km s⁻¹ in the previous step. By visual inspection we find that lines broader than 0.25 km s⁻¹ show profiles consistent with two velocity components on the line of sight. For the second time we check the Gaussian fit profiles and for those pixels that have residuals bigger than 2 × rms we do the fitting one more time with different initial guesses (60% of the pixels with S/N cut), which means they always have residuals less than two times the rms. As a final step, if the error on the velocity dispersion or velocity and on the amplitude is larger than 1 km s⁻¹ and 1 K, respectively, we remove the fit. By doing this, we remove fits with unreasonably large uncertainties, which indicates that they have been poorly fitted. Figure 3 shows the fit results overlapped with the data.

The grid of spectra of C¹⁸O (2–1) and DCO⁺ (3–2) lines for 40 positions at 18 arcsec intervals from each other around the core is shown in Fig. 3 (for more details about the position of each spectrum, see Appendi× A). The red histogram represents the DCO⁺ (3–2) spectrum and the black histogram is the C¹⁸O (2–1) spectrum. The blue curves represent the fit when the code uses a two-Gaussian and the green curves are the fit when the code used one Gaussian for the C¹⁸O (2–1) line. The vertical line is V_sys = 5.2 km s⁻¹, which is the systematic velocity computed from C¹⁸O (2–1) line which was reported in Yen et al. (2017, Vsys = 5.24 ± 0.03 km s⁻¹), which is consistent with the systematic velocity seen in our data.

The DCO⁺ (3–2) spectra always show a single velocity component with a Gaussian profile. This line traces only one velocity component that is associated with the high-density material. On the other hand, as already mentioned, C¹⁸O (2–1) shows two velocity components in their profiles. The brightest component is the one close to the velocity of the DCO⁺ (3–2) line, and it hence arises from the high-density material. The fainter component of C¹⁸O (2–1) can appear on the red or blue side of the stronger component depending on the location in the cloud. In total, there are three different velocity components, the main one at 5.2 km s⁻¹ (in all panels on all DCO⁺ (3–2) spectra locations), a lower velocity component seen in the northwestern panels, and a higher velocity component seen in the southwestern panels. In the first two rows of Fig. 3 the faint component appears on the blue side, and in the two bottom rows of the spectrum grid it appears on the blue side of the stronger component.

In Fig. 3 we can see the line profiles change across the grid. The multiple velocity components of C¹⁸O (2–1) along the line of sight appear to merge from east to west (moving toward the position of the protostar). C¹⁸O (2–1) is a lower density tracer than DCO⁺ (3–2); these fainter velocity components, seen only in the former tracer, are likely additional nearby low-density filamentary clouds along the line of sight. The core envelope then appears to be found in the correspondence of the merger of these structures. In the following subsections, we discuss the maps of the kinematics parameters for each tracer individually.

4.3 DCO⁺ (3–2) line

In Fig. 4 (top panel) we show the velocity dispersion map, which presents a clear increase toward the center of the protostar envelope, starting with very narrow lines in the outer region (around the filament). The mean velocity dispersion derived from DCO⁺ (3–2) of the gas in the filament is in fact 〈σ_v〉 = 0.12 ± 0.02 km s⁻¹, but it becomes broader toward the center; for the positions where N(H₂) > 3.2 × 10²² cm⁻² (inner contour in the Fig. 4), we derive 〈σ_v〉 = 0.175 ± 0.006 km s⁻¹. This increase is linked to the protostellar activity, injecting turbulence, and it could also be caused by the rotation of the core.

We also observe oscillatory motions in the velocity field (visible in the C¹⁸O data as well, see Fig. 7, the green dots), which have been seen before in the large-scale velocity patterns. This velocity pattern is consistent with core-forming motions (Hacar & Tafalla 2011). A small velocity gradient can be seen along the filament, which could be linked to the ongoing accretion material toward the central object.

The filament extends over ~0.1 km s⁻¹ in velocity toward the protostar from 5.1 km s⁻¹ in the west side to 5.0 km s⁻¹ at the eastern edge. These positions are shown with red plus signs in Fig. 4. We determine a conservative estimate of the length of the filament of 0.11 pc. This value is computed based on the H₂ column density from Herschel data, an area with a higher value than 1.5 × 10²² cm⁻². In this border we find almost all of the filament (second contour in Fig. 4). Thereby the velocity gradient is ∇V = ∆V/∆R = 0.91 ± 0.23 km s^–1pc⁻¹. To compute the mass of filament we employ this relation: $$M_{\text{fil}}={\displaystyle \sum _{i=1}N{\left({\text{H}}_2\right)}_i\times \mu \times m_{\text{H}}\times \text{A.}}$$(3)

In this case, we obtain M_fil = 1.04 M_⊙, where N(H₂) is gas column density of the i_th pixel and A is the area of the pixel; μ = 2.8 and m_H are gas mean molecular weight per hydrogen molecule and hydrogen mass, respectively. Therefore, we quantify the mass accretion rate along the filament, ${\dot{M}}_{\text{acc}}=M_{\text{fil}}\times \nabla V=9.7\times {10}^{-7}{M}_{\odot }{\text{yr}}^{-1}$. This value is affected by some sources of uncertainty. We consider a 2% error on the distance (Dzib et al. 2018), a 12% error on the calibration, and a 40% error due to the assumption of the dust opacity inde× (Benedettini et al. 2018; Roy et al. 2014). The uncertainty on the mass estimation with all these three errors is a total of 42%. With the 26% error on the velocity gradient, we get a final relative error of 49% for the value of the mass accretion rate. In addition, the inclination of the filament with respect to the plane of the sky is unknown, and it has an influence on the value of ${\dot{M}}_{\text{acc}}$ by a factor of tan i (see, e.g., Chen et al. 2019). The derived value of the accretion rate changes by up to 70% if the inclination varies between 30 and 60 degrees. Therefore, considering these factors, we assume that the derived ${\dot{M}}_{\text{acc}}$ value is accurate within a factor of two. This mass accretion rate is comparable to the value that Pineda et al. (2020) found for a streamer of material connecting a protostar, Per-emb-2 (IRAS 03292+3039), to the surrounding cloud $\left({\dot{M}}_{\text{acc}}={10}^{-6}{M}_{\odot }{\text{yr}}^{-1}\right)$. The accretion rate that we calculate here is for a filament, not a streamer. However, this filament is diffuse and not very dense. Furthermore, the value we found is comparable to the mass infall rate of protostellar envelopes estimated in other young objects ~3 × 10⁻⁶ (Evans et al. 2015).

We detect another velocity gradient around the protostar in the west–east direction, which could be due to the rotation of the core. In this scenario the rotation axis would lay in the north-south direction (PA = 16°; see last paragraph of the section), close to the direction of the detected bipolar outflows found by Bjerkeli et al. (2016), which are shown with the blue and red arrows in the bottom panel of Fig. 4 (PA = 35°). Redaelli et al. (2019a) also showed that this outflow direction is almost parallel to the mean magnetic field, which has the direction of PA = 45±7°. The black vectors in the bottom panel of Fig. 4 represent the polarization angles, tilted by 90 degrees to trace the magnetic field direction. The rotation axis and magnetic field lines are moderately well aligned (with an offset of 29°). According to the magnetohydrodynamic (MHD) collapse models, magnetic braking should be effective in this cloud (Joos et al. 2012; Li et al. 2013; Krumholz et al. 2013; Seifried et al. 2015), which is consistent with the absence of a resolved Keplerian disk (down to 30 AU, Yen et al. 2017). Okoda et al. (2018) observed IRAS 15398-3359 with better resolution (0.2" angular resolution). A disk of no more than 30 au in size was detected in their analysis with a mass between 0.006 and 0.001 M_⊙. This is a very small disk that is consistent with magnetic braking. Magnetic braking is in fact an efficient way to remove angular momentum from infalling and rotating material, suppressing envelope fragmentation and the formation of large disks (Li et al. 2014).

The velocity gradient over the protostar shows gas motions likely consistent with the rotation of the core. To analyze the gas motions further, we employ a two-dimensional representation of the velocity gradients. Figure 5 shows the centroid velocity of DCO⁺ (3–2), overlaid with the gradient velocity arrows. The arrow length represents the relative vector magnitude of the gradient and the direction of the arrows point in the direction of increasing velocity. The method used here followed the analysis of Goodman et al. (1993), and developed by Caselli et al. (2002), to find the velocity gradient based on a linear fit between offset declination and the velocity. By assuming that the core rotates as a solid body, would only depend on the coordinates in the sky and not on the distance along the line of sight. In this approximation, the centroid velocity of the line is a linear function of the offset on the plane of the sky $$V_{\text{LSR}}=V_0+a\text{Δ}\alpha +b\text{Δ}\delta \text{,}$$(4)

where Δα and Δδ represent offsets in right ascension and declination, and V₀ is the systemic velocity of the cloud with respect to the local standard of rest. The coefficients a and b, together with V₀, can be obtained by least-squares fitting. The velocity gradient magnitude is then (Goodman et al. 1993) $$\nabla V_{\text{LSR}}=\sqrt{\left(a^2+b^2\right)},$$(5)

and its direction toward increasing velocity Θ_∇V is given by $${\text{Θ}}_{\nabla V}={\tan }^{-1}\frac{a}{b}.$$(6)

We use this method and obtain the velocity gradient and the position angle with their uncertainties. The number of pixels used to carry out the fit is 9, which is appropriate for single-dish data that is Nyquist sampled (Caselli et al. 2002).

We only consider values with S/N > 3 in the final result for each velocity gradient value; the mean signal-to-noise ratio around the protostar position is 7. The mean velocity gradient magnitude distribution around the core peaks at 5.1 ± 0.7 km s⁻¹ pc⁻¹ and has a mean position angle of (106 ± 7)° (counterclockwise from the north toward the east), which is shown as a yellow arrow in Fig. 4. We can see that this mean velocity gradient is consistent at the 3 × σ level to be perpendicular to the direction of the bipolar outflow found by Bjerkeli et al. (2016, PA = 35°). We assume that the total gradient is the rotational direction of the core. This method is implemented in a Python code which is available via open access on GitHub⁵.

Fig. 3 Spectrum grid of C18O (2–1) (black histogram) and DCO+ (3–2) (red histogram) for 40 random positions around the core. Shown are the two-Gaussian fits (blue curves) and one-Gaussian fits (green curves) to the C18O (2–1) line. The vertical line is Vlsr = 5.2 km s−1. The blue star shows the spectrum for the position of the central protostar.

Fig. 4 Maps of kinematic parameters. Top: gas velocity dispersion of the cloud traced using the DCO+ (3-2). Bottom: centroid velocity map obtained fitting the observed DCO+ (3-2). The black vectors represent the polarization angles, tilted by 90 degrees to trace the magnetic field direction (Redaelli et al. 2019a). The blue and red arrows show the direction of the outflow (PA = 35°, from Bjerkeli et al. 2016), and the yellow arrow presents the mean velocity gradient direction around the core. The star shows the position of the protostar. The contours represent N(H2) column density levels, as derived from Herschel data: [1.0, 1.5, 2.0] 1022 cm−2.

Fig. 5 Centroid velocity map of the DCO+ (3–2) line overlaid with the gradient arrows (only vectors with S/N > 3 are shown). The arrow length represents the relative vector magnitude of the gradient, according to the scale shown in the top left corner, and the direction of the arrows points to the steepest velocity field change.

Fig. 6 Maps of kinematic parameters. Left panel: bright component of centroid velocity map of C18O (2–1). Right panel: centroid velocity map only showing the weak component of the Gaussian fitting. The black vectors represent the polarization angles, tilted by 90 degrees. The star shows the position of the protostar. The contours represent N(H2) column density levels, as derived from Herschel data: [1.0,1.5,2.0] × 1022 cm−2.

4.4 C¹⁸O (2–1) line

Figure 6 represents the centroid velocity map obtained fitting the observed C¹⁸O (2–1) spectra. The left panel of Fig. 6 represents the brightest components of the Gaussian fitting, and as we discussed in Sect. 4.2, in some locations of the cloud C¹⁸O (2–1) gas reveals two velocity components. The right panel of Fig. 6 shows the second component, the one with lower intensity. Depending on the location on the map, this less bright component has red or blue velocities compared to the main component. The right panel of Fig. 6 shows that by moving from north to south the faint component is going to the red side of the brightest components. The red velocities appear mostly on the south of the filament and the west side of the protostar position. This is similar to our discussion in Sect. 4.2, where there are three components along the line of sight: the bright one associated with the core and two fainter ones, one in the north and one in the south of the filament.

A small velocity gradient can also be seen in the brightest component of the C¹⁸O (2–1) emission around the protostar in the east-west direction, in agreement with what we discuss in Sect. 4.3 for the DCO⁺ (3–2) line, which is likely associated with the rotation of the core. For positions where N(H₂) > 2 × 10²² cm⁻², the mean value of the centroid velocity is 〈V_lsr〉 = 5.22 km s⁻¹ with an uncertainty 0.04 km s⁻¹. The V_lsr value at the west side of the core is 5.169 ± 0.005 km s⁻¹ and increases toward the east side of the core up to the value 5.290 ± 0.003 km s⁻¹. These two values are calculated at the edge of the core, where the N(H₂) is equal to 2 × 10²² cm⁻² (the inner contour in Fig. 6).

Frau et al. (2015) proposed for the Pipe nebula that the sharp changes in the magnetic field is produced by shocks between two clouds and, in comparison to the non-shocked gas, the column density and magnetic field strength double. Redaelli et al. (2019a) observed a sharp change in the magnetic field toward dust extension in the northwest direction of this cloud with respect to the core magnetic field. We predict that this is produced by merging two distinct clouds, two components of the C¹⁸O (2–1) line. These polarization vectors, which show different directions concerning the vectors in the core position, could result from the cloud collision.

4.5 Comparison between the C¹⁸O (2–1) and DCO⁺ (3–2) kinematic

The result of the fitting procedure is shown in Fig. 7, where we present an image depicting a 3D position-position-velocity (PPV) diagram highlighting the distributions of DCO⁺ (3–2) and C¹⁸O (2–1) gas throughout the cloud. Each data point illustrates the location and centroid velocity of an independent Gaussian component. The color of each data point relates to each spectral component (discussed in Sect. 4.2). Orange refers to DCO⁺ (3–2) emitting gas, and the others refer to C¹⁸O (2–1) emitting gas. The velocity structure of the C¹⁸O (2–1) emission is quite complex. Four velocity components are displayed in total. Overall, for C¹⁸O (2–1) we used only two-component fitting, but here we display it in different categories. The blue data correspond to the brightest components of the C¹⁸O (2–1) gas, and the red and cyan points instead are related to the secondary fainter component, when it is located at lower and higher velocities, respectively, with respect to the brightest one.

We observe a systematic difference between the centroid velocity of DCO⁺ (3–2) and the main component of C¹⁸O (2–1), suggesting that they are not tracing exactly the same gas. We speculate that this is related to the fact that one is an ion and the other is a neutral species, and they behave differently concerning the magnetic fields, or it could be because of the difference between their gas densities. It is important to note that if ions and neutrals behave differently at the same density, this indicates a violation of the flux-freezing assumption. We discuss this point in more detail in Sect. 5.

Figure 7 shows that the main velocity component of the C¹⁸O (2–1) data and the DCO⁺ (3–2) line present a systematic velocity shift. The C¹⁸O (2–1) data always appear at a higher velocity than the DCO⁺ (3–2) spectra. In order to investigate this point further, we show in Fig. 8 the velocity difference V_lsr(C¹⁸O) – V_lsr(DCO⁺) with the H₂ column density overlaid on top. In addition, we did not find any correlation between the velocity shift and the H₂ column density map. In this figure the shift between C¹⁸O (2–1) and DCO⁺ (3–2) is clearly visible. We report a mean velocity shift of 0.13 km s⁻¹ across the whole source. The largest velocity shift values are found on the west side of the protostellar core, where the infalling material from the envelope reaches the core, under the assumption that the small velocity gradient seen along the filament represents an accretion flow. On the southeast and northwest sides of the core there is a very low-velocity shift, equal to ~0.06 km s⁻¹ and these values are increasing toward the protostar position up to 0.10 km s⁻¹. To calculate these low-velocity shifts we use the pixels indicated with a black plus sign in Fig. 8.

Fig. 7 PPV image of C180 (2–1) and DCO+ (3–2) gas. Each data point denotes the location and centroid velocity of a Gaussian component and each color refers to a different Gaussian fit. The centroid velocity is shown for the DCO+ (3–2) line (in orange), and the green data points represent the brightest component of the C180 (2–1) line. The red and blue data points indicate the two lower intensity velocity components. The black circle represents the position of the protostar.

5 Discussion

Star-forming regions can be more completely understood by analyzing the distribution of different molecular species in the velocity. We compare the velocity shift between these two tracers in the cloud. The velocity shift between a neutral and an ionized species was observed in the past. For instance, Henshaw et al. (2013) studied the large-scale velocity field throughout the cloud in G035.39-00.33, and found a velocity shift between the two tracers of N₂H⁺(l–0) and C¹⁸O (1–0), in agreement with a model of collision between filaments that is still ongoing. It follows that the velocity structure of the core does not have intrinsic properties but is a product of large-scale motions on filamentary scales. They proposed that the velocity difference in the cloud occurs because of filament merging, implying that higher velocity filaments are interacting with a lower velocity, less massive filament, increasing the density of an intermediate velocity filament (Barnes et al. 2018).

Another scenario is that the velocity shift between the C¹⁸O (2–1) and DCO⁺ (3–2) reveals relative motions between the dense gas, traced by DCO⁺ (3–2), and the surrounding less dense envelope, traced by C¹⁸O (2–1). DCO⁺ (3–2) is simply tracing higher densities, which may not necessarily have the same velocities as the gas traced by C¹⁸0 (2–1), especially because the C¹⁸O (2–1) emission is much more extended, so the kinematics derived from the line profile is affected by lower density material not seen in DCO⁺ (3–2). According to Zhang et al. (2017), velocity shifts between high-density and low-density tracers within a cloud are indicative of gas expanding and contracting. It is based on the assumption that the higher critical density molecules trace the dense gas closer to the inner of a core, while lower critical density molecules trace the more extended gas in the outer envelope.

Magnetic braking might have had a great impact in this cloud. We determined that the rotation axis of the core and magnetic field lines are almost aligned. According to the mag-netohydrodynamic (MHD) collapse models, magnetic braking should be effective in this cloud, which is in agreement with the absence of a resolved Keplerian disk. Magnetic braking is an effective way to remove angular momentum from infalling and rotating material, suppressing envelope fragmentation and the formation of large disks.

Fig. 8 Map of velocity shift between C180 (2–1) and DCO+ (3–2). Overlaid in black contours is the H2 column density (levels: [1.0,1.5,2.0] × l022 cm−2). The black star represents the position of the protostar and the black plus signs are positions of the lowest velocity shifts, which are used for the mean value.

6 Conclusions

We studied the region around the young low-mass Class 0 source IRAS 15398 using APEX. The kinematic analysis performed with the C¹⁸0 (2–1) emission line as a low-density material tracer (extended gas) and with the DCO⁺ (3–2) line as a tracer of dense gas closer to the protostar. The measured kinematics parameters revealed several properties by performing Gaussian fitting, using the PYSPECKIT package. Our main conclusions can be summarized as follows:

• From the spectral line profiles we conclude that the two velocity components of C¹⁸0 (2–1) in the west side of the region have merged together toward the position of the protostar. C¹⁸0 (2–1) is a lower-density tracer than DCO⁺ (3–2); the fainter velocity component seen only in the former tracer is likely low-density filamentary material associated with the cloud. The core envelope appears to be located in correspondence of the merger of these filamentary structures.

• We see a velocity gradient along the filament in the DCO⁺ (3–2) gas. Therefore, we measured the ongoing accretion material toward the protostar core in this gas, Ṁ_acc = 9.7 10⁻⁷ M_⊙ yr⁻¹, where the accretion rate is expected to be accurate within a factor of 2.

• The mean velocity gradient is roughly 5.1 km s⁻¹ pc⁻¹ measured in DCO⁺ (3–2) around the protostar core, which is linked to the rotation of the core. This velocity gradient at the position of the protostar is in the east-west direction, oriented approximately perpendicular to the bipolar outflow previously found.

• Line widths of DCO⁺ (3–2) increase toward the position of the protostar, probably due to protostellar feedback.

• We observed a velocity shift between neutral and ionized species. A higher velocity is always present in the C¹⁸O (2–1) data compared to the DCO⁺ (3–2) data. The mean velocity difference, V_lsr(C¹⁸O) – V_lsr(DCO⁺), is equal to 0.13 km s⁻¹ across the full filament. This is consistent with a model of collision between filaments that is still ongoing. The velocity shift between the C¹⁸O (2–1) and DCO⁺ (3–2) illustrates the relative motion of the dense gas, traced by DCO⁺ (3–2), and the surrounding less dense envelope, traced by C¹⁸O (2–1).

Further observational investigations are needed to determine in more detail the connections within the kinematics and magnetic field in this source.

Acknowledgements

E.R. acknowledges the support from the Minerva Fast Track Program of the Max Planck Society. The authors would like to thank Jaime Pineda Fornerod for his support and discussion about the code to calculate the velocity gradient. This research has made use of data from the Herschel Gould Belt survey (HGBS) project (http://gouldbelt-herschel.cea.fr). The HGBS is an Herschel Key Programme jointly carried out by SPIRE Specialist Astronomy Group 3 (SAG 3), scientists of several institutes in the PACS Consortium (CEA Saclay, INAF-IFSI Rome, and INAF-Arcetri, KU Leuven, MPIA Heidelberg), and scientists of the Herschel Science Center (HSC; André et al. 2010).

Appendix A Positions of spectral grid

Figure A.1 shows the position of each spectra in Fig. 3. There are 40 black dots that show the exact location of each spectrum at an 18 arcsec interval from each other.

Fig. A.1 Map of centroid velocity od DCO+ (3–2). The black dots showing the position of each spectrum. Overlaid in black contours is the H2 column density (levels:[1.0, 1.5, 2.0] 1022 cm−2).

References

All Tables

All Figures

Fig. 1 Integrated intensity of C18O (1–2) (top) and DCO+ (3–2) (bottom) toward IRAS 15398. The white contour levels are 10, 20, and 30 times of mean rms value for the DCO+ line and 40, 50, and 60 times of mean rms value for the C18O (2–1) line. The beam size is shown in the bottom left corner. The yellow contours in both images show H2 column density (levels:[1.0, 1.5, 2.0] 1022 cm−2). The black star gives the position of the protostar.

In the text

	Fig. 2 Channel maps of C18O (2–1) emission. The black cross indicates the position of the protostar core. The velocity of each channel is shown above each panel.
In the text

	Fig. 3 Spectrum grid of C18O (2–1) (black histogram) and DCO+ (3–2) (red histogram) for 40 random positions around the core. Shown are the two-Gaussian fits (blue curves) and one-Gaussian fits (green curves) to the C18O (2–1) line. The vertical line is Vlsr = 5.2 km s−1. The blue star shows the spectrum for the position of the central protostar.
In the text

Fig. 4 Maps of kinematic parameters. Top: gas velocity dispersion of the cloud traced using the DCO+ (3-2). Bottom: centroid velocity map obtained fitting the observed DCO+ (3-2). The black vectors represent the polarization angles, tilted by 90 degrees to trace the magnetic field direction (Redaelli et al. 2019a). The blue and red arrows show the direction of the outflow (PA = 35°, from Bjerkeli et al. 2016), and the yellow arrow presents the mean velocity gradient direction around the core. The star shows the position of the protostar. The contours represent N(H2) column density levels, as derived from Herschel data: [1.0, 1.5, 2.0] 1022 cm−2.

In the text

	Fig. 5 Centroid velocity map of the DCO+ (3–2) line overlaid with the gradient arrows (only vectors with S/N > 3 are shown). The arrow length represents the relative vector magnitude of the gradient, according to the scale shown in the top left corner, and the direction of the arrows points to the steepest velocity field change.
In the text

Fig. 6 Maps of kinematic parameters. Left panel: bright component of centroid velocity map of C18O (2–1). Right panel: centroid velocity map only showing the weak component of the Gaussian fitting. The black vectors represent the polarization angles, tilted by 90 degrees. The star shows the position of the protostar. The contours represent N(H2) column density levels, as derived from Herschel data: [1.0,1.5,2.0] × 1022 cm−2.

In the text

Fig. 7 PPV image of C180 (2–1) and DCO+ (3–2) gas. Each data point denotes the location and centroid velocity of a Gaussian component and each color refers to a different Gaussian fit. The centroid velocity is shown for the DCO+ (3–2) line (in orange), and the green data points represent the brightest component of the C180 (2–1) line. The red and blue data points indicate the two lower intensity velocity components. The black circle represents the position of the protostar.

In the text

	Fig. 8 Map of velocity shift between C180 (2–1) and DCO+ (3–2). Overlaid in black contours is the H2 column density (levels: [1.0,1.5,2.0] × l022 cm−2). The black star represents the position of the protostar and the black plus signs are positions of the lowest velocity shifts, which are used for the mean value.
In the text

	Fig. A.1 Map of centroid velocity od DCO+ (3–2). The black dots showing the position of each spectrum. Overlaid in black contours is the H2 column density (levels:[1.0, 1.5, 2.0] 1022 cm−2).
In the text