© The Authors 2024

1 Introduction

Stars like our Sun form in the cold and dense parts of molecular clouds, which consist of filaments, clumps, and dense cores (Shu et al. 1987; Lada & Lada 2003; Kennicutt & Evans 2012). The formation of stars is regulated not only by gravity, but also by magnetic fields (B-fields), turbulence, and radiative and mechanical feedback from nearby stars (McKee & Ostriker 2007; Hennebelle & Falgarone 2012; Hennebelle & Inutsuka 2019). Turbulence counteracts gravity and prevents the collapse of the dense core, but it can also generate shocks that compress the gas and trigger the formation of new stars.

B-fields are thought to play a critical role in the star formation process by supporting the material in the cloud against self-gravitational collapse (see, e.g., Mouschovias 2001; Crutcher 2012; Pattle & Fissel 2019; Pattle et al. 2023). The collapse can occur when the dense core is magnetically super-critical, with the ratio of mass to magnetic flux (known as mass-to-flux ratio) greater than 1; that is, when B-fields are not strong enough to prevent the gravitational collapse. The collapse of super-critical cores will lead to the formation of protostars. Conversely, when the core is sub-critical, strong B-fields can act against gravitational collapse. However, several mechanisms might enhance the mass-to-flux ratio and cause the core to collapse, including ambipolar diffusion (see e.g., Mestel & Spitzer 1956; Mouschovias et al. 2006; Liu et al. 2022) and magnetic reconnection diffusion (Lazarian & Vishniac 1999; Santos-Lima et al. 2012).

Feedback from massive stars (O- or B-type) can influence the properties of the parental molecular cloud and the star formation processes (e.g., Krumholz et al. 2014; Pabst et al. 2019, 2020). Three-dimensional magnetohydrodynamic (3D MHD) simulations have suggested that B-fields also influence the structure of clouds along with radiation pressure feedback from massive stars (e.g., Henney et al. 2009; Mackey & Lim 2011a). It has been shown that when the B-field orientations are parallel to the ionizing radiation field, the cloud will tend to have a flattened shape, while the cloud will have a column-like or pillar structure if the B-fields are perpendicular to the radiation field (Mackey & Lim 2011b). These findings from simulations have been in agreement with observations and could explain various globular and elongated substructures in the IC 1396 region (Soam et al. 2018), the famous pillar structure in the M16 region (Pattle et al. 2018), and the Horsehead Nebula (Hwang et al. 2023).

Interstellar dust grains are known to have aspherical shapes and can be systematically aligned with a fixed direction in space, for instance, interstellar B-fields (Hall 1949; Hiltner 1949). The alignment of interstellar dust grains with B-fields involves the internal alignment of the grain axis of the maximum inertia moment (also the shortest axis) with its angular momentum (i.e., internal alignment) and the alignment of the angular momentum with B-fields (external alignment). Internal alignment is caused by internal relaxation processes, including Barnett relaxation and inelastic relaxation (Purcell 1979). The external alignment is caused by the radiative torque (RAT) alignment mechanism (Draine & Weingartner 1997; Lazarian & Hoang 2007) and by the more comprehensive magnetically enhanced RAT mechanism (MRAT, Hoang & Lazarian 2016). The internal and external alignment processes establish the alignment of the grain with its minor axis parallel to the B-fields (see, e.g., Hoang et al. 2022). Dust grains absorb or scatter starlight in the ultraviolet (UV), optical, or near-infrared (NIR) range. However, aligned dust grains tend to absorb and scatter the electric field component of the starlight along their major axis more efficiently. This results in the polarization of background starlight with the polarization vector parallel to the grain’s minor axis. Thus, the polarization of the starlight seen in the UV to optical or NIR wavelength has the polarization orientation parallel to the B-field orientation. Aligned dust grains also re-emit the energy absorbed from starlight in the longer wavelength from far-infrared (FIR) to sub-millimeter (sub-mm), the so-called thermal dust emission. As the electric field component of starlight is absorbed more along the grain’s major axis, the dust grain will re-emit more energy with the electric field component along this major axis than along the minor axis. This causes the thermal emission of aligned dust grains to be polarized with the polarization orientation along the grain’s major axis. Because dust grains have their minor axes aligned with the B-field orientation, the polarization orientation of thermal emission is perpendicular to the B-field orientation. Thus, if the position angle of the polarized thermal dust emission is known, one can rotate the polarization angle by 90° to infer the position angle of the B-fields in the medium. Linear polarized thermal dust emission has been routinely used to trace B-fields in various environments of molecular clouds in the Milky Way (see, e.g., Ward-Thompson et al. 2017; Chuss et al. 2019; Pattle et al. 2022) and extragalactic regions (see, e.g., Lopez-Rodriguez et al. 2021, 2022; Borlaff et al. 2023). The Davis-Chandrasekhar-Fermi method (DCF; Davis 1951; Chandrasekhar & Fermi 1953) has been widely used to measure the strengths of the B-field projected on the plane of the sky, B_pos, based on the assumption that the kinetic energy of turbulent gas balances to that of turbulent B-fields.

Adopting the mean value of B-field strengths along the line of sight (B_los) measured using the Zeeman technique over different molecular clouds, Crutcher et al. (2010) derived a relationship between the B-field strengths and gas volume densities (n_h). For low densities of n_H ≤ 300 cm⁻³, B-fields were found not to considerably change with densities. However, for high densities above n_H ~ 300 cm⁻³, the B-field strengths increase with densities as $n_{\mathrm{H}}^k$, where the power-law index k is ~0.65. Pattle et al. (2023) compiled a large number of measurements of the B-field strength projected on the plane of the sky (B_pos) estimated using the DCF method. These authors showed that these B_pos generally follow the results from Crutcher et al. (2010) until n_H ≥ 10⁷ cm⁻³. Beyond this n_H value, the B_pos distribution is more sparse (see Fig. 2a in Pattle et al. 2023). Theoretical studies have suggested that the power-law index varies between k ~ 2/3 for weak B-fields (Mestel 1966) and k ≲ 0.5 for strong B-fields (Mouschovias & Ciolek 1999). Therefore, the maps of the morphology and strength of the B-fields within a specific cloud are essential to understanding the role of the local B-fields on star formation within the cloud. Previous authors have attempted to derive the map of B-fields toward some well-known star-forming regions, for instance, OMC-1 (Guerra et al. 2021; Hwang et al. 2021), 30 Doradus (30 Dor, Tram et al. 2023), and the G11.11–0.12 filament (Ngoc et al. 2023). These studies revealed significant variations in the B-field morphology and strength across the clouds, demonstrating that B-fields strongly influence the cloud’s evolution and, thus, star formation activities.

In this work, we focus on a specific system of the ρ Ophiuchus-A (hereafter ρ Oph-A) molecular cloud to study the effects of stellar feedback and B-fields on the evolution of the cloud and star formation. The ρ Oph-A molecular cloud is located in the northern part of the L1688 active low-mass star-forming region in the ρ Ophiuchus dark cloud complex (Loren et al. 1990; Wilking et al. 2008; Esplin & Luhman 2020). The ρ Oph-A cloud was initially identified as a dense gas clump based on DCO⁺ observations performed by Loren et al. (1990), along with other clumps named ρ Oph-B to ρ Oph-F. This region has high extinctions (A_V ≳ 20 mag) corresponding to the high gas column densities compared to other regions, which are typically more quiescent (see Fig. 1). Being one of the closest low-mass star-forming regions in our Galaxy, located at a distance of ~137 pc (Ortiz-León et al. 2017), the gas and dust content in the Oph–A region has been studied extensively at almost every wavelength by both ground-based and space instruments for many years (see e.g., Wilking et al. 2008, for a review). The ρ Oph-A region hosts numerous prestellar and protostellar cores identified by the Ground Belt Survey (Ladjelate et al. 2020) and protostars in the early stages (i.e., Class 0/I Young Stellar Objects (YSOs), Enoch et al. 2009; Evans et al. 2009; Connelley & Greene 2010). The ρ Oph-A cloud is mostly influenced by a nearby embedded early high-mass B3/4 star (Brown & Zuckerman 1975; Mookerjea et al. 2018), Oph-S1, which has a mass of ~8 M_⊙ (Hamaguchi et al. 2003). The Oph-S1 star is surrounded by a compact H II region (Andre et al. 1988) and is believed to trigger the star formation activity in Oph-A, especially along the ridge of the dense cloud with indications of protostellar clumps (see Fig. 2a in Motte et al. 1998). The star formation activity in the ρ Oph-A is also believed to be affected by the Sco OB2 association located in the western part of the cloud at a distance of ~145 pc (de Zeeuw et al. 1999).

We investigated B-fields in ρ Oph-A using data observed in various wavelengths: NIR (JHK bands) using imaging polarimeter SIRPOL (Kandori et al. 2006) mounted on the 1.4-m Infrared Survey Facility (IRSF) telescope (Kwon et al. 2015), FIR by the High-resolution Airborne Wideband Camera-plus (HAWC+; Harper et al. 2018) on board the 2.7-m Stratospheric Observatory For Infrared Astronomy (SOFIA) (Santos et al. 2019), and sub-mm using the linear polarimeter POL-2 placed in front of the Submillimetre Common-User Bolometer Array 2 (SCUBA-2; Holland et al. 2013) camera mounted on the 15-m James Clerk Maxwell Telescope (JCMT) (Kwon et al. 2018). Santos et al. (2019) performed the analysis on polarimetric data from SOFIA/HAWC+ observations at 89 and 154 μm to study the change in the slope of the polarization spectra with respect to the column densities and dust temperatures in ρ Oph-A. Later, Tram et al. (2021) used a similar SOFIA/HAWC+ dataset to investigate the relation between the polarization fraction of thermal dust emission versus dust temperatures in the cloud and to understand the mechanism of dust grain alignment and disruption induced by RATs. Kwon et al. (2018) estimated the B_pos values in ρ Oph-A, using the JCMT/POL-2 data at 850 μm, toward a few of the densest parts of the cloud, yielding B_pos values in the range of ~0.5–5 mG. The highest strength was measured toward the densest part of the cloud associated with the starless core SM1 (see Fig. 1). These B_pos strengths in ρ Oph-A are comparable to values found in the OmC-1 region (Hwang et al. 2021). Kwon et al. (2018) also found that the B-field direction is generally perpendicular to the main structure of the cloud, which is generally seen in the dense cores or high-density structures of molecular clouds. The ρ Oph-A region provides an ideal isolated system for studying the impact of the B-fields in the star-forming region and investigating the interplay between B-fields and the stellar feedback into the molecular cloud. In this paper, we aim to: 1) determine the B-field morphology and strengths in the entire ρ Oph-A region using the SOFIA/HAWC+ dust polarization data; 2) investigate the role of B-fields relative to gravity and turbulence in the star formation activity in the cloud; and 3) quantify the combined effect of B-fields and radiation pressure from the nearby Oph-S1 star on the shape of the ρ Oph-A cloud.

The paper is organized as follows. We describe a set of archival observational data in Sec. 2. We present the results of mapping the B-field morphology and strength using various techniques in Sec. 3 and discuss the roles of B-fields relative to gravity, turbulence, and radiation pressure feedback from the Oph-S1 star in Sec. 4. Finally, we summarize our results in Sec. 5.

Fig. 1 Distribution of dust temperatures (colors) and H2 column densities (contours) in ρ Oph-A (Santos et al. 2019). Contour levels of $N_{{\mathrm{H}}_2}$ are shown at (1, 2, 4, 6, 8, and 10) × 1022 cm−2. The star symbol shows the position of Oph-S1. Small circles show the positions of YSOs (Enoch et al. 2009; Evans et al. 2009; Connelley & Greene 2010). Yellow squares and cyan diamonds indicate the positions of prestellar and protostellar cores, respectively, identified as part of the Herschel Gould Belt Survey (Ladjelate et al. 2020). The “×” symbol shows the position of the SM1 core.

2 Observational data

2.1 Dust polarization data from SOFIA/HAWC+

