Magnetic models of young stars have their peak magnetic field strength prescribed where $R = 0.5 R_{\star}$. This works well and permits models of young stars with strong surface magnetic fields ($\langle{\rm B}f\rangle$) to converge and evolve without issue. We can understand this by plotting how the peak magnetic field strength varies as a function of depth in the dipole magnetic field models. The magnetic field profile in this case is prescribed to be
\begin{equation} B(R) = B_{\rm surf} \left(\frac{R_{\star}}{R}\right)^3. \end{equation}We can therefore easily visualize the peak magnetic field strength as a function of the surface magnetic field strength and the cut-off radius where we define the peak strength.
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
In [2]:
radial_points = np.arange(0.01, 1.0, 0.01) # units of Rstar
bfield_scaling = radial_points**(-3.0) # see equation (1)
bfield_surface = np.arange(0.5, 4.1, 0.5) # units of kiloGauss
which leads to the following figure,
In [3]:
fig, ax = plt.subplots(1, 1, figsize=(8.0, 8.0))
# configure axis labels and limits
ax.set_xlabel('Radius Fraction, $R/R_{\\star}$', fontsize=20.)
ax.set_ylabel('Magnetic Field Strength (kG)', fontsize=20.)
ax.set_ylim(0.1, 1.0e3)
ax.tick_params(which='major', axis='both', length=10.0, labelsize=16.)
ax.grid(True)
for bf in bfield_surface:
ax.semilogy(radial_points, bf*bfield_scaling, '-', lw=3)
fig.tight_layout()
It is clear that interior magnetic field strengths remain below 100 kG for all values of the surface field strength when the peak magnetic field strength is assigned to a radial fraction greater than approximately $0.40 R_{\star}$. These values for magnetic field strengths ($< 100 \textrm{ kG}$) permit convergence as they produce values of the magnetic pressure and magnetic energy density that are well below values for the gas pressure and the internal energy of the ambient plasma.
This is important to understand once we introduce magnetic fields in young stars that later develop a radiative core during their pre-main-sequence evolution. Since the code was initially designed to treat the effects of magnetic fields in main-sequence stars, the development of a radiative core was not explicitly considered. Instead, radiative cores existed prior to the inclusion of a magnetic field. Since maost stars have radiative cores that extend no lower than $0.40 R_{\star}$, there was no need to be concerned about convergence outside of strong surface magnetic fields. As a result, we defined the peak magnetic field strength to occur at the base of the tachocline (radiative-convection zone interface).
Fully convective stars were assigned peak magnetic field strengths at $0.30 R_{\star}$, which is just beyond the maximum depth observed for the convection zone in partially convective stars and is also where the peak magnetic field strength appeared (a very slight peak) in 3D MHD models of fully convective dynamos (Browning et al. 2008).
Pre-main-sequence stars therefore pose an interesting problem: star that eventually become partially convective on the main sequence start off as fully convective stars along the Hyashi track. However, once a radiative core begins developing at the beginning of the Henyey track, the peak magnetic field strength jumps from $0.50 R_{\star}$ to $0.12 R_{\star}$, which causes an immense increase in the peak magnetic field strength. Such a deep convection zone is not seen in main-sequence low-mass stars (they'd become fully convective, instead). We therefore need to control for the evolutionary development of a radiative core and the receding convection zone boundary.
Possible solutions:
Solution 1 is now implemented in an experimental version that was quite easy to implement. Within the magnetic field module, the tachocline definition function was adjusted to select between the peak magnetic field strength location in a fully convection model (fc_tach) and the radius of the tachocline (r_tach) depending on which is larger. That is
r_tach = max(fc_tach, r_tach)
Test models with masses between $0.40 M_{\odot}$ and $0.65 M_{\odot}$ are now running. Models at the lower mass end of this range can be compared with previous converged models.
In [14]:
filename = 'm0400_GAS07_p000_p0_y26_mlt2.040_mag25kG.trk'
std_iso = np.genfromtxt('../../evolve/dmestar/trk/gas07/p000/a0/amlt2202/m0400_GAS07_p000_p0_y26_mlt2.202.trk')
old_iso = np.genfromtxt('../../evolve/dmestar/trk/gas07/p000/a0/amlt2040/mag25kG/{0}'.format(filename))
new_iso = np.genfromtxt('../../evolve/models/tmp/{0}'.format(filename))
In [18]:
fig , ax = plt.subplots(3, 1, figsize=(8, 12), sharex=True)
ax[0].set_title('$M = 0.40 M_{\\odot}$ with $\\langle{\\rm B}f\\rangle$ = 2.5 kG', family='serif', fontsize=20.)
ax[2].set_xlabel('Age (Myr)', fontsize=20.)
ax[0].set_ylabel('Effective Temperature (K)', fontsize=20.)
ax[1].set_ylabel('Luminosity ($L_{\\odot}$)', fontsize=20.)
ax[2].set_ylabel('Radius ($R_{\\odot}$)', fontsize=20.)
for axis in ax:
axis.set_xlim(1.0e-1, 1.0e1)
axis.tick_params(which='major', axis='both', length=10.0, labelsize=16.)
# Temperature
ax[0].semilogx(old_iso[:,0]/1.0e6, 10**old_iso[:,1], '-', lw=3, color='#800000')
ax[0].semilogx(new_iso[:,0]/1.0e6, 10**new_iso[:,1], '-', lw=2, color='#1e90ff')
ax[0].semilogx(std_iso[:,0]/1.0e6, 10**std_iso[:,1], '-', lw=4, color='#555555')
# Luminosity
ax[1].semilogx(old_iso[:,0]/1.0e6, 10**old_iso[:,3], '-', lw=3, color='#800000')
ax[1].semilogx(new_iso[:,0]/1.0e6, 10**new_iso[:,3], '-', lw=2, color='#1e90ff')
ax[1].semilogx(std_iso[:,0]/1.0e6, 10**std_iso[:,3], '-', lw=4, color='#555555')
# Radius
ax[2].semilogx(old_iso[:,0]/1.0e6, 10**old_iso[:,4], '-', lw=3, color='#800000')
ax[2].semilogx(new_iso[:,0]/1.0e6, 10**new_iso[:,4], '-', lw=2, color='#1e90ff')
ax[2].semilogx(std_iso[:,0]/1.0e6, 10**std_iso[:,4], '-', lw=4, color='#555555')
Out[18]:
In the figure above, the standard track (dark grey) is slightly cooler than the unperturbed magnetic tracks due to a difference in the depth at which the surface boundary conditions are defined. The above track is used only for illustration, as $\tau_{\rm ross} = 10$ standard models have not been computed, yet.
Now, at a higher mass where magnetic models have not previously been computed.
In [19]:
filename = 'm0600_GAS07_p000_p0_y26_mlt2.040_mag25kG.trk'
std_iso = np.genfromtxt('../../evolve/dmestar/trk/gas07/p000/a0/amlt2202/m0600_GAS07_p000_p0_y26_mlt2.202.trk')
new_iso = np.genfromtxt('../../evolve/models/tmp/{0}'.format(filename))
In [22]:
fig , ax = plt.subplots(3, 1, figsize=(8, 12), sharex=True)
ax[0].set_title('$M = 0.60 M_{\\odot}$ with $\\langle{\\rm B}f\\rangle$ = 2.5 kG', family='serif', fontsize=20.)
ax[2].set_xlabel('Age (Myr)', fontsize=20.)
ax[0].set_ylabel('Effective Temperature (K)', fontsize=20.)
ax[1].set_ylabel('Luminosity ($L_{\\odot}$)', fontsize=20.)
ax[2].set_ylabel('Radius ($R_{\\odot}$)', fontsize=20.)
ax[1].set_ylim(0.0, 5.0)
ax[2].set_ylim(0.0, 5.0)
for axis in ax:
axis.set_xlim(1.0e-1, 1.0e1)
axis.tick_params(which='major', axis='both', length=10.0, labelsize=16.)
# Temperature
ax[0].semilogx(new_iso[:,0]/1.0e6, 10**new_iso[:,1], '-', lw=2, color='#1e90ff')
ax[0].semilogx(std_iso[:,0]/1.0e6, 10**std_iso[:,1], '-', lw=4, color='#555555')
# Luminosity
ax[1].semilogx(new_iso[:,0]/1.0e6, 10**new_iso[:,3], '-', lw=2, color='#1e90ff')
ax[1].semilogx(std_iso[:,0]/1.0e6, 10**std_iso[:,3], '-', lw=4, color='#555555')
# Radius
ax[2].semilogx(new_iso[:,0]/1.0e6, 10**new_iso[:,4], '-', lw=2, color='#1e90ff')
ax[2].semilogx(std_iso[:,0]/1.0e6, 10**std_iso[:,4], '-', lw=4, color='#555555')
Out[22]:
Differences in the effeictive temperature frame are considerable; nearly 1000 K cooler in the case where magnetic fields are implemented! We can now push the models further with the new boundary definition. For example,
In [23]:
filename = 'm0900_GAS07_p000_p0_y26_mlt2.040_mag25kG.trk'
std_iso = np.genfromtxt('../../evolve/dmestar/trk/gas07/p000/a0/amlt2202/m0900_GAS07_p000_p0_y26_mlt2.202.trk')
new_iso = np.genfromtxt('../../evolve/models/tmp/{0}'.format(filename))
In [25]:
fig , ax = plt.subplots(3, 1, figsize=(8, 12), sharex=True)
ax[0].set_title('$M = 0.90 M_{\\odot}$ with $\\langle{\\rm B}f\\rangle$ = 2.5 kG', family='serif', fontsize=20.)
ax[2].set_xlabel('Age (Myr)', fontsize=20.)
ax[0].set_ylabel('Effective Temperature (K)', fontsize=20.)
ax[1].set_ylabel('Luminosity ($L_{\\odot}$)', fontsize=20.)
ax[2].set_ylabel('Radius ($R_{\\odot}$)', fontsize=20.)
ax[1].set_ylim(0.0, 7.0)
ax[2].set_ylim(1.0, 7.0)
for axis in ax:
axis.set_xlim(1.0e-1, 1.0e1)
axis.tick_params(which='major', axis='both', length=10.0, labelsize=16.)
# Temperature
ax[0].semilogx(new_iso[:,0]/1.0e6, 10**new_iso[:,1], '-', lw=2, color='#1e90ff')
ax[0].semilogx(std_iso[:,0]/1.0e6, 10**std_iso[:,1], '-', lw=4, color='#555555')
# Luminosity
ax[1].semilogx(new_iso[:,0]/1.0e6, 10**new_iso[:,3], '-', lw=2, color='#1e90ff')
ax[1].semilogx(std_iso[:,0]/1.0e6, 10**std_iso[:,3], '-', lw=4, color='#555555')
# Radius
ax[2].semilogx(new_iso[:,0]/1.0e6, 10**new_iso[:,4], '-', lw=2, color='#1e90ff')
ax[2].semilogx(std_iso[:,0]/1.0e6, 10**std_iso[:,4], '-', lw=4, color='#555555')
Out[25]:
With convergence issues largely mitigated with the revised definition of the peak magnetic field strength accomodating the transition from fully convective to partially convective interiors, I'm now computing a small grid of magnetic models for the generation of new isochrones.