Constance Emeriau-Viard

Allan Sacha Brun

ORIGIN AND EVOLUTION OF MAGNETIC FIELD IN PMS STARS: INFLUENCE OF ROTATION AND STRUCTURAL CHANGES

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

During stellar evolution, especially in the PMS phase, 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 energies 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.

STELLAR EVOLUTION AND MAGNETISM

Stellar rotation is known to significantly change over the course of stellar evolution (Gallet & Bouvier 2013). Once the star is formed, i.e. at the end of the protostellar phase, it enters the pre-main-sequence phase (PMS). This stage is characterized by a strong contraction of the star under the action of gravitation and the absence of long-term thermonuclear reactions. When H-burning occurs at the center of the star and stellar contraction stops, the young star arrives on the main sequence (MS) where it will spent the majority of its lifetime. The very beginning of this stage, when the H-burning starts, is called the zero-age-main-sequence (ZAMS). This stage also corresponds to a rapid overall increase of the rotation of the star.

At the very beginning of the PMS phase, 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 ZAMS. Later, on the main sequence, 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 imconstance.emeriau@cea.fr,sacha.brun@cea.fr pacted as the star evolves along the PMS phase. 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 stellar evolution, stars are meant to rotate quite fast since they contract and accrete angular momentum from the disk. However observations led for instance by Bouvier et al. (1986Bouvier et al. ( , 2014) ) show that they only rotate at one tenth of their breakup 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 windbraking 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 [START_REF] Charbonnel | EAS Publications Series[END_REF]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;[START_REF] Wright | Cambridge Workshop on Cool Stars, Stellar Systems[END_REF]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 mainsequence 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 L X /L bol 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. 2011Donati et al. , 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. M core > 0.4 M star , the internal magnetic field is complex, non axisymmetric and has weak dipole components. This behavior changes when the radiative core is smaller, e.g. 0 < M core < 0.4 M star , 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 phase, 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 [START_REF] Moffatt | Magnetic field generation in electrically conducting fluids[END_REF]. 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 P cyc increases as P rot increases since they found that P cyc ∝ P n rot , where n 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 P cyc vs P rot 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 evolution track, 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 helioand 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: E pol /E tot < 0.8. 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 [START_REF] Moffatt | Magnetic field generation in electrically conducting fluids[END_REF]. Its main ingredients are the Ω-effect, and the helical nature of small scales convective motions, called the α-effect (e.g. Parker 1955Parker , 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;[START_REF] Choudhuri | [END_REF][START_REF] Charbonneau | Solar and Stellar Dynamos, Saas-Fee Advanced Courses[END_REF]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;Browning 2008), the PENCIL code (Dobler et al. 2006;Warnecke et al. 2013;Käpylä et al. 2013;Guerrero et al. 2013), the EULAG code (Ghizaru et al. 2010;Racine et al. 2011;[START_REF] Charbonneau | Solar and Stellar Dynamos, Saas-Fee Advanced Courses[END_REF]) and the MAGIC code (Christensen & Aubert 2006;Gastine et al. 2014). 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.

MODEL SETUP

To study the evolution of magnetic field during the PMS phase, we compute 3D global magnetohydrodynamical models of one solar mass star. However the PMS stellar evolution 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 phase. 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).

1D evolution

We study the evolution of a 1 M solar-like star during the early phase of stellar evolution 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 phase, 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 phase, 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 R to 1 R . As the outer radius becomes smaller, temperature and density increase at the center of the star and the opacity of the core drops as : κ ∝ ρT -7/2 . 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.

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 4 Myr. Then the contraction of the star impacts its rotation rate that increases due to angular momentum conservation until the end of the ZAMS, at 50 Myr in our study. In the third and last phase, the star rotation decreases following the Skumanich trend: Ω 0 ∝ t -1/2 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.

3D numerical models

Computational methods

The 3D full sphere simulations of the evolution of solartype stars from the proto-stellar phase 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 X and a reference state X, which only depends on the radial coordinate and evolves slightly over time. We assume the linearized equation of state:

ρ ρ = P P - T T = P γ P - S c p , (1)
with the ideal gas law:

P = Rρ T , (2)
where ρ, P , T , S have their usual meaning, c p is the specific heat per unit of mass at constant pressure, γ is the adiabatic exponent and R is the ideal gas constant. The continuity equation is:

∇ • (ρv) = 0 (3) with v = (v r , v θ , v ϕ )
is the local velocity in spherical coordinates and (r, θ, ϕ) is the spherical frame rotating at a constant velocity Ω 0 êz . Under the LBR formulation, the momentum equation can be written as:

ρ ∂v ∂t + (v • ∇) v = -ρ∇ω - ρ S c P g -2ρΩ 0 × v -∇ • D + 1 4π (∇ × B) × B
(4) where ω = P/ρ is the reduced pressure that replaces pressure fluctuations in LBR formulation, g is the gravitational acceleration, B = (B r , B θ , B ϕ ) is the magnetic field and D is the viscous stress tensor given by:

D ij = -2ρν [e ij -1/3 (∇ • v) δ ij ]
(5)

with e ij = 1/2 (∂ i v j + ∂ j v i ) is the strain rate tensor and δ ij the Kronecker symbol. Since we study magnetohydrodynamical simulations, we need to consider the flux conservation equation for the magnetic field:

∇ • B = 0 (6)
and the induction equation

∂B ∂t = ∇ × (v × B) -∇ × (η∇ × B) (7)
with η the magnetic diffusivity. The magnetic field and mass flux are decomposed into

