Spectro-spatial evolution of the CMB. Part II. Generalised Boltzmann hierarchy

In this paper, we formulate a generalised photon Boltzmann hierarchy that allows us to model the evolution and creation of spectral distortion anisotropies caused by energy release in the early Universe. We directly build on our first paper in this series, extending the thermalisation Green's function treatment to the anisotropic case. We show that the problem can be described with the common Boltzmann hierarchy for the photon field extended by new spectral parameters — a step that reduces the complexity of the calculation by at least two orders of magnitude. Our formalism describes the effects of i) Doppler and potential driving, ii) spectral evolution by Compton scattering, iii) perturbed thermalisation and iv) anisotropic heating on the distortion anisotropies. We highlight some of the main physical properties of the equations and also outline the steps for computing CMB power spectra including distortion anisotropies. Limitations and extensions of the formulation are also briefly discussed. The novel Boltzmann hierarchy given here is the basis for a series of companion papers studying how distortion anisotropies evolve in the perturbed Universe and which physical processes could be constrained using future CMB imaging techniques.


JCAP11(2023)027 1 Introduction
The cosmic microwave background (CMB) has been a goldmine for furthering our understanding of the cosmos [1][2][3][4][5].We are now entering a new phase, in which novel cosmological observables are moving to the focus of our efforts [6][7][8][9].One of these observables is spectral distortions (SDs) of the CMB [10,11].In this work, we develop a novel Boltzmann formulation that can describe the evolution of SD anisotropies created by energy release in the early Universe.We base the treatment on our first paper in the series [12, henceforth referred to as paper I], which introduced a novel discretisation of the thermalisation Green's function to efficiently describe the evolution of average SDs.In paper I, we were limited to running one single thermalisation history at the time, effectively obtaining a solution for the average distortion intensity spectrum, ∆I ν (t), at a given time t and frequency ν.In simple words, we now deliver the tools to repeat the calculation along multiple lines of sight and including the effect of perturbations on the CMB signal evolution.This will open the way to studying SD physics using standard methods known from CMB imaging, without the need for absolute calibration, delivering new targets for experiments like Litebird [13], PICO [14] and CMB-S4 [6].
We refer the reader to paper I for more introduction and motivation to the topic.The steps taken here explain in detail how to extend the thermalisation Green's function treatment to the anisotropic case.We are in particular keen on including thermalisation effects to the evolution.This captures the full spectral evolution due to Compton scattering, double Compton (DC) and Bremsstrahlung (BR), which are so crucial in the formation of SDs at redshifts z 10 4 [15][16][17][18][19].These effects can change the type of the distortion and gradually rotate y-type distortion sources into µ and ultimately thermalise the distortion completely.The rotation is with respect to the energy in various spectral components, which by construction is conserved across the spectral basis.A perturbative formulation that includes some of the Doppler boosting effects has been given in [20] and [21,22]; however, this does not capture the effects of repeated Compton scattering at z > 10 5 and also does not describe the conversion of µ to T , thus having limited applicability, which we overcome here.
In section 2, we start by recapping the main outcomes of paper I.We then move on to describing some of the distortion effects in the Thomson limit of the kinetic equation.This only captures the effect of Doppler and potential driving as well as free-streaming mixing on the SD signals across the sky, but provides valuable insight preparing for the general case.In section 2.4, we add thermalisation terms to complete our treatment.Details of the derivation are giving in appendix A to C, which we recommend for in depth reading.The tools for computing the CMB power spectra for various combinations of temperature and SD parameters are given in section 2.8, and a discussion of the basic expectations is presented in section 2.9.Numerical solutions and first forecasts will be presented in paper III [23].In section 3, we present our conclusions and also highlight some of the limitation of the formulation and further work that may become important.

JCAP11(2023)027 2 Extended Boltzmann hierarchy with primordial spectral distortions
The goal of this section is the provide an approximate description of the spectro-spatial thermalisation problem at first order in perturbation theory.This greatly generalizes previous treatments of the problem allowing us to capture the main sources of distortions in the presence of perturbations.Unlike the usual split, here we have two small parameters to work with, one for the average energy release, γ = ∆ρ γ /ρ γ , and one for the primordial curvature perturbations, ζ.We keep terms up to O( γ ζ), thereby allowing to capture all linear order (spatial) perturbation effects.
We will use the approach presented in paper I to model the evolution of the local monopole across the chosen spectral basis, while the spatial evolution is treated using the standard multipole decomposition which includes the effect of Thomson scattering and free-streaming.We will make direct use of the main results from paper I and refer the interested reader to that paper for clarification of the notation and details of the Green's function discretisation and choice of basis functions.

Paper I and zeroth order problem
In paper I, we essentially solved the zeroth order photon equation for the evolution of the average CMB spectrum.This problem can be described with the kinetic equation [see section 4.1.1 of 24]: for the evolution of the photon occupation number distortion ∆n (0) ≡ ∆n (0) (t, x), defined with respect to the average blackbody, n bb (x = hν/k B T z ) = 1/(e x − 1) at a temperature T z .
Here, the dot denotes time derivatives and we introduced the Thomson optical depth, τ = σ T N e c dt, which is evaluated at the background level assuming the standard recombination history from CosmoRec [25].The temperature variables are presented as θ i = k B T i /m e c 2 , with θ z ∝ T z ∝ (1 + z).The Kompaneets operator is denoted by Kx = x −2 ∂ x x 4 ∂ x + x −2 ∂ x x 4 A(x), with A(x) = 1 + 2n bb (x), and Λ = Λ(x, θ z ) determines the photon production rate by double Compton (DC) and Bremsstrahlung (BR).These can be computed accurately using BRpack [26] and DCpack [27].The electron temperature perturbation, Θ (0) eq = ∆T (0) e /T z , and effective heating rate, Q(0) , are given by Θ (i)  eq = x 3 ∆n with y-weight factor w y = Y /G = xA(x) − 4 = x e x +1 e x −1 − 4 and where E f = x 3 f (x) dx is the energy density integral of f (x).The expression for Θ (i) eq can be obtained by balancing Compton heating and cooling (appendix C.1).The distortion sources from heating relate directly to Q(0)

JCAP11(2023)027
Usually, eq.(2.1) is solved on a frequency-grid for the spectral distortion, ∆n (0) (t, x), e.g., using CosmoTherm [28].This task can become very time-consuming and certainly is not easily extendable to anisotropic distortions.The main result of paper I was to demonstrate that the problem can be approximately written as a matrix equation.The derivation uses the Ansatz ∆n (0) (t, x) ≈ B(x) • y (0) (t) to discretise the average photon spectrum.Here, B = (G(x), Y (x), Y 1 (x), . . ., Y N (n), M (x)) T denotes the computation basis, based on the standard distortion shapes, G, Y and M as well as the boosted signals, Y k = (1/4) k Ôk x Y with boost generator, Ôx = −x∂ x .This then yields the evolution equation 4 , for the spectral parameters, y (0) (t).Here, x c is the critical frequency for photon emission, and the emission coefficients are γ T ≈ 0.1387 and γ N ≈ 0.7769.The Kompaneets mixing matrix, M K describes the (photon number-and energy-conserving) rotation of the spectral parameters in each scattering, and is pre-computed for different basis sizes.From eq. ( 2.3) one can already anticipate the main thermalisation terms; however, a few augmentations will become important as we show now.We also note, that for eq.( 2.3), the problem was linearised respect to the distortion.Possible sources of distortions from photon injection processes [e.g., 29,30] were also neglected and cannot be readily treated within the proposed framework.

