# MC-BLOS111software currently available at
https://github.com/MehrnooshTahani/MC-BLOS: Determination of the Line-of-Sight
Component of Magnetic Fields Associated with Molecular Clouds

Mehrnoosh Tahani Banting Fellow: Kavli Institute for Particle Astrophysics &
Cosmology (KIPAC), Stanford University, Stanford, CA 94305, USA  Department of
Physics, Stanford University, Stanford, CA 94305, USA  Dominion Radio
Astrophysical Observatory, Herzberg Astronomy and Astrophysics Research
Centre, National Research Council Canada, P. O. Box 248, Penticton, BC V2A 6J9
Canada  John Ming Ngo  Department of Physics, University of Alberta, Edmonton,
Alberta T6G 2E1, Canada  Dominion Radio Astrophysical Observatory, Herzberg
Astronomy and Astrophysics Research Centre, National Research Council Canada,
P. O. Box 248, Penticton, BC V2A 6J9 Canada  Jennifer Glover Department of
Physics, McGill University, Montréal, Quebec H3A 2T8, Canada  Dominion Radio
Astrophysical Observatory, Herzberg Astronomy and Astrophysics Research
Centre, National Research Council Canada, P. O. Box 248, Penticton, BC V2A 6J9
Canada  Ryan Clairmont  Department of Physics, Stanford University, Stanford,
CA 94305, USA  Gabriel M. Zarazua  Department of Physics and Astronomy, San
Francisco State University, 1600 Holloway Ave., San Francisco, CA 94132, USA
René Plume  Department of Physics and Astronomy, University of Calgary, 2500
University Drive NW, Calgary, Alberta T2N 1N4, Canada

###### Abstract

In recent years a number of surveys and telescopes have observed the plane-of-
sky component of magnetic fields associated with molecular clouds. However,
observations of their line-of-sight magnetic field remain limited. To address
this issue, Tahani et al. (2018) developed a technique based on Faraday
rotation. The technique incorporates an ON-OFF approach to identify the
rotation measure induced by the magnetic fields associated with the cloud. The
upcoming abundance of Faraday rotation observations from the Square Kilometer
Array and its pathfinders necessitates robustly-tested software to
automatically obtain line-of-sight magnetic fields of molecular clouds. We
developed software, called MC-BLOS (Molecular Cloud Line-of-Sight Magnetic
Field), to carry out the technique in an automated manner. The software’s
input are Faraday rotation of point sources (extra-galactic sources or
pulsars), extinction or column density maps, chemical evolution code results,
and a text/CSV file, which allows the user to specify the cloud name or other
parameters pertaining to the technique. For each cloud, the software invokes a
set of predefined initial parameters such as density, temperature, and
surrounding boundary, which the user can modify. The software then runs the
technique automatically, outputting line-of-sight magnetic field maps and
tables (including uncertainties) at the end of the process. This automated
approach significantly reduces analysis time compared to manual methods. We
have tested the software on previously-published clouds, and the results are
consistent within the reported uncertainty range. This software will
facilitate the analysis of forthcoming Faraday rotation observations, enabling
a better understanding of the role of magnetic fields in molecular cloud
dynamics and star formation.

††journal: ApJ††software: Astropy (Astropy Collaboration et al., 2013),
Matplotlib (Hunter, 2007), Math (Van Rossum, 2020), Numpy (Harris et al.,
2020), Pandas (pandas development team, 2020), Scikit-learn (Pedregosa et al.,
2011).

##  1 Introduction

Mapping the three-dimensional (3D) structure of magnetic fields is crucial for
understanding the role of magnetic fields in the formation and evolution of
molecular clouds and their subsequent star formation. While magnetic fields
can hinder star formation by providing additional support against
gravitational collapse (Seifried & Walch, 2015) or when combined with feedback
mechanism such as stellar outflows (Krumholz & Federrath, 2019), they are also
necessary for regulating core collapse through magnetic braking (Pudritz &
Ray, 2019). Recent studies have shown that the relative orientation between
magnetic fields and density structures such as molecular clouds can provide
insights into their formation and evolution (Planck Collaboration et al.,
2016; Arzoumanian et al., 2021; Pattle et al., 2023). However, to fully
understand the role of magnetic fields in molecular cloud dynamics and star
formation, it is essential to observe and analyze their 3D structure,
including both the plane-of-sky ($B_{\rm{POS}}$) and line-of-sight
($B_{\rm{LOS}}$) magnetic field components (Tahani et al., 2019, 2022a, 2022b;
Tahani, 2022).

Far-infrared (FIR) and sub-millimeter observations, such as those by Planck
(e.g., Planck Collaboration et al., 2015) or James Clerk Maxwell telescope
(e.g., Arzoumanian et al., 2021; Hwang et al., 2022; Ching et al., 2022;
Tahani et al., 2023), have successfully used dust emission polarization in
large surveys to map $B_{\rm{POS}}$. In this technique the measured
polarization is perpendicular to the magnetic field in the plane of the sky
(Pattle & Fissel, 2019), and $B_{\rm{POS}}$ is averaged along the lines sight
(mass weighted; Seifried et al., 2019). Observations using dust extinction
polarization at near-infrared (NIR) and optical wavelengths, referred to as
the starlight polarization technique, combined with stellar distances obtained
by Gaia, enable us to obtain a tomography of the two-dimensional magnetic
field along the line of sight (tomographic $B_{\rm{POS}}$; Doi et al., 2021;
Panopoulou et al., 2023; Doi et al., 2024). This technique involves measuring
the polarization of starlight that has passed through the interstellar medium,
where the interstellar dust grain have aligned their long axis perpendicular
to magnetic field lines (Lazarian & Hoang, 2007). In this case the measured
polarized light is parallel to the plane-of-sky magnetic field orientation. By
combining these measurements with stellar distances, it is possible to
construct a tomographic view of the $B_{\rm{POS}}$ orientation along the line
of sight. However, to study the true 3D magnetic fields, observing
$B_{\rm{LOS}}$ is crucial.

Zeeman splitting is sensitive to the line-of-sight component of magnetic
fields and is the most direct method for measuring the strength of magnetic
fields (Crutcher & Kemball, 2019), but it requires long telescope integration
times and strong magnetic fields (Crutcher, 2012; Robishaw, 2008). Faraday
rotation is a more widely applicable technique for determining $B_{\rm{LOS}}$,
as it does not have the same limitations as Zeeman splitting. However, it is
less direct (more model dependent) for inferring the field strengths.

Faraday rotation is a powerful technique for determining $B_{\rm{LOS}}$ and
refers to the rotation of the polarization plane of a linearly polarized
electromagnetic wave in a magneto-ionic region. The amount of rotation in the
medium is described by

$\Delta\Psi=\lambda^{2}\bigg{(}0.812\int n_{\text{e}}\mathbf{B\cdot
dl}\bigg{)}=\lambda^{2}\rm\text{RM}\rm~{}[rad],$ (1)