ρv = ∇ × ∇ × (W êr ) + ∇ × (Zê r ) B = ∇ × ∇ × (Cê r ) + ∇ × (Aê r ) (8)
to ensure that they remain divergenceless to machine precision throughout the simulation. Finally the internal energy conservation is:

ρ T ∂S ∂t + v • ∇ S + S = ρ + 4πη c 2 j 2 + 2ρν e ij e ij -1/3 (∇ • v) c/4π (∇ × B)
is the current density. The diffusivity κ 0 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 in the core of the star. We can fit the nuclear reaction rate by a simple model ε = ε 0 T n . For each model, we adjust both parameters ε 0 and n such as to have a heating source in agreement with the corresponding 1D STAREVOL model, see Table 1. During the early phase of stellar evolution, this energy comes from the stellar contraction (Iben 1965) instead of nuclear fusion. However as the star ages the contribution of stellar contraction ε grav to the total energy ε = ε grav + ε nucl decreases (from 99.95% in the C20 model to 3% in the C70 model). As the ε nucl becomes predominant, we notice that the exponent n is no longer 1 but grows to reach 9 for the model on the ZAMS, i.e. the C70 model.

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 κ 0 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 96% 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:

       v r | rtop = 0 ∂ ∂r v θ r rtop = ∂ ∂r v ϕ r rtop = 0 . ( 9
)
Fig. 2.-Radial luminosity balance for two models : C00 (purple) and C20 (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 also fix a constant heat flux

∂S ∂r rtop

= 0 and ∂ S

∂r rtop = cste. (10)
Finally, we want the surface magnetic field B to match an external potential field Φ that implies:

B = ∇Φ and ∇ 2 Φ| surface = 0. ( 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 d S/dr is negative and for d S/dr 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. 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 C20 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 C00 model (purple lines), the radiative flux is very low and the luminosity is mainly sustain by the enthalpy flux. In both models, the en- Note: Seven stellar parameters that characterize the five stars we choose to model (see Figure 1) : radius, thickness of the convective envelop (D = rtop -r bcz ), radius of the radiative core, luminosity, rotation rate and ε 0 and n characterizing the nuclear reaction rate : ε = ε 0 T n . In all our models, we choose to fix the outer radius of the simulation to 96% of the stellar radius. The C70 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.

Fig. 3.-Zoom, on the radiative core and the tachocline, of the radial evolution of the entropy gradient for the five hydrodynamical simulations. Left: 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. Right: 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.

tropy flux is confined to the surface layer and represent the flux carried by the unresolved motions.

For each model, we use the same numerical resolution N r × N θ × N ϕ = 500 × 768 × 1536 that gives a horizontal resolution with a maximum spherical harmonic degree l max = 512. 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:

ν = ν top c 1 + (1 -c 1 ) f (r) ρ ρ 0 m , (12)
with κ and η being calculated with the same type of formula. f is a step function:

f (r) = 1 2 tanh r -r t σ stiff + 1 , (13)
with the stiffness σ stiff = 0.09 × 10 10 cm, the density dependency m = -0.2 and c 1 = 0.01 remain the same for all the models. This profile of diffusivity is illustrated in Figure 4 for the C20 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 r t , that locates the jump in diffusivities, and ν top , that gives the value of the kinematic Note: Quantities caracterising diffusivity profiles. Radii rt 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. νtop, ν rbcz and νrz respectively represent the values of the diffusivity at the surface, at the base of the convective envelope and in the radiative core.

viscosity at the surface, are referenced in Table 2. We keep a constant Prandlt number, P r = ν/κ = 1. In the C00 simulation, 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 P m to choose the magnetic diffusivity. However the corresponding R m was not large enough to trigger dynamo action. Hence we decrease η to reach a sufficient R m which gives us P m = 2. 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 P m constant for all the models with a radiative core.

Table 3 shows us the characteristic numbers for our five models. The Elssasser numbers are small in all our simulations but it increases as the models rotate faster and get closer to a magnetostrophic state. The Rayleigh numbers in all models are greater than the critical value which is about 10 4 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.

HD PROGENITORS

As seen in Figure 1, we choose five models to represent the evolution of one solar mass star during the PMS phase. 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.

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.

In Figure 5, shell slices of radial velocity are shown at three different depths : 96%, 60% and 20% 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 phase (from (a) to (m)). These changes in the horizontal direction are due to the increase of the most instable modes m in the convective zone with the rotation rate Ω 0 (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 60% 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 60% 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 20% of the stellar radius, convection patterns are completely different depending on the model considered. In the C00 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 C60 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 C00 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. In the C70 simulation, we capture a global mode in the radiative zone. The amplitude of this mode is high enough to hide the weaker progressive internal gravity waves that can only be seen near the base of the convective zone. These waves are

Case

Ro Ra Ta Re Pr Re,m Pr,m Λ C00 4.4 × 10 -3 5.9 × 10 7 1.9 × 10 8 61.2 1 48.96 0.8 3.2 × 10 -3 C20 5.1 × 10 -3 4.0 × 10 7 1.2 × 10 8 55.6 1 111.2 2 9.1 × 10 -3 C40 6.4 × 10 -3 6.0 × 10 7 4.5 × 10 7 43.0 1 86.7 2 1.2 × 10 -2 C60 6.6 × 10 -3 2.3 × 10 8 2.9 × 10 7 35.5 1 71.0 2 3.6 × 10 -2 C70 1.7 × 10 -2 2.0 × 10 8 4.1 × 10 6 34.8 1 69.6 2 3.5 × 10 -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 Ro = vrms/(2Ω 0 D), the Rayleigh number Ra = (-∂ρ/∂S)(d S/dr)gD 4 /ν 2 , the Taylor Number Ta = 4Ω 2 0 D 4 /ν 2 , the Reynolds number Re = vrmsD/ν, the Prandtl number Pr = ν/κ, the magnetic Reynolds number Re,m = vrmsD/η, the magnetic Prandtl number Pr,m = ν/η and the Elsasser number Λ = B 2 rms /(8πρΩ 0 vrmsD) (here defined with vrms, not eta).

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. 2014Alvan et al. , 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. 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.

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.

Except for the C00 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 C20 model 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 : ∆Ω lat and ∆Ω r . ∆Ω lat is the contrast in differential rotation near the surface between two latitudes : 0°and 60°. ∆Ω r 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 C40 model up to the model C70, 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):

r sin θ ρv r = - 1 r ∂Ψ ∂θ and r sin θ ρv θ = ∂Ψ ∂r (14)
In Figure 8, cells in in the stellar interior are cylindrical and aligned with the rotation axis for all models. The amplitude of v θ is between 8 and 25 m.s -1 , depending on the model with no clear trend as the star evolves along the PMS phase. 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 convec-tive 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 momen- 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.

Kinetic energy

To further assess the dynamics in the five cases, we now turn to study their energetic content. The total kinetic energy contained in the convective zone decreases during the PMS phase 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:

DRKE = 1 2 ρ v φ 2 MCKE = 1 2 ρ v r 2 + v θ 2 CKE = 1 2 ρ (v r -v r ) 2 + (v θ -v θ ) 2 + (v φ -v φ ) 2
(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 C40 model to 92.7% in the C20 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: EK RZ/CZ is the ratio between the kinetic energy and KE RZ/CZ 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.

MAGNETIC FIELD PROPERTIES AND EVOLUTION DURING THE PMS PHASE

In the previous section, we saw how the evolution of the star along the PMS phase 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.

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 -Description of the procedure to study the PMS evolution of an stellar magnetic field. 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 the C70m model.

from this simulation into the C20 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.

Magnetic field generation and evolution

Magnetic and velocity fields are linked by dynamo action and Lorentz forces. Hence, as we injected the resulting magnetic field of a model n into the following hydrodynamical model n + 1, 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 C20m 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 Figure 12 for the C40 and C40m models. 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 evolution, except for the C70m model which has an amplitude of B r similar to the C40m model. B ϕ 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 phase. 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 B ϕ , 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 both C00m and C20m 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 C20m 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 B ϕ is higher than in the convective case, likely due to the larger radial shear of Ω in the C20m model. This can be explain as the tachocline is a zone of high shear where global dynamo might be seated. From the C40m model to the C70m one, the observed radius lies in the radiative core. The amplitude is much lower than in the tachocline and less structured than in the C00m model, i.e. we see not cycling patterns. The magnetic field relax in the radiative zone until the ZAMS.

Hence we see that B 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.

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 C20 and C20m models 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.

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 7, we show here the equatorial slices of normalized radial velocity for the C20 and C20m 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. 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 scale 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 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. Note: The first column gives the global kinetic energy (EK) of the convective zone (in erg). The third column represents the ratio between the kinetic energy density in the MHD simulations and the kinetic energy density in the HD simulations. 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 -3 . 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 -3 .

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 (see KE MHD /KE HD in Table 6), we note important variations. These are due to the quenching of the differential rotation by the magnetic field. Indeed we observe that variations of kinetic energies are linked to the importance of the contribution of DRKE in the kinetic energy of the hydrodynamical model. In C00 and C20 models, where DRKE contribution is the highest (92.7%), the ratio KE MHD /KE HD is weak (below 0.1) whereas in the C40 model, where DRKE contribution was the weakest (70.2%), the ratio KE MHD /KE HD is important (0.44). Contrary to the HD simulations, the kinetic 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. 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 phase. First of all, Table 7 shows that magnetic energy grows slightly as the star evolves through the PMS phase. 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.

TME = 1 8π B ϕ 2 PME = 1 8π B r 2 + B θ 2 FME = 1 8π (B r -B r ) 2 + (B θ -B θ ) 2 + (B ϕ -B ϕ ) 2 (
By looking at Table 7, we note that all models are in an equipartition state, and even in a superequipartition state for the faster case. M E/KE increases as the star ages. The ratio even almost reaches 5 for the C60m case. Only the C70m model behaves differently with M E/KE decreasing. As seen previously, the kinetic energy sity does not change much in the MHD models. The change observed in the ratio M E/KE is thus mostly due to the change in magnetic energy density. Indeed, we observe that M E increases as the star goes along the PMS phase except for the C70m 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 evolution of the star whereas the luminosity first decreases down to C60m model and then increases until the ZAMS. This can explain a different behavior for the C70m model compared to the other four cases. 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 phase. At the very beginning of the PMS phase, it represents 47% of the total energy quite close to the proportion of 49% found by a study led by Brown et al. (2011) (see also Dobler et al. 2006). As the star ages, this proportion decreases until 19% for the C60m 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 phase.

Topology of the magnetic fields

In the previous section, we have studied the evolution of the amplitude of the magnetic field during the PMS phase 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 ME surf as:

ME surf = ,m ME m Y ,m (17)
and we can define:

     ME = m ME m , ∀ ME m = ME m , ∀m. (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:

f dip ME 1 12 =1 ME
R axi = ME 0 m>0 ME m . ( 20
)
As in Schrinner et al. (2012), we also defined the local Rossby number of our simulations

R o,l = R o ¯ π (21) where ¯ = (v) • (v) v • v (22)
is the mean harmonic degree, with • the average over time and radius.

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 phase. 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 phase. By looking at the product Ω 0 D, 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 phase. As the mean harmonic degree and the Rossby number grow, the local Rossby number increases as the star is aging. By plotting f dip as a function of the local Rossby number, Schrinner et al. (2012) noted a transition between the dipolar and multipolar mode at R c o,l = 0.1. 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 phase. As in Schrinner et al. (2012), dipolar components of the magnetic field are weak when R o,l > 0.1. These components are bigger when R o,l < 0.1 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 phase. We see a transition between the fully convective model, with R axi 0.35, the two models with a small radiative core, with R axi 0.2, and the two simulations with a big radiative core R axi < 0.05. 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. Two different quantities P and M were proposed by Raedler et al. (1990) to study the parity and non-axisymmetry of the field

P = E S -E A E S + E A M = 1 - E 0 E ( 23
)
where E is the magnetic energy, E 0 the magnetic energy contained in the m = 0 modes, E S the magnetic energy in the even + m modes and E A the magnetic energy in the odd + m modes. M can easily be linked to the R axi quantity by M = 1 1+Raxi . P evaluates the dominant magnetic field parity in the dynamo solution. P fluctuates around the value of zero (indicating exact field parity distribution) with min/max value being around -0.2 to 0.2 most of the time. Our dynamo being non linear, it is no surprises that both family are equally excited (Tobias 1997;DeRosa et al. 2012).

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, advec-tion and compression for the magnetic field production:

∂ B ∂t = P MS +P FS +P MA +P FA +P MC +P FC +P MD (24)
with P MS representing the production of field by mean shear, P FS production by fluctuating shear, P MA advection by mean flows, P FA advection by fluctuating flows, P MC amplification due to the compressibility of mean flows, P FC amplification due to the compressibility of fluctuating flows and P MD the ohmic diffusion of the mean fields:

P MS = ( B • ∇) v , P FS = (B • ∇) v , P MA = -( v • ∇) B , P FA = (v • ∇) B , P MC = ( v r B ) ∂ ∂r lnρ
,

P FC = ( v r B ) ∂ ∂r lnρ
,

P MD = -∇ × (η∇ × B ). (25)
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

∆ B ϕ = t4 t1 dt (P MS + P FS + P MA + P FA + P MC + P FC + P MD ) | ϕ (26)
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 C00m model, where t 1 and t 4 are taken at the maximum and minimum of the dipole component of the magnetic field (see Figure 21). Therefore, the interval [t 1 ; t 4 ] captures one magnetic cycle and one magnetic polarity reversal (see Figure 20). In this figure, we notice that B ϕ 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:

dTME dt = V dV ∂ ∂t B 2 ϕ 8π = V dV B ϕ 4π • (P MS + P FS + P MA + P FA +P MC + P FC + P MD ) | ϕ . (27)
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 other plots are the generation terms integrated between t 1 and t 4 (see Figure 20). All the terms are reported in G/s. 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.

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

B P = ∇ × A ϕ φ. ( 28
)
By uncurling the induction equation 7 once, we obtain the evolution of the potential vector:

∂ A ϕ ∂t = v × B | ϕ -η∇ × B | ϕ . ( 29
)
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:

E MI = v r B θ -v θ B r E FI = v r B θ -v θ B r E MD = -η 1 r ∂ ∂r (r B θ ) - ∂ B r ∂θ . ( 30
)
This equation is illustrated for the C00m 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:

E FI = E = v × B .
This process can also be interpreted through the α-effect approximation which is a first order expansion of E around the mean magnetic field and its gradient:

E i = α ij B j + β ijk ∇ B + O ∂ B /∂t, ∇ 2 B (31)
with α ij a rank-two pseudo-vector and β ijk 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 [START_REF] Augustson | [END_REF]. In the following analysis, α has been decomposed into its symmetric and antisymmetric components

α B = α S B + γ × B (32)
with 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

α S = α (rr) α (rθ) α (rϕ) α (rθ) α (θθ) α (θϕ) α (rϕ) α (θϕ) α (ϕϕ) and γ = γ r γ θ γ ϕ . (33)
E α v rms = 3 2(r 3 2 -r 3 1 ) × i,j drdθr 2 sin θ α ij α ij {v • v } (34)
where {v • v } 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 34 to provide a measure of the importance of each component of α as

ε ij = E ij E 1 E α ij v rms = 3 2E(r 3 2 -r 3 1 ) drdθr 2 sin θ α ij α ij {v • v } = ε (rr) ε (rθ) ε (rϕ) ε γϕ ε (θθ) ε (θϕ) ε γ θ ε γr ε (ϕϕ) (35) with ε (xx) = α (xx) E and E γx = γ x E
. 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. γ r and γ θ have the same order of magnitude and are between 3 and 18 times smaller than γ ϕ . By looking at the symmetric part α S , we see the same trend. The predominant term is either α (rr) 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 α P /α ϕ :

α P α ϕ = 3 2(r 3 2 -r 3 1 ) × drdθr 2 sin θ B P • ∇ × E B ϕ φ • ∇ × E .
(36) With Table 8, we confirm the predominance of the Note: Results of the dynamo analysis on our PMS models. The first column represents the α tensor with its symmetric: αs (white background) and antisymmetric: γ (gray background) portions (see Eq 32). 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.

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 C20m model is smaller than in the C00m 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 phase. The C70m model behaves differently because of a change of trend in L * that impacts v rms .

As seen in equation 24, the P M S 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:

P M S = r sin θ B P • ∇Ω. (37)
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 B ϕ by the αeffect is small, compared to the one of B pol . 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:

Ω α ϕ = 3 2(r 3 2 -r 3 1 ) × drdθr 2 sin θ r sin θ B ϕ B P • ∇ Ω B ϕ φ • ∇ × E .
(38) 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 phase. The C70m 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 C70m model.

Time evolution and magnetic cycles

Time evolution of the magnetic field through the PMS phase, 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 C20m one, 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 Azimuthal averages of both toroidal and poloidal magnetic fields are chosen to show the evolution and reversal of magnetic field. At t 1 , B pol is almost completely positive. Reversal begins at t 2 with increasing negative values in B pol and decreasing patterns amplitude for Btor. At t 3 , negative patterns of Btor mostly vanish and B pol has opposite signs in the two hemispheres. At t 4 , the reversal is achieved with a fully negative B pol and patterns of Btor have the opposite signs of those observed at t 1 . 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 rss 2.5R * (see Schrijver et al. (2003)). t=0) in the figures correspond to the beginning of each MHD simulation. We note an increase of the cycle period length. No clear phasing relationship between the three components can be established, except for the C20m case that show that symmetric and anti-symmetric dynamo mode are in opposition of phase. We also notice that in the C70m model there is no dipole reversal, even if it presents a cycling behavior.

through the amplitudes of the m = 0 mode of the radial component of the magnetic field B r . These amplitudes are calculated by the expressions

D = 3 4π B r cos θ sin θdθdφ, Q = 1 2 5 4π B r 3 cos 2 θ -1 sin θdθdφ, O = 1 2 7 4π B r 5 cos 3 θ -3 cos θ sin θdθdφ, (39)
with the integral solid angle is at a fixed radius. The amplitudes D, Q and O are shown for both models in Figure 21. The measurements are done at 96% of r top .

As seen above, our simulations show equatorward propagation of the magnetic field (see Figure 13). In three of those models, C00m, C20m and C70m, 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)

S = -r sin θ ᾱê ϕ × ∇ Ω Ω 0 (40)
with ᾱ = -τ 0 v • /3, τ 0 is the convective overturning time and = ∇ × v the vorticity and • is the azimuthal average. To further illustrate the nature of dynamo action realized in our simulations, we calculated S θ for one of the models with a cycle dynamo : the fully convective one. In Figure 22, the Parker-Yoshimura rule is respected as S θ < 0 (resp. S θ > 0) 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 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 l 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. The fossil fields in early stellar evolution have been studied by [START_REF] Arlt | Cosmic Magnetic Fields: From Planets, to Stars and Galaxies[END_REF]. 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:

ME = ME pol + ME tor . (41)
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: ME pol RZ < 0.8ME RZ . 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 r bcz,n to r bcz,n+1 , the magnetic field "introduced" by this change in size already fulfill the stability criteria. Since the relaxation does not change much the partition of B 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 B pol /B tor all along our simulations to see if this stability criteria was always fulfilled. Since M E pol /M E must be lower than 0.8, we are looking if B pol /B tor < 2. Figure 23 shows the evolution of that ratio along the PMS phase. 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 B to relax in the radiative core. Relaxation takes place on several Alfvén times t a = v a /D RZ where D RZ MEtor CZ (10 32 erg) (10 4 erg.cm -3 ) (10 4 erg.cm -3 ) (10 4 erg.cm -3 ) (10 5 erg.cm -3 ) (10 5 erg.cm -3 ) (10 5 erg.cm Note: The first column gives the global magnetic energy (EM) in the radiative zone (in erg). The following columns show energy densities the magnetic energy (ME) divided into two its toroidal and poloidal part (ME tor , ME pol ). All energy densities are averaged over a period of 400 days and reported in erg cm -3 . is a characteristic length scale of the relaxation in the radiative zone. By looking at the internal structure of the model n + 1, we have three different areas. r > r bcz,n+1 defines the convective envelop of the model. r < r bcz,n+1 is the radiative core of the simulation. This zone can be split in two : r < r bcz,n that was already radiative in the previous model and r bcz,n < r < r bcz,n+1 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 n+1 model thus occurs in the portion limited by the radii : r bcz,n and r bcz,n+1 . That is why we choose to define D n as D n = r bcz,n+1 -r bcz,n rather than as the radius of the radiative zone of the model. The values of Alfvén times for our simulations are given in Table 10. We see that each MHD model is evolved over few Alfvén times to insure the relaxation of the magnetic field in the radiative zone. We conclude that interestingly dynamo action tends to generate mixed fields whose properties satisfy stability criteria in stratified radiative core. This result on the stability of the fossil field left over by dynamo action as the convective envelope becomes shallower along the PMS phase is a direct outcome of our set of simulations and could not be easily anticipated. It has direct consequences on the geometry of fossil field that can be expected in solar-like star's radiative core as stable mixed poloidal-toroidal configurations should be favored.

DISCUSSION AND CONCLUSION

During the PMS phase, between the protostellar phase and the ZAMS, the stellar radius decreases due to gravitational contraction. As the star contracts, internal temperature and pressure increase, opacity drops in the stellar core and a radiative zone appears and grows within the star. Moreover, the stellar contraction causes an increase of the rotation rate due to angular momentum conservation. Hence we expect the star's dynamical properties to vary significantly and these variations are key to characterize. In order to do so, we have developed a series of 3D MHD simulations of stellar convective dynamo. To make this study more realistic, we used to setup the spherically symmetric background state of our 3D models, radial profiles obtained from 1D stellar evolution model at various stages along the PMS evolutionary track. We choose five different models that represent the star at ages of the PMS phase with different rotation rates and radiative radii. At first, we run hydrodynamical 3D simulations of these models in order to equilibrate internal flows and coupling between the radiative core and the convective envelop. Then we inject a magnetic field into the fully convective model and once the MHD simulation is equilibrated we inject the resulting magnetic field into the following model. We compute all the MHD simulations by propagating the magnetic field into all HD models.

Our five MHD simulations show the mutual influence of the internal magnetic field and internal flows as the star evolves along the PMS phase. As seen in section 4, the introduction of the magnetic field in the hydrodynamical models leads to important modifications of the internal mean flows and the convective patterns. As the differential rotation profiles are quenched by the influence of the Maxwell stresses, the radial convective patterns are larger, since they are less sheared. The internal magnetic field also has an notable impact on the angular momentum transport since there are two additional contributions : the Maxwell stresses and the large scale magnetic torques. Indeed, in hydrodynamical models, the inward propagation of the angular momentum is due to viscous diffusion whereas in the MHD simulation this contribution becomes small, given the weak differential rotation present in these MHD simulations, and the inward propagation is sustained by the large scale magnetic torques and Maxwell stresses.

As the star ages along the PMS evolutionary track, we analyze the evolution of the magnetic energy. The proportion of mean field energy (TME + PME) decreases strongly from 47% to 6% of the total magnetic energy. In all models, the poloidal mean energy prevails on the toroidal mean energy. The decrease of the mean energy is coherent with results found by Gregory et al. (2012) in which the magnetic field is less axisymmetric and more complex as the radiative core is bigger.

As the magnetic field is propagated through the PMS models, we study the evolution of its topology and amplitude as well as its generation through the α -Ω effect. At first, we notice that in both zones, either convective or radiative, the magnetic energy increases as the star ages. In our five MHD simulations, we notice that the topology of the magnetic field changes strongly as the star ages. Since we follow the evolutionary path of a solar-like star, both rotation rate and aspect ratio of the convective envelop change as the star evolves along the PMS phase. The specific influence of each parameter is not always easy to disentangle in our study. The influence of internal structure was studied by Gregory et al. (2012). The results of this study show that as the radiative core becomes bigger the dipole components drop and the magnetic field becomes more and more complex. These properties are coherent with the results obtained with our MHD simulations with the ratio B axi /B non,axi dropping from 0.61 to 0.18 and the dipole field strength decreasing from 0.12 to 0.042. By plotting these two quantities as a function of the Rossby number, we want to analyze the influence of rotation on magnetic topology. Rossby numbers of our simulations are quite close to the transition value R o = 0.1 found in several studies (Pizzolato et al. 2003;[START_REF] Wright | Cambridge Workshop on Cool Stars, Stellar Systems[END_REF]Reiners et al. 2014;Schrinner et al. 2012). We notice that f dip is smaller when R o < 0.1. For small Rossby numbers, the dipole field strength increases slightly but we need more simulations to know if this value is the upper limit or if for lower Rossby number f dip continues to grow and shows that R o = 0.1 is a sharp transition. The generation of the mean magnetic field shows that as the convective zone becomes shallower and the rotation rate increases, the Ω effect becomes predominant in the generation of the mean toroidal magnetic field. Moreover the α-effect tends to generate more poloidal field than toroidal one. Hence, in each model, we see an α -Ω dynamo. Three out of our five MHD simulations display a magnetic cycle. In all cyclic cases, the time latitude diagram of the longitudinally averaged toroidal magnetic field shows a clear poleward branch starting from low latitude (see Figure 20). This magnetic field propagation is also compatible with the α -Ω dynamo concept as it satisfies Parker-Yoshimura rule (see Figure 22).

As the radiative zone grows in the star, we observe that, in all models, the magnetic field in the core, left over by the convective dynamo, is stable regarding the limit given by Braithwaite ( 2008): E pol /E tot < 0.8. The magnetic field in the radiative core of the star originates from the relaxation of the dynamo field coming from the previous stellar evolution phase in our sequence of models. By looking at this dynamo field, we notice, that in all convection zones, the magnetic field that comes from the dynamo action also fulfill the stability criteria E pol /E tot < 0.8 even if this has no obvious consequence in that zone. It seems that the relaxation of the field preserves this feature of the field and explains that all relaxed fields in radiative core are stable for the stability criteria.

The global properties of the magnetic fields we obtain in our study also have direct consequences for the coronae of PMS stars. We have computed the change of the Alfvén radius (i.e. the radius where the stellar wind decouples from the star) that such topological and ro-tation state implies, following the prescription described in Réville et al. (2015b). We find that the Alfvén radius shrinks from about 33 R to 10 R . This is in good qualitative agreement with the recent work of Réville et al. (2016) who computed realistic 3D stellar wind models along the evolutionary tracks of a solar mass star using spectro-polarimetric maps from Folsom et al. (2016). We intend in the near future to use the magnetic field coming out of our dynamo simulations to compute similar 3D wind solutions during the PMS phase. This will allow us to assess the loss of mass and angular momentum, that must also vary significantly given the large change of Alfvén radius we have identified. One can also guess that the braking due to the stellar wind would be less efficient for more complex magnetic topologies as in Brown (2014). However in the evolutionary phase studied in this article, i.e. earlier than the ZAMS, stellar contraction is the key driver of the angular momentum evolution, except, of course, during the disk-locking phase. For that earlier phase, various scenarii have been proposed (see Zanni & Ferreira 2013), and one can safely say that more compact magnetic geometry will lead to a different coupling to the disk than with a simple purely dipolar magnetic field geometry.

The 3D simulations studied in this paper, five HD progenitors and five MHD dynamo simulations, are an idealistic representation of the evolution of solar-like stars along the PMS phase. For instance, the turbulence degree in our numerical models is far from reaching that of a real star. We try to keep a comparable degree of turbulence in all our models while taking into account computational constraints (each simulation has required on average 1 Mh node hours). One could also consider a systematic parametric study of each model by varying for instance, in the simulations, their Reynolds and Prandlt numbers to better assess their sensitivity to parameter change. Still we are confident that trends presented in this work are robust as we have shown that there are in qualitative agreement with observations and compatible with previous numerical studies of stellar dynamos. In this study, we only propagate the magnetic field from one model to another. This leads to an important transient phase since this magnetic field and the global flows have to adapt to each other. One future study could be to propagate also global flows between the different models. Nevertheless, one challenge of this study is its serial aspect : before computing the model n, we need to reach an equilibrium state, in the statistically stationary sense, for the magnetohydrodynamical model n -1. Indeed, as already stated, we need the magnetic field of the simulation n -1 to initialize the model n. This serial aspect makes this study complex to compute as each model take several months to settle down. One improvement to this analysis would be to increase the number of models and make the gap between the different rotation rates and radiative radii smaller to have a smoother evolution of B and trends. This may be done by changing our way of simulating the stellar evolution.

It is important to notice that in this study, we choose our models to follow an astrophysical path along the PMS phase. A logical follow-up is therefore to apply this analysis to the evolution of solar-like stars along the following step of stellar evolution, i.e. the main sequence. In that study the main parameter will be the decrease of the rotation rate as the star is braked by the solar wind and the internal stellar structure of the star is fixed during this evolutionary phase. An additional study would be to study the impact of stellar structure with a fixed rotation rate. We have started doing such studies and their results will be reported in future communications.

Fig. 1 .

 1 Fig. 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 stellar evolution.Middle : Stellar rotation rate as a function of time. At the beginning of the PMS phase, 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.

Fig. 4 .

 4 Fig. 4.-Kinematic viscosity for the C20 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.

Fig. 5 .

 5 Fig. 5.-Radial velocity for the five HD simulations (rows) at three different depths (columns) : 96%, 60% and 20% 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.

Fig. 6 .

 6 Fig. 6.-Radial velocity for the five HD simulations timeaveraged over the last 400 days of each simulation. 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 C20 model, to 100 000 times smaller than the value observed in the convective zone.

Fig. 7 .

 7 Fig. 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 gravitoinertial waves that are excited by the downflows of the convective zone hitting the tachocline.

Fig. 8 .

 8 Fig. 8.-Contour plot-left side: Meridional circulation of HD simulations. Counter-clockwise circulation is represented in red and clockwise in blue. Contour plot-right side: Differential rotation of HD models: C00, C20, C40, C60 and C70. 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. All quantities have been time-averaged over the last 400 days of each simulation.

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 -3 . 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.

Fig

  Fig. 9.-Description of the procedure to study the PMS evolution of an stellar magnetic field. 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 the C70m model.

Fig. 10 .

 10 Fig. 10.-Evolution of energies when magnetic field resulting from the C00m model is introduced into the hydrodynamical simulation of the C20 model. The origin of time (t=0) in the figure corresponds to the introduction of the magnetic field in the C20 hydrodynamical model.

Fig. 11 .

 11 Fig. 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 Bϕ for each model.

Fig. 12

 12 Fig. 12.-3D views of the C40 and C40m models. 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.

Fig. 13 .

 13 Fig. 13.-Butterfly diagram for the longitudinal component of the magnetic field Bϕ .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). The origin of time (t=0) in the figure corresponds to the beginning of the C00m 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.

Fig. 14 .

 14 Fig. 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. All quantities have been time-averaged over the last 400 days of each simulation.

Fig. 15

 15 Fig. 15.-As in Figure7, we show here the equatorial slices of normalized radial velocity for the C20 and C20m 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.

Fig. 16 .

 16 Fig. 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, ¯ p 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 evolutionary track.

Fig. 17 .

 17 Fig. 17.-Balance between the variation of Bϕ and generation terms of the toroidal magnetic field for the C00m simulation. The first 2D plot shows the variation of Bϕ : ∆Bϕ = (B ϕ,2 -B ϕ,1). The other plots are the generation terms integrated between t 1 and t 4 (see Figure20). All the terms are reported in G/s.

Fig. 18 .

 18 Fig. 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.

Fig. 19 .

 19 Fig. 19.-Time evolution the vector potential Aϕ through one magnetic cycle of the C00m model. The three plot on the left show the vector potential at two different times t 1 and t 4 (see Figure 20), and their difference ∆Aϕ = A ϕ,2 -A ϕ,1 . 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 Bϕ : ∆Aϕ = (A ϕ,2 -A ϕ,1 ) /(t 4 -t 1 ) in G.cm/s.

Fig. 20 .

 20 Fig. 20.-Time evolution of the toroidal magnetic field at the surface of the C00m model. The origin of time (t=0) in the figure corresponds to the beginning of the C00m simulation. 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 t 1 , B pol is almost completely positive. Reversal begins at t 2 with increasing negative values in B pol and decreasing patterns amplitude for Btor. At t 3 , negative patterns of Btor mostly vanish and B pol has opposite signs in the two hemispheres. At t 4 , the reversal is achieved with a fully negative B pol and patterns of Btor have the opposite signs of those observed at t 1 . 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 rss 2.5R * (seeSchrijver et al. (2003)).

Fig. 21

 21 Fig. 21.-The amplitudes D, Q and O are shown for the cycling models : C00m, C20m and C70m. The origin of time (t=0) in the figures correspond to the beginning of each MHD simulation.We note an increase of the cycle period length. No clear phasing relationship between the three components can be established, except for the C20m case that show that symmetric and anti-symmetric dynamo mode are in opposition of phase. We also notice that in the C70m model there is no dipole reversal, even if it presents a cycling behavior.

Fig. 22 .

 22 Fig. 22.-Latitudinal component of the propagation direction of a dynamo wave S θ (see equation 40), for the C00m model. Red contours denote positive southward direction. The origin of time (t=0) in the figure corresponds to the beginning of the C00m model. Overplotted in white are isocontours of Bϕ 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.

Fig. 23 .

 23 Fig. 23.-Evolution of the ratio B pol /Btor over the PMS evolutionary track. The origin of time (t=0) in the figure corresponds to the injection the magnetic field in the C20m model.

Note:

  Relaxation time in the radiative zone. The first column gives the Alfvén speed in the core of our simulations. The second one shows the characteristic length scale of the relaxation. The Alfvén time is defined by ta = va/D RZ . The last columns shows the computational time of the MHD models in days and in Alfvén time.

TABLE 3 Characteristic numbers

 3numbers

  CZ KE RZ/CZ (10 39 erg) (10 6 erg.cm -3 ) (10 6 erg.cm -3 ) (erg.cm -3 ) (10 5 erg.cm -3 )

					TABLE 5
				Kinetic energies in hydrodynamical simulations
	Case EK RZ/C00 EK KE DRKE MCKE CKE 114 5.52 5.12 (92.7%) 379 (0.0068%) 4.05 (7.30%) -	-
	C20	62.9	6.89	6.39 (92.7%)	457 (0.0066%)	5.04 (7.30%)	4.0 × 10 -3	0.5
	C40	4.35	1.11	0.776 (70.2%)	68.2 (0.0062%)	3.29 (29.8%)	1.5 × 10 -2	0.22
	C60	4.04	3.11	2.63 (84.5%)	170 (0.0055%)	4.82 (15.5%)	1.2 × 10 -2	0.15
	C70	2.01	2.18	1.78 (81.1%)	454 (0.021%)	4.11 (18.8%)	2.7 × 10 -2	0.05

TABLE 8 α

 8 -Ω effect

			α tensor	Ω/αϕ α P /αϕ
		0.104	0.114 0.130
	C00m	0.114 0.127 0.147	4.20	12.5
		0.082	0.086 0.095
		0.133	0.129 0.115
	C20m	0.120 0.116 0.124	2.08	10.8
		0.084	0.087 0.073
		0.203	0.126 0.103
	C40m	0.149 0.097 0.098	3.42	19.1
		0.085	0.070 0.069
		0.279	0.137 0.061
	C60m	0.194 0.102 0.056	4.39	21.2
		0.077	0.053 0.040
		0.188	0.155 0.032
	C70m	0.099 0.084 0.024	1.54	9.19
		0.211	0.171 0.035

TABLE 9

 9 Magnetic energy decomposition

	Case	EM RZ	ME RZ	ME pol	RZ	MEtor RZ	ME CZ	ME pol	CZ

TABLE 10 Alfvén

 10 time for the radiative zone

	Case	va	Dn	ta	t model
		(cm.s -1 ) 10 10 cm (days) (days) (ta)
	C20m	181	2.6	1660	8062	4.87
	C40m	194	2	1190	7875	6.61
	C60m	202	1.47	838	2454	2.93
	C70m	461	0.69	1320	2555	1.94

+ ∇ • κ r ρc p ∇ T + T + κρ T ∇S + κ 0 ρ T ∇ Swhere ν and κ are effective eddy diffusivities that represent momentum and heat transport by subgrid-scale (SGS) motions, κ r is the radiative diffusivity and j =

First we wish to thank N. Featherstone for proposing the original idea that has lead to this project. We thank R.Arlt for helpful discussions. We wish also to thank A. Palacios and L. Amard for providing the 1D stellar evolution and structure models used as input to setup our 3D ASH simulations. Special thanks to K. Augustson for providing some of the analysis routines used in section 5.3. We acknowledge A. Strugarek for the 3D visualization of the C40 and C40m model. We also thank V. Réville for the analysis of the evolution of the Alfvén radius and loss of angular momentum in our MHD simulations. We acknowledge funding support by ERC STARS2 207430, CNES Solar Orbiter, Plato, INSU/PNST and ANR Toupies fundings.