JCAP11(2023)027
are added back as usual [31].Equation (2.5) will now be generalized to include the effect of a non-vanishing average distortion, spectral evolution, anisotropic heating and perturbed thermalisation.These effects lead to small corrections to the brightness temperature equations but give the leading order terms in the distortion hierarchy.
Accounting for all the thermalisation effects self-consistently is beyond the scope of this paper.However, if we assume that locally thermalisation is only mediated through the monopole part of the spectrum, the situation becomes more tractable.This is in fact wellmotivated since for the anisotropies the much faster Thomson scattering process dominates, while Thomson terms are absent in the monopole, making thermalisation terms dominant there [24,36].In particular the monopole here is defined in the local inertial frame.We note, however, that changing frame would introduce similar transformations in the Thomson and Compton scattering terms, and the arguments of timescales would still hold true for the relevant pairings of non-monopole scattering terms.We will furthermore not attempt to solve for higher order corrections to the CMB temperature anisotropy field, as this would entail a full treatment in second order perturbation theory [e.g., 35,38,39] including additional modifications to more correctly treat distortion terms.
To better appreciate the various new effects we will proceed in a step-by-step manner first only considering Thomson scattering terms, ∝ τ .These will reveal the optimal basis given an average distortion and also already show the main effects in terms of distortion anisotropy generation.The spectral evolution across the various distortion types is then included using the modified Green's function treatment with the spectral mixing occurring only in the local monopole.The scattering efficiency is ∝ τ θ z and hence suppressed relative to the Thomson terms.However, at z 10 4 , we expected spectral evolution to change the distortion signals in an observable way.
We give many of the derivation details in appendix A. In addition, we will consider various limiting solution in appendix B. In appendix C, we furthermore explicitly carry out some of the approximations that are required to simplify the problem to the expressions given in the main text.We refer the interested reader to these sections for an in depth understanding.However, we hope that the result presented in this section can be mostly understood without the detailed (and sometimes cumbersome) derivations.
We also mention that the picture and interpretation of the results we develop below could potentially be cast into another form using alternative variables and reference frames.For example, one could use the local baryon frame to overcome some of the kinematic complications in the scattering process [36].Similarly, log-temperature moments [40] might provide the means to separate some of the redshifting and scattering effects in a more transparent way.However, the methods applied here closely follow the standard approach for CMB anisotropy computations [31,41], extending them in a minimal way to account for the evolution of the average spectrum.This is deemed to be a valuable step in understanding the rich physics of CMB spectral distortion anisotropies.

Effect of Doppler and potential terms in the Thomson limit
As the first and simplest scenario, let us consider the case where heating happened well before the recombination era and the average distortion is long frozen into its final state.We shall JCAP11(2023)027 also assume that any thermalisation effects and energy exchange terms can be neglected.Thomson scattering isotropises the medium and leads to damping of the perturbations at small scales.This picture is valid for modes that become dynamic at z 10 4 .The photon evolution equation then reads ∂n (1)  ∂t + cγ a • ∇n (1) + Ôx n (0) ∂Φ (1)  ∂t 2 − n (1) + β (1) χ Ôx n (0) (2.6) using the Thomson collision term, eq.(A.1).
Here, we defined n (t, x, r, γ) = m n m (t, x, r)Y m (γ) using the spherical harmonic coefficients of the photon occupation number, n m (t, x, r).This definition implies n (0) ≡ n , where the solid angle is dΩ = d 2 γ = dφ dχ.We now assume that the average spectrum is spectrally frozen, n (0) (x) = n bb (x) + ∆n (0) (x), where ∆n (0) (x) is the departure from the blackbody, n bb (x) = 1/[e x − 1] at a temperature T z .We can then make the simple Ansatz n (1) ≈ Θ (1) G(x) + Σ (1)  Ôx ∆n (0) . 1 Assuming that Ôx ∆n (0) = G(x), this results in two identical photon Boltzmann hierarchies for Θ (1) and Σ (1) that are exactly like the standard equation for the temperature perturbations.These can be obtained by simply comparing coefficients of the two independent spectral functions or alternatively by taking the number and energy density moments of the Boltzmann equation (see appendix B.1).However, in contrast to the temperature perturbations, at a given wavenumber k we start with the initial condition Σ (1) = 0 such that distortion anisotropies are only sourced later when perturbations in the potentials and velocity field appear [24].The transfer functions and related power spectra can be computed directly using any Einstein-Boltzmann code.Cross-correlations of the distortion and temperature fields thus probe Doppler and integrated Sachs-Wolfe (ISW) terms.This can result in a novel noise floor for tests of primordial non-Gaussianity from µ/y × T correlations, as we highlight in paper III.

Preparing for more general spectral evolution
As explained above, in the Thomson limit the problem described by eq.(2.6) can in principle be solved with only two independent variables, Θ (1) and Σ (1) [see eq.(B.4) for the Boltzmann equations], and the new spectral shape Ôx ∆n (0) caused by boosts 2 of the average distortion.For example, if the average spectral distortion was a pure µ-distortion, the anisotropies would have the spectrum of Ôx M (x).We could then observationally search for G(x) and Ôx M (x) and thereby extract information on the perturbations, Θ (1) ∝ ζ and Σ (1) In more general situations, when spectral evolution is also included, using the computational spectral basis we can write ∆n (0) ≈ y (0) • B and hence Ôx M (x), (2.7) JCAP11(2023)027 where we used Ôx G = 3G+Y .The first two groups of terms directly fall back onto the original spectral basis, B; however, Ôx Y N (x) = 4Y N +1 (x) and Ôx M (x) = [M (x)/G(x)] Ôx G(x) − G(x)/x lie outside.As a simple fix, we could add these new spectral shapes to the basis and thereby keep the precision of the average distortion.However, once we consider Compton scattering effects, this approach becomes problematic, requiring a new truncation of the spectral hierarchy.
Instead of adding new spectra to the basis, we approximately represented Ôx Y N (x) and Ôx M (x) within the old spectral basis. 3The projection procedure is explained in detail in paper I. Just like for the Kompaneets operator, one has to ensure energy conservation in the calculation; photon number conservation is automatically built in, since Y N and M do not carry number and hence Ôx Y N and Ôx M do not either.It turns out that Ôx M (x) is extremely well represented as a simple sum of µ and y.Carrying out the projections, we find Ôx M (x) ≈ 0.3736 Y (x) + 1.9069 M (x); however, even better precision can be achieved when including additional terms in Y k (x) (see figure 1).As for the Kompaneets operator (see paper I), the largest error is always made when representing Y N +1 (x) with only terms up to Y N (x), as part of function space is omitted.However, less and less energy is carried in these contributions such that the precision does not suffer much once Y 5 or higher are included.
With these comments, schematically we can therefore always write the representations

