Time-dependent radiative transfer with PHOENIX

D. Jack; P. H. Hauschildt; E. Baron

doi:10.1051/0004-6361/200810982

Home

All issues

Volume 502 / No 3 (August II 2009)

A&A, 502 3 (2009) 1043-1049

Full HTML

Free Access

This article has an erratum: [erratum]

Issue		A&A Volume 502, Number 3, August II 2009


Page(s)		1043 - 1049
Section		Numerical methods and codes
DOI		https://doi.org/10.1051/0004-6361/200810982
Published online		15 June 2009

Time-dependent radiative transfer with `PHOENIX`

D. Jack¹ - P. H. Hauschildt¹ - E. Baron^1,2

1 - Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany
2 - Homer L. Dodge Department of Physics and Astronomy, University of Oklahoma, 440 W Brooks, Rm 100, Norman, OK 73019-2061, USA

Received 18 September 2008 / Accepted 3 April 2009

Abstract
Aims. We present first results and tests of a time-dependent extension to the general purpose model atmosphere code PHOENIX. We aim to produce light curves and spectra of hydro models for all types of supernovae.
Methods. We extend our model atmosphere code PHOENIX to solve time-dependent non-grey, NLTE, radiative transfer in a special relativistic framework. A simple hydrodynamics solver was implemented to keep track of the energy conservation of the atmosphere during free expansion.
Results. The correct operation of the new additions to PHOENIX were verified in test calculations.
Conclusions. We have shown the correct operation of our extension to time-dependent radiative transfer and will be able to calculate supernova light curves and spectra in future work.

Key words: stars: supernovae: general - radiative transfer - methods: numerical

1 Introduction

All types of supernovae are important for the role that they play in understanding stellar evolution and galactic nucleosynthesis and as cosmological probes. Type Ia supernovae are of particular cosmological interest, e.g., because the dark energy was discovered with type Ia supernovae (Perlmutter et al. 1999; Riess et al. 1998).

In dark energy studies, the goal now is to characterize the nature of the dark energy as a function of redshift. While there are other probes that will be used (gravitational lensing, baryon acoustic oscillations), a JDEM or Euclid mission will likely consider supernovae in some form. In planning for future dark energy studies both from space and from the ground, it is important to know whether the mission will require spectroscopy of modest resolution, or whether pure imaging or grism spectroscopy will be adequate. Several purely spectral indicators of peak luminosity have been proposed (Bongard et al. 2006; Nugent et al. 1995; Bronder et al. 2008; Le Du et al. 2009; Hachinger et al. 2006; Foley et al. 2008). What is required is an empirical and theoretical comparison of both light curve shape luminosity indicators (Pskovskii 1977; Phillips 1993; Goldhaber et al. 2001; Riess et al. 1996) and spectral indicators.

To make this comparison one needs to know more about the physics going on in a supernova explosion and to be able to calculate light curves and spectra self-consistently. Thus, we need to extend our code to time-dependent problems. While our primary focus is on type Ia supernovae, the time-dependent radiative transfer code is applicable to all types of supernovae, as well as to other objects, e.g., stellar pulsations.

In the following we present the methods we used to implement time dependence. First we focus on solving the energy equation (first law) and then turn our attention to time-dependent radiative transfer.

2 Equation of energy conservation

To compute light curves we extended PHOENIX with a time-dependent solution of the non-grey, NLTE radiative transfer problem in special relativistic environments and a simple hydrodynamic code to solve the free expanding case relevant to supernova light curves. The idea is to keep track of the conservation of energy for the gas and radiation together and to allow for different time scales for the gas and the radiation.

For the time-dependent radiative transfer problem we extended our existing radiative transfer code PHOENIX (Hauschildt & Baron 2004b,1999). The change in the energy density of a radiating material is given by Eq. (96.15) in Mihalas & Mihalas (1984)

$\begin{displaymath}\frac{D}{Dt}E=-\frac{\partial}{\partial M_{r}}\left(L_{r}+P_{r}\right)+\varepsilon, \end{displaymath}$

(1)

where E is the total energy density. All quantities are considered in the comoving frame. P_r is not the pressure, but rather mechanical power on the sphere of a radius r. Equation (1) is only valid to first order in v/c, and thus lacks the full special relativistic accuracy of PHOENIX. This is adequate for the velocities found in supernovae. The total energy density of a radiating fluid consists of the sum of the energy density of the material, the energy density of the radiation field, the kinetic energy density of the material, and the gravitational energy density:

$\begin{displaymath}E=E_{\rm gas}+\frac{E_{0}}{\rho}+E_{\rm kin}+E_{\rm grav}. \end{displaymath}$

(2)

For supernovae in the free expansion phase, the gravitational energy density $E_{\rm grav}$ is negligible since the potential is small in absolute value with the standard choice of zero at infinity. Homologous expansion is a reasonably good assumption for supernovae. The energy release by the decay of ⁵⁶Ni can influence the dynamics of the expansion (Pinto & Eastman 2000). Woosley et al. (2007) compared a study following this energy release to the results from assuming homologous expansion. Figure 2 in their paper shows the deviation and density variations can be as large as 10%. However, this is probably an upper limit due to the simple burning parameterization used in that study. Ultimately, when the deflagration to detonation transition is understood, it will be important to revisit this issue and replace 1-D calculations with full 3-D calculations, which include the effects of clumps, as well as nickel bubble expansion. For now the accuracy of homologous expansion should be adequate, given the other uncertainties in the problem.

With the assumption of homology, the velocity of a given matter element is then constant as is the kinetic energy density. Thus, we can neglect the kinetic energy term $\frac{D E_{\rm kin}}{Dt}$ . So for our approach, we only have to consider the energy densities of the radiation field and the material. For the material, this includes effects such as an energy deposition due to radioactive decay of nickel and cobalt in an SN Ia envelope.

The other possible cause of a change in the energy density is the structure term. This term is given by (Cooperstein et al. 1986)

$\begin{displaymath}\frac{\partial}{\partial M_{r}}\left(P_{r}+L_{r}\right)=\frac... ...ft\{4\pi r^{2}\left[u\left(p+P_{0}\right)+F_{0}\right]\right\} \end{displaymath}$

(3)

where p is the pressure of the material and P₀the radiation pressure, u the velocity of the expanding gas, the radiative flux is represented by F₀, and the mass inside of the radius r of a layer is given by M_r. The radiation pressure is a result of the solution of the detailed radiative transfer equation and given by

$\begin{displaymath}P_{0}=\frac{4\pi}{c}K, \end{displaymath}$

(4)

with K the second moment of the radiation field. The change of the energy density is given by the two quantities

$\begin{displaymath}L_{r}=4\pi r^{2}F_{0} \end{displaymath}$

(5)

and

$\begin{displaymath}P_{r}=4\pi r^{2}u\left(p+P_{0}\right). \end{displaymath}$

(6)

If the atmosphere is in radiative equilibrium, the structure term is zero and the energy density stays constant if there is no additional energy source and the atmosphere is not expanding.

All the quantities required for the structure term can be derived from thermodynamics or the solution of the radiative transfer problem. We need the energy density of the material $E_{\rm gas}$ and the energy density of the radiation field $\frac{E_{0}}{\rho}$ . For the latter, we have to solve the radiative transfer equation for the radiation field and the radiative energy. We use our radiative transfer code PHOENIX to solve the time-independent radiative transfer equation. The energy of the radiation field is given by

$\begin{displaymath}E_{0}=\frac{4\pi}{c}J, \end{displaymath}$

(7)

where J is the mean intensity and c the speed of light.

The energy density of the material is given by

$\begin{displaymath}E_{\rm gas}=\frac{3}{2}\frac{p}{\rho}=\frac{3}{2}\frac{R}{m_{\rm u}}T, \end{displaymath}$

(8)

with the mean molecular weigh $m_{\rm u}$ and the universal gas constant R. The gas pressure is represented by p and the density by $\rho$ . T stands for the temperature of the gas. The sum of the radiation and material energy density is then the total energy density

$\begin{displaymath}E_{\rm total}=E_{\rm gas}+\frac{E_{0}}{\rho}\cdot \end{displaymath}$

(9)

The change in this total energy density is given by Eq. (1). So the equation to calculate the new energy density $E_{\rm new}$ is given by

$\begin{displaymath}E_{\rm new}=E_{\rm old} -\frac{\partial}{\partial M_{r}}\left(L_{r}+P_{r}\right)\Delta t+\varepsilon\Delta t. \end{displaymath}$

(10)

