Preparation for AGB Group Meeting

Sofia requested that each participant bring a plot (of something!) for discussion at the AGB group meeting on 27 Aug 2015. A few ideas about which plot(s) to bring: (1) demonstration that more recent solar abundance distributions (Asplund et al. 2009; Caffau et al. 2011) are in better agreement with NIR broadband photometric colors and should be implemented in AGB atmosphere models. Lower oxygen abundance is much better for M dwarfs owing to decreased H20, whereas M giant wind driving mechanisms may be inhibited due to decreased oxygen available to capture into grains? (2) Magnetic fields have been tentatively detected at the surface of AGB stars with typical magnetic field strength estimates in the 1 – 5 G range (Ramstedt, priv. comm.). What magnetic field strengths are required to interact and inhibit convection in AGB star envelopes (neglecting large-scale dynamics)?


In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

Solar Abundances and NIR Broadband Photometric Colors

We'll begin by loading necessary data, here NIR photometric data for the Pleiades (and M67?).


In [36]:
Pleiades_NIR = np.genfromtxt('../Projects/pleiades_colors/data/Stauffer_Pleiades_nir.dat', 
                             usecols=(2, 3, 5, 6, 8, 9))  # J, errJ, H, errH, K, errK

dist_mod = 5.62 # distance modulus for the Pleiades

Now we'll need to load appropriate stellar evolution isochrones with bolometric corrections derived from MARCS model atmospheres. MARCS uses Grevesse et al. (2007), but differences between colors derived from GAS07, AGSS09, and CIFIST11 are negligible.


In [5]:
iso_gs98  = np.genfromtxt('../Projects/pleiades_colors/data/dmestar_00120.0myr_z+0.00_a+0.00_phx.iso')
iso_gas07 = np.genfromtxt('../Projects/pleiades_colors/data/dmestar_00120.0myr_z+0.00_a+0.00_marcs.iso')

We've loaded isochrones with both Grevesse & Sauval (1998) and Grevesse et al. (2007) solar abundance distributions. The former is similar to Anders & Grevesse (1989) and Grevesse & Noels (1993) distributions in that the solar surface $(Z/X) = 0.023$ – $0.025$. On the other hand, the Grevesse et al. (2007) distribution has a characteristic $(Z/X) = 0.014$.

Starting with CMDs against $M_K$.


In [33]:
fig, ax = plt.subplots(1, 3, figsize=(15., 8.), sharey=True)

ax[0].set_ylabel('$M_K$', fontsize=22.)
for axis in ax:
    axis.grid(True)
    axis.set_ylim(9., 3.)
    axis.tick_params(which='major', axis='both', length=15., labelsize=16.)

# (J-K)
ax[0].set_xlim(0.2, 1.2)
ax[0].set_xlabel('$(J - K)$', fontsize=22.)
ax[0].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,4] - 0.02, Pleiades_NIR[:,4] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[0].plot(iso_gs98[:,10] - iso_gs98[:,12], iso_gs98[:,12], dashes=(20., 5.), lw=3, c='#B22222')
ax[0].plot(iso_gas07[:,11] - iso_gas07[:,13], iso_gas07[:,13], lw=3, c='#0094b2')

# (J-H)
ax[1].set_xlim(0.0, 1.0)
ax[1].set_xlabel('$(J - H)$', fontsize=22.)
ax[1].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,2] - 0.01, Pleiades_NIR[:,4] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[1].plot(iso_gs98[:,10] - iso_gs98[:,11], iso_gs98[:,12], dashes=(20., 5.), lw=3, c='#B22222')
ax[1].plot(iso_gas07[:,11] - iso_gas07[:,12], iso_gas07[:,13], lw=3, c='#0094b2')