where $\Delta\Psi$ [rad], $\lambda$ [m], $\mathbf{B}$ [$\mu$G], $\mathbf{dl}$
[pc], and $n_{\text{e}}$ [cm-3] are the amount of rotation in polarization
angle, wavelength of the electromagnetic wave, external magnetic field, path
length through the magnetized region, and electron number density of the
region, respectively. The quantity in brackets is an integral of the product
of the electron density and the magnetic field along the path length, and it
is known as the rotation measure (RM).

The average $B_{\rm{LOS}}$ of a cloud or region can be determined by observing
its induced RM, where the RM’s sign indicates the $B_{\rm{LOS}}$ direction.
This is achieved by studying the observed polarization angle with respect to
$\lambda^{2}$ over a range of wavelengths (e.g., Brown & Taylor, 2001).
Assuming $B_{\rm{LOS}}$ is constant within the region and can be taken out of
the integral, only the electron column density of the region needs to be
determined to find $B_{\rm{LOS}}$:

$B_{\rm{LOS}}=\frac{RM}{\int n_{e}dl}=\frac{RM}{N_{e}},$ (2)

where $N_{e}$ is the electron column density of the region.

However, until relatively recently, despite previous attempts (Wolleben &
Reich, 2004), Faraday rotation had not been successfully applied to study
$B_{\rm{LOS}}$ in molecular clouds. Tahani et al. (2018) demonstrated the
feasibility of using Faraday rotation to determine $B_{\rm{LOS}}$ associated
with molecular clouds. Their technique (hereafter MC-BLOS method) uses the RM
value of point sources (such as unresolved pulsars, or radio galaxies) and an
ON-OFF approach to find the RM induced by the cloud and effectively separate
the RM contribution of the magnetized molecular cloud from that of the
foreground and background regions. This enables the determination of
$B_{\rm{LOS}}$ direction associated with the molecular cloud at each point. To
determine the electron densities and hence the strength of magnetic fields,
the technique uses a chemical evolution code and extinction ($A_{{V}}$ or
column density) maps to model the ionization state of the cloud.

The MC-BLOS technique has been applied to study and map the $B_{\rm{LOS}}$ of
several molecular clouds, and the results have been found to be consistent
with available Zeeman measurements (Tahani et al., 2018). These $B_{\rm{LOS}}$
maps have since been applied to study and reconstruct the 3D magnetic field
structure in the Orion A and Perseus molecular clouds (Tahani et al., 2022b,
a). These studies are the first to provide a complete picture of the 3D
magnetic field vectors in molecular clouds, including the $B_{\rm{LOS}}$
direction and the $B_{\rm{POS}}$ orientation and direction (i.e., complete
3D).

Using the 3D magnetic field structure of the Perseus molecular cloud, Tahani
et al. (2022b) proposed the cloud’s interaction with a previously-unidentified
structure, presence of which was later confirmed by Kounkel et al. (2022)
using kinematic observations. These studies demonstrate the power of the MC-
BLOS technique in revealing 3D magnetic fields and hence new insights into the
complex interplay between magnetic fields and cloud dynamics.

New and upcoming RM observations, such as those from POSSUM (Gaensler et al.,
2010; Vanderwoude et al., 2024), LOFAR (Van Eck et al., 2017; O’Sullivan et
al., 2023), and Square Kilometer Array (SKA; Heald et al., 2020), necessitate
a software and upgraded technique to determine the $B_{\rm{LOS}}$ associated
with molecular clouds in an automated manner.

In this paper, we introduce the python-based, open-source software, MC-BLOS
(v1.0), that performs the MC-BLOS technique automatically. Additionally, we
have made upgrades to the technique to include improvements for better
sampling of non-cloud contributions to the RM, such as an optimized algorithm
for selecting final OFF positions (referred to as reference points222we refer
to the the final selected OFF points as reference points. Reference and
“selected OFF” points may be used interchangeably.) and a more robust method
for estimating the background RM. Furthermore, the new software is capable of
performing the technique on multiple clouds in one run, enhancing its
efficiency and usability. These upgrades make MC-BLOS (v1.0) a valuable tool
for analyzing the upcoming wealth of RM data and studying the 3D magnetic
field structure of molecular clouds. We discuss the MC-BLOS technique, the
software (including the upgrades), the results, and a summary in sections 2,
3, 4, and 5, respectively.

##  2 General theory

The MC-BLOS technique consists of two main components: direction determination
and strength determination. The direction determination component isolates the
RM contribution of the molecular cloud from that of the foreground and
background regions (Galactic contribution), while the strength determination
component estimates the magnetic field strength using electron density
information derived from chemical evolution models and extinction maps. This
separation allows the MC-BLOS technique to effectively constrain both the
direction and strength of the line-of-sight magnetic field in molecular
clouds.

###  2.1 Determining the RM induced by the cloud

The RM induced by the cloud can be estimated by averaging the RM values of
reference points, which are located near the cloud but far enough to be
unaffected by its magnetic fields. OFF points are selected from regions with
low column densities (e.g., $A_{V}<1$ mag) that are not near high-column
density regions. Averaging OFF points from different areas around the cloud
ensures proper sampling of the Galactic contribution and that ON and OFF
points probe similar path lengths.

Tahani et al. (2018) define the average of reference positions as the
reference RM (RMref). They first identify all the potential OFF positions (N
points) by taking a threshold value for $A_{{V}}$ and finding the RM points
that have an associated $A_{{V}}$ of less than the threshold. To find the
optimal number of OFF positions (i.e., the number of points needed until the
$B_{\rm{LOS}}$ values stabilize), they introduce a “stability trend” analysis,
studying how derived magnetic field strength and direction change as the
number of OFF positions increases from 1 to N .

Initially, with few OFF positions, there is large variance in the results. As
the number increases, variations decrease, and magnetic field strengths and
directions stabilize. This approach ensures an accurate estimation of the
Galactic RM contribution using sufficient OFF positions, while retaining
enough RMs for $B_{\rm{LOS}}$ determination. In this work, we refer to the RM
of each individual OFF and ON point as RMOFF and RMON, respectively, while
RMref represents the final estimated Galactic contribution to the RM.

In the upgraded version (see Section 3), we improve upon the stability trend
by automating the selection of reference positions and by using a more robust
algorithm to identify the optimal number of OFF positions (see Section 3.2 for
details). Furthermore, we add additional steps to improve the selection of
reference points. Subtracting this RMref from the RMON at each point gives us
the RM contribution associated with the cloud’s magnetic field at each ON
position.

###  2.2 Determining the electron column density

After determining the RM contribution of the cloud at each position, the
electron column densities are needed to estimate the magnetic field strengths.
To this end a chemical evolution code and extinction (column density) maps are
required. Tahani et al. (2018) used a chemical evolution code previously
tested and utilized by Gibson et al. (2009).