We now have all the needed equations to calculate a simple light curve. One problem for the calculation is that we can only determine the change in the total energy density for the next time step. However, the total energy change is divided into a change in the gas energy density and the energy density of the radiation field. To obtain the correct distribution of the gas and radiative energy, one has to iterate for each time step by solving the radiative transfer equation to compute the correct new temperature at the next time step.

To get the correct new temperature we apply the following iteration scheme. The error in the energy density, $E_{\rm err}$ is given by

$\begin{displaymath}E_{\rm err}=\frac{E_{\rm current}-E_{\rm target}}{E_{\rm target}}\cdot \end{displaymath}$

(11)

Here, $E_{\rm target}$ is the desired new total energy density, which is known from Eq. (10), and $E_{\rm current}$ is the total energy density obtained by Eq. (9) with the current temperature guess and the radiative transfer solution. Tests have shown that the error is almost linearly proportional to the temperature T. Therefore, a new temperature guess can be calculated for the next iteration step. The new temperature guess $T_{\rm new}$ is obtained by

$\begin{displaymath}T_{\rm new}=\frac{E_{\rm err} T_{\rm old} - E_{\rm err_{\rm old}}T_{\rm cur}}{E_{\rm err}-E_{\rm err_{\rm old}}}, \end{displaymath}$

(12)

where $T_{\rm cur}$ and $E_{\rm err}$ are the current temperature guess and energy error. The variables $T_{\rm old}$ and $E_{\rm err_{\rm old}}$ are the temperature and energy error of the temperature iteration step before. With the new temperature guess we solve the radiative transfer equation again and check whether the total energy density is the desired one. It takes approximately four or five iteration steps to get the correct new temperature for a typical time step. The energy density is correct within an accuracy of 10^-5.

3 Test light curves

For the first test calculations we have solved the radiative transfer equation for a grey test atmosphere. Test model atmospheres are divided in 100 layers and we consider here time independent radiative transfer (in Sect. 4 below we discuss time dependent radiative transfer).

As a first test we applied our time evolution code to a static atmosphere. The test atmosphere is not expanding and no energy sources are present. Inside the test atmosphere we used a ``lightbulb'' radiating with a constant luminosity to simulate the internal energy flow from a star. We assumed an approximate temperature structure for t=0 and then let the atmosphere evolve towards radiative equilibrium. The atmosphere changed until it reaches steady state and the luminosity is constant in both space and time. The resulting final temperature structure should be identical to the structure for radiative equilibrium computed directly.

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f01.eps}} \end{figure}$	Figure 1: Three different light curves for the evolution to a stationary state. The three models have different luminosities produced by differentinner ``lightbulbs''.
Open with DEXTER

In Fig. 1 we show the light curves for three different static models. We can observe the model on its way towards radiative equilibrium. All calculations were started from the same temperature structure. After a certain time, the radiative relaxation time scale, each atmosphere has the (constant) luminosity of the lightbulb throughout the configuration.

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f02.eps}} \end{figure}$	Figure 2: The final temperature structures of the time evolution calculations.
Open with DEXTER

The final temperature structure of the model atmosphere with the lightbulb with a luminosity of 10³¹ erg/s is displayed in Fig. 2. Also shown is the result of a calculation using the temperature correction procedure of PHOENIX (Hauschildt et al. 2003) to directly compute the radiative equilibrium structure of the configuration. As one can clearly see, the resulting temperature structure of both methods agrees very well, and the maximum deviation of the temperature structure is less than one percent. Only the inner three layers have a deviation of up to 10 percent due to the implicit assumption of a diffusion approximation (and not a lightbulb) in the static calculation.

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f03.eps}} \end{figure}$	Figure 3: The flux in every layer at each time step.
Open with DEXTER

The next test is to look at time varying atmospheres. As an example we considered an atmosphere with a sinusoidally varying lightbulb inside. The resulting luminosity in each layer at each time point is plotted in Fig. 3. The luminosity in each layer is sinusoidal. It takes some time for the radiation to reach the outside boundary of the atmosphere and this results in a phase shift compared to the lightbulb.

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f04.eps}} \end{figure}$	Figure 4: The result for a sinusoidal varying lightbulb. Shown here is the luminosity in different layers. The phase shift between the lightbulb and the emergent flux is roughly $\pi$ .
Open with DEXTER

We now take a look at the luminosity in different layers. These are plotted in Fig. 4. One can see the phase shift of the sine in each layer. As expected, the luminosity from the central source needs time to move through the atmosphere.

