A simple model of time-dependent ionization in Type IIP supernova envelopes

TL;DR

This paper introduces a simple yet realistic model for time-dependent hydrogen ionization in Type IIP supernova envelopes, highlighting the persistent ionization freeze-out effect and providing criteria for non-equilibrium conditions.

Contribution

It presents a minimal kinetic model for hydrogen ionization in supernova atmospheres and generalizes criteria for non-equilibrium states, enhancing understanding of ionization dynamics during supernova evolution.

Findings

Ionization freeze-out occurs at large times in supernova envelopes.

Non-equilibrium ionization effects are inevitable if early conditions are non-equilibrium.

The model provides a simple criterion for assessing statistical equilibrium.

Abstract

We propose a model kinetic system of the hydrogen atom (two levels plus continuum) under the conditions typical for atmospheres of Type IIP supernovae in the plateau stage. Despite the simplicity of this system, it describes realistically the basic properties of the complete system. Analysis shows that the ionization "freeze-out" effect is always manifest at large times. We give a simple criterion for checking the statistical equilibrium of a system under the given conditions at any time. It is shown that if the system is in non-equilibrium at early times, the time-dependent effect of ionization necessarily exists. We also generalize this criterion to the case of arbitrary multilevel systems. We discuss various factors that affect the strength of the time-dependent effect in the kinetics during the photospheric phase of a supernova explosion.

Tables

Table 1. Table 1: Typical values of dimensionless constants for the system ( 23 , 24 ) obtained on the basis of STELLA calculations for SN 1999em (Baklanov et al., 2005 ; Potashov et al., 2017 ) , which is a standard SN IIP. The column labelled “Out” contains values that correspond to the outer, high-velocity layers located far from the photosphere. The column labelled “Ph” contains values that correspond to the near-photospheric layers. See Appendix  B for details.

	out	ph
$Q$	$3.5\cdot {10}^6$	$3.5\cdot {10}^6$
$A$	$3\cdot {10}^7$	$3\cdot {10}^4$
$B$	$6\cdot {10}^{-4}$	$6\cdot {10}^{-3}$
$P$	$7\cdot {10}^9$	$2\cdot {10}^{11}$
$R$	$50$	$5\cdot {10}^4$
$s_1$	$21$	$24$
$s_2$	$4$	$6$

Table 2. Table 2: Different kinetic systems.

Model	Elements	Levels	Collisional processes
$M_1$	H	2	no
$M_2$	H	2	yes
$M_3$	H	10	yes
$M_4$	H $+$ He $+$ metals	10	yes

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_{\mathrm{e}})(A_{2q}+A_{21})+n_1{B}_{12}{J}_{12}(t),$$

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

$$K_{\mathrm{e}}(t)=(N(t)-n_1-n_{\mathrm{e}})P_{2c}(t)-n_{\mathrm{e}}^2{R}_{c2}(t).$$

$$n_2=N(t)-n_1-n_{\mathrm{e}}\ll n_1,$$

$${\tau }_S(t)\simeq \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 }_S(t)}}{{\tau }_S(t)}\simeq \frac{1}{{\tau }_S(t)}.$$

$$S(t)\simeq \frac{2h{\nu }_{L\alpha }^3}{c^2}\left(\frac{g_1{n}_2}{g_2{n}_1}\right).$$

$${\dot{n}}_1=(N(t)-n_1-n_{\mathrm{e}})(A_{2q}+A_{21}\beta (t))-n_1{B}_{12}\beta (t)J_c({\nu }_{L\alpha },t)-\frac{3n_1}{t},$$

$${\dot{n}}_{\mathrm{e}}=(N(t)-n_1-n_{\mathrm{e}})P_{2c}(t)-n_{\mathrm{e}}^2{R}_{c2}(t)-\frac{3n_{\mathrm{e}}}{t}.$$

$$P_{2c}(t)=4\pi \underset{{\nu }_2}{\overset{\infty }{\int }}{J}_c(\nu ,t)\frac{{\alpha }_{2c}(\nu )}{h\nu }d\nu ,$$

$$R_{c2}=4\pi {\Phi }_{\mathrm{Saha}}(T_{\mathrm{e}})\underset{{\nu }_2}{\overset{\infty }{\int }}\frac{{\alpha }_{2c}(\nu )}{h\nu }\frac{2h{\nu }^3}{c^2}{e}^{-\frac{h\nu }{k{T}_{\mathrm{e}}}}d\nu \approx \frac{64{\pi }^{5}m{e}^{10}}{3\sqrt{3}{c}^3{h}^6}{g}_{\mathrm{II}}(2,{\nu }_2){\Phi }_{\mathrm{Saha}}(T_{\mathrm{e}})E_1\left(\frac{h{\nu }_2}{k{T}_{\mathrm{e}}}\right)$$

$$u_1(t)=\frac{n_1}{N(t)}=\frac{n_1}{N_0}\frac{t^3}{t_0^3},u_{\mathrm{e}}(t)=\frac{n_{\mathrm{e}}}{N(t)}=\frac{n_{\mathrm{e}}}{N_0}\frac{t^3}{t_0^3}$$