The chemical evolution code used in the MC-BLOS technique takes into account
various input parameters, such as gas density, temperature, UV field strength,
and cosmic ray ionization rate. These parameters affect the electron abundance
calculations and need to be carefully constrained for each molecular cloud.
The code uses the UMIST Rate 99 database for reaction rates (Le Teuff et al.,
2000) and assumes a constant density, temperature, UV field, and cosmic ray
ionization rate for each cloud, which can be set as input parameters. The
cloud is taken to be homogeneous and planar, divided into 100 layers of equal
width. The electron column density for a given position cannot be obtained
from the electron abundance of a single layer multiplied by the hydrogen
column density, as the electron abundance may vary in each layer.

MC-BLOS first determines the extinction of each ON position and then
calculates the electron column density of all cloud layers up to the point
where half the extinction is reached, assuming cloud symmetry along the line
of sight, as shown in Figure 1. The total electron column density of the ON
position is obtained by summing the electron column densities of all layers
and multiplying by two, as follows:

$\displaystyle N_{e}=$ $\displaystyle 2\times\Sigma N_{e,i}$ (3)
$\displaystyle=$ $\displaystyle 2\times\Sigma^{\frac{A_{\text{V,
MC}}}{2}}\big{(}(A_{V,i}-A_{V,i-1})\times X_{e,i}\big{)}\times\text{$A_{{V}}$\
to N conversion},$

where $N_{e,i}$ and $X_{e,i}$ are the electron column density and abundance at
layer $i$, respectively. $A_{V,i}$ is the extinction at the same layer.

Finally, it determines $B_{\rm{LOS}}$ using the total electron column density
and RMref. The upgraded version of the software (Section 3.3) allows users to
input their own chemical evolution code and provides default input parameters
for many commonly studied molecular clouds, such as those observed by
Herschel333http://www.herschel.fr/cea/gouldbelt/en/Phocea/Vie_des_labos/Ast/ast_visu.php?id_ast=66.

Figure 1: Cloud division into layers for accurate electron column density
estimation. Assuming a symmetrical cloud, the electron column density is
calculated for each layer up to the center of the cloud along the line of
sight, reaching $\frac{A_{V,\text{MC}}}{2}$.

###  2.3 Determining Uncertainties

To estimate the total uncertainties, Tahani et al. (2018) consider the
following contributions: a) RMON uncertainty and standard deviation of RMref
($\frac{\delta RM}{RM}$), b) electron density estimation using the chemical
evolution code input values, particularly cloud-specific ones such as
temperature ($\frac{\delta T}{T}$) and density ($\frac{\delta
n({\text{HI}+\text{H}_{2}})}{n({\text{HI}+\text{H}_{2}})}$), and c)
uncertainties due to overlapping RM values on extinction (or column density)
maps ($\frac{\delta A_{V}}{A_{V}}$).

They estimate uncertainties in the chemical evolution code input values by
varying each parameter within a reasonable range (e.g., $\pm 20\%$ for
temperature, $\pm 50\%$ for density) and calculating the corresponding change
in magnetic field. Uncertainties due to overlapping RM values on extinction
maps are estimated by comparing extinction values of nearby pixels and
propagating uncertainty caused by extinction differences through the magnetic
field strength calculation. The total uncertainty is given by:

$\displaystyle\delta B_{\text{LOS}}=B_{\text{LOS}}\bigg{(}(\frac{\delta
RM}{RM})^{2}+(\frac{\delta A_{V}}{A_{V}})^{2}$ (4)
$\displaystyle+(\frac{\delta
n({\text{HI}+\text{H}_{2}})}{n({\text{HI}+\text{H}_{2}})})^{2}+(\frac{\delta
T}{T})^{2}\bigg{)}^{1/2}$ .

To find the $B_{\rm{LOS}}$ uncertainty from RM at each point, we use:

$\displaystyle\Delta
B_{RM}=B_{\text{LOS}}\bigg{(}\frac{\delta(RM_{\text{ref}})+\delta(RM_{ON})}{RM_{ON}-RM_{\text{ref}}}\bigg{)},$
(5)

where $\Delta B_{RM}$ is the uncertainty in $B_{\rm{LOS}}$ caused by RM (at
each point and associated with RMref), $\delta(RM_{\text{ref}})$ is the
standard error (standard deviation/$\sqrt{N}$) of the selected OFF positions,
and $\delta(RM_{ON})$ is the tabulated uncertainty of the RM at each ON point
(from the RM catalog).

The presented software follows the same uncertainty determination approach as
described in Tahani et al. (2018) with new upgrades discussed in Section 3.4.
These upgrades include: 1) ensuring that while the uncertainties caused by the
RM or the RMref are allowed to result in a change of direction, the
uncertainties due to the chemical code or extinction can only influence the
strength of $B_{\rm{LOS}}$ and not its sign; 2) ensuring that the ON points
have an $A_{{V}}$ higher than $3\times$AV,ref.

##  3 Software

###  3.1 Data files

To run the software, column density or extinction maps of molecular clouds and
an RM catalog are needed. Upon initializing the software (running
Initialize.py), cloud maps and the Taylor et al.
(2009)444https://cdsarc.cds.unistra.fr/ftp/J/ApJ/702/1230/ RM catalog are
downloaded, as well as the Van Eck et al. (2023) consolidated RM
catalog555CIRADA-tools; https://github.com/CIRADA-Tools/RMTable. Any RM
catalog with a format consistent with Van Eck et al. (2023) or Taylor et al.
(2009) can be used in the software. The cloud maps are mainly downloaded from
publicly-available Herschel observations such as the Herschel Gould Belt
Survey666http://www.herschel.fr/cea/gouldbelt/en/Phocea/Vie_des_labos/Ast/ast_visu.php?id_ast=66.
Users can input their own cloud map (either extinction or column density) or
rotation measure catalogs by saving them in the Data directory and specifying
the file names in the relevant config files.

###  3.2 $B_{\rm{LOS}}$ Direction determination

This section describes how the software determines RMref, which sets the
direction of the magnetic field at each ON point associated with the cloud.
The rotation measure induced by the cloud at any ON position is determined by:

$RM_{MC}=RM_{ON}-RM_{ref}.$ (6)

Therefore, if RMMC is positive (negative), the direction of the magnetic field
associated with the cloud at that position is pointing toward (away from) us,
following the convention used in the Faraday rotation measure community. After
finding RMMC for each ON point, we can also find the strength of the magnetic
field at that point (see Section 3.3).

To find RMref, the following steps are taken as illustrated in Figure
(Flowchart) 2:

  1. 1.

An OFF extinction threshold (AV,th) is set by default or by the user. AV,th
determines the maximum $A_{{V}}$ value that can be considered as a potential
OFF position – higher $A_{{V}}$ values would indicate that their corresponding
RM is influenced by the cloud’s magnetic field. The software determines the
appropriate AV,th based on the Galactic coordinates of each cloud, using the
center location estimated from high extinction regions in the corresponding
fits map. The default AV,th is as follows:

    1. (a)

For high-latitude clouds (absolute central latitude $>15^{\circ}$), the
default AV,th is 1. These clouds are relatively isolated from the complex
Galactic plane, allowing for a lower threshold.

    2. (b)

