Exploring the properties of RHD model atmosphere inner boundaries for AGB stars. Liljegren finds that some models fail to converge due to a temperature instability. Temperatures at the inner boundary increase rapidly before finally reaching a temperature that is outside of the pre-computed opacity grid.
The question becomes, is the temperature instability the result of missing physics near the inner boundary? We suspect that the lack of a description of convection for the inner layers of the atmosphere model are causing a false build-up of heat, leading to the runaway temperature instability. To explore whether this is the case, we'll need to first take a look at the thermodynamic state of the gas in the envelope.
RHD models adopted by Liljegren are based on the RHD model atmosphere developed by Höfner et al. (2003, A&A, 399, 589). For simplicity, and to allow for better comparison with previous studies, the RHD models adopt a perfect gas equation of state with $\gamma = 5/3$ and $\mu = 1.26$. The mean molecular weight is that for a non-interacting (i.e., perfect) gas without consideration of ionization states, which effectively equates the mean molecular weight with the mean atomic molecular weight.
Liljegren provided the following properties of the gas near the inner boundary prior to the temperature instability: $P_{\rm gas} = 1680\ {\rm Ba}$ with $T = 5600\ {\rm K}$. Using the perfect gas equation of state, one can derive a density for the gas \begin{equation} \rho = \frac{\mu m_p}{k}\frac{P_{\rm gas}}{T_{\rm gas}} \end{equation} where $k$ is the Boltzmann constant and $m_p$ is the mass of a proton. We therefore find,
In [1]:
rho = (1.26*1.6726219e-24/1.3806488e-16)*(1680./5600.)
print "Density of the gas [g/cm**3] = {:11.5e}.".format(rho)
We can also easily estimate the thermodynamic properties of the gas under the assumption of a perfect gas. We know that $\gamma = c_P / c_V = 5/3$. Of importance for our hypothesis is knowing the various thermodynamic temperature gradients, $\nabla_{\rm ad}$ and $\nabla_{\rm rad}$. The former can be expressed generally as \begin{equation} \nabla_{\rm ad} = \frac{P_{\rm gas} \delta}{\rho T_{\rm gas} c_P}, \end{equation} where $\delta = - (\partial\ln\rho / \partial\ln T_{\rm gas})_P$ (coefficient of thermal expansion) and $c_P$ is the specific heat at constant pressure. The latter can be evaluated using the general thermodynamic relation, \begin{equation} c_P - c_V = \frac{P_{\rm gas} \delta^2}{\rho T_{\rm gas} \alpha}, \end{equation} with $\alpha = (\partial\ln\rho / \partial\ln P_{\rm gas})_T$ (compressibility coefficient). Under the assumption that the gas particles are non-interacting and neglecting ionization, we have $\alpha = \delta = 1$. Therefore we can write \begin{equation} c_P = \frac{5}{2} \frac{P_{\rm gas}}{\rho T_{\rm gas}}, \end{equation} meaning \begin{equation} \nabla_{\rm ad} = \frac{2}{5} = 0.4. \end{equation}
Now, to understand whether convection may or may not be important under these conditions, we can consider a simple comparison of the adiabatic and radiative temperature gradients. Adiabatic temperature gradient is a limiting case for the temperature gradient if the gas is undergoing convection, but should provide an estimate of the super-adiabaticity of the temperature gradient.
The radiative temperature gradient is \begin{equation} \nabla_{\rm rad} = \left(\frac{d\ln T_{\rm gas}}{d\ln P_{\rm gas}}\right)_{\rm rad} = \frac{3}{16\pi acG}\frac{\kappa P_{\rm gas}}{T_{\rm gas}^4}\frac{L_r}{M_r}, \end{equation} where $X_r$ are quantities defined at a radius $r$. For our purposes, we can assume $L_r = L_{\star}$ and $M_r = M_{\star}$. This is justified as there are no energy generation sources in the outer envelope of an AGB star and the outer envelope contains a negliglbe fraction of the total mass. We can esitmate the Rosseland mean opacity from low temperature opacity sources (e.g., Ferguson et al. 2005), which yield an opacity with the range of $10^{-2}$ - $10^{-1}\ {\rm cm}^2\, {\rm g}^{-1}$. Thus,
In [2]:
def del_rad(opac, P, T, L, M):
return 7.62586e9*opac*P*L/(T**4*M)
In [3]:
P = 1680.0 # Ba
T = 5600.0 # K
M = 1.5 # Msun
L = 6953.0 # Rstar = 412 Rsun; Teff = 2600 K
print "{:11.5e} < Del_rad < {:11.5e}".format(del_rad(1.0e-2, P, T, L, M), del_rad(1.0e-1, P, T, L, M))
In the case where the radiative opacity is close to $10^{-2}\ {\rm cm}^2\, {\rm g}^{-1}$, the radiative gradient may exceed the adiabatic temperature gradient, suggesting that layers may be unstable to convective instabilities. It is therefore quite possible that thermodynamic conditions developing at the inner boundary of the RHD atmosphere models are ammenable to convective energy transport.
Since the inner boundary ignores convection and ionization, a temperature purturbation may lead to a growing instability in the temperature as radiation is becomes less able to transport the flux (at a given temperature) and ionization stages of hydrogen do not exist to alter the specific heat when the temperature begins to rapidly increase toward $8\,000 - 10\,000\ {\rm K}$.
A temperature instability is all inevitable once conditions become ammenable to convection (at least within the framework of adiabaic convection). Radiation will be unable to efficiently transport flux and the neglect of physics related to convective flux transport will force the local temperature gradient (and thus the deep interior temperature) to force the flux through the gas using only radiation.
We may be able to estimate what temperature is required, assuming the pressure is defined by the condition for hydrostatic equilibrium. The approximate conditions required for radiation to carry all of the flux through an ideal gas (neglecting ionization) is \begin{equation} \nabla_{\rm rad} = \nabla_{\rm ad} = \frac{2}{5}. \end{equation} Therefore, the local temperature can be approximated at the inner boundary \begin{equation} T_{\rm gas} = \left(\frac{15}{32\pi acG} \kappa P_{\rm gas} \frac{L_r}{M_r}\right )^{1/4}. \end{equation} For conditions at the inner boundary, the temperature must exceed
In [4]:
def T_rad(opac, P, L, M):
return (1.906465e10*opac*P*L/M)**0.25
In [5]:
print "{:11.5e} < T < {:11.5e} K".format(T_rad(1.0e-2, 1680., L, M), T_rad(1.0e-1, 1680., L, M))
which produces a density
In [6]:
rho1 = (1.26*1.6726219e-24/1.3806488e-16)*(1680./6207.)
rho2 = (1.26*1.6726219e-24/1.3806488e-16)*(1680./11038.)
print "{:11.5e} < Density < {:11.5e} [g cm**-3].".format(rho1, rho2)
Thus, the density could be forced to decrease by a factor of two compared to predictions from a non-ideal gas if the local conditions become favorable to convective instability, but the models neglect physics related to convective flux transport.
We have seen that increasing the local temperature will decrease the local density, if the layer is to maintain a constant pressure. Opacity in cool, low density environments is an increasing function of temperature and a decreasing function of density. Therefore, an increase in the local temperature and decrease in the local density will cause a dramatic increase in the opacity.