Partial Ly\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}α thermalization in an analytic nonlinear diffusion model

During recombination, the cosmic background radiation is disturbed, in particular, by Lyman-alpha emissions from neutral hydrogen. It is proposed to account for the subsequent time-dependent partial thermalization of the Ly\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}α energy content in an analytically solvable nonlinear diffusion model. The amplitude of the partially thermalized and redshifted Ly\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}α line is found to be too low to be visible in the cosmic microwave spectrum, in accordance with previous numerical models and Planck observations.


Introduction
Following the predictions [1,2], and discovery of the cosmic microwave background (CMB) with a temperature of (3.5 ± 1.0) K at a frequency of 4080 MHz [3], its spectrum has been measured with ever increasing precision.The radiation has a Planck distribution, because the cosmic background radiation (CBR) had been thermalized essentially through Compton scattering and bremsstrahlung [4] at very early times corresponding to redshifts z > 10 7 , and expansion retains the thermal spectrum.Ground-, balloon-, and rocketbased observations confirmed the low-frequency Rayleigh-Jeans branch of the The nonlinear boson diffusion equation (NBDE) [13,14] accounts for timedependent (partial or complete) thermalization towards the Bose-Einstein stationary distribution that is reached in the limit t → ∞.In the general case of frequency-dependent transport coefficients, it can only be solved numerically, but analytical solutions exist for constant drift and diffusion coefficients.It is one of the few nonlinear partial differential equations with a clear physical meaning that has analytical solutions.In this work, these solutions are applied to the partial thermalization of the Lyα line that is emitted at recombination, corresponding to a redshift of the last-scattering surface z rec ≃ 1100, and an average recombination temperature of T rec ≃ 3000 K.
The focus is on the implementation of the model into the cosmological scenario using phenomenological values for the drift coefficient J and the associated diffusion coefficient D. As in a more general model with frequency-and time-dependent transport coefficients, these are related to the equilibrium temperature through a fluctuation-dissipation relation T = − lim t→∞ D(ν, t)/J(ν, t), thus constraining the value of the drift once the diffusion coefficient has been determined, and vice versa.The thermalization timescale is τ eq ∝ D/J 2 , and the proportionality factor will eventually have to be derived from astrophysical input, or could probably be measured in a laboratory.So far, no upper or lower limits are known.In this work, the drift coefficient J will be estimated on phenomenological grounds, and the diffusion coefficient is computed from the fluctuation-dissipation relation with the equilibrium temperature T.
The nonlinear boson diffusion model is adapted to the cosmological scenario in the next section, and the analytic solution in case of constant transport coefficients is reconsidered.In particular, the initial conditions for the specific case of Lyα thermalization in the early universe are incorporated into the analytical solution scheme.In Sect.3, the time-dependent results of the thermalization problem for the Lyα initial conditions are presented, and the relation of the transport coefficients to the equilibration timescale is discussed.As a model calculation that does not yet reflect the physically realistic situation in cosmology, the case of complete thermalization is investigated in Sect.4, where it is shown that the solutions of the NBDE correctly approach the Bose-Einstein limit for large times.In Sect.5, the case of incomplete thermalization is discussed that corresponds to the actual time evolution of the Lyα line from recombination in cosmology.An upper limit for the effect of the partially thermalized Lyα line from recombination on the CMB is calculated.The conclusions are drawn in Sect.6. Partial Lyα thermalization in cosmology 2 Nonlinear boson diffusion model in cosmology Planck's equilibrium spectrum for the specific intensity (spectral radiance) as function of frequency ν at temperature T is Due to the expansion of the universe the equilibrium temperature T decreases, but the thermal spectrum is maintained, because both temperature and frequency are reduced with redshift as (1 + z) such that ν/T is unchanged.The chemical potential µ may initially be smaller than zero, but is driven towards zero in the course of time and hence, does not appear in the thermal spectrum.
The cosmic background radiation at early times as well as the CMB radiation at present are therefore modelled as blackbody spectra with µ = 0, as in eq. ( 1).At the time of recombination τ rec ≃ 380 ky, the CBR intensity for an average recombination temperature of T = T rec ≃ 3000 K is shown in Fig. 1 together with the Lyman-α line at ν α = 2466 THz in the vacuum ultraviolet (VUV) region of the electromagnetic spectrum.This line is a doublet with transition frequencies 2466.071THz and 2466.060THz [18] corresponding to the respective 2p 3/2 → 1s 1/2 and 2p 1/2 → 1s 1/2 transitions in neutral hydrogen having Lorentzian line profiles.At the recombination temperature, however, the natural line widths are thermally broadened, such that a single Gaussian profile of line width Γ α ≃ 1.8 T rec ≃ 62.5 THz can be used to represent the occupation-number distribution of the Lyman-alpha line with normalization N α in THz, and standard deviation σ α = Γ α / √ 8 ln 2 = 1.8 T rec / √ 8 ln 2 ≃ 26.8 THz.The specific line intensity at recombination L α (ν, T rec ) = 2hν 3 /c 2 × n α (ν, T rec ) normalized to Rayleigh-Jeans is shown in Fig. 1.A physically realistic normalization will eventually have to rely on observational data, but in Sec. 5, an upper-limit estimate will be given.In the following, a system of units h/(2π) ≡ ℏ = c = k B = 1 will be used.
The position of a redshifted Lyα line in the CMB spectrum would be at about 2240 GHz, in the far Wien end.It is, however, likely that the line is partially thermalized until today through random scatterings with the remaining free electrons (few parts in 10 4 after recombination), resonance scattering, and other processes.If the timescale for thermalization turns out to be comparable to, or larger than the inverse expansion rate, the expansion during Partial Lyα thermalization in cosmology thermalization has to be taken into account.Since Lyα thermalization is likely not completed until the present time, a far-infrared nonthermal background (FIRB) could survive in the CMB above several hundred GHz.However, other nonthermal sources may also contribute there, such as light generated by star formation that is absorbed by interstellar dust and re-emitted in the far infrared.Indeed, FIRB radiation had been detected by COBE's FIRAS [19] and DIRBE [20] instruments, but it is difficult to disentangle its various possible sources.
Once the diffusion function D(ν, t) and the drift function J(ν, t) that account for thermalization are known, the time-dependence of the photonic single-particle occupation number distribution n(ν, t) from its initial distribution eq.( 2) towards the equilibrium distribution n ∞ (ν) can be calculated from solutions of the nonlinear boson diffusion equation [13,14] The mean occupation number n(ν, t) is equal to the mean energy, divided by the energy per photon.This nonlinear diffusion equation has been derived for bosonic systems in [13] to account for the fast thermalization of gluons in relativistic heavy-ion collisions, but it also accounts for thermalization in other nonequilibrium Bose systems, such as cold atoms, or photons.To obtain the NBDE, the quantum Boltzmann equation is first rewritten in form of a master equation, and the discrete transition probabilities between quantum states are expressed as integrals by introducing the corresponding densities of states.An approximation to the master equation is then obtained through a gradient expansion in energy space, and drift and diffusion coefficients are introduced as first and second moments of the transition probabilities, respectively, to finally arrive at the above nonlinear eq. ( 3).The drift term J(ν, t) in the NBDE is negative.It is mainly responsible for dissipative effects such as recoil [21] that drive the distribution towards lower frequencies and cause boson (photon) enhancement, the diffusion term D(ν, t) accounts via the fluctuation-dissipation theorem for the diffusion in frequency (energy) space.The derivative-term of the diffusion coefficient is required such that the stationary solution n ∞ (ν) becomes a Bose-Einstein equilibrium distribution.This can be seen by rewriting the equation, setting the time derivative to zero and solving for n ∞ .
With the condition that the ratio J(ν, t)/D(ν, t) must have no frequency (energy) dependence for t → ∞ such that lim t→∞ [−J(ν, t)/D(ν, t)] ≡ 1/T , it can be shown [15] that the stationary distribution n ∞ equals the Bose-Einstein equilibrium distribution n eq , respectively, Here, the chemical potential µ ≤ 0 appears as a parameter.In Fig. 2, the equilibrium distribution for µ = 0 is shown at recombination together with the Rayleigh-Jeans distribution and the thermally broadened Lyα line.
The nonlinear diffusion equation for the occupation-number distribution n(ν, t) becomes particularly simple for frequency-independent transport coefficients where the derivative-term of the diffusion coefficient is now absent, and the transport coefficients have been pulled in front of the derivatives.
The equation with constant transport coefficients differs significantly from a linear Fokker-Planck equation -which has the Maxwell-Boltzmann distribution as stationary solution -due to the nonlinear term: It preserves the essential features of Bose-Einstein statistics that are contained in the quantum Boltzmann equation.This refers especially to the Bose enhancement in the lowfrequency region that increases rapidly with time.Indeed, for ultracold bosonic atoms it has been shown in [14,16] that the simplified equation with constant transport coefficients -together with the requirement of particle-number conservation -already accounts for time-dependent condensate formation in agreement with recent Cambridge data [22].At much higher energies and temperatures, the NBDE has been used in [14,15] to account for the fast thermalization of gluons in relativistic heavy-ion collisions at energies reached at CERN's large hadron collider (LHC).
The diffusion equation with constant coefficients can be solved in closed form for any given initial condition n 0 (ν) using the nonlinear transformation outlined in [13,14].The resulting exact solution of the NBDE can be expressed as where the generalized (time-dependent) partition function Z(ν, t) obeys a linear diffusion (heat) equation The time-dependent partition function can be written as an integral over Green's function of the above eq.( 7) and an exponential function F (x) which depends on the initial occupation-number distribution n 0 according to x n 0 (y) dy .(8) Partial Lyα thermalization in cosmology Here, the integration constant can be omitted, because it will drop out once the logarithmic derivative of the partition function is taken.The time-dependent partition function with boundary conditions at the singularity becomes For sufficiently simple initial conditions, it can be calculated analytically, as has been done in [23] for ultracold atoms.In eq. ( 9), Green's function G bound accounts for the boundary conditions at the singularity 2πν = µ = 0.They can be expressed as lim ν↓0 n (ν, t) = ∞ ∀ t.One obtains a vanishing partition function -corresponding to an infinite occupation-number distribution -at the boundary, Z bound (ν = 0, t) = 0, and the energy range is restricted to ν ≥ 0. This requires a Green's function that equals zero at ν = 0 ∀ t.It can be written as with the free Green's function G ≡ G free of the linear diffusion eq. ( 7), Finally, the occupation-number distribution is obtained via the basic nonlinear transformation, eq. ( 6).

