Origin and evolution of magnetic field in PMS stars

Origin and evolution of magnetic field in PMS stars: influence of rotation and structural changes

Abstract

During stellar evolution, especially in the PMS, stellar structure and rotation evolve significantly causing major changes in the dynamics and global flows of the star. We wish to assess the consequences of these changes on stellar dynamo, internal magnetic field topology and activity level. To do so, we have performed a series of 3D HD and MHD simulations with the ASH code. We choose five different models characterized by the radius of their radiative zone following an evolutionary track computed by a 1D stellar evolution code. These models characterized stellar evolution from 1 Myr to 50 Myr. By introducing a seed magnetic field in the fully convective model and spreading its evolved state through all four remaining cases, we observe systematic variations in the dynamical properties and magnetic field amplitude and topology of the models. The five MHD simulations develop strong dynamo field that can reach equipartition state between the kinetic and magnetic energy and even super-equipartition levels in the faster rotating cases. We find that the magnetic field amplitude increases as it evolves toward the ZAMS. Moreover the magnetic field topology becomes more complex, with a decreasing axisymmetric component and a non-axisymmetric one becoming predominant. The dipolar components decrease as the rotation rate and the size of the radiative core increase. The magnetic fields possess a mixed poloidal-toroidal topology with no obvious dominant component. Moreover the relaxation of the vestige dynamo magnetic field within the radiative core is found to satisfy MHD stability criteria. Hence it does not experience a global reconfiguration but slowly relaxes by retaining its mixed stable poloidal-toroidal topology.

Subject headings:
Convection, Hydrodynamics, Magnetohydrodynamics, Stars: interiors, Sun: interior, dynamo, stellar magnetism

1. Stellar evolution and magnetism

Stellar rotation is known to significantly change over the pre main sequence (PMS) through angular momentum conservation as young stars contract under the action of gravitation. At the very beginning of the PMS, stellar rotation remains constant since the star is in a disk-locking phase until about 3 to 10 Myr, when it decouples from the vanishing disk. Then as the star contracts under the influence of gravitation, stellar rotation increases as a consequence of angular momentum conservation until it reaches the zero age main sequence (ZAMS). Later, on the main sequence (MS), stellar rotation decreases as contraction stops and magnetic winds start braking the star. This is not the only drastic evolution that young stars experience during this phase of stellar evolution as their luminosity also varies by a large factor. The internal structure is too strongly impacted as the star evolves along the PMS. Indeed starting from a fully convective state, their radiative zone grows continuously due to the ignition of thermonuclear reactions in their deep core, such as occupying most of their interior upon their arrival on the ZAMS. These major changes impact the star’s properties, especially their internal rotation and magnetic field.

Stellar rotation rate, internal rotation and magnetic field are strongly linked through complex physical processes. At the very beginning of the PMS, stars are meant to rotate quite fast since they contract and accrete angular momentum from the disk. However observations led by Bouvier et al. (1986, 2014) show that they only rotate at one tenth of their break-up velocity. Hence some physical process prevent stars from spinning up at the very beginning of their PMS evolution. Magnetic field is an likely candidate to explain this phenomena as it controls the interaction between the star and its disk. Even after the disk-locking phase, magnetic field has a strong link with rotation through wind-braking and core-envelop coupling. Magnetic field can also possibly modify the transport of angular momentum in stellar interiors through Maxwell stresses. For instance, it has been invoked to explain the flat rotation profile in the radiative interior of the Sun, along with other processes such as internal waves (Charbonnel & Talon, 2005; Eggenberger et al., 2005). It is quite clear that a feedback loop between rotation and magnetic field must exist. On one hand the rotation impacts the magnetic field through dynamo process, and especially through the shearing of magnetic field lines by the differential rotation of the convective envelop (e.g. the -effect). On the other hand, magnetic field topology and amplitude impact braking by the wind. Evidence of such an influence was studied for instance by (Pizzolato et al., 2003; Wright et al., 2011; Gondoin, 2012; Reiners et al., 2014; Matt et al., 2015; Blackman & Owen, 2016). Such analysis showed a correlation between coronal X-ray emission and stellar rotation in late-type main-sequence stars, revealing the existence of two regimes. In the first one, at Rossby number greater than 0.1-0.3, the X-ray emission is well correlated with the rotation period whereas in the second one, at low Rossby number, a constant saturated X-ray to bolometric luminosity ratio is attained. This implies that either the surface field or the stellar dynamo, or both, saturates at fast rotation rates, i.e. at low Rossby numbers.

Stellar magnetic field and internal rotation can also be influenced by internal structure changes. The correlation between the existence of the radiative core and X-ray emission was studied by Rebull et al. (2006). The results showed that stars with a radiative core have values that are systematically lower by a factor of 10 than those found for fully convective stars of similar mass. The flux reduction from fully convective stars to stars with a radiative core is likely related to structural changes that influence the efficiency of magnetic field generation and thus the amplitude and topology of magnetic field. A correlation between the growth of the radiative core and the reduction of the number of periodically variable T Tauri stars have been established by Saunders et al. (2009). Several surface magnetic maps of accreting T Tauri stars have been published (e.g. Donati et al., 2011, 2012; Hussain et al., 2009). These maps were used by Gregory et al. (2012) to study the influence on stellar magnetic fields of the apparition and growth of a radiative core. It has been found that for stars with a massive radiative core, e.g. , the internal magnetic field is complex, non axisymmetric and has weak dipole components. This behavior changes when the radiative core is smaller, e.g. , the field is less complex and more axisymmetric whereas the dipole component is still weak compared to higher order components. As young solar-like stars evolve along the PMS, the internal structure changes from fully convective stars to stars with a radiative core, one can expect similarities with low-mass stars behavior, in particular near the M3-M4 transition. Near the fully convective limit, most of the stars have axisymmetric fields with strong dipole components whereas in fully convective stars, the behavior of stellar magnetic fields might even be bistable with a mixture of different geometries and amplitudes (Morin et al., 2010). Observations have also shown that solar-like stars possess a magnetic field which is predominantly toroidal for fast rotators (Petit et al., 2008; See et al., 2016) and that a subset may even possess a magnetic cycle. To be more precise, some correlations between the period of these cycles and the stellar rotation were advocated during the last decades (Noyes et al., 1984; Soon et al., 1993; Baliunas et al., 1996; Saar & Brandenburg, 1999; Saar, 2001). These analysis show that increases as increases since they found that , where varies depending on the studies but remains positive. However these studies are now reconsidered since they are based on observations of the chromospheric cycle that may differ from the magnetic one See et al. (2016), as it is the case in the Sun (Shapiro et al., 2014). Recent non-linear simulations led by Strugarek et al. (submitted) seem to show that the vs relation may not be so straightforward. Indeed in some of dynamo models, the cycle period decreases while the rotation rate increases (see also Jouve et al., 2010).

Since during the stellar evolution along the PMS, the radiative zone of the star increases, we also wish to know how the magnetic field evolves in the radiative zone as convective dynamo action does not support it anymore. These magnetic fields are observed in massive stars, since their envelop is radiative, where they are often oblique dipoles. In solar-like stars, knowledge of magnetic field in radiative core is important even if it is buried under the dynamo field. Indeed it is a candidate for the transport of angular momentum in the stellar core that can explain the rotation profile observed by helio- and asterioseismology (Schou et al., 1998; García et al., 2007; Benomar et al., 2015). These magnetic fields left by a convective zone in a stably stratified zone are called fossil fields. Studies led by Tayler (1973), Markey & Tayler (1973), Braithwaite (2007) and Brun (2007) showed that purely poloidal or toroidal magnetic fields are unstable in such stably stratified zones. Tayler (1980) proposed that the field needs a mixed configuration to be stable in radiative regions. This statement was confirmed by numerical simulations and theoretical works (e.g. Braithwaite & Spruit, 2004; Duez et al., 2010). Braithwaite (2008) introduced a constraint on the relative amplitude of poloidal and toroidal field in a stable fossil field: . Hence it is interesting to assess if the left over magnetic field is stable or if it must relax to a different configuration.

