In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as scint
We need to import a file with stellar parameters and a table of stellar atmosphere properties in photospheric layers. We approximate the surface equipartition magnetic field strength as the equipartition field strength where the magnetic pressure is equal the gas pressure at a Rossland optical depth $\tau = 1$. This is roughly equivalent to the equipartition field strength in the optical photospheric layers. Stellar parameters will be taken from a non-magnetic stellar model isochrone at an age of 10 Myr.
Start by importing stellar parameters,
In [2]:
iso_10 = np.genfromtxt('../models/iso/std/dmestar_00010.0myr_z+0.00_a+0.00_phx.iso')
iso_20 = np.genfromtxt('../models/iso/std/dmestar_00020.0myr_z+0.00_a+0.00_phx.iso')
Now load a table with the atmospheric gas pressure and temperature at an optical depth $\tau = 1$ for a set of $\log(g)$s and $T_{\rm eff}$s. We'll adopt a solar metallicity to simplify things.
In [3]:
atm = np.genfromtxt('../models/atm/tab/Zp0d0.ap0d0_t001.dat')
Values for the effective temperature of the atmosphere model is tabulated in column 0, but we must define an array with $\log(g)$ values.
In [4]:
teffs = np.transpose(atm[:, 0])
loggs = np.arange(-0.5, 5.6, 0.5)
print teffs
print loggs
The atmosphere structure table isn't quite in the correct form for our purposes, as pressures and temperatures are intermingled. We should separate those properties into individual pressure and temperature tables.
In [5]:
temps = np.empty((len(teffs), len(loggs)))
press = np.empty((len(teffs), len(loggs)))
for i, teff in enumerate(atm[:, 1:]):
for j, prop in enumerate(teff):
if prop == 0.:
prop = np.nan
else:
pass
if j%2 == 0:
press[i, j/2] = prop
else:
temps[i, j/2] = prop
With the individual tables formed, we now need to construct interpolation surfaces using a 2D interpolation routine. Note that we only really care about the pressure table, as that sets the equipartition magnetic field strengths.
In [6]:
pres_surface = scint.interp2d(teffs, loggs, np.transpose(press), kind='linear')
We are in a position to compute surface pressures and, by extension, equipartition magnetic field strengths.
In [7]:
B_eq_10 = np.empty((len(iso_10[:62])))
B_eq_20 = np.empty((len(iso_20[:62])))
for i, star in enumerate(iso_10[:62]):
B_eq_10[i] = np.sqrt(8.0*np.pi*pres_surface(10**star[1], star[2]))
for i, star in enumerate(iso_20[:62]):
B_eq_20[i] = np.sqrt(8.0*np.pi*pres_surface(10**star[1], star[2]))
See what kind of values we obtain.
In [8]:
B_eq_20
Out[8]:
These values match estimates from convective energy equiparition and observational measurements (e.g., Saar). The general trend is that equipartition field strengths decrease toward higher masses/temperatures, which matches intuition since surface gas pressures decrease as stellar surface layers become more extended and fluffy.
We can visualize this:
In [9]:
fig, ax = plt.subplots(1, 1, figsize=(9., 6.))
ax.set_xlabel('${\\rm Mass}\ (M_{\\odot})$', fontsize=20)
ax.set_ylabel('$\\langle {\\rm B}f \\rangle_{\\rm eq}$', fontsize=20)
ax.grid()
ax.plot(iso_10[:62, 0], B_eq_10, 'o-', lw=1, color='#333333', alpha=0.7)
ax.plot(iso_20[:62, 0], B_eq_20, 'o-', lw=1, color='#1e90ff', alpha=0.7)
Out[9]:
and as a function of effective temperature:
In [10]:
fig, ax = plt.subplots(1, 1, figsize=(9., 6.))
ax.set_xlabel('$T_{\\rm eff}\ (K)$', fontsize=20)
ax.set_ylabel('$\\langle {\\rm B}f \\rangle_{\\rm eq} ({\\rm G})$', fontsize=20)
ax.grid()
ax.plot(10**iso_10[:62, 1], B_eq_10, 'o-', lw=1, color='#333333', alpha=0.7)
ax.plot(10**iso_20[:62, 1], B_eq_20, 'o-', lw=1, color='#1e90ff', alpha=0.7)
Out[10]:
For the manuscript, there is a table of stellar properties and the resulting equipartition surface magnetic field strengths. Here is that table,
In [11]:
for i, star in enumerate(iso_10[:62]):
if int(star[0]*100.)%10 == 0:
print '{:5.1f} & {:6.0f} & {:6.2f} & {:6.2f} \\\\'.format(star[0], 10**star[1], star[2], B_eq_10[i]/1.0e3)
else:
pass
When running models, an interpolation curve is constructed from the data in the table to find the equipartition value as a function of mass. Let's construct that curve.
In [12]:
B_eq_curve_data = np.empty((17, 2))
j = 0
for i, star in enumerate(iso_10[:62]):
if int(star[0]*100.)%10 == 0:
B_eq_curve_data[j, 0] = star[0]
B_eq_curve_data[j, 1] = B_eq_10[i]
j += 1
else:
pass
B_eq_interp_curve = scint.interp1d(B_eq_curve_data[:, 0], B_eq_curve_data[:, 1], kind='cubic', axis=0)
Now we can compare the accuracy of the interpolated curve.
In [13]:
fig, ax = plt.subplots(1, 1, figsize=(9., 6.))
ax.set_xlabel('${\\rm Mass}\ (M_{\\odot})$', fontsize=20)
ax.set_ylabel('$\\langle {\\rm B}f \\rangle_{\\rm eq}\ ({\\rm G})$', fontsize=20)
ax.grid()
ax.plot(iso_10[:62, 0], B_eq_10, 'o', lw=1, color='#333333', alpha=0.7)
ax.plot(np.arange(0.1, 1.7, 0.02), B_eq_interp_curve(np.arange(0.1, 1.7, 0.02)), '-', lw=3, color='#1e90ff')
Out[13]:
Surface magnetic field strengths are a over-estimated just prior to the ZAMS. Otherwise, the interpolated curve looks to reproduce the origianal data to a high degree.