Fitting an Analytic Magnetic Field to a Prestellar Core

The method to model the magnetic field using polarimetry data involves an approach that is slightly different from that done by Ewertowski & Basu (2013), who fit the model magnetic field components directly to simulation data. When dealing with polarimetry data we do not know the magnitude of each individual component of the magnetic field, but rather can estimate the ratio B_r /B_z . We begin with the expressions for the magnetic field components given in Ewertowski & Basu (2013):

$$\begin{eqnarray}\begin{array}{rcl}{B}_{r}(r,z) & = & \displaystyle \sum _{m=1}^{\infty }{k}_{m}\sqrt{{\lambda }_{m}}{J}_{1}(\sqrt{{\lambda }_{m}}r)\left[\mathrm{erfc}\left(\displaystyle \frac{\sqrt{{\lambda }_{m}}h}{2}-\displaystyle \frac{z}{h}\right)\right.\\ & & \times \left.{e}^{-z\sqrt{{\lambda }_{m}}}-\mathrm{erfc}\left(\displaystyle \frac{\sqrt{{\lambda }_{m}}h}{2}+\displaystyle \frac{z}{h}\right){e}^{z\sqrt{{\lambda }_{m}}}\right],\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}{B}_{z}(r,z) & = & \displaystyle \sum _{m=1}^{\infty }{k}_{m}\sqrt{{\lambda }_{m}}{J}_{0}(\sqrt{{\lambda }_{m}}r)\left[\mathrm{erfc}\left(\displaystyle \frac{\sqrt{{\lambda }_{m}}h}{2}+\displaystyle \frac{z}{h}\right)\right.\\ & & \left.\times {e}^{z\sqrt{{\lambda }_{m}}}+\mathrm{erfc}\left(\displaystyle \frac{\sqrt{{\lambda }_{m}}h}{2}-\displaystyle \frac{z}{h}\right){e}^{-z\sqrt{{\lambda }_{m}}}\right]+{B}_{0},\end{array}\end{eqnarray} \tag{ 2 }$$

where h is the scale height for the Gaussian z-distribution of electric current density, and $\lambda ={\left({a}_{{\rm{m}},1}/R\right)}^{2}$, where a_m,1 is the mth root of the Bessel function of the first kind. In this situation we can observe (measure) the relative angle that each polarization segment makes with the normal (the z-axis in our case) and use this offset angle to fit the ratio B_r /B_z , that is we minimize the quantity

$$\begin{eqnarray}&&{\rm{\Delta }}({r}_{i},{z}_{i})=\displaystyle \sum _{i}{\parallel \displaystyle \frac{{B}_{r}({r}_{i},{z}_{i})}{{B}_{z}({r}_{i},{z}_{i})}-\tan {\phi }_{i}\parallel }^{2},\end{eqnarray} \tag{ 3 }$$

where the angles ϕ are data quantities and the ratio B_r /B_z is taken from the model. We can express the magnetic field nondimensionally by defining the parameters

$$\begin{eqnarray*}&&\tilde{r}=r/R,\ \ \ \ \ \ \ \tilde{z}=z/R,\ \ \ \ \ \ \ \eta =h/R,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\tilde{B}}_{r}={B}_{r}/{B}_{0},\ \ \ \ \ \ \ {\tilde{B}}_{z}={B}_{z}/{B}_{0},\ \ \ \ \ \ \ {\beta }_{m}={k}_{m}\sqrt{{\lambda }_{m}}/{B}_{0},\end{eqnarray*}$$