We use the FIR polarimetric observations from the onboard SOFIA/HAWC+ toward ρ Oph-A. The observations were done on 2017 May 17, under the proposal ID 70-0511 (PI: Darren Dowell), using the nod-match-chop mode, with 154° chop angle, 480″ chop throw, and 4-points dithers. Here, we use the SOFIA/HAWC+ Level-4 data for two bands C and D, centered at 89 and 154 μm, respectively, which is fully calibrated and can be directly downloaded from the SOFIA archive¹. The final polarimetric maps cover an area of ~4.5′ × 4.5′ (corresponding to $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}18\times 0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}18{\mathrm{p}\mathrm{c}}^2$) toward the central region in ρ Oph-A with the angular resolutions of 7.8″ and 13.6″ full-width-half-maximum (FWHM) at 89 and 154 μm, respectively. The pixel size of HAWC+ maps is ~2.0″ at 89 μm and 3.4″ at 154 μm. The observations and data reduction are described in detail in Santos et al. (2019).

From the maps of the Stokes-I, Stokes-Q, and Stokes-U (along with their uncertainties σI, σ_Q, σ_U) the maps of polarization fraction (p₀) and its corresponding uncertainty (σ_p) were calculated. The de-biased polarization fraction, p, is calculated as: $p^2=p_0^2-{\sigma }_{\mathrm{p}}^2$ (Vaillancourt 2006). The maps of polarization angle (θ) and its corresponding uncertainty (σ_θ) in degrees can be derived from Stokes-Q and Stokes-U. We refer to Gordon et al. (2018) for a full description of estimates of the polarization fraction, polarization angle, and their corresponding uncertainties. To enhance the quality in measuring the maps of B-field orientation and strength in ρ Oph-A, we apply some critical conditions in choosing vectors (pixels) in SOFIA/HAWC+ maps. We mask out from SOFIA/HAWC+ maps the pixels for which I/σ_I < 54 for band C and I/σ_I < 214 for band D. In both bands, we further removed the pixels for which we have p/σ_p < 3 and p>50% from the HAWC+ maps (see Appendix A).

2.2 Spectroscopic data

We used the spectroscopic observations of HCO⁺(4–3) line at ~356.734 GHz obtained from the Heterodyne Array Receiver Program (HARP; Buckle et al. 2009) mounted on JCMT. The observations were carried out during June and July 2011 (Proposal ID: M11AU13), covering the ρ Oph-A region with a beam size (or FWHM) of ~14″. We used the reduced spectral map of HCO⁺(4–3) line directly downloaded from the JCMT science archive². The data cube had an original spectral resolution of 0.025 km s⁻¹ and was smoothed to have a final spectral resolution of 0.5 km s⁻¹. For comparison to the HCO⁺(4–3) data, we also used observations of N₂H⁺(3–2) line at ~279.515 GHz obtained from the APEX-2 instrument on the 12-m Atacama Pathfinder Experiment telescope (APEX; Güsten et al. 2006). The data cube has an original spectral resolution of 0.082 km s⁻¹. Here, we used the back-end data cube with a spectral resolution of 0.25 km s⁻¹, covering the ρ Oph-A region with the beam size of ~22″ (for details, see Liseau et al. 2015).

Fig. 2 Maps of the B-field orientation (white segments) on top of the total intensity (Stokes-I, color maps) toward ρ Oph-A, observed with SOFIA/HAWC+ in band C (left) and band D (right), centered at 89 and 154 μm, respectively. The gray contours in each map show the Stokes-I intensity at the corresponding wavelength, with contour levels at 0.12, 0.2, 0.4, 0.6, and 0.8 Jy/arcsec2. The star symbol marks the position of Oph-S1 star. Small white circle symbols indicate the positions of YSOs. The filled orange circles show the beam size of each color map.

2.3 Dust temperature and column density maps

We adopted the maps of dust temperature (T_dust) and gas column density ($N_{{\mathrm{H}}_2}$) from Santos et al. (2019). These maps were obtained from the fitting of a modified black body to the spectral energy distribution (SED) using the continuum flux maps from Herschel/PACS (Poglitsch et al. 2010) at 70, 100, and 160 μm, with an FWHM of 11.4″ (for details, see Sect. 3.2 in Santos et al. 2019). Figure 1 shows the spatial distribution T_dust toward ρ Oph-A, overlaid with contours of column density, $N_{{\mathrm{H}}_2}$. In this study, we use the map of T_dust as the map of the gas kinetic temperature, assuming that the gas is in the local thermodynamic equilibrium (LTE) and well-coupled with the dust.

3 Magnetic fields in ρ Oph-A

3.1 Magnetic field morphology

We construct the B-field morphology in ρ Oph-A by rotating polarization angles by 90°. Figure 2 shows the maps of B-field orientation (half-vectors) inferred from maps of the polarization angle at 89 and 154 μm. The background color maps show the Stokes-I intensity at the respective wavelength. Here, the position angle of B-fields (θ_B) is measured east of north, and in the range between 0° and 180°. The polarimetric data at both wavelengths can well reveal B-fields toward the dense region in ρ Oph-A, where Stokes-I>0.12 Jy/arcsec² (the area inside the lowest contour in each map). The data at 89 μm can further resolve the B-fields morphology extended toward the warmer dust region located at the east side of the Oph-S1 star.

The B-field vectors in ρ Oph-A are generally well-ordered and largely follow the material structure in low- to intermediate-densities regimes (see Fig. 2). B-field orientations seem to be bent toward the higher density region in the center of the map and are mostly perpendicular to the ridge of the cloud. On the eastern side of the map, the orientation of B-field vectors tends to vary in the counter-clockwise direction from north to south; thus, the position angle of B-fields increases. Particularly, in the northeast part of the map, B-field vectors run along the northeast to southwest direction with a median θ_B value of ~56° and ~49° at 89 and 154 μm, respectively. This feature is also seen in a small area in the northwest part of the map. Toward the east side, B-field vectors orient nearly to the east-west direction, with median θ_B values of ~83° at 89 μm and ~75° at 154 μm. On the other hand, B-field vectors in the southeast region tend to orient consistently to the southeast to northwest directions, with median θ_B values of ~92° and ~85° at 89 and 154 μm, respectively. Interestingly, the data at 154 μm shows B-field vectors running nearly vertically, that is, along the north-south direction toward a small area at the north part of the map, which is not seen in the data at 89 μm (see Fig. 2 and also Fig. B.1).

In addition, there are two outstanding components of B-fields toward the regions characterized by the highest H₂ column densities: one component associated with the peak of continuum emission at 154 μm runs along the northeast to southwest direction with the median B-field orientation of ~70°; the other is located at the north of the first component and associated with the peak of continuum emission at 89 μm, which exhibits the B-field vectors along the east-west direction with a median position angle of ~88° and −81° at 89 and 154 μm, respectively (see Fig. B.1). We will further discuss the B-field morphology together with B-field strengths in Sec. 4.1.

3.2 Measurement of magnetic field strength using the DCF method

We first used the classical DCF method (Davis 1951; Chandrasekhar & Fermi 1953) (hereafter we refer to as the DCF method) to calculate the strengths of the B-field projected in the plane-of-the-sky ($B_{\mathrm{pos}}^{\mathrm{DCF}}$) in the entire region of the ρ Oph-A cloud. We compared our obtained results with those obtained using the more advanced techniques described in Sec. 3.3.

According to the DCF method, the B-field strength can be estimated by Ostriker et al. (2001): $$B_{\mathrm{p}\mathrm{o}\mathrm{s}}^{\mathrm{D}\mathrm{C}\mathrm{F}}=f\sqrt{4\pi \rho }\frac{{\sigma }_{\mathrm{N}\mathrm{T}}}{\delta \theta },$$(1)

where f is the correction factor, $\rho =\mu {\mathrm{m}}_{\mathrm{H}}{n}_{{\mathrm{H}}_2}$ is the gas mass density in units of g cm⁻³, μ=2.8 is the mean molecular weight per hydrogen molecule (Kauffmann et al. 2008), and m_H is the mass of the hydrogen atom, $n_{{\mathrm{H}}_2}$ is the gas volume density in units of cm⁻³. Then, σ_NT is the one-dimensional (1D) nonthermal velocity dispersion along the line of sight (LOS) of the ideal gas species that traces the dense region, similarly to the dust-polarized emission in units of km s⁻¹, and δθ is the polarization angle dispersion in degrees. The correction factor, f, is uncertain, and numerical simulations by Ostriker et al. (2001) show f = 0.5 when the polarization angle dispersion does not exceed 25°. However, varying values of f will result in different B_pos strengths (see e.g., Skalidis et al. 2021; Liu et al. 2021, for more detailed discussions).

Here, we use the formula given in Crutcher (2004): $$B_{\mathrm{p}\mathrm{o}\mathrm{s}}\approx 9.3\sqrt{n_{{\mathrm{H}}_2}}\frac{\mathrm{\Delta }{v}_{\mathrm{N}\mathrm{T}}}{\delta \theta }(\mu \mathrm{G}),$$(2)

where $\mathrm{\Delta }{v}_{\mathrm{NT}}={\sigma }_{\mathrm{NT}}\sqrt{8\ln 2}$ is non-thermal FWHM velocity dispersion measured in units of km s⁻¹, and the correction factor f = 0.5 has been used. This formula is widely used to estimate B_pos via the DCF technique.

3.2.1 Polarization angle dispersion maps

It is essential to see how the B-field orientation changes locally within the molecular cloud. To do so, we calculated δθ across the ρ Oph-A region. We used the method presented in Hwang et al. (2021), for which δθ is calculated pixel-by-pixel throughout the map of ρ Oph-A. We summarize the procedure of calculating δθ toward the pixel i^th in the map as follows:

1. We used a box covering the pixel i^th with a size of 9×9 pixels, corresponding to two beam sizes of the SOFIA/HAWC+ observations (i.e., 15.6″ and 36.4″ in bands C and D, respectively).

2. We calculated the mean value of polarization angles within the small box, which implies the local mean polarization angle toward the i^th pixel.

3. For each pixel within the small box, we calculated the difference between the observed polarization angle and the mean value. We computed the root mean square (RMS) of the distribution of angle differences in the box. This RMS value indicates the polarization angle dispersion calculated toward the i^th pixel.

4. The uncertainty of the polarization angle dispersion toward the pixel i^th is propagated from the uncertainty of polarization angles in the small box and determined as the mean value of polarization angle uncertainties within the small box.

By repeating the process for all pixels in the entire region, we obtained the map of δθ in ρ Oph-A. We note that if there are more than 50% of pixels undefined (i.e., having NaN value of θ) in the computing box, the result of polarization angle dispersion calculation toward pixel i^th is ignored. We applied this criterion to keep the high quality of the δθ estimate. It is also necessary to de-bias the polarization angle dispersion due to the uncertainty of the observed polarized angles.

Figure 3 shows the maps of polarization angle dispersion toward ρ Oph-A using HAWC+ data in two bands C (left panel) and D (right panel), respectively. The polarization angle dispersion calculated using data in bands C and D cover ranges of 2.3°−32.2° and 1.3°−31.5°, which correspond to the median values of 5.7° and 5.8°, respectively. The differences of δθ in two wavelengths imply that the data at the shorter wavelength tends to trace warm dust surrounding the high-mass star S1 and toward the outer layers; whereas the data at the longer wavelength probes the colder dust in the denser layers of the cloud. We note, however, that in maps of δθ in both bands, there are a few pixels where δθ > 25° (27 pixels in band C and 45 pixels in band D). Thus, we do not use the δθ value in these pixels to estimate B_pos in Sect. 3.2.4.

3.2.2 Velocity dispersion map

