Hourglass Magnetic Field of a Protostellar System

An hourglass-shaped magnetic field pattern arises naturally from the gravitational collapse of a star-forming gas cloud. Most studies have focused on the prestellar collapse phase, when the structure has a smooth and monotonic radial profile. However, most observations target dense clouds that already contain a central protostar, and possibly a circumstellar disk. We utilize an analytic treatment of the magnetic field along with insights gained from simulations to develop a more realistic magnetic field model for the protostellar phase. Key elements of the model are a strong radial magnetic field in the region of rapid collapse, an off-center peak in the magnetic field strength (a consequence of magnetic field dissipation in the circumstellar disk), and a strong toroidal field that is generated in the region of rapid collapse and outflow generation. A model with a highly pinched and twisted magnetic field pattern in the inner collapse zone facilitates the interpretation of magnetic field patterns observed in protostellar clouds.


Introduction
The general picture of star formation involves the gravitational collapse of a dense subregion (a core) of a molecular cloud that is threaded by an ambient magnetic field.The near flux freezing at low densities leads to an hourglass magnetic field configuration in a star-forming core.While predicted theoretically many decades ago [1,2], the observational evidence for hourglass fields emerged more recently [3][4][5].These observations infer a projected magnetic field through dust polarization [5][6][7][8] or dichroic extinction [9].They tend to confirm that gravitational collapse begins with a relatively well-ordered magnetic field and that the turbulent energy does not dominate the magnetic energy within collapsing cloud cores.However, many of these sources are protostellar, which are expected to have a complex three-dimensional magnetic field structure.
The hourglass pattern is a good model for the prestellar phase of collapse, before a central star-disk system has formed.During the prestellar phase, the magnetic field can be modeled as essentially poloidal because the rotation rate and strength of the toroidal magnetic field B ϕ remain small, i.e., B ϕ ≪ B r < B z is the expected hierarchy.In the later protostellar phase, the rapid infall onto the protostar and the accompanying rotational spinup create significant values of B r and B ϕ in the inner regions.The dissipation of the magnetic field by ohmic dissipation and ambipolar diffusion also cause a distortion of the magnetic field pattern.In this phase, B ϕ ∼ B r ∼ B z in some inner regions.This transition can be seen in the simulations of [10] and subsequent works.
An analytic solution for a prestellar hourglass-shaped magnetic field was derived in cylindrical coordinates by [11].They calculated the magnetic field directly from Ampere's Law upon assuming an axisymmetric electric current density.This solution was used successfully to fit the magnetic field in simulations of collapsing cores [11] and the magnetic field pattern inferred from polarized starlight in the prestellar core FeSt-457 [12].In the latter case, ref. [12] found that the radial scale length R of the magnetic field is significantly larger than the observed radial scale of R gas of the gas distribution.This could be understood by the occurrence of significant flux redistribution by ambipolar diffusion in the formation phase of the core.See [12] for details.
In order to evaluate the observed magnetic field of protostellar systems, it is helpful to build an analytic or semi-analytic model for this phase, when the radial and toroidal field components are significant or even dominant in some regions.Prior to this work, the only analytic magnetic field model applicable to the protostellar phase was the magnetized toroid model of [13].That model applied to the pivotal moment when a protostar first forms, and does not include the effects of rotation and the disk and outflow that subsequently develop.Previous work to model dust polarization maps in the protostellar phase has utilized detailed simulation results [14,15], but analytic functions can simplify the analysis and reduce computation time.In this paper, we introduce an axisymmetric semi-analytic model for the protostellar phase that captures many features of simulation results that evolve well into the Class 0 phase and include a disk, outflow, and nonideal MHD [16].These features include an off-center peak in B z and significant amplitudes of B r and B ϕ .