where we normalize the magnetic field components to the background field and normalize the positional quantities r and z to the core radius R. These give rise to the nondimensionalized counterparts to Equations (1) and (2):

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{B}}_{r} & = & \displaystyle \sum _{m=1}^{\infty }{\beta }_{m}{J}_{1}({a}_{{\rm{m}},1}\tilde{r})\left[\mathrm{erfc}\left(\displaystyle \frac{{a}_{{\rm{m}},1}\eta }{2}-\displaystyle \frac{\tilde{z}}{\eta }\right){e}^{-{a}_{{\rm{m}},1}\tilde{z}}\right.\\ & & \left.-\mathrm{erfc}\left(\displaystyle \frac{{a}_{{\rm{m}},1}\eta }{2}+\displaystyle \frac{\tilde{z}}{\eta }\right){e}^{{a}_{{\rm{m}},1}\tilde{z}}\right],\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{B}}_{z} & = & \displaystyle \sum _{m=1}^{\infty }{\beta }_{m}{J}_{0}({a}_{{\rm{m}},1}\tilde{r})\left[\mathrm{erfc}\left(\displaystyle \frac{{a}_{{\rm{m}},1}\eta }{2}+\displaystyle \frac{\tilde{z}}{\eta }\right){e}^{{a}_{{\rm{m}},1}\tilde{z}}\right.\\ & & \left.+\mathrm{erfc}\left(\displaystyle \frac{{a}_{{\rm{m}},1}\eta }{2}-\displaystyle \frac{\tilde{z}}{\eta }\right){e}^{-{a}_{{\rm{m}},1}\tilde{z}}\right]+1.\end{array}\end{eqnarray} \tag{ 5 }$$

This formalism gives freedom to the model user to select a value of B₀ and have each β_m normalized to that value. To perform the optimization we use the Sequential Quadratic Programming algorithm, part of Matlab's fmincon function, to minimize Equation (3).

We can use our models to produce synthetic polarization maps that simulate the magnetic field as a three-dimensional integral through the entire core rather than taking a slice of a single plane as done in Figure 2. To perform the simulations we used the three-dimensional radiative transfer code POLARIS developed by Reissl et al. (2016). We include grain alignment in the simulated polarization maps and restrict the alignment to be perfectly perpendicular to the magnetic field. This allows one to directly infer the magnetic field morphology through the orientation of dust grains. We model the gas profile with a modified isothermal sphere:

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{\rho }_{c}}{1+{\left(r/a\right)}^{2}},\end{eqnarray} \tag{ 6 }$$

with central density ρ_c = m_n n_c , central number density n_c = 3.5 × 10⁵ cm⁻³ (Kandori et al. 2005), and neutral mass m_n = 2.3 m_H. The radial scale length a is defined by ${a}^{2}=A\,{c}_{s}^{2}/(2\pi G{\rho }_{c})$, where c_s is the isothermal sound speed. We adopt A = 2 so as to achieve a Bonnor–Ebert-like density profile with a gas mass M_gas ≈ 3.6 M_⊙ within a radius R_gas ≈ 0.1 pc, in agreement with the fit of the gas density profile for FeSt 1–457 obtained by Kandori et al. (2005). For a detailed discussion of the relation between a modified isothermal sphere and a Bonnor–Ebert sphere, see Dapp & Basu (2009).

The grain alignment by radiative torques requires an anisotropic radiation source (see Andersson et al. 2015). Therefore, we implement a weakly radiating protostar of temperature T = 300 K and radius R = 4R_⊙ as the radiation source for our simulations, mimicking a possible observationally undetected first hydrostatic core (see, e.g., Enoch et al. 2010). The grain size distribution is taken to follow an MRN (Mathis et al. 1977) power-law distribution in which the number density of grains n_g satisfies ${{dn}}_{{\rm{g}}}/{da}\propto {n}_{{\rm{g}}}^{-3.5}$ from a minimum grain size ${a}_{\min }=0.005\,\mu {\rm{m}}$ to a maximum grain size ${a}_{\max }=2\,\mu {\rm{m}}$. We adopt a dust-to-gas mass ratio of 0.01. These and other properties are listed in Table 1, where we also adopt the gas and dust temperatures from Kandori et al. (2017).

