Probing Three-Dimensional Magnetic Fields: II -- An Interpretable
  Convolutional Neural Network
Observing 3D magnetic fields, including orientation and strength, within the
interstellar medium is vital but notoriously difficult. However, recent
advances in our understanding of anisotropic magnetohydrodynamic (MHD)
turbulence demonstrate that MHD turbulence and 3D magnetic fields leave their
imprints on the intensity features of spectroscopic observations. Leveraging
these theoretical frameworks, we propose a novel Convolutional Neural Network
(CNN) model to extract this embedded information, enabling the probe of 3D
magnetic fields. This model examines not only the plane-of-the-sky magnetic
field orientation ($\phi$), but also the magnetic field's inclination angle
($\gamma$) relative to the line-of-sight, and the total magnetization level
(M$_A^{-1}$) of the cloud. We train the model using synthetic emission lines of
$^{13}$CO (J = 1 - 0) and C$^{18}$O (J = 1 - 0), generated from 3D MHD
simulations that span conditions from sub-Alfv\'enic to super-Alfv\'enic
molecular clouds. Our tests confirm that the CNN model effectively reconstructs
the 3D magnetic field topology and magnetization. The median uncertainties are
under $5^\circ$ for both $\phi$ and $\gamma$, and less than 0.2 for M$_A$ in
sub-Alfv\'enic conditions (M$_A\approx0.5$). In super-Alfv\'enic scenarios
(M$_A\approx2.0$), they are under $15^\circ$ for $\phi$ and $\gamma$, and 1.5
for M$_A$. We applied this trained CNN model to the L1478 molecular cloud.
Results show a strong agreement between the CNN-predicted magnetic field
orientation and that derived from Planck 353 GHz polarization data. The CNN
approach enabled us to construct the 3D magnetic field map for L1478, revealing
a global inclination angle of $\approx76^\circ$ and a global M$_A$ of
$\approx1.07$.
Yue Hu, A. Lazarian, Yan Wu, Chengcheng Fu
# Probing Three-Dimensional Magnetic Fields: II - An Interpretable Convolutional Neural Network

###### Abstract

Observing 3D magnetic fields, including orientation and strength, within the interstellar medium is vital but notoriously difficult. However, recent advances in our understanding of anisotropic magnetohydrodynamic (MHD) turbulence demonstrate that MHD turbulence and 3D magnetic fields leave their imprints on the intensity features of spectroscopic observations. Leveraging these theoretical frameworks, we propose a novel Convolutional Neural Network (CNN) model to extract this embedded information, enabling the probe of 3D magnetic fields. This model examines the plane-of-the-sky magnetic field orientation (\(\phi\)), the magnetic field's inclination angle (\(\gamma\)) relative to the line-of-sight, and the total magnetization level (M\({}_{A}^{-1}\)) of the cloud. We train the model using synthetic emission lines of \({}^{13}\)CO (J = 1 - 0) and C\({}^{18}\)O (J = 1 - 0), generated from 3D MHD simulations that span conditions from sub-Alfvenic to super-Alfvenic molecular clouds. Our tests confirm that the CNN model effectively reconstructs the 3D magnetic field topology and magnetization. The median uncertainties are under 5\({}^{\circ}\) for both \(\phi\) and \(\gamma\), and less than 0.2 for M\({}_{A}\) in sub-Alfvenic conditions (M\({}_{A}\approx\) 0.5). In super-Alfvenic scenarios (M\({}_{A}\approx\) 2.0), they are under 15\({}^{\circ}\) for \(\phi\) and \(\gamma\), and 1.5 for M\({}_{A}\). We applied this trained CNN model to the L1478 molecular cloud. Results show a strong agreement between the CNN-predicted magnetic field orientation and that derived from Planck 353 GHz polarization. The CNN approach enabled us to construct the 3D magnetic field map for L1478, revealing a global inclination angle of \(\approx\) 76\({}^{\circ}\) and a global M\({}_{A}\) of \(\approx\) 1.07.

keywords: ISM:general--ISM:structure--ISM:magnetic field--turbulence

## 1 Introduction

In the vast interstellar medium (ISM), magnetic fields are pervasive powers that significantly influence various astrophysical phenomena. These fields serve as invisible balancers against gravitational forces within the ISM, intricately maintaining its equilibrium (Wurster & Li, 2018; Abbate et al., 2020). They are instrumental in directing gas flows towards galactic nuclei, playing a crucial role in their sentence and the dynamic processes unfolding therein (Kim & Stone, 2012; Roche et al., 2018; Busquet, 2020; Whittingham et al., 2021; Hu et al., 2022). Magnetic fields also govern the trajectories of cosmic rays, affecting the energy distribution and overall dynamics of ISM (Fermi, 1949; Jokipii, 1966; Yan & Lazarian, 2002, 2004; Ferrand & Marcowith, 2010; Xu & Yan, 2013; Xu & Lazarian, 2020; Hopkins et al., 2021; Hu et al., 2022). Furthermore, they are deeply involved in the star formation processes within molecular clouds, influencing both the rate and nature of star births (Mestel, 1965; Mac Low & Klessen, 2004; McKee & Ostriker, 2007; Lazarian et al., 2012; Federrath & Klessen, 2012; Hu et al., 2021). Despite their pivotal roles, our understanding of these magnetic fields remains far from complete.

Our primary challenge lies in the formidable task of probing a three-dimensional (3D) magnetic field in 3D spatial space. Current approaches -- such as polarized dust emission (Lazarian, 2007; Andersson et al., 2015; Planck Collaboration et al., 2015, 2020; Fissel et al., 2016; Li et al., 2021; Liu et al., 2023), polarized synchrotron emission (Xiao et al., 2008; Planck Collaboration et al., 2016; Guan et al., 2021), provide 2D measurements of the plane-of-sky (POS) magnetic field direction, while Zeeman splitting (Crutcher, 2004, 2012), and Faraday rotation (Haverkorn, 2007; Taylor et al., 2009; Oppermann et al., 2012; Xu & Zhang, 2016; Tahani et al., 2019) provide line-of-sight (LOS) components of the magnetic field. Yielding valuable insights, these techniques probe into distinct and typically different regions of the multiphase ISM. Thus, despite their individual strengths, merging these insights into a coherent, full 3D magnetic field vector, which includes both the 3D orientation and total strength, presents a non-trivial task.

A significant advance in probing the 3D magnetic fields in molecular clouds has been made by leveraging polarized dust emission, drawing on the depolarization effect induced by different magnetic field orientations (see Chen et al., 2019) and by accounting for the properties of turbulent magnetic fields Hu & Lazarian 2023a,c). As a separate development, Tahani et al. (2019, 2022) has succeeded in employing the synergy of Faraday rotation and dust polarization to infer a helical 3D magnetic field topology across the Orion A, Orion B, Perseus, and California clouds. Subsequently, Hu et al. (2021)
and Hu et al. (2021c) proposed the use of anisotropic properties of magnetohydrodynamic (MHD) turbulence, inherited by young stellar objects (Ha et al., 2022) and spectroscopic lines (Lazarian & Pogosyan, 2000; Kandel et al., 2016; Hu et al., 2023), to obtain the LOS and POS components of the magnetic field's orientation and total magnetization simultaneously.

Importantly, the underlying theory of Hu et al. (2021c)'s approach demonstrates that spectroscopic observations embody the anisotropy of MHD turbulence (Lazarian & Pogosyan, 2000; Kandel et al., 2016; Hu et al., 2023), i.e., turbulent eddies elongate along the 3D direction of the magnetic field (Goldreich & Sridhar, 1995; Lazarian & Vishniac, 1999). The spatial features present in these observations imprint the anisotropy and thus carry detailed information about the magnetic fields. This implies that, given an extensive amount of training data, machine learning algorithms have the potential to capture these features and produce accurate measurements. This strategy has been employed to map the 2D POS magnetic field orientation using velocity channel maps from spectroscopic observations (Xu et al., 2023). The theoretical basis remains the anisotropy of the MHD turbulence, a principle previously utilized to trace magnetic fields via velocity gradients (Lazarian & Yuen, 2018; Hu et al., 2018, 2022c; Alina et al., 2022; Liu et al., 2022a; Schmaltz et al., 2023). However, Hu et al. (2021c) made the crucial discovery that anisotropy in velocity channel maps harbors not only information about the POS magnetic field orientation, but also the total magnetization and the magnetic field's inclination angle with respect to the LOS. This additional information paves the way for constructing the full 3D magnetic field vector from spectroscopic observations.