Magnetic Field
An axisymmetric poloidal magnetic field can be derived from Ampere's Law using cylindrical coordinates and by describing the electric current density j = j(r, z) φ.The poloidal field B pol = B r r + B z ẑ can be obtained from the vector potential A = A(r, z) φ.It was shown by [11] that, if the current density is separable, i.e., and where and in which a m,1 is the mth positive root of J 1 (x), and a m,1 < a m+1,1 .Here, R is the outer radius of the solution, and the magnetic field is assumed to equal B 0 = B 0 ẑ for r > R. See [11] for a detailed derivation and description of the boundary conditions adopted in this solution.
Given these results and the relation B = ∇ × A, along with the existence of a background magnetic field B 0 = B 0 ẑ, the axisymmetric poloidal magnetic field components are where B m ≡ k m √ λ m and carries units of the magnetic field.Equation (3) shows that, for any given radial current distribution function f (r), one can evaluate the coefficients k m that yield the poloidal magnetic field components.Note that j(r, z) is the induced current due to self-inductance and (at least partial) flux freezing during the gravitational contraction of a cloud, whereas a current density on a much larger scale is responsible for the background magnetic field B 0 .In [17], various forms of the radial current density function f (r) were used to derive hourglass magnetic field patterns.The adopted forms of f (r) were all centrally peaked, specifically a Gaussian, a Bessel function J 0 , and a smoothed power law.These were mostly applicable to prestellar cores, for the reasons described earlier.Here, we consider more complex forms that capture the physics of the protostar-disk-outflow system.In [16], it was shown through three-dimensional simulations that the early protostellar (Class 0) phase of evolution is characterized by an off-center peak in B z .This is due to the nonideal MHD effects that become significant at the typical density of protostellar disks.Furthermore, there are significant amplitudes of B r in the inner collapse zone due to rapid infall with magnetic field dragging and of B ϕ in the innermost regions of rapid rotation and strong outflow.
We let the radial current distribution function be defined by two peaks and a power law tail as follows: where A and B are scaling factors that carry units of current density, The step functions S 10 and S 01 switch from 1 to 0 or vice versa as r crosses r 1 , ensuring that the vicinity of the second peak is a Gaussian for r < r 1 and a power law for r > r 1 .Here, r 0 and r 1 are the locations of the two separate peaks.The first peak has a Gaussian standard deviation σ 0 .The second peak has a Gaussian inner side with standard deviation σ 1 and an outer side that is a smoothed power law with scale length a 0 and asymptotic dependence ∝ r −2 .
To produce an off-center magnetic field peak that reflects the numerical simulation results in [16], the scaling factors A and B usually take negative and positive values, respectively.While a single centrally peaked Gaussian in the radial current density f (r) yields a centrally peaked poloidal magnetic field [17], a double-peaked f (r) with A < 0 and B > 0 is capable of yielding an off-center peak in the poloidal magnetic field.A value of A < 0 is physically representative of an inner region of magnetic field dissipation; the magnetic field due to the inner negative (counterclockwise direction) electric current partially offsets the magnetic field due to the outer positive (clockwise direction) electric current.The second peak necessarily has a B > 0, and the outer power law decrease in f (r) is representative of an extended region of infall.The location r 1 of the second peak in f (r) sets the approximate location of the peak poloidal field strength on the midplane of the model.We note that, while the function f (r) is an ansatz, the solutions for the magnetic field satisfy Maxwell's equations such that they are self-consistent.
A nonzero toroidal magnetic field that has significant amplitude in the inner region and especially in the outflow zone is a special feature of the protostellar phase.In the earlier prestellar phase, B ϕ is expected to be small, and one can fit a projected magnetic field from polarization maps with a purely poloidal magnetic field model [12].In the protostellar phase, there is significant B ϕ in the rapidly rotating disk region, as well as in the outflow zone, where the inertia of the outflow can twist the ambient magnetic field.The projection of the toroidal field can significantly affect a polarization map of a star-forming region [14,15], so it is vital to include it in a model of the protostellar magnetic field.These can aid in the interpretation of the polarization maps of protostellar outflow systems [6,18].
To model the toroidal field B ϕ , we note that there is a physical correlation between B r and B ϕ within the region of rapid collapse (i.e., inside the expansion wave front-see [19]).The region of strongly pinched magnetic field is also a region of strong twist in the magnetic field due to inward collapse with at least partial angular momentum conservation.The components B ϕ and B r are both antisymmetric across z = 0. Hence, it is reasonable to expect that B ϕ is directly proportional to B r .However, simulations also show that the region of enhanced B ϕ is larger than that of enhanced B r due to twisting of the ambient magnetic field by an outflow [16].We achieve a suitable balance of accuracy and simplicity by introducing a multiplicative dependence on B z as well, and, preferentially in regions somewhat above and below the midplane z = 0, by introducing the vertical scale height H.This leads to the empirical formula where a 1 > 0 is a fitting coefficient.Altogether, our model can be cast in a normalized form for maximum generality and ease of computation.We define normalized quantities r = r/R, z = z/R, η = h/R, as well as Br = B r /B 0 , Bz = B z /B 0 , β m = B m /B 0 .This leads to normalized equations where H = H/R.In addition to picking B 0 as the unit of magnetic field and R as the unit of length, we choose c as the unit of speed.Hence, the unit of current density is B 0 c/R.Therefore, the normalized counterpart of Equation ( 7) is where f , Ã, and B are current densities normalized by B 0 c/R, and all length variables with a "tilde" are normalized by R. Equation ( 12) is used to obtain normalized coefficients for use in Equations ( 9) and (10).

