A simple model of time-dependent ionization in Type IIP supernova envelopes
M. Sh. Potashov A. V. Yudin


Multiple-copy state discrimination of noisy qubits
    Sarah Croke
    Received: date / Accepted: date
==================================================




  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.









§ 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)
  <cit.>,
  the spectral-fitting expanding atmosphere method (SEAM)
  <cit.>,
  and the dense shell dethod (DSM)
  <cit.>.
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)
  <cit.>.

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 <cit.>
  to explain the high Hα intensity in the spectra of SN 1970G,
  and by
  <cit.>
  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.

<cit.>
  found a strong time-dependent effect in the ionization kinetics
  and hydrogen lines in SNe IIP during the photospheric phase.
In the next paper,
  <cit.>
  also took into account the time-dependent effect in the energy equation.
In these papers, it was shown that
  the Hα 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 into the outer high-velocity layers.
<cit.> 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.
<cit.> still applied the steady-state approach
  that was implemented in the CMFGEN package,
  but the Hα line in the hydrogen-rich envelopes
  was weaker than the observed one during the recombination epoch.
In particular, the model did not reproduce Hα 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
  <cit.>,
  and for the latest version, that in the radiative transfer
  <cit.>.
This strengthened the Hα in the
  simulated spectrum and led to better agreement with observations.

On the other hand, based on the computations
  with the PHOENIX software package,
  <cit.>
  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.

<cit.>,
  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.

<cit.>
  showed the importance of the non-stationary kinetics
  for SN 1999em in the purely hydrogen case
  using the codes STELLA
  <cit.>
  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 <ref> describes the
  physical model of the problem.
In Section <ref>, it is concluded
  that the time-dependent effect is significant at least
  in the large time limit.
In Section <ref>, 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 <ref>,
  it is shown that the expression described in Section <ref> 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 <ref> 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.






§ 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
  2s and 2p
  are proportional to their statistical weights.
Thus, the second level is considered as a single so-called super-level
  <cit.>.




























The light-curve behaviour of a typical SN IIP
  can be divided into various characteristic stages
  <cit.>:

  
  * shock breakout;
  
  * adiabatic cooling phase;
  
  * photospheric phase (cooling and recombination
 wave);
  
  * phase of radiative diffusion cooling;
  
  * exhaustion of radiation energy;
  
  * plateau tail;
  
  * radioactive tail ().




















































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
  <cit.>
  this phase lasts from t_0 ∼ 20 d to ∼ 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 Rph,
  the radiation temperature T_c
  and the gas temperature Te
  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
  <cit.>,
  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 <ref>
  (see also Appendix <ref>).

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


    ∂ρ/∂ t = -▽(ρ v),


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


    Dρ/D t = -ρ(▽· v).



In the free homologous expansion stage (v ∝ r),
  equation (<ref>) is simplified:


    Dρ/D t + 3ρ/t = 0.


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


    Dn_i/D t + 3n_i/t = K_i(t),


  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


    K_1(t) = (N(t) - n_1 - ne) (A_2q + A_21) + n_1B_12J_12(t),


Here

    N(t) = N_0(t_0/t)^3

  is the hydrogen number density,
  and A_2q is the two-photon decay probability 2s→ 1s.
The reverse 1s→ 2s
  transition (two-photon absorption) rate
  is much lower than the 2s→ 1s rate,
  and we neglect this process
  <cit.>;
  we also neglect stimulated two-photon decays.
A_21 and B_12 are the Einstein coefficients for spontaneous and induced
  transition 1 ↔ 2,
  and J_12(t) is the mean intensity of radiation for the 2 → 1 transition
  averaged over the line profile.

Equation (<ref>) with the index
  e instead of i is valid
  for the number density of free electrons,
  with


    Ke(t) = (N(t) - n_1 - ne) P_2c(t) - ne^2 R_c2(t) .


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
  ν⩾ν_LyC
  is provided by free-bound processes
  <cit.>.
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 <ref>)
  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
  (<ref> and <ref>).

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

    n_2 = N(t) - n_1 - ne≪ n_1,

  and write down the standard formulas for the Sobolev approximation
  <cit.>.


The Sobolev optical depth in Lα is


    τ_S(t) ≃c^3/8πA_21/ν_Lα^3g_2/g_1n_1t.


The frequency-averaged mean intensity of the transition 1 ↔ 2 is


    J_12(t) = (1 - β(t))   S(t) + β(t)   J_c(ν_Lα, t),


  where J_c(ν_Lα,t) is the average intensity of the continuum
  on the Lα frequency.

Because the medium is optically thick for Lα photons, τ_S(t) ≫ 1.
Then the escape probability of Lα photons,
  summed over all directions and line frequencies, is


    β(t) = 1-e^-τ_S(t)/τ_S(t)≃1/τ_S(t).