By leveraging the capabilities of Convolutional Neural Networks (CNN; LeCun et al., 1998)--a type of deep learning model excelling in image and signal processing--we aspire to develop a novel method that can probe the 3D magnetic field. Earlier studies of the CNN explored the possibility to distinguish sub-Alfvenic and super-Alfvenic turbulence (Peek & Burkhart, 2019) and predict the POS magnetic field orientation (Xu et al., 2023). Our study, however, targets the simultaneous extraction of LOS and POS magnetic field orientations and the total magnetization. The foundation of our CNN model is the anisotropic MHD turbulence exhibited in spectroscopic observations (Lazarian & Pogosyan, 2000; Kandel et al., 2016; Hu et al., 2023), a theoretical underpinning that allows us to interpret the CNN model accurately. In other words, it enables us to discern the specific features that convey information about the magnetic field, the reasons why they are informative, and their underlying physical meanings. The effectiveness of training a CNN is highly dependent on the availability of comprehensive numerical simulations that accurately represent realistic ISM. In this research, we employ 3D MHD supersonic simulations that portray a range of ISM environments, spanning sub-Alfvenic scenarios (i.e., strong magnetic field), trans-Alfvenic, and super-Alfvenic conditions (i.e., weak magnetic field). We further post-process these simulations by incorporating the radiative transfer effect, which enables us to generate mock emission lines of \({}^{13}\)CO and C\({}^{18}\)O from diffuse molecular clouds. Through this trained CNN model, we present the 3D magnetic field map of the molecular cloud L1478.

This paper is organized as follows. In SS 2, we briefly review the basic concepts of MHD turbulence anisotropy in spectroscopic observations and their correlation with 3D magnetic field orientation and total magnetization. In SS 3, we give details of the 3D MHD simulations and mock observations used in this work, as well as our CNN model. We use mock observations to train the CNN model and present the results of numerical testing in SS 4. We further apply the trained CNN model to predict the 3D magnetic field in the molecular cloud L1478. In SS 5, we discuss the uncertainty and prospects of the machine learning approach, as well as implications for various astrophysical problems. We summarize our results in SS 6.

## 2 Theoretical consideration

### Anisotropy of MHD turbulence: revealing magnetic field orientation and magnetization

The earliest model of MHD turbulence was proposed to be isotropic (Iroshnikov, 1963; Kraichnan, 1965). However, this model underwent subsequent revisions through a series of theoretical and numerical studies, revealing that MHD turbulence exhibits anisotropy under sub-Alfvenic conditions and isotropy at large-scale, super-Alfvenic conditions (Montgomery & Turner, 1981; Shebalin et al., 1983; Higdon, 1984; Montgomery & Matthaeus, 1995).

A significant advance in this field was the introduction of the "critical balance" condition, i.e., equating the cascading time \((k_{\perp}v_{l})^{-1}\) and the wave periods \((k_{\parallel}v_{A})^{-1}\), proposed by Goldreich & Sridhar (1995), hereafter GS95. Here \(k_{\parallel}\) and \(k_{\perp}\) represent the components of the wavevector parallel and perpendicular to the magnetic field, respectively, while \(v_{l}\) denotes the turbulent velocity at scale \(l\), and \(v_{A}=B/\sqrt{4\pi\rho}\) represents the Alfven speed. Here \(B\) is the magnetic field strength and \(\rho\) is the gas mass density.

Taking into account Kolmogorov-type turbulence, i.e. \(v_{l}\propto l^{1/3}\), the GS95 anisotropy scaling can be straightforwardly derived.

\[k_{\parallel}\propto k_{\perp}^{2/3}, \tag{1}\]

which reveals the anisotropic nature of turbulence eddies, implying that the eddies are elongated along the magnetic fields. However, it should be noted that the considerations of GS95 are based on a global reference frame, where the direction of the wavevectors is defined relative to the mean magnetic field.

Scale-dependent anisotropy was later introduced via the study of fast turbulent reconnection by Lazarian & Vishniac 1999 (hereafter LV99), which proposed a local reference frame. This frame is defined relative to the magnetic field passing through an eddy at scale \(l\). According to LV99, the motion of eddies perpendicular to the direction of the local magnetic field adheres to the Kolmogorov law (i.e. \(v_{L,\perp}\propto l_{\perp}^{1/3}\)), since this is the direction in which the magnetic field offers minimal resistance. Applying the "critical balance" condition in the local reference frame: \(v_{L,\perp}l_{\perp}^{-1}\approx v_{A}l_{\parallel}^{-1}\), the scale-dependent anisotropy scaling is then given by:

\[l_{\parallel}=L_{\rm inj}(\frac{l_{\perp}}{L_{\rm inj}})^{\frac{2}{3}}{\rm M }_{A}^{-4/3},\ \ \ {\rm M}_{A}\leq 1, \tag{2}\]

where \(l_{\perp}\) and \(l_{\parallel}\) represent the perpendicular and parallel scales of eddies with respect to the local magnetic field, respectively. \(L_{\rm inj}\) denotes the turbulence injection scale and \({\rm M}_{A}=v_{\rm inj}/v_{A}\) is the Alfven Mach number. \({\rm M}_{A}^{-4}\) gives magnetization level of the medium.

Eq. 2 provides two critical insights: (1) **Turbulent eddies stretch along the local magnetic field (i.e., \(l_{\parallel}\gg l_{\perp}\))**, and (2) **the degree of anisotropy, defined as \(l_{\parallel}/l_{\perp}\), depends on the magnetization \({\rm M}_{A}^{-1}\). As we illustrated in Fig. 1, this indicates that eddies become increasingly anisotropic in a strongly magnetized medium. For the case where \({\rm M}_{A}\gg 1\), turbulence is essentially isotropic due to the predominance of hydrodynamic turbulence. However, the essence of turbulence lies in the cascading of energy from larger injection scales to smaller ones, which leads to a decrease in turbulent velocity. Eventually, at the transition scale \(l_{a}=L_{\rm inj}/M_{A}^{3}\), the strength of the magnetic field becomes comparable to the turbulence (i.e., the Alfven
Mach number at \(l_{a}\) is unity, see Lazarian 2006), and anisotropy starts to manifest.

Furthermore, **(3) changes in \(\mathrm{M}_{A}\) are distinctly reflected in the magnetic field topology.** Within a strongly magnetized medium, the magnetic field lines exhibit minimal variation due to the presence of weaker fluctuations, resulting in more straightened field lines. In contrast, in the context of a weaker magnetic field, which corresponds to a larger value of \(\mathrm{M}_{A}\), fluctuations in the magnetic field direction intensify significantly. This leads to the field lines adopting a more curved configuration (Yuen & Lazarian, 2020). As turbulent eddies extend along the local magnetic field, the topological changes induced by \(\mathrm{M}_{A}\) become evidently imprinted within these eddies.

### Obtaining velocity information from spectroscopic observation

The anisotropy outlined in Eq. 2 pertains to turbulent velocity fluctuations, and the turbulent eddy refers to velocity fluctuation contour. This suggests that anisotropy manifests in turbulent velocity fields. Such anisotropic velocity can be obtained from the velocity channel map of spectroscopic observations, due to the velocity caustics effect (Lazarian & Pogosyan, 2000). We briefly review this concept.

In position-position-velocity (PPV) space, the observed intensity distribution of a given spectral line is determined by both the density of emitters and their velocity distribution along the LOS. If coherent velocity shear -- for instance, from galactic rotation -- can be disregarded 1, the LOS velocity component, \(v\), becomes the sum of the turbulent velocity, \(v_{\mathrm{t}\mathrm{u}t}(x,y,z)\), and the residual component attributable to thermal motions. This residual thermal velocity, \(v-v_{\mathrm{t}\mathrm{u}t}(x,y,z)\), has a Maxwellian distribution, \(\phi(v,x,y,z)\). For emissivity proportional to density, it provides PPV emission density \(\rho_{s}(x,y,z)\) as (Lazarian & Pogosyan, 2004):

Footnote 1: The impact of galactic rotation on velocity caustics was explored by Lazarian & Pogosyan (2000). It demonstrated that its effects are insignificant (Hu et al., 2023).

\[\rho_{s}(x,y,v) =\kappa\int\rho(x,y,z)\phi(v,x,y,z)dz, \tag{3}\] \[\phi(v,x,y,z) \equiv\frac{1}{\sqrt{2\pi c_{s}^{2}}}\exp[-\frac{[v-v_{\mathrm{t }\mathrm{u}t}(x,y,z)]^{2}}{2c_{s}^{2}}], \tag{4}\]

where \(\kappa\) is a constant that correlates the number of emitters to the observed intensities. \(c_{s}=\sqrt{\gamma_{a}k_{\mathrm{B}}T/m}\) is the sound speed, with \(m\) being the mass of atoms or molecules, \(\gamma_{a}\) the adiabatic index, \(k_{\mathrm{B}}\) being the Boltzmann constant, and \(T\) the temperature, which can vary from point to point if the emitter is not isothermal. However, the variation of temperature has only a marginal contribution to the distribution of \(\rho_{s}(x,y,v)\)(see Hu et al., 2023). By integrating \(\rho_{s}(x,y,v)\) over a defined velocity range or channel width \(\Delta v\), we obtain a velocity channel:

\[p(x,y,v)=\int_{v-\Delta v/2}^{v+\Delta v/2}\rho_{s}(x,y,v^{\prime})dv^{\prime }. \tag{5}\]

By separating the 3D density into the mean density and zero-mean fluctuations, \(\rho(x,y,z)=\bar{\rho}+\bar{\rho}\phi(x,y,z)\), the channel intensity can be represented as the sum of two terms, \(p(x,y,v)=p_{vc}(x,y,v)+\bar{\rho}\phi(x,y,z)\).