Polarization Model
We use the three-dimensional radiative transfer code POLARIS [20] to produce synthetic polarization maps from our magnetic field model.POLARIS simulates the scattering, absorption, and emission processes of the dust grains within the domain of our model.Details of the radiative transfer technique, including the adopted dust emissivity, absorption, and scattering properties, can be found in [20].We assume a dust-to-gas mass ratio of 0.01, and other adopted dust grain properties are listed in Table 1 of [12].We assume perfect alignment of dust grains through the action of radiative anisotropic torques (RAT); therefore, the calculated polarization percentages may be taken to be upper limits.In this study, we focus on the spatial maps of the variations in polarization percentage rather than absolute values of the same.

Magnetic Field Model
We use the formulation described in Section 2 to find a model that is representative of the protostellar disk-outflow region in the early (Class 0) phase of star formation.In this phase, theoretical calculations [16] show that there is an off-center peak in the magnetic field strength approximately at the edge of the centrifugal disk.The peak magnetic field occurs at the interaction zone between the outwardly diffusing magnetic field of the disk and the inwardly advecting magnetic field of the pseudodisk.
Table 1 lists the parameters used in our magnetic field model.Figure 1 shows the radial dependence f (r) of the current density in our model, including a region of negative values at small radii that captures the effect of local magnetic field dissipation.Figure 2 shows the poloidal magnetic field structure and heat map on three different scales.Dimensional numbers are obtained with an adopted radius R = 600 au and background magnetic field B 0 = 1.62 mG.We choose r 1 = r1 R = 40 au, which is physically associated with the disk edge.The remaining parameters in Table 1 are chosen to give a good approximation of the magnetic field snapshot in Figures 3 and 10 of [16].For the parameter values listed in Table 1, we find that keeping the first 17 terms in the series expansions of Equations ( 9) and ( 10) gives accurate results.A much larger number of terms in the series is needed to be kept than in models of the prestellar phase [12] in order to resolve the small-scale structure of the off-axis magnetic field strength peak.Note that the model is meant to fit the region well within r = R and not the outer collapsing core envelope that is at r > R. With the adopted parameters, the magnetic field profile in the inner few hundred au resembles the magnetic field profiles of the three-dimensional nonideal MHD simulation presented by [16]-see their Figure 7.The left panel of Figure 2 shows that the poloidal magnetic field strength has an off-center peak, which is physically associated with the disk edge.The middle panel shows the transition to a flared magnetic field of the surrounding pseudodisk at r > 50 au.The right panel emphasizes the flared field lines of the infalling pseudodisk.Thus, this model captures the inner region with an outwardly diffused magnetic field as well as the surrounding pseudodisk containing an inwardly advected magnetic field that is dragged in with near flux freezing.12), and we use the parameter values listed in Table 1.The full picture of the magnetic field, affected by rapid rotation and an outflow, is revealed by including the toroidal magnetic field B ϕ .During the protostellar phase, the rapid infall within an outgoing expansion wave will spin up the gas due to near angular momentum conservation.A wide-angle outflow launched primarily from the diskpseudodisk interaction region [16] also transports angular momentum upward, which, in turn, generates significant magnitudes of B ϕ above and below the midplane.Figure 3 shows heat maps of the toroidal field strength and the ratio B ϕ /B pol .These figures capture some key features of nonideal MHD simulations of the disk, pseudodisk, and outflow regions [16]-see their Figure 10.This includes off-center peaks of |B ϕ | above and below the midplane and an increasing ratio |B ϕ /B pol | along the approximate direction of the outflow that emerges from just outside the disk.A visualization of the morphology of the total magnetic field is possible with a threedimensional visualization of the streamlines of the magnetic vectors.Figure 4 shows the magnetic field streamlines of our model when viewed either perpendicular or parallel to the magnetic axis, respectively.Figure 5 shows the corresponding views from angles tilted at 60 • and 45 • relative to the magnetic axis.

