Spectro-spatial evolution of the CMB. Part I. Discretisation of the thermalisation Green's function

Spectral distortions of the cosmic microwave background (CMB) have been recognized as an important future probe of the early Universe. Existing theoretical studies primarily focused on describing the evolution and creation of average distortions, ignoring spatial perturbations in the plasma. One of the main reasons for this choice is that a treatment of the spectro-spatial evolution of the photon field deep into the primordial Universe requires solving a radiative transfer problem for the distortion signals, which in full detail is computationally challenging. Here we provide the first crucial step towards tackling this problem by formulating a new spectral discretisation of the underlying average thermalisation Green's function. Our approach allows us to convert the high-dimensional partial differential equation system (≃ 103–104 equations) into and set of ordinary differential equations of much lower dimension (≃ 10 equations). We demonstrate the precision of the approach and highlight how it may be further improved in the future. We also clarify the link of the observable spectral distortion parameters (e.g., μ and y) to the computational spectral basis that we use in our frequency discretisation. This reveals how several basis-dependent ambiguities can be interpreted in future CMB analysis. Even if not exact, the new Green's function discretisation can be used to formulate a generalised photon Boltzmann-hierarchy, which can then be solved with methods that are familiar from theoretical studies of the CMB temperature and polarisation anisotropies. We will carry this program out in a series of companion papers, thereby opening the path to full spectro-spatial exploration of the CMB with future CMB imagers and spectrometers.


Introduction
Spectral distortions (SDs) of the cosmic microwave background (CMB) have now been recognized as an important future probe of early-universe and particle physics.In particular the ability of CMB SDs to constrain the primordial power spectrum at small scales [1][2][3][4] provides important motivation to push the observational frontier with the next generation CMB experiments [5,6].However, much of the recent theoretical work [e.g., [7][8][9][10][11] and experimental spectrometer concept studies [12][13][14][15][16][17] focused primarily on the science of average distortion signals.It is known however that an average distortion signal can form anisotropies through interaction with other perturbed quantities [4].While average distortion physics is often concerned with earlier times and thus sourced by larger k-modes, the anisotropies formed from this average distortion would be prevalent on similar scales to usual CMB anisotropies [explored further in 18, 19].From the theoretical point of view it is clear that distortion anisotropies should be smaller and thus harder to detect, however, one of the main reasons for preferring average distortion science is that the computations of distorted CMB anisotropies of primordial origin are difficult and currently beyond the possibilities of existing Boltzmann codes.

JCAP11(2023)026
To illustrate this statement we highlight that numerically solving the full thermalisation problem for SDs created on average by energy release now takes of order 30 seconds on a standard laptop using CosmoTherm [7,20].While this is already highly optimised, it will be difficult to extend this method to SD anisotropies, where in analogy to the standard CMB temperature fluctuations [21,22] one would have to solve the thermalisation problem for multiple k-modes.For each k-mode, a multipole hierarchy would furthermore be required, overall boosting the computation by a factor 10 3 .In addition, one would have to consider how to convert the final (frequency-dependent) signal transfer functions into CMB observables, which further increases the complexity of the problem over the standard CMB anisotropy computation, likely yielding single computations that would take O(10 5 -10 6 ) seconds.While not necessarily prohibitively expensive with modern computational resources, this brute force approach would be overly-complicated for exploratory calculations and not scalable in parameter forecasts and searches for new physics.
How could one make the problem more tractable?The most common approach is to simplify the problem by considering limiting cases.In particular, scenarios in which the evolution of distortions and primordial perturbations as well as thermalisation physics can be mostly separated come to mind.This brings us to the well-known Sunyaev-Zeldovich (SZ) effect [23,24], which is created by anisotropic heating effects in the late Universe, sourcing y-type distortion anisotropies that peak at several arcminute angular scales.This signal is highly non-Gaussian and requires an understanding of the non-linear large-scale structure evolution, but then analytically translates the statistical properties of the dark matter distribution into the y-field [25][26][27][28][29].The SZ effect is therefore an important probe for cosmology and cluster physics [30,31].
Another example is the sourcing of y-distortion anisotropies by the mixing of blackbodies in the perturbed universe [1,2].This second order effect leads to a fluctuating y-distortion sky [32,33] in addition to an average distortion [4] when perturbations dissipate by freestreaming and Thomson scattering effects.For the fluctuating part, no spectral evolution has to be considered at the late stages (redshift z 10 4 ), just like for the SZ effect -a linear perturbation description of the problem is furthermore possible, yielding y-parameter transfer functions that are excited by first order temperature perturbations [33,34].If the amplitude of the small-scale curvature perturbations is modulated by large-scale modes this can furthermore lead to correlated µ × T and y × T fluctuations [35][36][37][38][39][40][41][42], which can be directly constrained using CMB imagers [see [43][44][45], for most recent forecasts and constraints].Note that at the largest angular scales, the corresponding transfer problem was simplified by neglecting details of the distortion evolution in the perturbed Universe [35,40,46].
There are, however, a number of aspect to the thermalisation problem that have not been captured by any of these calculations.As explained in [4], if an average distortion is present during the pre-recombination era, the standard density perturbations at first order will source distortion anisotropies.Assuming the average SD is ∆n (0) ν in terms of the photon occupation number, the SD anisotropies will have a spectrum that follows ∆n (1) ν ∝ −ν∂ ν ∆n (0) ν [4].Even without any spectral evolution, the standard Doppler terms and potential perturbations therefore source distortion anisotropies, which have not been evaluated for a more general average spectrum [see 33, for a treatment of occupation number temperature and y-distortion].

JCAP11(2023)026
Assuming that the average distortion saturates the limits imposed by COBE/FIRAS [47,48], one can expect distortion anisotropies at the level of 10 −8 -10 −7 of the average CMB.This can exceed the signals expected from the aforementioned non-Gaussian signals and can also be directly constrained with existing and future CMB imaging data.In addition, the thermalisation efficiency should vary from patch to patch in the perturbed Universe.The required terms in the photon Boltzmann equation were already discussed in [4]; however, only recently has the effect been estimated using a separate universe approach [49].In particular for modes that cross the horizon at or after the recombination process completes this effect should be noticeable in the transfer function solutions, but has not been computed using a full Boltzmann treatment.
To fully capitalise on the potential of spectral distortion anisotropy studies, as a first step we need to formulate a generalized photon Boltzmann equation that goes beyond the standard temperature and polarisation anisotropies.The biggest bottleneck is due to the discretisation of the spectral evolution, which currently is done with 10 3 -10 4 bins in frequency, as explained above.In this work, we obtain a new discretisation for the average frequency evolution that reduces the computational burden by a factor of 10 3 (section 2).This allows us to model the thermalisation from y → µ → T with a small number ( 10) of new spectral parameters, that can represent the exact calculation from CosmoTherm to high precision.In contrast to other approximations, the solution is no longer limited to the three standard spectral shapes but allows one to capture the dominant contributions from the residual distortion [e.g., 50].We also explain how the computational distortion parametrisation can be mapped back onto the leading residual distortion spectra, which present the main spectral shapes that may be testable in future applications (section 3).we close with a brief discussion of the next steps and remaining limitations in section 4.
This paper is the first in a series of works that study the effect of spectro-spatial evolution of the CMB.In paper II [18], we will formulate the generalised Boltzmann equation for distortions cause by heating processes, strongly drawing on the results of this paper.In paper III [19], we will present a detailed discussion of the distortion transfer functions and power spectra, highlighting the importance of various physical effects and providing Fisher forecasts.We also plan subsequent papers that discuss how the dissipation of acoustic modes in the presence of primordial non-Gaussianty causes spectral distortion anisotropies, and which constraints on various scenarios can be expected.Overall we hope this will provide further motivation to study SDs in the future.

Approximate ODE representation of the thermalisation Green's function
In this section, we establish a novel way of modeling the spectral evolution of the average photon field under repeated Compton scattering and thermal photon emission processes.In terms of perturbation theory, this is akin to focusing on the background quantities only, which leads back to the Green's function approach for the thermalisation problem, as will be developed here.Anisotropic distortion evolution from heating processes will be considered in papers II and III of the series.

Brief recap of the thermalisation Green's function
The efficiency of photon production and Comptonisation in the primordial plasma dictate various eras with characteristic SD shapes which are defined below.At sufficiently early times, in the temperature or T -era (2 × 10 6 z), thermalisation processes are very efficient and any excess energy is rapidly converted into a temperature shift, G(x).Here x = hν/kT z where T z = T 0 (1 + z) is the background reference temperature, which is chosen to match today's CMB temperature T 0 = 2.7255 K [51]. 1 The subsequent µ-era (5 × 10 4 z 2 × 10 6 ) is characterised by a lack of photon production, leading to a chemical potential distortion, M (x).Finally, the y-era (z 5 × 10 4 ) renders photon energy redistribution inefficient, leading to a distortion, Y (x), related to the well-known SZ effect, albeit in this case of primordial origin.
The different characteristic spectra introduced above have the forms [e.g., 10, 52] with 1923.In section 2.4 we will describe how these can be obtained by boosts of the average blackbody spectrum.Central properties of these spectra are summarised by their dimensionless photon number density N f = x 2 f (x) dx and energy density E f = x 3 f (x) dx.The corresponding integrals can be carried out analytically in terms of Riemann ζ-functions.Also using the blackbody occupation number, n bb (x) = 1/(e x − 1), we then find: ) The absence of overall photon number for y and µ type distortions is by construction (and easily achieved by subtracting G from alternative definitions).This convention has already been commonplace in the literature, but will become a fundamental simplifying fact in the novel treatment introduced below.
While the heuristic decomposition into three distinct eras introduced above conveys the correct physics to relatively high precision, it is much more convenient to have a robust framework in which the results can be expanded and built upon.In [53] it was shown that the thermalisation problem can be expressed as a Green's function problem in the limit of small energy injection: where s ∈{Θ ≡ ∆T/T z , µ, y} gives the signal amplitude of the corresponding SD (f ∈{G, Y, M }), J s is a dimensionless energy branching ratio, and , where dQ c /dz directly follows from the photon collision term.For clarity we note that the three era picture of the early Universe would correspond to simple top-hat functions for J s [i.e., see 'Method A' in 54].Other approximations for the energy branching ratios of varying accuracy exist [54], including the addition of intermediate spectral shapes known as residual distortions [50] and perturbative SD approximations for moderate scattering y-parameter [55].A byproduct of this work is the ability to generate accurate Green's functions in a generalized spectral basis to a precision comparable with full numerical treatments (see section 2.2).
While the spectral shapes in eq.(2.1) are physically motivated -each characteristic of a limiting case for each phase in the early Universe -they are insufficient to model the general evolution of the spectrum caused by heating processes.In the following sections we introduce a method for extending this set of spectral functions and explain how this new spectral basis eventually allows for full spectro-spatial solutions of primordial perturbations in the photon field.