We used the HCO⁺(4–3) molecular data observed by JCMT/HARP to estimate the total velocity dispersion, σ_υ. This molecular line has a critical density of ~10⁶ cm⁻³ (Greve et al. 2009) and thus is considered a good tracer for probing the dense region in the ρ Oph-A cloud. The HCO⁺(4–3) emission is well-represented by a single Gaussian profile. Thus, in each pixel, we fit the line with a single Gaussian. The central velocity (υ₀) across ρ Oph-A covers the range from 2.25 to 3.76 km s⁻¹ with a median value of 3.12 km s⁻¹, which is in agreement with the systemic velocity of ρ Oph-A of ~3.44 km s⁻¹ (André et al. 2007). We used the line width of the Gaussian as the measure of total velocity dispersion toward each pixel.

To construct the map of the 1D non-thermal (turbulent) velocity dispersion (σ_NT), we subtracted the thermal component (σ_υ,th) from the map of the total velocity dispersion, given by: $${\sigma }_{\mathrm{N}\mathrm{T}}=\sqrt{{\sigma }_v^2-{\sigma }_{v,\mathrm{t}\mathrm{h}}^2}=\sqrt{{\sigma }_v^2-\frac{k_{\mathrm{B}}{T}_{\mathrm{g}\mathrm{a}\mathrm{s}}}{m_{\mathrm{m}\mathrm{o}\mathrm{l}}}}\left(\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}\right),$$(3)

where m_mol=29 amu is the mass of the HCO⁺ molecule, k_B is the Boltzmann constant, and T_gas is the kinetic temperature of the gas. Here, we assume that the gas is in LTE condition and well-coupled with the dust in the cloud. Thus, the gas kinetic temperature is equal to the excitation temperature in all levels (T_ex) and equal to the dust temperature: T_gas = T_ex = T_dust. We convolved the map of T_dust to the lower resolution of the molecular line map (14″) to obtain the same pixel size on both maps. We converted the velocity dispersion of the non-thermal component into the FWHM velocity dispersion, given as $\mathrm{\Delta }{v}_{\mathrm{NT}}=$ ${\sigma }_{\mathrm{NT}}\sqrt{8\ln 2}$. The uncertainty of non-thermal FWHMs is propagated from the uncertainty of the line width obtained from the Gaussian fit.

Figure 4 shows the map of the 1D non-thermal FWHM of the HCO⁺ molecular line toward ρ Oph-A. The συ_NT values are in the range from 0.62 to 2.98 km s⁻¹ with a median value of 1.19 km s⁻¹. We find a good agreement between our result obtained with HCO⁺(4–3) and those from other line tracers, namely, 0.95 km s⁻¹ using N₂H⁺(3–2) from APEX and 1.5 km s⁻¹ using C¹⁸O(3–2) from JCMT/HARP (see Appendix C for details). Thus, it is plausible to use results from the HCO⁺(4–3) line in the following analysis.

Fig. 3 Polarization angle dispersion maps using data in band C (left panel) and band D (right panel). The black contours are from the total intensity map at the corresponding wavelength, with the contour levels at 0.12, 0.2, 0.4, 0.6, and 0.8 Jy/arcsec2. The star symbol in each map indicates the position of the Oph-S1 star.

Fig. 4 Distribution of the non-thermal FWHM of the HCO+(4–3) molecular line toward ρ Oph-A. The star symbol marks the position of the Oph-S1 star. White contours show the continuum emission (Stokes-I) at 154 μm, with contours levels at 0.12, 0.2, 0.4, 0.6, and 0.8 Jy/arcsec2.

3.2.3 Gas volume density map

We estimated the gas volume, $n_{{\mathrm{H}}_2}$, across the dense cloud from $N_{{\mathrm{H}}_2}$. Here, we assume that the cloud has a spherical Plummer-like shape (Liseau et al. 2015; Santos et al. 2019), which is frequently used to describe the globular dense cores (Plummer 1911). We consider a sphere centered at the peak position of $N_{{\mathrm{H}}_2}$, located at (RA, Dec)= $({16}^{\mathrm{h}}{26}^{\mathrm{m}}{27}_{\cdot }^{\mathrm{s}}8,-{24}^{\circ }{24}^{\mathrm{\prime }}{00}_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}0)$, which is also associated with the position of the starless core SM1 (the separation is much less than 1″). We took the radius of the cloud R_cl ≈ 120″ (~0.074 pc, Santos et al. 2019). To describe the effect of stellar feedback from the high-mass Oph-S1 star, we further add a spherical shell with a radius of $R_{{\mathrm{S}}_1}\approx {80}^{\prime \prime }$ surrounding this star. We also assume that the positions of peak column density and the S1 star are located in the same plane perpendicular to the LOS. The top panels in Fig. 5 show the 3D schematic view of the cloud (left) and a cross-section on the xz-plane at a given y position (right). The x- and y-axes perpendicular to each other are in the plane of the sky and the z-axis is along the LOS. The center of the sphere is located at (x₀, y₀, z₀) = (0, 0, 0). The depth of the cloud at each position (x, y) along z-axis is calculated by: $$\mathcal{W}(x,y)=2(z_1-z_2)=2(\sqrt{R_y^2-x^2}-\sqrt{R_{y,S_1}^2-{(x-x_{{\mathrm{S}}_1})}^2}),$$(4)

where $R_y=\sqrt{R_{\mathrm{cl}}^2-y^2}$ and $R_{y,S1}=\sqrt{R_{\mathrm{S}1}^2-{(y-y_{S1})}^2}$ are the radius of horizontal circles perpendicular to the plane of sky at y position (see top right panel in Fig. 5). We obtain the map of W, shown in the bottom left panel in Fig. 5. The values of W are in the range of ~20″−222″. We then infer $n_{{\mathrm{H}}_2}$ from $N_{{\mathrm{H}}_2}$ using the formula: $$n_{{\mathrm{H}}_2}=\frac{N_{{\mathrm{H}}_2}}{\mathcal{W}}.$$(5)

The bottom right panel in Fig. 5 shows the spatial extent of $n_{{\mathrm{H}}_2}$ across ρ Oph-A, overlaid by the contours representing the continuum emission in band D (Stokes-I). The $n_{{\mathrm{H}}_2}$ map has an angular resolution of 11.4″ and a pixel size of ~3.2″, similar to that of the $N_{{\mathrm{H}}_2}$ map. The gas density spans from ~10⁴⁴ − 10⁶ cm⁻³ and is highest at the center of the cloud. The pattern of $n_{{\mathrm{H}}_2}$ distribution well follows the intensity map in band D, revealing the ridge shape of the dense cloud. Some pixels located at the circular boundary in the $n_{{\mathrm{H}}_2}$ map show anomalously high densities; this is due to the very low depths caused by the spherical model, and thus we exclude these pixels from our B_pos strength measurement.

Fig. 5 Cartoon of 3D structure of the ρ Oph-A cloud where we assume that the cloud has a spherical shape but is pushed by radiation feedback from a nearby high-mass star Oph-S1, shown at the top (left panel, see text) and a plane of the constant coordinate y with the LOS $\hat{z}$ (right panel). Bottom row displays the map of the depth of the cloud, W (left panel), and the map of the gas volume density $\left(n_{{\mathrm{H}}_2}\right)$ inferred from $N_{{\mathrm{H}}_2}$ (right panel). In each panel, the “×” and star symbols mark the peak position of $N_{{\mathrm{H}}_2}$ and location of Oph-S1 star, respectively. Toward maps of W and $n_{{\mathrm{H}}_2}$ in bottom panels, the black dashed circle has a radius Rcl=120″ and black solid contours show the continuum emission at 154 μm with levels similar as shown in Fig. 4.

3.2.4 Magnetic field strength maps

We used the maps of δθ, Δυ_nt, and $n_{{\mathrm{H}}_2}$ to measure $B_{\mathrm{pos}}^{\mathrm{DCF}}$ in the entire region of ρ Oph-A. All input maps were smoothed to the angular resolution of 14″, corresponding to that of the Δυ_nt map. We then reconstructed the maps of B_pos at 89 and 154 μm toward ρ Oph-A using Eq. (2), which are shown in Fig. 6.

The $B_{\mathrm{pos}}^{\mathrm{DCF}}$ values vary from 0.2–2.3 mG and from 0.2–2.5 mG for data at 89 and 154 μm, respectively. The highest strength is found at the center region of the cloud, associated with the emission peak at 154 μm and also the position of the SM1 core. The median $B_{\mathrm{pos}}^{\mathrm{DCF}}$ values are 0.7 ± 0.3 and 0.7 ± 0.4 mG³, as measured using 89 and 154 μm data, respectively. The uncertainty of B_pos is propagated based on the uncertainty of the velocity dispersion of the gas and polarization angle dispersion. Since $n_{{\mathrm{H}}_2}$ map is inferred from $N_{{\mathrm{H}}_2}$ resulted in SED fitting models (see Santos et al. 2019). In addition, the estimated depth of the cloud is based on the 3D geometrical shape of the cloud and also is uncertain. Thus, we did not take into account the uncertainty of the gas volume density in the total uncertainty of the B-field strengths. For this reason, the uncertainty of B_pos might be underestimated.

Fig. 6 Maps of B-field strengths projected in the plane-of-the-sky (Bpos) in ρ Oph-A using data in band C (left panel) and band D (right panel), respectively. The contours in each map are similar as in Fig. 2. The star symbol marks the position of the Oph-S1 star and small circle symbols show the position of YSOs.

3.3 Measurement of magnetic field strength using other methods

Several new techniques have been proposed to improve the DCF method for the last decades. For instance, Hildebrand et al. (2009) introduced the angle dispersion function (ADF) method where the dispersion in polarization angles δθ in Eq. (2) is replaced by the square root of the structure function of polarization angles. This provides an improved measurement of B-field strength but also gives us the ratio of the turbulent to large-scale B-field strengths. Alternatively, Skalidis & Tassis (2021) introduced a modification of the DCF method (hereafter ST), assuming that the mean (ordered) and fluctuation components of the total B-field energy, B₀ and B_t, respectively, are not independent. More recently, the differential measure approach (DMA; Lazarian et al. 2022) is considered promising to enhance the accuracy of B-field strength measurements at scales of ℓ smaller than the turbulence injection scale and takes into account the anisotropic properties of MHD turbulence (given by factor f′). For the sake of cross-checking, the B-field strength measurement in ρ Oph-A using the DCF method, we re-calculated the mean B-field strengths using the ADF, ST, and DMA methods at the two subregions (a) and (e) associated with the densest parts of the cloud (these subregions are named as in Kwon et al. (2018), see also in Fig. B.1).

3.3.1 ADF method

According to the ADF method, the $B_{\mathrm{pos}}^{\mathrm{ADF}}$ can be estimated using Eq. (2) but the polarization angle dispersion δθ is replaced by the structure function of the polarization angles, so-called angular structure function (Hildebrand et al. 2009). The angular structure function, $D_{\theta }^{1/2}$, at an angular scale ℓ is estimated as: $$D_{\theta }^{1/2}(\ell )\equiv {⟨\mathrm{\Delta }{\theta }^2(\ell )⟩}^{1/2}={\{\frac{1}{N(\ell )}{\displaystyle \sum _{i=1}^{N(\ell )}}[\theta (x)-\theta (x+\ell ){]}^2\}}^{1/2},$$(6)

where ⟨...⟩ denotes an average, Δθ(ℓ) = θ(x) − θ(x + ℓ) is the polarization angle difference between individual pairs of polarization vectors at position x and x + ℓ, where θ(x) and θ(x + ℓ) are the corresponding polarization angles at position x and x + ℓ, respectively, and N(ℓ) is number of pairs of vectors with a displacement of ℓ. At scales ℓ much smaller than the scale for a variation of the large-scale B-field structure (typically of ~5′ seen by Planck observations), the total angular dispersion function can be expressed as: $$D_{\theta }(\ell )\simeq b^2+m^2{\ell }^2+{\sigma }_{\mathbf{M}}^2(\ell ),$$(7)

where b represents the turbulent component of B-fields (B_t), mℓ describes the large-scale structure of B-fields (B), and σ_M(ℓ) is subject to measurement uncertainties. Also, the ratio between the turbulent and large-scale B-field strength indicates the angular dispersion, given as (Hildebrand et al. 2009): $$\delta \theta \simeq \frac{{⟨B_{\mathrm{t}}^2⟩}^{1/2}}{B}=\frac{b}{\sqrt{2-b^2}}.$$(8)