Polarization Map
To make further progress in understanding the possible observational realization of our magnetic field model, we generate a synthetic polarization map using the POLARIS radiative transfer code [20].The dust density can be obtained from a gas density profile assuming a dust-to-gas mass ratio of 0.01.
To perform this calculation, we adopt an analytic gas density function that fits the main features of the inner collapse zone and an outflow cavity: This equation is written in the spherical coordinate system (r s , θ, ϕ) and is axisymmetric (ϕ independent).Transformation to a cylindrical coordinate system is accomplished through r s = √ r 2 + z 2 , θ = tan −1 (r/z).We adopt β = −3/2, the typical power law index for gravitational collapse inside an outgoing expansion wave [19].The factor sin n θ accounts for the cloud flattening and outflow cavity, similar to in models of magnetized toroids [13], and the index n can be adapted to fit the outflow zone calculated in simulations.Here, we adopt n = 2.
Figure 6 shows segments of implied magnetic field direction (perpendicular to the emergent polarization direction) and a color map of the polarization percentage for the model with poloidal components only.We adopt R = 600 au and B 0 = 1.62 mG as in previous plots.The two limits of viewing perpendicular to the magnetic axis (edge-on) and along the magnetic axis (pole-on) are shown.These figures reveal the well-known hourglass pattern in the edge-on plot and purely radial segments in the pole-on plot, for which the background field is along the line of sight.We focus here on the morphology of the polarization segments and the relative strength of the polarization rather than on the degree of polarization.Figure 7 is the equivalent of Figure 6, which shows the polarization map of the full magnetic field, including both poloidal and toroidal components.The edge-on map shows a clear transition from an hourglass at outer radii to purely horizontal segments at inner radii, where the toroidal field dominates.There is also a depolarization at some heights above and below the midplane where there is a significant volume of poloidal field in the foreground that is oriented at a large angle relative to the toroidal field in the inner region-see [15] for more discussion of this effect.The pole-on map now emphasizes the circular pattern emerging from the toroidal field.Figure 8 shows the polarization maps of the full magnetic field from two additional viewing angles.From these intermediate viewing angles, the implied magnetic field directions and distribution of polarized intensity are far more complex than in Figure 7.An important feature is the pair of local polarization minima located at the top right and bottom left.There is also a pair of local polarization maxima at the top left and bottom right.The polarization map is not axisymmetric and also does not have a mirror symmetry about the midplane of the observation.This phenomenon was first demonstrated by [14], who explained it in terms of the relative angles of B z and B ϕ in the plane of sky when seen from different angles.These plots, in fact, have a point symmetry (symmetry under rotation by 180 • ).This arises because B z is symmetric about z = 0 while B ϕ is antisymmetric about z = 0. Observational evidence of such point symmetry in a polarization map would be strong evidence for the existence of significant toroidal magnetic fields [14,15,18].These plots also show a complex morphology that could be interpreted in terms of magnetic field line "streamers", i.e., radially directed magnetic field lines that point toward a central "hub".See [22] for an observational example of such a hub-filament system in a region of massive star formation.Qualitative features of such a pattern can, in principle, emerge from a tilted view of a highly twisted but axisymmetric magnetic field.