For anti-galactic-pointing (pointing away from the Galactic center) and low-
latitude clouds (absolute central latitude $<15^{\circ}$), the default AV,th
is 1.5. These clouds experience moderate Galactic interference, necessitating
a slightly higher threshold.

    3. (c)

For galactic-pointing (pointing towards the Galactic center) and low-latitude
clouds, the default AV,th is 2. These clouds are in regions of high Galactic
complexity, requiring the highest threshold to distinguish cloud effects from
the Galactic contribution.

  2. 2.

All the RM points are sorted based on their $A_{{V}}$ value, from lowest up to
AV,th. This set of RM points defines the potential OFF points.

  3. 3.

The software searches for anomalous OFF points in the potential OFF points
set. The anomalous points include any potential OFF position that is too close
to the cloud or has an RM value significantly different from the rest of the
OFF positions.

     * •

Too close to the cloud: If the user has set “use near high extinction
exclusion” = True (default) in the config file, then any point too close to
high extinction regions will be excluded from the potential OFF points. This
is because points too close to the cloud might be influenced by the cloud’s
magnetic field. The user can adjust the “too close” criterion by adjusting the
“near high extinction threshold multiplier” and “high extinction multiplier”
parameters in the config file. The “high extinction threshold multiplier” is
multiplied by AV,th to set a high $A_{{V}}$ value. Reference points should be
located far enough away from regions with extinction values above this
threshold. The “near high extinction multiplier” is multiplied by the
resolution or pixel size of the extinction or column density map to define the
proximity threshold. Points closer to high extinction regions than this
threshold are considered too close and will be excluded from the potential OFF
points. We discuss these choices in more details in Section 4.

     * •

Anomalous RM value: If the user has set “use anomalous value removal” = True
(default) in the config file, potential OFF points with RM values that are
significantly higher or lower compared to other OFF positions are excluded
from the final reference points. These outlier RM values are likely caused by
astrophysical phenomena other than the Galactic contribution to the RM, such
as intrinsic RM from the source itself. The user can adjust the definition of
an “anomalous RM” by modifying the “anomalous values iqr multiple (greater
than or equal to)” ($C_{IQR,RM}$) in the config file. The software calculates
the median ($RM_{med}$) and interquartile range ($RM_{IQR}$) of all potential
OFF points. The interquartile range is the difference between the third
quartile (Q3) and the first quartile (Q1). If a potential OFF point has an RM
value higher than $RM_{\text{Upper Limit}}=RM_{med}+C_{IQR,RM}\times RM_{IQR}$
or lower than $RM_{\text{Lower Limit}}=RM_{med}-C_{IQR,RM}\times RM_{IQR}$, it
is considered anomalous.

     * •

Too far from the cloud: This condition is similar to the “too-close-to-cloud”
condition. If the user has set “use far high extinction exclusion” = True in
the config file, then any point too far from the cloud or high extinction
regions will be excluded from the potential OFF points set. This is necessary
to ensure that we are sampling the Galactic contribution to the RM (for
subtraction from RMON), requiring a similar path and path length along the
line of sight. The software’s default setting is False, as the user already
identifies the desired region for the RMs (ON and OFF) to be considered in the
config file or through the cloud column density or extinction fits file. The
user can adjust the definition of “too far” by modifying the “far from high
extinction threshold multiplier”. The “far from high extinction multiplier”
parameter is multiplied by the resolution or the pixel size of the extinction
(or column density) map of the cloud to define the distance threshold. The
“high extinction threshold multiplier” is used to set an $A_{{V}}$ value, and
RM points should not be too far from this value. Figure 3 illustrates
anomalous potential OFF points identified using the default input parameters
for the Perseus molecular cloud.

  4. 4.

The stability trend is used to determine the optimal number of candidate OFF
points at which the trend of calculated $B_{\rm{LOS}}$ values stabilizes. This
is done by analyzing the differences in $B_{\rm{LOS}}$ values for different
numbers of OFF points and identifying the lowest number of OFF points needed
for a stable $B_{\rm{LOS}}$ behavior when more OFF points are added. The goal
is to find the lowest number of reference points that result in a stable trend
of $B_{\rm{LOS}}$ values, avoiding the use of too many OFF points in a limited
RM dataset and allowing determination of $B_{\rm{LOS}}$ associated with lower
density regions. Tahani et al. (2018) determined the optimal number of
reference points primarily through visual inspection of stability trend plots
(see their Figure 4), and this process is now automated. Figure 4 illustrates
the stability trend plot for the Perseus cloud. The software’s approach to the
stability trend is as follows:

     * $\blacksquare$

Initially, all non-anomalous potential OFF points are gathered. The software
calculates the $B_{\rm{LOS}}$ values using an increasing number of OFF points
ordered by extinction value from low to high. That is, the software starts
with a single OFF point with the lowest extinction value, then the two OFF
points with the lowest extinction values, and so on, up to including all N
potential OFF points. This process generates a set of $B_{\rm{LOS}}$ values
for each RM point, enabling the exploration of the effects of reference points
and the determination of the optimal number of OFF points.

     * $\blacksquare$

For each set of OFF points, the differences between adjacent $B_{\rm{LOS}}$
values (calculated using i and i+1 reference points) are computed for the
remaining RM points. The maximum and minimum $B_{\rm{LOS}}$ differences among
all the RM points are considered the stability maximum and minimum thresholds.
These thresholds enable the determination of when the $B_{\rm{LOS}}$ values
stabilize as OFF points are added. The software determines the number of OFF
points at which most RM points have adjacent $B_{\rm{LOS}}$ differences within
the stability thresholds, and this number is taken as the optimal number of
reference points.

  5. 5.

If the user has set the quadrant selection to True, the software divides the
cloud into four quadrants based on extinction values of the region777To divide
the cloud into quadrants, first a center point in the map is found weighted by
extinction distribution of the region. Then the software finds a line passing
through this center point and divides the entire region into two equally-
weighted regions. A perpendicular line passing through the first line and the
center point is then found, dividing the cloud into four quadrants.. The
software starts taking points one by one as sorted in the potential OFF points
set until it satisfies the minimum number of points per quadrant specified by
the user. If there are not enough potential OFF points in a quadrant, the
software continues taking OFF points from the sorted set until it takes all
the OFF points in that quadrant. Figure 3 illustrates the quadrants of the
Perseus cloud.

  6. 6.

After the final OFF (reference) points are selected, an average is calculated
to determine RMref. The software default is to simply average the RM value of
all the selected reference points to determine RMref.

$\displaystyle\text{RM}_{\text{ref}}=\sum_{i=1}^{N}\frac{\text{RM}_{\text{OFF,
i}}}{N}$ (7) $\displaystyle\text{A}_{\text{V,
ref}}=\sum_{i=1}^{N}\frac{\text{A}_{\text{{V~{}OFF, i}}}}{\text{N}}.$