Time-dependent calculations
To account for the time-dependent -partial or complete -thermalization of the Lyα energy content subsequent to recombination, the integral over the initial distribution of the Lyα line is obtained as which is inserted into the exponential function F (x) in eq. ( 8).With Green's function from eq. ( 10), the partition function is obtained from eq. ( 9), and the time-dependent occupation-number distributions can be computed from the nonlinear transformation, eq. ( 6).
For a realistic calculation of the time-dependent thermalization, the values of the transport coefficients are decisive.They are related to the equilibrium temperature through a fluctuation-dissipation relation T = −D/J.Moreover, the equilibration time scale is related to the transport coefficients [14] according to τ eq = a τ D/J 2 with a proportionality constant a τ .Hence, the transport coefficients are obtained as The proportionality factor a τ for thermalization in a Bose system depends significantly on the initial condition.So far, it has been calculated analytically only for a θ-function initial distribution that overlaps with the low-frequency branch of the thermal distribution [13].No derivation is available for the present case, where a narrow initial distribution at the UV side of the spectrum is far away from the thermal Bose enhancement at the opposite side of the spectrum, and hence, a τ and τ eq are treated as parameters that have to be determined in the cosmological context.Alternatively, the drift coefficient is taken as a parameter -see the next section -, and the diffusion coefficient is computed from the fluctuation-dissipation relation at temperature T .
With the above connections between the transport coefficients D, J, the equilibrium temperature T , and the equilibration time scale τ eq , the nonlinear diffusion eq. ( 5) in the limit of constant coefficients can also be expressed using the dimensionless time-like variable δ Here, the time variable t of eq. ( 5) has been replaced, t → t 0 + τ eq δ.In this dimensionless form of the nonlinear diffusion equation, the temperature T and the equilibration timescale appear as parameters instead of the transport coefficients J and D, together with the dimensionless constant a τ that characterises thermalization in a Bose system for a given initial condition.The solution method outlined above for the NBDE with constant transport coefficients can be applied to eq. ( 15) as well.

