Analytical Model of Time-Dependent Ionization in the Envelopes of Type II Supernovae at the Photospheric Phase

TL;DR

This paper presents an analytical model of time-dependent hydrogen ionization in Type IIP supernova envelopes, highlighting the importance of non-equilibrium effects during the photospheric phase.

Contribution

It introduces a simplified kinetic system and derives an analytical Lyapunov function to describe ionization freeze-out, emphasizing the need for time-dependent modeling.

Findings

Analytical expression for ionization freeze-out effect.

Recombination does not reach equilibrium during the photospheric phase.

Time-dependent effects are crucial for accurate supernova envelope modeling.

Abstract

We investigate a simplified kinetic system of the hydrogen atom (two levels plus continuum) under conditions of a type IIP supernova at the plateau phase that realistically describes the basic properties of the complete system. We have found the Lyapunov function for the reduced system using which we have analytically obtained the ionization freeze-out effect on long time scales. Since the system completely recombines in the equilibrium approximation on long time scales, which does not occur in reality, this result confirms the necessity of allowance for the time-dependent effect in the kinetics during the photospheric phase in a supernova explosion.

Figures

Equations

$$\frac{\partial \rho }{\partial t}=-▽(\rho v),$$

$$\frac{D\rho }{Dt}=-\rho (▽\cdot v)$$

$$\frac{D\rho }{Dt}+\frac{3\rho }{t}=0.$$

$$\frac{D{n}_i}{Dt}+\frac{3n_i}{t}=K_i(t),$$

$$K_1(t)=(N(t)-n_1-n_e)(Q+A_{21})+n_1{B}_{12}{J}_{12}(t)$$

$$K_e(t)=(N(t)-n_1-n_e)P_{2c}(t)-n_e^2{R}_{c2}(t).$$

$$N(t)=N_0\frac{t_0^3}{t^3}$$

$$\tau (t)\sim \frac{c^3}{8\pi }\frac{A_{21}}{{\nu }_{L\alpha }^3}\frac{g_2}{g_1}{n}_{1}t,$$

$$J_{12}(t)=(1-\beta (t))S(t)+\beta (t)J_c({\nu }_{L\alpha },t),$$

$$\beta (t)=\frac{1-e^{-\tau (t)}}{\tau (t)}\sim \frac{1}{\tau (t)}.$$

$$S_{12}(t)=\frac{2h{\nu }_{L\alpha }^3}{c^2},\left(\frac{g_2{n}_1}{g_1{n}_2}\right)$$

$$\dot{n_1}=(N(t)-n_1-n_e)(Q+A_{21}\beta (t))-n_1{B}_{12}\beta (t)J_c({\nu }_{L\alpha },t)-\frac{3n_1}{t}$$

$$\dot{n_e}=(N(t)-n_1-n_e)P_{2c}(t)-n_e^2{R}_{c2}(t)-\frac{3n_e}{t}$$

$$P_{2c}(t)=4\pi {\int }_{{\nu }_2}^{\infty }{\alpha }_{12}(\nu )\frac{J_c(\nu ,t)}{h\nu }d\nu .$$

$$\begin{array}{rl}{R}_{c2}=& 4\pi {\Phi }_{\mathrm{Saha}}(T_e){\int }_{{\nu }_2}^{\infty }{\alpha }_{12}(\nu )\frac{1}{h\nu }\frac{2h{\nu }^3}{c^2}{e}^{-\frac{h\nu }{k{T}_e}}d\nu \sim \\ & {\Phi }_{\mathrm{Saha}}(T_e)g_{II}(1,{\nu }_{L\alpha })\frac{\pi }{c^2}{E}_1\left(\frac{h{\nu }_2}{k{T}_e}\right)\end{array}$$

$$u_1=\frac{n_1}{N(t)}=\frac{n_1}{N_0}\frac{t^3}{t_0^3},u_e=\frac{n_e}{N(t)}=\frac{n_e}{N_0}\frac{t^3}{t_0^3},$$