# (H-K)
ax[2].set_xlim(0.0, 0.5)
ax[2].set_xlabel('$(H - K)$', fontsize=22.)
ax[2].plot(Pleiades_NIR[:,2] - Pleiades_NIR[:,4] - 0.01, Pleiades_NIR[:,4] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[2].plot(iso_gs98[:,11] - iso_gs98[:,12], iso_gs98[:,12], dashes=(20., 5.), lw=3, c='#B22222')
ax[2].plot(iso_gas07[:,12] - iso_gas07[:,13], iso_gas07[:,13], lw=3, c='#0094b2')


Out[33]:
[<matplotlib.lines.Line2D at 0x10efec750>]

Now CMDs against $M_J$


In [34]:
fig, ax = plt.subplots(1, 3, figsize=(15., 8.), sharey=True)

ax[0].set_ylabel('$M_J$', fontsize=22.)
for axis in ax:
    axis.grid(True)
    axis.set_ylim(10., 3.)
    axis.tick_params(which='major', axis='both', length=15., labelsize=16.)

# (J-K)
ax[0].set_xlim(0.2, 1.2)
ax[0].set_xlabel('$(J - K)$', fontsize=22.)
ax[0].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,4] - 0.02, Pleiades_NIR[:,0] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[0].plot(iso_gs98[:,10] - iso_gs98[:,12], iso_gs98[:,10], dashes=(20., 5.), lw=3, c='#B22222')
ax[0].plot(iso_gas07[:,11] - iso_gas07[:,13], iso_gas07[:,11], lw=3, c='#0094b2')

# (J-H)
ax[1].set_xlim(0.0, 1.0)
ax[1].set_xlabel('$(J - H)$', fontsize=22.)
ax[1].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,2] - 0.01, Pleiades_NIR[:,0] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[1].plot(iso_gs98[:,10] - iso_gs98[:,11], iso_gs98[:,10], dashes=(20., 5.), lw=3, c='#B22222')
ax[1].plot(iso_gas07[:,11] - iso_gas07[:,12], iso_gas07[:,11], lw=3, c='#0094b2')

# (H-K)
ax[2].set_xlim(0.0, 0.5)
ax[2].set_xlabel('$(H - K)$', fontsize=22.)
ax[2].plot(Pleiades_NIR[:,2] - Pleiades_NIR[:,4] - 0.01, Pleiades_NIR[:,0] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[2].plot(iso_gs98[:,11] - iso_gs98[:,12], iso_gs98[:,10], dashes=(20., 5.), lw=3, c='#B22222')
ax[2].plot(iso_gas07[:,12] - iso_gas07[:,13], iso_gas07[:,11], lw=3, c='#0094b2')


Out[34]:
[<matplotlib.lines.Line2D at 0x10f4aa610>]

And finally, against $M_H$


In [35]:
fig, ax = plt.subplots(1, 3, figsize=(15., 8.), sharey=True)

ax[0].set_ylabel('$M_H$', fontsize=22.)
for axis in ax:
    axis.grid(True)
    axis.set_ylim(10., 3.)
    axis.tick_params(which='major', axis='both', length=15., labelsize=16.)

# (J-K)
ax[0].set_xlim(0.2, 1.2)
ax[0].set_xlabel('$(J - K)$', fontsize=22.)
ax[0].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,4] - 0.02, Pleiades_NIR[:,2] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[0].plot(iso_gs98[:,10] - iso_gs98[:,12], iso_gs98[:,11], dashes=(20., 5.), lw=3, c='#B22222')
ax[0].plot(iso_gas07[:,11] - iso_gas07[:,13], iso_gas07[:,12], lw=3, c='#0094b2')

# (J-H)
ax[1].set_xlim(0.0, 1.0)
ax[1].set_xlabel('$(J - H)$', fontsize=22.)
ax[1].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,2] - 0.01, Pleiades_NIR[:,2] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[1].plot(iso_gs98[:,10] - iso_gs98[:,11], iso_gs98[:,11], dashes=(20., 5.), lw=3, c='#B22222')
ax[1].plot(iso_gas07[:,11] - iso_gas07[:,12], iso_gas07[:,12], lw=3, c='#0094b2')