When B_t ≫ B, i.e., b ≫ 1 rad, we can obtain $\delta \theta \sim b/\sqrt{2}$.

For each subregion (a) and (e) in ρ Oph-A, we calculate the angular structure function D_θ(ℓ) using HAWC+ data at 154 μm via Eq. (6) and then perform the fit to the quadratic function described in Eq. (7) using the Python curve_fit routine. We select data points that 13.6″ < ℓ ≤ 50″ to avoid fitting to structures higher than the beam size and to better fit to smaller-scale structures. Top panels in Fig. 7 show structure functions of polarization angles as a function of the distance of displacement ℓ toward subregions (a) and (e). From the fitted b value of the fits, we derive δθ values of 10.9° ± 0.2° and 17.5° ± 0.2° toward subregions (a) and (e), respectively. We find these values in good agreement with mean polarization angle dispersion values estimated using our approach described in Sect. 3.2.1, which are ~9.2° and ~14.3° in subregions (a) and (e), respectively. Toward subregions (a) and (e), mean gas volume densities are 2.1 × 10⁵ and 1.1 × 10⁵ cm⁻³ respectively, and their mean non-thermal FWHMs are 1.74 and 1.45 km s⁻¹, respectively. We substitute new δθ values together with the above mean values of $n_{{\mathrm{H}}_2}$ and συ_NT into Eq. (2) to calculate $B_{\mathrm{pos}}^{\mathrm{ADF}}$ strengths in subregions (a) and (e), which are 0.69 mG and 0.26 mG, respectively. These values are in good agreement with $B_{\text{pos}}^{\mathrm{DCF}}$ strengths obtained using the DCF method, where $B_{\text{pos}}^{\mathrm{DCF}}$ values are in ranges of 0.32–1.55 mG and 0.15–0.97 mG toward subregions (a) and (e), respectively, with corresponding mean $B_{\text{pos}}^{\mathrm{DCF}}$ values of ~0.84 and ~0.37 mG.

Fig. 7 Structure function of polarization angles (top panels) and of central velocity (middle panels) measured in the subregions (a) (left) and (e) (right). The bottom panels show the resulting $B_{\text{pos}}^{\text{DMA}}$ strength estimated via Eq. (10) in the two subregions (a) and (e). The HAWC+ data at 154 μm is used to calculate the structure-function (data points). The red solid line in the top panels shows the fitted line of the angular structure-function to the quadratic function, where we chose data points with 13.6″ < ℓ ≤ 50″. The blue dash-dotted line in each panel indicates the upper ℓ-limit of 50″ used for fitting the angular structure-function and for estimating the mean value of $B_{\mathrm{pos}}^{\mathrm{DMA}}$ in the two subregions (a) and (e) (see texts). The blue dashed line in each panel indicates the beam size of HAWC+ observation at 154 μm.

3.3.2 ST method

The ST method is applicable for the compressible B-fields mode and thus is expected to provide a more accurate estimate of B-field strengths (Skalidis et al. 2021). Another advantage of the ST method is that this method accounts for the correlation between the fluctuation and mean components of the total B-field energy not to be zeros. Here, we measure maps of B_pos strength $\left(B_{\mathrm{pos}}^{\mathrm{ST}}\right)$ in band C and band D using the formula as appeared in Skalidis & Tassis (2021): $$B_{\mathrm{p}\mathrm{o}\mathrm{s}}^{\mathrm{S}\mathrm{T}}=\sqrt{4\pi \rho }\frac{\mathrm{\Delta }{v}_{\mathrm{N}\mathrm{T}}}{\sqrt{2\delta \theta }}=\frac{1}{\sqrt{2\delta \theta }}{B}_{\mathrm{p}\mathrm{o}\mathrm{s}}^{\mathrm{D}\mathrm{C}\mathrm{F}}.$$(9)

We show the maps of $B_{\text{pos}}^{\mathrm{ST}}$ across ρ Oph-A in band C and band D in Fig. B.3. We find that $B_{\text{pos}}^{\mathrm{ST}}$ value ranges from 0.03 −0.7 mG at 89 μm and 0.02–1.3 mG at 154 μm, appeared to be lower than $B_{\text{pos}}^{\mathrm{DCF}}$. A median $B_{\mathrm{pos}}^{\mathrm{ST}}$ value of −0.2 ± 0.1 mG is seen in both bands. For a direct comparison with those using other methods, we estimate the mean value of $B_{\mathrm{pos}}^{\mathrm{ST}}$ at 154 μm toward two subregions (a) and (e), which are ~0.19 and ~0.07 mG, respectively. These values are about 4.5 times lower than those estimated using the DCF method.

3.3.3 DMA method

We also calculate B_pos trengths in two subregions (a) and (e) using the DMA method, $B_{\text{pos}}^{\text{DMA}}$, using Eq. (10) following Lazarian et al. (2022): $$B_{\mathrm{p}\mathrm{o}\mathrm{s}}=f^{\mathrm{\prime }}\sqrt{4\pi \rho }\frac{D_v^{1/2}(\ell )}{D_{\theta }^{1/2}(\ell )},$$(10)

where $D_{\theta }^{1/2}(\ell )$ and $D_v^{1/2}(\ell )$ are the second order structure function of polarization angles and central velocity at an angular scale ℓ, f′ is the correction factor depending on the anisotropic properties of MHD turbulence. Here, we adopt f′=1 related to the naive formalism as suggested in Lazarian et al. (2022).

We calculate the structure function of the central velocity of HCO⁺(4–3) line toward subregions (a) and (e), using the similar form in Eq. (6) but for the central velocity υ₀, shown in middle panels in Fig. 7. We then measure $B_{\text{pos}}^{\text{DMA}}$ strengths in two subregions (a) and (e) via Eq. (10), using the information from the $D_{\theta }^{1/2}(\ell )$ together with $D_v^{1/2}(\ell )$ and the mean values of $n_{{\mathrm{H}}_2}$ in the subregions. The $B_{\text{pos}}^{\mathrm{DMA}}$ as a function of displacement ℓ is shown in bottom panels of Fig. 7 for subregions (a) and (e). It can be seen that the resulting $B_{\mathrm{pos}}^{\mathrm{DMA}}$ in the range 13.6″ < ℓ < 50″ is likely flatten within their uncertainty. Thus, in each subregion we take the mean of $B_{\mathrm{pos}}^{\mathrm{DMA}}$ in this range to be the average value, which are ~0.34 mG and ~0.12 mG in subregions (a) and (e), respectively. These values are slightly higher than the resulting B_pos estimated using the ST method, but they are at the low-end of that measured using the DCF method (i.e., ~2–3 times less than the mean values estimated using the DCF method presented in Sect. 3.2). These differences are indeed expected since the DCF method somewhat overestimates the B_pos strength in the cloud. We discuss this issue further in Sect. 4.1.

3.4 Mass-to-flux ratio

The relative importance of B-fields compared to gravity in the cloud can be quantified by the mass-to-flux ratio, M/ϕ, which can be estimated in units of the critical ratio, ${\left(M/{\varphi }_B\right)}_{\text{crit}}=1/(2\pi \sqrt{G})$ (Nakano & Nakamura 1978): λ ≡ (M/ϕ)_observed / (M/ϕ)_crit. We measure λ in the entire region of ρ Oph-A, using the formula introduced by Crutcher (2004): $$\lambda =7.6\times {10}^{-21}\frac{N_{{\mathrm{H}}_2}}{B}\text{, }$$(11)

where $N_{{\mathrm{H}}_2}$ and B are in units of cm⁻² and μG, respectively.

Due to the lack of B_los component measurements in the ρ Oph-A cloud, we adopted the B_pos component as the total B-field strengths B: B = B_tot ≈ B_pos. For the analysis, we employed only the $B_{\mathrm{pos}}^{\mathrm{DCF}}$ map estimated using the DCF method and inferred from the polarimetric data in band D; this is because it covers a larger area of the higher density region in ρ Oph-A than that coming from the data in band C.

The top panel in Fig. 8 shows the spatial distribution of the mass-to-flux ratio in ρ Oph-A. It reveals an increasing gradient of λ from the central and eastern parts of the clump to the outer regions, spanning from 0.03–2.56, with a median value of ~0.26. The region exhibiting the lowest values of λ ≪ 1 is close to the Oph-S1 star, suggesting that the region is firmly magnetically sub-critical, namely, B-fields are strong enough to prevent the gravitational collapse. This can also be explained by the influence of the strong winds and radiation fields from the Oph-S1 star which tend to push the material away to to its vicinity at the west side. The mass-to-flux ratio increases toward the outskirt regions of the Oph-S1 star, especially toward regions along the ridge of the dense cloud (i.e. center, north, and southeast parts of the map) where the gas column density is ≳4 × 10²² cm⁻². These regions are magnetically super-critical and so, the self-gravitational collapse might have occurred there, resulting in several protostars and a single protostellar core identified in the region (Enoch et al. 2009; Evans et al. 2009; Connelley & Greene 2010).

Fig. 8 Maps of λ (top panel), Alfvén Mach number (middle panel), and plasma β parameter (bottom panel) in ρ Oph-A using the $B_{\mathrm{pos}}^{\mathrm{DCF}}$ in band D. In each map, yellow square and cyan diamond symbols indicate the positions of prestellar and protostellar cores, respectively. Green circles indicate the positions of YSOs (Enoch et al. 2009; Evans et al. 2009; Connelley & Greene 2010). The star and “×” symbols mark the positions of the Oph-S1 star and the SM1 core, respectively. Black contours show $N_{{\mathrm{H}}_2}$, similar as in Fig. 1.

3.5 Alfvénic Mach number

To investigate the interplay between B-field and turbulence due to the non-thermal motion of the gas in ρ Oph-A, we calculate the map of Alfvénic Mach number, given as: $${\mathcal{M}}_{\mathrm{A}}=\sqrt{3}\frac{{\sigma }_{\mathrm{N}\mathrm{T}}}{v_{\mathrm{A}}},$$(12)

where σ_NT is given by Eq. (3), and υ_A is the Alfvénic speed driven by the B-fields ($v_{\mathrm{A}}=B/\sqrt{4\pi \rho }$), both in units of km s⁻¹. The factor $\sqrt{3}$ in Eq. (12) assumes that turbulence is isotropic, and allows the calculation of the 3D Alfvénic Mach numbers. We note that adding this factor will provide upper limits to the estimated ℳ_A.

The middle panel in Fig. 8 shows the map of ℳ_A toward ρ Oph-A using results of polarimetric data at 154 μm and HCO⁺(4–3) as the gas tracer. The ℳ_A values are in the range from 0.1–1.5 across the cloud, with a median value of 0.38. Most of ρ Oph-A shows ℳ_A < 1 (magnetically sub-Alfvénic) suggesting the B-field dominates turbulence in the cloud. However, two small regions at the north and southwest of the map exhibit ℳ_A ≥ 1, indicating strong Alfvénic motions of the gas (magnetically super-Alfvénic). These regions are also associated with the positions of a few prestellar and protostellar cores, and YSO candidates (Enoch et al. 2009; Evans et al. 2009; Connelley & Greene 2010).

3.6 Plasma β parameter

To compare the magnetic presure (P_mag) and thermal pressure (P_th) in ρ Oph-A, we constructed the map of the plasma beta parameter (β) given as: $$\beta =\frac{P_{\mathrm{t}\mathrm{h}}}{P_{\mathrm{m}\mathrm{a}\mathrm{g}}}=2{\left(\frac{c_{\mathrm{s}}}{v_{\mathrm{A}}}\right)}^2,$$(13)

where $c_{\mathrm{s}}=\sqrt{k_{\mathrm{B}}{T}_{\text{gas}}/\left(\mu {\mathrm{m}}_{\mathrm{H}}\right)}$ is the thermal sound speed measured in units of km s⁻¹, assuming T_gas=T_dust (see Sect. 3.2.2).