The origin of stellar magnetic activity and regular cycle is supposed to be linked to a global scale dynamo acting in and at the bottom of the convective envelop Parker (1993). This dynamo is a complex dynamical process that can be understand using the mean field theory Moffatt (1978). Its main ingredients are the -effect, and the helical nature of small scales convective motions, called the -effect (e.g. Parker, 1955, 1977; Steenbeck & Krause, 1969). Alternative mechanisms based on the influence of surface magnetic fields have also been developed, as in the Babcock-Leighton dynamo process (Babcock, 1961; Leighton, 1969; Choudhuri et al., 1995; Charbonneau, 2005; Jouve & Brun, 2007; Miesch & Brown, 2012). These theories allow us to reproduce large scales behavior of the magnetic fields in solar-like stars, such as cycles, but lack an explicit treatment of turbulence and many non linear effects. Thus we cannot rely only on the mean field theory and 3D simulations are an ideal tool to perform such studies. The earliest non-linear, turbulent and self-consistent works on stellar convection and dynamo models in spherical geometry, first with the Boussinesq approximation then with the anelastic hypothesis, were done during the mid 80’s by P. Gilman and G. Glatzmaier (Gilman, 1983; Glatzmaier, 1985a, b). During the last three decades, several groups developed global or wedge-like MHD simulations of convective dynamo, including the ASH code (Brun et al., 2004; Browning et al., 2006; Dobler et al., 2006; Browning, 2008), the PENCIL code (Warnecke et al., 2013; Käpylä et al., 2013; Guerrero et al., 2013), the EULAG code (Ghizaru et al., 2010; Racine et al., 2011; Charbonneau, 2013) and the MAGIC code (Christensen & Aubert, 2006). Observations, theoretical models and 3D numerical simulations enable us to improve our understanding of the MHD processes in solar-like stars. Nowadays we believe that magnetic field in convective envelope of solar-like stars is due to a dynamo with two separate ranges of spatial and temporal scales. The global dynamo explains regular cycles and butterfly diagrams and might be seated at the base of the convective zone and in the tachocline. A local dynamo is likely to be at the origin of the rapidly varying and smaller scale magnetism. All these phenomena take place in the convective zone of the star as all the dynamo theories cited above need convection motions and differential rotation as essential ingredients for regenerating magnetic field.

In our study, we compute 3D global magnetohydrodynamical models of PMS stars at different stages of their early evolution to understand the impact of both structural change and rotation on the internal mean field flows and magnetic field. We choose five different models with specific stellar parameters as presented in section 2. In section 3, we present the hydrodynamical progenitors of our five models. We study the influence of both stellar rotation rate and internal structure on the internal flows and convective motions. Thereafter, we introduce magnetic fields in these HD progenitors. Thus, we can observe the resulting changes on the hydrodynamical flows and convection (see section 4). In section 5, we analyze the amplitude and topology of the magnetic field. We also look at the magnetic field generation in the convective zone and how it evolves as the models go through its evolution and as rotation rate and internal structure change. We finally follow the evolution of the magnetic field in the radiative zone by looking at its stability in the core and its relaxation along the stellar evolution. In section 6, we discuss the results and we conclude.

2. Model setup

To study the evolution of magnetic field during the PMS, we compute 3D global magnetohydrodynamical models of one solar mass star. However stellar evolution in the PMS lasts for several Myrs whereas 3D MHD simulations can only compute stellar evolution for several hundreds of years. Since we cannot compute our simulations for secular time-scales, we select specific models that represents the important stages of the PMS. To characterize these models, we need adequate values for luminosity, rotation rate, radii … These physical values are given by 1D stellar evolution models that were computed with the STAREVOL code (Amard et al., 2016).

2.1. 1D evolution

We study the evolution of a 1 solar-like star during PMS from a fully convective progenitor to the ZAMS. This evolution drastically changes the main stellar parameters: radius, size of the radiative core, rotation rate, luminosity and temperature, as represented in Figure 1. On the upper plot, we can observe that during this evolution the radial structure of the star evolves. At the very beginning of the PMS, no nuclear reactions occur in the core of the young star. The energy of PMS stars is due to the release of gravitational potential energy. Convection is efficient enough to transport energy in the stellar interior. Thus, at the very beginning of the PMS, young stars are fully convective. Since there is no internal process to counterbalance the gravitational contraction, the radius of the star decreases and we can see in Figure 1 that the stellar radius contracts from 2.5 to 1 . As the outer radius becomes smaller, temperature and density increase at the center of the star and the opacity of the core drops as : . When the opacity becomes small enough, the radiative zone appears at the center of the star. We see on the Figure that this radiative zone grows up to 70% of the outer radius as the star reaches the ZAMS and remains stable later on the main sequence.

Figure 1.— Choice of 3D ASH models based on stellar evolution computed with STAREVOL code. Top : Evolution of stellar radius and radius of the radiative core during the PMS and the MS. Middle : Stellar rotation rate as a function of time. At the beginning of the PMS it is constant since star is still in the disk-locking phase. Then it increases as the star contracts under the effect of gravity until the ZAMS. Stellar rotation rate starts decreasing as the stellar contraction stops and magnetic wind brakes the star. Bottom : Relation between luminosity and temperature through stellar evolution. Luminosity decreases, as star ages, for the 4 first models and increases until star arrives at the ZAMS.

As the star contracts, the rotation rate of the star also changes through angular momentum conservation (Gallet & Bouvier, 2013). This evolution shows three different main phases that can be observed in Figure 1 (middle panel). First of all, we see a locking phase with the protostellar disk where stellar rotation remains constant until around . Then the contraction of the star impacts its rotation rate that increases due to angular momentum conservation until the end of the ZAMS, at in our study. In the third and last phase, the star rotation decreases following the Skumanich trend: due to the influence of wind-braking (Réville et al. (2015a)). This modeling of the stellar rotation rate has free parameters such as the initial period at 1Myr, the time coupling between the core and the envelope, the lifetime of the disk and the scaling constant of the wind-braking law. Hence, slow, median and fast rotators can be modeled, as seen in Gallet & Bouvier (2013) and Gallet & Bouvier (2015). For our study, we choose an intermediate rotational evolution profile that comes from STAREVOL models (Amard et al., 2016).

1D simulations, such as STAREVOL models, can compute stellar evolution over secular time-scales. However, they are restricted in space since they do not consider the angular dependencies. On the contrary, 3D global simulations give a more precise understanding of the physical processes that take place in the star, but they do not yet enable the study of star’s life over several million years. To study stellar evolution, we use STAREVOL structures to perform relevant 3D models at key instants in the star’s evolution. We decided to select five models such that the radiative core radii (in stellar radius unit) are well distributed, almost every 20%, and the ratio between the rotation rate of two consecutive models is smaller than two. Such radial structures enable us to create reference states and thus initialize our 3D ASH simulations. We now turn to describe the 3D setup.

2.2. 3D numerical models

Computational methods

The 3D full sphere simulations of the evolution of solar-type stars during the PMS to the ZAMS presented here are computed with the ASH code (Clune et al., 1999; Brun et al., 2004; Alvan et al., 2014). This code evolves the Lantz-Braginski-Roberts (LBR) form of the anelastic MHD equations for a conductive plasma in a rotating sphere (Jones et al., 2011). The anelastic approximation filters fast magnetoacoustic waves but Alfvén and slow magnetoacoustic waves remain. In the ASH code, the equations are non-linear in velocity and magnetic field, and are linearized for thermodynamical variables. These variables are separated into fluctuations and a reference state , which only depends on the radial coordinate and evolves slightly over time. We assume the linearized equation of state:

(1)

with the ideal gas law:

(2)

where , , , have their usual meaning, is the specific heat per unit of mass at constant pressure, is the adiabatic exponent and is the ideal gas constant. The continuity equation is:

(3)

with is the local velocity in spherical coordinates and is the spherical frame rotating at a constant velocity . Under the LBR formulation, the momentum equation can be written as:

(4)

where is the reduced pressure that replaces pressure fluctuations in LBR formulation, is the gravitational acceleration, is the magnetic field and is the viscous stress tensor given by:

(5)

with is the strain rate tensor and the Kronecker symbol. Since we study magneto-hydrodynamical simulations, we need to consider the flux conservation equation for the magnetic field:

(6)

and the induction equation

(7)

with the magnetic diffusivity. The magnetic field and mass flux are decomposed into

(8)

to ensure that they remain divergenceless to machine precision throughout the simulation. Finally the internal energy conservation is:

where and are effective eddy diffusivities that represent momentum and heat transport by subgrid-scale (SGS) motions, is the radiative diffusivity and is the current density. The diffusivity also represents a subgrid process. It is fitted to have the unresolved eddy flux carrying the stellar flux outward the top of the domain. This flux drops exponentially with depth since it should play no role inside the radiative zone. In the energy equation, we have a volume heating term that represents the energy generation by nuclear fusion in the core of the star. We can fit the nuclear reaction rate by a simple model . For each model, we adjust both parameters and such as to have a heating source in agreement with the corresponding 1D STAREVOL model, see Table 1.

Case Stellar radius Radiative radius Luminosity Rotation rate
() () () () ()
FullConv
20 %
40 %
60 %
70 %

Note: Seven stellar parameters that characterize the five stars we choose to model (see Figure 1) : radius, thickness of the convective envelop (), radius of the radiative core, luminosity, rotation rate and and characterizing the nuclear reaction rate : . In all our models, we choose to fix the outer radius of the simulation to 96% of the stellar radius. The simulation has the same internal structure than the Sun. The main difference between this model and the Sun is its rotation rate which is 14 times greater in our model.