We use the data set published by Kandori et al. (2017). Figure 1 shows the polarimetry data for the core FeSt 1–457, extracted from the data image using the Plot digitizer. We observe that the data is not symmetric about x = 0 and should expect this asymmetry to hinder our model fit in some sense. It is also not exactly symmetric about z = 0 but this is a lesser discrepancy and we do not investigate it further. Because of the differences in the number and morphology of the polarization segments on the different sides of x = 0, we perform three different fits of the magnetic field. We fit the entire data set, the left half of the data set and the right half of the data set. The left and right half fits are performed in an effort to see if we can improve the fit by removing the asymmetry in the data, but conclude that both models demonstrate poor fits to modeling the core as a whole. We perform a fit to the entire data set, adopting R = 0.1 pc, yielding the nondimensionalized parameters

$$\begin{eqnarray}&&{\boldsymbol{\beta }}=\left[\begin{array}{c}15.358\\ 0.502\\ 23.464\end{array}\right],\ \ \ \ \ \ \ \eta =0.0366.\end{eqnarray} \tag{ 7 }$$

We also perform a fit by fixing the core radius to R = 0.2 pc, corresponding physically to a condensation of the core from a significantly larger volume than implied by its density structure. In this case the optimal dimensionless parameters turn out to be

$$\begin{eqnarray}&&{\boldsymbol{\beta }}=\left[\begin{array}{c}0.611\\ 0.000\\ 4.480\end{array}\right],\ \ \ \ \ \ \ \eta =7.360\times {10}^{-4}.\end{eqnarray} \tag{ 8 }$$

The resulting field lines for both models are shown in Figure 2. The models' boundary conditions require that B_r will approach zero for r = R, i.e., the field becomes vertical at this radial distance. So, physically the model with R = 0.2 pc allows the curved magnetic field lines to span a larger region than is associated with the observed column density (R_gas ≈ 0.1 pc), while the R = 0.1 pc model fits the curved magnetic field lines in the same region as R_gas. Figures 3 and 4 show the values of B_r and B_z for both models along various cuts along the r- and z-directions. It can be seen that B_r approaches zero at both r = R and at large heights above the midplane. The strong degree of curvature in the field lines for offsets x ≳ 0.4 pc on the right side plane (see Figure 1) suggests a weak relative strength of the background magnetic field B₀ compared to the core-generated local field. In fact, the dominance of the self-induced magnetic field creates a region of closed magnetic field lines with negative B_z near the midplane and in the region 0.05 ≲ r ≲ 0.1 pc in both model fits. The statistical plots assessing the model's goodness of fit are given in Figure 5. Figure 6 demonstrates the magnetic field lines for both model fits overlaid with normalized magnetic field strength ($\sqrt{{B}_{r}^{2}+{B}_{z}^{2}}/{B}_{0}$) contours. These plots indicate that the model with smaller radius R = 0.1 pc requires a greater relative central concentration of magnetic field strength in order to fit the curvature of the observed field lines. Figure 7 demonstrates the decomposition of the magnetic field waveform into the three fitted modes. We truncate the series at m = 3 as done by Ewertowski & Basu (2013). We did make fits with higher order terms, and observed a similar lack of improvement to the fit and, in some cases, a reduction in the model performance. The plots show that the m = 1 mode (half-wavelength) needs to have a significant amplitude in order to fit the overall decline of B_z from center to edge. The high curvature of field lines within the core region emphasizes the m = 3 mode with little need for much contribution from the m = 2 mode.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To assess the goodness of fit, we need to make the standard assumption that the model residuals are normally distributed. In order to test this hypothesis, we perform a χ² goodness-of-fit test and a residual whiteness test. Pearson's χ²-test defines the statistic

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i}\displaystyle \frac{{\left({O}_{i}-{E}_{i}\right)}^{2}}{{E}_{i}},\end{eqnarray} \tag{ 9 }$$