One can also see that the amplitude of the sine decreases with increasing radius. What one would expect is that the amplitude is the same in every layer. The radiation is moving outwards and the incoming varying luminosity from the lightbulb should move through every layer and therefore the amplitude is supposed to be the same in each layer. If we integrate the luminosity over a whole sine, it stays constant in every layer. In the plot one can see that the mean luminosity is at the same level, so the luminosity is preserved. But why does the amplitude decrease? A possible explanation is that the radiation moving through the atmosphere is going backwards when the sine is declining, so this smears out the amplitude. If the sine varying radiation is moving through a very thick atmosphere, the amplitude at the top of the atmosphere is finally flat. An observer sees the mean luminosity.

For the next test we consider an atmosphere with an internal energy source. The initial structure is the radiative equilibrium structure of the static model with the lightbulb with a luminosity of 10³¹ erg/s. We assume a constant energy deposition rate in each layer of the model atmosphere. The luminosity is expected to increase over time and towards the outside. Figure 5 shows a plot of the light curve of this test atmosphere. The luminosity increases in time because of the energy deposition.

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f05.eps}} \end{figure}$	Figure 5: The light curve of an lightbulb with an additional energy source. An energy source in each layer causes an increasing luminosity of the model atmosphere.
Open with DEXTER

3.1 Expansion

In the case of SN Ia, homologous expansion is a good assumption after the initial breakout. In our test model each layer has a constant velocity, $u\propto r$ , to simulate a freely expanding envelope. The new radius $r_{\rm new}$ of a layer is, therefore, determined by

$\begin{displaymath}r_{\rm new}=u\cdot\Delta t+r_{\rm old} \end{displaymath}$

(13)

for a time step $\Delta t$ , while the layer is expanding with a velocity u. The radius before the new time step is $r_{\rm old}$ . With the same assumption of homologous expansion the new density $\rho_{\rm new}$ of a layer after at the new time step is determined by

$\begin{displaymath}\rho_{\rm new}=\rho_{\rm old}\cdot\left(\frac{r_{\rm old}}{r_{\rm new}}\right)^{3} \end{displaymath}$

(14)

where $\rho_{\rm old}$ is the old density. With the new radius and density, we are now able to calculate the new thermodynamics of the atmosphere. We assume the expansion of the supernova is an adiabatic process. The internal energy of the material changes due to work ${\rm d}W$ done during the expansion

$\begin{displaymath}{\rm d}E ={\rm d}W = - p {\rm d}V, \end{displaymath}$

(15)

where p is the pressure and ${\rm d}V$ a volume change. The energy conservation equation considers the energy density. The change in the energy density for the adiabatic expansion is given by

$\begin{displaymath}{\rm d}E_{\rm adia} = \frac{{\rm d}W}{m}= - \frac{p}{m} {\rm d}V. \end{displaymath}$

(16)

The mass m of a layer is given by its volume and density

$\begin{displaymath}m=V\cdot \rho. \end{displaymath}$

(17)

In our approach to a homologous expanding supernova of the type Ia, we consider each layer has a constant mass, so the derivative ${\rm d}V$ is given by

$\begin{displaymath}{\rm d}V = -\frac{m}{\rho^{2}} {\rm d}\rho. \end{displaymath}$

(18)

Together with the equation of state

$\begin{displaymath}p=\frac{R}{\mu}T\rho, \end{displaymath}$

(19)

we obtain the work for the adiabatic expansion as

$\begin{displaymath}\frac{{\rm d}W}{m}=\frac{R}{\mu}T\frac{{\rm d}\rho}{\rho}\cdot \end{displaymath}$

(20)

For differences in a time step $\Delta t$ , the work is

$\begin{displaymath}W_{\rm adia}= \Delta \frac{W}{m}=\frac{R}{\mu}T \ln \left( \frac{\rho_{2}}{\rho_{1}} \right)\cdot \end{displaymath}$

(21)

For the first test of a light curve for a supernova with an expanding atmosphere, we neglect the interaction between the layers, meaning the structure term is equal to the work of the adiabatic expansion so that

$\begin{displaymath}\frac{\rm D}{{\rm D}t}E=\frac{\rm D}{{\rm D}t}\left(\frac{E_{0}}{\rho}+e\right)=W_{\rm adia}. \end{displaymath}$