Table 1Stellar parameters of our numerical simulations

Problem setup and boundary conditions

ASH is a large eddy simulation (LES) code with a SGS treatment for motions whose scales are below the grid resolution of our simulations. These unresolved scales are modeled with diffusivities , and that represent transport of moment, heat and magnetic field in those small scales. The eddy thermal diffusivity that drives the mean entropy gradient is computed separately and occupies a tiny region in the upper convection zone (dashed plot in Figure 2). This diffusivity transports heat through the outer surface where radial convective motions vanish.

The radial structure of velocity, magnetic and thermodynamical variables is computed with a fourth-order finite differences while angular structure is computed with a pseudo-spectral method with spherical harmonics expansion. Time evolution is solved by a Crank-Nicolson/Adams-Bashorth second-order technique, advection and Coriolis been computed thanks to Adams-Basforth part and diffusion and buoyancy terms thanks to the semi-implicit Cranck-Nicolson scheme (Clune et al., 1999).

The domain of our simulations goes from the center of the star to of the stellar radius for each case considered in this study. Indeed ASH code does not compute 3D simulations up to 100% of the stellar surface since more complex equation of state and very small convection scales would required extreme resolution. Moreover, except for the first one, all our models have two zones with a convective envelope and a radiative core. Since our models are full-sphere, boundary conditions only have to be imposed at the surface of the star and regularization of the solution is done at the center of the star, as described in Alvan et al. (2014). The velocity boundary conditions are impenetrable and torque-free:

(9)

We also fix a constant heat flux

(10)

Finally, we want the surface magnetic field to match an external potential field that implies:

(11)

From 1D stellar structure computed by STAREVOL, a 1D Lagrangian hydrodynamical stellar evolution code (Siess et al., 2000), we initialize the 3D ASH simulations. The gravitational acceleration is fitted from the STAREVOL models. The entropy gradient gives us the internal structure of the star: convection occurs where is negative and for positive, we have the radiative zone (see Figure 3). In our models we impose a small constant negative entropy gradient in the convection zone. The entropy gradient in the radiative zone is deduced from 1D structure as shown in the upper panel of the Figure 3. In this Figure, we see that as the star evolves along the evolutionary path, the radiative zone grows. Moreover it becomes more and more stratified, i.e. the values of entropy gradient are greater.

Figure 2.— Radial luminosity balance for two models : FC (purple) and 20% (blue). Luminosities and radii are normalized to their respective stellar values. In these balances, we show the contribution to the total luminosity (solid line) from radiative diffusion (short dashes), enthalpy (dot-dashed), kinetic energy (three-dot-dashed), modeled SGS processes (long dashes) and viscous diffusion (dot) averaged over 400 days. We notice that in the center of the star, flux luminosity is greater in the model which have a small radiative zone. Indeed it is sustain by the radiative flux (dashed) that is very low in the fully convective star.

We run the simulations over hundreds of convection overturning times and we obtain the flux balance given by the Figure 2. The luminosity flux can be decomposed in several contributions : radiative diffusion, enthalpy, kinetic energy, modeled SGS processes and viscous diffusion. By looking on the blue plots, corresponding to the 20% model, we notice that in the radiative core, the radiative flux is the main contributor to the total flux and the enthalpy flux is negligible whereas in the convective zone, we see that the radiative flux decreases with radius and the convective flux dominates. In the fully convective model (purple lines), the radiative flux is very low and the luminosity is mainly sustain by the enthalpy flux. In both models, the entropy flux is confined to the surface layer and represent the flux carried by the unresolved motions.

Figure 3.— Zoom, on the radiative core and the tachocline, of the radial evolution of the entropy gradient for the five hydrodynamical simulations. Top : In the radiative core, the entropy gradient is positive. As the radiative core grows in the different models, it also becomes more stratified since the entropy gradient amplitude becomes larger. Bottom : In the convective envelop, the entropy gradient is constant and fixed to a small negative value at initial time, then it evolves to become more superadiabatic.
Figure 4.— Kinematic viscosity for the 20% model. In the convective envelop, the diffusivity decreases as the depth increases. The interface between the convective and the radiative zone is characterized by a jump of two order of magnitude of the diffusivity. In the radiative core, we keep a constant diffusivity.
Case (cm) (cm.s) (cm.s)
20 %
40 %
60 %
70 %

Note: Quantities caracterising diffusivity profiles. Radii locate the jump of diffusivity for the radiative core. This jump prevents the spreading of the convective envelope into the radiative core as we compute the simulations. , and respectively represent the values of the diffusivity at the surface, at the base of the convective envelope and in the radiative core.

Table 2Characteristics of diffusivity profiles

For each model, we use the same numerical resolution that gives a horizontal resolution with a maximum spherical harmonic degree . We adapt the radial dependence of diffusivities to accommodate to the coexistence of turbulent convective envelope with stably stratified radiative interior as the star evolves. We use the following formula:

(12)

with and being calculated with the same type of formula. is a step function:

(13)

with the stiffness , the density dependency and remain the same for all the models. This profile of diffusivity is illustrated in Figure 4 for the 20% model. Diffusivity is constant in the radiative core, the interface between the two zones is characterized by a jump of two order of magnitude and the diffusivity slightly decreases into the convective envelop. The radii , that locates the jump in diffusivities, and , that gives the value of the kinematic viscosity at the surface, are referenced in Table 2. We keep a constant Prandlt number, . In the fully convective star, we choose a magnetic Prandlt number at 0.8, which allows us to get dynamo action. When computing the following model, with a small radiative core, we first keep this value of to choose the magnetic diffusivity. However the corresponding was not sufficient enough to trigger dynamo action. Hence we decrease to reach a sufficient which gives us . Since in the last three models this value enables us to have dynamo action in the convective envelop, we choose to keep it and to work at constant for all the models with a radiative core.

Case
FullConv 61.2 1 48.96 0.8
20 % 55.6 1 111.2 2
40 % 43.0 1 86.7 2
60 % 35.5 1 71.0 2
70 % 34.8 1 69.6 2

Note: Characteristic numbers for the MHD simulations. The characteristic numbers are evaluated at mid-depth of the convective zone. These numbers are defined as the Rossby number , the Rayleigh number , the Taylor Number , the Reynolds number , the Prandtl number , the magnetic Reynolds number , the magnetic Prandtl number and the Elsasser number .

Table 3Characteristic numbers
Figure 5.— Radial velocity for the five HD simulations (rows) at three different depths (columns) : , and of stellar radii. Up flows are shown in red/white and down flows in blue/black. At the same depth, the internal structure is not the same in the different models since radiative cores have different sizes. This change, linked with the evolution of stellar rotation rate, leads to important difference for the convection amplitude and patterns.

Table 3 shows us the characteristic numbers for our five models. The Rayleigh numbers in all models are greater than the critical value which is about for the values of Taylor numbers used her (see Jones et al., 2009). The Rossby numbers are significantly smaller than 1. Hence, according to Brun et al. (2017), we can expect the differential rotation profile of our simulations to be solar-like stars with fast equator and slow poles. The Elsasser number is smaller than 1 in all models. As the star ages, this number increases with the Lorentz begin to counterbalance the Coriolis ones.

3. HD progenitors

As seen in Figure 1, we choose five models to represent the evolution of one solar mass star during the PMS. The different parameters of these models are listed in Table 1. We then performed ASH 3D simulations of such model stars for which typically 600kh cpus are needed.

3.1. Internal flow fields

First of all, we want to analyze the impact of the evolution of internal structures and rotation rates in our five hydrodynamical simulations on the convection patterns and internal flow fields. Convection patterns of our five HD simulations evolve as the radiative core grows and the rotation rate increases.

Figure 6.— Radial velocity for the five HD simulations. At the base of the convective envelop, we see a jump in amplitude with values of velocity in the radiative zone that are 100 times, for the 20% model, to 100 000 times smaller than the value observed in the convective zone.
Figure 7.— Equatorial slices of locally normalized radial velocity, e.g. at each depth the velocity is normalized by the horizontally averaged rms velocity. Downflows are shown in blue/black and upflows are in red/white. The radiative core is delimited by a dotted black line. In that radiative core, we can see the gravito-inertial waves that are excited by the downflows of the convective zone hitting the tachocline.