where O_i and E_i are representative of the observed and expected counts, respectively. After removing the outliers ⁴ (we utilize Matlab's rmoutliers function), we achieve a χ² statistic of 6.51. This means that the test does not reject the null hypothesis that the model residuals are normally distributed. A histogram of our residuals are demonstrated in the left panel of Figure 5. We additionally perform a residual whiteness test using the residual autocorrelation function. If the autocorrelations between residuals at lag k are within a defined confidence interval, then the residuals are uncorrelated. Adopting a 99% confidence interval (2.58 Gaussian standard deviations) we find the majority of lags (∼95%) fall within the interval. Results are shown in the right panel of Figure 5. One may ignore the lag at k = 0, which is representative of the autocorrelation of a data point with itself.

POLARIS solves a three-dimensional radiative transfer equation for all Stokes parameters:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dl}}{\boldsymbol{S}}=-\widehat{{\boldsymbol{R}}}(\alpha )\widehat{{\boldsymbol{K}}}\widehat{{\boldsymbol{R}}}{\left(\alpha \right)}^{-1}{\boldsymbol{S}}+{\boldsymbol{J}},\end{eqnarray} \tag{ 10 }$$

where $\widehat{{\boldsymbol{R}}}(\alpha )$ is the rotation matrix, $\widehat{{\boldsymbol{K}}}$ is the Muller matrix describing the extinction, absorption, and scattering, respectively, and J is the energy transfer contribution due to emission. The vector ${\boldsymbol{S}}={\left[I,Q,U,V\right]}^{T}$ is the Stokes vector having components I and V representing the total intensity and circular polarization, respectively, and the components Q and U representing the states of linear polarization (for further details on the pipeline in POLARIS, refer to the software manual ⁵ ). We run simulations outputting the intensity of total polarization oriented at θ = 0 (magnetic axis in the plane of the sky) and at θ = π/2 (magnetic axis along the line of sight). We simulate both the R = 0.1 pc and R = 0.2 pc models and the results are given in Figure 8. Additionally, Figure 9 shows the polarization segments rotated by 90°, representative of the inferred magnetic field directions. In our simulations, we obtain polarizations as high as 20% whereas observed polarizations are typically in the range of 4%–6%. This discrepancy is expected since we have enforced perfect alignment of the dust grains, resulting in a more efficient polarization. Had we adopted an imperfect alignment mechanism, we would expect to see lower polarization, but our purpose here is to study the magnetic field morphology. In that case, the perfect alignment assumption provides the most straightforward path.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now use the model to study the stability properties of the core. To compute the normalized mass-to-flux ratio

$$\begin{eqnarray}&&\mu \equiv \displaystyle \frac{(M/{\rm{\Phi }})}{{\left(M/{\rm{\Phi }}\right)}_{\mathrm{crit}}}=2\pi \sqrt{G}\left(\displaystyle \frac{M}{{\rm{\Phi }}}\right),\end{eqnarray} \tag{ 11 }$$

the core's mass and flux must be estimated explicitly. The mass-to-flux ratio gives insight into the contraction mechanisms governing the core formation. We calculate the magnetic flux threading the core, which we can evaluate efficiently in the equatorial plane (z = 0) as

$$\begin{eqnarray}&&{\rm{\Phi }}={\int }_{S}{B}_{z}(r,0)\hat{{\boldsymbol{z}}}\cdot d{\boldsymbol{S}}.\end{eqnarray} \tag{ 12 }$$

Inserting the dimensionless expression introduced for the model magnetic field, we define the dimensionless flux counterpart to be

$$\begin{eqnarray}&&\widehat{{\rm{\Phi }}}=2\pi {\int }_{0}^{\tilde{r}}{\tilde{B}}_{z}(\tilde{r}^{\prime} ,0)\ \tilde{r}^{\prime} d\tilde{r}^{\prime} .\end{eqnarray} \tag{ 13 }$$

Inserting Equation (5) evaluated at $\tilde{z}=0$ we find