$$\dot{u_1}=(1-u_1-u_e)\left(Q+\frac{A}{u_1}{\left(\frac{t}{t_0}\right)}^2\right)-\tilde{B}{J}_c({\nu }_{L\alpha },t){\left(\frac{t}{t_0}\right)}^2$$

$$\dot{u_e}=(1-u_1-u_e)P_{2c}(t)-u_e^{2}R{\left(\frac{t_0}{t}\right)}^3$$

$$A=\frac{8\pi {\nu }_{L\alpha }^3}{c^3}\frac{g_1}{g_2}\cdot \frac{1}{N_0{t}_0},\tilde{B}=\frac{4\pi }{hc}\cdot \frac{1}{N_0{t}_0},R=N_0{R}_{c2}$$

$$W(t)\sim \frac{1}{4}{\left(\frac{R_{ph}}{Vt}\right)}^2$$

$$J_c({\nu }_{L\alpha },t)=J_c({\nu }_{L\alpha },t_0){\left(\frac{t_0}{t}\right)}^{s_1}$$

$$P_{2c}(t)=P_{2c}(t_0){\left(\frac{t_0}{t}\right)}^{s_2}=P{\left(\frac{t_0}{t}\right)}^{s_2}$$

$${\dot{u_1}}^{td}=(1-u_1^{td}-u_e^{td})\left(Q+\frac{A}{u_1^{td}}{\left(\frac{t}{t_0}\right)}^2\right)-B{\left(\frac{t_0}{t}\right)}^{s_1-2}$$

$${\dot{u_e}}^{td}=(1-u_1^{td}-u_e^{td})P{\left(\frac{t_0}{t}\right)}^{s_2}-(u_e^{td}{)}^{2}R{\left(\frac{t_0}{t}\right)}^3,$$

$$(1-u_1^{ss}-u_e^{ss})\left(Q+\frac{A}{u_1^{ss}}{\left(\frac{t}{t_0}\right)}^2\right)-B{\left(\frac{t_0}{t}\right)}^{s_1-2}=0$$

$$(1-u_1^{ss}-u_e^{ss})P{\left(\frac{t_0}{t}\right)}^{s_2}-(u_e^{ss}{)}^{2}R{\left(\frac{t_0}{t}\right)}^3=0.$$

$$\dot{x}=-(x+y)Q-\frac{(1-{\tilde{u}}_e)x+{\tilde{u}}_{1}y}{{\tilde{u}}_1({\tilde{u}}_1+x)}A{\left(\frac{t}{t_0}\right)}^2$$

$$\dot{y}=-(x+y)P{\left(\frac{t_0}{t}\right)}^{s_2}-y(2{\tilde{u}}_e+y)R{\left(\frac{t_0}{t}\right)}^3.$$

$$V(t,x,y)=x^2+2xy+y^2\left(2+\frac{B^2}{A}\frac{1}{{\tilde{u}}_1}\frac{t_0^2}{t}\right).$$

$$\dot{V}(t,x,y)=\frac{\partial V}{\partial t}+\frac{\partial V}{\partial x}\dot{x}+\frac{\partial V}{\partial y}\dot{y},$$

Full text

K****eywords supernovae  $\cdot$ atmospheres  $\cdot$ spectral line formation

INTRODUCTION AND FORMULATION OF THE PROBLEM

New data, the photometric distances to objects with known redshifts, are required to investigate the present-day structure of the Universe. Among the great variety of distance measurement techniques there are methods that do not rely on the cosmological distance ladder, for example, the expanding photosphere method (EPM) [1], the spectral-fitting expanding atmosphere method (SEAM) [2], or the dense shell method (DSM) [3, 4, 5] that use type IIP and IIn supernovae (SNe) as objects. Using such a method as the SEAM requires the construction of a complete physical model for a type II SN that reproduces in detail its spectrum.