In Figure 5, shell slices of radial velocity are shown at three different depths : , and of the stellar radius. At the surface of our models ((a),(d),(g),(j) and (m)), we observe two types of convective patterns in all simulations with, at low latitudes, elongated flows that are aligned with the rotation axis, the so-called banana cells, and smaller scales at high latitudes. These banana cells produce correlations in the velocity field that increase the Reynolds stresses (Gastine et al., 2014; Brun et al., 2017). The consequence is an acceleration of the equator and a slowdown of the poles that explain the differential rotation profiles of our simulations. Convective scales at the surface become smaller as the star evolves along the PMS (from (a) to (m)). These changes in the horizontal direction are due to the increase of the most instable modes in the convective zone with the rotation rate (Jones et al., 2009). After the surface, we look at a deeper shell in the star, i.e. at 60% of stellar radius. In the first three models, (b), (e) and (h), we are looking in the convective envelop, since the radiative zone is smaller than of their stellar radius. In these simulations, we see that the size of the convective patterns slightly narrows as the star evolves. As at the surface, we see banana cells near the equator. Moreover these patterns are modulated in amplitude. This phenomenon is called active nest and is linked to small Rossby and Prandtl numbers (see Ballot et al., 2007; Brown et al., 2008). The amplitude of convective patterns at that depth slightly decreases compared to the one at the surface. As the radiative zone reaches of the stellar radius (plot (k)), observations at this depth are not longer focused on the convective envelope but on the tachocline of the star. Thus we observe that convection patterns drastically change. Small structures disappear but there is still some persistence of the convective patterns, especially at the equator. Moreover the amplitude of the convection becomes almost four order of magnitude smaller than the one of the previous models (see Figure 6). In the last model (n), we look at a shell in the radiative core of the star and we can see that the last traces of the convective patterns vanish and we are only left with large-scale structures. At the lowest depth, near of the stellar radius, convection patterns are completely different depending on the model considered. In the fully convective model (c), banana cells can be seen. Their horizontal extents are larger than the ones at the surface and at 60%. As the flows go deeper into the stars, they merge and form larger structures. The amplitude of the velocity keep the same order of magnitude than in the other depths even if it decreases slightly. In the first model with a radiative zone, the observed shell slice (f) corresponds to the tachocline of the star. As in the mid-depth cut for the 60% case (k), small convective structures vanish and a large-scale structure emerge. The amplitude of convective patterns is smaller than the one observed at the same depth in the convective zone of the fully convective model. In the last three simulations ((i), (l) and (o)), the observed shell slice is inside the radiative zone. A large-scale structure is predominant in all these models. The radial velocity amplitude is much smaller than the one observed in the first two models.

Radial velocity gives a good picture of the convection flows and patterns that occur in our 3D stellar simulations. In Figure 7, we show an equatorial slice of the radial velocity for each hydrodynamical model. In all the models that have a radiative zone, gravito-inertial waves can be observed in that core. These waves are generated by the downflows of the convective envelop. These flows come from the surface to the base of the convective zone, where they hit the radiative zone and excite waves in there (Brun et al., 2011; Alvan et al., 2014, 2015, for a detailed discussion of such phenomena in 3D). In the convective envelop, changes in convective patterns are observed in the different models. They are due to the changes in rotation rate (Jones et al., 2009) and in aspect ratio size of the convective envelope. At first, we focus on the differences between the first two models since they have the same stellar rotation rate : the main difference between the models is the size and geometry of the convective envelop. By comparing these models, we see that the size of convective patterns do not change much through the change of internal structure and the appearance of a radiative core of 20%. In the following models, rotation rate increases as the star contracts. We notice that both radial and horizontal extents of the convective patterns decrease as the rotation rate increases. Horizontal changes are coherent with the ones observed on the shell slices and with the increase of rotation rate discussed by commenting Figure 5. Changes on the radial direction may be linked to the changes in geometry and size on the convection zone.

Figure 8.— Contour plot-left side: Meridional circulation of HD simulations. Clockwise circulation is represented in red and counter-clockwise in blue. Contour plot-right side: Differential rotation of HD models: FullConv, 20%, 40 %, 60% and 70%. These models have different rotation rates: 3.5, 3.5, 4.7, 8.7 and 14 (in solar units). Differential rotation in the five simulations is cylindrical. Moreover like in the Sun, the rotation profile is prograde, i.e. the equator rotates quicker than the poles. 1D plots: Radial cuts of the differential rotation profile. The radiative cores, when presents, are in solid-body rotation. There is a drop of the rotation rate in the tachocline with respect to the value observed in the radiative core.

As the rotation changes, differential rotation profiles also vary. Figure 8 shows the differential rotation profile of our models in the meridional plane, averaged over longitude and 400 days. All our models have a solar-like profile with a fast prograde equator and slow retrograde poles which is coherent with their Rossby number (see Gastine et al., 2014; Brun et al., 2017). They are also solar-like in the sense that their rotation rate monotically decreases from the equator to the poles, except for a small region around the poles where the average is not stable due to small level arm. Profiles are more cylindrical than the one provided by helioseismology for the Sun. This effect is expected for rapidly rotating stars and linked to Taylor-Proudman constraint (Brun & Toomre, 2002; Brown et al., 2008; Brun et al., 2017). In Figure 8, radial cuts show an important radial shear at low latitudes.

Case
FullConv 18.0 19.4
20 % 75.1 86.3
40 % 27.8 30.6
60 % 56.2 67.5
70 % 89.7 96.7

Note: Differential rotation in nHz with measured near the surface and measured at the equator between the surface and the base of the convective zone. represents the relative latitudinal shear measured near the surface. All the values have been averaged over 400 days.

Table 4Contrasts in differential rotation

Except for the fully convective one, all our models have two zones with a convective envelope and a radiative core. Another characteristic of the solar differential rotation profile is the solid body rotation of the radiative core. In the first model, with a small radiative core of 20% of the stellar radius, the core rotation is not constant. However is does depend only of the radius and no longer of the latitude. In models with a larger core, we see that the core is in solid rotation, except for the overshooting zone where the convective motions penetrate the radiative zone. There is an interface of shear between the differentially rotating convection zone and the radiative interior which is in solid body rotation. This interface is called a tachocline (Spiegel & Zahn, 1992) and plays an important role in the dynamo process (Browning et al., 2006).

A more quantitative analysis of the differential rotation can be achieved by calculating the following quantities : and . is the contrast in differential rotation near the surface between two latitudes : and . is the difference in differential rotation near the equator at two depths : near the surface and at the base of the convective zone. The values of these quantities are listed in Table 4. As the rotation rate of the star increases, from the 40% model up to the model 70%, angular velocity contrast increase whereas the relative shear remains quite constant.

As the Rossby number is small in all our simulations, we expect our meridional circulation to be multi-cells which is the case (Featherstone & Miesch, 2015). Meridional circulation can be represented by plotting the contours of the stream function , as defined by Miesch et al. (2000):

(14)

In Figure 8, cells in in the stellar interior are cylindrical and aligned with the rotation axis for all models. The amplitude of is between 8 and 25 , depending on the model with no clear trend as the star evolves along the PMS. In each hemisphere, there is an changeover of sign between the cylindrical cells with an anti-symmetry with respect to the equator. At the surface, the meridional circulation keeps the same behavior, regardless of the model. It is counter-clockwise in the northern hemisphere and clockwise in the southern, i.e. the flows go from the equator to the pole at the surface and come back to the equator deeper in the star’s interior.

Changes in rotation rate and in geometry of the convective envelop also impact the angular momentum transport. Since we choose stress-free and potential-field boundary conditions at the top of our simulations, no net external torque is applied, and thus angular momentum is conserved. Following previous studies (e.g. Brun & Toomre, 2002), we study the contribution of the different terms in the balance of angular momentum: viscous diffusion, meridional circulation and Reynolds stresses. In all models, Reynolds stresses redistribute the angular momentum outward while the viscous diffusion is inward, down the radial gradient of . The amplitude of the contribution of meridional circulation is smaller and its sign varies with radius following the numbers of dominant cells as seen in Figure 8. Overall the radial balance is well established in the five 3D hydrodynamical simulations. Considering the latitudinal flux balance, Reynolds stresses carry angular momentum towards the equator as they are positive (resp. negative) in the northern (resp. southern) hemisphere. Viscous diffusion is poleward since they tend to erase the differential rotation in the star. Meridional circulation is a response to the torque applied by the sum of the Reynolds and viscous stresses and its sign varies with the multiple cells seen in 8. Latitudinal balance is longer to attain than the radial one and is not fully established in our hydrodynamical simulations.

Case EK KE DRKE MCKE CKE
( erg) ( erg.cm) ( erg.cm) (erg.cm) ( erg.cm)
FC 114 5.52 5.12 (92.7%) 379 4.05 (7.30%) - -
20 % 62.9 6.89 6.39 (92.7%) 457 5.04 (7.30%) 0.5
40 % 4.35 1.11 0.776 (70.2%) 68.2 3.29 (29.8%) 0.22
60 % 4.04 3.11 2.63 (84.5%) 170 4.82 (15.5%) 0.15
70 % 2.01 2.18 1.78 (81.1%) 454 4.11 (18.8%) 0.05