Basic idea and lowest order solution of the thermalisation problem
As previously mentioned, in the µ-era all injected energy rapidly converts into the µdistortion [1,[56][57][58].The net µ-parameter is given by the evolution equation ∂µ ∂t ≈ γ ρ Q, where γ ρ ≡ α M ≈ 1.4007 and Q = dQ/dt.For a given Q, this equation can be solved with initial µ = 0.
Physically, the energy injection first leads to an increase in the distortion y-parameter by ẏ ≈ Here, 4y is the relative momentary energy density within the y-distortion part, θ z = k B T z /m e c2 is the dimensionless temperature, 2 and τ = dτ / dt = N e σ T c denotes the differential Thomson optical depth, all with the common choice of constants.The average energy exchange rate is ∆ν/ν 4θ z per scattering [59,60], which determines how quickly energy flows from y to µ.This identifies τ θ z as a fundamental timescale in the thermalisation problem which contrasts with the timescale of Thomson scattering τ -a fact which will become important for the generalised Boltzmann hierarchy (paper II).As we see in figure 1, the solution of this simple system roughly captures the transition between the µ and y-eras, yielding the y-distortion visibility J y ≈ e −4 yz , where we introduce the scattering y-parameter y z = τ θ z dt, and, from energy conservation, the µ-visibility J µ ≈ 1 − J y .
Following [61], the reduction of the chemical potential by the Bremsstrahlung (BR) and double Compton (DC) processes is approximately given by μ| em/abs ≈ −γ N x c τ θ z µ with γ N ≈ 0.7769.Here, x c is the critical frequency between DC/BR emission and Compton scattering, with approximations as a function of redshift given in section 3.3.1 of [11].Every absorption event then removes ∆ ργ /ρ γ ≈ α −1 M μ| em/abs of energy from the µ-distortion, which JCAP11(2023)026 with γ T = γ N 4 γρ ≈ 0.1387, a task that can be easily carried out numerically.Assuming that the photon production process only becomes important when y is already negligible, it is easy to show that J T ≈ 1−J bb with distortion visibility J bb ≈ e −(z/zµ) 2.5 and z µ = 1.98×10 6 [57,58].We then have As can be seen from figure 1, these simple approximations already capture the main dependence of the distortion visibility on the injection redshift.The question of the next section is now whether we can improve on this description to also include terms relating to the residual distortion.

Preliminaries
Neglecting photon production and heating terms, the relevant evolution equation for the distortion, ∆n = n − n bb , from the blackbody n bb = 1/(e x − 1) in the expanding Universe

JCAP11(2023)026
can be cast into the compact form [7,53] ∂∆n(x) where Θ e = ∆Te Tz is the relative electron temperature difference and Kx = Dx + D * x A is the Kompaneets operator, constructed from the diffusion and recoil operators, Dx = x −2 ∂ x x4 ∂ x and D * x = x −2 ∂ x x 4 , with A = 1 + 2n bb = (e x + 1)/(e x − 1).The time variable is the scattering y-parameter as defined above.The problem has been linearised in the distortion, an approximation that will be good unless very large distortions are encountered [11,20]. 3or the electron temperature correction, ∆T e = T e − T z , we assume that Compton equilibrium is reached at all times. 4In the absence of external heating, this means that x 3 ∂ yz ∆n dx ≈ 0, which implies Θ e ≈ Θ eq with [e.g., see 63] Kx ∆n dx and the y-weight factor w y = Y /G = xA(x) − 4 = x e x +1 e x −1 − 4. Since the integrals in eq.(2.8) will appear multiple times, for convenience we introduce The exact integrals that are encountered in our computations can all be given in terms of the Riemann ζ-functions.For the basic spectral shapes we have η G = 1, η Y ≈ 5.3996 and η M ≈ 0.4561, as well as G = 4, Y = 4 and M = 1/1.4007.For numerical applications we pre-compute all these integrals.

Spectral basis and approximate representation of the Kompaneets operator
The goal is to find an efficient spectral representation that captures the changes of the spectrum under repeated Compton scattering as described by eq.(2.7).The simplest decomposition considers the three main spectral types appearing in the thermalisation problem introduced in eq.(2.1).To build intuition, we discuss this case in some detail, but eventually find it is insufficient.The required refinements are presented right after.
It is instructive to understand the links of these basic spectra to that of the background blackbody spectrum.Both G and Y are generated by applications of the boost generator,5 Ôx = −x∂ x : (2.10)

JCAP11(2023)026
Making the Ansatz ∆n = Θ G(x)+y Y (x)+µM (x) and inserting back into eq.(2.7) we obtain where the prime indicates the derivative with respect to y z .The stationary solution of the Kompaneets equation, Kx G(x)/x = 0, defines the null-space.Hence, the number-conserving definition of the µ-distortion, which we are using here, will transform like the temperature shift.In fact both functions nicely map back onto Y (x).However, the function K Y (x) = Kx Y (x) has contributions that are not spanned by G(x), Y (x) and M (x).We can nevertheless enforce a representation of Here, we are mostly interested in intensity functions.
Knowing that K Y (x) does not carry photon number, 6 we understand that only Y (x) and M (x) can contribute in our current basis.Since with our finite basis, the representation will not be exact, we can demand that the contribution of M to follow from energy conservation to improve matters, as discussed below.This operator-function view of the mathematics allows us to make a choice of an inner product on the functions.We can choose to define the matrix elements of the operator X between the two function F (x) and J(x) as Although by no means rigorous, this inner product-like operation helps us to define the expansion coefficients heuristically.It is equivalent to taking the integrals over the two intensities x 3 F and x 3 X J, and will be helpful later in ensuring energy conservation of the entire system.
To decompose K Y = Kx Y in the most simple approach we remap back to the basis using the Ansatz | Kx Y ≈ a 0 |Y + a 1 |M and solve the system T and a = (a 0 , a 1 ) T .Carrying out the projection integrals and inverting the system we obtain K Y ≈ −8.8169 Y (x) + 40.409 M (x).However, since we used an incomplete basis, this approximation does not satisfy energy conservation.Carrying out the energy integrals, we find

JCAP11(2023)026
Since energy and photon number conservation are the most fundamental aspects of the thermalisation problem, this is not a solution we can work with.
To fix the problem, we replace the last equation in the system, eq.(2.14), with the energy conservation equation.This yields the augmented system (2.15a) We have now reformulated the problem once we also determine Θ e ≈ Θ eq .In vector notation, our Ansatz reads ∆n = B • y, where now we include G(x) in the basis, i.e., B = (G(x), Y (x), M (x)) T and y = (Θ, y, µ) T .By inserting this Ansatz for ∆n into eq.(2.8) for the Compton equilibrium temperature perturbation, and carrying out the energy exchange integrals one finds Inserting everything back into eq.( 2.11) and collecting terms, with eq. ( 2.12) we then obtain (2.17) We note that the terms in the Compton equilibrium temperature ∝ Θ and µ cancel identically due to the identities in eq.(2.12).We furthermore comment that eq.(2.17) can be also obtained by directly carrying out the projections onto the basis starting from eq. (2.7).We show this more formally in appendix B for the extended basis that is discussed in the next section.
Since the system in eq. ( 2.17) has to be fulfilled for any x and because the spectral basis is non-degenerate, by comparing coefficients, we obtain the ordinary differential equation (ODE) system with Y / M ≈ 5.6026.While these equations correctly represent the conservation of photon number [only G(x) carries photon number but Θ does not change] and also energy (the sum of the energies in µ and y does not change), they do not yield the correct overall evolution: for the y parameter, the solution is y(y z ) y(0) e 1.9403yz , while we saw in section 2.2 that it should be more close to y(y z ) ≈ y(0) e −4yz , which due to the sign ensures that energy correctly flows from y → µ.What has gone wrong?The approximate representation of K Y (x) ≈ −3.4593 Y (x) − 10.871 M (x) is insufficient, as could have been guessed.This can be appreciated in figure 2, where we compare the exact solution of K Y (x) with various approximations.In particular the high-frequency part of K Y (x) is not well-captured by this simplest approximation, a problem that we fix next.

Extension of the basis
To make progress, we need to extend the spectral basis, B. One of the natural selections is to use the boost operator Ôx to find the extensions. 7In principle other choices that could potentially even simplify the calculation either from the theoretical or the computational point of view can be imagined.Our choice is motivated by the fact that Ôx is one of the fundamental operators generating the Kompaneets operator, Kx .In fact, at x 1, the Kompaneets operator reduces to Klow x = 6 + Ôx ( Ôx − 5), which commutes with Ôx (i.e., [ Klow x , Ôx ] = 0) indicating a common basis.The boost operator also commutes with the diffusion operator, [ Dx , Ôx ] = 0, which further supports this choice also in more general cases (see appendix A).Finally, it appears in log-moment expansions of distortion spectra related to heating processes, which were shown to have useful properties in terms of gaugechoices [32,34,64]. 8For G(x), we have Ôx G(x) = Ô2 x n bb (x) = 3G(x) + Y (x), which directly maps back onto the old basis.For the boosts of M (x) and Y (x), new spectral shapes are generated.However, since K M = Kx M (x) = −η M Y (x) already maps back onto our basis, for now we only need to think about extensions based on the functions These functions can be readily computed using Mathematica or through the combinatoric sums given in appendix A, and are illustrated for a few cases in figure 3. The Y k are similar to those functions appearing in asymptotic expansions of the SZ effect [65][66][67][68] and can furthermore be found in perturbative expansions of the photon transfer problem [42,55,69,70].

JCAP11(2023)026
Note that Y 0 (x) ≡ Y (x).We also added the factor of (1/4) k to make each of the Y k more comparable in amplitude.This choice also ensures Y k ≡ 4.These functions all conserve photon number ( x 2 Y k (x) dx = 0) and hence provide a natural extension of the simple Y and M basis.As we will see in paper II, these also naturally appear once Doppler-driving in the perturbed Universe is included.

Generalization of the ODE system
In this section, we outline the basic approach for obtaining a generalized ODE system in the extended basis.By deciding about how many Y k we include in the Ansatz for ∆n, we have to determine the representations for each of the9 K Y k = Kx Y k within this basis.
To simplify the notation, let us again write the extended representation basis as a vector denoting the coefficients of each term in the representation basis, a = (a 0 , a 1 , . . ., a N , a N +1 ) T .As above, we now have to compute the projection of K Y k onto each of the R i .To ensure energy conservation, we will again determine µ using the energy integrals. 10This then yields the following system of equations that determines the representation vector a Y k : (2.20) The last equation is the energy conservation equation to determine the coefficient of M (x).We thus have a matrix equation of the form b K Y k = M R a K Y k , which we can solve for a K Y k given a finite representation basis.We again highlight the fact that the system was obtained using the energy conservation equation.The matrix M R is therefore again nearly equivalent to the full basis mixing matrix M R,ij = R i |R j .However, the last equation is replaced by the energy conservation equation, even if not explicitly distinguished in the notation.
As an example, if we choose R(x) = (Y (x), Y 1 (x), M (x)) T , we only have to determine the representations for Kx Y and Kx Y 1 .Solving the corresponding systems of equations then yields Looking at figure 2, we see that now the match with the exact result for K Y (x) is already very good.We note that adding more terms to the basis evidently changes all the coefficients of the solution, and also improves the result for the correspondence (see figure 2).

JCAP11(2023)026
Inserting eq.(2.21) back into eq.(2.7) and using η Y 1 ≈ 7.8246 in eq.(2.8), after collecting coefficients we then find In appendix B we give an alternative derivation that avoids the intermediate step of first representing the K Y k in terms of the basis.However, mathematically this is equivalent.By comparing the coefficients, one can again obtain a system for the evolution of Θ, y, y 1 and µ.The solution of this system now has the correct main properties.It conserves photon number and energy and leads to a solution11 y(y z ) y(0) e −3.8 yz .Indeed this is very close to the correct Green's function solution that neglects any residual distortion contributions.However, the precision can be improved by further extending the spectral basis (see figure 2).Below we will give the solutions for systems that include up to Y 15 (x) in the basis.This already provides a very accurate approximation for the exact themalisation Green's function.The related system can be readily generated using Mathematica following the procedure above.Schematically, we can then express the effect of the Kompaneets operator on the distortion in the form with y = (Θ, y, y 1 , . . ., y N , µ) T and where M K is the Kompaneets mixing matrix that directly depends on the chosen spectral basis. 12By construction, this matrix merely rotates y → y k → µ under energy conservation and thus has a zero first row and column.Even order systems are omitted, as they are found to be numerically unstable. 13We suspect this is due to the second order nature of the Kompaneets operator, but have no additional prove for this.The solution at any moment is then ∆n(x, y z ) ≈ B(x) • y(y z ) with the full spectral basis B(x) = (G(x), Y (x), Y 1 (x), . . ., Y N (x), M (x)) T .

Adding the effect of photon production and heating
To add the effect of photon production by double Compton (DC) and Bremsstrahlung (BR), we make use of the fact that once these become important, the y k will be extremely short-lived (i.e., decay quickly, y k → 0).In this case, we can neglect the role of the Y k 's for photon production and the analytic results for the µ-distortion evolution can be used [1,61].The net photon emission and absorption term has the explicit form [7,58,61]

JCAP11(2023)026
In the last step, we again linearised the problem with respect to the distortion [and Θ e O(∆n)].The DC and BR emissivities can be computed accurately using DCpack [71] and BRpack [72].As already explained in section 2.2, we can think of the effect that photon emission and absorption has on the distortion as a redistribution between µ and Θ. Overall this means as in eq.(2.5).This greatly simplifies the thermalisation problem, essentially converting the collision term into a source-sink term with built-in energy conservation.
To also add the effect of external heating, we assume that the distortions are generated through a y-distortion source, y = (1/4) Q , where Q = dQ/ dy z = ( τ θ z ) −1 dQ/ dt in this context.For energy release scenarios, this will be a very good approximation in the prerecombination era, since heat that is transferred to the baryons quickly reaches the photons through Compton scattering [e.g., see 7,59].The factor of α Y = 1/4 converts the change of the relative energy density into the y-parameter.Together we then have This equation now allows us to account for the effects of external heating and emission/absorption with the source vectors, Q and D, respectively.Refinements to the treatment of photon emission and absorption that include the effects of Y k (x) as well as other corrections to the leading order terms can in principle be added following the method of [61]; however, for now we stop with this simple description, emphasizing again that most of the thermalisation of distortions occurs deep into the µ-era, when these effects are expected to be small.It is worth noting that by formulating the heating as a vector in the extended spectral basis this assumption can be generalised and other spectral shapes can in principle be sourced at late times [see 73, for photon injection distortion evolution].This will not be explored in this paper -we adopt the standard thermalisation picture and study excitations of the y-distortion amplitude through energy injection as our benchmark in the following sections.

Solutions for the Green's function after single injection
The thermalisation Green's function has been successfully used to represent the spectral distortion shapes from continuous heating [50,53,74].With the above description we can reproduce the Green's function to high precision, as we show now.For this, we model the scenario of single injection in eq.(2.25) and introduce a narrow Gaussian heating rate at z injection (or alternatively set an effective initial condition for y at that redshift).Allowing this to evolve under successive scatterings we study the state of the system y(z f ) at the final redshift z f , and extract the corresponding spectral shape.
The result of this calculation for various basis sizes is shown in figure 4. Also shown is a comparison to the exact result of CosmoTherm, which performs an analogous calculation by directly binning the frequency space.The latter approach can be thought of as applying a "top JCAP11(2023)026 hat" basis in x to the same formalism discussed in section 2, and thus is more precise at the cost of tracking thousands of equations simultaneously.Despite the relative simplicity of the treatment derived here it is possible to capture the transition from temperature shift to y-distortion through the intermediate µ and residual eras accurately.The residual era in particular is captured by the expanded basis Y k , with N max = 9 already yielding very accurate results.
It is noteworthy that with the inclusion of the new spectral shapes Y k the definitions of y and µ have a degree of degeneracy.This is most notable in the recession of the µ-era with increasing N max and the existence of J y > 1 in the residual era (energy conservation is ensured by cancellation with the negative J y k ).This apparent arbitrary labelling of energy with different coefficients is not problematic, since the real physical observable that must converge is the spectrum, which indeed remains stable as seen in the right panels of figure 4.This physical observable will itself be projected onto some beneficial spectral shapes, as discussed in section 3, which depend on the characteristics of observing instrument [10,50] or other theoretical choices.
While the spectra in figure 4 show each snapshot being captured accurately, we note that the precise timings of the transition from one phase to the next appear slightly delayed (right) the different order of bases approach the residual distortion shape, but perform especially badly at low-frequencies.Also shown is a least-squares fit (red dotted line) using the Y 11 basis.
relative to the full CosmoTherm calculation. 14In figure 5 we show two time slices in the transition phases T → µ and µ → y, again with their respective CosmoTherm comparison and an optimised least squares fit to the full solution using the approximate basis of spectral functions.At z = 5 × 10 5 (left panel) we see the approximate solution having part of a temperature shift while the full solution is almost a pure µ distortion.Seeing that the optimised fit reproduces the CosmoTherm solution well, we conclude that the approximate treatment slightly overestimates the thermalisation timescale.As explained in [61], several additional aspects that are not captured by the simple treatment here do matter at the level of a few percent.By more carefully treating the DC and BR thermalisation rate, which will lead to a refined scaling of x c with time, one can probably improve the treatment; however, for our purpose the current approximation shall suffice, and refinements are left to future work.
In the right panel of figure 5, we also see a snapshot at z = 10 5 .It is apparent that the approximate basis approaches the full solution, but does not capture it fully.For comparison, an optimised fit (which allows one to smooth over any time-dependent mismatch) is shown but also fails to exactly reproduce the curve in this case.We can therefore conclude the basis only has enough freedom to capture some -but not all -of the nuances of the residual era.Departures from the solution are visible at low and intermediate frequencies (i.e., x 2) owing to the nature of the chosen basis (see figure 3) and our focus on energy conservation, which is driven by the high-frequency tail.Additional work on the optimal basis will likely remedy these limitations; however, we highlight that our treatment already greatly improves the modeling of the residual era, which is barely captured using a simple y + µ approximation.Hence, we again shall be content with the performance of the ODE treatment and focus on applications to anisotropic distortions as the main next step (paper II).

JCAP11(2023)026 3 Defining spectral distortion observables
In figure 4 we saw that the amplitude of the y distortion changes depending on the other amplitudes within the expanded basis while leaving the actual photon spectrum unchanged.In this section, we will formalize and further discuss this phenomenon in the context of changing basis.Heuristically, we can see the space of valid spectra as an abstract vector space, and as such choose a basis for this space.Provided the spectrum is a continuous function we expect a formal basis to be infinite, but computationally a finite basis can suffice, if chosen well.The bottom line statement we emphasise and highlight here then is that the underlying physics will be (and must be) independent of the choice of basis, where the physics here is captured only by the full photon spectrum and not any individual branching ratio or transfer function.
This gives the freedom to choose a basis which suits a given purpose most appropriately.In this light we will introduce two new bases, which are useful for packaging and exporting the results of the Y k basis or computation basis.At a given spectral sensitivity, only a finite number of spectral parameters will be directly measurable, and it is moot to attempt determining the amplitudes (or power and cross-power spectra) for all the spectral parameters inherent to the computation basis.The other two bases introduced here are guided by the principle to compress the information in the spectrum and prepare for easily extracting and interpreting the physics in observations.
Although from our discussion it is clear that a better computation basis which captures all the spectral complexity at low frequencies may exist, we are now interested in finding alternative representations for the space spanned by our Y k basis.As explained in [50], for a given experimental setting (e.g., frequency coverage and channel sensitivities) one can ask which spectral shapes are best constrained aside from the standard distortions.These spectral shapes can be determined using a principal component analysis.Mathematically, this can be thought of as an expansion of the spectrum into µ, y and Θ plus some additional spectral parameters, r i , to describe the residual distortion shapes. 15The residual distortion shapes are the principal spectral components spanning the residual distortion space and can be ranked by their observability, defining the observation basis.Denoting the residual distortion eigenspectra as S (k) , we find where ∆I ν = 2hν 3 /c 2 ∆n ν is the intensity corresponding to ∆n ν , which is integrated over the bandpass, B i (ν).In our computation we shall use a simple top-hat bandpass centered around frequency ν i with a width ∆ν i .Similarly, ∆I G i , ∆I Y i and ∆I µ i are the band-averaged versions of the corresponding G, Y and M intensities.The (band-averaged) residual distortion, ∆R i , space is orthogonal to M , Y and G, for the selected instrumental configuration.Since the binned spectral shapes can all be thought of as simple vectors, we can directly obtain the JCAP11(2023)026 µ, y, Θ and r k values for any distortion signal as This assumes that the covariance of the spectral bands is diagonal, but extensions can be readily given.The residual distortion parameters, by construction, will only receive contributions from the y i of our computation basis, while Θ o , µ o and y o will be a superposition of the Θ, y and µ values in the previous basis with extra contributions from the y i .The relevant rotation of the basis can be precomputed (see section 3.1).Given the observation basis S (k) we can therefore usually compress the information into fewer observational parameters, as we show below.
In figure 6, we show the first few S (k) used in our computations below.The basis was created assuming constant channel sensitivity and channel widths ∆ν = 1 GHz in the range ν min = 30 GHz to ν max = 1000 GHz, mainly for illustration.We normalized all of these to carry ∆ρ γ /ρ γ = 4 of energy.This choice makes them comparable in amplitude to the standard distortion shapes and the level of the corresponding residual distortion parameter gives away its relative importance.Creating the optimal distortion eigenmodes for more realistic experimental configurations is straightforward following the procedure outlined in [10,50].We can see that the distortion eigenmodes exhibit an increasing number of nodes, JCAP11(2023)026 As for figure but now with results into the "observation" basis.Notice how now the y and µ amplitudes are stable with increasing basis size.The spectra cover a smaller frequency range, as dictated by realistic observational scenarios, however they do not otherwise change compared to the computational basis.The residual era shows an effective negative temperature shift to achieve the correct spectral shape.reminiscent of other orthogonal functions sets.In applications, this will typically lead to the corresponding residual distortion parameter, r i , decreasing in amplitude.
In figure 7 we illustrate how this mapping to the observation basis modifies the appearance of the branching ratios (left panels) and photon spectra (right panels).We immediately see that the y distortion does not take on a relative energy contribution > 1 as it did in the computation basis (figure 4).Further to this point, the total amplitude of y and µ are more stable for increasing y n>0 , revealing that the modelling of spectral evolution is improving with basis size, but without the usual incurred coefficient ambiguity as a trade-off.The spectra cover a smaller frequency range, as discussed above, but otherwise show no significant departure from the result of the computational basis.Recall the statement that the bottom line physical result -the spectrum -is independent of the chosen basis.

Efficient change of the basis
To accelerate the calculation we can precompute all 'rotations' from one basis to the other given the distortion vectors (which depend on the experimental setting).Algorithmically, we have to bin all the involved spectra from the various bases and then compute the relevant JCAP11(2023)026  for various representations of the signal.The labels gives the maximal spectral component in the respective basis aside from the standard Θ, y and µ-description.The 'exact' result was obtained with the computation basis up to Y 15 .The simple Θ, y and µ-descriptions fails at the level of several tens of percent in particular at low frequencies.On the other hand, the distortion is extremely well represented once r 2 or r 3 are included.This is a compression of the information by a factor of more than 5.
mixing matrices and subsequently invert the problem.This then defines the mixing matrix L, which maps y = (Θ, y, y 1 , . . ., y N , µ) The dimension of the two spaces need not be the same, with the observation basis having a lower dimension given that the observability of various independent signal modes is usually reduced.
For our analysis, we pre-compute L for N = 15 and M = 6, but usually will only need r 1 , r 2 and r 3 to obtain a highly accurate representation of the full Y k basis result.Even for computations of power spectra, this significantly reduces the dimensionality of the problem (see paper III).As shown in figure 8, the residual distortion representation performs as well as the computation basis but with a lot fewer components (see discussion next section).

Performance and convergence
We are now in the position to compare the performance of the observation basis in representing the distortion solutions obtained by using the computational Y k basis.Two aspects are immediately worth noting: since the observation basis has a limited frequency coverage, it will not provide a description of the distortion solution outside this domain.This is analogous to having limited sky coverage, although there the properties of the spherical harmonic basis allows for some level of statistical deconvolution in CMB analyses [75].For the distortion spectra, this inversion problem will not be possible unless as many distortion parameters as basis parameters are observed accurately.
Second, the number of independently observable modes will depend on the frequency domain and frequency resolution as well as the sensitivity of the experiment.For example, it has been demonstrated that distinguishing µ-type distortion spectra benefits from having frequency channels below 30 GHz [44,76,77].However, a more comprehensive exploration of these dependencies aspects is beyond the scope of this paper, and for our illustrations we will stick to the modes shown in figure 6.
To illustrate the performance of the observation basis, we consider the distortion caused by a single energy injection at z h = 5 × 10 4 with ∆ρ γ /ρ γ = 10 −5 .In this regime, the residual distortion contributions are expected to be largest and hence the departures from the standard µ and y description are maximized.Looking at figure 8, we can immediately see that only the first few residual distortion spectra are needed to accurately represent the distortion shape at the level of a few percent of the dominant signal.This is a significant compression of the required information for the signal processing.However, it also implies that from the precise distortion shape not as much information can be directly extracted unless a very large distortion signal is present or extremely high sensitivity is achieved [50].

Caveats of the observation basis and alternative description
We would like to highlight a few important aspects of the observation basis.While by construction, Y , M and the Y k spectra do not carry photon number, the same is not true for the residual distortion spectra.This implies that in the new representation, not only the y and µ parameters change but also the temperature parameter is affected.Concretely, we have Θ ≈ 2.6 × 10 −12 , y ≈ 6.1 × 10 −6 and µ ≈ 3.0 × 10 −8 for the example shown in figure 8 in the Y 15 -representation.When projecting onto the observation basis we find Θ o ≈ −9.2 × 10 −7 , y o ≈ 1.3 × 10 −6 and µ o ≈ 1.3 × 10 −5 .We injected energy at a redshift where DC and BR are already very inefficient, such that Θ based on scattering physics alone should be negligibly small.After the change of basis, in particular µ o picks up a noticeable contribution and Θ o even drops below zero.This effect is known and originates from the fact that the residual distortion construction is based on intensity projections [see figure 2 of 50].The chosen procedure is most close to what would be obtained using standard component separation methods in future spectrometer analyses [e.g., 43,44,77].Although the total energetics of the problem and also the spectrum remain unchanged by the change of the representation, this behavior seems ambiguous.
An alternative observational procedure, without this apparent ambiguity, could be to fix the temperature parameter Θ based on the number density of the photon field.In this -21 -JCAP11(2023)026 case, one could fully orthogonalize G to the distortion space and construct a pure residual scattering distortion representation that is unaffected by the aforementioned effects.In figure 9, we show the result for the basis vectors in this alternative construction procedure.While generally very similar to the previous set of distortion modes (see figure 6), the alternative modes show a slightly differing pattern and overall trend.These modes can only be use in cases where the temperature contribution can be independently separated, as the modes no longer are orthogonal to G(x).
Figure 10 shows how this basis again stabilises the y and µ amplitudes across basis size while reproducing the same spectrum as the other bases.The difference however is that now the first residual mode constitutes a more dominant fraction of the energy, and the temperature shift never takes on its effective negative value.This is closer to the full scattering physics, since now number is conserved, but the description is not akin to a realistic observation of the sky.In paper III we will use these two bases wherever they are most illustrative, but always being careful and explicit.
In figure 11, we demonstrate that these alternative modes also represent the distortion shape very well, with the difference being that the contribution from G was fixed independently using a number density constraint.In the Y 15 representation, one has Θ ≈ 2.6 × 10 −12 , y ≈ 6.1 × 10 −6 and µ ≈ 3.0 × 10 −8 as before.Taking the full spectrum and imposing the photon number constraint to obtain the amplitude of G and then fitting for y and µ, we obtain Θ ≈ 2.6×10 −12 , y ≈ 1.6×10 −6 and µ ≈ 7.7×10 −6 .Just like before, we see a significant change in the values for y and µ, but this time no change to Θ.The total energy carried by G, Y and M is ∆ρ γ /ρ γ 4y + µ/1.4 ≈ 1.2 × 10 −5 , implying that the residual distortion contributes  4, but now results cast into the "scattering" basis.Notice that in contrast to figure 7 there is no production of a negative temperature shift, since here we enforce a strict number conservation of the residual modes, meaning no distortion shape can project onto G in the change of basis.Again the spectra show no change compared to the computation basis.∆ρ γ /ρ γ −0.2 × 10 −5 .In contrast, for the Y 15 representation we have ∆ρ γ /ρ γ ≈ 2.5 × 10 −5 stored in the G, Y and M components, implying that about ∆ρ γ /ρ γ ≈ −1.5 × 10 −5 is in the Y k>0 terms, which is no small total correction.If we compare all this to the lowest order computation using only Θ, y and µ in the ODE (i.e., a Y 0 -representation) we obtain Θ ≈ 1.1 × 10 −9 , y ≈ 1.6 × 10 −6 and µ ≈ 5.1 × 10 −6 .This demonstrates that the µ and y decomposition is well captured by the alternative distortion eigenmodes.
However, in the outlined alternative procedure an observer has to evaluate the number integral ∝ x 2 ∆n x dx of the photon field, which in experimental settings has several challenges.First, unless the distortion is measured in a sufficiently wide range of frequencies, this number integral would not evaluate accurately.Specifically, over a finite frequency domain even x 2 Y dx and x 2 M dx are no longer guaranteed to vanish, thereby breaking the 'photon number orthogonality'.Second, to carry out the integral, the frequency sampling has to be fine, which again usually runs into observational difficulties.Third, the estimation of errors will be non-standard since the observable is based on weighted sums of fluxes.Therefore, this approach is not expected to be realized in actual observations.Nevertheless, for theoretical calculations, we can use it for illustration when the focus is on the energetics of the problem.We will therefore refer to this alternative basis as scattering basis, given that it is constructed to focus on the spectral shapes that are introduced purely by Compton scattering terms, which conserve photon number.We further discuss the benefits and differences of changing the basis in paper III.

Discussion and conclusions
In this work we obtained an approximate ODE treatment for the thermalisation Green's function for heating processes, which captures most aspects of the full CosmoTherm calculation using an extended spectral basis to describe the residual distortion evolution (section 2).Instead of the expensive 'top-hat' frequency binning we use a spectral basis that is derived from boosts of the y-distortion spectrum.This reduces the computational burden by a factor of 10 3 , thereby providing one of the main steps towards formulating a generalised photon Boltzmann hierarchy, that will allow us to compute the evolution of distortion anisotropies at first order in perturbation theory (see papers II and III).We also clarify how the computational spectral basis can be compressed into fewer distortion shapes that -24 -

JCAP11(2023)026
can be distinguished with a given experimental configuration, introducing the observation and scattering basis (section 3).
The new ODE representation of the thermalisation Green's function given here is not perfect, because not all the spectral shapes that are excited by heating processes can be spanned by the basis functions we choose (see figure 5).However, we have demonstrated that this is not a severe problem for the average distortion evolution, which specifically relies on conversion of Y (x) [as the main distortion source] to M (x) and G(x).The representation of the full Green's function could probably be improved by studying the eigenfunctions of the Kompaneets and boost operators more carefully.In addition, weighting schemes and a modified truncation of the distortion basis could likely improve the performance.One could also refine the treatment of photon emission processes, including the effect of the new distortion shapes.This is expected to modify the thermalisation efficiency, a problem that may be solved perturbatively.Nevertheless, the novel ODE representation of the average thermalisation Green's function is sufficiently accurate for approximate applications to SD anisotropies caused by energy release, as we show in papers II and III.
It is very important to highlight that not all spectral distortion problems can be treated this way.First, if the distortions are large, a full non-linear thermalisation problem is encountered, which require several extensions [11,78].However, given the distortion limits by COBE/FIRAS, this situation is usually not relevant to the description of primordial distortions.Second, the Kompaneets equation itself does not describe the general scattering problem in an isotropic medium.At high energies, the full Compton kernel has to be applied, which can cause noticeable departures to the evolution [e.g., 63].Similarly, non-thermal distortions can in principle be excited by energy injection from particle cascades [78][79][80], which cannot be treated with the presented formalism.And finally, distortions created by direct photon injection can have a much richer spectral structure [73,81] that cannot be captured here.We leave these cases for future studies, noting that at least for linear problems extended spectral bases should provide a means forward.
Overall, this paper is the first in a series of works discussing the evolution of SD anisotropies generated by various physical mechanisms and how these might be constrained with future CMB spectrometers and imagers.The results from these works should open the path for more realististic SD anisotropy forecasts over a wide range of physics which previously were not possible.This will hopefully spur additional activity on CMB spectral distortions, uniting the efforts of CMB imaging and spectrometer approaches for probing the early Universe.

A Useful operator properties
We can somewhat reduce the complexity of the above calculations by studying the properties of -and relationships between -the operators Dx , D * x and Ôx .This will furthermore illustrate the suitability of the expanded basis {Y k (x)}.
We first note that two of the main operators commute with one another: Dx , Ôx = 0, (A.1) thus implying the potential existence of a shared eigenbasis.A logical step is to express the larger of the operators in terms of the other, revealing the identities The final equality can be easily found with a single application of the chain rule.However, it hints towards a more generic recurrence relation, which yields the following combinatoric sum: where the square brackets indicate Stirling set numbers, which counts partitions of an k-set into m nonempty subsets.This expansion of the boost operator reveals that the reverse operation is non-trivial -very specific weighted sums of x m ∂ m x terms are needed to make a boost operator with some power.Because of this, it is useful to be able to compose these expanded Ôx terms directly: As noted in [68], x k ∂ k x n bb has a recursion relation allowing for another combinatoric analytic solution where the angle brackets denote Eulerian numbers, defined as the number of permutations of the numbers 1 to m in which exactly k elements are greater than the previous element.This expression has very good convergence properties when summed starting from the highest power of e −mx .

JCAP11(2023)026
Defining H k (x) = (−x) k e −x /(1 − e −x ) k+1 , we are now in a position to write general expressions using the above formulae: showing that no such shared basis will exist, and thus for now we resort to the approximate numerical projections discussed in the main text (see especially figure 2).However, some more progress can be made my realising that D * x always appears in conjunction with the factor A = (1 + 2n bb ).It can be shown that A = This expression of the Kompaneets operator makes it clearer to see how certain results arise algebraically.Consider for example that Kx G = −Y , and similarly Kx M = −Y , where the latter result follows from the former together with Kx (G/x) = 0.The spectral shape, Y 1 , appears as intermediate step in these calculations, but ends up cancelling.These results may not be interesting in isolation, but they emphasise the fact that the cancellations only occur for simple shapes.Once you apply Kx to a distortion shape like Y you naturally get Y 1 and Y 2 that do not analytically cancel.

B Alternative derivation of the ODE system
To obtain the ODE system for the evolution of the spectrum, we can also directly project the evolution equation.Making the Ansatz ∆n ≈ B • y (with definitions as in the main section for a given basis) and then inserting this into the evolution equation, eq.where y = (y, y 1 , . . ., y k , µ) T and M R,ij = R i |R j is the full mixing matrix.We also have the source vector b Y,i = R i |R 0 = R i |Y and Kompaneets matrix As already explained in the main text, the system above will not yield a solution that correctly conserves energy (although it will become better and better the more Y k are included).We therefore replace the last row in the matrices M R and K and the last entry -28 -

JCAP11(2023)026
in b Y with the corresponding energy equation (as shown in the main text).The modified system has the same form as eq.(B.4), just with redefined matrices and vectors which we do not explicitly distinguish in the notation.The system can be solved for y to obtain the evolution equation for y, y 1 , . . ., y N , µ as The rows of the matrix M −1 R K are composed of the representation vectors for the operators K Y k .Note that the matrix K is an (N +2)×(N +1) matrix, while M R −1 is an (N +2)×(N +2) matrix, such that M −1 R K also is an (N +2)×(N +1) matrix.In addition we have M −1 R b Y = δ i0 , which simply follows from the fact that bY is the first column vector of the matrix M R .Since the matrix M −1 R K can be determined by independently solving for the representations of K Y k in terms of the representation basis R, this means we have proven the equivalence with the approach used in the main text.
.14b) This system is equivalent to the matrix equation b = M R a, where b i = R i |K Y for each function of the representation basis, i.e., R 0 = Y and R 1 = M in the considered case.Similarly, we have the basis mixing matrix M R,ij = R i |R j and the corresponding representation coefficients a 0 and a 1 .The solution is then a .15b) which can still be thought of as b = M R a, but with modified last rows in b and M R according to the energy conservation equation.By inverting the new system, this then yields the improved representation K Y (x) ≈ −3.4593 Y (x) − 10.871 M (x).Carrying out the energy integrals, we find E K Y ≈ −3.4593 × 4 − 10.871/1.4007= −21.598, in agreement with the direct integral result.