JCAP11(2023)027
where O N and O µ are the corresponding solution vectors in terms of y, y 1 , . . ., y N , µ (the projections of G(x) never matter as Ôx Y k and Ôx M do not carry photon number).This means we can express the boosted average spectrum as Ôx where we implicitly defined the boosting matrix, M B , noting that B has N + 3 entries starting with G(x) for k = 0 and ending with M (x) for k = N + 2. Inserting n (1) = y (1) • B into eq.(2.6) and rearranging terms we can then write the photon Boltzmann equation in the Thomson limit as ∂y (1)  ∂η + γ • ∇y (1) ≈ −b (0) 0 ∂Φ (1)  ∂η + γ • ∇Ψ (1) + τ y ( 2 − y (1) + β (1) χ b where we now converted to conformal time, η = t 0 c dt /a with τ = dτ / dη = N e σ T a, and also used y (t, r, γ) = m y m (t, r)Y m (γ) as before for n .We furthermore made explicit that b (0) only has a monopolar dependence, b (0) = b (0) 0 ≡ b (0) dΩ 4π .For the example above, b (0) 0 was time-independent.In more general situations it becomes a function of time, where ∆b 0 provides sources for distortion and temperature anisotropies.However, before we can treat this problem we need to include the effects of spectral evolution into the equations, as we discuss in section 2.4.

Changing the background temperature
One of the simplest problems that one can already consider with eq.(2.10) is a time-dependent change of the background temperature that leads to a departure from the standard ∝ (1 + z) scaling.This limit describes average injection of G, which physically is hard to establish without efficient thermalisation terms [42], but provides some intuition.Speculative scenarios could relate to the decay of dark energy into photons [43,44] or a treatment of energy release signals in the instantaneous thermalisation limit.It is also relevant to early treatments of temperature fluctuations at second order in perturbations [38], warranting the discussion presented here.
Using the zeroth order equation, eq.(2.1), with ∆n (0) = Θ (0) 0 G, and assuming that the change in the temperature is created by a time-dependent source term, G S T (t), we can easily -7 -

JCAP11(2023)027
verify that, as expected, this leaves the average spectrum unchanged: Here, we used Θ (0) eq = Θ (0) 0 and (1 − e −x )/x = e −x (e x − 1)/x = n bb /G and assumed that the external source has a spectrum G.We also used Kx G = −Y , which causes all scattering terms to vanish, and therefore no new spectral shapes other than G appear in the average evolution.
In the absence of thermalisation effects, this indeed causes distortion anisotropies, as we illustrate now.However, since even a pure blackbody temperature fluctuation has a y-type spectral contribution at second order in Θ [40,45,46], to demonstrate that a distortion anisotropy is created we have to show that the photon number and energy densities in the anisotropies do not simply obey the blackbody relations, with one effective temperature describing the full spectrum.
To make progress, we therefore first ask the question how a spatially-varying blackbody spectrum changes when the average temperature is varied in a way that departs from the standard ∝ (1 + z) scaling.For this, we write the expression where Θ describes the temperature anisotropies while Θ is only time-dependent.At zeroth (no spatial terms) and first order in perturbations (only up to terms ρ ζ) this implies where we neglected terms O( Θ) 2 and O(Θ (1) ) 2 .We comment immediately on the consistency of this limit: terms O( Θ) 2 do not add any spatial effects and thus are merely higher order corrections to the average spectrum.Terms relating to O(Θ (1) ) 2 O(ζ) 2 are second order in the primordial curvature perturbations, which we also do not consider here.

JCAP11(2023)027
The expressions in eq.(2.14) seems to suggest that a y-type distortion is sourced at first order, even if we started with a blackbody.However, the terms shown above are merely needed to precisely transition from a blackbody at the initial temperature to a new blackbody when Θ varies [e.g., 45,46].This begs the question if this Ansatz solves the evolution equations in eq.(2.12), which would show that no real distortion anisotropy is actually created?4Using D t [X] from eq. (B.2a), we have where Ôx Put together with the rest of the Liouville operator [see eq.(B.2) for relevant definitions], we then have where we equated with the collision term, formally including both the Thomson terms, C T [Θ (1) ], and thermalisation effects, C therm [n].We also used n (1) ≈ Θ (1)  Ôx n (0) to factor Ôx n (0) out of C (1) If we now only consider the Thomson terms and integrate photon number and energy density, with L[n (1) ] x 3 dx → L[Θ (1) ] + 4Θ (1) S T ≈ C T [Θ (1) ]. (2.17) Here, we used the number density, of a blackbody at a temperature T z and also neglected higher order terms in Θ.Since both equations have to be fulfilled, this shows that the spectrum cannot be a blackbody anymore without additional thermalisation processes.The reason is that Y (x) does not carry photon number and only contributes to the second equation, leading to an extra source of Θ (1) S T .
The CMB anisotropies are thus distorted.
Since the only source spectra that are present in the Ansatz are G and Y , in the Thomson limit we need an extra independent y-parameter to describe the full spectrum.Adding an extra term y d Y (x) to the Ansatz we previously used [i.e. using n (1) + y (1) d Y (x) as our input], we then have L[n (1) ] = L[Θ (1) ] (1 + 3 Θ) + 3Θ (1) S T G + ΘL[Θ (1) ] + Θ (1) (2.18)

JCAP11(2023)027
Here, we explicitly split the terms ∝ G and ∝ Y , since we now wish to keep track of distortion contributions.By comparing coefficients we can find the modified system L[Θ (1) ] + 3Θ (1) S T ≈ C T [Θ (1) ], (2.19a) ΘL[Θ (1) ] + Θ (1) d ] − β (1) χ, (2.19b) where we again neglected thermalisation effects.With ΘL[Θ (1) ] ≈ Θ C T [Θ (1) ], after dropping terms O( Θ) 2 , one then has the explicit evolution equations for Θ (1) and y (1) d ∂Θ (1)  ∂t + cγ a • ∇Θ (1) + ∂Φ (1)  ∂t 2 − Θ (1) + β (1) χ , (2.20a) These equations show that a non-vanishing y-parameter is created as perturbations mix through Thomson scattering and free streaming in the presence of average CMB temperature changes.As alluded to above, the reason is that the spectrum of the initial temperature perturbations disagrees with that of the evolving average blackbody.The change in the average blackbody temperature sources Θ (1) and y (1) d perturbations but on average no photons are created in the fluctuating part, implying that the perturbed photon field is not consistent with that of pure blackbody temperature fluctuations: one cannot simply transform all y-type terms away by redefining the temperature perturbations.
This discussion also shows that the Ansatz in eq.(2.14) can be recast in terms of the effective parameters Θ (1) eff = Θ (1) (1 + 3 Θ) and y d + ΘΘ (1) .Indeed, using this redefinition with eq.(2.20) we recover eq.(2.12).This clearly separates the origin of the distortion anisotropies, identifying temperature raising y-contributions, y (1) mix = ΘΘ (1) from contributions that cannot thermalise, y (1) d .These terms can in principle be distinguished from simple higher order corrections ∝ (Θ (1) ) 2 or uncertainty ∝ ∆T 0 in the exact average present-day blackbody temperature ∝ ∆T 0 Θ (1) .For the former, this statement is evident since the related terms exhibit a different correlation structure.For the latter, the effect would be time independent, which does change the correlation structure.However, a more in depth discussion of these effects is beyond the scope of this work.
We mention that corrections ∝ Θ ζ to the potentials and baryon velocity equations in principle also need to be added to the Boltzmann hierarchy.For the velocity terms, we will more explicitly discuss this in section 2.6.However, here we are not interested in computing the exact corrections to the temperature power spectra caused by changes of the background temperature, nor will we consider cases when these terms become very important.For computing the physical y-parameter caused by changes in the average CMB temperature, we can simply use the old set of equations to solve for Θ (1) , since all corrections ∝ Θ ζ recursively enter the evolution equation for y (1) at higher order in Θ.In the Thomson limit, we have thus completed the formulation of the problem and demonstrated that real y-type distortion anisotropies are created by changing the average temperature.These could potentially be used to test for departures from the standard CMB temperature-redshift relation; however, we leave a more detailed discussion to future work.