Note: The first column gives the global energy in the convective zone (in erg). The four following columns show the kinetic energy density (KE) split into three components : convection (CKE), differential rotation (DRKE) and meridional circulation (MCKE). The kinetic energy densities KE, DRKE, MCKE and CKE are reported in erg cm. They take into account the changes in size and geometry of the convective zone in the different models. All values are averaged over a period of 400 days.

Table 5Kinetic energies in hydrodynamical simulations

3.2. Kinetic energy

To further assess the dynamics in the five cases, we now turn to study their energetic content. The total energy contained in the convective zone decreases during the PMS as seen in Table 5. However this trend can be due to the decrease of the size of the convective envelop. Thus we look at the kinetic energy density that does not depend on the volume and we notice that there is no clear trend. Kinetic energy density can be split into different components linked to the fluctuating convection (CKE), the differential rotation (DRKE) and to the meridional circulation (MCKE) as done by Miesch et al. (2000). These energies are defined as:

(15)

with KE = DRKE + MCKE + CKE and is the longitudinal average. The values corresponding to these energy densities are reported in Table 5. In the convective zone, the MCKE is negligible compared to CKE and DRKE. The differential rotation has the more important contribution in all models even if this proportion can vary, from 70% in the 40% model to 92.7% in the 20% model. The fraction of the fluctuating convection is smaller but not negligible.

In the radiative zone, proportions are drastically different since the core is stably stratified. Convection and meridional circulation are quite reduced and hence the associated energy densities are negligible compared to the one due to differential rotation. DRKE represents more than 99% of the kinetic energy of all the models possessing a radiative core. We defined two quantities to study the evolution of kinetic energy in both zones: is the ratio between the kinetic energy and is the ratio between the kinetic energy densities. The values of these ratios are given in Table 5. Hence we notice two opposite trends. In total value, the energy in the radiative zone tends to increase more rapidly than the one in the convective zone, even if the energy stored in the convective zone is still much higher than the one in the core. However the sizes of both zones change since the radiative core grows. Thus we also have to look at the ratio of energy densities. At first we notice that the fraction of density energy stored in the radiative core is not negligible, especially when the core is small. Secondly, this fraction tends to decrease as the radiative core grows. To put it in a nutshell, there is more and more energy in the radiative core and its contribution to the total energy, even if it remains small, grows. But by looking at the energy densities, we notice that if the energy ratio grows, the energy density ratio decreases as the radiative core grows: there is more energy in the radiative zone but it is less concentrated.

4. Magnetic field properties and evolution during the PMS phase

In the previous section, we saw how the evolution of the star along the PMS modified its internal structure and flows. We now want to study the impact of this evolution on the resulting internal magnetic field of the star both in the convective and radiative zones.

Figure 9.— Description of the procedure to study the evolution of an stellar magnetic field through the PMS. Once we verified that the hydrodynamical models have equilibrated internal flows and coupling between the radiative core and the convective envelop, we introduce a seed magnetic field in the first model, here the fully convective one. This field is chosen to represent the internal magnetic field of the star after the proto-stellar phase. We run the computation of the fully convective model with this magnetic field until this model is equilibrated, we take its final magnetic field and put it in the following model. Then we re-do all the steps until we reach the end of the PMS with a star with a 70% radiative core.

4.1. The procedure

Figure 9 shows how we proceed to reproduce this evolution with ASH simulations. First of all, we inject a seed magnetic field in the fully convective hydrodynamical model. This weak seed confined dipole magnetic field represents the field left by the proto-stellar phase. We run the MHD simulation of the fully convective model until it reaches an equilibrium state with a dynamo generated field. Then we inject the magnetic field resulting from this simulation into the 20% hydrodynamical model. Hence, we can see how the change of internal structure affects the magnetic field. Once this simulation reaches an equilibrium state, in the statistically stationary sense (see Figure 10), and the magnetic field has relaxed in the radiative core, we introduce the resulting magnetic field in the following hydrodynamical model. By reproducing these operations with all the hydrodynamical models, we can see the influence on the magnetic field of the changes of internal structure and rotation rate caused by stellar evolution.

Figure 10.— Evolution of energies when magnetic field resulting from the fully convective MHD model is introduced into the hydrodynamical simulation of the 20% model.

4.2. Magnetic field generation and evolution

Figure 11.— The two first columns show the shell slices at the surface of radial velocity and radial magnetic field for the five MHD simulations. The last column shows the equatorial slice of for each model.
Figure 12.— 3D views of the 40% model. Top: Radial velocity of the hydrodynamical model. Bottom: Radial velocity of the MHD model with the potential extrapolation of magnetic field outside the star. Upflows in red and downflows in blue. Field lines are color coded with the radial component of the magnetic field.
Figure 13.— Butterfly diagram for the longitudinal component of the magnetic field across the PMS. The five cases are represented on this 2D plot in time and latitude at different radii of 96%, 60% and 20% of the stellar radius (for each model). As we go from one model to another, changing the internal structure but propagating the magnetic field, we can see that structures are conserved.

Magnetic and velocity fields are linked by dynamo action and Lorentz forces. Hence, as we injected the resulting magnetic field of a model into the following hydrodynamical model , this field has to adapt to the new internal structure and flows. Moreover, the injection of the magnetic field also has an influence on the internal structure and flows. Figure 10 shows the evolution of the kinetic and magnetic energies in the 20% model after the injection of the magnetic field resulting from the fully convective dynamo model as we now run it in MHD mode. Kinetic energy drops with the presence of the magnetic field. The decrease is due to a drastic change in the differential rotation profile which will be describe below (see section 4.3). Magnetic energy has a burst as it has to adapt to the new internal structure and flows. Both poloidal and toroidal energies grow before decreasing. After a transient phase, of roughly 2500 days, kinetic and magnetic energies stabilize and a genuine dynamo process occurs in the convective envelop.

Figure 11 shows us radial velocity at the surface of the stars. By comparing it with those observed in 3.1, we see similar patterns with prograde banana cells at the equator. Size of convective patterns do not seem to change much between the HD and MHD simulations. The 3D topology of the radial velocity is illustrated by for the 40% model Figure 12 both in HD and MHD cases. We notice the cylindrical patterns of convection linked to fast rotation in both cases. The tangent cylinder can be seen as velocity is smaller in it. These two figures, 11 and 12, also show the topology of the magnetic field inside and outside the star. The radial component of the magnetic field is shown at the surface of the star and we see well-defined patterns with a growing amplitude as the star grows along the PMS, except for the 70% which has an amplitude of similar to the 40% model. is shown with equatorial slices. In these slices we can notice that there is two different behaviors of magnetic field depending on the zone, convective or radiative, of the star. In the convective zone, the magnetic field, that comes from the dynamo process, is quite turbulent whereas in the radiative core the field relaxes and possesses smoother and larger structures. The 3D view shows us a potential extrapolation of the magnetic field outside the star which is complex, highly non-axisymmetric and exhibits as well extended transequatorial loops.

As we propagate the magnetic field from one model to another, we want to analyze its time evolution through the PMS. Thus we plot, in Figure 13, a butterfly diagram of our complete set of simulations. This diagram is a 2D plot in time and latitude at different radii (96%, 60% and 20% of the stellar radius) of , the longitudinal magnetic field averaged over . At 96% of the stellar radius, i.e. at the surface of the simulation, and at 60% of the stellar radius, we see cycles in the fully convective and 20% models (see the upper and middle panels). A finer study of these cycles will be done in section 5.5 on the cycle period and the sense of the dynamo wave. Sharp transitions occurs as the magnetic field goes from one model to the following, which is coherent with the burst observed in the volume integrated energy analysis. However some magnetic structures are preserved even during the transition between the different simulations. In the lower plot, at 20% of the stellar radius, we observe that the butterfly diagram has a different behavior depending on the models. In the fully convective star, we still see cycling patterns and a propagation of the field from the equator to the poles which is coherent since we are still in the convective zone of the star. When the magnetic field is injected into the 20% model, there is a major change. Indeed, in this model, the observed radius is no longer in the convective envelop, but in the tachocline of the star. We can notice that we do not see any propagation patterns and the amplitude of is higher than in the convective case, likely due to the larger radial shear of in the 20% model. This can be explain as the tachocline is a zone of high shear where global dynamo might be seated. From the 40% model to the 70% one, the observed radius lies in the radiative core. The amplitude is much lower than in the tachocline and less structured than in the fully convective model, i.e. we see not cycling patterns. The magnetic field relax in the radiative zone until the ZAMS.

Hence we see that evolves quite substantially from one model to the other. The source of these changes, and the type of field generation it produces, will be discussed in section 5.

Figure 14.— Contour plot-left side: Meridional circulation of the MHD simulations. Contour plot-right: Differential rotation of the MHD simulations. 1D plots: Radial cuts of the differential rotation profile. Black plots come from the HD simulations and the red plots results from the MHD models. In all models, the differential rotation is flatten by the introduction of magnetic field.

4.3. Mean flows HD vs. MHD