Figure 1: Illustration of how the observed intensity structures in channel map regulated by \(\mathrm{M}_{A}\) and \(\gamma\). Within all three panels, these intensity structures elongate along the POS magnetic field direction where \(l_{\parallel}>l_{\perp}\). Structures 1 and 2, depicted in panels (a) and (b), are projected onto the POS with identical inclination angles \(\gamma_{1}=\gamma_{2}\), yet exhibit different magnetizations with \(\mathrm{M}_{A,1}^{-1}>\mathrm{M}_{A,2}^{-1}\). Notably, the anisotropy observed, represented as \(l_{\parallel}/l_{\perp}\), in the weakly magnetized Structure 2 is less pronounced than in Structure 1. Structure 2 is less straightened because the weak magnetic field has more fluctuations. The curvature of the observed magnetic structures is suggested for magnetization studies by Yuen & Lazarian (2020). Comparatively, Structures 1 and 3—showcased in panels (a) and (c)—possess equivalent magnetizations \(\mathrm{M}_{A,1}^{-1}=\mathrm{M}_{A,3}^{-1}\), but divergent inclination angles with \(\gamma_{1}>\gamma_{3}\). The observed anisotropy decreases with smaller \(\gamma\), though it is crucial to note that the straightness of Structure 3 remains unaffected by this projection. It should be noted that here the projection effect is simplified. The intensity structures are predominantly created by the velocity caustics effect, due to MHD turbulence. The projection effect is applied to the velocity field and then subsequent intensity structures in velocity channels.
\[p_{dc}(x,y,v)\text{(Hu et al., 2023)}\] \[p_{vc} \equiv\int_{\gamma-\Delta v/2}^{v+\Delta v/2}\!dv^{\prime}\int\bar{ \rho}\phi(v^{\prime},x,y,z)dz, \tag{6}\] \[p_{dc} \equiv\int_{v-\Delta v/2}^{v+\Delta v/2}\!dv^{\prime}\int\bar{ \rho}\delta(x,y,z)\phi(v^{\prime},x,y,z)dz. \tag{7}\]

The first term, \(p_{vc}\), encompasses the mean intensity in the channel and carries fluctuations exclusively produced by velocity, called the velocity caustics effect (Lazarian & Pogosyan, 2000). The second term, \(p_{dc}\), reflects the inhomogeneities in the real 3D density.

The relative importance of \(p_{vc}\) and \(p_{dc}\) depends on the channel width (Lazarian & Pogosyan, 2000; Kandel et al., 2016; Hu et al., 2023). The narrower the channel width, the greater the contribution from \(p_{vc}\). When the channel width \(\Delta v\) is less than the velocity dispersion \(\sqrt{\delta(v^{2})}\) of the turbulent eddies under investigation, that is, \(\Delta v<\sqrt{\delta(v^{2})}\), the intensity fluctuation in such a thin channel is predominantly due to velocity fluctuation. Consequently, \(p(x,y,v)\) inherits the anisotropy of MHD turbulence. The intensity structures within \(p(x,y,v)\) elongate along the POS magnetic fields, and their corresponding anisotropy degree, as well as the topology, is correlated with the magnetization and inclination angle. On the other hand, the dominance of \(p_{vc}\) ensures that the morphology of intensity fluctuation within \(p(x,y,v)\) is less sensitive to M\({}_{s}\), because the anisotropy in MHD turbulence's velocity field is not affected by M\({}_{s}\)(Kowal & Lazarian, 2010).

It is important to note that Clark et al. (2019) questioned the validity of velocity caustics in the presence of thermal broadening in multiphase HI gas and suggested that the thin velocity channel is dominated by density fluctuations from cold filaments. The nature of the striations in channel maps was tested in Hu et al. (2023), by explicitly evaluating velocity and density contributions in velocity channels obtained from multi-phase HI simulations and GALFA-HI observations. This study confirmed that the velocity caustics were responsible for the observed striation.

### Anisotropy in thin velocity channels: dependence on the inclination angle of magnetic fields

The anisotropy of the observed intensity in a PPV channel, represented by \(p(x,y,v)\), is also affected by the inclination angle \(\gamma\) of the magnetic field with respect to the LOS, due to the projection effect (Hu et al., 2021). For example, as illustrated in Fig. 1, we consider two magnetized structures (or eddies), \(s_{1}\), and \(s_{3}\), both having identical magnetization. Although these unprojected structures have the same anisotropy degree, their projections differ. Specifically, a projection with a smaller inclination angle results in a lower anisotropy degree by reducing the scale parallel to the magnetic fields. When \(\gamma=0\), the parallel scale of the eddy aligns with the LOS, making the anisotropy unobservable on the POS.

However, as previously mentioned, the degree of anisotropy is also controlled by magnetization. As shown in Fig. 1, although two magnetized structures (\(s_{1}\) and \(s_{2}\)) share identical inclination angles, the projection of the weakly magnetized \(s_{2}\) shows less anisotropy. Importantly, the topology of \(s_{2}\) is further changed being less straightened. This is because a weak magnetic field has more deviations and exhibits significant curvature in terms of its POS orientation (Yuen & Lazarian, 2020). Consequently, the observed structure, as well as the structure's topology, in \(p(x,y,v)\) is governed by both M\({}_{\text{A}}\) and \(\gamma\)(Hu et al., 2021).

To summarize succinctly, the thin channel maps \(p(x,y,v)\) from spectroscopic observations capture the anisotropy of MHD turbulence. This leads to the following important implications:

1. The intensity structures in \(p(x,y,v)\) align with the POS magnetic field.
2. The degree of anisotropy observed in these intensity structures is influenced by two distinct factors: M\({}_{\text{A}}\) and \(\gamma\). These factors contribute to the anisotropy: (a) \(\gamma\) introduces a projection effect that consequently decreases the anisotropy. (b) M\({}_{\text{A}}\) defines the magnetization level of the medium. A larger M\({}_{\text{A}}\) represents a weaker magnetic field, resulting in less pronounced anisotropy.
3. Additionally, changes in M\({}_{\text{A}}\) alter the topology of the magnetic field lines, as well as the observed intensity structure, manifesting itself as significant curvature.

The interconnection between magnetic field topology and M\({}_{\text{A}}\) is vital to extracting accurate 3D magnetic fields. A subtle change in the degree of anisotropy responds sensitively to variations in both M\({}_{\text{A}}\) and \(\gamma\), leading to a degeneracy. This degeneracy necessitates the introduction of an additional feature that is sensitive to M\({}_{\text{A}}\) or \(\gamma\) to solve for these parameters, and the topology of the magnetic field conveniently provides this required information.

Additionally, it is crucial to acknowledge that relying solely on anisotropy does not offer a clear distinction regarding the magnetic field's orientation along the LOS, specifically whether the field is directed towards or away from our observation point. Consequently, the value of \(\gamma\) is inherently restricted to a limited range between 0 and 90\({}^{\circ}\).

## 3 Numerical method

### Convolutional neural network (CNN)

To construct a deep neural network for the purpose of tracing the 3D magnetic field from a spectroscopic map, we adopt a CNN-based (LeCun et al., 1998) architecture. CNNs have demonstrated significant success in processing multidimensional data. The typical CNN architecture, as illustrated in Fig. 2, consists of initial layers comprising a stack of convolutional layers followed by pooling layers. To facilitate faster convergence during the network training process using backpropagation of the loss and enhance the stability of learning, we introduce a batch normalization layer following each convolution layer. After several iterations of convolution and pooling layers, we extract a compressed image feature, which is then processed by the fully connected layers to predict the desired properties. In the following, we introduce the core modules in the CNN architecture as well as the training procedure for the CNN network.

**Convolutional Layer:** Serving as the fundamental component of a CNN, the convolutional layer processes input data to produce feature maps (LeCun et al., 1989). In this layer, each neuron connects to a local region of the input feature map. This connection is achieved by applying a 2D convolutional kernel \(w_{I}\) to the input feature map. This process can be mathematically described as follows:

\[a_{I}=\sigma(w_{I}*h_{I-1}+b_{I}), \tag{8}\]

where \(h_{I-1}\) and \(a_{I}\) are the input and output feature map for the \(I-\)th convolutional layer, respectively, and \(w_{I}\) is the learnable convolution kernel, and \(*\) indicates the convolution operation. In addition, a learnable bias \(b_{I}\) is applied to the input feature map. To be more concrete,

\[a_{I}(x,y)=\sigma(\sum_{i=-k}^{k}\sum_{j=-k}^{k}w_{I}(i,j)h_{I-1}(x-i,y-j)+b_{I}( x,y)). \tag{9}\]
By applying the 2D convolution kernel \(w_{I}\in\mathbb{R}^{(2k+1)\times(2k+1)}\) to the input feature map \(h_{I-1}\in\mathbb{R}^{d^{\text{in}}\times d^{\text{in}}}\), we yield the output feature map with size \((d^{\text{in}}-k-1)\times(d^{\text{in}}-k-1)\). Here, \(d^{\text{in}}\) denotes the size of the input feature map and \(2k+1\) is the size of the convolution kernel. The resulting locally-weighted sum, once added to the learned bias, undergoes a non-linear transformation via the ReLU activation function \(\sigma(\cdot)\).