The bottom panel in Fig. 8 shows the map of the β parameter toward ρ Oph-A in band D. We obtained β < 1 across the map with a median value of 0.02. This result implies that the cloud is supported by strong magnetic pressure, which is much more prominent than that induced by the thermal motions of the gas.

3.7 Virial mass analysis

To probe the dynamical state of the ρ Oph-A cloud, we calculated the virial parameter (α_vir), which is the ratio of the virial mass (M_vir) and the isothermal mass (M_iso) of the cloud. Here, we assume the cloud as a uniform density sphere with an effective radius of R_eff as a whole, having a mean mass density of ⟨ρ⟩. The isothermal mass, M_iso, considering the self-gravity of the cloud can be estimated as: $$M_{\mathrm{i}\mathrm{s}\mathrm{o}}=\frac{4}{3}\pi ⟨\rho ⟩R_{\mathrm{e}\mathrm{f}\mathrm{f}}^3.$$(14)

The virial mass, M_vir, determined from magnetic, thermal, and non-thermal kinetic energies can be estimated as (see e.g., Bertoldi & McKee 1992; Pillai et al. 2011): $$M_{\mathrm{v}\mathrm{i}\mathrm{r}}=\frac{5R_{\mathrm{e}\mathrm{f}\mathrm{f}}}{G}(\frac{{⟨v_{\mathrm{A}}⟩}^2}{6}+{⟨c_{\mathrm{s}}⟩}^2+{⟨{\sigma }_{\mathrm{N}\mathrm{T}}⟩}^2),$$(15)

where ⟨(υ_A), ⟨c_s⟩, and ⟨σ_NT⟩ are the mean values of the Alfvénic velocity, thermal sound speed, and of non-thermal velocity dispersion toward the cloud, respectively, while G is the gravitational constant.

We calculate R_eff by selecting an ellipse covering the area where B_pos was measured. This ellipse is centered at (RA, Dec) = $({16}^{\mathrm{h}}{26}^{\mathrm{m}}{29}_{\cdot }^{\mathrm{s}}5,-{24}^{\circ }{23}^{\mathrm{\prime }}{50}_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}0)$ with the semi-major and semi-minor axes are a=120″, and b=75″, respectively. The R_eff is then determined as $R_{\text{eff}}=\sqrt{ab}\approx {95}^{\prime \prime }$ (~0.063 pc assuming the distance to ρ Oph-A of 137 pc). Within this area, the mean gas volume density is $⟨n_{{\mathrm{H}}_2}⟩=2.2\times {10}^5{\mathrm{cm}}^{-3}$, corresponding to the mean mass density ⟨ρ⟩ = 1.01 × 10⁻¹⁸g cm⁻³. The corresponding mean values of velocities are: ⟨υ_A⟩ = 3.36 km s⁻¹, ⟨c_s)⟩ = 0.09 km s^−1, and ⟨σ_NT⟩ = 0.59 km⁻¹. We used these mean values of the velocities, non-thermal velocity dispersion, and the gas volume density to calculate M_iso and M_vir in Eqs. (14) and (15), yielding α_vir = M_vir/M_iso ≈ 15. This shows that the cloud is gravitationally unbound and that B-fields and turbulence can support the cloud against self-gravity. This high value of α obtained in the cloud overalls could explain the low efficiency of star formation in ρ Oph-A; however, in a few regions, the density is very high (up to 10⁵ or 10⁶ cm⁻³). This is consistent with the lack of prestellar cores and protostellar sources in the central region, which can be seen in the lower density parts of the cloud only (see e.g., Fig. 8). To study smaller scales across the clump, polarimetric observations with a higher spatial resolution would be required.

4 Discussion

4.1 Morphology and strengths of magnetic field in ρ Oph-A

The B-field orientation in ρ Oph-A is well-ordered and predominantly perpendicular to the ridge of the dense cloud (Sec. 3.1). The feature is illustrated most clearly by the data at the longer wavelength (154 μm). The data at 89 μm better traces warmer dust in the lower density region of ρ Oph-A, which surrounds the Oph-S1 star.

The B-field properties measured using SOFIA/HAWC+ polarization data are qualitatively consistent with previous measurements of B-field orientations in ρ Oph-A obtained at 850 μm using SCUBA-2/POL-2 (Kwon et al. 2018) and SCUBA/SCUPOL (Matthews et al. 2009) on JCMT and at NIR wavelengths (JHK bands) using IRSF/SIRPOL (Kwon et al. 2015). In the POL-2 dataset, ρ Oph-A was divided into ten sub-regions, identified based on their B-field position angle and the polarization fraction, covering the entire extent of HAWC+ maps (see Fig. B.1). The highest density region is characterized by an average B-field position angle of 50°, which is also consistent with the NIR data toward the L1688 main cloud (Kwon et al. 2015). Table 1 shows the median B-field position angle (θ_B) in each subregion, which has been covered by the FIR and sub-mm observations a well. The median B-field position angles from the 89, 154, and 850 μm datasets are in rather good agreement, with the difference in θ_B values ranging from 3° to 19° (≤30%). These differences could be explained by the fact that the sub-mm data probe colder dust layers of the cloud, while the FIR data trace warmer layers. Indeed, such foreground layers were identified through the detection of the self-absorption features in the spectra of the [O I] line at 63 μm (Liseau et al. 2006) and the [CII] line at 158 μm (Mookerjea et al. 2018).

Kwon et al. (2018) calculated B_pos in three subregions associated with the densest parts of ρ Oph-A using non-thermal velocity dispersion from the N₂H⁺(1–0) line (André et al. 2007) and the sub-mm data at 850 μm In our study, we calculated the pixel-by-pixel B_pos toward the entire ρ Oph-A region. To facilitate comparisons with results obtained from the sub-mm data, we calculated the mean value of B_pos in two out of these three subregions, covered by the HAWC+ data (see Table 1). Mean values of B_pos using FIR data are lower than that obtained with the sub-mm data, which traces colder dust layers in the cloud characterized also by the higher densities. The HAWC+ data, on the other hand, essentially traces the warmer layers in the lower density region. In addition, Kwon et al. (2018) adopted values of the polarization angle dispersion as the standard deviation of the polarization angles, which are much smaller than those using the spread of the dust polarization angles (see Sect. 3.2.1). As a consequence, the B_pos obtained with sub-mm data is expected to be stronger than the one obtained with FIR data. We find that our B_pos estimate for ρ Oph-A is typically higher than those in other molecular clouds, ranging from 50 to 400 μG (for e.g., Coudé et al. 2019; Ngoc et al. 2021; Pattle et al. 2021; Ward-Thompson et al. 2023), but similar to those found in the DR21 (OH) region (Poidevin et al. 2013) or the OMC-1 cloud in Orion A region (Pattle et al. 2017; Hwang et al. 2021) where B_pos of a few mG was measured.

In this study, we adopt B_pos component as the total B-field strength. If we consider the contribution of the LOS component, the total B-field strength will be stronger than what we measured. The mass-to-flux ratio, Mach Alfvénic number, and plasma β parameter will decrease, but this will not significantly change our conclusion that B-fields are generally dominant over gravity, turbulent, and gas kinetic energies, especially toward the center of the cloud. Furthermore, as presented in Sect. 3.3, we test that the B_pos strengths obtained with the DCF method in subregions (a) and (e) are rather in good agreement with values obtained with the ADF method but higher than those using the ST and DMA methods (see Table 2). We use new B_pos values estimated with ADF, ST, and DMA methods toward two subregions (a) and (e) to estimate the mass-to-flux ratio, Mach Alfvénic number, and β plasma, which are all higher than that measured with the DCF method (see Table 2). Subregion (a) changes from magnetically sub/trans-critical to super-critical but remains to be super-Alfvénic, and has β < 1 in overall. A similar behavior is seen in subregion (e) except that this subregion changes from trans- to super-Alfvénic (see Table 2).

We go on to assume that the B-field strength in the entire cloud decreases by a factor of ~4.5, as measured with the ST method. In this case, the majority of the cloud will become magnetically trans- or super-critical (λ ≳1) except for the regions in the vicinity of the Oph-S1 star (see the top panel of Fig. D.2). However, this trend of λ remains similar as obtained using $B_{\mathrm{pos}}^{\mathrm{DCF}}$. The cloud will become mostly super-Afvénic with ℳ_A>1, except for the southeast region of the cloud (see the middle panel of Fig. D.2). This assumption will also result in the Alfvénic velocity decreasing by a factor of ~4.5 and thus, M_vir ≈ 2M_iso. This implies that the cloud is in the magneto hydrostatic equilibrium and may undergo gravitational collapse due to turbulent fluctuations on smaller scales.

4.2 Magnetic field orientation versus density structure of ρ Oph-A

To investigate what is the relative orientation of the B-fields with respect to the gas structure in ρ Oph-A, we quantify the projected offset angle between the orientations of B-field vectors projected in the plane-of-the-sky and the main structure associated with the ridge of the cloud, given as σθ = θ_B − θ_ridge, where θ_ridge, is the position angle of the ridge. We identified the shape of the ridge using the Python package RadFil⁴ (Zucker & Chen 2018, see Appendix E).

Figure 9 shows the map of offset angles (σθ) in ρ Oph-A, overlaid by contours of continuum emission at 154 μm. Clearly, the offset angle appears to increase when moving from the lower density regions to the denser regions, indicating that the relative orientation of B-fields with respect to the main structure of the cloud turns from parallel to perpendicular toward the high-density region. This is in agreement with Planck polarization observations at 353 GHz (850 μm) toward the entire L1688 region where the relative orientation of B-fields, with respect to the gas structure changes gradually. It runs parallel in the low-density region to perpendicular toward the higher density region, with the transition occurring at $N_{{\mathrm{H}}_2}\sim 3\times {10}^{21}{\mathrm{cm}}^{-2}$ and corresponding to the visual extinction A_V of ~3.3 mag⁵ (Soler 2019). This transition has also been observed in other nearby molecular clouds in the Gould Belt using Planck 353 GHz polarization observations, typically occurring at $N_{{\mathrm{H}}_2}\sim 2.5\times {10}^{21}{\mathrm{cm}}^{-2}$ or A_V ~2.7 mag (see e.g., Planck Collaboration Int. XXXV 2016; Soler et al. 2017). This phenomenon can be explained by the motion of the high-density material during the self-gravitational collapse. When the contraction occurs, the material becomes denser and compressed, which could cause the B-field lines to be re-configured. Eventually, the B-field lines are dragged and reshaped to become more uniform and perpendicular to the high-density structure of the field.

The HAWC+ data traces B-fields in the dense region of ρ Oph-A with $N_{{\mathrm{H}}_2}\sim (3-240)\times {10}^{21}{\mathrm{cm}}^{-2}$ (see Fig. 1), which corresponds to A_V in range of ~3.2–256 mag. Thus, we further investigate in detail the relation between the relative orientation of B-fields with respect to the cloud structure and the A_V within ρ Oph-A, as shown in Fig. 10. Clearly, in the regime where A_V ≲ 30 mag the offset angle (Δθ) decreases with A_V and drops from ~60° at A_V=4 mag down to 0° at A_V ~ 20 mag. This indicates that the relative orientation of B-fields changes from perpendicular back to parallel at A_V ~ 18–33 mag. Similar trends of this transition were found in the dense region of the Serpens South filament (Pillai et al. 2020) and in the Serpens Main dense cloud (Kwon et al. 2022), occurring at A_V ~ 21 mag and ~50 mag, respectively. This is likely because the B-fields are already dragged into the flow of very dense material during the gravitational contraction. Interestingly, we observe that the offset angle then progressively increases with A_V ≳ 40 mag in our region, indicating the relative orientation of B-field with respect to the ridge changes back from parallel to perpendicular (with the transition at A_V ~ 42–60 mag). We note, however, that this feature was not seen in the Serpens South filamentary cloud, whose maximum is A_V ~ 60 mag (see Fig. 3 in Pillai et al. 2020); however, it was seen in the Serpens Main cloud region, where A_V ≳ 170 mag (see Fig. 6 in Kwon et al. 2022).