The introduction of magnetic field in the hydrodynamical models strongly impacts the internal flows we studied in the previous section 3.1. One major change is the profile of differential rotation due to the influence of Maxwell stresses (Brun et al., 2004). Comparison between hydrodynamical and MHD simulations are shown in Figure 14. We can see that the presence of magnetic field tends to reduce the latitudinal variation of the differential rotation in the convective envelop. Solid rotation in the radiative core is also altered by magnetic fields. This change is certainly due to diffusion and spreading of the differential rotation of the convective envelope into the radiative zone. The rotation profiles of the MHD simulations remain prograde but less monotonic and contrasted. Changes on structures of the internal flows can be observed through the equatorial slice of the 20% model in Figure 15. In this figure, the left side of the slice comes from the progenitor HD simulation and the right side shows the result of the MHD model. The horizontal extent of the convective patterns do not drastically change when the magnetic field is introduced. On the contrary, we notice that the radial extent of these structures is larger in the MHD models that in the hydrodynamical ones. This property is coherent with the changes in differential rotation profiles shown in Figure 14 and can be seen in each simulation. Indeed, as the differential rotation profiles are flatter in the MHD models, the shearing is smaller and thus the radial extent of the convective structures are more elongated.

Figure 15.— As in Figure 7, we show here the equatorial slices of normalized radial velocity for the 20% HD and MHD simulations, with the downflows represented in blue/dark and upflows in red/white. The HD simulation is shown on the left side while the MHD one is plotted on the right side of the slice. We see that the internal magnetic field does not drastically change the horizontal extent of the convective patterns whereas the radial patterns are larger in the MHD simulation than in the HD one.

Meridional circulation is also impacted by the introduction of the magnetic field. The cells are still cylindrical and aligned with the rotation axis but there is no anti-symmetry at the equator, as it was observed in the hydrodynamical case. A large clockwise cell spread both sides of the equator for the first two models. In the simulations with a larger core, this cell breaks and we see two smaller cells with opposite signs on both sides of the equator.

The change in differential rotation, observed in the MHD simulations, can be understood by the presence of two additional contributions to the way angular momentum is redistributed in the convective shell: e.g. Maxwell stresses and large scaled magnetic torques (Brun et al., 2004). As seen in section 3.1, Reynolds stresses carry angular momentum outward whereas the viscous diffusion is inward. The introduction of the magnetic field modifies this balance as the inward transport in no longer supported by the viscous diffusion, since the differential rotation has been quenched, but by the Maxwell stresses and large scale magnetic torques. In the MHD simulations, the radial balance of angular momentum transport is well established. However the latitudinal balance is not fully achieved yet when we propagate the magnetic field from one model to the following. We choose to stop our simulations when DRKE is mostly constant in time to inject the resulting magnetic field in the following model. The latitudinal balance is quite long to established and for computional resources issues we cannot compute each model until it reaches this balance. However we notice some trends in the latitudinal transport of angular momentum. The Reynolds stresses are still equatorward and the viscous diffusion poleward. The transport linked to Maxwell stresses and large scale magnetic torques is not completely stable but is mainly poleward. The meridional circulation that helps to established the balance varies in signs as it has multiple cells in each hemisphere as seen in Figure 14. Overall, the action of the magnetic field is to quench the angular velocity by reducing both the radial and latitudinal contrast.

5. Explaining the dynamics

5.1. Energy content and radial flux balance

Kinetic energies are also impacted by the injection of the magnetic field in the model. By comparing the energy densities between the hydrodynamical models and the MHD ones, we note that there is small variation, between 3 and 12% (see Table 6). Contrary to the HD simulations, the energy densities are quite the same in all MHD models. Since differential rotation is flatten, the relative influence of the components of the kinetic energy in the convective envelope is changed, as shown in Table 6.

Case EK KE KE DRKE MCKE CKE
() () (%) () () ()
FC 101 4.93 12 4.0 (8.2%) 2.2 () 4.5 (91.8%)
20% 45.2 4.95 4 2.5 (5.1%) 5.1 () 4.7 (94.8%)
40% 19.4 4.95 8 2.6 (5.2%) 4.6 () 4.7 (94.7%)
60% 5.74 4.42 3 0.44 (1.0%) 2.3 () 4.4 (98.9%)
70% 4.17 4.51 5 1.2 (2.6%) 2.8 () 4.4 (97.3%)

Note: The first column gives the global kinetic energy (EK) of the convective zone (in erg). The third column represents the variations of kinetic energy density (KE) between the hydrodynamical and MHD simulations: KE = KEKEKE The kinetic energy density is split into three contributions: convection (CKE), differential rotation (DRKE) and meridional circulation (MCKE) for MHD models. All energy densities are averaged over a period of 400 days and reported in erg cm.

Table 6Kinetic energies of the MHD simulations
Case EM ME TME PME FME
() () () () ()
FC 116 5.62 1.14 11 (19.9%) 1.5 (27.0%) 3.0 (53.1%)
20% 68.4 7.49 1.51 6.5 (8.74%) 1.1 (14.4%) 5.8 (76.9%)
40% 42.0 10.7 2.16 9.5 (8.86%) 2.5 (22.9%) 7.3 (68.2%)
60% 27.3 21.0 4.75 7.3 (3.48%) 3.2 (15.5%) 17 (81.0%)
70% 8.13 8.8 1.95 1.2 (1.37%) 0.5 (5.23%) 8.2 (93.4%)

Note: The first column gives the global magnetic energy (EM) of the convective zone (in erg). The third column represents the ratio of magnetic to kinetic energy. The magnetic energy density (ME) is split into three contributions: poloidal mean energy (PME), toroidal mean energy (TME) and fluctuating energy (FME). All energy densities are averaged over a period of 400 days and reported in erg cm.

Table 7Magnetic energies

In the convective envelop of the MHD simulations, the DRKE decreases compared to the hydrodynamical case and the dominant term becomes the convective one (CKE). As in the hydrodynamical models, the contribution of MCKE remains negligible.

As seen by plotting the temporal evolution of the energy densities, in Figure 10, the magnetic energy varies when the magnetic field is propagated from one simulation to the following one. Thus our MHD models enable us to study how magnetic energy evolves as the magnetic field is propagated along the PMS. First of all, Table 7 shows that magnetic energy grows slightly as the star evolves through the PMS. For a finer analysis, as for the kinetic energy, we split the magnetic energy (ME) into three different components linked to the mean toroidal magnetic energy (TME), to the mean poloidal magnetic energy (PME) and to the fluctuating energy (FME):

(16)

where is the longitudinal mean and ME = TME + PME + FME.

By looking at Table 7, we note that all models are in a superequipartition state with . The ratio even almost reaches 5 for the 60%. Hence in all our models is in an equirepartition state, and even in a superequipartition state for the faster cases. increases as the star ages. Only the 70% model behaves differently with decreasing. As seen previously, the kinetic energy density does not change much in the MHD models. The change observed in the ratio is thus mostly due to the change in magnetic energy density. Indeed, we observe that increases as the star goes along the PMS except for the 70% model. The specificity of this simulation is a change in the stellar luminosity evolution. The size of the radiative core and the rotation rate monotically change during the PMS whereas the luminosity first decreases down to 60% model and then increases until the ZAMS. This can explain a different behavior for the 70% model compared to the other four cases. Indeed, increases with the size of the radiative core and the rotation rate. By looking at Table 7, we note that, in all simulations, magnetic energy in the convective zone is mainly contained in the fluctuating part (FME). The proportion of mean field energy (TME and PME) varies strongly as the star evolves along the PMS. At the very beginning of the PMS, it represents 47% of the total energy quite close to the proportion of 49% found by a study led by Brown et al. (2011). As the star ages, this proportion decreases until 19% for the 60% model (similar to results obtained by Browning (2008)). Finally as the stars arrives on the ZAMS, mean field energy only represents a few percent of the total energy as in the study led by Brun et al. (2004). However, in all these simulations, the mean toroidal energy prevails in the mean energy, whereas in our models, the predominant term is the mean poloidal energy. In summary, the decrease of the mean field energy is coherent with results obtained by Gregory et al. (2012) in which the magnetic field becomes less axisymmetric and more complex as the star ages along the PMS.

5.2. Topology of the magnetic fields

In the previous section, we have studied the evolution of the amplitude of the magnetic field during the PMS and found that ME mostly grows. We will now focus on the topology of these fields. We can express the energy of the magnetic field at the surface as:

(17)

and we can define:

(18)

We define two ratios that characterize the topology (Christensen & Aubert, 2006): the dipole field strength:

(19)

and the ratio between the axisymmetric and non axisymmetric field:

(20)

As in Schrinner et al. (2012), we also defined the local Rossby number of our simulations

(21)

where

(22)

is the mean harmonic degree, with the average over time and radius.