However, if the user has set the quadrant weighting scheme to True, a
weighting scheme is applied to ensure that the effects of each cloud quadrant
are the same in finding RMref. This ensures that if, for example, different
quadrants have different numbers of OFF points, they both contribute equally
to RMref.

  7. 7.

After RMref is found, the contribution of the molecular cloud’s magnetic field
to each RMON can be estimated as follows:

$\displaystyle\text{RM}_{\text{MC}}=\text{RM}_{\text{ON}}-\text{RM}_{\text{ref}}=\bigg{(}0.812\int
n_{\text{e}}B_{\text{LOS}}dl\bigg{)}_{\text{MC}},$ (8)

and the extinction value of the cloud at each ON point is:

$A_{\text{V, MC}}=A_{\text{V, ON}}-A_{\text{V, ref}}.$ (9)

In summary, this process for determining RMref is designed to be both robust
and flexible. It accounts for various factors that could influence the
accuracy of the reference RM, including proximity to high-extinction regions,
anomalous RM values, and the spatial distribution of reference points across
the cloud. The stability trend analysis ensures that an optimal number of
reference points is used, balancing the need for statistical reliability with
the preservation of data for $B_{\rm{LOS}}$ determination. This comprehensive
approach to RMref determination forms a crucial foundation for accurate
$B_{\rm{LOS}}$ mapping in molecular clouds using the MC-BLOS technique.

RMref Determination Identify potential OFF points: sort all the RMs up to
$A_{\rm{v}}<A_{\rm{V,~{}th}}$ Too close to the cloud? Exclude the point and
continue with the remaining potential OFF points Anomalous RM points Exclude
the point and continue with the remaining potential OFF points Run stability
trend to determine optimal number of reference points Quadrant selection
Continue selecting potential OFF points until the specified minimum number of
OFF points per quadrant is satisfied Rerun the stability trend for user’s
study Weighting Scheme Average reference points for RMref determination
Determine the $RM_{ref}$ weighted equally by each quadrant
TRUEFALSETRUEFALSETRUEFALSEFALSETRUE

Figure 2: Flowchart showing the direction determination steps. The software
finds RMref by identifying potential OFF points, excluding anomalous ones,
finding the optimal number of OFF points, and finally averaging them. Figure
3: Perseus quadrants and anomalous (rejected) reference points (green circles)
identified using the default input parameters for the Perseus Molecular Cloud.
The background image shows the extinction map of the region. Points 63 and 76
are both too close to high extinction regions. No anomalous reference RM
values were found. The cloud is divided into four regions for a thorough
selection of reference points. If selected, the software ensures that each
quadrant contributes equally to the determination of RMref. Figure 4:
Stability trend plot for the Perseus cloud, showing the number of reference
points (x-axis) versus the calculated $B_{\rm{LOS}}$ values (y-axis) for each
RM point. The $B_{\rm{LOS}}$ of each RM point is calculated by taking an
increasing number of OFF points to find the optimal number.

###  3.3 Strength determination

After RMref and AV,ref are determined, the $B_{\rm{LOS}}$ of each remaining ON
point can be obtained using $A_{{V}}$ (or column density) maps and a chemical
evolution code. The software reads and analyzes chemical results from output
files generated by a previously tested chemical evolution code (Gibson et al.,
2009). These results cover various cloud parameters: density ($n$),
temperature ($T$), cosmic ionization rate ($IR_{CR}$), and ultra-violet
ionization rate ($IR_{UV}$). Users can also employ their own chemical
evolution code (see the user manual for details). The current chemical code
divides the cloud into layers and calculates the extinction and element
abundances (including electrons) for each layer.

Concurrently, the $A_{{V}}$ (or column density) of each point is read from the
maps. Assuming a symmetrical cloud, the electron column density is calculated
for each layer up to the cloud’s center along the line of sight, reaching
$\frac{A_{V,MC}}{2}$, as illustrated in Figure 1. The total electron column
density for each ON point is then obtained by doubling this value, as
explained in Section 2.2. Finally, the $B_{\rm{LOS}}$ strength is determined
using equations 2, 3, and 6.

$N_{e}$ Determination Retrieve chemical results using corresponding input
parameters ($IR_{CR}$, $IR_{UV}$, $T$, $n$) Find corresponding $A_{V}$ (or
$N$) at each $RM_{ON}$ Divide cloud into layers at each ON point to
$\frac{A_{V}}{2}$ and find $X_{e,i}$ in each layer Determine N$e,i$ based on
$X_{e,i}$ and $A_{v,i}$ Find $2\times\Sigma_{i}N_{e,i}$

Figure 5: Flowchart illustrating the strength determination steps. In this
process, a chemical evolution code along with extinction maps are used to
estimate the electron column densities at each layer of the cloud. The output
results of a chemical code are included in the software for different cloud
density ($n$), temperature ($T$), cosmic ionization rate ($IR_{CR}$), and
ultra-violet ionization rate ($IR_{UV}$).

###  3.4 Finding total uncertainty

As discussed in Section 2.3, the uncertainties of each point are obtained by
considering the uncertainties in RM, RMref, electron column densities, and the
overlap of RM and $A_{{V}}$. The electron column density uncertainties stem
from uncertainties in extinction and the chemical evolution code. We calculate
the $N_{e}$ uncertainty by varying the chemical code’s input temperature by
$\pm 20\%$ and input volume density by $\pm 50\%$, as well as finding the
maximum and minimum $A_{{V}}$ values within a Jeans length radius around the
RM value.

The uncertainties in RM and RMref can cause a change in the direction of the
magnetic field, as these parameters determine the direction of each
$B_{\rm{LOS}}$ point. Conversely, the $N_{e}$ uncertainties only affect the
strength of $B_{\rm{LOS}}$, not its direction. The direction is determined
solely by RM and RMref, as the sign of RMON determines the direction. The
upgraded technique now accounts for this distinction; if the $N_{e}$
uncertainties exceed the $B_{\rm{LOS}}$ value, the software adjusts these
uncertainties to ensure they do not affect the sign of $B_{\rm{LOS}}$.

Additionally, the software ensures all ON points have an $A_{{V}}$ value that
is a user-specified multiple higher than AV,ref (i.e.,
$A_{{\rm{V,ON}}}\geq\text{multiple}\times A_{\rm{V,ref}}$). This criterion
enhances the fidelity of magnetic field strength determination, as the
$A_{{V}}$ value is crucial for strength determination. The multiplication
factor for setting the minimum ON point extinction can be adjusted in the
input parameter file by changing the “on point extinction multiple of off
point average multiplier”. For example, if set to three, the software ensures
all ON points have $A_{V}\geq 3\times A_{\rm{V,ref}}$. While ideally this
factor should be greater than one, the current default is set to one due to
the low source density in available RM catalogs. Future RM catalogs with
higher source densities will enable better selection of ON points.

Regardless of the chosen multiplication factor, the software always ensures
that any ON point has an $A_{{V}}$ value higher than each of the reference
points. This criterion does not influence the $B_{\rm{LOS}}$ direction. Even
with a multiplication factor of one, we obtain $B_{\rm{LOS}}$ directions that
match previously determined directions from Zeeman observations (Heiles, 1987,
1997).