Fig. 9 Map of the offset angle between B-field orientations and main cloud structure in ρ Oph-A. White contours show the continuum emission at 154 μm with similar levels as in Fig. 4. The dashed circles centered at the position of the Oph-S1 star have a radius of 60″ (black), 90″ (blue), and 120″ (green), corresponding to the levels of AV of ~6, 30, and 200 mag. The star symbol marks the position of the Oph-S1 star.

Fig. 10 Distribution of the offset angles between the orientations of B-fields and the ridge in ρ Oph-A, as a function of AV (circle symbols). The data in gray color represent the offset angle collected within the black circle “I” centered at the position of the Oph-S1 star with a radius of 60″ in Fig. 9. The data in blue and green additionally represent the offset angle collected within the two rings “II” and “III” in Fig. 9, respectively. The red diamonds represent the median value of the offset angles in each interval of AV=3 mag. The number of red diamonds is reduced to improve the clarity of the figure. The red curve connecting all the red diamonds shows the trend of the offset angles vs. AV in ρ Oph-A. The cyan box indicates the range of AV=42–60 mag where relative orientation between B-fields and the ridge changes from parallel to perpendicular. The gray horizontal dashed lines indicate perpendicular, 45°, and parallel relative orientations of B-fields with respect to the ridge.

4.3 Magnetic field versus gas density

The B-field strength (B_tot or B) can be expressed as a function of $n_{\mathrm{H}}(n_{\mathrm{H}}=2n_{{\mathrm{H}}_2}):B\propto n_{\mathrm{H}}^k$. Theoretical models predict k ~ 2/3 in the weak B-field regions (Mestel 1966) and k ≲ 0.5 in the strong B-field regions, where magnetic energy dominates over the gravity (Mouschovias & Ciolek 1999). Consequently, in the latter case, it suggests that n_H would increase faster than B if there is any collapse in the dense cloud.

Here, we show in Fig. 11 the relation log B − log n_h toward the ρ Oph-A cloud, using B-fields measured with the DCF method in band D. We find that for densities ≲10^5.2 cm⁻³, the B-fields distribution tends to be flattened. In a regime with densities higher than 10^5.2 cm⁻³ we start seeing the increasing trend of B-fields with respect to gas densities. Thus, we decided to test the log B − log n_h relation at the densest part of the cloud only. We chose this region as the area within the third contour of Stokes-I at band D (see in the bottom right panel of Fig. 5), where Stokes-I ≥ 0.4 Jy/arcsec². Within this region, it reveals a relatively weak correlation between the B-fields and densities (with the Pearson coefficient of the correlation is found r ~ 0.32, corresponding to a significance of ~5σ). We performed the power-law fit using the Python curve_fit routine (log B = k × log(n_H) + b), using mean values of B estimated within every interval of log n_H=0.08⁶. This gives a power-law index of k = 0.47 ± 0.18. In Fig. 11, we also overplot the results of B_pos using other methods in subregions (a) and (e). We find that the fitted line well covers the data point representing the resulting B_pos value using ADF method and shows good agreement with those estimated with ST and DMA in subregion (a) within and uncertainty of 1σ. The resulting k = 0.47 ± 0.18 is found consistent with theoretical predictions for strong B-fields (Mouschovias & Ciolek 1999). This is also consistent with β ≪ 1, representing strong B-fields in ρ Oph-A (see bottom panel in Fig. 8 and Sect. 3.6). We then performed the fit to $B_{\text{pos}}^{\mathrm{ST}}$ values measured with the ST method from Sect. 3.3.2 and find k = 0.49 ± 0.27. We found a similar power-law index k for the dataset at band C, namely, 0.43 ± 0.12 and 0.39 ± 0.08, using result of $B_{\text{pos}}^{\mathrm{DCF}}$ and $B_{\text{pos}}^{\mathrm{ST}}$, respectively (see Appendix C).

Our resulting k value is smaller than k ~ 0.65 found by Crutcher et al. (2010). However, their result was based on the mean values of B and n_h in different molecular clouds, which contain significant observational uncertainty and has been questioned by later studies (for e.g., Tritsis et al. 2015; Li 2021).

Fig. 11 Correlation between log B and log nH in ρ Oph-A. Small circles represent B strength estimated with the DCF method at band D. Small circles in blue further represent data within the region where Stokes-I ≥ 0.4 Jy/arcsec2, corresponding to the densities log nH ≥ 5.2. The bigger circles in green represent median values of B binned in each interval of log nH of ~0.08 (see texts). Mean values of B estimated with ADF (red triangles), DMA (blue diamonds), and ST (orange squares) methods toward subregion (a) and (e) are overplotted. The solid green line represents the median value log B ≈ 2.81 at low densities. The solid blue line shows the power-law fit toward the binned data of B using the DCF method toward more denser region (covered by the orange shaded box), resulting k = 0.47 ± 0.18. The black dash-dotted line indicates the critical power-law $B\propto n_{\mathrm{H}}^{0.5}$ for strong B-fields, predicted from the theoretical models (Mouschovias & Ciolek 1999). The black dashed line indicates the relation found by Crutcher et al. (2010) using the Zeeman observational data.

4.4 Magnetic fields versus stellar feedback

High-mass stars irradiate, heat up, and ionize the surrounding gas and create the H II region (see e.g., Zinnecker & Yorke 2007). The ionized, hot gas also expands and interacts with the vicinity environment. The radiation pressure acting on both dust and gas (assuming dust is well-coupled to the gas in the cloud) can drive the material away into the denser region. The lower density matter in the outer regions of ρ Oph-A is more exposed to the ionization front and pushed away faster than the higher density matter in the inner regions. This results in dust and gas being compressed and forming the dense shell structure at the edge of the cloud. Here, we investigate a possible impact of the radiation field from the high-mass star Oph-S1 on the ρ Oph-A cloud. The ρ Oph-A molecular cloud is directly exposed to the strong UV radiation field emitted by the Oph-S1 star from the east. Figure 12 shows the map of the continuum emission at 154 μm from SOFIA/HAWC+, superimposed by the white contours showing the H₂ column density.

We calculated the external radiation pressure of the Oph-S1 star acting on dust and gas in the densest part of the cloud, given by: $$P_{\mathrm{e}\mathrm{x}\mathrm{t},\mathrm{r}\mathrm{a}\mathrm{d}}=\frac{L_{\mathrm{b}\mathrm{o}\mathrm{l}}}{4\pi {\mathcal{R}}^{2}c},$$(16)

where L_bol=1100 L_⊙ is the bolometric luminosity of the Oph-S1 star adopted from Lada & Wilking (1984) and ℛ ≈ 0.06 pc is the projected distance from the cloud to the Oph-S1 star. We obtained P_ext,rad/k_B ≈ 2.4 × 10⁶ K cm⁻³, which is in agreement with the result calculated from the ionizing photons emitted by the Oph-S1 star (~2×10⁶ K cm⁻³, Pattle et al. 2015). Here, we omit the fact that the UV radiation from the Oph-S1 star can be absorbed by intervening dust in between the star and the ρ Oph-A region. Thus, the actual P_ext,rad acting on dust and gas in ρ Oph-A is overestimated in the Eq. (16). To make a comparison with the external radiation pressure, we calculated the internal thermal pressure due to the hot gas in ρ Oph-A, $P_{\mathrm{th}}=n_{{\mathrm{H}}_2}{k}_{\mathrm{B}}{T}_{\text{gas}}$. We considered the hot gas within the densest part of ρ Oph-A, where Stokes-I ≥0.4 Jy/arcsec² at 154 μm only, which corresponds to the area inside the third contour shown in bottom right panel of Fig. 5. Within this region, we adopted $n_{{\mathrm{H}}_2}$ in the range of 10⁴−10⁶ cm⁻³ and the average temperature of the hot gas T_gas ≈ T_dust ≈ 20 K, which yields P_th/k_B in the range of (0.2–20)×10⁶ K cm⁻³. With this assumption, we note that P_th/k_B could be underestimated, since there is a gas heating process in the cloud dominated by the UV radiation from the Oph-S1 star and by the photoelectric emission from dust exposed to the UV radiation (i.e., T_gas ≳ T_dust). The highest internal thermal pressure of ~20×10⁶ K cm⁻³ is seen toward the densest part of the cloud, associated with the position of the SM1 core. This suggests that the internal thermal pressure is strong enough to help the material refrain from external pressure. We note that, however, in the vicinity region away from the center (e.g., the tails of the ridge shown in Fig. 12), the internal thermal pressure is weaker, at least one order of magnitude less than the external one. In this case, the radiation pressure from the Oph-S1 star on the east side is dominant and able to push away all the matter into the surrounding environment toward the west side. However, it is clear that the material in the low-density region still survives, despite owning a low internal thermal pressure (see Fig. 12). This indicates that there must be another mechanism at play, for example, the presence of strong B-fields of a few mG in the cloud. We roughly estimated the magnetic pressure (P_mag = B²/(8π)) to be 10⁷−10⁸ K cm⁻³ for B-field strengths between 1–2 mG. This indicates that the magnetic pressure is at least two orders of magnitude higher than the external pressure of ~10⁶ K cm⁻³, sufficient to support the matter against the external pressure along with the gas thermal pressure.

We further investigated the interplay between B-fields and radiation directions on the evolution of the cloud. We determined the main radiation direction from the Oph-S1 star toward the densest region associated with the position of the starless core SM1 as the sky-projected direction from the Oph-S1 star and SM1 core. We found that the main radiation direction is ~70°, measured east of north. The main B-field orientation toward the densest region of ρ Oph-A of ~74° is almost parallel to the radiation field. The strong radiation pressure from the Oph-S1 star pushes all the material away and creates a dense shell toward the west side. To prevent dust and gas compression in the dense region, the B-fields exert a magnetic pressure perpendicular to the B-field orientation, seen perpendicular to the direction of the radiation pressure. This suggests that the magnetic pressure does not help the matter resist the external pressure, but tends to stretch the structure out along the northwest-southeast direction, resulting in an elongated shape of the dense shell. In the west side, the external thermal pressure is blocked by the gas thermal pressure due to the high-density material, resulting in the ridge shape in the cloud with very high densities. These results are consistent with results from 3D MHD simulation, which predicted that if the B-field orientation is parallel or along with the radiation direction, the cloud will tend to have a flattened shape Henney et al. (2009). However, in the northeast and southeast regions, the offset angles between the B-field vectors and radiation field are larger; namely, the B-field orientation becomes perpendicular to the radiation direction. In this case, the magnetic pressure and tension induced by B-fields run parallel with the radiation pressure, suggesting magnetic pressure can support the matter to confront the radiative compression and prevent the dust and gas from being swept away. This results in bending the material structure and forming a dense ridge with an umbrella-like shape. Interestingly, in the case study of ρ Oph-A the ionization source is very close to the cloud (<0.1 pc) and the offset position angle between the radiation and B-fields is not homogeneous; instead, it varies from parallel to orthogonal. This allows us to closely examine how the interplay between the B-fields and radiation fields could influence the material structure of the cloud.

Fig. 12 Image of the continuum emission at 154 μm toward Oph-A. B-field vectors obtained from SOFIA/HAWC+ at 154 μm are shown in magenta solid segments. The green dashed arrows indicate the radiation directions from the Oph-S1 star toward the cloud. The white star marks the position of the Oph-S1 star. The red curve illustrates the elongated shape of the ridge in the dense cloud. The gray contours are the H2 column density, with contour levels at (1, 2, 5, 10, 20)×l022 cm−2.

5 Summary

In this paper, we present our analysis of the Ophiuchus dark cloud complex, using the SOFIA/HAWC+ polarimetric observations toward the ρ Oph-A molecular cloud in the L1688 active star-forming region. Our main results are as follows:

• Using the dust polarization maps at 89 and 154 μm, we derived the map of B-field morphology in ρ Oph-A. We found that B-fields are generally well-ordered and primarily perpendicular to the ridge of the cloud toward the higher density regime. We estimated B_pos in the entire region of ρ Oph-A using the DCF method, ranging from 0.2 to 2.3 mG at 89 and to 2.5 mG at 154 μm. We found that B_pos is prominent toward the densest region of the cloud associated with the starless core SM1 and decreases toward the outskirts. This is the first time we have successfully estimated the distribution of the B-field strengths in ρ Oph-A, which is essential to understanding how B-fields vary in the molecular cloud globally. Our derived B-field morphology and strengths well agree with previous measurements using the sub-mm dust polarization (Kwon et al. 2018). We also estimated the B_pos strength using more advanced techniques, including the ADF, ST, and DMA, toward two subregions associated with the densest part of the cloud. While the B_pos strengths estimated with the ADF method were in agreement with those from DCF, the ST and DMA methods provided lower B_pos values of ~4.5 and 2–3 times, respectively;

• We studied the role of B-fields relative to gravity and turbulence by calculating the maps of the mass-to-flux ratio, Alfvénic Mach number, and plasma β parameter toward ρ Oph-A, using the inferred $B_{\mathrm{pos}}^{\mathrm{DCF}}$ from SOFIA/HAWC+ data at 154 μm. We found that toward the central dense parts the cloud is magnetically sub-critical and sub-Alfvénic where B-fields dominate over the gravity and turbulence. We obtained a result of β ≪ 1 in ρ Oph-A, suggesting the cloud is supported by the strong magnetic pressure and much greater than the thermal gas pressure. The lower B_pos strengths obtained with the ST and DMA methods do not significantly change our conclusions on the roles of B-fields relative to gravity and turbulence on star formation. Our virial analysis in the entire cloud showed the virial parameter α_vir ~ 15, suggesting that the cloud is gravitationally unbound and able to against the self-gravitational collapse. The cloud may be in the magneto hydrostatic equilibrium and may undergo gravitational collapse due to turbulent fluctuations on smaller scales. We should notice that our rough estimation of α_vir is a global parameter, applied for the entire cloud. In contrast, ρ Oph-A is known as one of the densest parts of the active star-forming region L1688 hosting a few pre- and proto-stellar cores and embedded protostars. Higher spatial resolution observations toward the smaller scale region in the cloud will be necessary to investigate the dynamical state of the cloud in greater detail;

• We examined the role of B-fields in supporting the cloud against the radiation pressure from the nearby high-mass star Oph-S1. We found that at the center part of the cloud, the magnetic pressure indeed does not play a role in helping the material to avoid the onset of radiation pressure. However, toward the northeast and southeast regions of the cloud, the magnetic pressure positively supports the cloud to prevent both the effects of radiative pressure and gas thermal pressure.

Acknowledgements

We thank our anonymous referee for their careful reading of our manuscript and valuable suggestions, which improved the quality of this work, particularly with regard to the realistic estimate of the gas volume density in the ρ Oph-A cloud. We thank Bhaswati Mookerjea for the valuable comments and discussion at the early stage of the project and for kindly sharing the final data cube of the HCO⁺(4–3) transition line. We thank Fabio Santos for sharing the maps of dust temperatures and gas column density. We also thank Alex Lazarian for the constructive discussions concerning the DCF and DMA methods. N.L. thanks Nguyen Tat Thang for his valuable comments as part of the 3D structure of the ρ Oph-A. N.L. and A.K. acknowledge support from the First TEAM grant of the Foundation for Polish Science No. POIR.04.04.00-00-5D21/18-00 (PI: A. Karska). A.K. acknowledges also support from the Polish National Agency for Academic Exchange grant No. BPN/BEK/2021/1/00319/DEC/1. N.B.N. was funded by the Master, Ph.D. Scholarship Programme of Vingroup Innovation Foundation (VINIF), code VINIF.2023.TS.077. M.H. acknowledges support by the National Science Centre through the OPUS grant no. 2015/19/B/ST9/02959. N.L. also acknowledges support by the grant from the Simons Foundation to IFIRSE, ICISE (916424, PI: N.H.). P.N.D. and N.T.P. were funded by the Vietnam Academy of Science and Technology under project code THTETN.03/24-25. This work is based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 2002 to the University of Stuttgart. This work made use of data from the JCMT/HARP and APEX. JCMT is operated by the East Asian Observatory on behalf of National Astronomical Observatory of Japan; the UK STFC under the auspices of grant number ST/N005856/1; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; the Operation, Maintenance and Upgrading Fund for Astronomical Telescopes and Facility Instruments, budgeted from the Ministry of Finance of China. This work used the facilities of the Canadian Astronomy Data Centre (CADC) operated by the National Research Council of Canada.

Appendix A Histogram of SOFIA/HAWC+ dataset

Figures A.1 and A.2 show histograms of Stokes-I, σ_I, SNR_I, p, σ_p, and SNR_p of the SOFIA/HAWC+ data in band C and band D. Table A.1 shows the summary of the mean and RMS of the distribution of Stokes-I, p, their uncertainties and SNRs.

Fig. A.1 Histograms of Stokes-I (top panels), σI (middle panels), and SNRI (bottom panels) of the SOFIA/HAWC+ data. Left panels represent data in band C at 89 μm and right panels represent data in band D at 154 μm. The blue vertical dashed lines in each panel represent the mean value of the distribution. The red vertical dashed lines show the position of mean value plus 1σ, 2σ, and 3σ of the distribution from left to right, respectively.

Fig. A.2 Histograms of p (top panels), σp (middle panels), and SNRp (bottoms panels) of the SOFIA/HAWC+ data. Left panels represent data in band C at 89 μm and right panels represent data in band D at 154 μm. The blue vertical dashed lines in each panel represent the mean value of the distribution. The red vertical dashed lines show the position of mean value plus 1σ, 2σ, and 3σ of the distribution from left to right, respectively.

Appendix B Magnetic field orientations and strengths at different wavelengths

In Fig. B.1, we compare the B-fields vectors obtained from HAWC+ data at 89 and 154 μm (this study), and from POL-2 data at 850 μm (Kwon et al. 2018). Vectors show the mean B-field position angles at the subregions in ρ Oph-A.

Fig. B.1 Same as the right panel in Fig. 2 but with B-field vectors from POL-2 data added (red dash-dotted vectors, Kwon et al. 2018). These vectors represent the mean B-field position angles in subregions (see texts). White and black dashed vectors show the corresponding mean B-field position angles at 89 and 154 μm.

Figure B.2 shows the histogram distribution of the estimated $B_{\mathrm{pos}}^{\mathrm{DCF}}$ obtained with HAWC+ data at 89 and 154 μm. Figure B.3 shows the maps of B-field strength estimated using the ST method $\left(B_{\mathrm{pos}}^{\mathrm{ST}}\right)$ at 89 and 154 μm.

Fig. B.2 Histogram distribution of the estimated Bpos using the DCF method, obtained with HAWC+ data at 89 μm (in blue color with hatches) and 154 (in orange color) μm.

Fig. B.3 Maps of the B-field strength projected in the plane-of-the-sky in ρ Oph-A obtained by using the ST method, $B_{\text{pos}}^{\mathrm{ST}}$, in band C (left panel) and band D (right panel). Other notations are similar as in Fig. 6.

Appendix C Velocity dispersion traced by other molecular lines

To test whether the FWHM map of HCO⁺(4–3) line is a good choice, we further use the N₂H⁺(3–2) molecular data observed by APEX to estimate the total velocity dispersion, σ_υ. The majority of the ρ Oph-A cloud exhibits the N₂H⁺(3–2) emission showing 3 single Gaussian profiles; the middle peak is much stronger compared with the other two components. This strongest hyperfine line indicates the central velocity which is in the range of 2.54–3.97 km s⁻¹, with an average velocity of 3.42 km s⁻¹, in agreement with the systemic velocity of ρ Oph-A of ~3.44 km s⁻¹ estimated by André et al. (2007). To each pixel we fit the spectrum with the sum of 3 single Gaussian profiles and then use the line width of the primary peak as the σ_υ toward the computing pixel. We construct the map of the FWHM velocity dispersion, following the procedure presented in Sect. 3.2.2 for the case of the HCO⁺(4–3) molecular line.

Figure C.l shows the map of Δυ_nt using the N₂H⁺(3–2) molecular line toward ρ Oph-A. The Δυ_nt values ranges from 0.18 to 1.73 km s⁻¹ with a median value of 0.95 km s⁻¹, in good agreement with those obtained by HCO⁺(4–3) line. We note, however, that the spatial extent of the HCO⁺(4–3) emission line is largeer than that of the N₂H⁺(3–2). Thus, we use the map of FWHM of the HCO⁺(4–3) line to calculate the $B_{\mathrm{pos}}^{\mathrm{DCF}}$ strengths in ρ Oph-A. Similarly, White et al. (2015) reported a mean FWHM velocity dispersion of 1.5 km s⁻¹ using the C¹⁸O(3–2) data from JCMT/HARP in overall entire ρ Ophiuchus cloud. Using the N₂H⁺(1–0) data from IRAM, André et al. (2007) estimated an average value of the 1D nonthermal velocity dispersion toward ρ Oph-A, σ_NT=0.2 km s⁻¹, corresponding to Δυ_NT ~0.47 km s⁻¹. Observations of the NH₃ hyperfine emission from the Green Bank Telescope also returned a similar 1D velocity dispersion with the average value of 0.3 km s⁻¹, corresponding to an average of FWHM of 0.65 km s⁻¹ (Chen et al. 2019).

Fig. C.1 Distribution of the non-thermal velocity FWHM of the N2H+(3–2) molecular line toward ρ Oph-A. White star marks the position of the Oph-S1 star. White contours show the continuum emission (Stokes-I) at 154 μm, with contours levels at 0.12, 0.2, 0.4, 0.6, and 0.8 Jy/arcsec2.

Appendix D Mass-to-flux ratio, Alfvénic Mach number, beta plasma, and B − n_H relation results obtained by using $B_{\text{pos}}^{ST}$

Figure D.1 shows the maps of mass-to-flux ratio (top panel), Alfvénic Mach number (middle panel), and beta plasma (bottom panel) toward ρ Oph-A using $B_{\mathrm{pos}}^{\mathrm{DCF}}$ in band C. Figure D.2 show maps of mass-to-flux ratio (top), Alfvénic Mach number (middle), and beta plasma (bottom) toward ρ Oph-A using $B_{\text{pos}}^{\mathrm{ST}}$ in band C (left panels) and D (right panels).

Figure D.3 shows the relation between log B and log n_H using $B_{\text{pos}}^{\mathrm{ST}}$ in band D.

Fig. D.1 Similar to Fig. 8 but using $B_{\text{pos}}^{\text{DCF}}$ in band C.

Fig. D.2 Similar to Fig. 8 but using $B_{\text{pos}}^{\mathrm{ST}}$ in band C (left panels) and D (right panels).