The importance of direct cosmological distance measurement methods is particularly topical in light of the problem of an uncertainty in measuring the Hubble parameter (Hubble tension) [6, 7, 8].

To completely model the physical processes occurring in a SN, it is necessary to simultaneously take into account the envelope expansion hydrodynamics, the matter–radiation field interaction, the radiative transfer in lines and continuum, and the kinetics of level populations in the atoms of a multiply charged plasma. This gives a system of integro-differential equations of radiation hydrodynamics that cannot yet be completely solved numerically even in the one-dimensional case. One has to resort to unavoidable simplifications in this complete system. One of such simplifications is the steady-state approximation of the kinetic system of level populations, when the system is assumed to be in statistical equilibrium.

The time-dependent hydrogen ionization effect in the envelopes of type II SNe at the photospheric phase was used by Kirshner and Kwan [9] to explain the high H$\alpha$ intensity in the spectra of SN 1970G and by Chugai [10] to explain the high degree of hydrogen excitation in the outer atmospheric layers ($v>7000$ km s*-1*) of SN 1987A in the first 40 days after its explosion.

Utrobin and Chugai [11] found a strong time-dependent effect in the ionization kinetics and hydrogen lines in type IIP SNe during the photospheric phase. In their next paper Utrobin and Chugai [12] also took into account the time-dependent effect in the energy equation. An important consequence of these papers was the conclusion that including the time-dependent ionization allowed the spectra of peculiar SN 1987A with a stronger H$\alpha$ line to be obtained, which could not be done previously without mixing radioactive 56Ni into the outer high-velocity layers in the steady-state approximation. In the next paper [13] the importance of this effect was also shown for normal SN 1999em.

The conclusions reached by Utrobin and Chugai were confirmed by Dessart and Hillier using the CMFGEN software package. In Dessart et al. [14] the applied approach was still the steady-state one and it was implemented in the CMFGEN package. Modeling revealed a problem: the H$\alpha$ line in hydrogen-rich envelopes was weaker than that observed at the recombination epoch. In particular, the model did not reproduce the line for times later than four days for SN 1987A and later than 20 days for SN 1999em. Next, Dessart and Hillier improved the code by including the time dependence in the kinetic system and the energy equation [15] and then in the radiative transfer [16, 17]. This allowed the H$\alpha$ line to be strengthened in the resulting spectrum, which led to better agreement with observations.

On the other hand, based on their computations with the PHOENIX software package, De et al. [18] found the time-dependent kinetics to be important only in the first days after SN explosion. Moreover, they argue that the role of the time-dependent effect is not very strong even in these first days by illustrating this with the models of SN 1987A and SN 1999em as an example. Using the open TARDIS code and without negating the importance of the time-dependent effect in the kinetics, Vogl et al. [19] nevertheless neglect it when modeling the spectra of SN 1999em and obtain good agreement of them with the observed ones. The overwhelming majority of Monte Carlo simulation codes also neglect the time-dependent effect in the kinetics. Thus, the conclusions of various research groups disagree and the importance of this effect is still called into question.

Using the STELLA and LEVELS codes, Potashov et al. [20] showed the importance of the time-dependent kinetics in the purely hydrogen case for SN 1999em as an example. The influence of metal admixtures on the intensity of the effect was also investigated: an increase in the concentration of metals in the envelope led to a weakening of the time-dependent effect in the kinetics.

In this paper we will attempt to answer the question of whether the time-dependent ionization effect is important or not, at least on long time scales, within a simple analytical model.

MODELING

Let us describe the construction of a simple analytical model for the behavior of multiply charged plasma electron populations in a SN envelope. We consider a purely hydrogen envelope, where the hydrogen atom is represented by a “two levels + continuum” system. We assume an $l$-equilibrium for the second atomic level. This means that the populations of sublevels including the fine structure of $2s$ and $2p$ are proportional to their statistical weights. Thus, the second level is considered as a single superlevel [21].