(22)

As we now consider expanding atmospheres, we use more supernova-like structures for the tests. We set the maximum velocity to 30 000 km s^-1and use a radius of $5\times 10^{15}~$ cm. With this setup, we calculated a light curve for the expanding atmosphere.

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f06.eps}} \end{figure}$	Figure 6: Light curve of an atmosphere that is just expanding.
Open with DEXTER

Figure 6 shows a plot of the light curve of supernova test atmosphere that is simply expanding. We considered time-independent radiative transfer and a grey atmosphere. As one can see, the luminosity is decreasing, because the atmosphere is cooling down adiabatically.

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f07.eps}} \end{figure}$	Figure 7: Light curve of an atmosphere that is expanding and has an energy source. It has the typical shape of a light curve of a SN Ia. The luminosity rises due to the energy deposition. After the maximum we see the decline resulting from the expansion.
Open with DEXTER

Now we test a setup more closely resembling a real supernova light curve. Therefore, we take an initial atmosphere structure and add an energy source (radioactive decay) in each layer. The energy source exponentially decreases to simulate declining activity of the radioactive species. Figure 7 shows the plot of the light curve of this test, resulting in a light curve with a supernova-like shape. We have a rising part of the light curve at the beginning because of the energy deposition. After the maximum, the luminosity decreases due to the ongoing expansion and decreasing energy deposition. Of course this light curve is far from correct because the assumption of a grey atmosphere is a bad assumption for an SN Ia. But the tests show that the code behaves as expected.

3.2 Entropy

We calculated the entropy to test the code. Because energy is moving through the atmosphere, entropy is not conserved. Nevertheless, testing entropy conservation is an important test of the correctness of the code. Therefore, we calculated the entropy change in the case of pure adiabatic expansion. For testing, we neglected interaction between layers and gamma ray deposition and calculated the entropy change for a time step for the just expanding case.

To be consistent with our hydrodynamics equations, we deduced the entropy from the first law of thermodynamics. A change in the entropy during a time step is therefore given by

$\begin{displaymath}\frac{\Delta S}{mR}=\frac{3}{2}\frac{1}{\mu}\ln\left( \frac{T... ...\right)-\frac{1}{\mu}\ln\left(\frac{\rho_{2}}{\rho_{1}}\right) \end{displaymath}$

(23)

where 2 is the index of quantities at the new time and 1 the one of the old. For the integration of the temperature and the density term, we kept $\mu$ fixed. This is not correct, but it is simpler to solve and the resulting differences for the entropy are small.

The entropy given here is only the entropy of the gas, so to check that the expansion worked right we had to neglect the radiation energy density in our energy density conservation equation. Furthermore we ignored the interaction between the layers and neglected the structure term. We now have only a change in the energy due to the adiabatic expansion.

With this setup we calculated a few time steps and checked the entropy. Even for a long time step of 1000 s the entropy stays almost constant. The relative change of the entropy was $\approx$ 10^-6.

4 Time-dependent radiative transfer

So far, our radiative transfer code only solves the time-independent radiative transfer equation. For an implementation of the time dependence in the radiative transfer itself, the spherical symmetric special relativistic radiative transfer equation (SSRTE) for expanding atmospheres (Hauschildt & Baron 1999) is extended so that the additional time dependent term is given by

$\begin{displaymath}\frac{\gamma}{c}\left(1+\beta\mu\right)\frac{\partial I}{\partial t} \end{displaymath}$

(24)

where $\beta=\frac{v}{c}$ is the velocity in units of the speed of light c and $\gamma=(1-\beta^{2})^{-1/2}$ is the usual Lorentz factor. Here, I is the intensity, $\mu$ the cosine of the angle between the radial direction and the propagation vector of the light. Using the notation of Hauschildt & Baron (2004a), the comoving frame SSRTE with the additional time dependent term is given by

$\begin{displaymath}a_{t}\frac{\partial I}{\partial t} +a_{r}\frac{\partial I}{\p... ...artial \lambda I}{\partial \lambda} +4a_{\lambda}I=\eta-\chi I \end{displaymath}$

(25)

where $\eta$ is the emissivity and $\chi$ the extinction coefficient. The wavelength is represented by $\lambda$ . The additional time dependent coefficient is given by

$\begin{displaymath}a_{t}=\frac{\gamma}{c}\left(1+\beta\mu\right). \end{displaymath}$

(26)

Along the characteristics the equation has the form

$\begin{displaymath}\frac{{\rm d}I_{l}}{{\rm d}s}+a_{t}\frac{\partial I}{\partial... ... \lambda I}{\partial \lambda} =\eta_{l}-(\chi_{l}+4a_{l})I_{l} \end{displaymath}$

(27)

where ${\rm d}s$ is a line element along a (curved) characteristic and I_l() is the specific intensity along the characteristic at point $s \ge 0$ (s=0 denotes the beginning of the characteristic). For a definition of the other coefficients see Hauschildt & Baron (2004a).

4.1 Discretization of the time derivative

We used the first discretization as described in Hauschildt & Baron (2004a) and added the time dependent term. We discretized the time-dependent, as well as the wavelength, derivative in the SSRTE with an fully implicit method. The discretization of the time dependent term is given by

$\begin{displaymath}\left.\frac{\partial I}{\partial t}\right\vert _{t=t_{j}} =\frac{I_{t_{j}}-I_{t_{j-1}}} {t_{j}-t_{j-1}}, \end{displaymath}$

(28)

so the SSRTE including the time discretization can be written as

$\begin{displaymath}\frac{{\rm d}I}{{\rm d} s}+a_{\lambda}\frac{\lambda_{l}I-\lam... ..._{j-1}}=\eta_{\lambda_{l}}-(\chi_{\lambda_{l}}+4a_{\lambda})I, \end{displaymath}$

(29)

where I is the intensity at wavelength point $\lambda_l$ and time point t_j. We define the optical depth scale along the ray as

$\begin{displaymath}{\rm d}\tau=\chi+a_{\lambda}\left(4+\frac{\lambda_{l}}{\lambd... ...mbda_{l-1}}\right)+\frac{a_{t}}{\Delta t}=-\hat{\chi}{\rm d}s. \end{displaymath}$

(30)

Introducing the source function $S=\eta/\chi$ , we get

$\begin{displaymath}\frac{{\rm d} I}{{\rm d}\tau}=I-\frac{\chi}{\hat{\chi}}\left(... ...}{\chi}\frac{1}{\Delta t}I_{t_{j-1}} \right) \equiv I-\hat{S}. \end{displaymath}$