The source function is


    S(t) ≃2 h ν_Lα^3/c^2( g_1 n_2/g_2 n_1).


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α
  (a τ_S ≫ 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
  <cit.>.
However, <cit.> and <cit.>
  have shown that the estimation (<ref>) remains true
  for the conservative scattering of Lα 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.
<cit.> showed
  that the frequency-weighted mean intensity J_12(t) defined in (<ref>) is
  weakly dependent on the mechanism of frequencies redistribution.
Expressions (<ref>) and (<ref>) remain correct
  in the case of non-conservative scattering with partial frequencies redistribution
  and partial non-coherence owing to the Stark effect
  <cit.>.

When absorption in the continuum in the region of a line is significant
  (for example, Lα photons can ionize Ca II from the second level),
  one must take into account corrections to the Sobolev approximation
  <cit.>.
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α
  <cit.>.
These processes increase the value of the effective escape probability
  (<ref>) (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
  <cit.>.
In the current analytical approach, we use a simple model and
  do not take into account these processes,
  using hereafter equations (<ref>) and (<ref>).
However, in Section <ref>,
  these modifications of the Sobolev approximation
  will be taken into account numerically.

Combining
(<ref>–<ref>,
  <ref>–<ref>),
  we obtain the system


    ṅ_1 =
        (N(t) - n_1 - ne) (A_2q + A_21β(t))
        -
        n_1B_12β(t) J_c(ν_Lα, t) - 3n_1/t,



    ṅe = (N(t) - n_1 - ne) P_2c(t) - ne^2 R_c2(t) - 3ne/t.


In accordance with
  <cit.>
  and
  <cit.>,
  the total photoionization coefficient for the second level is the integral


    P_2c(t) = 4π∫_ν_2^∞ J_c(ν, t)α_2c(ν)/hνdν,

  where ν_2 is a frequency of the 2 ↔ c transition,
  c means the continuum,
  and α_2c is a photoionization cross-section of the second level
  at the frequency ν.
The total radiative recombination coefficient
  of the second level for a purely hydrogen plasma
  looks like


    R_c2
        = 4π Φ_Saha(Te) ∫_ν_2^∞α_2c(ν)/hν2hν^3/c^2e^-hν/kTedν≈64 π^5 m e^10/3 √(3) c^3 h^6 g_ II(2,ν_2) Φ_Saha(Te) E_1(hν_2/kTe)

  if the processes of induced emission are neglected.




Here
  Φ_Saha(Te) is the Saha-Boltzmann factor;
  g_ II(2,ν_2) is a bound-free Gaunt factor for the 2 ↔ c transition;
  and E_1 is the exponential integral.
Note that R_c2 is constant over time if Te 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 (<ref>).
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
  <cit.>
  for the recombination coefficient calculation.
In the opposite Case A
  <cit.>,
  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
  <cit.>.
We choose the average value of R_c2 from <cit.>
  for our simple system.

Let us now introduce the dimensionless variables


    u_1(t) = n_1/N(t) = n_1/N_0t^3/t_0^3,
          ue(t) = ne/N(t) = ne/N_0t^3/t_0^3


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


    u̇_1 = (1-u_1-ue) [A_2q + Ã/u_1(t/t_0)^2
        ]
        -
        B̃ J_c(ν_Lα, t) (t/t_0)^2,



    u̇e = (1-u_1-ue) P_2c(t)-ue^2R̃(t_0/t)^3.


Here we also introduce a new notation for constants:


    Ã=8 πν_Lα^3/c^3g_1/g_21/N_0 t_0,  B̃=4π/h c1/N_0 t_0,  R̃= N_0 R_c2.


Especially important for further simplification of the system
  (<ref>, <ref>)
  is an investigation of the behaviour of J_c(ν_Lα, 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),


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


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


    W(t) ≃1/4(Rph/Vt)^2.


Then the continuum intensity J_c(ν_Lα, t)
  and the photoionization coefficient P_2c(t)
  decrease as ∼ 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,


    J_c(ν_Lα, t) ≃ J_c(ν_Lα, t_0) (t_0/t)^s_1,



    P_2c(t) ≃ P_2c(t_0) (t_0/t)^s_2 = P̃(t_0/t)^s_2.


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 ⩾ 2 and s_2 ⩾ 2.
In the optically thick case, the
  typical values of these exponents can be significant
  (see Table <ref>, Appendix <ref>).

Now let us introduce the normalized time


    τ = t/t_0.


Taking into account
   (<ref>), (<ref>) and (<ref>),
  we can rewrite the system (<ref>, <ref>) as


    u̇_1 = (1-u_1-ue)(Q + A/u_1τ^2)
      -B/τ^s_1-2,



    u̇e = (1-u_1-ue)P/τ^s_2
        -ue^2R/τ^3.


  Now u_1 and ue are functions of τ, and

    Q = A_2q t_0,  
      A = Ã t_0,    
      B = B̃J_c(ν_Lα, t_0) t_0,    
      P = P̃ t_0,    
      R = R̃ t_0.


Table <ref> 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 <ref> for details).

The initial conditions for the problem are:


    0 < u_1(1) ⩽ 1,
          
        0 ⩽ue(1) ⩽ 1,
       
        u_2(1) = 1 - u_1(1) - ue(1) ≪ u_1(1)
          (<ref>).







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


    (1-u_1^ss-ue^ss)(Q + A/u_1^ssτ^2)
      -B/τ^s_1-2 = 0,



    (1-u_1^ss-ue^ss)P/τ^s_2 - (ue^ss)^2R/τ^3 = 0,


  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
  (<ref>, <ref>)
  and (<ref>, <ref>).
In the next section, it will be shown that at large time
  the deviation of u_1, ue from u_1^ss, ue^ss
  persists regardless of the initial conditions.






§ THE SYSTEM AT LARGE TIME

  

We will investigate the solutions of systems
  (<ref>, <ref>) and
  (<ref>, <ref>)
  at large time τ→∞.
In this case, we will assume that the photospheric phase lasts for an indefinitely
  long time.

Let us transform the original system
  (<ref>, <ref>)
  by introducing the function u_2
  according to the relationship u_1+u_2+ue=1.
We should keep only terms linear in u_2
  because the smallness of this function was already assumed in the derivation
  (see <ref> and <ref>).
Thus, this is a reduction of the source system to a “normal” form:


    u̇_1   =u_2 g_1(τ,u_1)-g_2(τ),
    u̇_2   =-u_2 g_3(τ,u_1)+g_4(τ,u_1),


  where we have introduced the following functions:


    g_1(τ,u_1) = Q + A/u_1τ^2,
       g_2(τ) = B/τ^s_1-2,
       g_3(τ,u_1) = Q + A/u_1τ^2+P/τ^s_2
        +2 (1-u_1)R/τ^3,
       g_4(τ,u_1) = B/τ^s_1-2
        +(1-u_1)^2R/τ^3.


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


    u_2(τ) = 
        e^-G_3(τ)[u_2(1)+∫_1^τg_4(τ',u_1(τ'))e^G_3(τ')dτ'],


  where


    G_3(τ)≡∫_1^τg_3(τ',u_1(τ'))dτ'.


Equation (<ref>) together with (<ref>)
  represents a formal solution for u_2, because u_1 is still unknown.

Analysis of the system (<ref>, <ref>)
  shows that the solutions that satisfy the initial conditions
  0 < u_1(1) ⩽ 1,
  0 ⩽ u_2(1) ⩽ 1
  remain bounded
  and stable in the sense of Lyapunov
  Demidovich1967book, p. 66; Khalil2002book, p. 111,
  see also Appendix <ref>.
The functions g_1(t) and g_3(t)
  increase rapidly,
  g_1,3∼𝒪(τ^2),
  if u_1 is bounded and finite.
In contrast, the functions g_2,4 are decreasing
  as some power of τ.
We can use the fact that
  for any bounded behaviour of u_1(τ),
  G_3 is a rapidly increasing function of time,
  at least as ∼𝒪(τ^3).
Therefore, terms with exp(G_3) are most important at large time.
The term with u_2(1) in (<ref>)
  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 τ→∞, we obtain


    u_2≃g_4(τ,u_1)/g_3(τ,u_1)-1/g_3(τ,u_1)d/dτ(g_4(τ,u_1)/g_3(τ,u_1))+…


The first term of this decomposition is actually a steady-state approximation
  of equation (<ref>).
Substituting it into (<ref>), we obtain the equation for u_1:


    u̇_1 = f(τ,u_1) = g_1(τ,u_1)g_4(τ,u_1)/g_3(τ,u_1)
        -g_2(τ).


Equation (<ref>)
  is a special case of the Appell equation <cit.>.
That is, it is a generalized Abel equation of the second kind
  <cit.>
  and unfortunately it is generally not integrable by quadratures.

Usually, approximate analytical methods, such as
  the small parameter method <cit.>,
  Chaplygin method <cit.>
  or power series method,
  are used to solve nonlinear differential equations.
Case (<ref>) 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 (<ref>).
It is sufficient to prove that any solution of this equation
  that satisfies the initial conditions 0 < u_1(1) ⩽ 1
  is bounded in the interval (0, 1).
To do this, consider two functions


    f_l(τ,u) = - |u| P B/A τ^s_1+s_2,



    f_u(τ,u) = (1-u)^2R/τ^3.


Let us take an arbitrary moment of time τ_1 ≫ 1.
In the domain (τ∈ [1, τ_1]; u ∈ℝ),
  the inequalities f_l(τ,u) ⩽ f(τ,u) ⩽ f_u(τ,u) are obviously satisfied.
Because the function f(τ,u) is continuously differentiable
  in this domain and satisfies the Lipschitz condition,
  according to the Chaplygin theorem on differential inequalities
  BerezinZhidkov1959book, p. 260; Khalil2002book, p. 102,
  one can write u_l(τ) ⩽ u_1(τ) ⩽ u_u(τ)
  for any time in [1, τ_1],
  where u_l, u_u are solutions of the equations
  u̇_l = f_l(τ,u_l),
  u̇_u = f_l(τ,u_u), respectively.
Their initial conditions must match the initial conditions of
  (<ref>)
  u_1(1) = u_l(1) = u_u(1) = u_0.
These solutions provide bounds for the function u_1(τ)
  (Fig. <ref>),
  according to Chaplygin they are called “barrier” solutions.
The following expressions complete the proof:


    u_l(τ) = u_0 exp[-B P (1-τ^-(s_1+s_2-1))/A (s_1+s_2-1)] > 0,




    u_u(τ) = 1-
        [
        1/1-u_0-R/2(1-τ^-2)
        ]^-1 < 1.







In addition, we give another proof that points
  to an interesting property of equation (<ref>).
The function f(u_1,τ) in the considered range of u_1
  decreases strictly monotonically, and


    f(τ,0)   =R/τ^3 > 0, 
    
      f(τ,1)   =-B P τ^-(s_1+s_2-2)/Q + A τ^2+P/τ^s_2 < 0.


These boundary properties do not allow the function u_1(τ)
  to go beyond the interval (0, 1).
In addition, f(τ,u_m)<f(τ,u_n)
  holds for any u_m>u_n
  owing to the monotonicity of f(τ,u).
Therefore,


    d(u_m-u_n)/d τ = f(τ,u_m) - f(τ,u_n) < 0.


Equation (<ref>) 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 (<ref>)
  cannot intersect according to Cauchy's theorem on existence and uniqueness.
  Fedoruk1985book, p. 10; Khalil2002book, p. 88.
Thus, any solution of equation (<ref>)
  starting from arbitrary time
  eventually lies inside a cylinder of non-zero radius
  and never goes out of it (Fig. <ref>).