Complete thermalization
The analytic solutions of the nonlinear diffusion equation are first applied to an idealized situation of complete thermalization of the Lyα line during the time evolution.This is, of course, not realistic when accounting for the physics of recombination and the subsequent partial frequency redistribution [12] of the photons that are emitted following hydrogen and helium recombination: Scattering is known to hardly be able to thermalize the Ly-alpha distortion even at earlier times.Nevertheless, this calculation serves to demonstrate the method, and it will subsequently be adapted to the actual physical situation of incomplete thermalization of the Lyα-line in the next section.The value of the drift coefficient in this schematic calculation is chosen as J = −1 THz/ky.The frequency shift of the initial Lyα line at short time intervals ∆t ≪ τ eq is approximately ∆ν α ≃ J∆t THz.It later becomes a nonlinear function of time, especially when Bose enhancement sets in at smaller frequencies.Physically, the drift is a consequence of several effects that cause a frequency redistribution towards lower energies such as atomic recoil [21] and electron scattering during recombination.Neglecting expansion and cooling for the moment (it will be discussed in the next section), the equilibrium temperature is kept at T = T rec , and due to the fluctuation-dissipation relation D = −T /J the corresponding diffusion coefficient becomes D ≃ 63 THz 2 /ky.
Results of the time-dependent thermalization with the above parameters are shown for eight timesteps in Fig. 3, lower frame.The Lyα emissions are first broadened and shifted to the low-frequency region, until Bose's enhancement sets in.Equilibrium is reached at t = τ eq ≃ 4 × 10 3 ky in this particular example.Dividing both transport coefficients by a factor of 4 (upper frame) retains the same equilibrium temperature, but thermalization occurs more slowly: It would take 4 × τ eq to reach equilibrium, because τ eq ∝ D/J 2 .In case of complete thermalization, no remainder of the Lyα emissions from recombination would survive in today's CMB spectrum.With the consideration of expansion and cooling, and more realistic values for the transport parameters from astrophysical arguments as will be discussed below, the time dependence will differ, and the occupation-number distribution will remain far from equilibrium, but the principal effects of the approach to equilibrium shall persevere.The result of the lower frame in Fig. 3 is shown again for eight timesteps in Fig. 4 in a double-log plot, which emphasizes the approach to Bose-Einstein equilibrium in the near-infrared region.It can also be seen that even at short times, the solutions of the NBDE generate a low-frequency branch (shortdashed curves on the left) that thermalizes very quickly.It is due to the Bose enhancement that is contained in the NBDE, and leads to a Rayleigh-Jeanslike slope in the infrared when calculating the specific intensity, see the next section.Beyond the equilibration time scale, the full distribution function of the specific intensity becomes again a Planck spectrum, eq. ( 1).  3 (increasing dash lengths), but in a double-log plot to illustrate the approach towards the thermal occupation-number distribution in the visible and infrared region of the spectrum.