To constrain the number of parameters that need to be learned in our network, we generally use small kernel sizes. While each layer has a limited receptive field focusing on local features through the utilization of small convolutional kernels, stacking multiple layers allows for the gradual expansion of this receptive field. Consequently, the network becomes capable of capturing global features within the image as the depth increases.

**Batch Normalization Layer:** it is a technique frequently utilized in neural networks, playing a pivotal role in stabilizing them and hastening the convergence of the training loss during the backpropagation process (Ioffe and Szegedy, 2015). During each training iteration, it functions on a mini-batch of data. The layer normalizes each feature within the input data by centering its values around the mean and scaling based on the feature's standard deviation within the given batch. This normalization process is instrumental in mitigating the internal covariate shift -- a phenomenon where the distribution of inputs at each layer undergoes changes during training -- facilitating a more stable and efficient training process.

Following the normalization, batch normalization introduces two learnable parameters per feature: a scaling parameter and a shifting parameter. These parameters allow the network to learn the optimal scale and shift for the normalized values autonomously, providing the model with the flexibility to modify the normalization if it learns that such reversal or adjustment is beneficial for its predictive performance. These dynamic adjustments, enabled by the introduced parameters, imbue the network with a degree of adaptability, allowing it to fine-tune the transformations applied to the features as needed during the training.

**Pooling Layer:** following the detection of local features in the input feature maps by the convolution layer, a pooling layer is typically employed to merge similar local features into a singular feature (Sermanet et al., 2013). One common variant of the pooling layer is the _Max Pooling Layer_. This layer works by calculating the maximum value within a local patch of neurons and then outputting this maximum value as a single neuron. Importantly, the patches of input neurons for adjacent pooling units are shifted by more than one row or column, which effectively reduces the dimensionality of the feature representation. This process imparts the network with a degree of invariance to minor shifts and distortions in the input data, as it condenses the information in the feature maps while retaining the most salient features. This reduction not only helps in making the detection of features invariant to scale and orientation changes but also enhances computational efficiency by reducing the number of parameters and computations in the network.

**Fully Connected Layer:** After sequential operations that involve multiple convolutional layers and aggregation, the network derives a lower-dimensional compressed image feature map. Subsequently, this 2D feature map undergoes a transformation, being flattened into a 1D vector. The fully connected layer then processes this vector (Goodfellow et al., 2016). The role of the layer is critical, as it integrates the high-level reasoning of the features extracted and flattened previously. The mechanism involves applying learned weights and biases to this flattened vector to predict the final output. Mathematically, this operation can be represented as:

\[\mathbf{y}=\sigma(\mathbf{Wh+b}), \tag{10}\]

In this equation, \(\mathbf{h}\in\mathbb{R}^{d_{\text{in}}}\) represents the flattened, compressed image feature vector, and \(\mathbf{y}\in\mathbb{R}^{d_{\text{out}}}\) symbolizes the predicted result. Here, \(\mathbf{W}\in\mathbb{R}^{d_{\text{out}}\times d_{\text{in}}}\) and \(\mathbf{b}\in\mathbb{R}^{d_{\text{out}}}\) denote the learnable weights and biases for the fully connected layer, respectively. \(d^{\text{out}}\) represents the size of the output feature map. These weights and biases are integral to the layer's functionality, providing the means for it to learn and adapt during the training phase, ultimately allowing for the accurate prediction of the desired output from the input images.

**Network Training:** The trainable parameters within the CNN are optimized by adhering to a conventional neural network training methodology, where the mean-squared error of the 3D magnetic field prediction serves as the training loss for backpropagation, as outlined in the seminal work by Rumelhart et al. (1986). During the training process, we implement a strategy designed to enrich the diversity of the training dataset and consequently enhance the generalization capabilities of the deep neural network. Specifically, this involves augmenting the input images by subjecting them to random cropping operations, resulting in smaller patches of size \(22\times 22\) cells. Such augmentation introduces variability and randomness into the training data, which is instrumental in refining the network's ability to generalize from the training data to unseen data, thereby bolstering its predictive accuracy and robustness. In total, we generated \(\approx 1.7\times 10^{7}\) input \(22\times 22\)-cell maps, with 20% of them serving as a validation set, for each molecular species.

Figure 2: Architecture of the CNN-model. The input image is a \(22\times 22\) pixel map cropped from the thin velocity channel map. The network outputs the prediction of \(\phi\), \(\gamma\), or M\({}_{A}\).
### MHD simulations

The numerical simulations used in this study were executed using the ZEUS-MP/3D code (Hayes et al., 2006). We performed an isothermal simulation of a 10 pc cloud by solving the ideal MHD equations in an Eulerian frame under periodic boundary conditions:

\[\begin{split}&\partial\rho/\partial t+\nabla\cdot(\rho\mathbf{v})=0, \\ &\partial(\rho\mathbf{v})/\partial t+\nabla\cdot\left[\rho\mathbf{v}\mathbf{v}^ {T}+(c_{s}^{2}\rho_{+}\frac{B^{2}}{8\pi})\mathbf{I}-\frac{\mathbf{B}\mathbf{B}^{T}}{4\pi} \right]=\mathbf{f},\\ &\partial\mathbf{B}/\partial t-\nabla\times(\mathbf{v}\times\mathbf{B})=0,\\ &\nabla\cdot\mathbf{B}=0,\end{split} \tag{11}\]

where \(\mathbf{f}\) represents the stochastic forcing term used to drive turbulence. \(\rho\), \(\mathbf{v}\), and \(\mathbf{B}\) are mass density, velocity, and magnetic field, respectively. Given the isothermal equation of state, the sound speed \(c_{s}\) was held constant at approximately 187 m/s, corresponding to a gas temperature of 10 K. Purely turbulent scenarios were also considered, excluding the impact of self-gravity. Kinetic energy was solenoidoidally (i.e., the forcing term is divergence-free) injected at the wavenumber \(k=2\pi/l\approx 2\) (in the unit of \(2\pi/L_{\rm box}\), where \(L_{\rm box}\) is the length of simulation box) in Fourier space, where \(l\) is the length scale in real space, producing a Kolmogorov-like power spectrum. Turbulence was continuously stimulated until it reached a state of statistical saturation. The simulation was solved on a regular grid of \(792^{3}\) cells and the turbulence was numerically dissipated at scales of approximately 10 - 20 cells.

The simulations were initialized with a uniform density field and a magnetic field, with the initial mean magnetic field oriented along the y-axis. Furthermore, we rotated the simulation cubes so that the mean angle of inclination with respect to the LOS (or z-axis) reached \(90^{\circ}\), \(60^{\circ}\), and \(30^{\circ}\). The sonic Mach number, \(\mathrm{M}_{\mathrm{s}}=v_{\rm inj}/c_{s}\), and the Alfvenic Mach number, \(\mathrm{M}_{\mathrm{A}}=v_{\rm inj}/v_{A}\), characterize MHD turbulence simulations. To model different ISM conditions, we used a typical mean number density of 300 cm\({}^{-3}\) and varied the initial uniform magnetic field and the injected kinetic energy to obtain a range of \(\mathrm{M}_{\mathrm{A}}\) and \(\mathrm{M}_{\mathrm{s}}\) values. In this paper, we refer to the simulations in Tab. 1 by their model name or key parameters.

### Emission lines of \({}^{13}\)CO and C\({}^{18}\)O

We generate synthetic emission lines for two CO isotopologues: \({}^{13}\)CO (1-0) and C\({}^{18}\)O (1-0), following the procedures used in Hu & Lazarian (2021). This was achieved using the SPARX radiative transfer code (Hsieh et al., 2019). SPARX solves the radiative transfer equation (RTE) for finite cells, which means that it considers the emission from a homogeneous finite element. The equation of statistical equilibrium for molecular levels takes into account molecular self-emission, stimulated emission, and collisions with gas particles. Information on the distribution of molecular gas density with mean density \(\sim 300\) cm\({}^{-3}\) and LOS velocity was extracted from the MHD simulations mentioned above.

The fractional abundances of the CO isotopologues \({}^{13}\)CO(1-0) and C\({}^{18}\)O(1-0) were set at \(2\times 10^{-6}\) and \(1.7\times 10^{-7}\), respectively. We derive the \({}^{12}\)CO-to-H\({}_{2}\) ratio of \(1\times 10^{-4}\) from the cosmic value of C/H = \(3\times 10^{-4}\) and the assumption that 15% of C is in molecular form. The abundance of \({}^{13}\) CO is determined using a \({}^{13}\)CO/\({}^{12}\)CO ratio of 1/69, as indicated by Wilson (1999), giving a \({}^{13}\) CO / H\({}_{2}\) ratio of approximately \(2\times 10^{-6}\). Using a \({}^{12}\)CO/C\({}^{18}\)O ratio of 500, as given by Wilson et al. (2013), we obtained a C\({}^{18}\)O-to-H\({}_{2}\) ratio of \(1.7\times 10^{-7}\). When generating these synthetic emission lines, we specifically focused on the lowest-transition J = 1-0 of the CO isotopologues, with the Local Thermodynamic Equilibrium (LTE) satisfied.