While the $B_{\rm{LOS}}$ strengths of points with $A_{{V}}$ not significantly
higher than AV,ref may exceed those from Zeeman observations, they show good
correspondence with Zeeman-obtained strengths. For example, Tahani et al.
(2018) found that the error-weighted average $B_{\rm{LOS}}$ on the western
side of Orion A is half that of its eastern side, consistent with Zeeman
observations (Tahani et al., 2022a).

##  4 Results and Discussion

In this section, we present the MC-BLOS final results for the Orion A,
Perseus, California, and Orion B molecular clouds and compare them with the
previously-published results by Tahani et al. (2018) and nearby molecular
Zeeman effect observations. Although users can influence the selection of
reference points, the software-suggested reference points produce results
consistent with Tahani et al. (2018), with differences in the final
$B_{\rm{LOS}}$ values falling within the uncertainty ranges for each cloud.

To compare our results with those of Tahani et al. (2018), we use the RM
catalog of Taylor et al. (2009). As the same RM catalogs are used in both
studies, the main differences in final $B_{\rm{LOS}}$ values should be due to
differences in RMref values. We examine the software-selected RMref using
three input parameter sets:

  1. 1.

Default input parameters for each cloud (“Default” parameter set).

  2. 2.

Default input with at least one reference point in each cloud quadrant
(“Min1PerQ” parameter set).

  3. 3.

Default input with the “high extinction threshold multiplier” increased from 2
to 5, positioning OFF points further from high extinction regions (“NearExt5”
parameter set).

Table 1 compares the RMref values generated using these parameter sets with
those of Tahani et al. (2018). The software-selected reference points in all
three parameter sets result in RMref values similar to those of Tahani et al.
(2018), falling within their uncertainty ranges. For Orion A and California,
while the RMref values deviate from the original Tahani et al. (2018), they
remain consistent within the uncertainties.

Although the RMref values of Orion A slightly deviate from the original ones,
the direction of the reversal seen in the Orion A cloud maintains the same
orientation as Tahani et al. (2018) and Zeeman-obtained directions (Heiles,
1997). In both Orion A and Orion B clouds, we find points that closely align
with previously published molecular Zeeman observations.

In the Orion A region near $\alpha(J2000)\simeq
83.81^{\circ},\delta(J2000)\simeq-5.37^{\circ}$, different molecular Zeeman
values have been reported, such as $+360\pm 80,\mu$G (Falgarone et al., 2008;
Crutcher, 1999; Crutcher et al., 2010), $-79\pm 99,\mu$G (Crutcher et al.,
1996), and $-40\pm 240,\mu$G (Crutcher et al., 1999b, 2010). Considering the
uncertainties, these may indicate magnetic field strengths ranging from
$+440,\mu$G to $-280,\mu$G. Near this point, similar to the Zeeman
measurements, we find high uncertainties and magnetic field values of $-45\pm
36,\mu$G and $-6\pm 35,\mu$G at
$\alpha(J2000)=83.90^{\circ},\delta(J2000)=-5.39^{\circ}$ and
$\alpha(J2000)=83.90^{\circ},\delta(J2000)=-5.38^{\circ}$, respectively.

In Orion B, the observations are near $\alpha(J2000)\simeq
85.44^{\circ},\delta(J2000)\simeq 1.93^{\circ}$ with the values of $-270\pm
330,\mu$G (Crutcher et al., 1999b) and $-87\pm 5.5,\mu$G (Crutcher et al.,
1999a). Near this location, at
$\alpha(J2000)=85.42^{\circ},\delta(J2000)=1.92^{\circ}$ and
$\alpha(J2000)=85.42^{\circ},\delta(J2000)=1.90^{\circ}$, we find values of
$-110_{-24}^{+24}\,\mu$G and $-119_{-27}^{+28}\,\mu$G, respectively, with the
same direction.

We note that in high extinction regions, our $B_{\rm{LOS}}$ values may be
smaller than expected. This is because our method considers all cloud layers
along the line of sight, including potential depolarization effects.
Depolarization can occur due to variations in magnetic field direction across
different cloud layers. Changes in the sign of the magnetic field along the
line of sight can lead to a reduction in the net observed $B_{\rm{LOS}}$
value. This integrated approach may result in lower apparent field strengths
in regions where the magnetic field structure is more complex or tangled.

Cloud Model | Default | Min1PerQ | NearExt5 | Tahani+18
---|---|---|---|---
Orion A | $15.2\pm 12.6$ | $15.2\pm 12.6$ | $8.7\pm 13.0$ | $1.4\pm 13.7$
Perseus | $34.1\pm 10.6$ | $34.1\pm 10.6$ | $38.7\pm 12.7$ | $31.1\pm 11$
California | $13.4\pm 12.9$ | $13.4\pm 12.9$ | $13.4\pm 12.9$ | $4.0\pm 14$
Orion B | $26.6\pm 9.6$ | $31.3\pm 11.7$ | $26.6\pm 9.6$ | $32.3\pm 10$
Table 1: RMref values for each cloud and model using the upgraded technique
and RMref from Tahani et al. (2018). The differences between the models for
each cloud are within RM uncertainties.

Figure 6: Final $B_{\rm{LOS}}$ maps of the Orion A, Orion B, California, and
Perseus molecular clouds using Default parameter set. The background color
image shows extinction, and the circles represent $B_{\rm{LOS}}$ strengths
(without uncertainties). Red (blue) circles indicate magnetic fields pointing
away from (toward) the observer. Figure 7: Final $B_{\rm{LOS}}$ maps of the
Perseus molecular clouds using NearExt5 parameter set, positioning the
reference point (green circles) further from the cloud. The background color
image shows extinction, and the circles represent $B_{\rm{LOS}}$ strengths
(without uncertainties). Red (blue) circles indicate magnetic fields pointing
away from (toward) the observer.

Figure 6 illustrates the final $B_{\rm{LOS}}$ maps of the four clouds. The
background image shows extinction, and the circles depict $B_{\rm{LOS}}$
direction and strength (without uncertainty values). Red, blue, and green
circles illustrate magnetic fields pointing away from us, toward us, and the
reference points, respectively.

As mentioned earlier, the software’s default multiplication factor for setting
the minimum accepted AV,ON value is currently set to one. This can lead to ON
and reference points being positioned near each other. Despite this, we still
obtain acceptable $B_{\rm{LOS}}$ results because having a low $A_{{V}}$ value
is the most important criterion for setting OFF positions due to the
clumpiness of the clouds. The abundance of RM sources in upcoming catalogs
will result in better separation of ON and OFF points. As shown in Figure 7,
this can also be addressed by using the NearExt5 parameter set, which
positions OFF points further from the cloud. We note that both the NearExt5
and Default parameter sets produce very similar results.