# (H-K)
ax[2].set_xlim(0.0, 0.5)
ax[2].set_xlabel('$(H - K)$', fontsize=22.)
ax[2].plot(Pleiades_NIR[:,2] - Pleiades_NIR[:,4] - 0.01, Pleiades_NIR[:,2] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[2].plot(iso_gs98[:,11] - iso_gs98[:,12], iso_gs98[:,11], dashes=(20., 5.), lw=3, c='#B22222')
ax[2].plot(iso_gas07[:,12] - iso_gas07[:,13], iso_gas07[:,12], lw=3, c='#0094b2')


Out[35]:
[<matplotlib.lines.Line2D at 0x10f941a90>]

Spot Temperatures on AGB Stars

Using fairly simple energy equipartition arguments with a polytropic equation of state, I was able to estimate spot temperature contrasts for low mass stars with convective envelopes (G – M) with reasonble accuracy. Can we do this for AGB stars?

We'll use MARCS model hydrostatic model atmospheres to derive approximations for atmospheric properties where $\tau \approx 1$.


In [47]:
# log(g) = -0.5; [m/H] = 0.0
Teff = np.array([2500., 2700., 2800.])
Ttau = np.array([2850., 3000., 3200.])
Peff = np.array([5.0e1, 1.3e2, 1.7e2])
Ptau = np.array([3.8e2, 4.7e2, 5.5e2])
Rtau = np.array([1.41e-9, 2.33e-9, 2.61e-9])

gammas = (np.log(Ptau) - np.log(1.0e13))/np.log(Rtau)


[ 0.87109289  0.86036103  0.86095866]

Energy equipartition dictates that the temperature ratio of a spot to its surroundings is roughly equal to $0.4^{(\gamma - 1)/\gamma}$. This leads to the following spot temperatures,


In [48]:
fig, ax = plt.subplots(1, 1, figsize=(8,4))

ax.set_xlabel('Effective Temperature (K)', fontsize=20.)
ax.set_ylabel('Spot Temperature (K)', fontsize=20.)
ax.tick_params(which='major', axis='both', length=10., labelsize=16.)
ax.grid(True)

ax.plot(Teff, Ttau*0.4**((gammas - 1)/gammas), '--o', c="#b22222", markersize=10.)


Out[48]:
[<matplotlib.lines.Line2D at 0x1093abbd0>]

Temperature contrasts are on the order of 86 – 87%, meaning spots are only slightly cooler than the surrounding photosphere.

Estimating magnetic field strengths required to generate spots is more difficult to assess. One could compare the local Alfven velocity to the convective velocity. Alternatively, following an analogy with the Sun, umbral magnetic field strengths often reach equipartition with the gas pressure outside of a spot. Therefore,


In [50]:
fig, ax = plt.subplots(1, 1, figsize=(8,4))

ax.set_xlabel('Effective Temperature (K)', fontsize=20.)
ax.set_ylabel('Equipartition Strength (G)', fontsize=20.)
ax.tick_params(which='major', axis='both', length=10., labelsize=16.)
ax.grid(True)

ax.plot(Teff, np.sqrt(8.0*np.pi*Ptau), '--o', c="#b22222", markersize=10.)


Out[50]:
[<matplotlib.lines.Line2D at 0x11047d390>]

These values of the magnetic field strength are two orders of magnitude larger than observed magnetic field strengths. Are they dynamically important? At the risk of over-relying on simple theories, I'll adopt the convective velocities from mixing length theory.


In [54]:
# arbitrary velocity values; rho ~ 10**-9
velocities = np.arange(1.0, 5.0, 0.5) # log-scale

# B-fields to get V_A = 0.1 V_C
B_010 = np.sqrt(4.0e-9*np.pi*(0.1*10**velocities)**2)
B_100 = np.sqrt(4.0e-9*np.pi*(10**velocities)**2)

print B_100


[  1.12099824e-03   3.54490770e-03   1.12099824e-02   3.54490770e-02
   1.12099824e-01   3.54490770e-01   1.12099824e+00   3.54490770e+00]

From these magnetic field strengths, we see that magnetic fields may be dynamically important once they reach a strength of order 1 – 10 G, which corresponds to the range observed among AGB stars.