This property is called dissipativity, and systems that
  demonstrate it, for example the original system
  (<ref>, <ref>),
  are called a dissipative systems
  Demidovich1967book, p. 287; Khalil2002book, p. 168.
Fig. <ref> 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 <ref> for details.





From (<ref>) it follows that
  lim_τ→∞u̇_1 = 0.
The function u_1(τ) becomes constant
  u_1(τ=∞) at large time τ=∞.
Because u_l(τ) and u_u(τ) bound the function u_1(τ),
  it follows from (<ref>, <ref>)
  that 0 < u_1(τ=∞) < 1.
From (<ref>), it also follows
  that u_2 → 0 at τ→∞.
This means that the real normalized electron number density
  becomes constant 0 <ue(τ=∞) < 1
  at large time, regardless of the initial conditions.
But by solving (<ref>, <ref>)
  it can be shown that the equilibrium normalized electron number density
  approaches zero as


    ue^ss∼τ^-(s_1 + s_2 - 3)/2.


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
  Raizer1959; ZeldovichRajzer2008book and
  in the early Universe under cosmological conditions
  with a slowing down of the recombination of the primordial plasma
  <cit.>.
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. <ref> and <ref>.
The first figure presents a numerical calculation for the case of the outer, high-velocity layers
  far from the photosphere (Table <ref> column 1),
  while the second one is for physical parameters typical
  for the near-photospheric layers of SNe IIP
  (Table <ref> column 2).
The calculations were carried out
  under the assumption that the system is initially in equilibrium.

The value of ue(τ=∞) may be small in some cases.
The deviation of true number densities from the equilibrium ones
  in Fig. <ref> at large time is not significant.
The solutions u_1^ss and u_1 saturate to 1,
  while ue, ue^ss are close to zero
  (which corresponds to almost full recombination).
However, there is a significant increase of deviation
  at the t-t_0∼ 10 d
  (Fig. <ref>).
For the outer, high-velocity layers far from the photosphere
  (Fig. <ref>)
  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.










§ THE SYSTEM AT EARLY TIME

  

In Section <ref>,
  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
  <cit.>.
In contrast to Section <ref>, 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 <ref>).

Let us rewrite the system (<ref>, <ref>) in vector notation.


    𝐮̇ = 𝐟(τ,𝐮),


  where 𝐮=(u_1,ue)^T is the vector of normalized populations
  and T denotes the transpose of a row-vector into a column-vector.
We construct a system of equations in deviations
  𝐱^ss=(x_1^ss,x_e^ss)^T
  of the true normalized populations 𝐮
  from the steady-state solution of the system
  (<ref>, <ref>)
  𝐮^ss=(u_1^ss,ue^ss)^T,
  that is, 𝐱^ss≡𝐮 - 𝐮^ss.
Therefore, in general, one can write


    d/dτ(𝐱^ss+𝐮^ss) = 𝐟(τ,𝐱^ss+𝐮^ss),


and hence


    𝐱̇^ss = 𝐟(τ,𝐱^ss+𝐮^ss) - 𝐮̇^ss.


Because by definition 𝐟(τ,0+𝐮^ss) = 0,
  the value 𝐱^ss = 0 is an equilibrium point.
Therefore, one can write


    𝐱̇^ss = 𝐉(τ) 𝐱^ss + 𝐡(τ,𝐱^ss) - 𝐮̇^ss,


  where
  
    𝐉(τ) = .∂ (f_1,fe)/∂ (u_1^ss,ue^ss)|_u_i=u_i^ss

  is a Jacobi matrix and
  𝐡(τ,𝐱^ss) represents higher-order terms.


    𝐉(τ) =
      [ -Q - A/u_1^ss(1-ue^ss)/u_1^ssτ^2                 -Q - A/u_1^ssτ^2;                         -P/τ^s_2            -P/τ^s_2-2 ue^ssR/τ^3;                                  ].


Detailed analysis of 𝐡
  shows that these terms vanish rapidly relative to 𝐱
  as the origin is approached, and the relationship
  
    lim_𝐱→ 0sup_τ ⩾ 1𝐡(τ,𝐱)/𝐱 = 0

  holds.
Then the first-order approximation is accurate and defined as
  the linearization of the system
  <cit.>.
The final Cauchy problem can thus be formulated as


    𝐱̇^ss = 𝐉(τ) 𝐱^ss-𝐮̇^ss,
          𝐱^ss(1) = 0,


  where we assume the initial conditions to be zero.

The linear time-varying system (<ref>)
  describes the behaviour of the system
  (<ref>) in a small segment of initial time
  while deviations 𝐱^ss are small.
We will be able to understand how the deviations of 𝐱^ss
  from the steady-state solution grow by analysing the linearized case.

The general solution of (<ref>) in vector form is written as
  Demidovich1967book, p. 77; Kuijstermans2003thesis, p. 41


    𝐱^ss(τ) =
        -∫_1^τ𝐊(τ,τ') 𝐮̇^ss(τ') dτ',


  where 𝐊(τ, τ') = 𝐗(τ)𝐗(τ')^-1
  is a normalized fundamental matrix at point τ',
  and 𝐗(τ) is fundamental matrix.
