The Initial Magnetic Field Distribution in AB Stars

Stars are born with magnetic fields, but the distribution of their initial field strengths remains uncertain. We combine observations with theoretical models of magnetic field evolution to infer the initial distribution of magnetic fields for AB stars in the mass range of 1.6 - 3.4 M$_{\odot}$. We tested a variety of distributions with different shapes and found that a distribution with a mean of $\sim$800 G and a full width of $\sim$600 G is most consistent with the observed fraction of strongly magnetized stars as a function of mass. Our most-favored distribution is a Gaussian with a mean of $\mu$ = 770 G and standard deviation of $\sigma$ = 146 G. Independent approaches to measure the typical field strength suggest values closer to 2 - 3 kG, a discrepancy that could suggest a mass-dependent and bimodal initial field distribution, or an alternative theoretical picture for the origin of these magnetic fields.

Stars are expected to form with a range of initial magnetic field strengths. Unfortunately this Initial Field Distribution (IFD) is difficult to directly observe, as magnetic fields interact with convective regions and can be erased by subsurface convection layers as a star evolves (e.g. Gough & Tayler 1966;Jermyn & Cantiello 2020). Nonetheless, observations of early-type stars (OBA stars) can probe the present-day distribution of magnetic field strengths, and have revealed that it is bimodal. Stars are observed either with strong fossil fields, in excess of ∼ 300G, or ultra-weak fields, with amplitudes below a ∼ 3G, and these two modes appear to be massdependent. Few or no stars are observed with field strengths in between (Aurière et al. 2007;Grunhut et al. 2017;Fossati et al. 2015).
Corresponding author: Eoin Farrell The bimodal field distribution suggests that there are at least two different origin stories for early-type stellar magnetism. Strong magnetic fields may be the relic of fossil field generation during star formation (e.g. Donati & Landstreet 2009) or during a stellar merger (Schneider et al. 2016(Schneider et al. , 2019a, while weak magnetic fields could be the result of ongoing dynamo processes (Cantiello & Braithwaite 2011Jermyn & Cantiello 2021).
Exploring this scenario, Jermyn & Cantiello (2020) suggested that the difference between the two modes could be explained by the strength of the initial magnetic field in a star. Strong enough initial magnetic fields can suppress nearsurface convective motions and evolve just under flux conservation, producing the strong-field mode of the present-day distribution. Below some critical field strength, however, convection begins and can act to erase the near-surface magnetic field, producing the observed weak-field mode. This naturally explains both modes of the distribution as well as the approximate (mass-dependent) field strengths of the magnetic desert between the modes.
If this suggestion holds then whether or not a star will be observed in the strong-field mode can be calculated given its initial magnetic field strength, mass and evolutionary state. Hence, the fraction of strongly magnetized stars as a function of mass should depend only on the IFD. Importantly, this fraction has recently been observed in a volume-limited survey of AB stars by Sikora et al. (2019a,b). Our aim is to infer the Initial Field Distribution. We begin in Section 2 by introducing the observational data (Section 2.1), our forward model of the magnetic field evolution (Section 2.2) and the Approximate Bayesian Computation approach (Section 2.3) that we use to infer the IFD from the observations. We then apply these methods and obtain the IFD in Section 3. We find that a distribution with a mean of ∼ 800 G and a full-width of ∼ 600 G is most consistent with the observations, though we cannot precisely infer the shape of the distribution. We explore uncertainties in our analysis in Section 4 and discuss the implications of this IFD in Section 5. Sikora et al. (2019a) recently produced a volume-limited survey of all identified intermediate mass MS stars within a heliocentric distance of 100 pc as determined using Hipparcos parallaxes. Sikora et al. (2019b) used spectropolarimetry to measure the magnetic field strengths of the magnetic chemically peculiar stars (mCP) in the sample. Figure 1 shows their sample in the Hertzsprung-Russell (HR) diagram.

Observations
We then choose to divide the data into five mass bins. Using less than five bins limits the power of our inference method, while more than five leads to limited sample sizes of just 1 -2 stars in the upper mass bins. In each bin we compute the fraction of strongly-magnetized stars. We identify objects as strongly magnetic if Sikora et al. (2019b) were able to measure a magnetic field and as weak otherwise, associating these with the two modes of the present-day field distribution. Table 1 summarises the number of magnetic stars n B , non-magnetic stars n nonB and fraction of magnetic stars f B = n B /n nonB in the five mass bins. These five values of f B , along with the mass bins, form the data which we use to infer the IFD.

Theory
begin by noting that there is a critical vertical magnetic field strength above which systems are stable to convection.

MacDonald & Petit (2019) compute this as
where is the sound speed, Γ 1 is the first adiabatic index, ∇ ad is the adiabatic temperature gradient, ∇ rad is the radiative temperature gradient, is the density, and Here gas is the gas pressure and is the total pressure. Using this criterion, we model the evolution of the magnetic fields assuming the following: 1. The magnetic field is uniform near the surface of a star.
2. So long as > crit , convection is suppressed and the field evolves according to flux conservation ( ∝ −2 , where is the stellar radius).
3. If the magnetic field is ever < crit , near-surface convective motions begin and quickly erase the surface magnetic field.
We calculated stellar evolution tracks for stars ranging from 1.6 − 3.4 using revision 15140 of the Modules for Experiments in Stellar Astrophysics (MESA, Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019 software instrument. Details on the MESA microphysics inputs are provided in Appendix A. We swept the mass range in increments of 0.2 and used solar metallicity ( = 0.02). The evolutionary tracks of our models are plotted in the HR diagram in Fig. 1. In these calculations we forced the temperature gradient to be ∇ rad in the outer envelope where subsurface convective layers would otherwise form. This allows us to simulate the inhibition of convection by a strong magnetic field in those regions. This choice has a very small effect on the surface properties, as shown by Jermyn & Cantiello (2020). However, it does impact the thermodynamic properties of the subsurface convection regions, which are relevant for computing crit .
We evaluate crit for each model at each point in time using equation (1) and taking the maximum value inside the subsurface convective layers, because in our scenario even a small convective region is enough to erase the surface magnetic field. For the stars observed by (Sikora et al. 2019a) the relevant subsurface convective regions are driven by hydrogen (H CZ) and helium ionization (HeI CZ, HeII CZ) (Cantiello & Braithwaite 2019). We note that the HeI CZ is frequently stabilized by viscosity and thermal diffusion ), so formally we should exclude its contribution to crit . However, our models indicate that the HeI CZ never produces the largest crit among available convection zones. Fig. 2   (2020), we find that the critical field strength decreases with increasing stellar mass. We also find that forcing the temperature gradient to the radiative gradient in the subsurface convective layers typically decreases crit by about 1 % relative to what Jermyn & Cantiello (2020) calculated, though this can increase to 10% for low-mass models at towards the late main-sequence. Stars expand during the main-sequence (MS). Assuming conservation of magnetic flux, the surface magnetic field decrease with radius as ∝ 1/ 2 . Therefore the initial field strength required to produce a strong (super-critical) magnetic field for a given mass and MS fractional age is the cumulative maximum of crit from the upper panel in Fig. 2 scaled by the square of the increase of the radius. This initial crit is plotted in lower panel of Fig. 2, where the increase as a function of MS age is due to the expansion of the star. As in the upper panel, the critical field strength decreases with increasing mass.

Approximate Bayesian Computation
We use an Approximate Bayesian Computation (ABC) method (e.g. Alsing et al. 2018;Sunnåker et al. 2013) to perform likelihood-free inference of the initial field distribution. The ABC method allows us to reconstruct the IFD as a probability density function ( ) = / from the observed fraction of strongly magnetised stars f B ( ★ ) using only forward simulations, free from any likelihood assumptions or approximations. The advantage of using ABC is that it is relatively easy to to simulate mock data based on a given IFD, but it is not feasible to analytically propagate an IFD to an expected distribution of observed fields. We use the elfi python package (Lintusaari et al. 2018) for ABC to perform our simulations, which take place as follows.
First, we choose a parameterised, functional form for ( ). We write this as , ,... ( ), where the subscripts denote the parameterization. For simplicity, we assume that the fossil field distribution is mass independent throughout. We discuss the implications of this assumption later on.
Next, we choose an appropriate prior distribution ( , , ...) over these free parameters. For example, one of our parameterizations is a Gaussian distribution described by a mean and a standard deviation , and we choose flat priors over and .
Finally, we enter a loop: • We sample values for the free parameters of our distribution from the prior.
• We construct a population of stars with initial magnetic field strengths sampling the resulting IFD ( ... ( )). These stars have masses and ages chosen to match those in the observed volume-limited sample described above.
• For each star, we use our stellar evolution models and frozen flux evolution to determine if the magnetic field is weaker than the critical field at any point between the pre-main-sequence and the age of the star inferred from observations. If this occurs, we mark that star as 'not magnetic', assuming that convection will erase the fossil field. Otherwise we mark the star as 'magnetic'.
• We bin the simulated stars by mass and compute the magnetic fraction B in each bin. The critical initial magnetic field required to suppress the formation of subsurface convective layers and maintain a strong magnetic field until a given age as a function of normalised main sequence age. The trend in age is due to expansion during the main sequence.
• We compare B to that reported from observations by computing the sum of the absolute differences between the fractions in each mass bin.
• If this difference is less than some threshold distance, the values of the parameters chosen at the beginning are saved, otherwise they are rejected.
We perform this loop 10 5 times for each functional form of ( ). The accepted samples form our posterior distribution over the free parameters of each functional form. We verified that the accepted posterior distributions did not change substantially when we increased the number of iterations in the loop or the threshold distance between the simulations and observations

INFERRED INITIAL FIELD DISTRIBUTION
Because magnetohydrodynamics is complicated, we do not have a strong physical prior for the functional form of the IFD. For simplicity, we begin by assuming that all stars are born with some non-zero magnetic field and that the distribution is independent of stellar mass. We discuss the impact of relaxing these assumptions later on. We tested the following different forms for the IFD, all of which are summarized in the upper half of To properly compare these different forms we performed our ABC procedure over 10 5 distributions of each form, sampling the distribution parameters from their priors. We chose our acceptance tolerance separately for each form so as to accept the best 10 3 samples, so that the posterior would be well-sampled, and used these to construct the posterior distributions for each parameter. We verified that the posterior distributions we find are not sensitive to this acceptance threshold. Corner plots of these posteriors are provided in Appendix B.
The left panel of Fig. 3 shows examples of all four forms, with their respective parameters set to the means of their respective posterior distributions. The middle panel then compares the distributions from the 100 best simulations for each functional form with the observations. Finally, the right panel shows the distribution of the "distance" between the observations and simulations for each distribution type. The 'distance' is the sum of the distance between each of the observational data points in the middle panel and the corresponding simulated points.
We find that a linear distribution from 300 G to 1500 G is inconsistent with the observed data (Fig. 3, middle). The Table 2. Details of a range of distributions that we tested, including the free parameters, prior distributions, means and standard deviations ( ) of the posterior distributions and the 'distance' between the threshold 1000th best simulation and the observations (D 1000 ).  distribution of distances in the lower panel of Fig. 3 confirms that linear distributions perform poorly compared with the other distributions. The reason for this is that the linear distribution produces many more strongly-magnetized ∼ 1.8M stars than are actually observed. By contrast, the trapezoidal, Gaussian and triangular forms are all able to reproduce the observations within a much tighter tolerance, and visually reproduce the observed trend in B much more closely (Fig. 3,  middle).
The trapezoidal form is a useful one for building intuition about our inference procedure. This form favours a relatively uniform distribution from 500 G to 1100 G. Because the critical magnetic field strength decreases with increasing stellar mass, the lower bound of the distribution is sensitive to higher mass stars and the upper bound of the distribution is sensitive to lower mass stars. The upper bound of the trapezoidal distribution is relatively strongly constrained due to the very low fraction of strongly magnetised ∼ 1.8M stars. The lower bound and slope are less constrained because the fraction of strongly magnetised stars at higher masses just tells us how many stars have field strengths below the lowest crit , not how those stars are distributed.
F . The Gaussian and triangular forms are very similar in shape, and both rather different from the trapezoidal form. They both favour a peak around 800G and a similar width as that of the trapezoidal distribution (e.g. ∼ 130G). This suggests that the initial field distribution could have a peak at 800 G and a full-width of around 600 G, with very few stars born with magnetic fields of less than 500 G or greater than 1100 G.
Given that Gaussian, triangular, and trapezoidal all fit the data similarly well we cannot tell much else about the shape of the distribution. At the same time, the fact that all three forms obtain similar IFDs is encouraging, and tells us that features which all three agree on, like the mean and the width of the IFD, are likely robust.
The trapezoidal and triangular forms, with their sharp kinks and discontinuities, seem less likely to be found in nature than the Gaussian Therefore, given the similarity between their performance and that of the Gaussian, we choose to explore the Gaussian form in more detail below.
The posterior distribution for the Gaussian form has mean = 770 ± 44 G and standard deviation = 145 ± 77 G. These are anti-correlated, such that we see that distributions with higher mean tend to be narrower, and those with higher variance tend to have lower means. Fig. 4 shows samples of distributions drawn from this posterior, as well as the average over these samples, which broadly tells the same story as above: the initial field distribution is peaked around 800 G and has a full-width of around 600 G, with very few stars born below 500 G or above 1100 G.

UNCERTAINTIES
We now consider various uncertainties in our approach. We have assumed that all stars were born with some nonzero initial magnetic field. While this is almost certainly the case, it could well be that the IFD is itself bimodal, and observations of young stars suggest the possibility of a second mode at very weak field strengths (Villebrun et al. 2019).

Bimodal IFDs
To explore this possibility we investigated distributions in which some fraction of stars are born with an initial magnetic field of init = 0, while the rest have initial fields distributed according to a Gaussian form. We refer to this as the Gaussian+ 0 form. Fig. 6 shows the posterior distribution over 0 , , and for this form. We see that while is mostly uncorrelated with 0 , and 0 are strongly correlated. This makes sense: as the fraction of stars with zero initial field strength increases, more stars with strong magnetic fields are needed to match the observed magnetic fraction B . For this reason, regardless of the precise form we choose, any functional form allowing an additional mode at weak field strengths will exhibit such a degeneracy. We cannot exclude such forms, and this degeneracy introduces an additional uncertainty to our models.

Possible Mass Dependence
We investigated the possibility of a dependence of the distribution of initial fossil field on mass. To do this, we tested a Gaussian IFD with a mean that scales linearly with the stellar mass such that = 0 + * . We found that 0 and are strongly degenerate (Fig. 7).
To understand this, note that our approach is primarily using the fact that crit varies with mass, and so differences in the magnetic fraction between different mass bins tell us about Figure 6. Posterior distribution of the mean and standard deviation and fraction of stars born with no strong magnetic field ( 0 ) for a Gaussian initial field distribution. the IFD. As a result if we permit our parameterized IFD to vary with mass we effectively lose our ability to constrain the distribution. What little inferential power we retain comes from the evolution of crit across the main sequence, but we do not have a large enough sample of stars to also bin by age and so this is only minimally constraining. As a result, we emphasize that our results are strongly contingent on the assumption of a mass-independent IFD over the range from 1.6 − 3.4M .

Sensitivity to overshooting
To study the sensitivity of our results to the details of our stellar evolution calculations, we computed a grid of models with a higher value of core convective overshooting ( ov = 0.020 compared to ov = 0.014). Convective overshooting affects the size of stellar cores, and consequently F . the stellar luminosity and lifetime. We find that the posterior distributions of the parameters inferred using stellar models with larger overshooting is consistent with the distribution inferred using the original models to within 1%.

Magnetic Field Evolution
We have assumed that, so long as the magnetic field is strong enough to shut off convection, the field strength evolves according to flux conservation from the beginning of the main sequence phase onward. This neglects a variety of different possible phenomena. On the main sequence the most important are magnetic diffusion, differential rotation, and stellar interactions.
We think it is safe to neglect magnetic diffusion, as the diffusion time across a star is long compared with its mainsequence lifetime (e.g. Braithwaite & Spruit 2017). We likewise suspect that differential rotation can be neglected, as M F 9 super-critical magnetic fields are likely strong enough to inhibit any significant shears from developing on the mainsequence. Stellar interactions, on the other hand, do pose a challenge to our analysis. Mass accretion can potentially enhance surface stellar magnetic fields (Grunhut et al. 2013), and very strongly magnetized stars are believed to be the outcome of stellar mergers (Schneider et al. 2019b). Depending on the incidence of stellar interactions on the main sequence, these processes could impact our ability to infer the IFD at the ZAMS. While the expected fraction of massive stars undergoing interactions during the main sequence is fairly high (∼ 30%, de Mink et al. 2014), it is much smaller for the mass range explored here. This is because the fraction of short period binaries is smaller in these later-type stars (Moe & Di Stefano 2017). In the future it would be useful to carefully quantify the impact of stellar interactions on the magnetization of main-sequence, intermediate mass stars.

DISCUSSION
We have inferred the initial distribution of magnetic field strengths (IFD) from a volume-limited sample of AB stars with magnetic field measurements. To do this we have assumed that the IFD is independent of stellar mass and that convection erases near-surface subcritical magnetic fields. Our favoured distribution is a Gaussian with mean = 770 G and standard deviation = 145 G. We now turn to the astrophysical implications of this distribution. Jermyn & Cantiello (2020) proposed that the observed distribution of magnetic fields -and the "magnetic desert"can be explained in a scenario where some fraction of stars are born with a smooth distribution of initial fields and these fields are then modified by convection. The fact that our preferred IFD reproduces the observed magnetic fraction in AB stars is consistent with this hypothesis.

Magnetic Desert
It is important to note that this did not need to be the case: we could not have successfully fit arbitrary data. For instance this story would have trouble reproducing a non-monotone magnetic fraction as a function of mass, and would need to invoke more complicated assumptions such as a massdependent IFD.

Connection to Star Formation
It is unclear what physical mechanism sets the IFD. It is possible that the IFD reflects some frozen-in flux present in the molecular cloud a star formed from. However, that story runs into trouble explaining the mean of the IFD. Typical molecular clouds are have strong enough magnetic fields that they cannot collapse without shedding magnetic energy (Troland & Crutcher 2008). If just enough magnetic field is shed to allow the cloud to collapse then we should expect initial magnetic fields in approximate equipartition with gas pressure, predicting field strengths which are orders of magnitude stronger than what is observed.
On the other hand if convection indeed erases any magnetic fields that came before then we would expect the magnetic field on the ZAMS to be generated by a convective dynamo as the star descends the Hayashi track. Villebrun et al. (2019) compiled measurements of magnetic fields on the pre-mainsequence in their figure 1. For fully convective stars at the base of the Hayashi track they see a broad distribution of magnetic field strengths ranging from ∼ 100 G to ∼ 800 G. These field strengths are in equipartition with the kinetic energy of convection in the outer layers of a star carrying luminosity ∼ , and occur for stars with luminosities of that order, so it is at least plausible that the these fields are dynamogenerated.
From that point on though the story becomes unclear. These relatively strong fields persist through the development of a radiative interior, though they mostly vanish by the time early-type stars develop convective cores and radiative envelopes (Villebrun et al. 2019), so there could well be nontrivial field evolution happening in the immediate evolution prior to the ZAMS.

Solar Magnetism
If the IFD we have inferred holds at lower masses, we should expect the Sun to have been born with a ∼ 800 G magnetic field, and that field could well have survived below the solar convection zone to this day (Cowling 1953).
So far as we know, the fossil field in the interior of the Sun has not been measured. There are, however, various upper bounds which have been obtained. Reasoning from measurements of the solar oblateness, Friedland & Gruzinov (2004) found an upper bound of 7 × 10 6 G, which is certainly consistent with our prediction. By contrast Boruta (1996) inferred an upper bound of ∼ 30 G in the radiative zone based on the lack of bias in the solar cycle, contingent on some assumptions including that the dynamo is localized to the tachocline. If their assumptions hold then this tight limit is evidence against our IFD, and favours more complicated models like the bimodal Gaussian+ 0 forms we considered in Section 4.
The results of Stello et al. (2016) also support the idea that the interior of stars like the Sun is not strongly magnetized. They used asteroseismology to infer the presence of strong magnetic fields in evolving red giant cores (Fuller et al. 2019a), and claimed that the incidence of core magnetization goes to zero for stars below 1.1 . These stars are the descendants of main sequence stars with radiative cores, while our inference used stars with convective cores on the main sequence. The interaction of convective cores with a F .
large scale magnetic field could be an important piece of the puzzle (e.g. Featherstone et al. 2009).

Compact Remnants
Our results predict that some early-type stars with subcritical surface magnetic fields may have strong fossil fields hiding just beneath their subsurface convection zones (Jermyn & Cantiello 2021). Assuming simple flux conservation , a field of ∼ 800 G would give rise to fields of order 10 6 G in red giant cores and white dwarfs (WD) and 10 15 G in neutron stars (NS) (e.g. Woltjer 1964;Duncan & Thompson 1992). Asteroseismology has revealed magnetic fields of order 10 6 G in about 20% of red giant cores (Fuller et al. 2015;Stello et al. 2016), and roughly 20% of white dwarfs have 10 5−6 G fields (Landstreet & Bagnulo 2019). This does not necessarily mean that all strongly magnetized NS (magnetars) and magnetic WD inherit their fields directly from the IFD. Indeed there are other mechanisms for generating or destroying magnetic fields in the course of stellar evolution (Duncan & Thompson 1992;Spruit 2002b;Fuller et al. 2019b), but the scenario discussed here provides one possible, simple mechanism for the generation of (strongly) magnetized compact remnants.