Figure 2 .
Figure 2. Distortion shapes x3 K Y (x) and x 3 K Y3 (x) for various approximations.Here the order refers to the largest term in Y k that is included [i.e., 0. order only Y (x) and M (x); 5. order includes M (x) and all Y k (x) up to Y 5 (x)].The representations become increasingly accurate the more terms we add to the basis.Typically only poor representation is obtained upon acting on the largest function in the basis -a problem which is mitigated by the fact that less energy occupy these higher modes in numerical solutions.

Figure 3 .
Figure 3. First few basis function Y k in comparison with the standard distortion shapes (dashed lines, with the positive peak from left to the right relating to G, M , Y , respectively).For increasing k the functions Y k occupy more and more of the high frequency part of the spectrum.

Figure 4 .
Figure 4.A figure showing the iterative improvements of augmenting the Y k basis.The rows show the branching ratios across redshifts (left) and final spectrum at various energy injection redshifts (right) for N max = 1, 5, 9. Gray lines in the branching ratio plots correspond to the Y k>0 coefficients, and dotted lines are negative values.Dotted black lines in the spectrum plots show the full results performed with CosmoTherm.

Figure 5 .
Figure 5.A figure showing two characteristic time slices where the full solution is poorly approximated.At z = 5 × 10 5 (left) the approximate basis still has an excessive contribution from G. At z = 10 5(right) the different order of bases approach the residual distortion shape, but perform especially badly at low-frequencies.Also shown is a least-squares fit (red dotted line) using the Y 11 basis.