𝐊(τ, τ') is also called a matriciant or Cauchy matrix.
Our goal is to obtain an expression for 𝐊(τ, τ').
To find the fundamental matrix 𝐗(τ),
  it is convenient to introduce the concept of
  dynamic eigenvalues and eigenvectors
  Wu1980; NeerhoffVanderKloet2001; Kuijstermans2003thesis, p. 47.
If there exists a scalar function λ(τ)
  and a non-zero differentiable vector function 𝐯(τ),
  so that they satisfy the following condition


    [𝐉(τ)-λ(τ)𝐈] 𝐯(τ) = 𝐯̇(τ),


  where 𝐈 is an identity matrix,
  then λ(τ) is called a dynamic eigenvalue
  of matrix 𝐉(τ) associated with
  a dynamic eigenvector 𝐯(τ).
Quasi-static λ^qs(τ) and 𝐯^qs(τ)
  are the classical eigenvalues and eigenvectors
  obtained by solving equation
  (<ref>) with the right-hand side equal to zero.
Similar to the classical algebraic case,
  the number of distinct dynamic eigenvalues λ_i(τ)
  of the matrix is equal to or less than the size of the matrix.
The number of all dynamic eigenvalues of 𝐉
  is equal to the size of the matrix.

Furthermore, we change the variables
  using the transformation
  𝐱^ss = 𝐋(τ) 𝐲,
  where 𝐲 = 𝐲(τ)
  is a vector of new unknown variables.
The 𝐉 matrix is converted to a new matrix,


    𝐁 = 𝐋^-1𝐉𝐋 - 𝐋^-1𝐋̇.


In the classical algebraic case
  for constant 𝐉 and 𝐋,
  𝐁 is Jordan matrix
  whenever 𝐉 has repeated eigenvalues.
In the case with time-varying matrices, one can choose
  an algebraic transformation 𝐋
  so that 𝐁 transforms
  into a diagonal Λ matrix
  <cit.>.
This transformation belongs to the class of Lyapunov transformations
  <cit.>.
Because the coordinate transformation 𝐋 preserves
  the dynamic eigenvalues
  <cit.>,
  Λ(τ) = diag[λ_i(τ)].
The final expression for the fundamental matrix
  𝐗(τ) is


    𝐗(τ) = 𝐋(τ) diag[e^γ_i(τ)],


  where


    γ_i(τ) ≡∫_1^τλ_i(ζ') d ζ'.


Thus, the problem of finding solutions of (<ref>)
  is reduced to the search of its eigenvalues λ_i(τ)
  and the parameters for a Lyapunov transformation
  𝐋 that obtains the diagonal matrix Λ.

One way to find 𝐋(τ) has been proposed by
  <cit.>;
  <cit.>.
This method is based on an iterative algorithm


    Λ̅_j = 𝐐_j^-1(Λ̅_j-1
        -𝐐_j-1^-1𝐐̇_j-1)𝐐_j
          
        (j = 1,2,…),


  with the conditions


    Λ̅_0 = 𝐉,
          𝐐_0 = 𝐈,


  where 𝐈 is an identity matrix,
  Λ̅_j is a diagonal matrix,
  and 𝐐_j is a transformation matrix that consists
  of quasi-static eigenvectors of the matrix
  Λ̅_j-1-𝐐_j-1^-1𝐐̇_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
  𝐐_j-1^-1𝐐̇_j-1.
<cit.> have proved
  that the iterative process converges:


    lim_j →∞Λ̅_j(τ) = Λ(τ),
      
    lim_j →∞𝐐_1(τ)𝐐_2(τ)…𝐐_j(τ) = 𝐋(τ).



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 <ref>),
  the following inequality already holds for the first iterative step:


    ‖Λ̅_1(τ) ‖ =
        ‖𝐐_1^-1(τ)𝐉(τ)𝐐_1(τ) ‖≫‖𝐐_1^-1(τ)𝐐̇_1(τ) ‖.


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


    ‖Λ̅_1(1) ‖ > P
      ≫ 1 >
      ‖𝐐_1^-1(1)𝐐̇_1(1) ‖,

  where P is defined in (<ref>).


The inequality is strengthened as
  ‖Λ̅_1(τ) ‖∼𝒪(τ^2)
  with time τ→∞ and
  ‖𝐐_1^-1(τ)𝐐̇_1(τ) ‖∼𝒪(1/τ).
Therefore, we can conclude
  that in our case the following approximation is perfectly suitable:


    Λ(τ) ≈diag[λ_i^qs(τ)],
          𝐋(τ) ≈𝐕_qs(τ) = 𝐐_1(τ),


  where 𝐕_qs= (𝐯_1^qs, 𝐯_2^qs)
  is the matrix with quasi-static eigenvector columns
  𝐯_1^qs, 𝐯_2^qs
  of the matrix 𝐉
  corresponding to the quasi-static eigenvalues of the Jacobian matrix.
Therefore, expression (<ref>)
  for the fundamental matrix 𝐗
  can be rewritten as


    𝐗(τ) = 𝐕_qs(τ) diag[e^γ^qs_i(τ)]
        =
        (
          𝐯_1^qs(τ)  e^γ_1^qs(τ),
          𝐯_2^qs(τ)  e^γ_2^qs(τ)),


  where γ_i^qs(τ) is calculated from (<ref>)
  for λ_i^qs.
The normalized fundamental matrix looks like


    𝐊(τ, τ')
      =
      𝐕_qs(τ)diag[exp(γ_i^qs(τ)-γ_i^qs(τ'))]𝐕_qs^-1(τ').



The quasi-static eigenvalues and eigenvectors
  for the two-dimensional matrix 𝐉 can be written explicitly as


    λ_1,2^qs(τ) =
        1/2(tr(𝐉) ±√(tr(𝐉)^2 - 4 Δ)),




    𝐯_1,2^qs(τ) = (λ_1,2^qs(τ) - 𝐉[2,2]/𝐉[2,1], 1)^T.


Here, tr(𝐉) is the trace of the Jacobi matrix
  and Δ is its determinant.
It should be noted that
  tr(𝐉) = λ_1^qs+λ_2^qs
  and Δ = λ_1^qsλ_2^qs.

From (<ref>), it follows that
  the function ue^ss/τ^3 decreases as
  ∼τ^-(s_1 + s_2 + 3)/2,
  which is steeper than τ^-s_2.
Because of P ≫ R (see Appendix <ref>),
  P/τ^s_2≫ 2 R ue^ss/τ^3
  for any time τ⩾ 1.
Therefore, 𝐉[2,2] ≃𝐉[2,1]
  and


    𝐯_1,2^qs(τ) ≃(-1 - λ_1,2^qs(τ)/Pτ^s_2, 1)^T


  according to (<ref>).

Consider now the first eigenvalue of the matrix 𝐉
  defined in (<ref>).
  It is equal to


    λ_1^qs(τ)
        =
        -ue^ssR/τ^3(2 Q u_1^ss + (1/u_1^ss+1) Aτ^2)/Aτ^2+u_1^ss P/τ^s_2,

  which takes into account (<ref>).

The expression (<ref>) implies that λ_1^qs(τ)
  is a negative time function.
One can ensure that initially
  |λ_1^qs(1)| ≪ P
  (see Appendix <ref>).
As time increases, the module of the function |λ_1^qs(τ)|
  decreases monotonically faster than
  ue^ss R/τ^3 ∼τ^-(s_1 + s_2 + 3)/2,
  which is steeper than τ^-s_2.
This means that |λ_1^qs(τ)| ≪ P/τ^s_2
  for any τ⩾ 1.
Thus, the first eigenvector is simplified to


    𝐯_1^qs(τ) ≃ (-1, 1)^T.



Consider now the behaviour of the function tr(𝐉(τ)).
This is a negative function of τ.
For its module, the inequality |λ_1^qs(1)| ≪ P < |tr(𝐉(1))| holds.
Because P ≫ A ≫ R
  (see Appendix <ref>),
  the module |tr(𝐉(τ))| initially
  decreases as P/τ^s_2 with time
  not steeper than |λ_1^qs(τ)|
  and then increases as A τ^2.
Hence, |λ_1^qs(τ)| ≪ |tr(𝐉(τ))|
  for any τ⩾ 1.
And finally, from


    λ_2^qs(τ) = tr(𝐉(τ)) - λ_1^qs(τ)


  it follows that |λ_1^qs(τ)| ≪ |λ_2^qs(τ)|
  and λ_2^qs(τ) ≃tr(𝐉(τ))
  for any τ⩾ 1.
Thus, the second eigenvector can be rewritten as


    𝐯_2^qs(τ) ≃(-1-tr(𝐉(τ))/Pτ^s_2, 1)^T.



For small time,
  as long as |tr(𝐉(τ))|
  is dominated by the photoionization term P/τ^s_2, we have


    ∫_τ'^τλ_2^qs(ζ) d ζ≃
        -∫_τ'^τP/ζ^s_2 d ζ≃
        -P/τ^s_2 (τ-τ') + 𝒪(τ-τ')^2,


but at large time, when the spontaneous emission A τ^2 dominates, we can write


    ∫_τ'^τλ_2^qs(ζ) d ζ≃
        -∫_τ'^τ A ζ^2 d τ≃
        -A τ^2 (τ-τ') + 𝒪(τ-τ')^2.


Because


    min(|λ_2^qs(τ)|) > A(P/A)^2/2+s_2≫ 1,

  one has that
  exp(γ_2^qs(τ)-γ_2^qs(τ'))
  is an exponentially decreasing function of
  (τ-τ'), for any τ⩾ 1.
One can estimate the normalized fundamental matrix
  by summarizing all the above observations.
From (<ref>),
  it follows that when τ' and τ are nearly close,


    (τ-τ') ≲min(|λ_1^qs(τ)|)^-1≪ 1,


  𝐊(τ, τ') ∼𝐈.
In the opposite case, we assume that
  exp(γ_2^qs(τ)-γ_2^qs(τ')) ≃ 0
and find


    𝐊(τ, τ') ≃exp(γ_1^qs(τ)-γ_1^qs(τ'))
          /1+𝐯_2^qs[1](τ')[             -1 -𝐯_2^qs[1](τ');             -1 -𝐯_2^qs[1](τ');                ].



Finally, we consider the integrand function in (<ref>).
If (<ref>) is true
  then 𝐊(τ,τ') 𝐮̇^ss(τ') ∼𝐮̇^ss(τ),
  and the contribution of this term can be neglected.
When (<ref>) is not satisfied, the estimate can be written as


    𝐊(τ,τ') 𝐮̇^ss(τ') ≃exp(γ_1^qs(τ)-γ_1^qs(τ')) 
        M(τ')
        [ -1; -1;    ],

  where


    M(τ) =
          u̇_1^ss(τ) - u̇e^ss(τ) 𝐯_2^qs[1](τ)/1+𝐯_2^qs[1](τ) =
          
    =u̇_1^ss(τ) -(
            1+P/tr(𝐉(τ)) τ^s_2)
          (u̇_1^ss(τ) +u̇e^ss(τ)).


The analysis of the behaviour of the function tr(𝐉(τ))
  presented above implies that
  
    (1+P/tr(𝐉(τ)) τ^s_2) ≲ 2.

And because |(u̇_1^ss(τ) + u̇e^ss(τ))| ≪ |u̇_1^ss(τ)|
  for any τ⩾ 1,
  we obtain the final answer


    𝐱^ss(τ) ≃∫_1^τexp(γ_1^qs(τ)-γ_1^qs(τ'))
         u̇_1^ss(τ') dτ'
        [ -1; -1;    ].


The solution (<ref>) includes only the larger
  of the two eigenvalues λ_1^qs(τ),
  which is the smaller in absolute value
  because they are negative.
The smaller the |λ_1^qs(τ)| 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 (<ref>)
  of the system (<ref>)
  is applicable only when the deviations
  𝐱^ss are small.
However, we have previously seen that the deviations can be significant.
Fig. <ref>(a) presents the solution x_1^ss(t)
  of the system (<ref>)
  and its linearization x_1^ss,lin(t)
  derived from (<ref>)
  for physical time for the case of outer, high-velocity layers
  far from the photosphere
  (Table <ref> column 1).
  while Fig. <ref>(b) presents
  solutions of the same systems
  for the case of near-photospheric layers
  (Table <ref> 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. <ref>,<ref>).
Thus, the solution (<ref>) describes
  the evolution of the system only at early time.





However, in the next section,
  we show that the final expression for |λ_1^qs(τ)|
  determines whether the kinetic system is in equilibrium
  at any moment of time.






§ 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
  Khalil2002book, p. 365; Vidyasagar2002book, p. 248.

Let us rewrite the system of equations
  (<ref>, <ref>)
  where all the variable coefficients are fixed at some moment τ_1.
In fact, a non-autonomous system have transformed into an autonomous one.


    u̇_1^a = (1-u_1^a-ue^a) (Q + A/u_1^aτ_1^2)
        -B/τ_1^s_1-2,



    u̇e^a = (1-u_1^a-ue^a)P/τ_1^s_2
        - (ue^a)^2R/τ_1^3.


Unknown u_1^a(τ_a), ue^a(τ_a)
  are the functions of a new time variable τ_a.
Obviously, the values of expressions u_1^ss(τ_1), ue^ss(τ_1),
  are the solutions of the system
  (<ref>, <ref>).
The right-hand side of
  (<ref>, <ref>)
  vanishes for those values.

Let us prove that for any solution of the system
  (<ref>, <ref>)
  for τ_a →∞
  under initial conditions 0 < u_1,e^a(τ_a=1) ⩽ 1
  one obtains u_1,e^a(τ_a) → u_1,e^ss(τ_1).

The change of variables
  x_1(τ_a) = u_1^a(τ_a) - u_1^ss(τ_1)
  and xe(τ_a) = ue^a(τ_a) - ue^ss(τ_1)
  transforms the original system to the reduced one
  Demidovich1967book, p. 234; Khalil2002book, p. 147:


    ẋ_1 = -(x_1+xe) Q
        -
        (1-ue^ss(τ_1)) x_1+u_1^ss(τ_1) xe/u_1^ss(τ_1)(u_1^ss(τ_1)+x_1)
        Aτ_1^2,



    ẋe = -(x_1+xe)P/τ_1^s_2
        -xe (2 ue^ss(τ_1)+xe)R/τ_1^3.


The Jacobi matrix of this system
  coincides with the matrix (<ref>)
  taken at time τ_1: 𝐉 = 𝐉(τ_1).
From the analysis carried out in Section <ref>,
  it follows that the two eigenvalues λ_1, λ_2
  of the Jacobi matrix 𝐉 are negative and different.
Consequently, by Lyapunov's indirect method
  Fedoruk1985book, p. 289; Khalil2002book, p. 133, it appears that
  the system (<ref>, <ref>)
  is exponentially stable,
  and its solution is written as
  Fedoruk1985book, p. 61; Khalil2002book, p. 37


    𝐱(τ_a) = (x_1, xe)^T = C_1 e^λ_1 τ_a𝐯_1 + C_2 e^λ_2 τ_a𝐯_2,


  where 𝐯_i are the eigenvectors of the Jacobi matrix 𝐉
  and C_i are arbitrary constants.
Obviously, lim_τ→∞𝐱(τ_a) = 0,
  which shows the convergence property of the system
  (<ref>, <ref>)
  <cit.>.
Thus, the search for the limit of solutions of
  (<ref>, <ref>)
  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
  τ_ss is needed for the solution u_1^a(τ_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 (<ref>), it follows that
  the period of equilibration
  is determined by the smaller of the two modules of the eigenvalues
  τ_ss^-1 = min(|λ_1(τ_1)|, |λ_2(τ_1)|).
From Section <ref>, we know that
  |λ_1(τ_1)| ≪ |λ_2(τ_1)|,
  and thus τ_ss^-1 = |λ_1(τ_1)| holds.
The relaxation time τ_ss for any τ_1 from (<ref>) is


    τ_ss(τ)
        =
        Aτ^2+u_1^ssP/τ^s_2/ue^ssR/τ^3(2 Q u_1^ss + (1/u_1^ss+1) Aτ^2).


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 (<ref>)
  with the typical time-scale of change of SN envelope's physical parameters.
The calculations show (Fig. <ref>)
  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 = τ_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 τ_ss
  increases with time (Fig. <ref>),
  so at τ→∞ the time-dependent effect
  will always be important,
  which has already been proved in Section <ref>.
However, we have seen (Fig. <ref>) 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 (<ref>)
  for the dense near-photospheric layers
  can be estimated as


    τ_ss∼1/ue^ss RP/Q
        =
        t_recP/Q,

  where t_rec is the classical recombination time
  <cit.>.
This formula is consistent with the estimates of the effective recombination time-scale
  UtrobinChugai2005, eq. 4; DessartHillier2007, eq. 21.

However, a more general conclusion can be obtained from (<ref>)
  for layers of any density located both
  near and far from the photosphere.
As the constant B is small
  (see Table <ref>, Appendix <ref>),
  τ_ss in equation (<ref>) can be expand as


    τ_ss(τ)
        ≈
        f(τ,A,P,s_2,R)√(τ^s_1/B)
        +
        𝒪(√(B)).

Here, f(τ,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/τ^s_2
  dominates the spontaneous emission Aτ^2,
  and an expression for f can be obtained explicitly.
Then (<ref>) will be rewritten as


    τ_ss(τ)
        ≈√(A τ^3/R)/2(Q + A τ^2)√(P/τ^s_2)√(τ^s_1/B) =
        
    =√(A τ^3/(R B̃))/2(Q + A τ^2)√(P_2c(t_0)/J_c(ν_Lα, t_0))(t/t_0)^s_1-s_2/2,

  where we used equations
  (<ref>), (<ref>), (<ref>), (<ref>).
Let us rewrite this expression for the relaxation time
  using physical dimensional values:


    t_ss(t)
        ≈f̃(t,N_0,Te)
        √(P_2c(t_0)/J_c(ν_Lα, t_0))(t/t_0)^s_1-s_2/2.

The minimum and maximum of the function f̃
  for the domain of N_0, Te (see Appendix <ref>)
  at any given time t ⩾ 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α 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.










§ 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 <cit.>.
Calculations are carried out using the original code LEVELS
  and based on the model  for SN 1999em
  <cit.>.
The physical parameters for the calculation are taken from
  the Lagrangian layer with the velocity
  v ≈ 5.8 · 10^3  km s^-1.
This layer is close to the photosphere
  at the beginning of the plateau epoch
  and its characteristics were used in the
  Table <ref>, column 2.
Eigendecomposition of the system Jacobi matrix is done
  by the FEAST eigensolver package
  <cit.>.





We considered four systems presented in Table <ref>.
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
  <cit.>.
For model M_4,
the selective absorption of Lα photons
  in metal lines (for example, Fe II)
  <cit.>
  and overlapping lines
  <cit.>
  is taken into account.
The metal abundances are taken from the
   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
  <cit.>.

Fig. <ref> shows the results of the calculation.
The thin solid line “ph” corresponds
  to the analytical estimation
  obtained with the formula (<ref>).
This is the same line as that one designated by “ph” in Fig. <ref>.
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 <cit.>,
  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 <ref>,
  we can conclude the following.
The concentration of Lα photons that come from outside,
  as well as the change of this concentration
  and the rate of the selective absorption of Lα
  photons in metal lines
  are the most significant factors affecting the non-stationarity effect.
The greater the number of external Lα 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
  <cit.>.
On one hand,
  this means that the low metallicity
  leads to an increase of the number
  of external Lα 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. <ref>, model M_4).






§ 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. So­ro­ki­na for
  thorough reading the manuscript before publishing.
M.Sh. Po­ta­shov 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 <ref>
  were funded by the Russian Science Foundation, grant 19–12–00229.
The work of AY was supported by the RFBR, grant 18–29–21019.




@urlcharsothermakeother $&#_%
@doi@urlcharsother ifnextchar [ @doi@
  @doi@[]
@doi@[#1]#2tempa#1tempaempty http://dx.doi.org/#2 doi:#2http://dx.doi.org/#2 #1
@eprint#1#2@eprint@#1:#2::nil
@eprint@arXiv#1http://arxiv.org/abs/#1 arXiv:#1
@eprint@dblp#1http://dblp.uni-trier.de/rec/bibtex/#1.xml
  dblp:#1
@eprint@#1:#2:#3:#4niltempa #1tempb #2tempc
  #3tempc empty tempc tempb tempb tempa tempb empty tempb arXivifundefined
  mn@eprint@tempbtempb:tempcmn@eprint@tempbtempc

[AndronovaAndronova1991]Andronova1990
Andronova A. A.,  1991, @doi [] 10.1007/BF01005505, 32, 235

[AppellAppell1889]Appell1889
Appell P.,  1889, Journal de Mathématiques Pures et Appliquées, 5,
  361

[Baker & MenzelBaker &
  Menzel1938]BakerMenzel1938
Baker J. G.,  Menzel D. H.,  1938, @doi [] 10.1086/143959, 88, 52

[Baklanov, Blinnikov  & PavlyukBaklanov
  et al.2005]BaklanovBlinnikovPavlyuk2005
Baklanov P. V.,  Blinnikov S. I.,   Pavlyuk N. N.,  2005, @doi []
  10.1134/1.1958107, 31, 429

[Baklanov, Blinnikov, Potashov  &
  DolgovBaklanov et al.2013]BaklanovBlinnikovPotashovEtal2013a
Baklanov P. V.,  Blinnikov S. I.,  Potashov M. S.,   Dolgov A. D.,  2013,
  @doi [JETP Letters] 10.1134/S0021364013200034, 98, 432

[Baron, Nugent, Branch  & HauschildtBaron
  et al.2004]BaronNugentBranchEtal2004
Baron E.,  Nugent P. E.,  Branch D.,   Hauschildt P. H.,  2004, @doi []
  10.1086/426506, 616, L91

[Bartunov & MozgovoiBartunov &
  Mozgovoi1987]BartunovMozgovoi1987
Bartunov O. S.,  Mozgovoi A. L.,  1987, @doi [] 10.1007/BF01005551,
  26, 136

[Berezin & ZhidkovBerezin &
  Zhidkov1959]BerezinZhidkov1959book
Berezin I.,  Zhidkov N.,  1959, Computing Methods 2 (in Russian).
GIFML, Moscow

[Blinnikov, Eastman, Bartunov, Popolitov  &
  WoosleyBlinnikov et al.1998]BlinnikovEastmanBartunovEtal1998
Blinnikov S. I.,  Eastman R. G.,  Bartunov O. S.,  Popolitov V. A.,   Woosley
  S. E.,  1998, @doi [] 10.1086/305375, 496, 454

[Blinnikov, Lundqvist, Bartunov, Nomoto  &
  IwamotoBlinnikov et al.2000]BlinnikovLundqvistBartunovEtal2000
Blinnikov S. I.,  Lundqvist P.,  Bartunov O. S.,  Nomoto K.,   Iwamoto K.,
  2000, @doi [] 10.1086/308588, 532, 1132

[Blinnikov, Röpke, Sorokina,
  Gieseler, Reinecke, Travaglio, Hillebrandt  &
  StritzingerBlinnikov et al.2006]BlinnikovRopkeSorokinaEtal2006
Blinnikov S. I.,  Röpke F. K.,  Sorokina E. I.,  Gieseler M.,
  Reinecke M.,  Travaglio C.,  Hillebrandt W.,   Stritzinger M. D.,
  2006, @doi [] 10.1051/0004-6361:20054594, 453, 229

[Blinnikov, Potashov, Baklanov  &
  DolgovBlinnikov et al.2012]BlinnikovPotashovBaklanovEtal2012
Blinnikov S. I.,  Potashov M. S.,  Baklanov P. V.,   Dolgov A. D.,  2012,
  @doi [JETP Letters] 10.1134/S0021364012150039, 96, 153

[CastorCastor1970]Castor1970
Castor J. I.,  1970, @doi [] 10.1093/mnras/149.2.111, 149, 111

[ChugaiChugai1980]Chugai1980a
Chugai N. N.,  1980, Soviet Astronomy Letters, 6

[ChugaiChugai1987]Chugai1987
Chugai N. N.,  1987, , 26, 89

[ChugaiChugai1988a]Chugai1988b
Chugai N. N.,  1988a, , 29, 74

[ChugaiChugai1988b]Chugai1988a
Chugai N. N.,  1988b, Nauchnye Informatsii, 65

[ChugaiChugai1991]Chugai1991
Chugai N. N.,  1991, in , Supernovae.
Springer New York, New York, NY, pp 286–290,
  @doi10.1007/978-1-4612-2988-9_40

[ChugaiChugai1998]Chugai1998
Chugai N. N.,  1998, , 24

[De, Baron  & HauschildtDe
  et al.2010]DeBaronHauschildt2010a
De S.,  Baron E.,   Hauschildt P. H.,  2010, @doi []
  10.1111/j.1365-2966.2009.15808.x, 401, 2081

[DemidovichDemidovich1967]Demidovich1967book
Demidovich B.,  1967, Lectures on the Mathematical Stability Theory (in
  Russian).
Nauka, Moscow

[Dessart & HillierDessart &
  Hillier2007]DessartHillier2007
Dessart L.,  Hillier D. J.,  2007, @doi []
  10.1111/j.1365-2966.2007.12538.x, 383, 57

[Dessart & HillierDessart &
  Hillier2010]DessartHillier2010a
Dessart L.,  Hillier D. J.,  2010, @doi []
  10.1111/j.1365-2966.2010.16611.x, 405, 23

[Dessart et al.,Dessart
  et al.2008]DessartBlondinBrownEtal2008
Dessart L.,  et al., 2008, @doi [] 10.1086/526451, 675, 644

[Dessart, Hillier, Waldman  & LivneDessart
  et al.2013]DessartHillierWaldmanEtal2013
Dessart L.,  Hillier D. J.,  Waldman R.,   Livne E.,  2013, @doi []
  10.1093/mnras/stt861, 433, 1745

[Ezquiaga & ZumalacárreguiEzquiaga &
  Zumalacárregui2018]EzquiagaZumalacarregui2018
Ezquiaga J. M.,  Zumalacárregui M.,  2018, @doi [Frontiers in
  Astronomy and Space Sciences] 10.3389/fspas.2018.00044, 5

[FedorukFedoruk1985]Fedoruk1985book
Fedoruk M.,  1985, Ordinary differential equations (in Russian).
Nauka, Moscow

[GrachevGrachev1988]Grachev1988
Grachev S. I.,  1988, @doi [] 10.1007/BF01014860, 28, 119

[GrachevGrachev1989]Grachev1989
Grachev S. I.,  1989, @doi [] 10.1007/BF01003949, 30, 211

[Hillier & DessartHillier &
  Dessart2012]HillierDessart2012
Hillier D. J.,  Dessart L.,  2012, @doi []
  10.1111/j.1365-2966.2012.21192.x, 424, 252

[Hubeny & LanzHubeny &
  Lanz1995]HubenyLanz1995
Hubeny I.,  Lanz T.,  1995, @doi [] 10.1086/175226, 439, 875

[Hubeny & MihalasHubeny &
  Mihalas2014]HubenyMihalas2014book
Hubeny I.,  Mihalas D.,  2014, Theory of Stellar Atmospheres.
Princeton University Press

[HummerHummer1994]Hummer1994
Hummer D. G.,  1994, @doi [] 10.1093/mnras/268.1.109, 268, 109

[Hummer & RybickiHummer &
  Rybicki1985]HummerRybicki1985
Hummer D. G.,  Rybicki G. B.,  1985, @doi [] 10.1086/163232, 293, 258

[Hummer & RybickiHummer &
  Rybicki1992]HummerRybicki1992
Hummer D. G.,  Rybicki G. B.,  1992, @doi [] 10.1086/171076, 387, 248

[Hummer & StoreyHummer &
  Storey1987]HummerStorey1987
Hummer D. G.,  Storey P. J.,  1987, @doi [] 10.1093/mnras/224.3.801,
  224, 801

[Hummer & StoreyHummer &
  Storey1992]HummerStorey1992
Hummer D. G.,  Storey P. J.,  1992, @doi [] 10.1093/mnras/254.2.277,
  254, 277

[Kestyn, Polizzi  & Peter TangKestyn
  et al.2016]KestynPolizziTang2016
Kestyn J.,  Polizzi E.,   Peter Tang P. T.,  2016, @doi [SIAM Journal on
  Scientific Computing] 10.1137/15M1026572, 38, S772

[KhalilKhalil2002]Khalil2002book
Khalil H. K.,  2002, Nonlinear Systems.
Pearson Education, Prentice Hall

[Kirshner & KwanKirshner &
  Kwan1974]KirshnerKwan1974
Kirshner R. P.,  Kwan J.,  1974, @doi [] 10.1086/153123, 193, 27

[Kirshner & KwanKirshner &
  Kwan1975]KirshnerKwan1975
Kirshner R. P.,  Kwan J.,  1975, @doi [] 10.1086/153527, 197, 415

[KuijstermansKuijstermans2003]Kuijstermans2003thesis
Kuijstermans F. C.,  2003, PhD thesis, Enschede, @doi10.6100/ir807325

[Kuntsevich & LychakKuntsevich &
  Lychak1977]KuncevichLychak1977book
Kuntsevich V.,  Lychak M.,  1977, Synthesis of Automatic Control Systems with
  the Help of Lyapunov Functions (in Russian).
Nauka, Moscow

[Kurt & ShakhvorostovaKurt &
  Shakhvorostova2014]KurtShakhvorostova2014a
Kurt V. G.,  Shakhvorostova N. N.,  2014, @doi [Physics-Uspekhi]
  10.3367/UFNe.0184.201404f.0423, 57, 389

[MihalasMihalas1978]Mihalas1978book
Mihalas D.,  1978, Stellar atmospheres /2nd edition/.
W. H. Freeman and Co.

[Mörtsell & DhawanMörtsell &
  Dhawan2018]MortsellDhawan2018
Mörtsell E.,  Dhawan S.,  2018, @doi [Journal of Cosmology and
  Astroparticle Physics] 10.1088/1475-7516/2018/09/025, 2018

[Neerhoff & van der KloetNeerhoff & van der
  Kloet2001]NeerhoffVanderKloet2001
Neerhoff F.,  van der Kloet P.,  2001, in ISCAS 2001. The 2001 IEEE
  International Symposium on Circuits and Systems (Cat. No.01CH37196). IEEE, pp
  779–782, @doi10.1109/ISCAS.2001.921448

[Nussbaumer & SchmutzNussbaumer &
  Schmutz1984]NussbaumerSchmutz1984
Nussbaumer H.,  Schmutz W.,  1984, , 138

[OlsonOlson1982]Olson1982
Olson G. L.,  1982, @doi [] 10.1086/159825, 255, 267

[Osterbrock & FerlandOsterbrock &
  Ferland2006]OsterbrockFerland2006book
Osterbrock D. E.,  Ferland G. J.,  2006, Astrophysics of gaseous nebulae and
  active galactic nuclei.
University Science Books, Sausalito

[PeeblesPeebles1968]Peebles1968
Peebles P. J. E.,  1968, @doi [] 10.1086/149628, 153, 1

[Peter Tang & PolizziPeter Tang &
  Polizzi2014]TangPolizzi2014
Peter Tang P. T.,  Polizzi E.,  2014, @doi [SIAM Journal on Matrix
  Analysis and Applications] 10.1137/13090866X, 35, 354

[PolizziPolizzi2009]Polizzi2009
Polizzi E.,  2009, @doi [] 10.1103/PhysRevB.79.115112, 79, 115112

[Polyanin & ZaitsevPolyanin &
  Zaitsev2002]PolyaninZaitsev2002book
Polyanin A. D.,  Zaitsev V. F.,  2002, Handbook of exact solutions for ODE.
Chapman & Hall/CRC

[Potashov, Blinnikov, Baklanov  &
  DolgovPotashov et al.2013]PotashovBlinnikovBaklanovEtal2013
Potashov M. S.,  Blinnikov S. I.,  Baklanov P. V.,   Dolgov A. D.,  2013,
  @doi [] 10.1093/mnrasl/slt016, 431, L98

[Potashov, Blinnikov  & UtrobinPotashov
  et al.2017]PotashovBlinnikovUtrobin2017a
Potashov M. S.,  Blinnikov S. I.,   Utrobin V. P.,  2017, @doi []
  10.1134/S1063773717010030, 43, 36

[RaizerRaizer1959]Raizer1959
Raizer Y. P.,  1959, JETP, Vol. 10, 411

[Riess et al.,Riess
  et al.2018]RiessCasertanoYuanEtal2018
Riess A. G.,  et al., 2018, @doi [] 10.3847/1538-4357/aac82e, 861, 126

[SemenovSemenov2014]Semenov2014
Semenov E.,  2014, The Bulletin of Irkutsk State University. Series
  Mathematics, 7, 124

[SobolevSobolev1960]Sobolev1960book
Sobolev V. V.,  1960, Moving envelopes of stars.
Harvard University Press

[UtrobinUtrobin2007]Utrobin2007
Utrobin V. P.,  2007, @doi [] 10.1051/0004-6361:20066078, 461, 233

[Utrobin & ChugaiUtrobin &
  Chugai2002]UtrobinChugai2002
Utrobin V. P.,  Chugai N. N.,  2002, @doi [] 10.1134/1.1484138, 28,
  386

[Utrobin & ChugaiUtrobin &
  Chugai2005]UtrobinChugai2005
Utrobin V. P.,  Chugai N. N.,  2005, @doi []
  10.1051/0004-6361:20042599, 441, 271

[Van Der Kloet & NeerhoffVan Der Kloet
  & Neerhoff2000]VanderKloetNeerhoff2000
Van Der Kloet P.,  Neerhoff F. L.,  2000, @doi [International Journal of
  Systems Science] 10.1080/002077200412195, 31, 1053

[VidyasagarVidyasagar2002]Vidyasagar2002book
Vidyasagar M.,  2002, Nonlinear Systems Analysis.
Classics in Applied Mathematics, Society for Industrial and Applied Mathematics

[Vogl, Sim, Noebauer, Kerzendorf  &
  HillebrandtVogl et al.2019]VoglSimNoebauerEtal2019
Vogl C.,  Sim S. A.,  Noebauer U. M.,  Kerzendorf W. E.,   Hillebrandt W.,
  2019, @doi [] 10.1051/0004-6361/201833701, 621, A29

[WuWu1980]Wu1980
Wu M.-Y.,  1980, @doi [IEEE Transactions on Automatic Control]
  10.1109/TAC.1980.1102448, 25, 824

[Zeldovich & RaizerZeldovich &
  Raizer2002]ZeldovichRajzer2008book
Zeldovich Y. B.,  Raizer Y. P.,  2002, Physics of Shock Waves and
  High-Temperature Hydrodynamic Phenomena.
Dover Publications

[Zeldovich, Kurt  & SyunyaevZeldovich
  et al.1969]ZeldovichKurtSunyaev1969
Zeldovich Y. B.,  Kurt V. G.,   Syunyaev R. A.,  1969, Soviet Journal of
  Experimental and Theoretical Physics, 28, 146












§ 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
  ν⩾ν_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
  <cit.>:


    χ_ν = (n_1 - n_1^* e^-h νk Te) α_1c(ν),



    η_ν = n_1^* (1 - e^-h νk Te) α_1c(ν)B_ν(Te).


Here α_1c(ν) is the photoionization cross section from the first level,
  B_ν(Te) is the intensity
  of blackbody radiation for the electron temperature Te, and


    n_1^* = nenpΦ_Saha(Te)


  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 ne and np.

For an optically thick medium,
  the transport equation in the continuum is solved as
  J_c(ν, t) = S_c(ν, t),
  where S_c(ν, t) is a source function,
  defined as


    S_c(ν, t) = η_ν/χ_ν.


Substituting (<ref>) and (<ref>) into (<ref>)
  we have


    J_c(ν, t) = (1 - e^-h νk Te)/n_1/n_1^* - e^-h νk TeB_ν(Te).



In accordance with
  <cit.>
  and
  <cit.>,
  the photoionization rate from the ground level is the integral


    n_1 P_1c(t) = 4π n_1-10mu∫_ν_LyC^∞-10muα_1c(ν)/hν J_c(ν, t) dν.


The photorecombination rate looks like


    nenp R_c1(t)
        = 4πnenp Φ_Saha(Te)
        ×-10mu∫_ν_LyC^∞-10muα_1c(ν)/hν(2hν^3/c^2
        +
        J_c(ν, t))e^-hν/kTedν.


Substituting equation (<ref>) into (<ref>) and (<ref>),
  and taking into account (<ref>) gives n_1 P_1c(t) = nenp 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.






§ 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
  <cit.>,
  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 “ph”
  and to the outer, high-velocity layers with “out”.

The initial time is t_0 ≈ 20 d after the explosion
  in the core.

The two-photon 2s→ 1s decay rate is taken from
  <cit.>.
We consider two l-sublevels 2s and 2 p
  as a single super-level 2
  <cit.>
  and assume that the populations of the sublevels are proportional to their
  statistical weights g_2s, g_2p.
This leads to


    Q ≈ 8.2249 g_2s/g_2s+g_2p t_0 ≈
      2.74 t_0
      ≈ 3.5 · 10^6.



The initial number densities of  gas are
  N_0(ph) ≈ 10^11 cm^-3
  and
  N_0(out) ≈ 10^8  cm^-3.
It follows that for the constant A
  defined in (<ref>, <ref>) we have corresponding values


    A(ph) ≈ 3·10^4,  
      A(out) ≈ 3·10^7.


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


    J_c(ν_Lα, t_0, ph) ≈ 10^-8 erg/cm^2 s Hz.
      
    
      J_c(ν_Lα, t_0, out) ≈ 10^-12 erg/cm^2 s Hz,


Thus, for B defined in (<ref>, <ref>) we have


    B(ph) ≈ 6·10^-3,  
      B(out) ≈ 6·10^-4.



The photoionization coefficient from the second level P_2c
  at the initial moment is also
  taken from STELLA calculations:


    P_2c(t_0,ph) ≈ 10^5 s^-1,  
      P_2c(t_0,out) ≈ 4·10^3 s^-1.


Then we get from (<ref>):


    P(ph) ≈ 2·10^11,  
      P(out) ≈ 7·10^9.



The effective photorecombination coefficient to the second level
  at the initial moment can be estimated from
  <cit.>
  using the material temperatures typical for a supernova envelope
  Te≈ 3000 - 7000 K as
  R_2c≈ 3 · 10^-13 cm^3/s.
Then we will have for R from (<ref>, <ref>):


    R(ph) ≈ 5·10^4,  
      R(out) ≈ 50.


Values for the power exponents taken from STELLA calculations are


    s_1(ph) ≈ 24,  
      s_1(out) ≈ 21,
      
    
      s_2(ph) ≈ 6,  
      s_2(out) ≈ 4.


Because in the optically thin atmosphere
  s_1 = 2 and s_2 = 2,
  in general case, the domains of these exponents are
  s_1 ⩾ 2 and s_2 ⩾ 2.






§ THE TUBE
In general case, we limit the domain of these exponents to
  s_1 ⩾ 2 and s_2 ⩾ 2.

  

Let us prove the dissipativity of the system
  (<ref>, <ref>)
  in another way.
Suppose that we know some particular bound solution
  0 <ũ_1 ⩽ 1 and 0 ⩽ũe⩽ 1
  (unperturbed motion)
  of a nonautonomous nonlinear differential system
  (<ref>, <ref>).
The change of variables
  x_1 = u_1 - ũ_1 and xe = ue - ũe
  transforms the system into the form of perturbed motion system
  Demidovich1967book, p. 234; Khalil2002book, p. 147:


    ẋ_1 = -(x_1+xe) Q
        - (1-ũe)x_1+ũ_1 xe/ũ_1(ũ_1+x_1)Aτ^2,




    ẋe = - (x_1+xe)P/τ^s_2
        -xe (2ũe+xe)R/τ^3.


The trivial solution x_1 = 0, xe = 0
  is an equilibrium point of the transformed system for any τ⩾ 1.

Consider the Lyapunov function candidate


    V(τ,x_1,xe) =
        x_1^2 (1-ũe)/ũ_1
        + 2 x_1 xe
        + xe^2 (1+d).


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
  (<ref>, <ref>) is given by


    V̇(τ, x_1, xe) = ∂ V/∂τ
        + ∂ V/∂ x_1ẋ_1 + ∂ V/∂xeẋe,


  where ẋ_1 and ẋe
    are (<ref>) and (<ref>), respectively.
At large time and at small x_1 and xe it expresses as


    V̇(τ, x_1, xe) =
        -2 A τ^2/ũ_1[(1-ũe)(1-ũe-ũ_1) x_1^2/2 ũ_1^2
        +
        ((1-ũe)/ũ_1x_1+xe)^2] + 𝒪(τ).


This function is negative semidefinite in the domain
  Ω = (0 < ũ_1 + x_1 ⩽ 1, 0 < ũe
  +xe⩽ 1).

From the Lyapunov direct method and (<ref>, <ref>)
  it follows that the trivial solution of system (<ref>, <ref>)
  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
  τ > τ_1
  with τ_1 ⩾ 1.
The stability on a half-axis τ > 1 is obtained taking into account
  the theorem of continuous dependence on the parameter
  <cit.>
  for the solution on a finite interval 1 ⩽τ⩽τ_1.

The analysis of the function (<ref>) reveals that the system
  (<ref>, <ref>)
  is  dissipative.
Indeed, Section <ref> shows
  that ũ_1 ≈ 1-ũe
  at large time.
Then (<ref>) can be rewritten as


    V̇(τ, x_1, xe) = -2 A τ^2/ũ_1(x_1+xe)^2 + 𝒪(τ).


This function is a parabolic cylinder
  with the zero values along the line l
  defined by the equation x_1 = -xe.
This means that one can specify a value μ
  so that the set l does not belong
  to the Ω̃ domain
  defined as the difference between the Ω domain
  and the domain specified by the inequality


    sup_τ⩾ 1 V(τ,x_1,xe) ⩽μ.


  for all moments of time τ⩾ 1.
Thus, the function V̇(τ, x_1, xe)
  is negative definite in Ω̃.
Therefore, according to the theorem of Yoshizawa
  Demidovich1967book, p. 290; KuncevichLychak1977book, p. 47-48, 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 xe.
Let us prove by contradiction that dissipativity is true for any
  x_1, xe in Ω.
Note additionally that the solutions of the initial value problem
  (<ref>, <ref>, <ref>)
  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
  (<ref>, <ref>) with the initial conditions
  x_1(1) ≈ũ_1(1),
  xe(1) = 1 - ũe(1).
It was shown in Section <ref>
  that the initial condition of u_2 does not affect the solution.
Because u_2 = 1- u_1 - ue,
  then from (<ref>) it follows that the
  expression -(x_1+xe)
  decreases exponentially as e^-G_3(τ).
Therefore, without a loss of generality,
  only the behaviour of x_1 is considered below.
Fig. <ref> 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 ∼ 10 d.
The situation changes dramatically
  if one takes into account collisional processes.
For example, in Fig. <ref> 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^col≪ 1 at time greater than 1 d)
  one can say that the system “forgets” the initial conditions!