Incomplete thermalization
Whereas the above results show the general viability of the nonlinear diffusion model to account for thermalization in a bosonic (here: photonic) system, the actual Lyα physics in the course of recombination remains far from equilibrium.One important reason is the rapidly falling free-electron number density during recombination, thus diminishing scattering processes that are required for thermalization: At the end of the recombination era at redshift z ≃ 500 corresponding to an equilibrium temperature of T = T f ≃ 0.45 T rec , the relative free-electron density has dropped below 10 −3 [12], the values of the transport coefficients in the NBDE must be reduced accordingly, and their time dependence during recombination has to be considered.
The partial frequency redistribution of Lyα photons has been treated in [24] based on a Fokker-Planck approximation -which may, however, not be sufficient towards the end of hydrogen recombination.That numerical approach had been proposed by Rybicki [25], who had also discussed a correspondence to Kompaneets' equation [26] when written in terms of the photon occupation number.This equation concentrates, in particular, on the role of the Compton effect in the establishment of equilibrium between quanta and electrons in a nonrelativistic approximation.
As a complement, the nonlinear diffusion equation offers a related solution to the problem of partial thermalization that properly accounts for Bose statistics, takes into account the boundary conditions, and provides a transparent analytical solution through a suitable nonlinear transformation.
First-principles calculations of drift and diffusion coefficients in the NBDE based on the relevant physical processes electron scattering, and resonance scattering off moving atoms -both including recoil, Doppler broadening and induced scatterings -are not yet available in the cosmological context for the proposed NBDE-model.The transport coefficients are instead estimated here on phenomenological grounds.We start from the equilibrium temperature at the beginning of the recombination era, T i = T rec ≃ 3000 K at redshift z ≃ 1100.The end of recombination is taken to be reached at T f ≃ 0.45 T rec ≃ 1350 K at redshift z ≃ 500.Using the Planck cosmological parameters H 0 = (67.4±0.5)km/s/Mpc, Ω m = 0.315±0.007[9], these values correspond to time scales [27] in the ΛCDM model.In the redshift region 1100 ≳ z ≳ 500 where partial thermalization is acting in the schematic nonlinear model, the relation between redshift and time can approximately be expressed as an exponential function.
Since the temperature depends linearly on the redshift, T (z) = T CMB (1+z), a related redshift-or time dependence is taken for the transport coefficients.For an equilibrium temperature at the end of recombination T f ≃ 1350 K at redshift z ≃ 500 and an equilibration time scale τ eq ≃ 1.5 × 10 3 ky, this corresponds to exponential time dependencies, with J 0 = −1 THz/ky as in the previous model calculation, and D 0 = −0.45T rec /J 0 = 28 THz 2 /ky.The ratio −D(t)/J(t) = T f remains time independent, as required.The value of the drift may turn out to be somewhat overestimated when compared with detailed numerical simulations of the frequency redistribution such as in [24], but could be adapted accordingly.The value of D(t) is, however, computed from the fluctuation-dissipation relation, which is inherent to the present model.It could only be modified by permitting frequency-dependent transport coefficients, as in eq. ( 3).The result of the time-dependent calculation with the above parameter set is shown in Fig. 5 in four timesteps up to t = 1.5 × 10 3 ky, corresponding to redshift z ≃ 500.The distribution functions become already asymmetric (note the log-scale), indicating that a linear Fokker-Planck approximation may not be justified.They remain, however, far from equilibrium -except for the lowfrequency branch, which thermalizes quickly even for short times, and leads to a Rayleigh-Jeans slope in the specific intensity of the partially thermalized Lyα line, see below.Increasing the equilibration time τ eq in eqs.( 16), (17) artificially by a factor of two causes larger values for drift and diffusion at any given time and thus, somewhat faster partial thermalization (the distribution for t = 1.5 × 10 3 ky peaks at 1290 THz rather than 1700 THz), thus decreasing the probability that a remainder of the Lyα-line could be visible in the CMB.Assuming that thermalization indeed terminates at the end of the recombination epoch, the effect of the Lyα photons from hydrogen recombinationwhich represent most of the photons from both, hydrogen and helium (≃ 24%) recombination -on the CMB distribution based on the NBDE evolution can be estimated by propagating the distribution function at z = 500 taken from Fig. 5 (solid yellow curve) to z = 0.The normalization is taken according to the ratio of photons in the isotropic blackbody radiation to baryons -mostly protons -, which is approximately 1.6 × 10 9 .This yields an upper limit to the effect of the partially thermalized hydrogen Lyα line on distorting today's CMB spectrum.
The result can be seen in Fig. 6, where the blackbody CMB (solid curve) including the COBE-FIRAS data [6] is shown together with the renormalized solution of the NBDE at z = 500 from Fig. 5, propagated to z = 0 (dashed curve).Its peak resides on the Wien side of the CMB, and the amplitude is more than seven orders of magnitude below the one of the CMB.To test the sensitivity of the nonlinear model, I have reduced the amplitudes J 0 of the drift coefficient and D 0 of the diffusion coefficient by factors of two, thus keeping the equilibrium temperature at the same value as before.The result is shown as a dotted curve in Fig. 6, which is still below the CMB signal by almost seven orders of magnitude.This largely analytical calculation gives a clear hint why no Lyα signal from recombination -or oscillatory signal when taking into account all other radiative transitions in hydrogen and helium -is visible at the present level of precision in the CMB, although the distortions can lead to biases to several cosmological parameters [9].
As shown in Fig. 6, the NBDE result yields a Rayleigh-Jeans slope in the intensity of the Lyα line at low frequencies.This is a consequence of the Bose enhancement, which would not occur in a linear Fokker-Planck type approach to partial thermalization.To reach the proper Wien limit also at large frequencies ν > 2000 GHz, energy-dependent transport coefficients would be required -which is beyond the scope of an analytical model.