(31)

We now have a modification of the source function and a new definition of the optical depth scale along a ray. With this redefinition of the optical depth and the source function, one can now proceed with the formal solution as described in Hauschildt & Baron (2004a).

4.2 Test calculations: time-dependent radiative transfer

For the test of the time-dependent radiative transfer, we used a static atmosphere structure to see the direct effects of the time dependence of the radiative transfer. Therefore, the temperature, radius, and density are all constant in time. For the test we then changed the inner boundary condition for the radiation (the ``lightbulb'') to initiate a perturbation of the radiation field, which then moves through the atmosphere via the time-dependent radiation transfer.

For the first test we switched on an additional light source inside the atmosphere. This light source has a luminosity of 10⁹ times higher then the inner lightbulb.

$\begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{10982f08.eps}} \end{figure}$	Figure 8: Luminosity of each layer and time step of our test atmosphere.
Open with DEXTER

In Fig. 8 we show the result of the time-dependent radiative transfer calculation. One can see that the light propagates outwards through the atmosphere. It takes a while before the additional light has propagated everywhere throughout the atmosphere. In Fig. 9 the light curves of a few layers are shown. One can see that the radiation needs time to get to the outer layers.

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f09.eps}} \end{figure}$	Figure 9: The light curves of a few layers. In this case the atmosphere has an additional light source switching on inside. The radiation needs time to get to the outer layers.
Open with DEXTER

We compared our time scale to the radiation diffusion time scale. Assuming a random-walk process the mean free path for a photon is given by $\lambda_{\rm p}=\frac{1}{\bar{\chi}}$ , where $\bar{\chi}$ is the mean opacity. For a travel distance l the time $t_{\rm p}$ a photon needs is given by

$\begin{displaymath}t_{\rm p}\approx \frac{1}{3}\frac{l^{2}}{c\lambda_{\rm p}}=\frac{1}{3}\frac{l^{2}}{c}\bar{\chi}, \end{displaymath}$

(32)