Discussion
We have demonstrated the ability of an analytic model to capture essential features that emerge from complex three-dimensional nonideal MHD simulations.These include an off-center peak in all three magnetic field components, a significant radial magnetic field near the pseudodisk, and an increasing relative strength of the toroidal field in regions associated with an outflow.A strength of the model is that the field components follow analytic expressions, with adjustable parameters that can be tuned to capture the extent of magnetic dissipation (through a negative current density at small radii near the midplane) and the amount of twisting of the magnetic field in the outflow zone.Here, we have modeled the inner few hundred au squared region as the outflow zone around the protostar since this is where the greatest values of |B ϕ | occur [16].However, the full extent of outflowing gas can reach to at least thousands of au.Future work can explore models with a much larger outer radius R than we adopted here.Such models will need to capture a larger dynamical range of scales and require a larger number of terms in the series solution of Equations ( 9) and (10) than we have adopted.
The analytic forms of the magnetic field and density profile ease the generation of synthetic polarization maps.The POLARIS code can take analytic functions as input, and runs much faster than when working with grid-based data.Analytic functions can also be tuned easily to reflect objects of different sizes and evolutionary stages.
The synthetic polarization maps show the complexity of interpreting observations.The toroidal field makes the projected magnetic field transition from an hourglass at outer radii to a series of horizontal segments at smaller radii.Therefore, the inferred magnetic field direction in these regions is perpendicular to the outflow direction.Observations that show the inferred magnetic field in perpendicular or random directions relative to the outflow direction [23] have been used to suggest that the magnetic field is not dynamically important on scales smaller than ∼ 1000 au.However, the dominance of the toroidal field at small radii means that there is a preference for the projected field to be perpendicular to the outflow when observed from 90 • relative to the magnetic axis.Furthermore, when viewing from an intermediate angle, the relationship is complex, and unusual patterns emerge.
The toroidal field used in our model is an empirical one that captures some of the essential features seen in simulations.Future work can refine the approach or start from the first principles of the current density, as performed for the poloidal field.The model has a peak in B ϕ that is off-center in the radial direction, is not on the midplane z = 0, and has a region of increasing B ϕ /B pol where an outflow is expected to exist.These features are advantageous compared to a previous implementation of B ϕ in the literature for fitting observational data [24] where B ϕ was chosen to be curl free and force free.The model of [24] results in B ϕ having an r −1 dependence, which does not have an off-axis peak and does not match the corresponding profiles emerging from simulations.We also expect that B ϕ is not force free in that it should exert forces that drive the outflow.
Future work can leverage the simplicity of using the analytic forms of the magnetic field and density to make useful comparisons with observed polarization maps.Many of these maps are quite complex but may be probing systems that still have an underlying symmetry.

Conclusions
We introduced an analytic magnetic field model for the protostellar phase.It captures essential physical features that define a protostellar disk-pseudodisk-outflow system.The model agrees well with snapshots of the Class 0 stage as calculated by [16], including an off-center peak in the field strength and strong radial and toroidal field components.The analytic model facilitates the creation of synthetic polarization maps that can be compared with observations.The cases presented here show that observed polarization maps can appear quite complex, leading to difficulty in interpretation.This is the case even though the underlying three-dimensional pattern is quite ordered and contains a significant toroidal field component.

Figure 1 .
Figure 1.The normalized radial current density f (r) versus normalized radius r in our model.The function f (r) comes from Equation (12), and we use the parameter values listed in Table1.

Figure 2 .
Figure 2. Lines of poloidal magnetic field at three different scales, along with a color table of the field strength in units of mG.

Figure 3 .
Figure 3. Left : heat map of the toroidal magnetic field B ϕ .The color bar is in units of mG.Right: heat map of the ratio B ϕ /B pol .

Figure 4 .
Figure 4. Magnetic field streamlines.Left: a view from 90 • relative to the magnetic axis.Right: a view from 0 • relative to the magnetic axis.The colors represent the relative field strength, with red the strongest and dark blue the weakest.The streamlines are plotted using the software package MAYAVI [21].

Figure 5 .
Figure 5. Magnetic field streamlines.The plots are made in the same manner as in Figure 4. Left: a view from 60 • relative to the magnetic axis.Right: a view from 45 • relative to the magnetic axis.The colors represent the relative field strength, with red the strongest and dark blue the weakest.

Figure 6 .
Figure 6.Synthetic polarization maps of the poloidal magnetic field model.The vectors show the implied direction of the magnetic field (rotated by 90 • from the electric field polarization) and the polarization percentage in a color heat map.Left: a view from 90 • relative to the magnetic axis.Right: a view from 0 • relative to the magnetic axis.

Figure 7 .
Figure 7. Synthetic polarization maps of the total (poloidal and toroidal) magnetic field model.The vectors and heat map have the same meaning as in Figure 6.Left: a view from 90 • relative to the magnetic axis.Right: a view from 0 • relative to the magnetic axis.

Figure 8 .
Figure 8. Same as Figure 7 but with different viewing angles.Left: a view from 60 • relative to the magnetic axis.Right: a view from 45 • relative to the magnetic axis.

Table 1 .
Parameter values used in the magnetic field model.