$$\begin{eqnarray}\begin{array}{rcl}\widehat{{\rm{\Phi }}} & = & 2\pi {\displaystyle \int }_{0}^{\tilde{r}}\left(\displaystyle \sum _{m=1}^{3}{\beta }_{m}{J}_{0}({a}_{{\rm{m}},1}\tilde{r}^{\prime} )\left[\mathrm{erfc}\left(\displaystyle \frac{{a}_{{\rm{m}},1}\eta }{2}\right)\right.\right.\\ & & \left.\left.+\mathrm{erfc}\left(\displaystyle \frac{{a}_{{\rm{m}},1}\eta }{2}\right)\right]+1\right)\ \tilde{r}^{\prime} d\tilde{r}^{\prime} ,\end{array}\end{eqnarray} \tag{ 14 }$$

which has a closed form solution

$$\begin{eqnarray}&&\widehat{{\rm{\Phi }}}=\pi {\tilde{r}}^{2}\left(\displaystyle \frac{4}{\tilde{r}}\displaystyle \sum _{m=1}^{3}\displaystyle \frac{{\beta }_{m}}{{a}_{{\rm{m}},1}}{J}_{1}({a}_{{\rm{m}},1}\tilde{r})\mathrm{erfc}\left(\displaystyle \frac{{a}_{{\rm{m}},1}\eta }{2}\right)+1\right).\end{eqnarray} \tag{ 15 }$$

Inserting $\tilde{r}=1$ yields the total flux

$$\begin{eqnarray}&&{\widehat{{\rm{\Phi }}}}_{\mathrm{core}}=\pi \left(4\displaystyle \sum _{m=1}^{3}\displaystyle \frac{{\beta }_{m}}{{a}_{{\rm{m}},1}}{J}_{1}({a}_{{\rm{m}},1})\mathrm{erfc}\left(\displaystyle \frac{{a}_{{\rm{m}},1}\eta }{2}\right)+1\right).\end{eqnarray} \tag{ 16 }$$

Since the a_m,1 are roots of J₁, this evaluates to ${\widehat{{\rm{\Phi }}}}_{\mathrm{core}}=\pi$, hence

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{core}}=\pi {R}^{2}{B}_{0}.\end{eqnarray} \tag{ 17 }$$

This result shows that our model is fitting the magnetic field structure within r = R and assuming that the cloud condensed from a uniform background field B₀ with the fixed magnetic flux shown above. We note that the model allows the field lines to become parallel to the z-axis at r = R, i.e., B_r goes to zero at r = R, however, B_z may not converge to B₀ at r = R for z-values near the midplane (see Ewertowski & Basu 2013, Section 4 for further discussion). The normalized mass-to-flux ratio of the model with R = 0.1 pc can be written as

$$\begin{eqnarray}&&\mu =0.8\left(\displaystyle \frac{{M}_{\mathrm{core}}}{3.6\,{M}_{\odot }}\right){\left(\displaystyle \frac{R}{0.1\mathrm{pc}}\right)}^{-2}{\left(\displaystyle \frac{{B}_{0}}{15\mu {\rm{G}}}\right)}^{-1}.\end{eqnarray} \tag{ 18 }$$

For the adopted core mass and radius this means that the core is somewhat subcritical (μ₀ ≲ 1) if we use the estimated 15 μG magnetic field strength for the intercore medium of molecular clouds using the OH Zeeman effect (Thompson et al. 2019). This field strength is also similar to an estimate made by Kandori et al. (2018) using the DCF method and extrapolating a B − ρ relation to the edge of FeSt 1–457. If alternatively we take B₀ = 5 μG, representative of the general interstellar medium (Ferrière 2015), as an estimate for the background field, then the core is somewhat supercritical (μ₀ ≈ 2). Given that the mass ${M}_{\mathrm{core}}$ has an estimated ≈20% uncertainty (Kandori et al. 2005), the most general conclusion is that the core is transcritical (μ₀ ≈ 1) within a factor of about 2. We note that an overall value μ = 1.41 ± 0.38 was estimated by Kandori et al. (2018) using the DCF method, by measuring deviations of the polarization segments from a parabolic fit to the magnetic field lines.