$${\dot{u}}_1=(1-u_1-u_{\mathrm{e}})\left[A_{2q}+\frac{\tilde{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}}_{\mathrm{e}}=(1-u_1-u_{\mathrm{e}})P_{2c}(t)-u_{\mathrm{e}}^2\tilde{R}{\left(\frac{t_0}{t}\right)}^3.$$

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

$$J_c(t)=W(t)B(T_c),$$

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

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

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

$$\tau =\frac{t}{t_0}.$$

$$\dot{u}_{1}=(1-u_{1}-u_{\mathrm{e}})\Big{(}Q+\frac{A}{u_{1}}\tau^{2}\Big{)}-\frac{B}{\tau^{s_{1}-2}},$$

$${\dot{u}}_{\mathrm{e}}=(1-u_1-u_{\mathrm{e}})\frac{P}{{\tau }^{s_2}}-u_{\mathrm{e}}^2\frac{R}{{\tau }^3}.$$

$$Q=A_{2q}{t}_0,A=\tilde{A}{t}_0,B=\tilde{B}{J}_c({\nu }_{L\alpha },t_0)t_0,P=\tilde{P}{t}_0,R=\tilde{R}{t}_0.$$

$$\begin{array}{rl}& 0<u_1(1)⩽1,0⩽u_{\mathrm{e}}(1)⩽1,\\ & u_2(1)=1-u_1(1)-u_{\mathrm{e}}(1)\ll u_1(1)(\text{mbox}seeeq.\text{leavevmode}\text{ref}eq:n{2}_{s}mall).\end{array}$$

$$(1-u_{1}^{ss}-u_{\mathrm{e}}^{ss})\Big{(}Q+\frac{A}{u_{1}^{ss}}\tau^{2}\Big{)}-\frac{B}{\tau^{s_{1}-2}}=0,$$

Full text

K****eywords atomic processes  $\cdot$ methods: analytical  $\cdot$ stars: atmospheres  $\cdot$ supernovae: general.

1 Introduction

The investigation of the structure of the Universe involves the measurement of photometric distances to objects with known redshifts. There are a large variety of methods for measuring distances. Among them, there are methods that do not rely on the cosmological distance ladder, such as the expanding photosphere method (EPM) (Kirshner & Kwan, 1974), the spectral-fitting expanding atmosphere method (SEAM) (Baron et al., 2004), and the dense shell dethod (DSM) (Blinnikov et al., 2012; Potashov et al., 2013; Baklanov et al., 2013). These methods are based on the properties of Type IIP and Type IIn supernovae (SNe IIP and SNe IIn). Some of them (e.g., SEAM) are very complicated and require full physical modelling of the SN with a detailed reproduction of its spectrum. Direct cosmological methods for distance measurement are especially important owing to the problem of the uncertainty in the Hubble parameter (Hubble tension) (Riess et al., 2018; Mörtsell & Dhawan, 2018; Ezquiaga & Zumalacárregui, 2018).

In order to model the physical processes occurring within SN ejecta, it is necessary to solve a system of partial integro-differential equations of radiation hydrodynamics, which includes the envelope expansion hydrodynamics, the interaction of the radiation field with matter, the radiation transfer in lines and in the continuum, and the kinetics of level populations in atoms of multiply charged plasma. The complete numerical solution of this system is still an impossible task even in the one-dimensional case, and one has to resort to unavoidable simplifications. An important and frequently used simplification is the steady-state approximation for the kinetic system of level populations, when the system is assumed to be in statistical equilibrium. The effect of time-dependence is manifest in the deviation of the actual occupation numbers of the atomic levels from their steady-state (equilibrium) values.

The time-dependent effect of hydrogen ionization in the envelopes of SNe II in the photospheric phase was used by Kirshner & Kwan (1975) to explain the high H$\alpha$ intensity in the spectra of SN 1970G, and by Chugai (1991) to explain the high degree of hydrogen excitation in the outer atmospheric layers ($v>7000$ km s*-1*) of SN 1987A for the first 40 d after the explosion.

Utrobin & Chugai (2002) found a strong time-dependent effect in the ionization kinetics and hydrogen lines in SNe IIP during the photospheric phase. In the next paper, Utrobin & Chugai (2005) also took into account the time-dependent effect in the energy equation. In these papers, it was shown that the H$\alpha$ line was enhanced in the spectrum of the peculiar SN 1987A owing to the time-dependent ionization. In the steady-state approximation, this effect has been achieved only by mixing radioactive 56Ni into the outer high-velocity layers. Utrobin (2007) showed the importance of non-stationarity also for a normal SN IIP, SN 1999em.

The conclusion of Utrobin and Chugai was confirmed by Dessart and Hillier using the CMFGEN software package. Dessart et al. (2008) still applied the steady-state approach that was implemented in the CMFGEN package, but the H$\alpha$ line in the hydrogen-rich envelopes was weaker than the observed one during the recombination epoch. In particular, the model did not reproduce H$\alpha$ after the SN age of 4 d for SN 1987A and 20 d for SN 1999em. Another version of CMFGEN was improved by including the time dependence in the kinetic system and the energy equation (Dessart & Hillier, 2007), and for the latest version, that in the radiative transfer (Dessart & Hillier, 2010; Hillier & Dessart, 2012). This strengthened the H$\alpha$ in the simulated spectrum and led to better agreement with observations.

On the other hand, based on the computations with the PHOENIX software package, De et al. (2010) found that the time-dependent kinetics is important only during the first few days after the SN explosion. Moreover, they argued that the role of the time-dependent effect is not very strong even in the first few days, and illustrated this with models of SN 1987A and SN 1999em.

Vogl et al. (2019), recognized the importance of time-dependent effect for the ionization kinetics, but neglected it. Nevertheless, they obtained a good agreement of the spectra of the SN 1999em simulated with the open source code TARDIS with the observed ones. The vast majority of Monte Carlo simulation codes also neglect the time-dependent effect in kinetics. Thus, the conclusions of the various research groups disagree, and the importance of the effect is still being debated.

Potashov et al. (2017) showed the importance of the non-stationary kinetics for SN 1999em in the purely hydrogen case using the codes STELLA (Blinnikov et al., 1998, 2000, 2006) and LEVELS. The influence of metal admixtures on the non-stationarity was also studied. The increase of the metal abundance in the envelope led to a weakening of the time-dependent effect.

In this paper, we use a simple analytical model to determine whether and when the time-dependent ionization effect is important. We give a simple criterion for checking the statistical equilibrium of a system. We distinguish two important periods in the evolution of the system under study. The first period starts from the beginning of the photospheric phase and lasts not longer than the time-scale for the change of the photospheric parameters, We refer to this period as “early time”. Typically, it lasts for a few days. Much longer periods will be called “large time”.

A brief outline of the paper is as follows. Section 2 describes the physical model of the problem. In Section 3, it is concluded that the time-dependent effect is significant at least in the large time limit. In Section 4, the evolution of deviation from stationarity with time is studied, assuming that initially it is small. A formula for this evolution is derived. This formula depends on an expression that defines the strength of the time-dependent effect. In Section 5, it is shown that the expression described in Section 4 gives a simple criterion for checking the statistical equilibrium of the system under the given conditions at any time. For the system in equilibrium, there will be no time-dependent effect. In this section, we also conclude that the key factor affecting the relaxation time of the entire system is the intensity and rate of change of the external hard continuum radiation between the Lyman and Balmer ionization thresholds. Finally, in Section 6 we investigate various kinetic systems using a generalized equilibrium criterion and find additional factors that influence the effect of non-stationary ionization. We also discuss how this effect is affected by the metallicity of the stellar envelope.

2 Physical model

We consider a reasonably simple analytical model for the behaviour of level populations of a purely hydrogen plasma in the supernova envelope. The hydrogen atom is represented by the system “two levels plus continuum”. We assume an $l$-equilibrium in the kinetic system for the second level. This means that the populations of the fine-structure sublevels $2\mathrm{s}$ and $2\mathrm{p}$ are proportional to their statistical weights. Thus, the second level is considered as a single so-called super-level (Hubeny & Lanz, 1995).

The light-curve behaviour of a typical SN IIP can be divided into various characteristic stages (Utrobin, 2007):

• •

shock breakout;

• •

adiabatic cooling phase;

• •

photospheric phase (cooling and recombination

wave);

• •

phase of radiative diffusion cooling;

• •

exhaustion of radiation energy;

• •

plateau tail;

• •

radioactive tail (56Ni $\Rightarrow$ 56Co $\Rightarrow$ 56Fe).

To study the time-dependent ionization of the hydrogen plasma, we will consider the behaviour of the system only in the photospheric phase. For the typical supernova SN 1999em (Baklanov et al., 2005; Utrobin, 2007) this phase lasts from ${t_{0}\sim 20}$ d to ${\sim 100}$ d. As the envelope expands, a cooling and recombination wave is formed. The bolometric luminosity of the SN is equal to the luminosity at the outer edge of this wave. The photosphere is located at the same level. It is important that during the photospheric phase, the photospheric radius $R_{\mathrm{ph}}$, the radiation temperature $T_{c}$ and the gas temperature $T_{\mathrm{e}}$ are nearly constant. Consequently, the luminosity of the SN does not change in time, and one observes a plateau in the light-curve.

A radiation-hydrodynamical simulation of the SN 1999em envelope with the code STELLA shows that the transition to homologous expansion (with high accuracy) is completed by about 15 d after the explosion (Baklanov et al., 2005), which is before the beginning of the photospheric phase $t_{0}$. We assume that the gas expands isotropically (i.e. a one-dimensional spherically symmetric approximation is used) and do not take collisional processes into account. The validity of the latter assumption is discussed in Section 6 (see also Appendix C).

Let us select a small area of the envelope above the photosphere. The continuity equation in Eulerian coordinates for the gas 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 free homologous expansion stage (${v\propto r}$), equation (2) is simplified:

[TABLE]

In this case, the rate of transitions to any discrete bound or free level $i$ of neutral or ionized hydrogen can be written as

[TABLE]

where $n_{i}$ is the population of level $i$ for an atom or ion. Neglecting the processes of induced emission for the ground level of hydrogen, we define the function $K_{1}(t)$ as

[TABLE]

Here

[TABLE]

is the hydrogen number density, and $A_{2q}$ is the two-photon decay probability ${2\mathrm{s}\to 1\mathrm{s}}$. The reverse ${1\mathrm{s}\to 2\mathrm{s}}$ transition (two-photon absorption) rate is much lower than the ${2\mathrm{s}\to 1\mathrm{s}}$ rate, and we neglect this process (Potashov et al., 2017); we also neglect stimulated two-photon decays. $A_{21}$ and $B_{12}$ are the Einstein coefficients for spontaneous and induced transition ${1\leftrightarrow 2}$, and $J_{12}(t)$ is the mean intensity of radiation for the ${2\to 1}$ transition averaged over the line profile.

Equation (4) with the index $\mathrm{e}$ instead of $i$ is valid for the number density of free electrons, with

[TABLE]

Here $P_{2c}(t)$ is the photoionization coefficient for the second level, and $R_{c2}(t)$ is the radiative recombination coefficient for the second level.

In our model, we will use the fact that the main contribution to the opacity in the frequency band of the Lyman continuum ${\nu\geqslant\nu_{LyC}}$ is provided by free-bound processes (Potashov et al., 2017). We neglect the relatively small contribution from bound-bound processes (expansion opacity) and free-free processes in the emission and absorption coefficients. The absorption in this band is mainly due to neutral hydrogen, and the optical depth is very large. Therefore, there is virtually no photospheric radiation, and the radiation field for the regions above the photosphere is determined by the diffusive radiation. In this case, it can be shown (see Appendix A) that the rates of photoionization transitions from the ground level of hydrogen and recombination to the ground level completely coincide (even if there is no pure hydrogen in the SN envelope). Thus, the ground level of hydrogen is in detailed balance with the continuum, and the related processes are not included in the system of equations (5 and 7).

Let us assume that the population of the second level is relatively small,

[TABLE]

and write down the standard formulas for the Sobolev approximation (Sobolev, 1960; Castor, 1970).

The Sobolev optical depth in L$\alpha$ is

[TABLE]

The frequency-averaged mean intensity of the transition ${1\leftrightarrow 2}$ is

[TABLE]

where $J_{c}(\nu_{L\alpha},t)$ is the average intensity of the continuum on the L$\alpha$ frequency.

Because the medium is optically thick for L$\alpha$ photons, ${\tau_{S}(t)\gg 1}$. Then the escape probability of L$\alpha$ photons, summed over all directions and line frequencies, is

[TABLE]

The source function is

[TABLE]

All other notations are standard.

It is known that the optical depth of the SN II envelope in the Lorentz wings is very large for L$\alpha$ (${a\tau_{S}\gg 1}$, where $a$ is the Voigt damping parameter), and the profile cannot be considered as the Doppler one. In this case, under the hypothesis of complete frequency redistribution, the formal criterion of the applicability of Sobolev theory is violated (Chugai, 1980). However, Chugai (1980) and Grachev (1989) have shown that the estimation (11) remains true for the conservative scattering of L$\alpha$ photons if one assumes partial frequency redistribution and neglects the effects of recoil. In these works, the Fokker-Planck approximation was used for the redistribution function. Hummer & Rybicki (1992) showed that the frequency-weighted mean intensity $J_{12}(t)$ defined in (10) is weakly dependent on the mechanism of frequencies redistribution. Expressions (10) and (11) remain correct in the case of non-conservative scattering with partial frequencies redistribution and partial non-coherence owing to the Stark effect (Chugai, 1988b).

When absorption in the continuum in the region of a line is significant (for example, L$\alpha$ photons can ionize Ca II from the second level), one must take into account corrections to the Sobolev approximation (Hummer & Rybicki, 1985; Chugai, 1987; Grachev, 1988). The additional selective absorption in lines of metal admixtures may also play a role, because a large number of lines of Fe II and Cr II are in the vicinity of L$\alpha$ (Chugai, 1988a, 1998). These processes increase the value of the effective escape probability (11) (the so-called loss probability of a photon in flight). The influence of absorption in the continuum and in the lines of admixtures on the time-dependence effect was investigated by Potashov et al. (2017). In the current analytical approach, we use a simple model and do not take into account these processes, using hereafter equations (10) and (11). However, in Section 6, these modifications of the Sobolev approximation will be taken into account numerically.

Combining (4–7, 10–12), we obtain the system

[TABLE]

[TABLE]

In accordance with Mihalas (1978) and Hubeny & Mihalas (2014, p. 273), the total photoionization coefficient for the second level is the integral

[TABLE]

where $\nu_{2}$ is a frequency of the ${2\leftrightarrow c}$ transition, $c$ means the continuum, and $\alpha_{2c}$ is a photoionization cross-section of the second level at the frequency $\nu$. The total radiative recombination coefficient of the second level for a purely hydrogen plasma looks like

[TABLE]

if the processes of induced emission are neglected. Here $\Phi_{\mathrm{Saha}}(T_{\mathrm{e}})$ is the Saha-Boltzmann factor; $g_{\,\mathrm{II}}(2,\nu_{2})$ is a bound-free Gaunt factor for the ${2\leftrightarrow c}$ transition; and $E_{1}$ is the exponential integral. Note that $R_{c2}$ is constant over time if $T_{\mathrm{e}}$ is constant. It would be more realistic to use the effective coefficient of recombination to the second level calculated for a multilevel atom taking into account cascade transitions instead of the direct coefficient of recombination to the same level (14). Usually, Lyman lines are supposed to be optically thick, and the escape probabilities for the Lyman line photons are zero, while all other lines are optically thin. This is the so-called Case B (Baker & Menzel, 1938; Osterbrock & Ferland, 2006) for the recombination coefficient calculation. In the opposite Case A (Baker & Menzel, 1938), full transparency is assumed and the escape probabilities for the Lyman line photons are equal to unity. In reality, the escape probability of a photon in Lyman line differs from both zero and unity. Therefore, neither Case A nor Case B is fully implemented, and the real situation is intermediate (Hummer & Storey, 1987, 1992). We choose the average value of $R_{c2}$ from Hummer (1994) for our simple system.

Let us now introduce the dimensionless variables

[TABLE]

that are normalized to the full current number density. By rewriting the system we obtain

[TABLE]

[TABLE]

Here we also introduce a new notation for constants:

[TABLE]

Especially important for further simplification of the system (15, 16) is an investigation of the behaviour of $J_{c}(\nu_{L\alpha},t)$ and $P_{2c}(t)$ over time. The continuum radiation between the Lyman and Balmer edges is assumed to be given, and not determined from the self-consistent calculation.

In the optically thin case we can write, introducing a dilution factor $W(t)$,

[TABLE]

If we additionally assume that the region under consideration is sufficiently far from the photosphere, the dilution factor is

[TABLE]

Then the continuum intensity $J_{c}(\nu_{L\alpha},t)$ and the photoionization coefficient $P_{2c}(t)$ decrease as ${\sim 1/t^{2}}$. This case is called the free-streaming approximation.

In a real SN, the medium at the considered frequencies in the continuum is optically thick. Numerous metal lines between the Lyman and Balmer ionization thresholds contribute to the large expansion opacity for SN envelopes with non-zero metallicities. The averaged intensity of such a quasi-continuum is lower than in the optically thin limit. Even in this case, numerical simulation (for example using the STELLA code) shows the power dependence of the intensity and the photoionization coefficient on time. Namely,

[TABLE]

[TABLE]

The values of $s_{1}$ and $s_{2}$ depend on the distance from the photosphere, but they are always greater than 2. In the general case, we limit the domain of these exponents to ${s_{1}\geqslant 2}$ and ${s_{2}\geqslant 2}$. In the optically thick case, the typical values of these exponents can be significant (see Table 1, Appendix B).

Now let us introduce the normalized time

[TABLE]

Taking into account (20), (21) and (22), we can rewrite the system (15, 16) as

[TABLE]

[TABLE]

Now $u_{1}$ and $u_{\mathrm{e}}$ are functions of $\tau$, and

[TABLE]

Table 1 shows typical values of the constants $Q$, $A$, $B$, $P$, $R$, $s_{1}$, $s_{2}$ that correspond to the physical conditions for the outer, high-velocity layers located far from the photosphere (column 1), and the near-photospheric layers (column 2) of SNe IIP, (see Appendix B for details).

The initial conditions for the problem are:

[TABLE]

Equilibrium populations in the same approximation can be found by solving a system of algebraic equations:

[TABLE]

[TABLE]

where the superscript $ss$ denotes steady state. Thus, the answer to the question regarding the importance of the time-dependence effect in kinetics can be found by comparing the solutions of the systems (23, 24) and (27, 28). In the next section, it will be shown that at large time the deviation of $u_{1}$, $u_{\mathrm{e}}$ from $u_{1}^{ss}$, $u_{\mathrm{e}}^{ss}$ persists regardless of the initial conditions.

3 The system at large time

We will investigate the solutions of systems (23, 24) and (27, 28) at large time ${\tau\to\infty}$. In this case, we will assume that the photospheric phase lasts for an indefinitely long time.

Let us transform the original system (23, 24) by introducing the function $u_{2}$ according to the relationship ${u_{1}{+}u_{2}{+}u_{\mathrm{e}}{=}1}$. We should keep only terms linear in $u_{2}$ because the smallness of this function was already assumed in the derivation (see 8 and 26). Thus, this is a reduction of the source system to a “normal” form:

[TABLE]

where we have introduced the following functions:

[TABLE]

It is important that, because of its linearity, equation (30) can be integrated explicitly and its solution can be written as

[TABLE]

where

[TABLE]

Equation (31) together with (32) represents a formal solution for $u_{2}$, because $u_{1}$ is still unknown.

Analysis of the system (29, 30) shows that the solutions that satisfy the initial conditions $0{<}u_{1}(1){\leqslant}1$, $0{\leqslant}u_{2}(1){\leqslant}1$ remain bounded and stable in the sense of Lyapunov \NAT@openDemidovich, 1967, p. 66; Khalil, 2002, p. 111\NAT@close, see also Appendix C. The functions $g_{1}(t)$ and $g_{3}(t)$ increase rapidly, ${g_{1,3}\sim\mathcal{O}(\tau^{2})}$, if $u_{1}$ is bounded and finite. In contrast, the functions $g_{2,4}$ are decreasing as some power of $\tau$. We can use the fact that for any bounded behaviour of $u_{1}(\tau)$, $G_{3}$ is a rapidly increasing function of time, at least as ${\sim\mathcal{O}(\tau^{3})}$. Therefore, terms with $\exp(G_{3})$ are most important at large time. The term with $u_{2}(1)$ in (31) is exponentially small at large time; that is, $u_{2}$ “forgets” the initial conditions. After multiple integration of the second term by parts and neglecting the exponentially small terms at ${\tau\to\infty}$, we obtain

[TABLE]

The first term of this decomposition is actually a steady-state approximation of equation (30). Substituting it into (29), we obtain the equation for $u_{1}$:

[TABLE]

Equation (34) is a special case of the Appell equation (Appell, 1889). That is, it is a generalized Abel equation of the second kind (Polyanin & Zaitsev, 2002; Semenov, 2014) and unfortunately it is generally not integrable by quadratures.

Usually, approximate analytical methods, such as the small parameter method (Fedoruk, 1985, p. 405), Chaplygin method (Berezin & Zhidkov, 1959, p. 260) or power series method, are used to solve nonlinear differential equations. Case (34) requires a large number of steps in each particular method or many summands in the decomposition. The analytical approach becomes cumbersome and useless. But we do not need to solve (34). It is sufficient to prove that any solution of this equation that satisfies the initial conditions ${0{<}u_{1}(1){\leqslant}1}$ is bounded in the interval ${(0,1)}$. To do this, consider two functions

[TABLE]

[TABLE]

Let us take an arbitrary moment of time ${\tau_{1}\gg 1}$. In the domain (${\tau\in[1,\tau_{1}];\;u\in\mathbb{R}}$), the inequalities $f_{l}(\tau,u)\leqslant f(\tau,u)\leqslant f_{u}(\tau,u)$ are obviously satisfied. Because the function $f(\tau,u)$ is continuously differentiable in this domain and satisfies the Lipschitz condition, according to the Chaplygin theorem on differential inequalities \NAT@openBerezin & Zhidkov, 1959, p. 260; Khalil, 2002, p. 102\NAT@close, one can write ${u_{l}(\tau)\leqslant u_{1}(\tau)\leqslant u_{u}(\tau)}$ for any time in ${[1,\tau_{1}]}$, where $u_{l}$, $u_{u}$ are solutions of the equations ${\dot{u}_{l}=f_{l}(\tau,u_{l})}$, ${\dot{u}_{u}=f_{l}(\tau,u_{u})}$, respectively. Their initial conditions must match the initial conditions of (34) ${u_{1}(1)=u_{l}(1)=u_{u}(1)=u_{0}}$. These solutions provide bounds for the function $u_{1}(\tau)$ (Fig. 1), according to Chaplygin they are called “barrier” solutions. The following expressions complete the proof:

[TABLE]

[TABLE]

In addition, we give another proof that points to an interesting property of equation (34). The function $f(u_{1},\tau)$ in the considered range of $u_{1}$ decreases strictly monotonically, and

[TABLE]

These boundary properties do not allow the function $u_{1}(\tau)$ to go beyond the interval ${(0,1)}$. In addition, ${f(\tau,u_{m})<f(\tau,u_{n})}$ holds for any ${u_{m}>u_{n}}$ owing to the monotonicity of $f(\tau,u)$. Therefore,

[TABLE]

Equation (41) guarantees the convergence of any two solutions $u_{m}$ and $u_{n}$. Their difference is continuously decreasing and never changes sign because two different solutions of equation (34) cannot intersect according to Cauchy’s theorem on existence and uniqueness. \NAT@openFedoruk, 1985, p. 10; Khalil, 2002, p. 88\NAT@close. Thus, any solution of equation (34) starting from arbitrary time eventually lies inside a cylinder of non-zero radius and never goes out of it (Fig. 2). This property is called dissipativity, and systems that demonstrate it, for example the original system (23, 24), are called a dissipative systems \NAT@openDemidovich, 1967, p. 287; Khalil, 2002, p. 168\NAT@close. Fig. 2 shows that the convergence time of the solutions of the dissipative system to the cylinder-tube is large (of the order of 10 d for the typical SN IIP conditions). It turns out that if one takes collisional processes into account, the convergence time decreases to only thousands of seconds. Moreover, if one neglects the width of the cylinder-tube, it is possible to say that the system “forgets” the initial conditions! See Appendix C for details.

From (34) it follows that ${\lim_{\tau\to\infty}\dot{u}_{1}=0}$. The function $u_{1}(\tau)$ becomes constant ${u_{1}(\tau{=}\infty)}$ at large time $\tau{=}\infty$. Because $u_{l}(\tau)$ and $u_{u}(\tau)$ bound the function $u_{1}(\tau)$, it follows from (37, 38) that ${0<u_{1}(\tau{=}\infty)<1}$. From (33), it also follows that ${u_{2}\to 0}$ at ${\tau\to\infty}$. This means that the real normalized electron number density becomes constant ${0<u_{\mathrm{e}}(\tau{=}\infty)<1}$ at large time, regardless of the initial conditions. But by solving (27, 28) it can be shown that the equilibrium normalized electron number density approaches zero as

[TABLE]

Its clear that in the time-dependent case the envelope expands with a greater degree of ionization compared with the steady-state solution. A similar phenomenon is observed in atmospheric explosions \NAT@openRaizer, 1959; Zeldovich & Raizer, 2002\NAT@close and in the early Universe under cosmological conditions with a slowing down of the recombination of the primordial plasma (Zeldovich et al., 1969; Peebles, 1968; Kurt & Shakhvorostova, 2014). The number density of free electrons experiences a “freeze-out” in this case. Unlike the “freeze-out” effect in terrestrial atmospheric explosions, the effect of time-dependence in SNe remains important even when the temperatures of material and radiation are constant. This is true, for example, for the optically thin case ${s_{1}=2}$, ${s_{2}=2}$ under the free-streaming approximation.

The “freeze-out” effect can be seen in Fig. 4 and 4. The first figure presents a numerical calculation for the case of the outer, high-velocity layers far from the photosphere (Table 1 column 1), while the second one is for physical parameters typical for the near-photospheric layers of SNe IIP (Table 1 column 2). The calculations were carried out under the assumption that the system is initially in equilibrium.

The value of $u_{\mathrm{e}}(\tau{=}\infty)$ may be small in some cases. The deviation of true number densities from the equilibrium ones in Fig. 4 at large time is not significant. The solutions $u_{1}^{ss}$ and $u_{1}$ saturate to 1, while $u_{\mathrm{e}}$, $u_{\mathrm{e}}^{ss}$ are close to zero (which corresponds to almost full recombination). However, there is a significant increase of deviation at the ${t-t_{0}\sim 10}$ d (Fig. 4). For the outer, high-velocity layers far from the photosphere (Fig. 4) the deviation of number densities also increases over these days.

In the next section, we will investigate the factors that influence the evolution of the non-stationarity at early time, when the system is still close to equilibrium.

4 The system at early time

In Section 3, the manifestation of violation of the steady-state approximation in kinetics was proved at large time. However, its investigation by analytical methods is difficult because the Abel equation is not generally integrable by quadratures. In this section, we are interested in deviation from the steady-state level populations at early time. Numerical calculations for SN 1999em, which is a typical SN IIP, show that at the initial time $t_{0}$ the deviation of real number densities from the equilibrium ones is negligible (Potashov et al., 2017). In contrast to Section 3, where the qualitative picture was sufficient to solve the problem at large time, here more quantitative estimates are required. Hence, the problem will be solved semi-analytically using typical values of the coefficients of the system (see Appendix B).

Let us rewrite the system (23, 24) in vector notation.

[TABLE]

where $\mathbf{u}{=}{(u_{1},u_{\mathrm{e}})^{\mathrm{T}}}$ is the vector of normalized populations and $\mathrm{T}$ denotes the transpose of a row-vector into a column-vector. We construct a system of equations in deviations $\mathbf{x}^{ss}{=}{(x_{1}^{ss},x_{e}^{ss})^{\mathrm{T}}}$ of the true normalized populations $\mathbf{u}$ from the steady-state solution of the system (27, 28) $\mathbf{u}^{ss}{=}{(u_{1}^{ss},u_{\mathrm{e}}^{ss})^{\mathrm{T}}}$, that is, $\mathbf{x}^{ss}\equiv\mathbf{u}-\mathbf{u}^{ss}$. Therefore, in general, one can write

[TABLE]

and hence

[TABLE]

Because by definition ${\mathbf{f}(\tau,0+\mathbf{u}^{ss})=0}$, the value ${\mathbf{x}^{ss}=0}$ is an equilibrium point. Therefore, one can write

[TABLE]

where

[TABLE]

is a Jacobi matrix and $\mathbf{h}(\tau,\mathbf{x}^{ss})$ represents higher-order terms.

[TABLE]

Detailed analysis of $\mathbf{h}$ shows that these terms vanish rapidly relative to ${\|\mathbf{x}\|}$ as the origin is approached, and the relationship

[TABLE]

holds. Then the first-order approximation is accurate and defined as the linearization of the system (see, for example, Vidyasagar, 2002, p. 210). The final Cauchy problem can thus be formulated as

[TABLE]

where we assume the initial conditions to be zero.

The linear time-varying system (46) describes the behaviour of the system (43) in a small segment of initial time while deviations $\mathbf{x}^{ss}$ are small. We will be able to understand how the deviations of $\mathbf{x}^{ss}$ from the steady-state solution grow by analysing the linearized case.

The general solution of (46) in vector form is written as \NAT@openDemidovich, 1967, p. 77; Kuijstermans, 2003, p. 41\NAT@close

[TABLE]

where ${\mathbf{K}(\tau,\tau^{\prime})=\mathbf{X}(\tau)\mathbf{X}(\tau^{\prime})^{-1}}$ is a normalized fundamental matrix at point $\tau^{\prime}$, and $\mathbf{X}(\tau)$ is fundamental matrix. $\mathbf{K}(\tau,\tau^{\prime})$ is also called a matriciant or Cauchy matrix. Our goal is to obtain an expression for $\mathbf{K}(\tau,\tau^{\prime})$. To find the fundamental matrix $\mathbf{X}(\tau)$, it is convenient to introduce the concept of dynamic eigenvalues and eigenvectors \NAT@openWu, 1980; Neerhoff & van der Kloet, 2001; Kuijstermans, 2003, p. 47\NAT@close. If there exists a scalar function $\lambda(\tau)$ and a non-zero differentiable vector function $\mathbf{v}(\tau)$, so that they satisfy the following condition

[TABLE]

where $\mathbf{I}$ is an identity matrix, then $\lambda(\tau)$ is called a dynamic eigenvalue of matrix $\mathbf{J}(\tau)$ associated with a dynamic eigenvector $\mathbf{v}(\tau)$. Quasi-static $\lambda^{qs}(\tau)$ and $\mathbf{v}^{qs}(\tau)$ are the classical eigenvalues and eigenvectors obtained by solving equation (48) with the right-hand side equal to zero. Similar to the classical algebraic case, the number of distinct dynamic eigenvalues $\lambda_{i}(\tau)$ of the matrix is equal to or less than the size of the matrix. The number of all dynamic eigenvalues of $\mathbf{J}$ is equal to the size of the matrix.

Furthermore, we change the variables using the transformation ${\mathbf{x}^{ss}=\mathbf{L}(\tau)\mathbf{y}}$, where $\mathbf{y}=\mathbf{y}(\tau)$ is a vector of new unknown variables. The $\mathbf{J}$ matrix is converted to a new matrix,

[TABLE]

In the classical algebraic case for constant $\mathbf{J}$ and $\mathbf{L}$, $\mathbf{B}$ is Jordan matrix whenever $\mathbf{J}$ has repeated eigenvalues. In the case with time-varying matrices, one can choose an algebraic transformation $\mathbf{L}$ so that $\mathbf{B}$ transforms into a diagonal $\mathbf{\Lambda}$ matrix (Wu, 1980, Theorem 2). This transformation belongs to the class of Lyapunov transformations (Kuijstermans, 2003, p. 133). Because the coordinate transformation $\mathbf{L}$ preserves the dynamic eigenvalues (Wu, 1980; Neerhoff & van der Kloet, 2001), ${\mathbf{\Lambda(\tau)}=\mathrm{diag}[\lambda_{i}(\tau)]}$. The final expression for the fundamental matrix $\mathbf{X}(\tau)$ is

[TABLE]

where

[TABLE]

Thus, the problem of finding solutions of (46) is reduced to the search of its eigenvalues $\lambda_{i}(\tau)$ and the parameters for a Lyapunov transformation $\mathbf{L}$ that obtains the diagonal matrix ${\mathbf{\Lambda}}$.

One way to find $\mathbf{L}(\tau)$ has been proposed by Wu (1980); Van Der Kloet & Neerhoff (2000); Kuijstermans (2003, p. 137). This method is based on an iterative algorithm

[TABLE]

with the conditions

[TABLE]

where $\mathbf{I}$ is an identity matrix, $\mathbf{\bar{\Lambda}}_{j}$ is a diagonal matrix, and $\mathbf{Q}_{j}$ is a transformation matrix that consists of quasi-static eigenvectors of the matrix ${\mathbf{\bar{\Lambda}}_{j-1}-\mathbf{Q}_{j-1}^{-1}\dot{\mathbf{Q}}_{j-1}}$, calculated at every moment of time. Thus, at each step of the iteration, the matrix from the previous step is diagonalized taking into account the error ${\mathbf{Q}_{j-1}^{-1}\dot{\mathbf{Q}}_{j-1}}$. Van Der Kloet & Neerhoff (2000) have proved that the iterative process converges:

[TABLE]

Let us compare the norms of the diagonal matrix and the error at each step of the iteration. For parameters typical for the problem that we are solving (see Appendix B), the following inequality already holds for the first iterative step:

[TABLE]

Indeed, for the initial moment, the norms can be estimated as

[TABLE]

where $P$ is defined in (25).

The inequality is strengthened as ${\lVert\mathbf{\bar{\Lambda}}_{1}(\tau)\rVert\sim\mathcal{O}(\tau^{2})}$ with time ${\tau\to\infty}$ and ${\lVert\mathbf{Q}_{1}^{-1}(\tau)\dot{\mathbf{Q}}_{1}(\tau)\rVert\sim\mathcal{O}(1/\tau)}$. Therefore, we can conclude that in our case the following approximation is perfectly suitable:

[TABLE]

where $\mathbf{V}_{qs}{=}(\mathbf{v}_{1}^{qs},\mathbf{v}_{2}^{qs})$ is the matrix with quasi-static eigenvector columns $\mathbf{v}_{1}^{qs}$, $\mathbf{v}_{2}^{qs}$ of the matrix $\mathbf{J}$ corresponding to the quasi-static eigenvalues of the Jacobian matrix. Therefore, expression (49) for the fundamental matrix $\mathbf{X}$ can be rewritten as

[TABLE]

where $\gamma_{i}^{qs}(\tau)$ is calculated from (50) for $\lambda_{i}^{qs}$. The normalized fundamental matrix looks like

[TABLE]

The quasi-static eigenvalues and eigenvectors for the two-dimensional matrix $\mathbf{J}$ can be written explicitly as

[TABLE]

[TABLE]

Here, $\mathrm{tr}(\mathbf{J})$ is the trace of the Jacobi matrix and $\Delta$ is its determinant. It should be noted that $\mathrm{tr}(\mathbf{J})=\lambda_{1}^{qs}+\lambda_{2}^{qs}$ and $\Delta=\lambda_{1}^{qs}\lambda_{2}^{qs}$.

From (42), it follows that the function $u_{\mathrm{e}}^{ss}/\tau^{3}$ decreases as ${\sim\tau^{-(s_{1}+s_{2}+3)/2}}$, which is steeper than $\tau^{-s_{2}}$. Because of $P\gg R$ (see Appendix B), ${P/\tau^{s_{2}}\gg 2R\,u_{\mathrm{e}}^{ss}/\tau^{3}}$ for any time ${\tau\geqslant 1}$. Therefore, $\mathbf{J}[2,2]\simeq\mathbf{J}[2,1]$ and

[TABLE]

according to (45).

Consider now the first eigenvalue of the matrix $\mathbf{J}$ defined in (53). It is equal to

[TABLE]

which takes into account (45). The expression (56) implies that $\lambda_{1}^{qs}(\tau)$ is a negative time function. One can ensure that initially $|\lambda_{1}^{qs}(1)|\ll P$ (see Appendix B). As time increases, the module of the function $|\lambda_{1}^{qs}(\tau)|$ decreases monotonically faster than $u_{\mathrm{e}}^{ss}R/\tau^{3}\sim\tau^{-(s_{1}+s_{2}+3)/2}$, which is steeper than $\tau^{-s_{2}}$. This means that $|\lambda_{1}^{qs}(\tau)|\ll P/\tau^{s_{2}}$ for any ${\tau\geqslant 1}$. Thus, the first eigenvector is simplified to

[TABLE]

Consider now the behaviour of the function $\mathrm{tr}(\mathbf{J}(\tau))$. This is a negative function of $\tau$. For its module, the inequality $|\lambda_{1}^{qs}(1)|\ll P<|\mathrm{tr}(\mathbf{J}(1))|$ holds. Because $P\gg A\gg R$ (see Appendix B), the module $|\mathrm{tr}(\mathbf{J}(\tau))|$ initially decreases as $P/\tau^{s_{2}}$ with time not steeper than $|\lambda_{1}^{qs}(\tau)|$ and then increases as $A\tau^{2}$. Hence, $|\lambda_{1}^{qs}(\tau)|\ll|\mathrm{tr}(\mathbf{J}(\tau))|$ for any ${\tau\geqslant 1}$. And finally, from

[TABLE]

it follows that $|\lambda_{1}^{qs}(\tau)|\ll|\lambda_{2}^{qs}(\tau)|$ and $\lambda_{2}^{qs}(\tau)\simeq\mathrm{tr}(\mathbf{J}(\tau))$ for any ${\tau\geqslant 1}$. Thus, the second eigenvector can be rewritten as

[TABLE]

For small time, as long as $|\mathrm{tr}(\mathbf{J}(\tau))|$ is dominated by the photoionization term $P/\tau^{s_{2}}$, we have

[TABLE]

but at large time, when the spontaneous emission $A\,\tau^{2}$ dominates, we can write

[TABLE]

Because

[TABLE]

one has that ${\exp(\gamma_{2}^{qs}(\tau)-\gamma_{2}^{qs}(\tau^{\prime}))}$ is an exponentially decreasing function of ${(\tau-\tau^{\prime})}$, for any ${\tau\geqslant 1}$. One can estimate the normalized fundamental matrix by summarizing all the above observations. From (52), it follows that when $\tau^{\prime}$ and $\tau$ are nearly close,

[TABLE]

$\mathbf{K}(\tau,\tau^{\prime})\sim\mathbf{I}$. In the opposite case, we assume that ${\exp(\gamma_{2}^{qs}(\tau)-\gamma_{2}^{qs}(\tau^{\prime}))\simeq 0}$ and find

[TABLE]

Finally, we consider the integrand function in (47). If (60) is true then $\mathbf{K}(\tau,\tau^{\prime})\,\dot{\mathbf{u}}^{ss}(\tau^{\prime})\sim\dot{\mathbf{u}}^{ss}(\tau)$, and the contribution of this term can be neglected. When (60) is not satisfied, the estimate can be written as

[TABLE]

where

[TABLE]

The analysis of the behaviour of the function $\mathrm{tr}(\mathbf{J}(\tau))$ presented above implies that

[TABLE]

And because ${|(\dot{u}_{1}^{ss}(\tau)+\dot{u}_{\mathrm{e}}^{ss}(\tau))|\ll|\dot{u}_{1}^{ss}(\tau)|}$ for any ${\tau\geqslant 1}$, we obtain the final answer

[TABLE]

The solution (62) includes only the larger of the two eigenvalues $\lambda_{1}^{qs}(\tau)$, which is the smaller in absolute value because they are negative. The smaller the $|\lambda_{1}^{qs}(\tau)|$ module, the faster the growth of deviation from the steady state, and the greater the role of time dependence. Thus, the value of this eigenvalue is a key parameter influencing the time-dependent effect.

The linearization (46) of the system (43) is applicable only when the deviations ${\mathbf{x}^{ss}}$ are small. However, we have previously seen that the deviations can be significant. Fig. 5(a) presents the solution $x_{1}^{ss}(t)$ of the system (46) and its linearization $x_{1}^{ss,\mathrm{lin}}(t)$ derived from (62) for physical time for the case of outer, high-velocity layers far from the photosphere (Table 1 column 1). while Fig. 5(b) presents solutions of the same systems for the case of near-photospheric layers (Table 1 column 2). The strong nonlinearity of the solutions $x_{1}^{ss}(t)$ is clearly seen when the normalized number densities $u_{1}^{ss}$ and $u_{1}$ are saturated to 1, which corresponds to almost full recombination (see Fig. 4,4). Thus, the solution (62) describes the evolution of the system only at early time.

However, in the next section, we show that the final expression for $|\lambda_{1}^{qs}(\tau)|$ determines whether the kinetic system is in equilibrium at any moment of time.

5 Frozen system

The question remains wether it is possible to specify a simple way to check the importance of the time-dependence effect under given conditions. The solution of the system with “frozen” coefficients would help us to find the answer \NAT@openKhalil, 2002, p. 365; Vidyasagar, 2002, p. 248\NAT@close.

Let us rewrite the system of equations (23, 24) where all the variable coefficients are fixed at some moment $\tau_{1}$. In fact, a non-autonomous system have transformed into an autonomous one.

[TABLE]

[TABLE]

Unknown ${u_{1}^{a}(\tau_{a})}$, ${u_{\mathrm{e}}^{a}(\tau_{a})}$ are the functions of a new time variable $\tau_{a}$. Obviously, the values of expressions $u_{1}^{ss}(\tau_{1})$, $u_{\mathrm{e}}^{ss}(\tau_{1})$, are the solutions of the system (27, 28). The right-hand side of (63, 64) vanishes for those values.

Let us prove that for any solution of the system (63, 64) for $\tau_{a}\to\infty$ under initial conditions $0<u_{1,e}^{a}(\tau_{a}=1)\leqslant 1$ one obtains $u_{1,e}^{a}(\tau_{a})\to u_{1,e}^{ss}(\tau_{1})$.

The change of variables $x_{1}(\tau_{a})=u_{1}^{a}(\tau_{a})-u_{1}^{ss}(\tau_{1})$ and $x_{\mathrm{e}}(\tau_{a})=u_{\mathrm{e}}^{a}(\tau_{a})-u_{\mathrm{e}}^{ss}(\tau_{1})$ transforms the original system to the reduced one \NAT@openDemidovich, 1967, p. 234; Khalil, 2002, p. 147\NAT@close:

[TABLE]

[TABLE]

The Jacobi matrix of this system coincides with the matrix (45) taken at time $\tau_{1}$: ${\mathbf{J}=\mathbf{J}(\tau_{1})}$. From the analysis carried out in Section 4, it follows that the two eigenvalues $\lambda_{1}$, $\lambda_{2}$ of the Jacobi matrix $\mathbf{J}$ are negative and different. Consequently, by Lyapunov’s indirect method \NAT@openFedoruk, 1985, p. 289; Khalil, 2002, p. 133\NAT@close, it appears that the system (65, 66) is exponentially stable, and its solution is written as \NAT@openFedoruk, 1985, p. 61; Khalil, 2002, p. 37\NAT@close

[TABLE]

where $\mathbf{v}_{i}$ are the eigenvectors of the Jacobi matrix $\mathbf{J}$ and $C_{i}$ are arbitrary constants. Obviously, $\lim_{\tau\to\infty}\mathbf{x}(\tau_{a})=0$, which shows the convergence property of the system (63, 64) (Demidovich, 1967, p. 281). Thus, the search for the limit of solutions of (63, 64) provides another way to find the algebraic approximation of $u_{1}^{ss}$, $u_{e}^{ss}$.

The next question that arises is how long a period $\tau_{ss}$ is needed for the solution $u_{1}^{a}(\tau_{a})$ to obtain the stationary value $u_{1}^{ss}$. If this period is significantly longer than the time-scale of change of the envelope’s physical parameters, the system does not have time to relax to equilibrium. From (67), it follows that the period of equilibration is determined by the smaller of the two modules of the eigenvalues ${\tau_{ss}^{-1}=\mathrm{min}(|\lambda_{1}(\tau_{1})|,|\lambda_{2}(\tau_{1})|)}$. From Section 4, we know that ${|\lambda_{1}(\tau_{1})|\ll|\lambda_{2}(\tau_{1})|}$, and thus ${\tau_{ss}^{-1}=|\lambda_{1}(\tau_{1})|}$ holds. The relaxation time $\tau_{ss}$ for any $\tau_{1}$ from (56) is

[TABLE]

A simple criterion for checking the equilibrium of the system during the photospheric phase at any time arises from a comparison of the value of expression (68) with the typical time-scale of change of SN envelope’s physical parameters. The calculations show (Fig. 7) that in the case of an optically thick medium, for the outer, high-velocity layers (“out” label) and for the near-photospheric layers (“ph” label), the time ${t_{ss}=\tau_{ss}t_{0}}$ substantially exceeds the duration of the photospheric phase itself. This means that the system is always in the non-equilibrium state for the whole envelope. The time-dependence effect is more significant for the outer layers.

The value of $\tau_{ss}$ increases with time (Fig. 7), so at ${\tau\to\infty}$ the time-dependent effect will always be important, which has already been proved in Section 3. However, we have seen (Fig. 4) that, for example, in the near-photospheric layers the deviations of the number densities from the steady-state values can decrease when the atoms are almost fully recombined. It can be concluded that the normalized populations of the equilibrium and non-equilibrium systems are close only when these populations are saturated.

Let us determine the physical quantities that affect the relaxation time more than others. At the very beginning of the plateau stage, equation (68) for the dense near-photospheric layers can be estimated as

[TABLE]

where ${t_{\mathrm{rec}}}$ is the classical recombination time (Osterbrock & Ferland, 2006, p. 22). This formula is consistent with the estimates of the effective recombination time-scale \NAT@openUtrobin & Chugai, 2005, eq. 4; Dessart & Hillier, 2007, eq. 21\NAT@close.

However, a more general conclusion can be obtained from (68) for layers of any density located both near and far from the photosphere. As the constant $B$ is small (see Table 1, Appendix B), $\tau_{ss}$ in equation (68) can be expand as

[TABLE]

Here, $f(\tau,A,P,s_{2},R)$ is a function that depends only on the parameters $A$, $P$, $s_{2}$, $R$. We limit ourselves to calculations for the first half of the photospheric phase, when the photoionization term $P/\tau^{s_{2}}$ dominates the spontaneous emission $A\tau^{2}$, and an expression for $f$ can be obtained explicitly. Then (70) will be rewritten as

[TABLE]

where we used equations (17), (20), (21), (25). Let us rewrite this expression for the relaxation time using physical dimensional values:

[TABLE]

The minimum and maximum of the function $\tilde{f}$ for the domain of $N_{0}$, $T_{\mathrm{e}}$ (see Appendix B) at any given time $t\geqslant t_{0}$ differ by less than an order of magnitude. It can be concluded that the key factor that affects the value of the relaxation time for the whole system is the value and the rate of change of the external hard continuous radiation between the Lyman and Balmer ionization thresholds. Because usually ${s_{1}>s_{2}}$, the rate of change of the concentration of L$\alpha$ photons coming from outside plays the most important role for the existence of the time-dependence effect. If there are a lot of hard photons, and their concentration decreases very slowly, then there will be no time-dependence effect. Similar consideration can be made for the second half of the photospheric phase, and the conclusion will be the same.

6 Full kinetic system

The equilibrium criterion can be expanded. The reciprocal of the smallest by modulus of the Jacobi matrix eigenvalue can be calculated for a complete kinetic system that takes into account all kinds of processes.

Let us fix the external field of continuous radiation and find the time evolution of relaxation times for different kinetic systems.

A description of the modelling for various kinetic systems is given in Potashov et al. (2017). Calculations are carried out using the original code LEVELS and based on the model R450_M15_Ni004_E7 for SN 1999em (Baklanov et al., 2005; Potashov et al., 2017). The physical parameters for the calculation are taken from the Lagrangian layer with the velocity $v\approx 5.8\cdot 10^{3}\,\textnormal{km}\;\textnormal{s}^{-1}$. This layer is close to the photosphere at the beginning of the plateau epoch and its characteristics were used in the Table 1, column 2. Eigendecomposition of the system Jacobi matrix is done by the FEAST eigensolver package (Polizzi, 2009; Peter Tang & Polizzi, 2014; Kestyn et al., 2016).

We considered four systems presented in Table 2. System $M_{1}$ is the simplest one. It contains pure hydrogen with “two levels plus continuum” and does not include any collisional processes. This system has been considered analytically in previous sections. System $M_{2}$ is similar to $M_{1}$, but with collisional processes included. System $M_{3}$ is also a pure hydrogen one, with collisional processes, but the model of the hydrogen atom contains 10 levels. System $M_{4}$ is similar to $M_{3}$, but with metal admixtures. For model $M_{3}$, we use a modified Sobolev approximation that takes into account the absorption in the continuum in the region of a line (Hummer & Rybicki, 1985; Chugai, 1987; Grachev, 1988). For model $M_{4}$, the selective absorption of L$\alpha$ photons in metal lines (for example, Fe II) (Chugai, 1988a, 1998) and overlapping lines (Olson, 1982; Andronova, 1991; Bartunov & Mozgovoi, 1987) is taken into account. The metal abundances are taken from the R450_M15_Ni004_E7 model. Moreover, for model $M_{4}$ we do not assume a detailed balance of the ground level with the continuum, and calculate the diffusive average intensity of the Lyman continuum in a self-consistent manner (Utrobin & Chugai, 2005, eq. 5).

Fig. 7 shows the results of the calculation. The thin solid line “ph” corresponds to the analytical estimation obtained with the formula (68). This is the same line as that one designated by “ph” in Fig. 7. Numerical calculation of the model $M_{1}$ (thick dashed line) agrees well with it. The thick solid curve $M_{2}$ deviates from $M_{1}$ during the first 10 d after $t_{0}$. Collisional processes decrease the relaxation time only at the very beginning of the photospheric phase. Additional levels of $M_{3}$ (dotted line) slightly decrease the relaxation time. Taking metals into account in $M_{4}$ (bold dashed line) leads to a significant decrease of the relaxation time and hence to a weakening of the time-dependence effect. Because the relaxation time is still much longer than the duration of the photospheric phase, this is not enough to completely cancel the non-stationarity. This result was already obtained earlier in Potashov et al. (2017), where it was shown that the increase of the metal concentration in the envelope leads to a weakening of the time-dependence effect in kinetics when the continuous radiation is fixed. This allows us to say that the metal admixtures play the most important role in the non-stationarity among other effects studied in this section.

Summarizing the statements above with the findings of Section 5, we can conclude the following. The concentration of L$\alpha$ photons that come from outside, as well as the change of this concentration and the rate of the selective absorption of L$\alpha$ photons in metal lines are the most significant factors affecting the non-stationarity effect. The greater the number of external L$\alpha$ photons and the faster their absorption by metals, the weaker the effect of time-dependent ionization.

The metallicity of the envelope strongly affects the flux in the U band, and this complicates the situation. The emergent UV flux is significantly enhanced at low metallicity, owing to the reduction of line blanketing (Goldstein & Blinnikov in preparation; Dessart et al., 2013). On one hand, this means that the low metallicity leads to an increase of the number of external L$\alpha$ photons and, consequently, decreases the relaxation time of the system. On the other hand, there is an opposite effect of the relaxation time increase owing to the reduced number of absorbers. Therefore, the final answer to the question of whether the time-dependent ionization effect is important requires checking the equilibrium criterion described here. At least for SN 1999em, our calculation shows the importance of accounting for the time-dependence in kinetics (see Fig. 7, model $M_{4}$).

7 Discussion and Conclusions

There are two characteristic time-scales in the SN ejecta. One defines the time of relaxation to equilibrium conditions. This time equals the reciprocal of the smallest by modulus of Jacobi matrix eigenvalue for a complete kinetic system that takes into account all kinds of processes at any time (not only during the photospheric phase). The second one is the typical time-scale for the change of SN envelope parameters.

We come to the key conclusion of the paper. If the ratio of the two time-scales is small the system is in equilibrium. In this case one can investigate the steady-state algebraic approximation instead of the real time-dependent system. If this criterion of equilibrium is not satisfied, then a time-dependent effect is important, and the real atomic level populations deviate from the equilibrium values. In this case, the replacement of the original time-dependent system with the algebraic steady-state one is incorrect.

Acknowledgements

The authors are grateful to S.I. Blinnikov, V.P. Utrobin and P.V. Baklanov for useful discussions. The authors are also thankful to E.I. Sorokina for thorough reading the manuscript before publishing. M.Sh. Potashov thanks S.G. Moiseenko and G.S. Bisnovaty-Kogan for repeated invitations to the school-seminar organized by them in the town of Tarusa, where the results of this paper were presented. We are grateful to the anonymous referee for valuable comments. The work of MP was partially supported by the grant of the Russian Foundation for Basic Research (RFBR), grant 19–02–00567; the numerical calculations in the Section 6 were funded by the Russian Science Foundation, grant 19–12–00229. The work of AY was supported by the RFBR, grant 18–29–21019.

Appendix A Detailed balance of the ground level population of hydrogen

Let us consider the photo-processes that influence the electronic population of the ground level of the hydrogen atom. We assume that there are admixtures of other elements in the envelope. Let us introduce two assumptions. First, the main contribution to the opacity in the frequency band of the Lyman continuum ${\nu\geqslant\nu_{LyC}}$ is due to the free-bound processes and caused mainly by neutral hydrogen. Secondly, the optical depth in this band is very large.

Let us write down the coefficients of emission and true absorption, corrected for induced emission using the Einstein-Milne relations for continua (Mihalas, 1978; Hubeny & Mihalas, 2014):

[TABLE]

[TABLE]

Here $\alpha_{1c}(\nu)$ is the photoionization cross section from the first level, $B_{\nu}(T_{\mathrm{e}})$ is the intensity of blackbody radiation for the electron temperature $T_{\mathrm{e}}$, and

[TABLE]

is the local thermodynamic equilibrium (LTE) value of $n_{1}$ computed with the Saha equation using the actual values of the electron and proton number densities $n_{\mathrm{e}}$ and $n_{\mathrm{p}}$.

For an optically thick medium, the transport equation in the continuum is solved as $J_{c}(\nu,t)=S_{c}(\nu,t)$, where $S_{c}(\nu,t)$ is a source function, defined as

[TABLE]

Substituting (73) and (74) into (76) we have

[TABLE]

In accordance with Mihalas (1978) and Hubeny & Mihalas (2014, p. 273), the photoionization rate from the ground level is the integral

[TABLE]

The photorecombination rate looks like

[TABLE]

Substituting equation (77) into (78) and (79), and taking into account (75) gives $n_{1}P_{1c}(t)=n_{\mathrm{e}}n_{\mathrm{p}}R_{c1}(t)$. Thus, in the general case, with admixtures, the photoionization rate from the ground level of hydrogen and the radiative recombination rate to the ground level completely coincide and the ground level of hydrogen is in detailed balance with the continuum.

Appendix B Typical values

Let us write down the typical values of the physical parameters of the SN envelope and its radiation, based on calculations for SN 1999em (Baklanov et al., 2005; Potashov et al., 2017), which is a standard SN IIP. The computations were performed with the radiation-hydrodynamical code STELLA. We will mark values that correspond to the near-photospheric layers with “$\mathrm{ph}$” and to the outer, high-velocity layers with “$\mathrm{out}$”.

The initial time is ${t_{0}\approx 20}$ d after the explosion in the core.

The two-photon ${2\mathrm{s}\to 1\mathrm{s}}$ decay rate is taken from Nussbaumer & Schmutz (1984). We consider two $l$-sublevels ${2\mathrm{s}}$ and ${2\,\mathrm{p}}$ as a single super-level 2 (Hubeny & Lanz, 1995) and assume that the populations of the sublevels are proportional to their statistical weights ${g_{2\mathrm{s}}}$, ${g_{2\mathrm{p}}}$. This leads to

[TABLE]

The initial number densities of gas are ${N_{0}(\mathrm{ph})\approx 10^{11}\,\mathrm{cm}^{-3}}$ and ${N_{0}(\mathrm{out})\approx 10^{8}\,\mathrm{cm}^{-3}}$. It follows that for the constant $A$ defined in (17, 25) we have corresponding values

[TABLE]

Typical values of the intensity of the continuum radiation at the initial moment at the L$\alpha$ frequency is taken from the calculation of STELLA:

[TABLE]

Thus, for $B$ defined in (17, 25) we have

[TABLE]

The photoionization coefficient from the second level $P_{2c}$ at the initial moment is also taken from STELLA calculations:

[TABLE]

Then we get from (25):

[TABLE]

The effective photorecombination coefficient to the second level at the initial moment can be estimated from Hummer (1994, Table 1) using the material temperatures typical for a supernova envelope ${T_{\mathrm{e}}\approx 3000-7000\,\mathrm{K}}$ as ${R_{2c}\approx 3\cdot 10^{-13}\,\mathrm{\frac{cm^{3}}{s}}}$. Then we will have for $R$ from (17, 25):

[TABLE]

Values for the power exponents taken from STELLA calculations are

[TABLE]

Because in the optically thin atmosphere ${s_{1}=2}$ and ${s_{2}=2}$, in general case, the domains of these exponents are ${s_{1}\geqslant 2}$ and ${s_{2}\geqslant 2}$.

Appendix C The tube

In general case, we limit the domain of these exponents to ${s_{1}\geqslant 2}$ and ${s_{2}\geqslant 2}$.

Let us prove the dissipativity of the system (23, 24) in another way. Suppose that we know some particular bound solution $0{<}\tilde{u}_{1}{\leqslant}1$ and $0{\leqslant}\tilde{u}_{\mathrm{e}}{\leqslant}1$ (unperturbed motion) of a nonautonomous nonlinear differential system (23, 24). The change of variables ${x_{1}=u_{1}-\tilde{u}_{1}}$ and ${x_{\mathrm{e}}=u_{\mathrm{e}}-\tilde{u}_{\mathrm{e}}}$ transforms the system into the form of perturbed motion system \NAT@openDemidovich, 1967, p. 234; Khalil, 2002, p. 147\NAT@close:

[TABLE]

[TABLE]

The trivial solution ${x_{1}=0}$, ${x_{\mathrm{e}}=0}$ is an equilibrium point of the transformed system for any $\tau\geqslant 1$.

Consider the Lyapunov function candidate

[TABLE]

The function $V$ is positively defined at $d>0$ and it is an elliptic paraboloid. The derivative of $V$ along the trajectories of the system (80, 81) is given by

[TABLE]

where $\dot{x}_{1}$ and $\dot{x}_{\mathrm{e}}$ are (80) and (81), respectively. At large time and at small $x_{1}$ and $x_{\mathrm{e}}$ it expresses as

[TABLE]

This function is negative semidefinite in the domain ${\Omega=(0<\tilde{u}_{1}{+}x_{1}\leqslant 1,0<\tilde{u}_{\mathrm{e}}{+}x_{\mathrm{e}}\leqslant 1)}$.

From the Lyapunov direct method and (82, 83) it follows that the trivial solution of system (80, 81) is stable for small perturbations of the initial conditions for any positive $d$. The stability problem for this system can be solved in this way on the half-axis ${\tau>\tau_{1}}$ with $\tau_{1}\geqslant 1$. The stability on a half-axis ${\tau>1}$ is obtained taking into account the theorem of continuous dependence on the parameter (Khalil, 2002, p. 95) for the solution on a finite interval ${1\leqslant\tau\leqslant\tau_{1}}$.

The analysis of the function (82) reveals that the system (80, 81) is dissipative. Indeed, Section 3 shows that ${\tilde{u}_{1}\approx 1-\tilde{u}_{\mathrm{e}}}$ at large time. Then (83) can be rewritten as

[TABLE]

This function is a parabolic cylinder with the zero values along the line $l$ defined by the equation ${x_{1}=-x_{\mathrm{e}}}$. This means that one can specify a value $\mu$ so that the set $l$ does not belong to the $\tilde{\Omega}$ domain defined as the difference between the $\Omega$ domain and the domain specified by the inequality

[TABLE]

for all moments of time ${\tau\geqslant 1}$. Thus, the function $\dot{V}(\tau,x_{1},x_{\mathrm{e}})$ is negative definite in $\tilde{\Omega}$. Therefore, according to the theorem of Yoshizawa \NAT@openDemidovich, 1967, p. 290; Kuntsevich & Lychak, 1977, p. 47-48\NAT@close, the perturbed motion of the system is dissipative. Hence, any solution of the equation, starting from arbitrary time, eventually will lie inside a cylinder of a non-zero radius and will never go out of it.

The given proof of dissipativity assumes small perturbations $x_{1}$ and $x_{\mathrm{e}}$. Let us prove by contradiction that dissipativity is true for any $x_{1}$, $x_{\mathrm{e}}$ in $\Omega$. Note additionally that the solutions of the initial value problem (23, 24, 26) are bounded. If one suggests the presence of two cylinders then one inevitably obtains two close-lying bounded solutions that belong to different cylinders. Contra this, we showed earlier that all close-lying solutions must lie in one cylinder owing to dissipativity. This contradiction yields the proof of uniqueness of the cylinder-tube.

The theorem of Yoshizawa does not give practically useful estimates of the convergence time of the solutions of a dissipative system into a semi-infinite cylinder-tube. By restricting our analysis to the numerical estimates of this time, we can solve the system (80, 81) with the initial conditions ${x_{1}(1)\approx\tilde{u}_{1}(1)}$, ${x_{\mathrm{e}}(1)=1-\tilde{u}_{\mathrm{e}}(1)}$. It was shown in Section 3 that the initial condition of $u_{2}$ does not affect the solution. Because ${u_{2}=1-u_{1}-u_{\mathrm{e}}}$, then from (31) it follows that the expression ${-(x_{1}+x_{\mathrm{e}})}$ decreases exponentially as $e^{-G_{3}(\tau)}$. Therefore, without a loss of generality, only the behaviour of $x_{1}$ is considered below. Fig. 8 shows that the convergence time of the solutions of the dissipative system to the cylinder-tube is large for an optically thick atmosphere and physical parameters typical for the near-photospheric layers of SN IIP and equal to about ${\sim 10}$ d. The situation changes dramatically if one takes into account collisional processes. For example, in Fig. 8 we present the results of calculation where collisional processes have been taken into account under the same conditions. The convergence time in this case is only thousands of seconds. If one neglects the width of the cylinder-tube (${x_{1}^{\mathrm{col}}\ll 1}$ at time greater than 1 d) one can say that the system “forgets” the initial conditions!