Figure 6 exhibits the same $B_{\rm{LOS}}$ morphology as Tahani et al. (2018).
The Orion A, Perseus, and California molecular clouds show a reversal of
$B_{\rm{LOS}}$ from one side of the cloud to the other (along the cloud’s
minor axis). The RM points associated with Orion B are insufficient for a
field reversal conclusion. This reversal in Orion A with the same orientation
was previously observed using atomic Zeeman observations (Heiles, 1997; Heiles
& Robishaw, 2009).

Tahani et al. (2022b) reconstructed the complete 3D magnetic field of the
Perseus cloud (ignoring sub-parsec fluctuations along the field lines). The
reconstructed 3D magnetic field of the Perseus cloud by Tahani et al. (2022b)
suggested the presence of an interacting structure near the cloud, which was
responsible for the observed $B_{\rm{LOS}}$ reversal. The presence of proposed
structure was later supported by kinematic observations presented in Kounkel
et al. (2022), providing independent confirmation of the proposed interaction.
The agreement between the implications of the reconstructed 3D magnetic field
structure and subsequent kinematic observations strengthens the reliability
and validity of both the 3D fields (obtained using the MC-BLOS results) and
the MC-BLOS method itself.

We note that the uncertainties in RMref in our study and Tahani et al. (2018)
are expected to be significantly reduced by upcoming observations, such as
those with the Square Kilometer Array Observatory, which will provide more
accurate and abundant RM data. This should improve the overall accuracy and
precision of the MC-BLOS results.

##  5 Summary and Plans for maintaining MC-BLOS

###  5.1 Summary

We have upgraded the Faraday-based technique of Tahani et al. (2018) and
present MC-BLOS (v1.0), a Python software package for determining the line-of-
sight component of magnetic fields associated with molecular clouds. The
software automates the process of determining reference points, which
previously required manual analysis, and incorporates upgrades in the
selection and handling of OFF points and $B_{\rm{LOS}}$ uncertainty values.
MC-BLOS utilizes available rotation measure catalogs (point sources),
extinction (or column density) maps, and output results from chemical
evolution modeling to produce $B_{\rm{LOS}}$ maps and tables associated with
molecular clouds.

The software employs an ON-OFF approach to find the RM contribution of the
cloud. By using a number of OFF points that do not probe the magnetic field of
the clouds, the contribution of everything along the line of sight except the
cloud and its magnetic field can be determined. The sign of RMs at each point
associated with the cloud determines the $B_{\rm{LOS}}$ direction. To
determine the strength of magnetic fields, the software uses the RM induced by
the cloud at each point (determined using RMref and the RM value at each
point) and electron column densities, which are obtained for each point using
the extinction (column density) maps and the chemical evolution code.

The results presented in this paper incorporated the previously available RM
catalog of Taylor et al. (2009). This catalog used the VLA Sky Survey data
(NVSS; Condon et al., 1998), which had only two observing frequencies,
resulting in high uncertainties when inferring RM values. Upcoming RM
observations by POSSUM, VLASS, and SKA are expected to significantly improve
the RM uncertainties and increase the RM density by a factor of 10 (e.g.,
Vanderwoude et al., 2024), resulting in improved RMref uncertainties. Using
the MC-BLOS software along with these future observations, the line-of-sight
magnetic field component of many molecular clouds can be determined with
greater accuracy and precision.

###  5.2 Plans for maintaining and upgrading MC-BLOS software

As this paper is the first official publication of MC-BLOS software, we are
committed to its continuous improvement and upgrade. Our development roadmap
is guided by new and improved observations and feedback from the scientific
community. For instance, recent VLA observations (observation IDs VLA/19B-053
and VLA/24A-376) not only enable detailed $B_{\rm{LOS}}$ and 3D magnetic field
maps of the Perseus cloud (Hajizadeh et al., in prep.) but also provide an
opportunity to refine the RMref determination in MC-BLOS, thanks to their
increased source density compared to the NVSS catalog.

we envision several key areas for improvement:

  * •

Improved RMref determination: With higher source densities and more OFF
points, we plan to replace the current stability trend method with more
advanced techniques. This includes:

    * –

Implementing sophisticated statistical methods for robust outlier detection
and handling.

    * –

Exploring the potential of machine learning algorithms to recognize patterns
in RM distributions and optimize reference point selection.

    * –

Utilizing 3D dust maps.

These enhancements aim to improve the accuracy and reliability of RMref
determinations, which are crucial for precise $B_{\rm{LOS}}$ mapping.

  * •

Algorithmic Improvements: We plan to refine our algorithms as more RM data
become available, particularly in areas such as anomaly detection, cloud
boundary determination, and the handling of complex cloud morphologies.

To facilitate these advancements, we encourage feedback, bug reports,
suggestions, and contributions from the scientific community. We plan to
maintain and update the software regularly, releasing new versions as
significant improvements are made.

## acknowledgments

M.Tahani is supported by the Banting Fellowship (Natural Sciences and
Engineering Research Council Canada) hosted at Stanford University and the
Kavli Institute for Particle Astrophysics and Cosmology (KIPAC) Fellowship.
Claude.AI and QuillBot888https://quillbot.com/ were used for proofreading and
language refinement. We have used PyCharm.