Instantaneous thermalisation
We can extend our discussion to the case where thermalisation is always instantaneous.At the background level we then have a changing temperature according to Θ = S T relative to T z like before.In the anisotropies, photon production and Compton scattering would always ensure that the spectrum of the fluctuating part thermalises under energy conservation.From the arguments leading up to eq. ( 2.17), we can therefore directly write ∂Θ (1)  ∂t + cγ a • ∇Θ (1) + ∂Φ (1)  ∂t 2 − Θ (1) + β (1) χ , ( which with the parameter Θ (1) eff = Θ (1) (1 + 4 Θ) can be cast into the equivalent form eff + b ∂Φ (1)  ∂t eff,0 + bβ (1) χ , The factor of 4 instead of 3 originates from the fact that the energy carried by the y-distortion in the previous case now appears in the blackbody part due to the assumption of instantaneous thermalisation.Equation (2.22) indicates that one can simplify the computation by scaling all variables (also Φ, Ψ and β) to absorb the factor of b = (1 + 4 Θ).However, assuming a general timedependence of Θ, new effects on the CMB power spectra would appear if this equation would be valid throughout the recombination era: the time-dependence of the modes would be modified by the evolution of the background temperature and a time-dependent rescaling of the temperature contrast would be folded into the shape of the CMB power spectra, an effect that can principally be captured with the above equations.Physically, this of course is an academic example to highlight how a coupling between modes and the background can be formulated.Indeed, this effect is an apparent 'super-horizon' effect, which becomes even more obvious when we think about times before BBN, where multiple phase transitions do increase the average blackbody temperature [e.g., 33,47].Due to quasi-instantaneous thermalisation, this leads to a rescaling of the temperature variables at all scales, leaving the adiabatic nature of the initial perturbations totally unchanged and merely modifying the initial conditions to account for this modification.Clearly, for modes entering the horizon during BBN or before, there is no tracer of the time-dependent effects that could be observed today, since all modes that witnessed this effect have dissipated away by Silk damping.At the normal CMB scales, we are left with the standard CMB anisotropy evolution (relative to the present-day higher temperature).However, in the primordial gravitational wave background, as an example, the traces of the phase transitions are in principle still visible [48,49].

Effect of Compton scattering and photon production
Now that we understand how to account for the effect of boosting and Thomson scattering, we will next include Compton scattering and photon production.These only affect the spectral evolution of the local monopole and can be treated using our ODE representation of -11 -JCAP11(2023)027 the Green's function.We will start by a simplistic treatment (section 2.4.1) that just uses a perturbed version of the zeroth order thermalisation terms, eq.(2.3).A more rigorous derivation (section 2.4.2) shows that a few additional terms appear; however, the overall picture does not change crucially.

Simplistic inclusion of thermalisation terms
Starting from the description in eq. ( 2.3), one can obtain a simple version for the thermalisation terms at first order in perturbation theory.Perturbing τ and including the potential perturbation due to the local inertial frame transformation [e.g., following 35], one readily finds ∂y (1)  ∂t therm ≈ τ θ z M K y (1) + D (1) 4 , While the equation above has been obtained in a simplistic manner, it actually turns out that even a more rigorous derivation does not change the result that much.The first group of terms simply describes the spectral evolution of the first order distortion parameters.The second accounts for corrections with respect to the average evolution from perturbations in the electron density and potentials.Here, we used τ (1) / τ ≈ δ (1) b , assuming that we are far from the recombination era, such that perturbed recombination effects [e.g., 35,50,51] can be omitted.The perturbed heating rate similarly includes perturbations in the first order heating term from collisions in the local inertial frame, Q(1) c , but also the effect of potentials on the zeroth order term.
However, the term 3 2 τ θ z Θ (1) 0 D (0) deserves a bit more explanation.It stems from perturbing the critical frequency, x c , with respect to the local photon temperature.We assumed that only DC is relevant for the emission processes, such that the photon emissivity is 0 .This term thus relates to perturbed emission effects and modulates the zeroth order term, D (0) .However, as we show below the coefficient of this term indeed changes to unity when a more careful account for modifications to the local Compton and DC rate is carried out.

More rigorous treatment
In this section, we now obtain the thermalisation terms in a more rigorous manner.For this we have to follow a few steps as outlined in appendix A and C.After writing the full photon collision term in appendix A, one has to obtain the electron temperature in the given distorted radiation field.For this we make use of the result of [35] and include the effect of Compton scattering (appendix C.1).While thermalisation terms are relevant, the electron temperature will always be extremely close to the Compton equilibrium temperature, which JCAP11(2023)027 greatly simplifies the problem and leads to an effective heating term in the photon field as we demonstrate in appendix C.2.The main result of that section is eq.(C.11) for the Compton terms and heating sources.This expression only neglects one stimulated scattering correction, as explained in that section, but otherwise is consistent at order O( ρ ζ).
Translating eq.(C.11) into matrix form using our spectral basis, we find the generalized Kompaneets terms at first order in perturbations ∂y (1)  ∂t where is the average electron temperature.We also used δ 0 for the photon energy density perturbation.Aside from the last group of terms, all the other are already present in eq.(2.23).As explained in appendix C.2, the latter do not add any energy to the system, but merely lead to a more minor change in the spectral distortion evolution.
In appendix C.3 we carry out a careful derivation of the photon production terms.It turns out that at first order in perturbations, a correction to the DC emissivity causes a modification of the related temperature correction term, changing 3  2 τ θ z Θ (1) 0 D (0) (see appendix C.3.2).The cause of this change is the precise balance between Compton and DC scattering terms, as explained there.It is expected that the coefficient can change depending on which approximation for the DC emissivity is actually used.In paper III, we will see that this modification is not expected to be as severe, but additional work may be needed.Put together we then have the thermalisation terms ∂y (1)  ∂t therm ≈ τ θ z M K y (1) + D (1) As explained in appendix A.2.1, we neglected kinematic corrections to the perturbed thermalisation terms, as these contribute at order β τ θ z and are therefore smaller than the boosting terms ∝ β τ in the Thomson contributions.We do not expect this to modify the main conclusions although some details might differ during the µ-era.Alternatively, one could perform the computation in the baryon frame to avoid this complication of the calculation [see 36, for discussion], however, this is beyond the scope of this work.

Final evolution equation
We are now in the position to write the full evolution equation in terms of the spectral basis y (0) and y (1) .Together with eq.(2.9), (2.10) and eq.(2.25) we find ∂y (1)  ∂η + γ • ∇y (1) = −b (0) 0 ∂Φ (1)  ∂η + γ • ∇Ψ (1) + τ y 2 − y (1) + β (1) χ b where we made it explicit that the zeroth order solution is only a monopole and converted to conformal time, η.The first line accounts for all terms in the Thomson limit and the effect of anisotropic heating.The last row describes the spectral evolution, with M K y (1) 0 + D (1) 0 determining the main terms and the other being related to perturbed thermalisation effects.Before we convert this into a photon distortion parameter hierarchy in Fourier space, let us briefly discuss the expected effect of distortions on the other Einstein-Boltzmann equations.

Effect of distortion on the other perturbation equations
While we have completed our reformulation of the photon evolution equation in the presence of spectral distortions, we still have to consider how the other perturbation equations might be modified.As an example, let us update the momentum exchange equation with the baryons.For this we have to compute the integral 1 2 x 3 χC (1) [n] dx dχ.We shall neglect any small corrections from photon emission and absorption terms.Since the Compton terms all only act on the monopole,5 these also do not contribute and one is only left with the Thomson scattering terms.The change of the baryon momentum by scattering with photons at first order in perturbations then is [33] ∂β (1)  ∂η T ≈ −τ ρ z ρ b 1 2 2 − n (1) + β (1) χ Ôx n with β(1) ≡ iβ (1) and f (η, x, r) = i 2 P (χ)f (η, x, χ, r) dχ = i f 0 (η, x, r)/ 4π (2 + 1) based on the Legendre polynomials, P (χ).Also, ρ z and ρ b are the background level quantities (without any distortion effects included).

JCAP11(2023)027
The presence of distortions modifies the momentum exchange by a small amount.For consistency, these terms should be included in the hierarchies presented below.Here, ρ z (1 + δ (0) γ,0 ) is the total energy density of the average CMB with δ Similarly, the momentum carried by the photon dipole is modified to ρ z ( Θ(1) γ,1 ) with δ 1 .In the tight coupling limit, one then has 1 ∆n (0) ). (2.28) However, in practice the terms δ γ,0 and δ γ,1 are small corrections to the evolution equations of the standard variables unless the average distortion amplitude becomes close to unity, which is ruled out by COBE/FIRAS.In a similar way, the equations for the potentials and neutrinos as well as temperature polarisation states see negligible effects from the presence of distortions unless a full second order treatment is attempted.For our computation below we shall leave all the first order equations unchanged and only add the new equations describing the distortion anisotropies.Unless anisotropic heating effects are explicitly considered, this means that distortion perturbations are directly driven by the standard perturbations known from the evolution of the CMB temperature with no direct distortion anisotropy sources from the distortions themselves, the latter being O(∆n (0) ) 2 .