Fig. D.3 Similar to Fig. 11 but using B-field strength estimated with the ST method in band D. Small squares represent B strength estimated with the ST method at band D. Blue small squares further show the data at the denser region (log nH ≥ 5.2) with Stokes-I ≥ 0.4 Jy/arcsec2. Green bigger squares represent median value of log B binned in each interval of log nH of ~0.08 in the denser region. Mean values of log B estimated with ADF (red triangles) and DMA (blue diamonds) methods toward subregions (a) and (e) are overplotted. The solid green line shows the median data at lower densities region where log nH ≥5.2. The solid blue line shows the power-law fit toward the binned data of B using the DCF method toward denser region (covered by the orange shaded box, resulting k = 0.49 ± 0.27. The black dash-dotted line indicates the critical power-law $B\propto n_{\mathrm{H}}^{0.5}$ for strong B-fields, predicted from the theoretical models (Mouschovias & Ciolek 1999). The black dashed line indicates the relation found by Crutcher et al. (2010) using the Zeeman observational data.

Appendix E Main structure of the ridge identification

We identified the elongated shape of the ridge in the densest region of ρ Oph-A using the Python package RadFil. We first applied a threshold of Stokes-I ≥0.5 Jy/arcsec² in the map of Stokes-I at 154 μm to build a mask for the main shape of the densest region (the ridge). The spine of the ridge was identified by using the Python package FilFinder⁷ from Koch & Rosolowsky (2015). We then smoothed the ridge spine and define a tangent line at every position on the ridge, where the orientation of the ridge (θ_ridge) at every position was identified. Figure E.1 shows the schematic view of the ridge and the perpendicular cuts at each position on the ridge. The position angle of the tangent line at each position located on the ridge is defined as the orientation of the ridge at that position, measured east of north, and ranges from 0° to 180°.

Fig. E.1 Schematic view of the ridge in ρ Oph-A. The background color map shows the Stokes-I intensity at 154 μm. The thick red line indicates the main shape of the ridge. Thin red lines show the lines perpendicular to the tangent line at each position on the ridge.

References

All Tables

All Figures

Fig. 1 Distribution of dust temperatures (colors) and H2 column densities (contours) in ρ Oph-A (Santos et al. 2019). Contour levels of $N_{{\mathrm{H}}_2}$ are shown at (1, 2, 4, 6, 8, and 10) × 1022 cm−2. The star symbol shows the position of Oph-S1. Small circles show the positions of YSOs (Enoch et al. 2009; Evans et al. 2009; Connelley & Greene 2010). Yellow squares and cyan diamonds indicate the positions of prestellar and protostellar cores, respectively, identified as part of the Herschel Gould Belt Survey (Ladjelate et al. 2020). The “×” symbol shows the position of the SM1 core.

In the text

Fig. 2 Maps of the B-field orientation (white segments) on top of the total intensity (Stokes-I, color maps) toward ρ Oph-A, observed with SOFIA/HAWC+ in band C (left) and band D (right), centered at 89 and 154 μm, respectively. The gray contours in each map show the Stokes-I intensity at the corresponding wavelength, with contour levels at 0.12, 0.2, 0.4, 0.6, and 0.8 Jy/arcsec2. The star symbol marks the position of Oph-S1 star. Small white circle symbols indicate the positions of YSOs. The filled orange circles show the beam size of each color map.

In the text

	Fig. 3 Polarization angle dispersion maps using data in band C (left panel) and band D (right panel). The black contours are from the total intensity map at the corresponding wavelength, with the contour levels at 0.12, 0.2, 0.4, 0.6, and 0.8 Jy/arcsec2. The star symbol in each map indicates the position of the Oph-S1 star.
In the text

	Fig. 4 Distribution of the non-thermal FWHM of the HCO+(4–3) molecular line toward ρ Oph-A. The star symbol marks the position of the Oph-S1 star. White contours show the continuum emission (Stokes-I) at 154 μm, with contours levels at 0.12, 0.2, 0.4, 0.6, and 0.8 Jy/arcsec2.
In the text

Fig. 5 Cartoon of 3D structure of the ρ Oph-A cloud where we assume that the cloud has a spherical shape but is pushed by radiation feedback from a nearby high-mass star Oph-S1, shown at the top (left panel, see text) and a plane of the constant coordinate y with the LOS $\hat{z}$ (right panel). Bottom row displays the map of the depth of the cloud, W (left panel), and the map of the gas volume density $\left(n_{{\mathrm{H}}_2}\right)$ inferred from $N_{{\mathrm{H}}_2}$ (right panel). In each panel, the “×” and star symbols mark the peak position of $N_{{\mathrm{H}}_2}$ and location of Oph-S1 star, respectively. Toward maps of W and $n_{{\mathrm{H}}_2}$ in bottom panels, the black dashed circle has a radius Rcl=120″ and black solid contours show the continuum emission at 154 μm with levels similar as shown in Fig. 4.

In the text

	Fig. 6 Maps of B-field strengths projected in the plane-of-the-sky (Bpos) in ρ Oph-A using data in band C (left panel) and band D (right panel), respectively. The contours in each map are similar as in Fig. 2. The star symbol marks the position of the Oph-S1 star and small circle symbols show the position of YSOs.
In the text

Fig. 7 Structure function of polarization angles (top panels) and of central velocity (middle panels) measured in the subregions (a) (left) and (e) (right). The bottom panels show the resulting $B_{\text{pos}}^{\text{DMA}}$ strength estimated via Eq. (10) in the two subregions (a) and (e). The HAWC+ data at 154 μm is used to calculate the structure-function (data points). The red solid line in the top panels shows the fitted line of the angular structure-function to the quadratic function, where we chose data points with 13.6″ < ℓ ≤ 50″. The blue dash-dotted line in each panel indicates the upper ℓ-limit of 50″ used for fitting the angular structure-function and for estimating the mean value of $B_{\mathrm{pos}}^{\mathrm{DMA}}$ in the two subregions (a) and (e) (see texts). The blue dashed line in each panel indicates the beam size of HAWC+ observation at 154 μm.

In the text

Fig. 8 Maps of λ (top panel), Alfvén Mach number (middle panel), and plasma β parameter (bottom panel) in ρ Oph-A using the $B_{\mathrm{pos}}^{\mathrm{DCF}}$ in band D. In each map, yellow square and cyan diamond symbols indicate the positions of prestellar and protostellar cores, respectively. Green circles indicate the positions of YSOs (Enoch et al. 2009; Evans et al. 2009; Connelley & Greene 2010). The star and “×” symbols mark the positions of the Oph-S1 star and the SM1 core, respectively. Black contours show $N_{{\mathrm{H}}_2}$, similar as in Fig. 1.

In the text

Fig. 9 Map of the offset angle between B-field orientations and main cloud structure in ρ Oph-A. White contours show the continuum emission at 154 μm with similar levels as in Fig. 4. The dashed circles centered at the position of the Oph-S1 star have a radius of 60″ (black), 90″ (blue), and 120″ (green), corresponding to the levels of AV of ~6, 30, and 200 mag. The star symbol marks the position of the Oph-S1 star.

In the text

Fig. 10 Distribution of the offset angles between the orientations of B-fields and the ridge in ρ Oph-A, as a function of AV (circle symbols). The data in gray color represent the offset angle collected within the black circle “I” centered at the position of the Oph-S1 star with a radius of 60″ in Fig. 9. The data in blue and green additionally represent the offset angle collected within the two rings “II” and “III” in Fig. 9, respectively. The red diamonds represent the median value of the offset angles in each interval of AV=3 mag. The number of red diamonds is reduced to improve the clarity of the figure. The red curve connecting all the red diamonds shows the trend of the offset angles vs. AV in ρ Oph-A. The cyan box indicates the range of AV=42–60 mag where relative orientation between B-fields and the ridge changes from parallel to perpendicular. The gray horizontal dashed lines indicate perpendicular, 45°, and parallel relative orientations of B-fields with respect to the ridge.

In the text

Fig. 11 Correlation between log B and log nH in ρ Oph-A. Small circles represent B strength estimated with the DCF method at band D. Small circles in blue further represent data within the region where Stokes-I ≥ 0.4 Jy/arcsec2, corresponding to the densities log nH ≥ 5.2. The bigger circles in green represent median values of B binned in each interval of log nH of ~0.08 (see texts). Mean values of B estimated with ADF (red triangles), DMA (blue diamonds), and ST (orange squares) methods toward subregion (a) and (e) are overplotted. The solid green line represents the median value log B ≈ 2.81 at low densities. The solid blue line shows the power-law fit toward the binned data of B using the DCF method toward more denser region (covered by the orange shaded box), resulting k = 0.47 ± 0.18. The black dash-dotted line indicates the critical power-law $B\propto n_{\mathrm{H}}^{0.5}$ for strong B-fields, predicted from the theoretical models (Mouschovias & Ciolek 1999). The black dashed line indicates the relation found by Crutcher et al. (2010) using the Zeeman observational data.

In the text

Fig. 12 Image of the continuum emission at 154 μm toward Oph-A. B-field vectors obtained from SOFIA/HAWC+ at 154 μm are shown in magenta solid segments. The green dashed arrows indicate the radiation directions from the Oph-S1 star toward the cloud. The white star marks the position of the Oph-S1 star. The red curve illustrates the elongated shape of the ridge in the dense cloud. The gray contours are the H2 column density, with contour levels at (1, 2, 5, 10, 20)×l022 cm−2.

In the text

Fig. A.1 Histograms of Stokes-I (top panels), σI (middle panels), and SNRI (bottom panels) of the SOFIA/HAWC+ data. Left panels represent data in band C at 89 μm and right panels represent data in band D at 154 μm. The blue vertical dashed lines in each panel represent the mean value of the distribution. The red vertical dashed lines show the position of mean value plus 1σ, 2σ, and 3σ of the distribution from left to right, respectively.

In the text

Fig. A.2 Histograms of p (top panels), σp (middle panels), and SNRp (bottoms panels) of the SOFIA/HAWC+ data. Left panels represent data in band C at 89 μm and right panels represent data in band D at 154 μm. The blue vertical dashed lines in each panel represent the mean value of the distribution. The red vertical dashed lines show the position of mean value plus 1σ, 2σ, and 3σ of the distribution from left to right, respectively.

In the text

	Fig. B.1 Same as the right panel in Fig. 2 but with B-field vectors from POL-2 data added (red dash-dotted vectors, Kwon et al. 2018). These vectors represent the mean B-field position angles in subregions (see texts). White and black dashed vectors show the corresponding mean B-field position angles at 89 and 154 μm.
In the text

	Fig. B.2 Histogram distribution of the estimated Bpos using the DCF method, obtained with HAWC+ data at 89 μm (in blue color with hatches) and 154 (in orange color) μm.
In the text

	Fig. B.3 Maps of the B-field strength projected in the plane-of-the-sky in ρ Oph-A obtained by using the ST method, $B_{\text{pos}}^{\mathrm{ST}}$, in band C (left panel) and band D (right panel). Other notations are similar as in Fig. 6.
In the text

	Fig. C.1 Distribution of the non-thermal velocity FWHM of the N2H+(3–2) molecular line toward ρ Oph-A. White star marks the position of the Oph-S1 star. White contours show the continuum emission (Stokes-I) at 154 μm, with contours levels at 0.12, 0.2, 0.4, 0.6, and 0.8 Jy/arcsec2.
In the text

	Fig. D.1 Similar to Fig. 8 but using $B_{\text{pos}}^{\text{DCF}}$ in band C.
In the text

	Fig. D.2 Similar to Fig. 8 but using $B_{\text{pos}}^{\mathrm{ST}}$ in band C (left panels) and D (right panels).
In the text

Fig. D.3 Similar to Fig. 11 but using B-field strength estimated with the ST method in band D. Small squares represent B strength estimated with the ST method at band D. Blue small squares further show the data at the denser region (log nH ≥ 5.2) with Stokes-I ≥ 0.4 Jy/arcsec2. Green bigger squares represent median value of log B binned in each interval of log nH of ~0.08 in the denser region. Mean values of log B estimated with ADF (red triangles) and DMA (blue diamonds) methods toward subregions (a) and (e) are overplotted. The solid green line shows the median data at lower densities region where log nH ≥5.2. The solid blue line shows the power-law fit toward the binned data of B using the DCF method toward denser region (covered by the orange shaded box, resulting k = 0.49 ± 0.27. The black dash-dotted line indicates the critical power-law $B\propto n_{\mathrm{H}}^{0.5}$ for strong B-fields, predicted from the theoretical models (Mouschovias & Ciolek 1999). The black dashed line indicates the relation found by Crutcher et al. (2010) using the Zeeman observational data.

In the text

	Fig. E.1 Schematic view of the ridge in ρ Oph-A. The background color map shows the Stokes-I intensity at 154 μm. The thick red line indicates the main shape of the ridge. Thin red lines show the lines perpendicular to the tangent line at each position on the ridge.
In the text