The characteristic stages of the behavior of the light curve for a typical type IIP SN can be written as follows [13]:

• •

shock breakout;

• •

the adiabatic cooling phase;

• •

the photospheric phase (a cooling and recombination wave is formed);

• •

the radiative diffusion cooling (the radiative diffusion time is less than the characteristic envelope expansion time);

• •

the beginning of thermal energy exhaustion;

• •

the end of thermal energy exhaustion (plateau tail phase);

• •

the nonthermal emission due to the decays 56Ni $\rightarrow$ 56Co $\rightarrow$ 56Fe (radioactive tail).

To study the time-dependent hydrogen ionization in the envelope plasma, we will consider the behavior of the system only at the photospheric phase. For typical SN 1999em [22, 13] this phase lasts from $t_{0}\sim 20$ to $t_{1}\sim 100$ days. A cooling and hydrogen recombination wave is formed as the envelope expands. The bolometric luminosity at the outer boundary of this wave is equal to the luminosity of the entire star. The photosphere is also located at the same boundary. An important fact is that in this period the photospheric radius $R_{ph}$, the radiation temperature $T_{c}$, and the matter temperature $T_{e}$ remain almost constant. Consequently, the total luminosity of the star does not change with time and the light curve reaches a plateau.

Our subsequent description suggests that $t\geqslant t_{0}$. Our modeling with the STELLA code shows that the transition to a homologous (with a high accuracy) expansion of the SN 1999em envelope ends approximately by day 15 after explosion [22], i.e., earlier than the beginning of the photospheric phase $t_{0}$. We assume an isotropic spherically symmetric expansion. In our initial analysis we also disregard the collisional excitation and ionization processes.

Let us select some small region of the envelope above the photosphere. The continuity equation in Eulerian coordinates for the matter in this region is

[TABLE]

where $\rho$ is the density of the envelope expanding with a velocity $v$. In the Lagrangian formalism in the comoving frame we obtain

[TABLE]

In the period of a free homologous expansion Eq. (2) is simplified to

[TABLE]

The rate of transitions to any discrete bound or free level $i$ of a hydrogen atom or ion can then be written as

[TABLE]

where $n_{i}$ is the population of level $i$ in the atom or ion. In turn, neglecting the stimulated emission processes, we define the function $K_{i}(t)$ as

[TABLE]

[TABLE]

Here,

[TABLE]

is the hydrogen number density; $Q$ is the two-photon $2s\to 1s$ decay rate; The reverse $1s\to 2s$ transition (two-photon absorption) rate is much lower than the $2s\to 1s$ rate and we disregard this process [20]; $A_{21}$ and $B_{12}$ are the Einstein coefficients for the spontaneous emission and photoexcitation of the $1\leftrightarrow 2$ transition; $J_{12}(t)$ is the intensity of radiation in the $2\to 1$ transition averaged over the line profile; $P_{2c}(t)$ is the total photoionization coefficient from the second level; $R_{c2}(t)$ is the total radiative recombination coefficient to the second level.

In our model we use the fact that the bound–free processes make a major contribution to the opacity in the Lyman continuum frequency band $\nu\geqslant\nu_{LyC}$ [20]. We neglect the relatively small contributions of the bound-bound processes in lines (the so-called expansion opacity) and free–free processes to the emission and absorption coefficients. The absorption in this band is caused mainly by neutral hydrogen and the optical depth is very large. Therefore, there is virtually no photospheric radiation here and the radiation field is determined for the regions above the photosphere by diffusive radiation. In this case, it can be shown that the photoionization rate from the ground hydrogen level and the recombination rate to the ground level closely coincide (even if there is not only hydrogen in the envelope). Thus, the first hydrogen level is in detailed balance with the continuum and the corresponding processes do not enter into the system of equations (5), (6).