Expression in Fourier space
To obtain the SD transfer functions we now carry out the final step of going to Fourier space.This allows us to obtain the photon transfer functions using a treatment that is similar to that of the standard perturbation equations [52,53].All the variables appearing above are functions of X(η, χ, r).Going to Fourier space and using the expansion, X(η, χ, r) = (2 + 1) (−i) X (η, r) P (χ) in terms of the Legendre polynomials, we then have [33,37] where k denotes the wavenumber of the mode.The background evolution equation in eq.(2.26) needs no further modification.Carrying out the Fourier and Legendre transformations of all remaining terms, from eq. (2.26) we then find the final photon hierarchy Aside from the new terms in the monopole equations and the generalization of the potential and Doppler sources in terms of the distortion parameters this system is identical to the standard brightness equations [31,33].Setting the terms τ θ z and Q to zero and using b (0) 0 = (1, 0, . . ., 0) these are identically recovered. 6e note, however, that here we neglected corrections from polarisation terms, which affect the quadrupole equation [37,54,55].For the standard temperature terms we can include these as usual, but no attempt to correct polarisation effects for the distortion evolution are made here.This means that the damping scales of the temperature and distortion anisotropies will differ slightly.For the distortions, the Thomson scattering terms will be identical [56]; however, additional effects in the thermalisation terms will have to be studied more carefully, a task that is left to future work but should be possible, e.g., by starting from [36,39].

Line of sight integration and power spectra
To close our discussion of the theoretical aspects, we briefly present the line of sight approach for obtaining the final signal power spectra.In principle the solutions to the photon hierarchy in eq.(2.30) are enough to compute the power spectra and cross-power spectra of two observables X and Y : once the corresponding transfer functions, in this context indicated by the hat, are obtained.However, this brute force approach becomes numerically challenging, as known from the standard CMB anisotropies [57].To simplify matters, we start by Fourier transforming eq.(2.26b), which yields ∂y (1)  ∂η +ikχ y (1) +τ y (1) = S LOS , (2.32a) where now all variables are functions f (1) = f (1) (η, χ, k).Realizing that the only real differences with respect to the temperature only case are that we now are dealing with a solution vector and a vector for the sources, following the standard steps (see appendix D)
Concretely we give j (1,0) (x) = j (x) = ∂ x j (x) and j (2,0) (x) = 1 2 [3j (x) + j (x)] as in [37].In comparison to the standard line of sight approach, we have the extra source terms coming firstly from the spectral mixing by Compton scattering and perturbed thermalisation effects (∝ θ z ) and secondly from anisotropic heating (∝ Q (1) ).In addition, we have a more general Doppler and potential source vector (∝ b (0) 0 ), which in general is time-dependent.

Reduction to the experimental basis
As explained in section 3 of paper I, with an experiment in mind we can reduce the number of spectral parameters by going from the computation basis to the residual distortion representation.From the results for the transfer functions, this can be achieved by applying a matrix L to the solution vector y (1) .We can then plot the transfer functions for a more limited number of parameters, e.g., Θ, y, r 1 and µ, without loosing much information.In a similar manner we can simplify the computation of power spectra to those of the parameters in the experimental basis.Here is it best to directly apply the rotation to the source vector in the line-of-sight solution, eq.(2.33).This reduces the computational burden, since the power spectra need to be computed for a reduced number of variables.We comment that directly applying the projection on the limited observation basis before even computing the transfer function is not an optimal choice as in this case energy conservation is not guaranteed to the same level of precision as with the computation basis.Even if this would further reduce the computational burden, we thus do not recommend such an approach.

Basic expectations for the evolution of distortion anisotropies
While it is difficult to anticipate the detailed behavior of the distortion transfer functions by looking at the system in eq.(2.30), we can already understand some of the most important features.Firstly, without average distortions, non-standard evolution of the average temperature (see section 2.3.2) or anisotropic heating no distortion anisotropies are generated JCAP11(2023)027 at first order even in the perturbed universe [24].This limit is evident without further explanation, as in this case the system simply becomes identical to the standard temperature anisotropies at first order in perturbation theory.
Secondly, photon emission and Comptonisation effects ∝ θ z 4.6 × 10 −10 (1 + z) can only be important at z 10 4 , like for the average distortion evolution.At later times, the spectral evolution is mostly frozen and only direct distortion anisotropy sources from boosting and potentials (∝ Ôx ∆n (0) ) or anisotropic heating (∝ Q(1) Y ) can create additional distortion anisotropies at significant levels.We note that if isotropic heating is present (e.g., by particle decay) it is naturally expected that anisotropic heating occurs unless the heating mechanism is detached from perturbations in the standard cosmic fluid (e.g., independent of the dark matter density or local clock speed).In paper III, we will discuss how distortion anisotropies from decaying or annihilating particles will manifest, and then provide simple constraints on these cases.
Thirdly, the effect of distortions (and thermalisation effects) on the evolution of the standard perturbation variables is negligible unless non-standard sources of temperature perturbations are considered.This means that in many scenarios with distortion anisotropies, the transfer functions of the distortion variables will behave like a driven oscillator, following closely the corresponding temperature variables.This approximation is possible unless we look at the generation and thermalisation of large distortions during the early phases (z 5 × 10 6 ), where momentarily one could still imagine average distortions µ (0) 0 10 −2 [59,60] and hence efficient conversion into temperature perturbations of noticeable level.However, we leave a more detailed discussion of this regime to future work.
Moving away from the main time-dependent aspects, let us consider the scale-dependence of the distortion evolution.In the absence of anisotropic heating, distortion anisotropies are only generated once average distortions are present.This means that there will be a difference in the transfer functions depending on when during the cosmic evolution of the mode the injection occurs.Crucial moments are horizon crossing, diffusion damping and free-streaming, all known from the standard CMB temperature anisotropies [e.g., 33,61].
If we pick a mode of wavenumber k, we can distinguish three main regimes: i) the average distortion is generated after the mode has fully damped away by Silkdamping (k k D , where k D 4.0 × 10 −6 (1 + z) 3/2 Mpc −1 is the Silk-damping scale [54,62]).In this case, no distortion anisotropies are formed because fluctuations in the plasma are no longer present at the corresponding scale.
ii) the average distortion is generated when the mode is well within the horizon but at k k D .In this regime, no sources of distortions from potential variations arise, and distortion anisotropies are generated solely by Doppler terms and perturbed thermalisation effects while at z 10 4 .
iii) the average distortion is generated while the mode is still super-horizon.In this case, both Doppler and potential perturbations affect the distortion fluctuations.We will furthermore be sensitive to perturbed thermalisation effects when considering the evolution at z 10 4 .

JCAP11(2023)027
We will quantify these general expectations in paper III.We also note that once anisotropic heating is included, the statistical properties of the source of heat could be highly non-Gaussian in which case a more general Boltzmann system may have to be solved.In addition, direct sourcing of distortions in the regime k k D can still lead to significant distortions at very small scales, e.g., by phase transitions or non-linear baryonic physics.