### Training images

Our training input is a thin velocity channel map, \(p(x,y,v_{0})\), derived from either the \({}^{13}\)CO (1-0) or C\({}^{18}\)O (1-0) line, calculated from:

\[p(x,y,v_{0})=\int_{v_{0}-\Delta v/2}^{v_{0}+\Delta v/2}T_{\mathrm{e}}(x,y,v) dv, \tag{12}\]

where \(v_{0}\) is the velocity associated with the line's central peak, \(T_{\mathrm{e}}\) is the emission line's intensity, and \(\Delta v=\sqrt{\delta(v^{2})}\). Here \(\sqrt{\delta(v^{2})}\) is the velocity dispersion derived from the moment-1 map (velocity centroid map). The \({}^{12}\)CO line, a common diffuse cloud tracer, is not used in this work due to numerical limitations related to the saturation of the intensity of \({}^{12}\)CO in the channel centering at \(v_{0}\), which obliterates the spatial features of that channel (Hsieh et al., 2019). However, the CNN method could be extended to include wing channels centering at \(|v|<v_{0}\) to bypass this numerical saturation, a possibility we might explore in future work.2

Footnote 2: The use of wing channels has its own advantages through increasing the ratio of velocity to density fluctuations (Yuen et al., 2021; Hu et al., 2023).

We generate \(p(x,y,v_{0})\) for the full cloud, a region of \(792\times 792\) cells, then randomly segment \(p(x,y,v_{0})\) into \(22\times 22\)-cell subfields for input into the CNN model. The choice of \(22\times 22\)-cell avoids that the features fall into the numerical dissipation range, in which the

\begin{table}
\begin{tabular}{c c c c c c c} \hline Run & \(\mathrm{M}_{\mathrm{s}}\) & \(\mathrm{M}_{\mathrm{A}}\) & \(\mathrm{min}(\mathrm{M}_{\mathrm{A}}^{\mathrm{sub}})\) & \(\mathrm{max}\{\mathrm{M}_{\mathrm{A}}^{\mathrm{sub}}\}\) & \(\mathrm{min}[\mathrm{M}_{\mathrm{s}}^{\mathrm{sub}}]\) & \(\mathrm{max}\{\mathrm{M}_{\mathrm{s}}^{\mathrm{sub}}\}\) \\ \hline \hline A0 & 5.33 & 0.20 & 0.03 & 0.28 & 2.97 & 7.84 \\ A1 & 5.38 & 0.41 & 0.10 & 0.81 & 2.90 & 7.24 \\ A2 & 5.40 & 0.61 & 0.21 & 1.00 & 3.15 & 7.33 \\ A3 & 5.20 & 0.79 & 0.29 & 1.37 & 3.10 & 6.55 \\ A4 & 5.23 & 0.95 & 0.30 & 1.99 & 3.00 & 7.18 \\ A5 & 5.12 & 1.13 & 0.32 & 2.49 & 3.17 & 6.80 \\ A6 & 5.38 & 1.09 & 0.41 & 3.37 & 3.13 & 6.96 \\ A7 & 5.23 & 1.39 & 0.40 & 4.13 & 3.19 & 7.41 \\ A8 & 5.16 & 1.46 & 0.39 & 4.94 & 3.21 & 6.76 \\ A9 & 5.08 & 1.43 & 0.48 & 6.06 & 2.87 & 7.10 \\ \hline \end{tabular}
\end{table}
Table 1: \(\mathrm{M}_{\mathrm{s}}\) and \(\mathrm{M}_{\mathrm{A}}\) are the sonic Mach number and the Alfvénic Mach number calculated from the global injection velocity, respectively. \(\mathrm{M}_{\mathrm{A}}^{\mathrm{sub}}\) and \(\mathrm{M}_{\mathrm{s}}^{\mathrm{sub}}\) are determined using the local velocity dispersion calculated along each LOS in a \(22\times 22\) cell sub-field. The expressions ”\(\mathrm{min}[...]\)” and ”\(\mathrm{max}[...]\)” denote the minimum and maximum value averaged over each \(22\times 22\) cell sub-field within the corresponding simulation.
anisotropy of MHD turbulence is distorted by numerical diffusivity. In observation, the inertial range of MHD turbulence is much longer and the velocity channel map is not affected by the dissipation. The size of the sub-field, thus, could be smaller to achieve higher resolution. For each subfield, we also generate corresponding projected maps of \(\phi^{\rm sub}\), \(\gamma^{\rm sub}\), \(M_{A}^{\rm sub}\), and \(M_{s}^{\rm sub}\) as per the following:

\[\phi^{\rm sub}(x,y) =\arctan(\frac{\int B_{y}(x,y,z)dz}{\int B_{x}(x,y,z)dz}), \tag{13}\] \[\gamma^{\rm sub}(x,y) =\arccos(\frac{\int B_{z}(x,y,z)dz}{\int B(x,y,z)dz}),\] \[\mathrm{M}_{A}^{\rm sub} =\frac{v_{\rm inj}^{\rm los}\sqrt{4\pi\left\langle\rho\right\rangle _{\rm los}}}{\left\langle B\right\rangle_{\rm los}},\] \[\mathrm{M}_{\rm s}^{\rm sub} =\frac{v_{\rm inj}^{\rm los}}{c_{s}},\]

where \(B=\sqrt{B_{x}^{2}+B_{y}^{2}+B_{z}^{2}}\) is the total magnetic field strength, and \(B_{x}\), \(B_{y}\), and \(B_{z}\) are its \(x\), \(y\), and \(z\) components. \(\left\langle\rho\right\rangle_{\rm los}\) and \(\left\langle B\right\rangle_{\rm los}\) are the gas mass density and magnetic field strength averaged along the LOS. \(\mathrm{M}_{A}^{\rm sub}\) and \(\mathrm{M}_{\rm s}^{\rm sub}\) are defined using the local velocity dispersion for each LOS (i.e., \(v_{\rm inj}^{\rm los}\)), rather than the global turbulent injection velocity \(v_{\rm inj}\) used to characterize the full simulation. The ranges of \(\mathrm{M}_{A}^{\rm sub}\) and \(\mathrm{M}_{\rm s}^{\rm sub}\) averaged over the subfield in each simulation with different \(\gamma\) are listed in Tab. 1, while \(\gamma^{\rm sub}\) spans from 0 to 90\({}^{\circ}\). These values of \(\mathrm{M}_{A}^{\rm sub}\), \(\mathrm{M}_{\rm s}^{\rm sub}\), and \(\gamma^{\rm sub}\) cover typical physical conditions of diffuse molecular clouds (Hu & Lazarian, 2023c).

## 4 Results

### Numerical training and tests

Fig. 3 provides a visualization detailing the influence of \(\mathrm{M}_{A}\) and \(\gamma\) on the anisotropy of intensity structures within thin velocity channels. In scenarios where both \(\mathrm{M}_{A}\) and \(\gamma\) values are small, the intensity structures distinctly manifest as slender strips, extending in alignment with the POS magnetic fields. These structures are produced predominantly by the turbulent velocity (Lazarian & Pogosyan, 2000), as demonstrated in Hu et al. (2023). As \(\mathrm{M}_{A}\) increases, representing a weakening in the magnetic field, the MHD turbulence begins to more closely resemble isotropic hydrodynamical turbulence. This shift brings about a marked change in the topology of intensity structures, making them less anisotropic. Alternatively, when dealing with smaller values of \(\gamma\), which imply that magnetic fields are oriented more proximally to the LOS, the inherent anisotropy is subdued due to the projection effect. Comparing \({}^{13}\)CO and \(\mathrm{C}^{18}\mathrm{O}\), \(\mathrm{C}^{18}\mathrm{O}\) is more sensitive to denser gas, so its associated intensity structures exhibit distinct characteristics. Despite these differences, the underlying physical principle of anisotropic MHD turbulence remains the same, suggesting \(\mathrm{M}_{A}\) and \(\gamma\) continue to shape the observed structural formations.

Fig. 4 provides a comparative visualization between the actual 3D magnetic fields and those predicted through the utilization of the trained CNN model with \({}^{13}\)CO. This comparison is framed within two distinct conditions: sub-Alfvenic (simulation with \(\langle\mathrm{M}_{A}\rangle\approx 0.5\) and \(\left\langle\gamma\right\rangle\approx 90^{\circ}\)) and super-Alfvenic (simulation with \(\langle\mathrm{M}_{A}\rangle\approx 2.0\) and \(\left\langle\gamma\right\rangle\approx 30^{\circ}\)). Within these settings, the mean projected total Alfven Mach number on the POS is given as \(\langle\mathrm{M}_{A}\rangle\approx 0.5\) for sub-Alfvenic conditions and \(\langle\mathrm{M}_{A}\rangle\approx 2.0\) for super-Alfvenic ones.