where c is the speed of light. For the mean opacity $\bar{\chi}$ we used the Rosseland mean. In our test case the atmosphere was divided into 128 layers. The mean opacity $\bar{\chi}$ ranges from $3\times 10^{-19}$ cm^-1 in the outer parts to $3\times 10^{-14}$ cm^-1 in the inner parts of the atmosphere. The distance l is the thickness of each layer, and it ranges between $6\times10^{12}$ cm and $3\times 10^{13}$ cm. The result is that the diffusion time for a photon through the whole model atmosphere is about $3\times 10^{3}$ s. As one can see in our plots our time scale is roughly 10⁵. The assumption of a diffusion through the atmosphere is only valid for optically thick regions. Another problem is the choice of the right mean opacity. Considering this our estimate of the time scale is reasonable.

An important test is also to check that the results of the time-dependent radiative transfer calculation depend on the size of the time step. We tested this with a model that has a small perturbation inside that is moving outwards. This was calculated with two different time steps. In Fig. 10 the results of calculations with two different time steps are shown. As one can see, the result does not depend on the size of the time step.

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f10.eps}} \end{figure}$	Figure 10: The light curve of an atmosphere where a small perturbation moves outwards. The two plotted light curves of layer 110 were calculated with different time steps.
Open with DEXTER

In the next test we put a sinusoidally varying light bulb inside of our test atmosphere. In Fig. 11, we show a plot of the luminosity for each time step. After a while, the luminosity of the whole atmosphere varies sinusoidally and steady state is reached.

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f11.eps}} \end{figure}$	Figure 11: Luminosity in each layer and time step. The light source inside the atmosphere varies sinusoidally.
Open with DEXTER

The luminosities in different layers are shown in Fig. 12. One can see a phase shift of the sine curve, due to the time required for the radiation field to propagate through the atmosphere.

5 Conclusion

We have implemented a time evolution code into our general-purpose model atmosphere code PHOENIX, which keeps track of the energy conservation. Because a homologous expansion is a good assumption for supernovae in general and particularly for type Ia supernovae, we considered adiabatic free expansion for our code. With first test calculations, we reproduced the expected behavior of the test cases. Static atmospheres adopted the luminosity of an internal lightbulb and perturbations of the lightbulb moved outwards in time. We also calculated a light curve that had the shape of a typical supernovae Ia light curve.

$\begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{10982f12.eps}} \end{figure}$	Figure 12: Flux at different radii of a sinusoidally varying atmosphere.
Open with DEXTER

To complete the physics of time-dependent processes in a supernova atmosphere, we extended our code even further. PHOENIX now solves the time dependent radiative transfer equation. We presented the new discretization scheme of the time dependent $\frac{\partial}{\partial t}$ term. We checked our code with a series of tests. A perturbation can be followed on its way through the atmosphere. We ran a test model with a sinusoidal source. In steady state the whole atmosphere was varying sinusoidally, responding to the source. A phase shift between the inner and outer layers could also be observed. All tests indicate that our extended code works fine.

Future work is to calculate realistic light curves for type Ia supernova hydro models. As a first step we can consider the atmosphere to be in LTE, but including detailed opacities and treating full line blanketing.

Also with our time-dependent radiative transfer, we can address a longstanding debate about the importance of time dependence in calculating the spectra of type Ia supernovae (Eastman & Pinto 1993; Baron et al. 1996). One challenge is the computation time needed for a whole light curve with more complex NLTE radiative transfer.

Acknowledgements

This work was supported in part by the Deutsche Forschungsgemeinschaft (DFG) via the SFB 676, NASA grant NAG5-12127, NSF grant AST-0707704, and US DOE Grant DE-FG02-07ER41517. This research used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, and the Höchstleistungs Rechenzentrum Nord (HLRN). We thank all these institutions for generous allocations of computer time.

References