It should be noted that in the Lyman continuum frequency band the intensity of continuum diffusive radiation $J_{c}(\nu)$ coincides with the equilibrium one $B_{\nu}(T_{e})$ only in a purely hydrogen envelope, suggesting a matter–radiation equilibrium. In the general case, with admixtures, $J_{c}(\nu)\neq B_{\nu}(T_{e})$.

Let us write out the standard formulas of Sobolev’s approximation [23, 24], but in a simplified form, using the condition that the population of the second level is relatively small, $N(t)-n_{1}-n_{e}\ll n_{1}$.

The Sobolev optical depth is

[TABLE]

the intensity of radiation averaged over the profile and the angles is

[TABLE]

where $J_{c}(\nu_{L\alpha})$ is the intensity of continuum radiation at the L$\alpha$ frequency.

We assume that $\tau(t)\gg 1$. The local escape probability of an L$\alpha$ photon without scattering integrated over the directions and over the line frequencies is then

[TABLE]

The source function is

[TABLE]

all the remaining notation is the standard one.

For an optically thick (in the L$\alpha$ line) SN envelope the estimate (10) breaks down. The absorption of photons in continuum should be taken into account [25, 26]. However, applying these corrections will not change qualitatively the final result of this paper.

Combining (4)–(6) and (9)–(11), we obtain the system

[TABLE]

[TABLE]

According to Mihalas [27] the photoionization rate is the integral

[TABLE]

while the photorecombination rate for a purely hydrogen plasma in the case where the stimulated emission is neglected will appear as

[TABLE]

Here $\Phi_{\mathrm{Saha}}(T_{e})$ is the Saha factor; $g_{II}(1,\nu_{L\alpha})$ is the Gaunt factor for the bound–free transition $1\leftrightarrow 2$, and $E_{1}$ is a modified exponential integral. It can be noted that the constancy of $R_{c2}$ in time follows from the constancy of $T_{e}$.

Let us now introduce dimensionless variables:

[TABLE]

which are the level populations normalized to the total current number density.

Rewriting the system in them, we obtain

[TABLE]

[TABLE]

New notation is also introduced here:

[TABLE]

Investigating the behavior of $J_{c}(\nu_{L\alpha},t)$ and $P_{2c}(t)$ with time is fundamentally important for a further simplification of system (16), (17). In the optically thin case, we can write $J_{c}(t)=W(t)B(T_{c})$. If, in addition, we assume that the region under consideration is sufficiently far from the photosphere, then the dilution factor will change with time as

[TABLE]

The continuum intensity $J_{c}(\nu_{L\alpha},t)$ and the photoionization rate $P_{2c}(t)$ will then drop as $\propto 1/t^{2}$.

In the observed SN the medium at the frequencies under consideration in continuum is optically thick. A large amount of metal admixtures changes the behavior of the radiation intensity in the hard band by reducing it. However, numerical simulations (for example, using the STELLA code) show a power-law time dependence of the intensity and the photoionization rate even in this case. Specifically,

[TABLE]

[TABLE]

The exponents $s_{1}$ and $s_{2}$ depend on the distance from the photosphere, but they are always greater than 2. Thus, generally, we restrict the domain of definition of the powers as $s_{1}\geqslant 2$ and $s_{2}\geqslant 2$.

System (16), (17), given (19), (20), will appear as

[TABLE]

[TABLE]

where $B=\tilde{B}J_{c}(\nu_{L\alpha},t_{0})$.

The equilibrium populations in the same approximation can be found by solving the following system of algebraic equations:

[TABLE]

[TABLE]

Thus, an answer to the question about the importance of allowance for the time-dependent effect in the kinetics should be sought by comparing the solutions of systems (21, 22) and (23, 24) or (equivalently) $u_{1}^{td}$, $u_{e}^{td}$ and $u_{1}^{ss}$, $u_{e}^{ss}$, respectively.

Below we show that any physically reasonable bounded solution of (21, 22) is Lyapunov stable “in the small”, i.e., stability is guaranteed at sufficiently small deviations.