The visual segment displayed in Fig. 4 is constructed from the POS
magnetic field's position angle, \(\phi\), and the inclination angle, \(\gamma\), with a superimposed color representation signifying the projected M\({}_{A}\). Upon comparison with the intrinsic magnetic field embedded within the simulation, a noteworthy observation is the alignment between the orientations of the CNN-predicted 3D magnetic field and the actual field, evident under both sub-Alfvenic and super-Alfvenic conditions. In the sub-Alfvenic case, the CNN-predicted M\({}_{A}\) is slightly larger (by \(\approx 0.1-0.2\)) than the actual values. Conversely, in the super-Alfvenic scenario, the predicted value is somewhat smaller, with a deviation ranging from \(\approx 0.5-1.0\). Another example with \(\langle\)M\({}_{A}\rangle\approx 0.15\) and \(\langle\gamma\rangle\approx 60^{\circ}\), is presented in Appendix A. Although this simulation shows an anisotropy degree similar to the case with \(\langle\)M\({}_{A}\rangle\approx 0.5\) and \(\langle\gamma\rangle\approx 90^{\circ}\), the CNN model effectively resolves the degeneracy in the correlation of the anisotropy degree with \(\gamma\) and M\({}_{A}\) (see SS 2), successfully recovering the 3D magnetic field (see Fig. A1). It should be noted that the predicted M\({}_{A}\) is still overestimated by approximately 0.1 - 0.2.

Fig. 5 offers a similar visual comparison but focuses on the C\({}^{18}\)O line. This line is generally recognized as denser tracers compared to \({}^{13}\)CO. Despite these differences in tracer density, the CNN predictions for C\({}^{18}\)O lines maintain a general alignment with the actual 3D magnetic fields observed within the simulations. Moreover, there is less significant overestimation and underestimation in the CNN-predicted M\({}_{A}\).

Figs. 6 and 7 present 2D histograms illustrating the correspondence between CNN predictions--\(\phi^{\rm CNN}\), \(\gamma^{\rm CNN}\), and M\({}_{A}^{\rm CNN}\)--and actual values obtained from two test simulations, A2 and A6. In sub-Alfvenic cases for both \({}^{13}\)CO and C\({}^{18}\)O molecules, we observe a close alignment between the CNN predictions and the real values. The scatter of the predictions, which includes \(\phi^{\rm CNN}\), \(\gamma^{\rm CNN}\), and M\({}_{A}^{\rm CNN}\), demonstrates a small deviation from the actual values, tightly cornegating near the one-to-one reference line. This minimal deviation suggests that the CNN model offers a high degree of accuracy and reliability when operating under sub-Alfvenic conditions.

However, the scenario is a bit different in super-Alfvenic cases. Here, the scatter is noticeably more widespread, indicating that deviations from the real values increase in these conditions. The \(\phi^{\rm CNN}\) predictions, in particular, show a tendency for both overestimation and underestimation. In contrast, the \(\gamma^{\rm CNN}\) predictions are primarily characterized by overestimations, a trend that is especially prominent in cases involving C\({}^{18}\)O molecules. Meanwhile, the scatter related to the M\({}_{A}^{\rm CNN}\) predictions is distributed more uniformly around the reference line.

This suggests predicting the 3D magnetic field under super-Alfvenic conditions is more challenging with higher uncertainty.

Figure 4: An comparison of the CNN-predicted 3D magnetic fields using \({}^{13}\)CO in sub-Alfven (top, \(\langle\)M\({}_{A}\rangle\approx 0.5\) and \(\langle\gamma\rangle\approx 90^{\circ}\)) and super-Alfven (bottom, \(\langle\)M\({}_{A}\rangle\approx 2.0\) and \(\langle\gamma\rangle\approx 30^{\circ}\)) conditions. Each magnetic field segment is constructed by the POS magnetic field’s position angle (i.e., \(\phi\)) and the inclination angle \(\gamma\). Note that the magnetic field obtained is the projection along the LOS and averaged over 132\(\times\)132 pixels for visualization purposes. The third axis of the LOS is for 3D visualization purposes and does not provide distance information here. The total intensity map \(I\) is placed on the POS, i.e., the \(x-y\) plane.
In these environments, the magnetic field exerts a weaker influence, and the turbulence observed more closely resembles that of hydrodynamic turbulence, thereby complicating the prediction process. Enhancing prediction accuracy is feasible through two strategies. First, it is possible to further refine and optimize the CNN model to improve its adaptability and responsiveness to the unique features of super-Alfvenic MHD turbulence. For instance, Peek & Burkhart (2019) put forth a CNN model designed specifically to differentiate between sub-Alfvenic and super-Alfvenic turbulence. This model, with its specialized focus, offers a promising avenue for enhancing the accuracy of predictions in super-Alfvenic environments. Second, enrich the data set to train the CNN model. By incorporating a broader and more diverse range of images, the model can be exposed to a wider array of scenarios and conditions, thereby reducing uncertainty and improving its ability to make accurate predictions across different environments and conditions.

Figs. 8 and 9 plot the histograms of the deviation between the CNN-predicted and the actual 3D magnetic field. We calculate the absolute difference between \(\phi^{\text{CNN}}\) and \(\phi\), between \(\gamma^{\text{CNN}}\) and \(\gamma\), and between \(\text{M}^{\text{CNN}}_{\text{A}}\) and \(\text{M}_{\text{A}}\), respectively. These differences are denoted as \(\sigma_{\phi}\), \(\sigma_{\gamma}\), and \(\sigma_{\text{M}_{\text{A}}}\). In the sub-Alfvenic scenarios, we observed that the distributions of \(\sigma_{\phi}\) and \(\sigma_{\gamma}\) are relatively condensed, primarily falling within the 0 to 20\({}^{\circ}\) range. This concentration indicates a close alignment between the CNN predictions and the actual values in sub-Alfvenic environments, suggesting that the CNN model performs with high precision in these conditions. However, as \(\langle\text{M}_{\text{A}}\rangle\) increases, the distributions of \(\sigma_{\phi}\) and \(\sigma_{\gamma}\) broaden, spanning a more extensive range from 0 to 60\({}^{\circ}\). This dispersion is indicative of larger deviations between predicted and actual values under these conditions, implying that the CNN model may face challenges in accurately capturing the magnetic field dynamics when \(\langle\text{M}_{\text{A}}\rangle\) increases.

Examining specific molecules, for \({}^{13}\)CO under sub-Alfvenic conditions, the median deviation values are relatively low: \(\sigma_{\phi}=3.26^{\circ}\), \(\sigma_{\gamma}=2.98^{\circ}\), and \(\sigma_{\text{M}_{\text{A}}}=0.16\). In contrast, under super-Alfvenic conditions, these values increase to \(12.32^{\circ}\), \(9.08^{\circ}\), and \(1.1\), respectively, highlighting an increase in prediction deviation as the environment transitions from sub- to super-Alfvenic. Similarly, for C\({}^{18}\)O, the median deviation values are \(2.22^{\circ}\), \(3.20^{\circ}\), and \(0.16\) under sub-Alfvenic conditions and \(12.08^{\circ}\), \(13.60^{\circ}\), and \(1.36\) under super-Alfvenic scenarios, underlining a consistent trend of increased deviation in super-Alfvenic environments across different molecules.

### Observational prediction

For the observational tests, our target is the nearby L1478 cloud. We utilized \({}^{13}\)CO spectral line from a previous study Lewis et al. (2021). The data has a beam resolution of 38\({}^{\text{n}}\) and was regrid to a pixel resolution of 10\({}^{\text{n}}\), while achieving a velocity resolution of 0.3 km s\({}^{-1}\). The 1D velocity dispersion \(\sigma_{\text{V}}\) of the \({}^{13}\)CO line was reported within the range of 0.40 - 0.70 km s\({}^{-1}\)(Lewis et al., 2021). Assuming an isotropic velocity dispersion in 3D and uniform temperature of 10 K (corresponding to an isothermal sound speed of \(c_{s}\sim 0.187\) km s\({}^{-1}\),

Figure 5: Same as Fig. 4, but for C\({}^{18}\)O.
see Hu et al. 2021b), we find the sonic Mach number \(\mathrm{M}_{s}=\sqrt{3}\sigma_{\nu}/c_{s}\) ranges from 3.69 to 6.45, falling into the parameter regimes in our numerical simulations. With these refined data, we applied our aderply trained CNN model to the \({}^{13}\)CO channel map, aiming to predict the key 3D magnetic field parameters, denoted as \(\phi^{\mathrm{CNN}}\), \(\gamma^{\mathrm{CNN}}\), \(\mathrm{M}_{A}^{\mathrm{CNN}}\).

For the purpose of validating the results yielded through our CNN application, we engaged in a comparative analysis with POS magnetic field orientations as predicted through Planck 353 GHz polarization data. The data harnessed for this comparative process was drawn from the third Public Data Release (DR3), provided by Planck's

Figure 6: 2D histogram of the \({}^{13}\)CO CNN-predictions, i.e., \(\phi^{\mathrm{CNN}}\) (left), \(\gamma^{\mathrm{CNN}}\) (middle), and \(\mathrm{M}_{A}^{\mathrm{CNN}}\) (right) and the corresponding actual values in simulation (Top: sub-Alfvén, \((\mathrm{M}_{A})\approx 0.5\) and \(\langle\gamma\rangle\approx 90^{\circ}\). Bottom: super-Alfvén, \((\mathrm{M}_{A})\approx 2.0\) and \(\langle\gamma\rangle\approx 30^{\circ}\)). The dashed reference line represents the ideal scenario, where the predicted values and actual values match perfectly.