Figure 6 .
Figure 6.First three residual distortion eigenmodes obtained for ∆ν = 1 GHz in the range ν min = 30 GHz to ν max = 1000 GHz.These signals are orthogonal to standard G, Y and M spectra and also among each other.They have all been normalized to carry an energy of ∆ρ γ /ρ γ = 4.

Figure 7 .
Figure 7.As for figure but now with results into the "observation" basis.Notice how now the y and µ amplitudes are stable with increasing basis size.The spectra cover a smaller frequency range, as dictated by realistic observational scenarios, however they do not otherwise change compared to the computational basis.The residual era shows an effective negative temperature shift to achieve the correct spectral shape.

1 2x 10 6 Θ, y and µ only r 1 r 2 r 3 Figure 8 .
Figure 8. Distortion (i.e., x 3 ∆n x ) after a single injection at z h = 5 × 10 4 with ∆ρ γ /ρ γ = 10 −5for various representations of the signal.The labels gives the maximal spectral component in the respective basis aside from the standard Θ, y and µ-description.The 'exact' result was obtained with the computation basis up to Y 15 .The simple Θ, y and µ-descriptions fails at the level of several tens of percent in particular at low frequencies.On the other hand, the distortion is extremely well represented once r 2 or r 3 are included.This is a compression of the information by a factor of more than 5.