Conclusions
The time-dependent incomplete thermalization of the Lyα line that is emitted from neutral hydrogen atoms during recombination has been accounted for in a nonlinear diffusion model.This approach is complementary to the available detailed numerical treatments of the release of photons during the recombination epoch, their partial frequency redistribution, its effect on the recombination history, and possible observable distortions of the CMB.
In the analytical model, the thermally broadened Lyα emission line at an average recombination temperature of 3000 K provides the initial condition.With the proper boundary condition that causes the low-frequency Bose enhancement, the diffusion equation is solved through a nonlinear transformation in the limit of constant transport coefficients.The stationary solution is equal to the Bose-Einstein equilibrium distribution.
As is well known, the system remains far from equilibrium, because the interaction of the radiation field with the electrons can not transform a non-Planckian spectrum into a Planckian one in the course of, or after, recombination.However, the low-frequency Rayleigh-Jeans slope in the specific intensity indicative for thermalization correctly emerges already at very short times from the analytical solutions of the nonlinear diffusion equation -which would not be the case in a linear theory for the frequency redistribution of Lyα, or other recombination lines.Variations of the transport coefficients by a factor of two have shown that the model gives robust results.Moreover, it is not overly sensitive to modifications in the time-dependence of the transport coefficients, thus enhancing its reliability in the astrophysical and cosmological context.
In the present work, I have investigated partial Lyα thermalization in a redshift range 1100 ≳ z ≳ 500.Additional thermalization may result at later times, in particular, during the epoch of reionisation, when the ultraviolet light from the first stars at t ≃ 370 My and redshift z ≃ 12 re-ionizes hydrogen and helium until about z ≃ 6.This may be another interesting topic for the application of the nonlinear diffusion model.
Regarding possible signatures of the redshifted Lyα recombination line in today's CMB, an upper limit of the specific intensity following partial frequency redistribution as calculated from the nonlinear diffusion model is estimated to be about seven orders of magnitude smaller than the CMB signal, and therefore unlikely to be directly detectable at present.This is in line with other numerical calculations, and also with Planck observations of the CMB, which at the present accuracy do not exhibit a frequency-modulated signal from the recombination spectrum.Here this result has been obtained in a novel nonlinear diffusion model that is analytically solvable and offers a transparent description of the partial thermalization process.The model is also expected to be useful regarding other astrophysical or cosmological equilibration processes in the non-thermal universe.An example is the reheating following inflation, which could be modelled using the nonlinear diffusion equation in case the timescale turned out to be finite.