Figure 16.— Evolution of the magnetic field topology with respect to the local Rossby number of our simulations. Top: Mean harmonic degree evolution, as the local Rossby number grows, grows. Middle: Dipole field strength decreases as the local Rossby number increases. Bottom: The amplitude of the axisymmetric field decreases with respect to the non axisymmetric one as the star evolves along the PMS

Increasing radiative core and rotation rate lead to a larger mean harmonic degree (see Figure 16). This result can be understood by considering the geometrical change of the convective zone, e.g. the smaller aspect ratio that favors larger . As seen in section 2.2.2, the Rossby number also grows as the star evolves along the PMS. This result may seem counterintuitive since as faster rotation should lead to smaller Rossby number. However in our cases, rotation rate is not the only stellar parameter to change: the thickness of the convective envelop decreases along the PMS. By looking at the product , we see that it decreases as the star ages. Hence it is logical that the Rossby number increases as the star evolves along the PMS. As the mean harmonic degree and the Rossby number grow, the local Rossby number increases as the star is aging. By plotting as a function of the local Rossby number, Schrinner et al. (2012) noted a transition between the dipolar and multipolar mode at . The local Rossby number of our simulations are around this transition value. We also observe a transition in the topology evolution of the magnetic field along the PMS. As in Schrinner et al. (2012), dipolar components of the magnetic field are weak when . These components are bigger when even if the transition is less strong than the one observed by Schrinner. The amplitude of the axisymmetric part of the magnetic field also decreases as the star evolves along the PMS. We see a transition between the fully convective model, with , the two models with a small radiative core, with , and the two simulations with a big radiative core . These observations are coherent with the different behavior studied by Gregory et al. (2012). By looking to the dipolar field strength, we also find similarities with this study since dipolar components of the fully convective models are greater than in the other simulations. However it has to be tempered since they are still quite weak contrary to what was observed in the study of Gregory et al. (2012). Moreover the dipolar field strength is quite different between the two models with a small core even if they are both supposed to have weak dipolar components.

Figure 17.— Balance between the variation of and generation terms of the toroidal magnetic field for the MHD simulation of the fully convective model. The first 2D plot shows the variation of : . The other plots are the generation terms integrated between and . All the terms are reported in G/s.
Figure 18.— Balance of time-averaged generation terms of the toroidal magnetic energy reported in %. The toroidal field is maintained thanks to the fluctuating compression and to the shear action, both mean and fluctuating. It is destroyed by the fluctuating advection and mean diffusion terms. The mean advection and mean compression terms are not represented since they are negligible compared to the other terms.

5.3. Mean field generation

In order to better understand the dynamics of the magnetic fields, we now turn to the generation of mean field. To study the generation of both mean toroidal and poloidal magnetic fields, we use the decomposition of the induction equation 7 developed by Brown et al. (2011). Thus we get the different contributions of shear, advection and compression for the magnetic field production:

(23)

with representing the production of field by mean shear, production by fluctuating shear, advection by mean flows, advection by fluctuating flows, amplification due to the compressibility of mean flows, amplification due to the compressibility of fluctuating flows and the ohmic diffusion of the mean fields:

(24)

The expression can be directly used to study the importance of the different terms in the generation of the toroidal magnetic field. The time integral of this equation for the longitudinal component gives

(25)

In these models, over the seven physical processes that contribute to the toroidal mean field generation the mean compression has a negligible role and the mean advection and the fluctuating compression have small contributions. Thus in the following analysis we will neglect the mean advection and the mean compression terms since they are negligible in all our simulations. The result of this calculation is shown in Figure 17 for the fully convective model, where and are taken at the maximum and minimum of the dipole component of the magnetic field (see Figure 21). Therefore, the interval captures one magnetic cycle and one magnetic polarity reversal (see Figure 20). In this figure, we notice that is mostly created by the mean shear and destroyed by the mean diffusion.

To have a better understanding of the generation of the toroidal field and to compare the different models, these generation terms are integrated over radius and latitude and we look at the toroidal mean energy:

(26)

The balance of time-averaged generation terms for the toroidal energy is shown in Figure 18 for all cases. We notice that the toroidal magnetic energy is sustained by the fluctuating compression term and by the action of both mean and fluctuating shear. It is annihilate by the mean diffusion and fluctuating advection terms. The mean advection and mean compression terms are negligible compared to the other contributions. The mean diffusion has the main contribution to the destruction of the toroidal field. The fluctuating advection also counteracts the generation of ME, with a smaller contribution. The generation of toroidal field can be split in two quasi-constant contributions : the fluctuating compression (FC) and the shear (FS + MS). As the star ages, the contribution of the mean shear increases while the one of the fluctuating shear decreases, except for the last simulation where it is almost the same. This tells us about the nature of the mean field generation, shear remains critical in these models.

Figure 19.— Time evolution the vector potential through one magnetic cycle of the fully convective model. The three plot on the left show the vector potential at two different times and (see Figure 20), and their difference . The mid plot stands for the sum of the right-hand side. The components of the sum are shown individually with the mean EMF, the fluctuating EMF and the resistive diffusion. e left 2D plot shows the variation of : in G.cm/s.

The production of mean poloidal field is simpler to understand if we represent by its vector potential :

(27)

By uncurling the induction equation 7 once, we obtain the evolution of the potential vector:

(28)

The first term of the equation is the electromotive force (EMF) coming from the coupling of internal flows and magnetic fields and the second is the ohmic diffusion. These terms can also be decomposed into mean and fluctuating components:

(29)

This equation is illustrated for the fully convective case by Figure 19. The main contribution to the generation of poloidal field is the fluctuating EMF and it will be studied in detail in the following section through the analysis of the effect.

5.4. Assessing the relative contribution of dynamo effects

The generation of poloidal magnetic field is dominated by the action of the fluctuating EMF: . This process can also be interpreted through the -effect approximation which is a first order expansion of around the mean magnetic field and its gradient:

(30)

with a rank-two pseudo-vector and a rank-three tensor. In the following, we will neglect the term. However, this will increase the systematic error when estimating the term. Thus a single-value decomposition (SVD) including the -effect has been calculated in order to provide a lower-bound on the systematic error (Augustson et al., 2015). In the following analysis, has been decomposed into its symmetric and antisymmetric components

(31)

with

(32)

In our study, we will focus on the efficiency of the -effect and on the characterization of our dynamo through the relative influence of its regenerating terms. At first, one interesting measure of these dynamo is to quantify the capacity of the convective flows to regenerate mean magnetic fields. This can be evaluated by finding the average magnitude of an estimated -effect relative to the rms value of the non axisymmetric velocity field

(33)

where is the sum of the diagonal elements of the Reynolds stress tensor averaged over time and over all longitudes. If we want to refine the analysis, we can use the equation 33 to provide a measure of the importance of each component of as

(34)

with and . By calculating this matrix, see Table 8, we notice that for the antisymmetric part , the predominant term is always that impacts the poloidal component of the magnetic field. and have the same order of magnitude and are between 3 and 18 times smaller than . By looking at the symmetric part , we see the same trend. The predominant term is either or which both act on the poloidal component of the magnetic field and the smallest term is always which is at least two times smaller than the predominant term. Thus we may conclude that the relative influence of the poloidal field regeneration by the EMF is more important than the toroidal field regeneration. We can quantify this relative influence :

(35)
tensor
0.104 0.114 0.130
FullConv 0.114 0.127 0.147 4.20 12.5
0.082 0.086 0.095
0.133 0.129 0.115
0.120 0.116 0.124 2.08 10.8
0.084 0.087 0.073
0.203 0.126 0.103
0.149 0.097 0.098 3.42 19.1
0.085 0.070 0.069
0.279 0.137 0.061
0.194 0.102 0.056 4.39 21.2
0.077 0.053 0.040
0.188 0.155 0.032
0.099 0.084 0.024 1.54 9.19
0.211 0.171 0.035

Note: Results of the dynamo analysis on our PMS models. The first column represents the tensor with its symmetric: (white background) and antisymmetric: (gray background) portions (see Eq 31). The middle column gives the relative importance of the -effect to the -effect for the regeneration of the toroidal field. The last column quantifies the ratio of the -effect used for the regeneration of the poloidal magnetic field to the one used for the regeneration of the toroidal field.

Table 8 effect

With Table 8, we confirm the predominance of the poloidal field regeneration over the toroidal field regeneration for all our models. This quantity also enables us to see the evolution of this relative influence over all our models. The influence of the poloidal part of in the 20% model is smaller than in the fully convective model. This decrease cannot be due to the rotation rate since the two models have the same stellar rotation rate. The main difference between these models is the size and geometry of the convective envelope. In the following models, with a decrease in size of the convective envelope and an increase of the rotation rate, a trend appears with an increase of the weight of the poloidal part in as the star evolves along the PMS. The 70% model behaves differently because of a change of trend in that impacts .