Discussion and conclusions
This work provides the main formulation of a new photon Boltzmann hierarchy, eq.(2.26) and eq.(2.30), that allows us to compute the evolution of distortion anisotropies at first order in perturbation theory (section 2).Distortion anisotropy sources from Doppler and potential terms as well as anisotropic heating and perturbed thermalisation are accounted for.We also account for the spectral evolution of the distortion based on an approximate ODE treatment for the thermalisation Green's function, which captures most aspects of the full calculation using an extended spectral basis to describe the residual distortion evolution (see paper I).We furthermore demonstrate that the line-of-sight approach can be generalized to simplify the computation of the CMB signal power spectra (section 2.8).We briefly explain how the computation can be simplified by converting to a spectral basis that is optimised for the experiment using a basis rotation as explained in paper I (section 2.8.1) Overall, this paper is the second step in a series of works discussing the production and evolution of SD anisotropies generated by various physical mechanisms and how these might be constrained with future CMB spectrometers and imagers, enabling more realististic SD anisotropy forecasts over a wide range of physics.The new formulation furthermore is a first step towards a more general and precise treatment of CMB temperature anisotropies at second order in perturbation theory.In the standard approach [35,38,39], an 'instantaneous thermalisation' approximation is applied which ensures full energy conservation for the photons but without allowing a rigorous separation of distortion and temperature terms.With our new ODE representation of the distortion Green's function, one should be able to overcome this problem, even if additional generalisations will be needed.
We highlight that Doppler boosting and potential driving distortion source terms were omitted in most previous discussions of primordial SD anisotropies, although the importance of these terms was recognized earlier [24].These terms are indeed small with respect to the standard temperature perturbation (hence not affecting the their evolution notably); however, for the SD anisotropies they provide the leading order source terms once average distortions are present.Since SDs can be spectrally distinguished, these terms remain relevant, leading to the independent distortion parameter hierarchies we presented here.The resultant distortion anisotropies are expected to be significant as long as the average distortion is at the level that COBE/FIRAS allows.A more quantitative discussion will be given in paper III, and subsequent works.Importantly, in paper III we will demonstrate that this effect allows one to derive limits on the average heating rate using existing and future CMB anisotropy measurements.
Finally, the generalised Boltzmann system can still be improved.We did not consider the effect of polarised distortions nor thermalisation/Compton scattering in the anisotropies, which could further modify the result.We furthermore neglected kinematic corrections as

JCAP11(2023)027
well as small non-linear Comptonisation terms and Comptonisation effects in the anisotropies.A brute force treatment of the problem for single modes, even if computationally demanding, might reveal additional augmentations to the problem that may require extra attention.It would also be extremely important to formulate the problem in alternative gauges/frames, to check the consistency of the equations and also explore possibilities for simplifying the computations, finding analytic approximations, and clarifying the physical interpretation of the signals.We leave all these improvements of our method to the future; however, our main results should not be affected.We therefore conclude that one of the main steps towards quasi-exact computations of anisotropic SDs in the perturbed Universe is taken.

A.1 Thomson terms
Starting from the standard physics of CMB temperature anisotropies, the Thomson terms and first order Doppler boosts carry over trivially, leading to [e.g., 38] τ −1 C (1) [n] T = n 2 − n (1) + β (1) χ Ôx n (0) (A.1) without further ado.Here, we introduced the direction cosine χ = γ • β between the velocity and photon direction.We also defined n (t, x, r, γ) = m n m (t, x, r)Y m (γ) using the spherical harmonic coefficients of the photon occupation number, n m (t, x, r).It is important that, in contrast to the usual treatment, we now include the average distortion in Ôx n (0) = G + Ôx ∆n (0) , as this leads to distortion anisotropies, which would not arise otherwise.This is the leading order source term at the level ∝ τ β (1)  Ôx ∆n (0) 0 , which becomes noticeable once the mode enters the horizon.Without average distortions, this term vanished and consequently only temperature fluctuations are sourced.
In eq.(A.1), we neglected small Klein-Nishina corrections ∝ θ e = k B T e /m e c 2 , as these do not change the spectral shape of the photon field but merely modify the Thomson scattering rates of the dipole, quadrupole and octupole [24,64].Some of these corrections have been considered in [22].

A.2 Kompaneets terms
Let us next consider the Comptonisation terms from the scattering of electrons and photon with non-zero energy exchange, which leads to spectral evolution.We shall include these effects only for the local monopole as Thomson terms are not suppressed for > 0 and hence dominate the evolution there.The standard Kompaneets equation is [34,65] where we used Dx From the first group of terms we see that the local difference between the electron and photon temperature is the main source of distortions, leading to injection of Y with a source correction from Dx ∆n 0 .We have not yet included the effect of transforming to the local inertial frame, which causes an overall factor of (1 + Ψ) that usually only becomes important at second order in perturbation theory [35,63].However, here we are dealing with a non-vanishing zeroth order term, which then required inclusion of this factor already at first order in perturbation once non-zero distortions are present.
Adding this factor and collecting terms, we then have the Compton collision terms at zeroth and first order in perturbations,
At first order, the leading source term is again due to differences in the local electron temperature with respect to T z .However, since Compton equilibrium is reached quickly, ∆θ (1) e will be comparable to that of the local photon temperature perturbation with corrections from the local distortion and heating rate (see appendix C.2).This source term therefore mainly captures the modulation of the local thermal equilibrium by the variation of the ambient blackbody temperature.
The term ∆θ e Dx ∆n 0 .This is dominated by the blackbody part, ∝ D * x A ∆n (1) 0 , within the Kompaneets operator, so that we neglect this correction.It is clear that this term in principle is of similar order as terms that we indeed keep; however, since it merely modifies the exact timing of leading order scattering terms, we believe that its omission does not change the main conclusions significantly.A full numerical solution will be required to check the error this simplification introduces.For this, we can in principle start by mapping this term back onto the spectral basis just like for the Kompaneets operator, Kx .This will yield a representation in matrix form, creating a rotation of the spectral vector like for Kx ; however, we leave an exploration of this effect to the future.
The term ∝ τ (1) / τ are also due to scattering time-scale modulations but this time from perturbations in the electron density, which we include as τ (1) / τ = N b .We shall neglect perturbed recombination effects here [e.g., 35,50,51].These should never become important in the regimes we are interested in, when θ z is not too small already.The terms ∝ Ψ (1) account for the transformation to the local JCAP11(2023)027 inertial frame [35], which were omitted in [24].These terms lead to perturbed scattering effects once an average distortion is present.
To close the discussion of all the effects from Compton scattering, we mention that additional velocity corrections appear that are related to dipole scattering and kinematic effects.However, to simplify the problem we neglect these terms.From [24], it is clear that these terms can only be relevant at high redshifts, where significant energy exchange occurs, as we briefly discuss now.

A.2.1 Kinematic corrections to the Kompaneets term
As already mentioned in the main text, we do not include any energy exchange effects mediated by Compton scattering on the anisotropies in the spectrum.This means that in the rest frame of the moving thermal electron distribution one has where the primes on variables denote that they are evaluated in the moving frame.After a boost, the average occupation number in the moving frame is Here we neglected possible corrections O(βn 1 ) from aberration effects on the restframe dipole spectrum to the monopole in the moving frame.Using the invariance of Dx ≡ Dx and D * x ≈ (1 − βχ) D * x together with the transformation of the scattering y-parameter we then have ∂n(x, γ) where in the last step we used x ≈ x(1 − βχ) and n 0 (x ) ≈ n 0 (x) + βχ T (x) with T (x) = Ôx n 0 (x).This result is consistent with eq.(C37) in [24] after omitting dipole scattering effects.Assuming small departures from the average blackbody spectrum, ∆n 0 1, we can insert n 0 = n bb + ∆n 0 and then linearize in ∆n 0 , which with T = G + Ôx ∆n 0 yields: The first two terms are just the standard Kompaneets terms, while the last group of terms accounts for first order velocity corrections.In computations, the terms in braces would be -23 -

JCAP11(2023)027
evaluated at background level and are of order βχ θ z per Thomson event.However, given that from Thomson terms we have a source βχ T (x), these temperature correction terms are merely a small modification to the main source term (suppressed by a factor of θ z 1).It is therefore well-justified to neglect these contribution unless high precision is needed, as already mentioned above.A more rigorous study of these effects is deemed important for obtaining a self-consistent (gauge-independent) formulation of the problem, in particular once scattering corrections from the anisotropies [24] are included.