Figure 9 .
Figure 9. First three residual distortion eigenmodes obtained for ∆ν = 1 GHz in the range ν min = 30 GHz to ν max = 1000 GHz and with an explicit photon number constraint to deproject G.These signals are orthogonal to the standard Y and M spectra and also among each other, but no longer leave the amplitudes of G unaltered.They have all been normalized to carry an energy of ∆ρ γ /ρ γ = 4.

Figure 10 .
Figure10.As for figure4, but now results cast into the "scattering" basis.Notice that in contrast to figure7there is no production of a negative temperature shift, since here we enforce a strict number conservation of the residual modes, meaning no distortion shape can project onto G in the change of basis.Again the spectra show no change compared to the computation basis.

6 Θ, y and µ only r 1 r 2 r 3 Figure 11 .
Figure 11.Same as figure 8 but using a photon number constraint to obtain the value for Θ.

1 x 4
+ Y G .This loose factor of 1/x combines nicely with the commutators noted above.Specifically we can then write D *x A = −( Ôx − 3) 4 + the misbehaving part of the Kompaneets operator to the previously named y-weight factor w y = Y /G.
(2.7), we haveΘ G(x) + N k=0 y k Y k (x) + µ M (x) = Θ e Y (x) − Θ Y (x) + N k=0 y k K Y k (x) − µ η M Y (x).(B.1)Here, Y 0 ≡ Y , K Y k = Kx Y k and we used the identities in eq.(2.12).Since only G(x) carries number we immediately obtain Θ = 0 by carrying out the number integral x 2 dx over this equation.Since we know that −Θ Y (x) on the right hand side of eq.(B.1) cancels the corresponding term in the Compton equilibrium temperatureΘ e ≈ x 3 w y (x) B dx 4E n bb • y ≈ Θ + N k=0 η Y k y k + η M µ, (B.2)and because there also is no term ∝ G(x), we only have to worry about the reduced problemN k=0 y k Y k (x) + µ M (x) = (Θ e − Θ − η M µ) Y (x) + N k=0 y k K Y k (x).(B.3)By performing the projections onto all function of the representation basis R = (Y, Y 1 , . .., Y k , M ) T , we obtain the system M R y = (Θ e − Θ − µ η M ) b Y + K y, (B.4) 21.598 by direct integration of the exact function and E K Y ≈ −8.8169 × 4 + 40.409/1.4007= −6.4178from the approximation.
JCAP11(2023)026 grant ST/T506291/1.JC was furthermore supported by the Royal Society as a Royal Society University Research Fellow at the University of Manchester, U.K. (No.URF/R/191023).AR acknowledges support by the project "Combining Cosmic Microwave Background and Large Scale Structure data: an Integrated Approach for Addressing Fundamental Questions in Cosmology", funded by the MIUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2017 -grant 2017YJYZAH.