Next, we evaluate Equation (14) for the R = 0.2 pc model, and use $\tilde{r}=0.5$ to evaluate the flux within the observed dense core. This leads to $\widehat{{\rm{\Phi }}}=0.44$ for the region r ≤ 0.5R = 0.1 pc. Therefore, the normalized mass-to-flux ratio within the half-radius (0.1 pc) is

$$\begin{eqnarray}&&\mu =4.6\left(\displaystyle \frac{{M}_{\mathrm{core}}}{3.6\,{M}_{\odot }}\right){\left(\displaystyle \frac{R}{0.2\mathrm{pc}}\right)}^{-2}{\left(\displaystyle \frac{{B}_{0}}{15\mu {\rm{G}}}\right)}^{-1}.\end{eqnarray} \tag{ 19 }$$

The model fit with R = 0.2 pc yields a more decidedly supercritical core. The estimated value of μ is high enough to imply a supercritical core even accounting for a moderate inclination angle of the magnetic axis relative to the plane of the sky and the uncertainty in the core mass. Furthermore, the value of μ would be mildly supercritical even if the background magnetic field strength at 0.2 pc is as low as the 5 μG value (Ferrière 2015) of the general interstellar medium. The magnetic field of the outer region (0.1 pc < r < 0.2 pc) of this model fit may imply a subcritical mass-to-flux ratio; however, the mass in that region is not well constrained observationally. The model we use is free to fit the magnetic field morphology without making any a priori assumptions about the spatial distribution of the mass-to-flux ratio.

The POLARIS simulations yield the degree of polarization as well as its orientation at different locations on the map. These emerge from the full three-dimensional magnetic field structure of our model as well as the adopted density structure and various assumptions about the dust properties listed in Section 2.2 and Table 1. Since the output from POLARIS is a grid of directions for the polarization of dust emission, we identify this with the direction of the elongation of the dust grains in a subsequent discussion. On the other hand, the polarization segments due to dichroic extinction measured by Kandori et al. (2017) are perpendicular to the grain elongation. The left-hand column of Figure 9 shows a comparison of the polarimetry data with the inferred magnetic field from the POLARIS simulation, which are expected to be parallel to one another. This provides a consistency check of whether the three-dimensional magnetic field model that is obtained from fitting the midplane field to the observed polarization segments can also yield a synthetic polarization map that would match the observations. There is a general tendency for the red polarization segments to be close to parallel to the inferred magnetic field direction that is found by rotating the dust grain directions by 90°. We further perform a residual estimate of the relative angles made with the z-axis for the magnetic field directions inferred from the POLARIS output as well as from the best-fit midplane magnetic field model. Because there are more dust grain vectors in the POLARIS output than polarization segments in the Kandori et al. (2017) data set, we perform the calculations using only the closest dust grain vector to each polarization segment. For the polarization segments, we rotate the vectors by 90° to get the inferred magnetic field direction. We calculate the difference each set of vectors make with the z-axis. For the midplane magnetic field model (Equations (4) and (5)), we estimate the angle using the relation

$$\begin{eqnarray*}&&\varphi ={\tan }^{-1}\displaystyle \frac{{B}_{r}}{{B}_{z}},\end{eqnarray*}$$

evaluated at the midpoints of the observed polarimetry data. We find that the mean residual angle of the dust grain directions from the synthetic POLARIS map to be $\overline{\varphi }=0.4081\ \mathrm{rad}$ and the mean residual angle using the model magnetic field to be $\overline{\varphi }=0.4064\ \mathrm{rad}$. The good agreement gives us confidence that we can find good fits to the observed polarization segments using a midplane fit of the magnetic model.