Baron, E., Hauschildt, P. H., & Mezzacappa, A. 1996, MNRAS, 278, 763 [NASA ADS]
Bongard, S., Baron, E., Smadja, G., Branch, D., & Hauschildt, P. 2006, ApJ, 647, 480 [NASA ADS] [CrossRef]
Bronder, T. J., Hook, I. M., Astier, P., et al. 2008, A&A, 477, 717 [NASA ADS] [CrossRef] [EDP Sciences]
Cooperstein, J., van den Horn, L. J., & Baron, E. 1986, ApJ, 309, 653 [NASA ADS] [CrossRef] (In the text)
Eastman, R., & Pinto, P. 1993, ApJ, 412, 731 [NASA ADS] [CrossRef]
Foley, R. J., Filippenko, A. V., & Jha, S. W. 2008, A&A [ArXiv:0803.1181]
Goldhaber, G., Groom, D. E., Kim, A., et al. 2001, ApJ, 558, 359 [NASA ADS] [CrossRef]
Hachinger, S., Mazzali, P. A., & Benetti, S. 2006, MNRAS, 370, 299 [NASA ADS]
Hauschildt, P. H., & Baron, E. 1999, J. Comp. Appl. Math., 109, 41 [NASA ADS] [CrossRef]
Hauschildt, P. H., & Baron, E. 2004a, A&A, 417, 317 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
Hauschildt, P. H., & Baron, E. 2004b, Mitteilungen der Mathematischen Gesellschaft in Hamburg, 24, 1
Hauschildt, P. H., Barman, T. S., Baron, E., & Allard, F. 2003, in Stellar Atmospheric Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner (San Francisco: ASP), 227 (In the text)
Le Du, J., et al. 2009, preprint
Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (Oxford: Oxford University) (In the text)
Nugent, P., Phillips, M., Baron, E., Branch, D., & Hauschildt, P. 1995, ApJ, 455, L147 [NASA ADS] [CrossRef]
Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565 [NASA ADS] [CrossRef]
Phillips, M. M. 1993, ApJ, 413, L105 [NASA ADS] [CrossRef]
Pinto, P. A., & Eastman, R. G. 2000, ApJ, 530, 744 [NASA ADS] [CrossRef] (In the text)
Pskovskii, I. P. 1977, Sov. Astron., 21, 675 [NASA ADS]
Riess, A., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009 [NASA ADS] [CrossRef]
Riess, A. G., Press, W. H., & Kirshner, R. P. 1996, ApJ, 473, 88 [NASA ADS] [CrossRef]
Woosley, S. E., Kasen, D., Blinnikov, S., & Sorokina, E. 2007, ApJ, 662, 487 [NASA ADS] [CrossRef] (In the text)

All Figures

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f01.eps}} \end{figure}$	Figure 1: Three different light curves for the evolution to a stationary state. The three models have different luminosities produced by differentinner ``lightbulbs''.
Open with DEXTER
In the text

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f02.eps}} \end{figure}$	Figure 2: The final temperature structures of the time evolution calculations.
Open with DEXTER
In the text

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f03.eps}} \end{figure}$	Figure 3: The flux in every layer at each time step.
Open with DEXTER
In the text

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f04.eps}} \end{figure}$	Figure 4: The result for a sinusoidal varying lightbulb. Shown here is the luminosity in different layers. The phase shift between the lightbulb and the emergent flux is roughly $\pi$ .
Open with DEXTER
In the text

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f05.eps}} \end{figure}$	Figure 5: The light curve of an lightbulb with an additional energy source. An energy source in each layer causes an increasing luminosity of the model atmosphere.
Open with DEXTER
In the text

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f06.eps}} \end{figure}$	Figure 6: Light curve of an atmosphere that is just expanding.
Open with DEXTER
In the text

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f07.eps}} \end{figure}$	Figure 7: Light curve of an atmosphere that is expanding and has an energy source. It has the typical shape of a light curve of a SN Ia. The luminosity rises due to the energy deposition. After the maximum we see the decline resulting from the expansion.
Open with DEXTER
In the text

$\begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{10982f08.eps}} \end{figure}$	Figure 8: Luminosity of each layer and time step of our test atmosphere.
Open with DEXTER
In the text

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f09.eps}} \end{figure}$	Figure 9: The light curves of a few layers. In this case the atmosphere has an additional light source switching on inside. The radiation needs time to get to the outer layers.
Open with DEXTER
In the text

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f10.eps}} \end{figure}$	Figure 10: The light curve of an atmosphere where a small perturbation moves outwards. The two plotted light curves of layer 110 were calculated with different time steps.
Open with DEXTER
In the text

$\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{10982f11.eps}} \end{figure}$	Figure 11: Luminosity in each layer and time step. The light source inside the atmosphere varies sinusoidally.
Open with DEXTER
In the text

$\begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics{10982f12.eps}} \end{figure}$	Figure 12: Flux at different radii of a sinusoidally varying atmosphere.
Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.