Not applicable • Availability of data and materials
The COBE-FIRAS data that are compared with my model calculations in Fig. 6

Fig. 1 Fig. 2
Fig. 1 Planck spectrum of blackbody radiation at recombination with an average temperature of Trec = 3000 K at redshift zrec ≃ 1100, and the thermally broadened Lyman-α line of neutral hydrogen at να = 2466 THz in the VUV region of the spectrum.The dot-dashed line is the Rayleigh-Jeans distribution, L RJ (ν, T ) ∝ ν 2 Trec.

Fig. 4
Fig.4Time-dependent evolution of the Lyα line in the ultraviolet part of the spectrum (blue, right at 2466 THz) towards equilibrium (solid, left) assuming complete thermalization towards T = Trec and time-independent D, J. Eight timesteps are shown as in Fig.3(increasing dash lengths), but in a double-log plot to illustrate the approach towards the thermal occupation-number distribution in the visible and infrared region of the spectrum.

Fig. 5
Fig.5Nonequilibrium evolution of the Lyα line (blue, right at 2466 THz): Incomplete thermalization with time-dependent transport coefficients D(t), J(t), see text.Timesteps t/10 3 ky= 0.05; 0.2; 0.5; 1.5 are shown, with no significant thermalization expected after the last step, which corresponds to redshift z ≃ 500.The curve at the lower left is the equilibrium distribution for T = 0.45 Trec ≃ 1350 K, which is reached only at low frequencies.

Fig. 6
Fig.6Today's CMB spectrum of blackbody radiation with T=2.725 K (solid curve), the Rayleigh-Jeans distribution (green dot-dashed line), and the partially thermalized NBDEsolution for the redistributed, and redshifted Lyα line featuring a low-frequency R-J branch due to Bose enhancement (dashed curve).The dotted curve is the NBDE-result for a 50% reduction in the drift and diffusion coefficients.See text for the cutoff at ν ≃ 2000 GHz.COBE-FIRAS data[6] are shown as blue dots, error bars are smaller than the symbol size.
were provided by John C. Mather.The model results shown in Figs.1-6 (Mathematica-plots) are available as tables from the corresponding author upon reasonable request.• Code availability The Mathematica-codes for the calculations shown in Figs.1-6 can be made available upon request.• Authors' contributions Single author GW of the manuscript.