Let we know one of the bounded solutions $0<\tilde{u}_{1}\leqslant 1$ and $0\leqslant\tilde{u}_{e}\leqslant 1$ (unperturbed motion) of the non-autonomous nonlinear differential system (21, 22). Let $x=\tilde{u}_{1}-u_{1}$ and $y=\tilde{u}_{e}-u_{e}$, i.e., $x$ and $y$ are the deviations of the solutions $u_{1}$, $u_{e}$ from $\tilde{u}_{1}$ and $\tilde{u}_{e}$, respectively.

For $x$ and $y$ we then obtain a reduced system of differential equations (it is called the Lyapunov system of equations of perturbed motion) [28, 29]:

[TABLE]

[TABLE]

It is important to note that the trivial solution $x=0$, $y=0$ is an equilibrium. Thus, investigating the Lyapunov stability of the solution $\tilde{u}_{1}$, $\tilde{u}_{e}$ is reduced to investigating the Lyapunov stability of the trivial solution (equilibrium position) $x=0$, $y=0$.

Let us next consider a scalar Lyapunov function of the following form:

[TABLE]

Obviously, (27) is positive definite for all instants of time. In view of the linearized system (25, 26), its time derivative can be written as

[TABLE]

where $\dot{x}$ and $\dot{y}$ are (25) and (26), respectively. It is a negative-definite quadratic form according to Sylvester’s criterion, because its corner minors when $t\rightarrow\infty$ are

[TABLE]

The derivative (29) itself is negative in sign on long time scales:

[TABLE]

Hence, according to Lyapunov’s first stability theorem [28, 29], the trivial solution of system (25, 26) is Lyapunov stable “in the small”. It should be noted that for the system under consideration the stability problem can be solved by this method on the semiaxis $t>\tilde{t}$ with a sufficiently distant boundary $\tilde{t}\geqslant t_{0}$. Stability on the pre-specified semiaxis $t>t_{0}$ is obtained by taking into account the known results on continuity in parameter [30] for the solution on the finite interval $t_{0}\leqslant t\leqslant\tilde{t}$.

It can be shown that system (25) and (26) is dissipative using Yoshizawa’s theorem [28, 31], i.e., the system is also stable “in the large”. Consequently, all solutions of system (21), (22) with physically reasonable initial conditions are bounded always.

Due to the boundedness in $u_{e}^{td}$, it follows from (22) that

[TABLE]

In turn, by solving (23), (24), it can be shown that

[TABLE]

Hence it follows that on long time scales the true relative electron number density reaches a constant, $u_{e}^{td}\sim c_{1}$, while the equilibrium relative electron number density approaches zero as $u_{e}^{ss}\sim t^{-(s_{1}+s_{2}-3)/2}$.

It can be seen that in the unsteady-state case the envelope expands with a higher degree of ionization than in the steady-state approximation. This should be taken into account when modeling the SN kinetics.

CONCLUSIONS

Such a phenomenon is also observed in atmospheric explosions [32, 33] and during the “protraction” of primordial plasma recombination in the early Universe under cosmological conditions [34, 35]. The number density of free electrons is commonly said to experience “freeze-out.” However, in contrast to the classical freeze-out in atmospheric explosions, the time-dependent effect in SNe remains even when the temperatures of both matter and radiation are constant.

In this paper we considered the breakdown of the equilibrium steady-state approximation itself in the kinetics. The magnitude of this breakdown and its evolution with time will be described in subsequent publications.

ACKNOWLEDGMENTS

We are grateful to N.N. Shakhvorostova, V.P. Utrobin, and A.V. Yudin for the stimulating discussions.

FUNDING

The work of M.Sh. Potashov was supported in part by RFBR grant no. 19–02–00567, while the work of S.I. Blinnikov on numerical SN models was supported by RSF grant no. 18–12–00522.