As seen in equation 23, the term plays an important contribution for the generation of toroidal magnetic field. This term, called the -effect, represents the action of mean shearing, i.e. differential rotation, on the poloidal magnetic field:

(36)

In mean-field theory, the regeneration of the toroidal field can both be due to the -effect, coming from the fluctuating EMF, and to the effect that acts on the poloidal field through differential rotation. In all our models, we note that the regeneration of by the -effect is small, compared to the one of . Therefore, we now want to measure the relative influence of the -effect to that of the -effect, since the toroidal magnetic field can be regenerated through both effects:

(37)

As the rotation rate remains constant and the convective envelope size decreases, the influence of the -effect decreases with respect to the -effect. In the following models, which have a growing radiative core and an increasing rotation rate, the influence of -effect increases with respect to . The poloidal effect also becomes more and more predominant with respect to . This can be understood as the -effect is strongly link to the differential rotation and as we see in section 5.1 that the contrast in differential rotation grows as the star evolves along the PMS. The 70% model behaves differently. Stellar luminosity increases in this model, whereas it decreases in the others. As luminosity increases, we increases both viscous and magnetic diffusivities to keep consistent Reynolds numbers. This increasing magnetic diffusivity possibly explains the behavior of the 70% model.

Figure 20.— Time evolution of the toroidal magnetic field at the surface of the fully convective model. Cycling period is almost 11 years so a magnetic energy cycle of 5.5 years. Azimuthal averages of both toroidal and poloidal magnetic fields are chosen to show the evolution and reversal of magnetic field. At , is almost completely positive. Reversal begins at with increasing negative values in and decreasing patterns amplitude for . At , negative patterns of mostly vanish and has opposite signs in the two hemispheres. At , the reversal is achieved with a fully negative and patterns of have the opposite signs of those observed at . The toroidal field in shown in color, with red for positive values and blue for negative ones. Contours of poloidal field are over plotted with solid line for positive values and dashed line for negative ones. This poloidal field is extrapolated outside the star using the PFSS model which neglects electrical currents in the corona. The lower boundary is given by our MHD simulations and the upper boundary, modeling the effect on the field of the outflowing solar wind, is characterized by an electric current source surface where the field lines are forced to be radial. In this extrapolation we used a source surface radius of (see Schrijver et al. (2003)).

5.5. Time evolution and magnetic cycles

Time evolution of the magnetic field through the PMS, as shown in the time-latitude plot in Figure 13, clearly possesses a cyclic behavior for three out of five cases. For illustrative purposes, we will now discuss one of these cyclic dynamo cases, namely the fully convective one. In Figure 20, we display a zoomed in version of its time-latitude diagram along with meridional cut for specific times samples. We notice that we have a dynamo cycle of almost 11 years, hence commensurable with what we know about the typical length of magnetic cycles in solar-like stars. We also notice a beginning of reversal at the end of the evolution of the second model, the one with a 20% radiative zone and we decide to continue the computation of this model and we see that a cycling dynamo appears.

The temporal evolution of the dipole, quadrupole and octupole moments of the magnetic field can by measured through the amplitudes of the mode of the radial component of the magnetic field . These amplitudes are calculated by the expressions

(38)

with the integral solid angle is at a fixed radius. The amplitudes , and are shown for both models in Figure 21. The measurements are done at of .

Figure 21.— The amplitudes , and are shown for the cycling models : FC, 20% and 70%. We note an increase of the cycle period length. No clear phasing relationship between the three components can be established, except for the 20% case that show that symmetric and anti-symmetric dynamo mode are in opposition of phase. We also notice that in the 70% model there is no dipole reversal, even if it presents a cycling behavior.

As seen above, our simulations show equatorward propagation of the magnetic field (see Figure 13). In three of those models, FC, 20% and 70%, we even see dynamo cycles. The equatorward propagation in kinematic dynamo is generally attributed to the propagation of a dynamo wave. In this theory, the propagation direction of the dynamo wave is given by the Parker-Yoshimura rule (Parker, 1955; Yoshimura, 1975)

(39)

with , is the convective overturning time and the vorticity and is the azimuthal average.

Figure 22.— Latitudinal component of the propagation direction of a dynamo wave (see equation 39), for the fully convective model. Red contours denote positive southward direction. Overplotted in white are isocontours of at 1 and 2 kG, with solid contours being of positive polarity and dashed of negative polarity. We note the good agreement of the Parker-Yoshimura rule with dynamo branch shown in the butterfly diagram of Figure 20.

To further illustrate the nature of dynamo action realized in our simulations, we calculated for one of the models with a cycle dynamo : the fully convective one. In Figure 22, the Parker-Yoshimura rule is respected as (resp. ) in the northern (resp. southern) hemisphere implies that the dynamo wave propagate from the equator to the poles at low latitudes, as in the butterfly diagram displayed in Figure 20 for the same fully convective case.

5.6. Magnetic field evolution in the radiative zone

Case EM ME ME ME ME ME ME
( erg) ( erg.cm) ( erg.cm) ( erg.cm) ( erg.cm) ( erg.cm) ( erg.cm)
FC 5.62 3.05 (54.3%) 2.57 (45.7%)
20 % 2.50 3.40 1.87 (55.15%) 1.52 (44.86%) 7.49 3.59 (47.8%) 3.91 (52.2%)
40 % 22.7 8.46 2.72 (32.13%) 5.74 (67.87%) 10.7 5.60 (52.1%) 5.14 (47.9%)
60 % 121 33.7 19.7 (58.48%) 14.0 (41.52%) 21.0 11.6 (55.5%) 9.29 (44.5%)
70 % 41.3 8.6 6.34 (73.5%) 2.29 (26.5%) 8.8 4.68 (53.2%) 4.12 (46.8%)

Note: The first column gives the global magnetic energy (EM) in the radiative zone (in erg). The following columns show energy densities with the magnetic energy (ME) divided into two its toroidal and poloidal part (ME, ME). All energy densities are averaged over a period of 400 days and reported in erg cm.

Table 9Magnetic energy decomposition

Observations conducted on the activity of massive stars have shown that some of these stars possess a strong surface magnetic fields (Babcock, 1947; Mathys et al., 2001; Donati et al., 2006; Wickramasinghe & Ferrario, 2005; Beuermann et al., 2007; Becker et al., 2003; Aurière et al., 2007). However these fields are drastically different of those observed in convective stars. None of these have very small scales fields, they all presented strong low components (dipole, quadrupole or octupole). Moreover none of these non convective stars have differential rotation which is mostly the case of our radiative cores. Hence the lack of ingredients for a dynamo generation leads us to infer, as Braithwaite (2008), that these magnetic fields in radiative zone are fossil fields, i.e. stable fields that evolves on diffusive time scales. These fields are very sensitive to instabilities and many of them are unstable and disappears quickly. This is coherent with observations since many massive stars do not seem to possess surface magnetic field. These instabilities were studied by Tayler (1973) and Markey & Tayler (1973) which showed that both purely poloidal and toroidal magnetic fields cannot be stable : we need a mixed configuration poloidal-toroidal. Braithwaite (2009) gave a quantitative upper limit to that stable mixed configuration : the poloidal magnetic energy must be less that 80% of the total magnetic energy.

Thus we use the decomposition of the magnetic field:

(40)

This decomposition can be analyzed in both convective and radiative zones (as shown in Table 9). These values are taken averaged over the last 400 days of each simulations. Hence each final magnetic field fulfill the Braithwaite criteria: MEME. We also notice that the magnetic energy is mostly in equipartition between the poloidal and the toroidal parts in the radiative core, but also in the convective envelop. Hence, each time that the radiative core grows from a radius to , the magnetic field “introduced” by this change in size already fulfill the stability criteria. Since the relaxation does not change much the partition of between poloidal and toroidal, the magnetic field logically fulfill the stability criteria. By looking at the table of magnetic energies (table 9), we notice that this limit is actually never reaches (when our models are relaxed). We decided to track the ratio all along our simulations to see if this stability criteria was always fulfilled. Since must be lower than 0.8, we are looking if . Figure 23 shows the evolution of that ratio along the PMS.

Figure 23.— Evolution of the ratio over the PMS.

Hence we notice that even during the transient time between two different internal structures, the stability criteria defined by Braithwaite is fulfilled. As we propagate the magnetic field from one stellar structure to the following one, we need to check that we compute the MHD simulation long enough to enable to relax in the radiative core. Relaxation takes place on several Alfvén times where is a characteristic length scale of the relaxation in the radiative zone. By looking at the internal structure of the model , we have three different areas. defines the convective envelop of the model. is the radiative core of the simulation. This zone can be split in two : that was already radiative in the previous model and that was convective in the previous model and becomes radiative in this model. In the first area, the relaxation process occurred in previous models as it was already radiative. The relaxation process that we look at in the model thus occurs in the portion limited by the radii : and . That is why we choose to define as