A.3 Photon number changing processes
Aside from scattering, which conserves photon number, we also have to treat the conversion of distortions into pure temperature shifts.This thermalisation process is mediated by the combined action of Compton scattering and photon emission processes (DC and BR).Modelling the exact evolution of the photon field when DC and BR are included is difficult; however, it is known that the evolution of the high frequency spectrum (ν 1 GHz) is not affected by these processes once the µ-era ends [15,18,19].We can therefore obtain an approximate description that simply leads to a redistribution of energy between the µ parameter and local photon temperature.The net photon emission and absorption term takes the explicit form [18,19,28] in the local inertial frame.We introduced the photon emissivity, Λ(x, θ e , θ γ ), which for DC scales as Λ(x, θ e , θ γ ) ∝ θ 2 γ being driven by the high-frequency blackbody photons [16,27,66].For convenience, we shall use the shorthand notation Λ(x, θ z ) ≡ Λ(x, θ z , θ z ) below.We neglect kinematic correction to the emission process and also do not consider the energy exchange corrections from this term.These are expected to be negligible and also require a more careful study of the emission process.Some first steps have been outline in [24].
Photon production processes are in equilibrium if the CMB occupation number is given by a blackbody at the electron temperature.Thus, defining the distortion with respect to n bb (x θ z /θ e ) would simplify several aspects of the computation; however, the invariance of the spectrum under redshifting is no longer guaranteed, such that we will not use this alternative description.
Perturbing the equation and neglecting terms that are higher order in the average distortions, we can find the emission term at zeroth and first order in perturbations as

JCAP11(2023)027
where we used (1 − e −x )/x = n bb /G.For the first order equation, we can see that the temperature derivatives of the emission coefficient modify the thermalisation effect.In the DC era, one has ∂ ln Λ/∂ ln θ e ≈ 0 and ∂ ln Λ/∂ ln θ γ ≈ 2, where the latter signifies that most of the DC emission is driven by the blackbody part of the CMB spectrum [27,66,67].
Corrections to the picture from BR will be neglected here, but can in principle be added using BRpack [26].However, in the BR era, photon production is already nearly frozen, such that our approximations should not make a significant difference to the final distortion evolution.

C.1 Compton electron temperature and Compton energy exchange
The energy exchange between electrons (and matter) and photons is controlled by the Compton scattering collision term.Smaller corrections due to BR and DC emission processes can be neglected.From eq. (A.3), we can directly write the zeroth order energy exchange term as8 C ≡ Θ (0) eq ≡ η ∆n (0) where E n bb is the energy density of a blackbody at a temperature T z and Θ e = ∆T e /T z denotes the relative temperature difference of T e and T z .For convenience we JCAP11(2023)027 also introduced the background Comptonisation rate, κ, which will be used frequently below.We furthermore used the Compton integral, η f , and the energy integral, E f = x 3 f (x) dx as defined in paper I.
With eq.(A.4), but setting Ψ (1) = 0 to evaluate in the local inertia frame, we can also directly write the first order Compton energy exchange rate as (1)  e − x 3  Kx ∆n e Dx ∆n C − Θ (1)   e + δ (1) b + 4Θ (1) 0 0 Θ (0) eq ≈ Θ (1)  eq − Θ (1) 0 , which accounts for a small correction to the energy density of the zeroth order photon field.We will see that this term ensures energy conservation with respect to the scattering correction, ∆θ e Dx ∆n (0) 0 in eq.(A.4).We comment that the Compton equilibrium temperature corrections, Θ C and Θ (1) C , could also have been obtained directly using the well-known expression [68,69] both at zeroth and first order in perturbation when non-linear distortion terms (∝ ∆n 2 0 ) are neglected.The direct path using the Compton collision term highlights the consistency of the result.

C.2 Electron temperature equation and effective photon heating rate
In this section, we consider solutions to the local electron temperature including perturbations in the medium.Due to the presence of Compton scattering, the electron temperature is always pushed extremely close to the Compton equilibrium temperature, with corrections from the heating terms that appear in the electron temperature equation.This leads to an effective heating term in the photon equation that is obtained here for conditions in the pre-recombination era.
To obtain the correct terms, we adapt the discussion of [35].Denoting time-derivatives in this context with 'dot', during the pre-recombination era and consistent up to first order in perturbations one has9 for the electron temperature.Here, C V = 3 2 k B N tot is the non-relativistic heat capacity of the baryons, with N tot = N H + N He + N e being the total number density of baryonic particles.

JCAP11(2023)027
Armed with these ingredients we can now write down the evolution equation for the electron temperature at background and perturbed level.This yields e +2HT (1)  e +2 β (1)  3a + Φ(1) T (0) e = Ψ (1) Assuming the plasma to be fully ionised one finds C (1) b .The background level equation is well-known in connection with the standard recombination and thermalisation problems [e.g., 18,70].In the pre-recombination era, T e will follow a sequence of quasi-stationary states due to rapid Compton interactions, pushing T where in the second step we neglected the correction from the adiabatic cooling effect, which leads to a tiny average distortion µ cool −3 × 10 −9 [28,71,72].We could have directly obtained the solution by setting Λ (0) + Q(0) c ≈ 0, a fact that we will use when computing the perturbed temperature solution.Starting with eq.(C.5) for the electron temperature, we then obtain for the Kompaneets terms.The combination Θ eq Y + Kx ∆n (0) 0 can be replaced with the Kompaneets matrix description → M K y (0) .Here, Θ 0 once ∆n 2 0 terms are neglected.From this we can also identify the relative effective heating rate which then leads to the formulation for the background distortion evolution given in the main text.We note that when exact Compton equilibrium is no longer reached, we can directly solve the zeroth order temperature equation given above to compute the main distortion source term in eq.(C.6).In this regime, we can neglect the terms τ θ z and simply find direct sources of y-distortions without any spectral evolution.This only is expected to become important at z 200, which is a regime we do not consider at this point, such that the quasi-stationary solution is valid.
To obtain the perturbed temperature solution, we take the anticipated short-cut assuming quasi-stationary conditions (or more accurately, C V H/κ 1).This then yields 0 ≈ Ψ (1)

JCAP11(2023)027
This expression neglects the perturbed corrections to the adiabatic cooling process, which are sub-dominant for our purposes, but should lead to a minimal distortion anisotropy in ΛCDM.In our limit, Λ C − Θ (1)   e − δ b + 4Θ (1) 0 By solving for Θ e and neglecting higher order distortion terms, we then finally find Θ (1)  e ≈ Θ where the term ∝ 4Θ (1) 0 accounts for variations of the local photon heat capacity.Putting everything together, from eq. (A.4) we find ∂∆n (1) 0 ∂t K ≈ τ θ z Θ (1)  e Y + Kx ∆n 0 + Θ (1)  e Dx ∆n 0 Dx ∆n 1)  eq Y + Kx ∆n (1) 0 for the Kompaneets terms at first order.We can identify Θ eq Y + Kx ∆n (1) 0 = M K y (1) like for the zeroth order.We note that in the second line we used the equilibrium temperature differences to isolate the effect of external heating terms.We also set Θ

JCAP11(2023)027
The only terms in eq.(C.11) leading to addition of energy are those ∝ Q(1) c /ρ z + Ψ (1) Q(0) c /ρ z .We thus identify them as which is used in the main text to formulate the first order heating rate.We highlight that the terms Q(0) c simply follow from the collision terms in the local inertial frame.One would have naturally guessed this perturbed heating term when thinking about Q ≈ (1 + Ψ) Qc /ρ z .

C.2.1 Adding extra Doppler terms to the system
Since the terms in the last line of eq.(C.11) do not add any extra energy (or photons) to the system, we can in principle neglect them in our computation, without changing the consistency of the system.However, it is fairly easy to include them once Dx ∆n (0) 0 is evaluated.Since they are expected to change the exact distortion evolution across the residual era, this could provide extra time-dependent information.Inserting our distortion representation, with Dx = −3 Ôx + Ô2 x we have where we separated the distortion only terms (i.e., those without G).We can in principle again use our projection method to obtain descriptions of Dx Y N −1 , Dx Y N and Dx M in our computation basis, but it is much easier to just apply the boost matrix, M B , defined in eq.(2.9) for this purpose.It turns out that we then have Dx ∆n where I is the identity matrix.One can precompute the Doppler matrix M D = (M B − 3I)M B for numerical applications to ease the computations.