## References

  * Arzoumanian et al. (2021) Arzoumanian, D., Furuya, R. S., Hasegawa, T., et al. 2021, A&A, 647, A78, doi: 10.1051/0004-6361/202038624
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Brown & Taylor (2001) Brown, J. C., & Taylor, A. R. 2001, ApJ, 563, L31, doi: 10.1086/338358
  * Ching et al. (2022) Ching, T.-C., Qiu, K., Li, D., et al. 2022, ApJ, 941, 122, doi: 10.3847/1538-4357/ac9dfb
  * Condon et al. (1998) Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693, doi: 10.1086/300337
  * Crutcher (1999) Crutcher, R. M. 1999, ApJ, 520, 706, doi: 10.1086/307483
  * Crutcher (2012) —. 2012, ARA&A, 50, 29, doi: 10.1146/annurev-astro-081811-125514
  * Crutcher & Kemball (2019) Crutcher, R. M., & Kemball, A. J. 2019, Frontiers in Astronomy and Space Sciences, 6, 66, doi: 10.3389/fspas.2019.00066
  * Crutcher et al. (1999a) Crutcher, R. M., Roberts, D. A., Troland, T. H., & Goss, W. M. 1999a, ApJ, 515, 275, doi: 10.1086/307026
  * Crutcher et al. (1996) Crutcher, R. M., Troland, T. H., Lazareff, B., & Kazes, I. 1996, ApJ, 456, 217, doi: 10.1086/176642
  * Crutcher et al. (1999b) Crutcher, R. M., Troland, T. H., Lazareff, B., Paubert, G., & Kazès, I. 1999b, ApJ, 514, L121, doi: 10.1086/311952
  * Crutcher et al. (2010) Crutcher, R. M., Wandelt, B., Heiles, C., Falgarone, E., & Troland, T. H. 2010, ApJ, 725, 466, doi: 10.1088/0004-637X/725/1/466
  * Doi et al. (2021) Doi, Y., Hasegawa, T., Bastien, P., et al. 2021, ApJ, 914, 122, doi: 10.3847/1538-4357/abfcc5
  * Doi et al. (2024) Doi, Y., Nakamura, K., Kawabata, K. S., et al. 2024, ApJ, 961, 13, doi: 10.3847/1538-4357/ad0fe2
  * Falgarone et al. (2008) Falgarone, E., Troland, T. H., Crutcher, R. M., & Paubert, G. 2008, A&A, 487, 247, doi: 10.1051/0004-6361:200809577
  * Gaensler et al. (2010) Gaensler, B. M., Landecker, T. L., Taylor, A. R., & POSSUM Collaboration. 2010, in American Astronomical Society Meeting Abstracts, Vol. 215, American Astronomical Society Meeting Abstracts #215, 470.13
  * Gibson et al. (2009) Gibson, D., Plume, R., Bergin, E., Ragan, S., & Evans, N. 2009, ApJ, 705, 123, doi: 10.1088/0004-637X/705/1/123
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Heald et al. (2020) Heald, G., Mao, S., Vacca, V., et al. 2020, Galaxies, 8, 53, doi: 10.3390/galaxies8030053
  * Heiles (1987) Heiles, C. 1987, in Astrophysics and Space Science Library, Vol. 134, Interstellar Processes, ed. D. J. Hollenbach & H. A. Thronson, Jr., 171–194, doi: 10.1007/978-94-009-3861-8_8
  * Heiles (1997) Heiles, C. 1997, ApJS, 111, 245, doi: 10.1086/313010
  * Heiles & Robishaw (2009) Heiles, C., & Robishaw, T. 2009, in IAU Symposium, Vol. 259, Cosmic Magnetic Fields: From Planets, to Stars and Galaxies, ed. K. G. Strassmeier, A. G. Kosovichev, & J. E. Beckman, 579–590, doi: 10.1017/S174392130903141X
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hwang et al. (2022) Hwang, J., Kim, J., Pattle, K., et al. 2022, ApJ, 941, 51, doi: 10.3847/1538-4357/ac99e0
  * Kounkel et al. (2022) Kounkel, M., Deng, T., & Stassun, K. G. 2022, arXiv e-prints, arXiv:2206.04703.  https://arxiv.org/abs/2206.04703
  * Krumholz & Federrath (2019) Krumholz, M. R., & Federrath, C. 2019, Frontiers in Astronomy and Space Sciences, 6, 7, doi: 10.3389/fspas.2019.00007
  * Lazarian & Hoang (2007) Lazarian, A., & Hoang, T. 2007, MNRAS, 378, 910, doi: 10.1111/j.1365-2966.2007.11817.x
  * Le Teuff et al. (2000) Le Teuff, Y. H., Millar, T. J., & Markwick, A. J. 2000, A&AS, 146, 157, doi: 10.1051/aas:2000265
  * O’Sullivan et al. (2023) O’Sullivan, S. P., Shimwell, T. W., Hardcastle, M. J., et al. 2023, MNRAS, 519, 5723, doi: 10.1093/mnras/stac3820
  * pandas development team (2020) pandas development team, T. 2020, pandas-dev/pandas: Pandas, latest, Zenodo, doi: 10.5281/zenodo.3509134
  * Panopoulou et al. (2023) Panopoulou, G. V., Markopoulioti, L., Bouzelou, F., et al. 2023, arXiv e-prints, arXiv:2307.05752, doi: 10.48550/arXiv.2307.05752
  * Pattle & Fissel (2019) Pattle, K., & Fissel, L. 2019, Frontiers in Astronomy and Space Sciences, 6, 15, doi: 10.3389/fspas.2019.00015
  * Pattle et al. (2023) Pattle, K., Fissel, L., Tahani, M., Liu, T., & Ntormousi, E. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 193, doi: 10.48550/arXiv.2203.11179
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of machine learning research, 12, 2825
  * Planck Collaboration et al. (2015) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2015, A&A, 576, A104, doi: 10.1051/0004-6361/201424082
  * Planck Collaboration et al. (2016) —. 2016, A&A, 586, A138, doi: 10.1051/0004-6361/201525896
  * Pudritz & Ray (2019) Pudritz, R. E., & Ray, T. P. 2019, Frontiers in Astronomy and Space Sciences, 6, 54, doi: 10.3389/fspas.2019.00054
  * Robishaw (2008) Robishaw, T. 2008, PhD thesis, University of California, Berkeley
  * Seifried & Walch (2015) Seifried, D., & Walch, S. 2015, MNRAS, 452, 2410, doi: 10.1093/mnras/stv1458
  * Seifried et al. (2019) Seifried, D., Walch, S., Reissl, S., & Ibáñez-Mejía, J. C. 2019, MNRAS, 482, 2697, doi: 10.1093/mnras/sty2831
  * Tahani (2022) Tahani, M. 2022, Frontiers in Astronomy and Space Sciences, 9, 940027, doi: 10.3389/fspas.2022.940027
  * Tahani et al. (2018) Tahani, M., Plume, R., Brown, J. C., & Kainulainen, J. 2018, A&A, 614, A100, doi: 10.1051/0004-6361/201732219
  * Tahani et al. (2019) Tahani, M., Plume, R., Brown, J. C., Soler, J. D., & Kainulainen, J. 2019, A&A, 632, A68, doi: 10.1051/0004-6361/201936280
  * Tahani et al. (2022a) Tahani, M., Glover, J., Lupypciw, W., et al. 2022a, A&A, 660, L7, doi: 10.1051/0004-6361/202243322
  * Tahani et al. (2022b) Tahani, M., Lupypciw, W., Glover, J., et al. 2022b, A&A, 660, A97, doi: 10.1051/0004-6361/202141170
  * Tahani et al. (2023) Tahani, M., Bastien, P., Furuya, R. S., et al. 2023, ApJ, 944, 139, doi: 10.3847/1538-4357/acac81
  * Taylor et al. (2009) Taylor, A. R., Stil, J. M., & Sunstrum, C. 2009, ApJ, 702, 1230, doi: 10.1088/0004-637X/702/2/1230
  * Van Eck et al. (2017) Van Eck, C. L., Haverkorn, M., Alves, M. I. R., et al. 2017, A&A, 597, A98, doi: 10.1051/0004-6361/201629707
  * Van Eck et al. (2023) Van Eck, C. L., Gaensler, B. M., Hutschenreuter, S., et al. 2023, ApJS, 267, 28, doi: 10.3847/1538-4365/acda24
  * Van Rossum (2020) Van Rossum, G. 2020, The Python Library Reference, release 3.8.2 (Python Software Foundation)
  * Vanderwoude et al. (2024) Vanderwoude, S., West, J. L., Gaensler, B. M., et al. 2024, AJ, 167, 226, doi: 10.3847/1538-3881/ad2fc8
  * Wolleben & Reich (2004) Wolleben, M., & Reich, W. 2004, A&A, 427, 537, doi: 10.1051/0004-6361:20040561