Figure 7: Same as Fig. 6, but for \(\mathrm{C}^{18}\mathrm{O}\).
High-Frequency Instrument (Planck Collaboration et al., 2020). The POS magnetic field orientation was inferred from Stokes parameters \(Q\) and \(U\) converted to IAU convention from HEALPix using the equation: \(\phi^{\rm Planck}=\frac{1}{2}\tan^{-1}(-U,Q)+\pi/2\). To enhance the signal-to-noise ratio, we smoothed the Stokes parameter maps from an angular resolution of \(5^{\prime}\) to \(10^{\prime}\) using a Gaussian kernel.

As presented in Fig. 10, a remarkable alignment between the magnetic field orientations as predicted by both the CNN model and the Planck polarization data is observed, while we notice the difference is apparent in the northeast clump (see the zoom-in plot in Fig. 10). To quantify the agreement between CNN-prediction and polarization, we utilize the Alignment Measure (AM; Gonzalez-Casanova &

Figure 8: Histograms of difference in CNN-predicted \(\phi^{\rm CNN}\) (left), \(\gamma^{\rm CNN}\) (middle), and \(\rm M_{A}^{\rm CNN}\) (right) and the actual values in simulations using \({}^{13}\)CO.

Figure 10: Comparison of the POS magnetic fields predicted by CNN-\({}^{13}\)CO (red segment) for the L1478 cloud and inferred from Planck polarization (blue segment). The background image is the integrated \({}^{13}\)CO intensity map.

Figure 9: Same as Fig. 8, but for \(\rm C^{18}\)O.
Lazarian, 2017), expressed as:

\[\mathrm{AM}=\langle\cos(2\theta_{t})\rangle, \tag{14}\]

here \(\theta_{t}\) is the relative angle between the two measurements. An AM value of \(\approx\) 0.94 confirms the CNN-prediction has an excellent agreement with Planck polarization 3, corresponding to an overall deviation of \(\approx\) 10\({}^{\circ}\).

Footnote 3: AM = 1 implies a perfect parallel alignment, while -1 indicates perpendicularity.

A noteworthy advantage of our CNN model over traditional polarization methodologies is its ability to trace the 3D magnetic fields. This is achieved through the model's predictions regarding \(\gamma\) and M\({}_{A}\). These predictions are summarised in histograms within Fig.11. According to the histograms, the median \(\gamma\) and M\({}_{A}\) of the L1478 cloud are estimated at \(\approx\) 76\({}^{\circ}\) and \(\approx\) 1.07, respectively. These measurements suggest that the L1478 is a trans-Alfvenic could. In this state, there is an equilibrium between magnetic and turbulent kinetic energies within the cloud. The parameters derived from the CNN application have been instrumental in creating the first-ever 3D magnetic field map for L1478, which can be viewed in Fig. 12.

## 5 Discussion

### Comparison with earlier studies

The realm of exploring magnetic fields within the ISM through CNN is experiencing swift advancements. As a pilot study presented by Xu et al. (2023), the Convolutional Approach to Structure Identification-3D (CASI-3D) model was employed to map the 2D POS magnetic field orientation. This is achieved similarly by using the velocity channel maps obtained from spectroscopic observations. The underlying physics principle is still founded on the anisotropic MHD turbulence. The training process underpinning this approach uses the emission lines of \({}^{12}\)CO and \({}^{13}\)CO (J = 1 - 0), generated through the RADMC-3D code (Dullemond et al., 2012).

In this study, we introduce a new CNN model. This advanced model is designed with the aim of predicting not merely the orientation \(\phi\) of the POS magnetic field but extends to encompass the angle of field inclination, \(\gamma\), as well as the total Alfven Mach number M\({}_{A}\). This approach allows the construction of 3D magnetic field vectors. For training the CNN model, we have utilized emission lines from \({}^{13}\)CO and C\({}^{18}\)O (J = 1 - 0), with data generated from the SPARX code (Hsieh et al., 2019).

We quantify the uncertainty of our CNN-predicted \(\phi\) and \(\gamma\). We found that the median value and the dispersion of uncertainty for

Figure 11: Histograms of CNN-predicted (as well as Planck measured) \(\phi^{\mathrm{CNN}}\) (left), defined east from the north, \(\gamma^{\mathrm{CNN}}\) (middle), and M\({}_{A}^{\mathrm{CNN}}\) (right).

Figure 12: An visilaration of the CNN-predicted 3D magnetic fields using \({}^{13}\)CO for the L1478 cloud. Each magnetic field segment is constructed by the position angle of the POS magnetic field (i.e., \(\phi\)) and the inclination angle \(\gamma\). Note that the magnetic field obtained is the projection along the LOS and averaged over 12\(\times\)12 pixels for visualization purposes. The third axis of the LOS is for 3D visualization purposes and does not provide distance information here. The total intensity map \(I\) is placed on the POS, i.e., the \(l-b\) plane.
C\({}^{18}\)O are approximately \(\sim 2.22^{\circ}\) and \(\sim 3.20^{\circ}\) under sub-Alfv'enic conditions (\(\langle\mathrm{M_{A}}\rangle\approx 0.5\)). These values shift to \(\sim 12.08^{\circ}\) and \(\sim 13.60^{\circ}\) undersuper-Alfvenic conditions (\(\langle\mathrm{M_{A}}\rangle\approx 2.0\)). When compared to the CASI-3D model, our CNN model demonstrates higher accuracy, as CASI-3D exhibits a median uncertainty of \(\sim 6.2^{\circ}\) and \(\sim 18.4^{\circ}\) under comparable sub-Alfvenic and super-Alfvenic conditions, respectively. Through the application of our CNN model to the L1478 molecular cloud, we successfully constructed the first 3D magnetic field map. The corresponding CNN-predicted POS magnetic field orientation shows remarkable alignment with that inferred from Planck 353 GHz polarization data.

It is crucial to acknowledge that despite the differences inherent between the CNN models used by Xu et al. (2023) and this study, the fundamental concept of utilizing spectroscopic channel maps for magnetic field investigation remains the same: (1) the intensity distribution observable in thin channel maps is predominantly influenced by turbulent velocity statistics, as outlined in (Lazarian and Pogosyan, 2000; Kandel et al., 2016; Hu et al., 2023); and (2) these channel maps capture the anisotropy intrinsic to MHD turbulence, thereby revealing the orientation of the POS magnetic field (Lazarian and Yuen, 2018; Hu et al., 2023). A crucial insight was provided by Hu et al. (2021c), highlighting that the degree of anisotropy in channel maps, as well as the magnetic field topology, is regulated by both the \(\gamma\) and the M\({}_{A}\). These are parameters that can be extracted efficiently using the CNN approach 4. Thus, drawing upon these foundational theoretical studies, we propose the use of the CNN model as an efficient tool for tracing 3D magnetic fields, providing convincing physical reasons for interpreting its feasibility.

Footnote 4: Note that while the anisotropy and magnetic field topology, that are sensitive to \(\gamma\) and the M\({}_{A}\), are the most apparent features in channel maps, it is also possible the CNN extracts additional features to facilitate the prediction.

### Synergy with other methods

Our newly proposed CNN model stands as a powerful complement to existing methodologies in the field. One notable technique, which involves utilizing polarized dust emission, has proven effective in tracing the 3D magnetic field orientation within diffuse clouds, where dust grains are perfectly aligned with magnetic fields (Chen et al., 2019; Hu and Lazarian, 2023a, c). However, this technique may encounter limitations within dense cloud environments - for example, those observed through tracing by C\({}^{18}\)O - where dust grains might not maintain perfect alignment (Lazarian, 2007; Andersson et al., 2015). This loss of alignment, resulting in a phenomenon known as the polarization hole (Seifried et al., 2019; Pattle et al., 2019; Hoang et al., 2021), introduces uncertainties when tracing 3D magnetic fields through polarized dust emission techniques.

Unlike these traditional approaches, the CNN approach remains immune to the effects of the polarization hole. When the CNN model is supplied with emission lines from dense tracers like C\({}^{18}\)O, HNC, and NH\({}_{3}\), it proves highly adept at probing the 3D magnetic fields present within dense clouds effectively. Nonetheless, it's important to consider that within these dense cloud environments, the forces of self-gravity can become a significant factor. This gravitational influence might induce alterations in the anisotropy observed within channel maps (Hu et al., 2020b). Therefore, it becomes imperative to input the CNN model with carefully selected numerical simulations before applying it to observational data to ensure accurate and reliable results.

Furthermore, it should be noted that the inclination angle predicted by the CNN model is inherently limited to the range of [0, 90\({}^{\circ}\)]. This limitation arises because the anisotropy within channel maps alone cannot definitively discern whether the magnetic field is oriented towards or away from the observer. However, recent advancements in the field, particularly in Faraday rotation measurements within molecular clouds (Tahani et al., 2019, 2022), offer promising avenues to resolve this degeneracy.