C.2.2 Compton equilibrium spectrum
An important property of the Compton collision term operator is the associated equilibrium spectrum, which leads to a stationary spectrum under repeated scatterings.At zeroth order in perturbations, we trivially have ∆n 0 , as can be verified by using Kx n 0 , which leads to cancellation of the Kompaneets term at zeroth order.
How does this work for the first order Comptonisation terms?If we assume no average distortion or temperature shift and no heating, then from eq. (C.11) we naturally have ∂∆n (1) 0 ∂t K ≈ τ θ z Θ (1)  eq Y + Kx ∆n

JCAP11(2023)027
determined (at leading order) by the implicit equation Λ(x c , θ z ) ≈ θ z x 2 c , which assumes that the emission coefficient, Λ(x c , θ z ), is a slowly varying function of x.In the DC-era, we have x c ≈ 8.6 × 10 −4 (1 + z)/2 × 10 6 , assuming the standard CMB temperature for T z .The idea is now to write the Ansatz ∆n 0 m(x) with frequency-dependent correction functions g(x) and m(x).At x x c one has g(x) m(x) 1.It is therefore useful to consider g(x) and m(x) as functions of ξ = x c /x instead of x.Starting from eq. (C.17), multiplying it by −x 2 c /[θ z ξ 4 ] (this is the leading order dependence) and replacing x → x c /ξ, with the Ansatz ∆n up to first order in x c 1. Neglecting the terms ∝ x c , we have the physical solutions m(ξ) = exp(−ξ) and g(ξ) = 1, which together gives for the modified temperature and µ-distortion spectrum.The goal is now use this to obtain an expression for the photon production rate at lowest order in x c .Integrating the collision term in terms of photon number, 10 we find 0 .Indeed, the reduction of the chemical potential amplitude is given by ∂ t µ (0) 0 em/abs ≈ −γ N τ θ z x c µ (0) 0 [19], which we use in main text.We highlight again that without the low-frequency modification to the spectrum by the factor e −xc/x one would not have obtained a finite photon production rate.At low frequencies, the distortion vanishes and the spectrum is in full equilibrium with the electrons, stopping any emission and absorption.
We note that, while for the temperature term ∝ Θ (0) 0 an exact cancellation occurred, in the final steps we had to drop the term η M µ (0) 0 , since it still diverges logarithmically towards x → 0. This is because we did not include higher orders in x c .We can in principle use eq.(C.18) to obtain correction to the distortion spectra at first order in x c .However, the JCAP11(2023)027 procedure becomes quite involved requiring additional normalisation conditions and also the including of time-dependent corrections [see 19].This is beyond the scope of this work and left to a future publication.
We also note that one can in principle evaluate the exact integral x 2 c M (x)m(x) dx numerically.This yields the correction x c → x c (1 − 1.65x 0.88 c ), which indicates that photon production is a little slower than in the soft photon limit.However, in terms of perturbations in x c this is not a consistent treatment and other corrections are expected to have a similar level.We therefore do not recommend adding these corrections until a more complete treatment of the problem is attempted.

C.3.2 First order treatment of photon production
We can now repeat the same computation but at first order in perturbations.The emission terms at first order in perturbations is given in eq.(A.10b).Remembering that ∆n (0) 0 ≈ −µ (0) 0 m(x)/x 2 gave the leading order term at zeroth order and checking the coefficients of all functions, keeping only leading order terms, we then find

D *
x ∆n (0) 0 G(x), which naturally cancels some of the logarithmic corrections at low frequencies.However, in the overall evolution at high frequencies it should not be as important.
From the first order computation it is clear that the term ∝ δ (1) b + Ψ (1) will be canceled by the corresponding zeroth order emission term.Inserting ∆n 0 .This suggests that the naive estimate does not capture the full effect of changing the DC emissivity with temperature.
While this is somewhat satisfying, we note that the required effect stems from a higher order correction caused by the absorption term ∝ −∆n 0 (e x θz/θe − 1) in eq.(A.9), which does not contribute at zeroth order, and hence can also not be captured by perturbing the zeroth order equation.In reality, additional corrections can be expected since the frequency dependence of the DC emissivity, Λ, may also modify matters.This was highlighted before as part of the analytic computations for the distortion visibility function [19,71].Overall this shows that the precise scaling may depend on which approximation for the photon emissivity is indeed used.However, a more detailed discussion is beyond the scope of this work.

Figure 1 .
Figure 1.Representation of Ôx M (x) using varying elements in the basis.Even the lowest order representation in terms of µ and y is highly accurate.

0
leads to a correction to the Comptonisation time-scale 7 of ∆n (1) 0 , with the dominant term given by θ z Kx ∆n the last step of eq.(A.4).We can also rewrite the term Dx ∆n described in the main text to simplify the correction ∆θ

B. 2 2 Θ2
Changing the average temperature at second order in ΘIn section 2.3.2,we showed that at first order in Θ, a photon source term ∝ G(x) leaves the spectrum invariant [confirm eq.(2.11)].Here we now ask what photon source do we need to obtain a blackbody up to second order in Θ 1. Writing the Taylor series for a blackbody as n ≈ n bb (x) + [ Θ + Θ2 ] G(x) + 1 Y (x) and taking the time derivative (denoted with dot at fixed x) we haveṅ ≈ G(x)[ Θ + 2 Θ Θ] + Θ Θ Y (x) ≡ S G [G(x) + f Y Y (x)].(B.5)If we demand that the average blackbody temperature should change like Θ = S T , then we can trivially writeS G = (1 + 2 Θ) S T .Similarly, we then have Θ Θ ≡ S G f Y , which implies f Y = Θ/[1 + 2 Θ].The overall evolution equation is then given by ṅ≈ S T [(1 + 2 Θ)G(x) + Θ Y (x)], which leaves a blackbody spectrum with changing temperature unchanged at second order in Θ.Note that the 'blackbody source' is no longer independent of the solution for the blackbody temperature and requires a y-type contribution to shift the maximum of n towards higher frequencies.At second order, one thus cannot create a blackbody with a temperature-independent source.However, we have an alternative way of describing the same situation.Defining the effective temperature and y-parameters as Θeff = Θ[1 + Θ] and ȳeff = Θ2 /2 ≈ Θ2eff /2, we can similarly determine the required source functions for G and Y.For the source of G, we then simply need S G ≡ Θeff = (1 + 2 Θ) S T = 1 + 4 Θeff S T ≈ (1 + 2 Θeff ) S T .Similarly, we have S Y = ẏeff = ΘS T = ( 1 + 4 Θeff − 1)S T /2 ≈ Θeff S T .Since up to second order in Θ we have Θ2 /2 ≈ Θ2 eff /2, this also directly follows from S Y ≈ Θeff Θeff ≈ Θeff S T .We therefore find the same source term as above, but with Θ → Θeff .The solution for Θeff is then Θeff = [exp(2 S t dt) − 1]/2 = [exp 2 Θ − 1]/2 ≈ Θ[1 + Θ], which as expected shows the equivalence of the two approaches.

3
that the other corrections to Θ(1) e are higher order in the distortion.It is important to highlight that the terms in the last line of eq.(C.11) do not add any extra energy to the system, but merely redistribute the energy spectrally, leading to a correction to the perturbed scattering effect.For the terms proportional to ∝ Y 1 − Y this is trivially seen when computing x 3 Y dx = x 3 Y 1 dx = 4E n bb .For the other term we havex Dx .With x 3 Y dx = 4E n bb this identically chancels the remaining term in the last line and also highlights how the term ∝ δ (0) γ,0 is crucial for conserving energy.-29-

x 2 e
−xc/x dx = 3 4 γ ρ γ N θ z x c µ term.This remarkable result is the key to understanding the rate of conversion from µ ∂n