Another relevant method worth discussing is the Velocity Gradient Technique (VGT; Gonzalez-Casanova and Lazarian, 2017; Lazarian and Yuen, 2018; Hu et al., 2018). Like our proposed CNN approach, the VGT is a technique that traces magnetic fields using spectroscopic observations. Importantly, both the CNN approach and VGT share a foundational physical principle: they rely on the anisotropy of MHD turbulence observed within thin channel maps. With VGT having undergone extensive and rigorous testing (Hu et al., 2019; Lu et al., 2020; Hu et al., 2021; Liu et al., 2022; Alina et al., 2022; Liu et al., 2023; Hu and Lazarian, 2023b; Tram et al., 2023; Schmaltz et al., 2023), it is established as an excellent benchmark for evaluating the accuracy of CNN models, especially in situations where polarization measurements are not readily available. This benchmarking is crucial when CNNs are deployed for tracing 3D Galactic Magnetic Fields, highlighting the important comparative and complementary roles these techniques play in advancing our understanding of magnetic fields in various astrophysical contexts.

### Prospects of the CNN method

In the present study, we introduced a CNN model adept at predicting 3D magnetic fields within molecular clouds, utilizing spectroscopic observations of molecular gas. However, the potential applications of this CNN method extend far beyond, encompassing various astrophysical environments and contexts, including neutral hydrogen (HI) regions, ionized gas, the Central Molecular Zone (CMZ), external galaxies, and supernova remnants. In the following sections, we outline several promising applications of this methodology.

#### 5.3.1 3D Galactic Magnetic Fields

A deep and comprehensive understanding of the 3D Galactic Magnetic Field (GMF; Jansson and Farrar, 2012) is paramount for addressing a host of astrophysical inequities. These include identifying the origins of ultra-high energy cosmic rays (Farrar, 2014; Farrar and Sutherland, 2019) and refining models of Galactic foreground polarization (Kovetz and Kamionkowski, 2015; Planck Collaboration et al., 2016).

Recent research indicates that thin channel maps of HI successfully capture the anisotropy inherent in MHD turbulence (Lazarian and Pogosyan, 2000; Lazarian and Yuen, 2018; Lu et al., 2020; Hu et al., 2023). Consequently, applying the CNN to HI channel maps constitutes a viable strategy for mapping 3D GMFs. Past efforts aimed at modeling the foreground polarization with HI primarily focused on mapping the POS magnetic field orientation (Clark and Hensley, 2019; Lu et al., 2020; Hu et al., 2020). These endeavors largely neglected the crucial depolarization factor, the inclination angle. However, the advent of sophisticated multi-phase HI simulations (Ho et al., 2021) has made it possible to train the CNN model for accurate predictions of 3D GMFs, yielding more realistic models of the foreground polarization.

Our primary goal in this paper is to explore the magnetic fields of molecular clouds, for which the isothermal approximation is applicable. Multi-phase HI requires separate training of the neural network.
For multi-phase HI where cooling and heating play a significant role, our general approach remains valid: intensity features/striations within channel maps continue to elongate along the POS magnetic field orientation. This is supported by several studies Lazarian and Yuen (2018); Clark and Hensley (2019); Lu et al. (2020); Hu et al. (2020, 2023). These intensity features/striations are also regulated by the Alfven Mach number (M\({}_{A}\)) and the projection effect associated with the inclination angle. However, additional physics, such as thermal instability, could modify the observed anisotropy, for instance, potentially leading to a smaller aspect ratio (Ho et al., 2023). The corresponding study employing our approach for multiphase HI will be provided elsewhere.

#### 5.3.2 3D Magnetic Fields in CMZ and external galaxies

Understanding the magnetic fields within cold molecular gas is essential for deciphering the processes of formation and fueling of Seyfert nuclei. Recent measurements of magnetic fields within the CMZ and in other Seyfert galaxies have been conducted using various techniques. These include far-infrared polarization observations from instruments like SOFIA/HAWC+ (Lopez-Rodriguez et al., 2021), JCMT (Pattle et al., 2021), and ALMA (Lopez-Rodriguez et al., 2020), as well as employing the VGT (Hu et al., 2022, 2023). However, these approaches primarily yield the POS magnetic field orientation, falling short of providing a comprehensive 3D perspective. Nevertheless, the successful application of VGT confirms the viability of using anisotropy in molecular emission channel maps as a tracer for magnetic fields in these environments. For instance, Hu et al. (2022) derived a POS magnetic field map surrounding Sgr A* using the [Ne II] emission line and Paschen-\(\alpha\) image observed with the Hubble Space Telescope (HST). Given these advances, extending the CNN methodology to incorporate optical/near-infrared observations from instruments like the HST and the James Webb Space Telescope is a feasible and promising approach for predicting 3D magnetic fields in both the Galactic center and external galaxies.

### Obtaining the full 3D magnetic field vector

3D magnetic fields, encompassing both orientation and strength, play a pivotal role in comprehending key astrophysical phenomena. These include processes such as star formation (Mestel, 1965; Mac Low and Klessen, 2004; McKee and Ostriker, 2007; Lazarian et al., 2012; Federrath and Klessen, 2012; Hu et al., 2021), the effects of stellar feedback (Pattle et al., 2022; Liu et al., 2023), as well as the acceleration and propagation of cosmic rays (Fermi, 1949; Jokipii, 1966; Yan and Lazarian, 2002; Xu and Yan, 2013; Xu and Lazarian, 2020; Hu et al., 2022; Beattie et al., 2022; Lazarian and Xu, 2023). Traditionally, to obtain the strength of these fields, the Davis-Chandrasekhar-Fermi (DCF) method is employed, which typically combines dust polarimetry with spectroscopic observations (see Davis, 1951; Chandrasekhar and Fermi, 1953). However, this often proves insufficient. the DCF method gives only the POS magnetic field strength, while the component along the LOS is missing. Other limitations of the DCF method have also been thoroughly dissected in the literature (Skadilidis et al., 2021; Lazarian et al., 2022; Chen et al., 2022; Liu et al., 2022).

In light of this, an alternative approach has been proposed: the use of the Alfven Mach number M\({}_{A}\) with the sonic Mach number M\({}_{s}\) to derive the magnetic field's strength (Lazarian et al., 2020). This method, aptly termed MM2, can be used to obtain the total strength, particularly since the vital term M\({}_{A}\) is readily available with the CNN approach proposed in this study. The sonic Mach number M\({}_{s}\) can be procured either directly via spectroscopic line broadening or by leveraging a CNN approach similar to our current study. Coupled with the 3D magnetic field orientation, this equips us with the necessary tools to construct a 3D magnetic field vector.

## 6 Summary

In this study, a CNN model was designed for the intricate task of probing 3D magnetic fields within molecular clouds. This model is not confined to determining the POS magnetic field orientation but extends its capabilities to accurately ascertain the field's inclination angle and the total Alfven Mach number, offering a more comprehensive understanding of the magnetic field in the observed regions. We summarize our major results below:

1. We developed a CNN model for probing the 3D magnetic fields, including the POS magnetic field orientation, inclination angle, and total Alfven Mach number.
2. The CNN model was trained using synthetic \({}^{13}\)CO and C\({}^{18}\)O (J = 1 - 0) emission lines, encompassing a range of conditions from sub-Alfvenic to super-Alfvenic. We quantified the uncertainties associated with the trained CNN model's predictions. Our findings revealed that the uncertainties are less than 5\({}^{\circ}\) for both \(\phi\) and \(\gamma\), and are smaller than 0.2 for M\({}_{A}\) under sub-Alfvenic conditions (with M\({}_{A}\approx 0.5\)). Under super-Alfvenic conditions (with M\({}_{A}\approx 2.0\)), the uncertainties increased slightly but remained below 15\({}^{\circ}\) for \(\phi\) and \(\gamma\), and were around 1.5 for M\({}_{A}\).
3. We implemented our trained CNN model to analyze the molecular cloud L1478. The CNN-predicted POS magnetic field orientation exhibited remarkable agreement with orientations inferred from Planck 353 GHz polarization data, with a marginal global difference of approximately 10\({}^{\circ}\).
4. This study facilitated the construction of the first 3D magnetic field map for the L1478 cloud. Through our analysis, we found that the cloud's global inclination angle is approximately 76\({}^{\circ}\), while the global total Alfven Mach number is close to 1.07.
5. We discussed the potential applications and future prospects of the CNN approach. Particularly, we discussed the feasibility and potential of utilizing the CNN model for predicting 3D Galactic Magnetic Fields. We also considered its application for understanding 3D magnetic fields in the CMZ and external galaxies.

## Acknowledgements

Y.H. and A.L. acknowledge the support of NASA ATP AAH7546, NSF grants AST 2307840, and ALMA SOSPADA-016. Financial support for this work was provided by NASA through award 09_0231 issued by the Universities Space Research Association, Inc. (USRA). This work used SDSC Expansive CPU at SDSC through allocations PHY230032, PHY230033, PHY230091, and PHY230105 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. Y.H. acknowledges the very kind computational and technical support of Bowen Cao.

## Data Availability

The data underlying this article will be shared on reasonable request to the corresponding author.