January 01, 1970
In this work, we present a novel and more efficient approach to constructing the relativistic BKW (Bobylev, Krook, and Wu) solution. By introducing a class of ansatz functions for the distribution function, we demonstrate that within this specific ansatz space, only the equilibrium and BKW-type forms yield exact solutions to the nonlinear Boltzmann equation. Furthermore, guided by physical insight and drawing upon the framework of relativistic kinetic theory, we derive the Lie algebra of invariant transformations admitted by the relativistic Boltzmann equation. From this algebra, the corresponding symmetry group transformations can be systematically constructed.
The relativistic Boltzmann equation serves as a fundamental framework for describing the kinetic behavior of dilute relativistic gases in non-equilibrium regimes and occupies a central role in relativistic kinetic theory. It provides a statistical description of how the distribution function of particles evolves under the influence of free streaming and particle collisions, making it particularly suitable for systems where deviations from local thermal equilibrium are significant. Due to its solid theoretical foundation, the equation has found broad applications across diverse domains of modern physics, ranging from high-energy heavy-ion collisions [1]–[7] to astrophysical phenomena and cosmology [8], [9]. In these contexts, it enables the study of non-equilibrium transport processes under extreme conditions, such as the quark-gluon plasma produced in ultra-relativistic heavy-ion collisions and in the early universe, thereby offering a crucial link between microscopic dynamics and macroscopic observables.
In 1872, Boltzmann first formulated a kinetic equation describing the time evolution of the single-particle distribution function \(f\) for dilute gas molecules [10], establishing the H-theorem to demonstrate the monotonic decrease of the H function in isolated systems, which provided a microscopic foundation for the macroscopic law of entropy production. Due to the foundational role of the Boltzmann equation in statistical physics, the search for its analytical solutions has remained a central task for theoretical physicists since then. However, the mathematical complexity of the collision integral and the strong nonlinearity inherent in the equation make the derivation of exact analytical solutions extremely difficult. For nearly a century thereafter, the only well-known exact solution of the Boltzmann equation was the equilibrium-state solution. From the perspective of solving the equation, this solution is undoubtedly trivial, as it can be directly constructed from the moment the equation is written, based on the principle of collisional invariance. This situation remained unchanged until the discovery of the BKW solution, the first non-trivial exact solution to the Boltzmann equation.
In 1976, Bobylev discovered a set of particular solutions from the self-similar solutions of the Fourier-transformed Boltzmann equation [11]. Almost simultaneously, Krook and Wu also found this particular solution using different methods [12], [13], leading to its designation as the BKW solution 1. This solution describes the nonlinear relaxation dynamics of a spatially homogeneous, non-relativistic gas, offering a nice interpretation of the behavior of the one-particle distribution function at the high-momentum tail of the spectrum [12]. The BKW solution is not only a paradigmatic example of non-equilibrium relaxation but also provides deep insight into the interplay between collisional dissipation and the restoration of thermodynamic equilibrium. Due to its exact analytical form, it serves as an indispensable benchmark for rigorously testing theoretical approximations — such as moment methods and model kinetic equations — as well as numerical schemes like the DSMC (Direct Simulation Monte Carlo) method [16] and other numerical Boltzmann solvers [6], particularly in far-from-equilibrium regimes. Furthermore, the BKW solution has been further extended to relativistic regimes recently, including expanding solutions in FLRW (Friedmann-Lemaitre-Robertson-Walker) spacetime [17], [18] and solutions with anisotropic scattering cross sections [19]. In this work, we employ a novel approach to reconstruct the relativistic BKW solution.
Given the known exact solutions to the Boltzmann equation, a natural question emerges: whether it is possible to construct or discover additional exact solutions of physical significance by building upon these known solutions. Historically, this question is intimately tied to the Lie group analysis of PDE (partial differential equation). From a mathematical perspective, a central theme is the classification of PDE solutions, where it is often found that solutions with apparently different forms are connected via nontrivial symmetry transformations. This idea is intuitively compelling and can be rigorously formulated using Lie group and Lie algebra methods: any transformation that leaves the form of the original PDE invariant will necessarily map one solution to another. Let us now return to a concrete discussion of the Boltzmann equation. As a complex integro-differential equation, it can still be subjected to Lie group and Lie algebra analysis. Such a symmetry-based approach offers a promising pathway to discovering new physically relevant exact solutions. Even if new desired solutions are not immediately found, this framework may help establish connections between existing ones, thereby uncovering deeper physical relationships hidden within the structure of the equation hopefully. Although the Lie algebra of the Boltzmann equation is of significant importance and could offer valuable insights, its systematic and comprehensive exploration within the relativistic kinetic theory has been lacking. Within this work, we aim to close this gap.
The paper is structured as follows: In Sec. 2, we provide a brief introduction to the relativistic Boltzmann equation, and reconstruct the BKW solution using an alternative yet straightforward approach. Section 3 is devoted to determining the Lie algebras admitted by the equation. Summary and outlook are given in Sec. 4. Natural units \(k_{\mathrm{B}} = c = \hbar = 1\) are used. The metric tensor is represented by \(g_{\mu \nu} = \mathrm{diag}(1,-1,-1,-1)\). We also use the abbreviation \(\mathrm{d}P\) to denote the Lorentz invariant integral measure for on-shell massless particles \(\int \mathrm{d}P \equiv \frac{2}{(2 \pi)^3} \int \mathrm{d}^4 p \theta (p^0) \delta(p^2)\).
In the absence of external fields, the spacetime evolution of a non-equilibrium Bose or Fermi gas is described by the relativistic Boltzmann equation \[\begin{align} \begin{aligned} p \cdot \partial f(x^{\mu} ,p^{\mu}) = C[f], \label{21} \end{aligned} \end{align}\tag{1}\] with the collision kernel defined as \[\begin{align} \begin{aligned} C[f] \equiv &\frac{1}{2} \int \mathrm{d}K \mathrm{d}P_f \mathrm{d}K_f \left(f(x^{\mu}, p_f^{\mu}) f(x^{\mu}, k_f^{\mu}) (1 + a f(x^{\mu}, p^{\mu})) (1 + a f(x^{\mu}, k^{\mu})) \right.\\ &\left. - f(x^{\mu}, p^{\mu}) f(x^{\mu}, k^{\mu}) (1 + a f(x^{\mu}, p_f^{\mu})) (1 + a f(x^{\mu}, k_f^{\mu}))\right) W_{p,k \to p_f,k_f}, \label{22} \end{aligned} \end{align}\tag{2}\] where \(f(x^{\mu}, p^{\mu})\) is the single-particle distribution function in phase space and only binary elastic collisions are taken into account. The parameter \(a\) is a coefficient labeling different statistics with \(a = 1 (-1)\) corresponding to the effect of Bose enhancement (Fermi blocking), while \(a = 0\) follows from the classical Maxwell-Boltzmann statistics. The transition rate \(W_{p,k \to p_f,k_f}\) can be expressed as \(W_{p,k \to p_f,k_f} = (2 \pi)^6 s \sigma(s,\Theta) \delta ^4 (p^\mu + k^\mu - p_f^\mu - k_f^\mu)\) with the differential cross-section \(\sigma(s,\Theta)\) depending on the scattering angle \(\Theta\) and the Mandelstam variable \(s\equiv\left(p+k\right)^2\). For convenience, in the following we assume that the angular dependence of the differential cross-section is solely through the cosine of the scattering angle, i.e., \(\cos \Theta \equiv \mu\). Accordingly, the differential cross-section is expressed as \(\sigma(s,\mu)\). Note that Eq.@eq:21 applies specifically to single-component systems; for extensions to multi-component systems, see [20], [21]. For the majority of this work, we focus on the case \(a=0\) , i.e., the Boltzmann equation under classical Maxwell–Boltzmann statistics, including a re-derivation of the BKW solution using a novel method. In the following sections, we will see that, when analyzing the symmetries of the equation, the classical Boltzmann equation exhibits a richer symmetry structure compared to its quantum counterparts. Informally speaking, it is “more symmetric”.
When reconstructing exact analytical solutions of the Boltzmann equation, we proceed by imposing the assumptions of spatial homogeneity, i.e., \(f = f(x^0, p^{\mu})\). Subsequently, Eq.(1 ) becomes \[\begin{align} \begin{aligned} p^0 \partial_0 f(x^0 ,p^{\mu}) = C[f], \label{23} \end{aligned} \end{align}\tag{3}\] where \(C[f]\) can be rewritten as \[\begin{align} \begin{aligned} C[f] & = \frac{(2 \pi)^6}{2} \int \mathrm{d}K \mathrm{d}P_f \mathrm{d}K_f s \sigma(s,\mu) \delta ^4 (p^\mu + k^\mu - p_f^\mu - k_f^\mu) \left(f(x^0, p_f^{\mu}) f(x^0, k_f^{\mu}) - f(x^0, p^{\mu}) f(x^0, k^{\mu})\right), \label{24} \end{aligned} \end{align}\tag{4}\] where the particles are taken to be massless so that \(p^0=p=|\boldsymbol{p}|\). Unless otherwise specified, we assume throughout that the particle mass is zero.
For later convenience, we present here a frequently used parametrization of the cross section in our following discussion \[\begin{align} \begin{aligned} \sigma(s, \mu) = s^{l - 1} g_l(\mu) \label{25} \end{aligned}. \end{align}\tag{5}\] We will see that this separable form of the cross-section—depending on the Mandelstam variable \(s\) and the scattering angle independently—plays a crucial role in analyzing the scaling behavior of the Boltzmann equation. In fact, our consideration extends beyond this specific form: any cross-section that respects a scale transformation of the form \(\sigma \rightarrow \lambda^w \sigma\) (in the case of Eq.@eq:25 , \(w=2(l-1)\)) under momentum rescaling \(p^\mu\rightarrow \lambda p^\mu\) falls within the scope of our interest. A more comprehensive discussion of such scale-covariant cross-sections will be provided in the following sections. It is rather straightforward to see that this parametrization includes two widely encountered scattering interactions as special cases: when \(g_l(\mu) = const\), \(l = 1\) corresponds to the hard-sphere interaction, while \(l = 0\) corresponds to the leading order scalar \(\phi^4\) theory.
In general, directly solving for the generators and their corresponding transformations from the definition is a cumbersome mathematical task. The standard Lie symmetry analysis typically leads to a set of highly involved, nonlinear PDEs, which are often prohibitively difficult to solve explicitly (see Section 3). Rather than pursuing brute-force computation, we adopt a physics-informed strategy, guided by accumulated knowledge or physical insights in kinetic theory, to determine as much as possible about the relevant Lie algebra. For instance, knowledge of the various closed-form solutions of the Boltzmann equation offers helpful insights that guide the construction of the Lie algebra. Below, we briefly review and collect the known exact solutions for reference and further analysis.
When the Boltzmann equation is formulated, it is already known that the equilibrium solution constitutes a special class of exact solutions. These solutions are of fundamental physical significance, as they form the basis of equilibrium statistical mechanics. Based on the principle of collisional invariance in microscopic scattering processes, they can be systematically constructed [20]. Depending on the value of \(a\) , these equilibrium distributions take the following forms: \[\begin{align} f(p^\mu)& = e^{- \frac{E_p-\mu_c}{T}},\quad a=0, \tag{6} \\ f(p^\mu)& = \frac{1}{e^{\frac{E_p - \mu_c}{T}} + 1} ,\quad a=-1, \tag{7} \\ f(p^\mu)& = \frac{1}{e^{\frac{Ep - \mu_c}{T}} - 1},\quad a=1, \tag{8} \end{align}\] where \(\mu_c\) is the chemical potential and \(E_p\equiv u\cdot p\) with a timelike four-vector \(u^\mu\) which is often interpreted as the fluid velocity applied to a fluid system. In the above expressions, we have suppressed the normalization factors for simplicity.
The second solution is the non-thermal fixed point (NTFP) solution to the spatially homogeneous and momentum-isotropic Boltzmann equation, which has attracted significant research interest in recent years [22], [23]. In the over-occupied regime of a Bose gas, where the distribution function \(f\gg 1\), the statistical factors in the collision term of the Boltzmann equation become \[\begin{align} &f(t, p_f^0) f(t, k_f^0) \left(1 + f(t, p^0)\right) \left(1 + f(t, k^0)\right) - f( t,p^0) f(t, k^0) \left(1 + f(t, p_f^0)\right) \left(1 + f(t, k_f^0)\right)\nonumber \\ \approx\,\, & f(t, p_f^0) f(t, k_f^0) \left(f(t, p^0)+ f(t, k^0)\right) - f( t,p^0) f(t, k^0) \left( f(t, p_f^0)+ f(t, k_f^0)\right).\label{29} \end{align}\tag{9}\] The above simplification is crucial to the analysis of the scaling behavior of NTFP, and we shall come to this point later. The Boltzmann equation with the replacement of Eq.@eq:29 admits a universal self-similar solution, i.e., the transport in isolated systems is well described in terms of a self-similar evolution, \[\begin{align} \begin{aligned} f(t, p^0) & = t^{\tilde{a}} f_s (\xi\equiv t^{\tilde{b}} p^0), \label{210} \end{aligned} \end{align}\tag{10}\] in a given scaling regime. The scaling exponents \(\tilde{a}\) and \(\tilde{b}\), as well as the functional form of the non-thermal fixed point distribution \(f_s(\xi)\), are universal [24]. The above solution exhibits a typical inverse particle cascade, which leads to Bose condensation in the highly occupied low-momentum regime. This particle transport toward low momenta is part of a dual cascade, in which energy is simultaneously transferred to higher momenta via a direct cascade. In this regime, distinct scaling exponents \({\tilde{a}}^\prime,{\tilde{b}}^\prime\) and a different scaling function \(f^\prime_s\) as compared to the infrared regime are observed. We refer interested readers to [22]–[24] for further technical details. Below, we would like to give a few remarks:
First, strictly speaking, the NTFP is not an analytical solution, as the explicit form of \(f_s\) cannot be determined analytically. Nor can it be considered an exact solution of the full Boltzmann equation, as it relies on simplifications or approximations such as Eq.@eq:29 . Nevertheless, it captures an approximately universal dynamical behavior, making it physically compelling. Notably, this self-similar solution is relevant for a wide range of applications from ultracold quantum gases to high-energy particle physics, and the existence of nonthermal fixed points acts as a nonequilibrium attractor for isolated many-body systems far from equilibrium [25].
Second, the NTFP is a fixed point distinct from thermal equilibrium. Its existence implies that the fate of a weakly coupled system’s evolution need not be thermalization—within physically relevant time scales. For instance, in the infrared regime, the inverse particle cascade eventually leads to the onset of Bose condensation.
Third, the NTFP may share a deep connection with another exact analytical solution of the Boltzmann equation with classical statistics to be discussed later. At first glance, both are self-similar solutions; however, a more profound relationship may exist, which goes beyond superficial similarity and warrants further investigation. However, it should be noted that only the behavior of the NTFP solution in the ultraviolet regime is truly comparable, as in this regime the scaling requires \(f\ll 1\), causing the transport equation to reduce to the classical Boltzmann equation.
The last exact analytical solution is the relativistic BKW solution of the Boltzmann equation. Notably, this non-equilibrium particular solution possesses a fixed point that corresponds to the equilibrium solution of the Boltzmann equation [19].
In this subsection, we will construct the relativistic BKW solution using a novel approach that is simpler and more efficient than the conventional moment method [17]–[19], as it only requires evaluating the collision integral and solving a system of ordinary differential equations. It should be noted that the relativistic BKW solution applies only to the case where \(\sigma(s, \mu) = g_l(\mu)\) [19], corresponding to \(l=1\) in Eq.@eq:25 . For simplicity, we consider the hard-sphere interaction, i.e., \(\sigma(s, \mu) \equiv \sigma = const\). Before proceeding to the explicit calculation, let us first simplify the collision integral as much as possible [26], [27].
In deriving the BKW solution, we further assume isotropy in momentum space. The collision term in Eq.(4 ) can be written as \[\begin{align} \begin{aligned} C[f] & = \frac{\sigma_\mathrm{T}}{2(2 \pi)^4} \int \frac{\mathrm{d}^3 \boldsymbol{k}}{k^0} \frac{\mathrm{d}^3 \boldsymbol{p}_f}{p_f^0} \frac{\mathrm{d}^3 \boldsymbol{k}_f}{k_f^0} s \delta ^4 (p^\mu + k^\mu - p_f^\mu - k_f^\mu)(f(x^0 ,p_f^0) f(x^0 ,k_f^0) - f(x^0 ,p^0) f(x^0 ,k^0)), \label{211} \end{aligned} \end{align}\tag{11}\] where \(\sigma_T \equiv 2 \pi \sigma\) is the total cross section. We can re-express the phase-space integral by introducing the definitions of the three-momentum difference \(\boldsymbol{q} \equiv \boldsymbol{p} - \boldsymbol{p}_f = \boldsymbol{k}_f - \boldsymbol{k}\) and the energy difference \(\omega \equiv p^0 - p_f^0 = k_f^0 - k^0\). Then the integral over the delta function can be cast into \[\begin{align} \begin{aligned} & \int \frac{\mathrm{d}^3 \boldsymbol{k}}{k^0} \frac{\mathrm{d}^3 \boldsymbol{p}_f}{p_f^0} \frac{\mathrm{d}^3 \boldsymbol{k}_f}{k_f^0} \delta ^4 (p^\mu + k^\mu - p_f^\mu - k_f^\mu) \\ & = \int \frac{\mathrm{d}^3 \boldsymbol{k}}{k^0} \frac{\mathrm{d}^3 \boldsymbol{p}_f}{p_f^0} \frac{\mathrm{d}^3 \boldsymbol{k}_f}{k_f^0} \mathrm{d} \omega\delta ^3 (\boldsymbol{p} + \boldsymbol{k} - \boldsymbol{p}_f - \boldsymbol{k}_f) \delta(p^0 - \omega - p_f^0) \delta(k^0 + \omega - k_f^0) \\ & = \int \mathrm{d}^3 \boldsymbol{k} \mathrm{d}^3 \boldsymbol{q} \mathrm{d} \omega \frac{1}{|\boldsymbol{k}| |\boldsymbol{p}_f| |\boldsymbol{q+k}|} \delta(|\boldsymbol{p}| - \omega - |\boldsymbol{p} - \boldsymbol{q}|) \delta(|\boldsymbol{k}| + \omega - |\boldsymbol{k} + \boldsymbol{q}|). \label{212} \end{aligned} \end{align}\tag{12}\]
Subsequently, the dependence on the vector modulus within the delta function is transformed into the angular dependence of the vectors \(\boldsymbol{p}\), \(\boldsymbol{q}\) and of \(\boldsymbol{k}\), \(\boldsymbol{q}\) \[\begin{align} \begin{aligned} \delta(|\boldsymbol{p}| - \omega - |\boldsymbol{p} - \boldsymbol{q}|) & = \delta\left(|\boldsymbol{p}_f| - (|\boldsymbol{p}| - \omega)\right) = \frac{|\boldsymbol{p}_f|}{|\boldsymbol{p}||\boldsymbol{q}|} \delta(\cos \theta_{\boldsymbol{pq}} - \frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|}) \theta(|\boldsymbol{p}| - \omega), \label{213} \end{aligned} \end{align}\tag{13}\] and \[\begin{align} \begin{aligned} \delta(|\boldsymbol{k}| + \omega - |\boldsymbol{k} + \boldsymbol{q}|) & = \delta\left(|\boldsymbol{k}_f| - (\omega + |\boldsymbol{k}|)\right) = \frac{|\boldsymbol{k}_f|}{|\boldsymbol{k}||\boldsymbol{q}|} \delta(\cos \theta_{\boldsymbol{kq}} - \frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|}) \theta(\omega + |\boldsymbol{k}|). \label{214} \end{aligned} \end{align}\tag{14}\]
Since \[\begin{align} |\boldsymbol{p}_f| &= |\boldsymbol{p}| - \omega = \sqrt{|\boldsymbol{p} - \boldsymbol{q}|^2} = \sqrt{\boldsymbol{p}^2 + \boldsymbol{q}^2 - 2|\boldsymbol{p}||\boldsymbol{q}| \cos \theta_{\boldsymbol{pq}}}, \tag{15} \\ |\boldsymbol{k}_f| &= \omega + |\boldsymbol{k}| = \sqrt{|\boldsymbol{k} + \boldsymbol{q}|^2} = \sqrt{\boldsymbol{k}^2 + \boldsymbol{q}^2 + 2|\boldsymbol{k}||\boldsymbol{q}| \cos \theta_{\boldsymbol{kq}}}, \tag{16} \end{align}\] we can further simplify the expression as \[\begin{align} \begin{aligned} & \quad \delta(|\boldsymbol{p}| - \omega - |\boldsymbol{p} - \boldsymbol{q}|) \delta(|\boldsymbol{k}| + \omega - |\boldsymbol{k} + \boldsymbol{q}|)\nonumber \\ & = \frac{|\boldsymbol{p}_f| |\boldsymbol{k}_f|}{|\boldsymbol{p}| |\boldsymbol{k}| |\boldsymbol{q}|^2} \delta(\cos \theta_{\boldsymbol{pq}} - \frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|}) \nonumber \delta(\cos \theta_{\boldsymbol{kq}} - \frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|}) \theta(|\boldsymbol{p}| - \omega) \theta(\omega + |\boldsymbol{k}|). \label{217} \end{aligned} \end{align}\tag{17}\]
The requirement that the two angles \(\theta_{\boldsymbol{pq}}\) and \(\theta_{\boldsymbol{kq}}\) fall within \([0, \pi]\), together with the spacelike nature of the momentum transfer, is incorporated into the collision term, giving rise to a piecewise structure governed by three-step theta functions, \[\begin{align} \begin{aligned} & \quad\quad\,\omega^2 - |\boldsymbol{q}|^2 \le 0 \to \theta(|\boldsymbol{q}| - |\omega|), \\ & -1 \le \cos \theta_{\boldsymbol{pq}} \le 1 \to \theta(|\boldsymbol{p}| - \frac{|\boldsymbol{q}| + \omega}{2}), \\ & -1 \le \cos \theta_{\boldsymbol{kq}} \le 1 \to \theta(|\boldsymbol{k}| - \frac{|\boldsymbol{q}| - \omega}{2}). \end{aligned} \end{align}\]
Ultimately, we obtain the desired form of the collision term, which can be explicitly written down as \[\begin{align} \begin{aligned} C[f] = &\frac{\sigma_\mathrm{T}}{2(2 \pi)^4} \int \mathrm{d}^3 \boldsymbol{k} \mathrm{d}^3 \boldsymbol{q} \mathrm{d} \omega \frac{1}{|\boldsymbol{p}| |\boldsymbol{k}|^2 |\boldsymbol{q}|^2} s \delta(\cos \theta_{\boldsymbol{pq}} - \frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|}) \delta(\cos \theta_{\boldsymbol{kq}} - \frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|}) \\ & \times \theta(|\boldsymbol{p}| - \omega) \theta(\omega + |\boldsymbol{k}|) \theta(|\boldsymbol{q}| - |\omega|) \theta(|\boldsymbol{p}| - \frac{|\boldsymbol{q}| + \omega}{2}) \theta(|\boldsymbol{k}| - \frac{|\boldsymbol{q}| - \omega}{2}) (f(x^0 ,p_f^0) f(x^0 ,k_f^0) - f(x^0 ,p^0) f(x^0 ,k^0)). \label{219} \end{aligned} \end{align}\tag{18}\]
The preparatory steps are now complete. In the following, we will demonstrate how to construct the BKW solution using the new method. We begin by introducing two dimensionless notations \(\tau \equiv T^3 \sigma_\mathrm{T} x^0\) and \(\hat{p}^ \mu \equiv \frac{p^ \mu}{T}\) where \(T\) is a typical energy scale and can be identified with the temperature in certain physical scenarios, and we will omit the hat notation from \(\hat{p}^\mu\) when nothing confusing occurs. Then the dimensionless Boltzmann equation takes the form of \[\begin{align} p^0 \partial_{\tau} f(\tau ,p^0) = C[f], \label{220} \end{align}\tag{19}\] with \[\begin{align} \begin{aligned} C[f] = & \frac{1}{2(2 \pi)^4} \int \mathrm{d}^3 \boldsymbol{k} \mathrm{d}^3 \boldsymbol{q} \mathrm{d} \omega \frac{1}{|\boldsymbol{p}| |\boldsymbol{k}|^2 |\boldsymbol{q}|^2} s \delta(\cos \theta_{\boldsymbol{pq}} - \frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|}) \delta(\cos \theta_{\boldsymbol{kq}} - \frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|}) \\ & \times \theta(|\boldsymbol{p}| - \omega) \theta(\omega + |\boldsymbol{k}|) \theta(|\boldsymbol{q}| - |\omega|) \theta(|\boldsymbol{p}| - \frac{|\boldsymbol{q}| + \omega}{2}) \theta(|\boldsymbol{k}| - \frac{|\boldsymbol{q}| - \omega}{2}) (f(\tau, p_f^0) f(\tau, k_f^0) - f(\tau, p^0) f(\tau, k^0)). \label{221} \end{aligned} \end{align}\tag{20}\]
We propose a class of ansatz distribution function with the following parametrization: \[\begin{align} \begin{aligned} f(\tau, p^0) = e^{- \frac{p^0}{\alpha(\tau)}} \sum_{i=0}^n A_i(\tau) \cdot (p^0)^i, \label{222} \end{aligned} \end{align}\tag{21}\] where \(A_i(\tau)\) are unknown functions of time \(\tau\) to be determined. Although this trial function does not encompass all possible parametrizations — hence it is not fully general — it is sufficiently flexible for our purposes. The proposed ansatz is not arbitrary, but guided by three physical and mathematical principles. First, to ensure analytical tractability, the distribution function is assumed to be a summation over functions separable in its variables, i.e., their dependence on \(\tau\) and \(p^0\) is factorized. This restricts its functional form to a series expansion in functions of \(p^0\), with coefficients that depend only on \(\tau\). Second, to guarantee regularity, ensuring that all kinetic moments are finite, a natural choice is to introduce an extra regulating function that suppresses potential ultraviolet divergences. Third, we invoke the principle of minimal complexity: in the absence of additional symmetry or dynamical constraints, the most economical choice is to expand the \(p^0\)-dependent part as a power series in \(p^0\), while taking the regulating function to be exponentially decaying (e.g., \(e^{-\frac{p^0}{\alpha(\tau)}}\)). Taken together, these considerations lead to the form adopted in this work. Substituting the above ansatz into both sides of the Boltzmann equation, we can determine the coefficients \(A_i\) by the method of undetermined coefficients—provided that the powers of \(p^0\) on the left- and right-hand sides match term by term. It is evident that not all values of \(n\) can satisfy this condition. One can count the powers of \(p^0\) on both sides of the equation: in the absence of subtle cancellations, the degrees no longer match for \(n>1\) (we will discuss these cases in Appendix.5). Only two exceptions exist: \(n=0\) corresponds to the equilibrium solution, while \(n=1\) yields a solution of BKW type. This may partly explain why the BKW solution is so special. Focusing on the case of \(n = 1\), Eq.(21 ) takes the following form \[\begin{align} \begin{aligned} f(\tau ,p^0) = e^{- \frac{p^0}{\alpha(\tau)}} \left(A_0(\tau) + A_1(\tau) p^0\right). \label{223} \end{aligned} \end{align}\tag{22}\]
Substituting it into the left-hand side of Eq.(19 ) leads us to \[\begin{align} \begin{aligned} e^{- \frac{p^0}{\alpha(\tau)}} \left(A_0'(\tau) p^0 + (A_1'(\tau) + \frac{A_0(\tau) \alpha'(\tau)}{\alpha^2(\tau)}) (p^0)^2 + \frac{A_1(\tau) \alpha'(\tau)}{\alpha^2(\tau)} (p^0)^3\right), \label{224} \end{aligned} \end{align}\tag{23}\] where the prime denotes the derivative of \(A(\tau)\) with respect to \(\tau\). Then the substitution of Eq.(22 ) into Eq.(20 ) gives \[\begin{align} \begin{aligned} C[f] = & \frac{1}{2(2 \pi)^4} \int \mathrm{d}^3 \boldsymbol{k} \mathrm{d}^3 \boldsymbol{q} \mathrm{d} \omega \frac{1}{|\boldsymbol{p}| |\boldsymbol{k}|^2 |\boldsymbol{q}|^2}e^{- \frac{|\boldsymbol{p}| + |\boldsymbol{k}|}{\alpha(\tau)}} A_1^2(\tau) \delta(\cos \theta_{\boldsymbol{p}\boldsymbol{q}} - \frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|}) \delta(\cos \theta_{\boldsymbol{k}\boldsymbol{q}} - \frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|}) \\ & \times \theta(|\boldsymbol{p}| - \omega) \theta(\omega + |\boldsymbol{k}|) \theta(|\boldsymbol{q}| - |\omega|) \theta(|\boldsymbol{p}| - \frac{|\boldsymbol{q}| + \omega}{2}) \theta(|\boldsymbol{k}| - \frac{|\boldsymbol{q}| - \omega}{2}) \\ & \times \left(\omega(|\boldsymbol{p}| - |\boldsymbol{k}|) - \omega^2\right) 2 |\boldsymbol{p}| |\boldsymbol{k}| \left(1 - (\cos \theta_{\boldsymbol{p}\boldsymbol{q}} \cos \theta_{\boldsymbol{k}\boldsymbol{q}} + \cos(\phi_{\boldsymbol{p}} - \phi_{\boldsymbol{k}}) \sin \theta_{\boldsymbol{p}\boldsymbol{q}} \sin \theta_{\boldsymbol{kq}})\right). \label{225} \end{aligned} \end{align}\tag{24}\] In obtaining the above expression, we use \(s = 2 |\boldsymbol{p}| |\boldsymbol{k}| (1 - \cos \theta_{\boldsymbol{p}\boldsymbol{k}})\) and \(\cos \theta_{\boldsymbol{p}\boldsymbol{k}} = \cos \theta_{\boldsymbol{p}\boldsymbol{q}} \cos \theta_{\boldsymbol{k}\boldsymbol{q}} + \cos(\phi_{\boldsymbol{p}} - \phi_{\boldsymbol{k}}) \sin \theta_{\boldsymbol{p}\boldsymbol{q}} \sin \theta_{\boldsymbol{kq}}\).
By transforming the integral measure to spherical coordinate, Eq.(24 ) can be rearranged to yield \[\begin{align} \begin{aligned} C[f] = & \frac{1}{(2 \pi)^4} e^{- \frac{|\boldsymbol{p}|}{\alpha(\tau)}} A_1^2(\tau) \int \mathrm{d}|\boldsymbol{k}| \mathrm{d} \cos \theta_{\boldsymbol{kq}} \mathrm{d} \phi_{\boldsymbol{k}} \mathrm{d}|\boldsymbol{q}| \mathrm{d} \cos \theta_{\boldsymbol{pq}} \mathrm{d} \phi_{\boldsymbol{q}} \mathrm{d} \omega \\ & \times |\boldsymbol{k}| \delta(\cos \theta_{\boldsymbol{pq}} - \frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|}) \delta(\cos \theta_{\boldsymbol{kq}} - \frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|}) \\ & \times \theta(|\boldsymbol{p}| - \omega) \theta(\omega + |\boldsymbol{k}|) \theta(|\boldsymbol{q}| - |\omega|) \theta(|\boldsymbol{p}| - \frac{|\boldsymbol{q}| + \omega}{2}) \theta(|\boldsymbol{k}| - \frac{|\boldsymbol{q}| - \omega}{2}) \\ & \times e^{- \frac{|\boldsymbol{k}|}{\alpha(\tau)}} (\omega(|\boldsymbol{p}| - |\boldsymbol{k}|) - \omega^2) (1 - (\cos \theta_{\boldsymbol{pq}} \cos \theta_{\boldsymbol{kq}} + \cos(\phi_{\boldsymbol{p}} - \phi_{\boldsymbol{k}}) \sin \theta_{\boldsymbol{pq}} \sin \theta_{\boldsymbol{kq}})) \\ = & \frac{1}{(2 \pi)^4} e^{- \frac{|\boldsymbol{p}|}{\alpha(\tau)}} A_1^2(\tau) \int \mathrm{d}|\boldsymbol{k}| \mathrm{d} \phi_{\boldsymbol{k}} \mathrm{d}|\boldsymbol{q}| \mathrm{d} \phi_{\boldsymbol{q}} \mathrm{d} \omega e^{- \frac{|\boldsymbol{k}|}{\alpha(\tau)}} |\boldsymbol{k}| (\omega(|\boldsymbol{p}| - |\boldsymbol{k}|) - \omega^2) \\ & \times \theta(|\boldsymbol{p}| - \omega) \theta(\omega + |\boldsymbol{k}|) \theta(|\boldsymbol{q}| - |\omega|) \theta(|\boldsymbol{p}| - \frac{|\boldsymbol{q}| + \omega}{2}) \theta(|\boldsymbol{k}| - \frac{|\boldsymbol{q}| - \omega}{2}) \\ & \times \left[1 - \left((\frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|}) (\frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|}) \right.\right.\\ &\left.\left. + \cos(\phi_{\boldsymbol{p}} - \phi_{\boldsymbol{k}}) \sqrt{1 - (\frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|})^2} \sqrt{1 - (\frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|})^2}\right)\right]. \label{226} \end{aligned} \end{align}\tag{25}\] Finally, by translating all step functions into constraints on the upper and lower limits of integration, we arrive at \[\begin{align} \begin{aligned} C[f] =& \frac{1}{(2 \pi)^2} e^{- \frac{|\boldsymbol{p}|}{\alpha(\tau)}} A_1^2(\tau) \left[ \int_{\frac{|\boldsymbol{q}| - \omega}{2}}^{+ \infty}\mathrm{d}|\boldsymbol{k}|\int_{-|\boldsymbol{q}|}^{+|\boldsymbol{q}|} \mathrm{d} \omega \int_0^{|\boldsymbol{p}|}\mathrm{d}|\boldsymbol{q}| e^{- \frac{|\boldsymbol{k}|}{\alpha(\tau)}} |\boldsymbol{k}| \left(\omega(|\boldsymbol{p}| - |\boldsymbol{k}|) - \omega^2\right) \right.\\ &\left. \times \left(1 - (\frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|}) (\frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|})\right) \right.\\ &\left. + \int_{\frac{|\boldsymbol{q}| - \omega}{2}}^{+ \infty}\mathrm{d}|\boldsymbol{k}| \int_{-|\boldsymbol{q}|}^{2|\boldsymbol{p}| - |\boldsymbol{q}|}\mathrm{d} \omega \int_{|\boldsymbol{p}|}^{+ \infty}\mathrm{d}|\boldsymbol{q}| e^{- \frac{|\boldsymbol{k}|}{\alpha(\tau)}} |\boldsymbol{k}| \left(\omega(|\boldsymbol{p}| - |\boldsymbol{k}|) - \omega^2\right) \right.\\ &\left. \times \left(1 - (\frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|}) (\frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|})\right)\right] \\ =& \frac{1}{6 \pi^2}\left[e^{- \frac{|\boldsymbol{p}|}{\alpha(\tau)}} A_1^2(\tau) |\boldsymbol{p}| \left(|\boldsymbol{p}| - 6 \alpha(\tau)\right) \left(|\boldsymbol{p}| - 2 \alpha(\tau)\right) \alpha^3(\tau)\right] \\ =& e^{- \frac{p^0}{\alpha(\tau)}} \left[\frac{2 A_1^2(\tau) \alpha^5(\tau)}{\pi^2} p^0 - \frac{4 A_1^2(\tau) \alpha^4(\tau)}{3 \pi^2} (p^0)^2 + \frac{A_1^2(\tau) \alpha^3(\tau)}{6 \pi^2} (p^0)^3\right]. \label{227} \end{aligned} \end{align}\tag{26}\] Remarkably, the collision integral — typically intractable in closed form — can be worked out analytically! This is undoubtedly a nontrivial result, as even a slight increase in interaction complexity—such as the inclusion of angular dependence in the scattering cross-section—would render the integral analytically intractable. After completing this nontrivial collision integral, we compare the powers of \(p^0\) on both sides of the Boltzmann equation, which allows us to obtain \[\begin{align} \begin{aligned} A_0'(\tau) & = \frac{2 A_1^2(\tau) \alpha^5(\tau)}{\pi^2}, \\ A_1'(\tau) + \frac{A_0(\tau) \alpha'(\tau)}{\alpha^2(\tau)} & = - \frac{4 A_1^2(\tau) \alpha^4(\tau)}{3 \pi^2}, \\ \frac{A_1(\tau) \alpha'(\tau)}{\alpha^2(\tau)} & = \frac{A_1^2(\tau) \alpha^3(\tau)}{6 \pi^2}. \label{228} \end{aligned} \end{align}\tag{27}\] Solving this system of ordinary differential equations directly is rather difficult. It should be noted that the particle number density \(n_0\) and the energy density \(e_0\) are the conserved quantities in the homogeneous case (see also [19] for details) with their definitions \[\begin{align} \begin{aligned} n_0 = \int \frac{\mathrm{d}^3 \boldsymbol{p}}{(2 \pi)^3} f(\tau,p^0), \qquad e_0 = \int \frac{\mathrm{d}^3 \boldsymbol{p}}{(2 \pi)^3} p^0 f(\tau,p^0). \label{229} \end{aligned} \end{align}\tag{28}\] These two conserved quantities provide the initial conditions necessary for solving the equations. Combining Eq.(28 ) with Eq.(22 ), \(A_0(\tau)\), \(A_1(\tau)\) can be expressed as a combination of \(n_0\), \(e_0\), and \(\alpha(\tau)\), \[\begin{align} \begin{aligned} A_0(\tau) = \frac{\pi^2 (-e_0 + 4 n_0 \alpha(\tau))}{\alpha^4(\tau)}, \quad A_1(\tau) = \frac{\pi^2 (e_0 - 3 n_0 \alpha(\tau))}{3 \alpha^5(\tau)}. \label{230} \end{aligned} \end{align}\tag{29}\]
Substituting Eq.(29 ) back into the first equation in Eq.(27 ) yields the following solution \[\begin{align} \begin{aligned} \alpha(\tau) = \frac{e_0}{3 n_0} + c e^{- \frac{n_0 \tau}{6}}, \label{231} \end{aligned} \end{align}\tag{30}\] where \(c\) is a constant independent of \(p^0\) and \(\tau\) acting as one free initial data. Note Eqs.@eq:230 and 30 are also consistent with the other two equations in Eq.@eq:228 . We also verify that our solutions reproduce the result given in [19] in the case of hard-sphere interaction. For comparison with [17], [18], we specify the initial condition by choosing the value of \(c\) as \(c = -\frac{1}{4}\), then the resulting expression reproduces the solution given in Ref.[17], [18] in Minkowski spacetime. A similar discussion on how the choice of initial conditions affects the final physical results can be found in Ref.[19]. It is also worthwhile to note that the present calculations can be readily extended to isotropically expanding FLRW spacetime [28].
Let’s first consider a general group transformation \[\begin{align} \begin{aligned} & f \to f' = \mathcal{A}(x^{\mu}, p^{\mu}, f, \theta), \\ & x^{\mu} \to x'^{\mu} = \mathcal{B}^{\mu}(x^{\mu}, p^{\mu}, f, \theta), \\ & p^{\mu} \to p'^{\mu} = \mathcal{C}^{\mu}(x^{\mu}, p^{\mu}, f, \theta), \label{31} \end{aligned} \end{align}\tag{31}\] with initial conditions \[\begin{align} \begin{aligned} & \mathcal{A}(x^{\mu}, p^{\mu}, f, 0) = f, \\ & \mathcal{B}^{\mu}(x^{\mu}, p^{\mu}, f, 0) = x^{\mu}, \\ & \mathcal{C}^{\mu}(x^{\mu}, p^{\mu}, f, 0) = p^{\mu}, \label{32} \end{aligned} \end{align}\tag{32}\] where \(\theta\) is the group parameter. Expanding the functions \(\mathcal{A}\), \(\mathcal{B}^{\mu}\), and \(\mathcal{C}^{\mu}\) in a Taylor series around \(\theta = 0\), while taking into account the initial conditions 32 , one can obtain infinitesimal transformations \[\begin{align} \begin{aligned} & f' \approx f + \theta \Omega(x^{\mu}, p^{\mu}, f), \\ & x'^{\mu} \approx x^{\mu} + \theta X^{\mu}(x^{\mu}, p^{\mu}, f), \\ & p'^{\mu} \approx p^{\mu} + \theta P^{\mu}(x^{\mu}, p^{\mu}, f). \label{33} \end{aligned} \end{align}\tag{33}\] where \[\begin{align} \begin{aligned} & \Omega(x^{\mu}, p^{\mu}, f) = \frac{\partial \mathcal{A}(x^{\mu}, p^{\mu}, f, \theta)}{\partial \theta} \big{|}_{\theta = 0}, \\ & X^{\mu}(x^{\mu}, p^{\mu}, f) = \frac{\partial \mathcal{B}^{\mu}(x^{\mu}, p^{\mu}, f, \theta)}{\partial \theta} \big{|}_{\theta = 0}, \\ & P^{\mu}(x^{\mu}, p^{\mu}, f) = \frac{\partial \mathcal{C}^{\mu}(x^{\mu}, p^{\mu}, f, \theta)}{\partial \theta} \big{|}_{\theta = 0}. \label{34} \end{aligned} \end{align}\tag{34}\] Subsequently, we are able to obtain the generator of the group \[\begin{align} \begin{aligned} \hat{L} = \Omega (x^{\mu}, p^{\mu}, f) \frac{\partial}{\partial f} + X^{\mu}(x^{\mu}, p^{\mu}, f) \frac{\partial}{\partial x^{\mu}} + P^{\mu}(x^{\mu}, p^{\mu}, f) \frac{\partial}{\partial p^{\mu}}. \label{35} \end{aligned} \end{align}\tag{35}\] If we have infinitesimal transformations 33 or generator 35 , the transformation 31 can be defined using the following Lie equations \[\begin{align} \begin{aligned} & \frac{\mathrm{d} \mathcal{A}}{\mathrm{d} \theta} = \Omega(\mathcal{A}, \mathcal{B}^{\mu}, \mathcal{C}^{\mu}), \qquad \mathcal{A} \big{|}_{\theta = 0} = f, \\ & \frac{\mathrm{d} \mathcal{B}^{\mu}}{\mathrm{d} \theta} = X^{\mu}(\mathcal{A}, \mathcal{B}^{\mu}, \mathcal{C}^{\mu}), \qquad \mathcal{B}^{\mu} \big{|}_{\theta = 0} = x^{\mu}, \\ & \frac{\mathrm{d} \mathcal{C}^{\mu}}{\mathrm{d} \theta} = P^{\mu}(\mathcal{A}, \mathcal{B}^{\mu}, \mathcal{C}^{\mu}), \qquad \mathcal{C}^{\mu} \big{|}_{\theta = 0} = p^{\mu}. \label{36} \end{aligned} \end{align}\tag{36}\]
We aim to obtain certain symmetry transformations, and then use them to map any solution of the Boltzmann equation back onto another solution of the same equation. That is, let Eq.@eq:31 be a symmetry transformation of Eq.@eq:21 , and let the function \[\begin{align} f = \Phi(x^{\mu}, p^{\mu}) \end{align}\] be a solution of Eq.@eq:21 . Since Eq.@eq:31 defines a symmetry transformation, the above solution can equivalently be expressed in terms of the transformed variables as \[\begin{align} f' = \Phi(x'^{\mu}, p'^{\mu}), \end{align}\] where \(f'\) satisfies the Boltzmann equation formulated with respect to the transformed coordinates and momenta \((x'^{\mu}, p'^{\mu})\). Substituting the transformation relations from Eq. 31 into the expression above, we finally obtain \[\begin{align} \begin{aligned} \mathcal{A}(x^{\mu}, p^{\mu}, f, \theta) = \Phi(\mathcal{B}^{\mu}(x^{\mu}, p^{\mu}, f, \theta), \mathcal{C}^{\mu}(x^{\mu}, p^{\mu}, f, \theta)). \label{37} \end{aligned} \end{align}\tag{37}\] If Eq.@eq:37 with respect to \(f\) is solved, we can then obtain a new solution of Eq.@eq:21 . Because the Boltzmann equation is a nonlinear integro-differential equation, its determining equations for Lie symmetries become highly nonlinear and analytically intractable. This difficulty motivates a shift in perspective: rather than solving the full system directly, we seek the corresponding Lie algebra following a problem-oriented and physics-informed approach. With these preparations in place, the following subsection will demonstrate how to achieve this.
The relativistic Boltzmann equation (1 ) is formally invariant under Poincaré transformations by construction. Therefore, it is natural to show that the following generators \[\begin{align} \begin{aligned} & \hat{L}^{(0)}_{ij} = (x_i \frac{\partial}{\partial x^j} - x_j \frac{\partial}{\partial x^i}) + (p_i \frac{\partial}{\partial p^j} - p_j \frac{\partial}{\partial p^i}), \\ & \hat{L}^{(1)}_i = \frac{\partial}{\partial x^i}, \\ & \hat{L}^{(2)}_i = (x_0 \frac{\partial}{\partial x^i} - x_i \frac{\partial}{\partial x^0}) + (p_0 \frac{\partial}{\partial p^i} - p_i \frac{\partial}{\partial p^0}), \\ & \hat{L}^{(3)} = \frac{\partial}{\partial x^0}. \label{38} \end{aligned} \end{align}\tag{38}\] These generators form a basic and closed Lie algebra. However, they still differ from the conventional Poincaré algebra. Specifically, the generators \(\hat{L}^{(0)}_{ij}\) and \(\hat{L}^{(2)}_i\) generate rotations and boost transformations in phase space, whereas spacetime translations remain confined to coordinate space. Since this algebra is trivially closed, no new elements can be generated through Lie brackets. This also implies that only by incorporating new physical insights can we obtain additional generators to construct a larger Lie algebra.
Recalling that in Sec.2.2, we have presented several analytical solutions of the Boltzmann equation, which could offer rich physical insights. If a particular solution exhibits a certain symmetry, that symmetry should be consistent with the structure of the Boltzmann equation. In other words, by examining the symmetries manifested in these physical solutions, we can uncover the underlying symmetries of the Boltzmann equation itself. The solution described by Eq.@eq:210 is precisely such a candidate that may provide valuable insight, as it respects the following scaling symmetry: If we perform the transformation \[t\rightarrow \lambda t,\quad p^0\rightarrow p^0\lambda^{-\tilde{b}}, \label{39}\tag{39}\] then the distribution function transforms accordingly as \(f\rightarrow \lambda^{\tilde{a}} f\). Inspired by the scaling symmetry of the NTFP, it is natural to ask whether there is a similar scaling symmetry respected by the Boltzmann equation. It is easy to draw a conclusion that the following scaling transformation generator \[\begin{align} \begin{aligned} \hat{L}^{(4)} = x^{\mu} \frac{\partial}{\partial x^{\mu}} - f \frac{\partial}{\partial f} \label{310} \end{aligned} \end{align}\tag{40}\] together with the resulting scale transformation \[\begin{align} x^\mu\rightarrow x^{\prime\mu}=e^\theta x^\mu,\quad p^\mu\rightarrow p^{\prime\mu}=p^\mu, \quad f\rightarrow f^\prime=e^{-\theta}f \label{311} \end{align}\tag{41}\] is the qualified one.
It is evident that under the scaling transformation generated by \(\hat{L}^{(4)}\), the momentum remains invariant, in contrast to the scaling transformation in Eq.@eq:39 . If the cross section transforms as \(\sigma \to \lambda^{2(r - 1)} \sigma\) under \(p^{\mu} \to \lambda p^{\mu}\), its specific form coincides with Eq. (5 ) when \(r = l\), which in turn allows us to identify another independent generator and its corresponding transformation \[\begin{align} \begin{aligned} & \hat{L}^{(5)} = v x^{\mu} \frac{\partial}{\partial x^{\mu}} - \frac{1}{2r + 1} p^{\mu} \frac{\partial}{\partial p^{\mu}} + (1 - v) f \frac{\partial}{\partial f}, \\ & x^\mu\rightarrow x^{\prime\mu}=e^{v\theta} x^\mu,\quad p^\mu\rightarrow p^{\prime\mu}=e^{-\frac{1}{2r+1}\theta}p^\mu, \quad f\rightarrow f^\prime=e^{(1-v)\theta}f. \label{313} \end{aligned} \end{align}\tag{42}\] Although \(\hat{L}^{(4)}\) and \(\hat{L}^{(5)}\) both appear to generate scaling transformations, they are distinct: \(\hat{L}^{(5)}\) cannot be reduced to \(\hat{L}^{(4)}\). As in the case of \(\hat{L}^{(4)}\), the inclusion of \(\hat{L}^{(5)}\) within the algebra generated by \(\hat{L}^{(0)}\) – \(\hat{L}^{(4)}\) results in a closed structure.
In the previous subsection, we derived the Lie algebra of invariant transformations for the relativistic Boltzmann equation. Here, we summarize these results and present a physical discussion. The complete set of generators forming the Lie algebra is presented below \[\begin{align} \begin{aligned} & \hat{L}^{(0)}_{ij} = (x_i \frac{\partial}{\partial x^j} - x_j \frac{\partial}{\partial x^i}) + (p_i \frac{\partial}{\partial p^j} - p_j \frac{\partial}{\partial p^i}), \\ & \hat{L}^{(1)}_i = \frac{\partial}{\partial x^i}, \\ & \hat{L}^{(2)}_i = (x_0 \frac{\partial}{\partial x^i} - x_i \frac{\partial}{\partial x^0}) + (p_0 \frac{\partial}{\partial p^i} - p_i \frac{\partial}{\partial p^0}), \\ & \hat{L}^{(3)} = \frac{\partial}{\partial x^0}, \\ & \hat{L}^{(4)} = x^{\mu} \frac{\partial}{\partial x^{\mu}} - f \frac{\partial}{\partial f}, \\ & \hat{L}^{(5)} = v x^{\mu} \frac{\partial}{\partial x^{\mu}} - \frac{1}{2r + 1} p^{\mu} \frac{\partial}{\partial p^{\mu}} + (1 - v) f \frac{\partial}{\partial f}. \label{314} \end{aligned} \end{align}\tag{43}\] The range of applicability of the transformations generated by these generators varies. The generators \(\hat{L}^{(0)}\) – \(\hat{L}^{(3)}\) form a Poincaré algebra, which applies to Eq.@eq:21 , regardless of whether quantum statistics are included or whether the particle mass is zero. In contrast, \(\hat{L}^{(4)}\) is applicable exclusively to classical systems and remains valid independently of the particle mass. On the other hand, \(\hat{L}^{(5)}\) is valid only for massless classical systems, where the cross-section exhibits scaling behavior \(\sigma \to \lambda^{2(r - 1)} \sigma\) under the transformation \(p^{\mu} \to \lambda p^{\mu}\).2 Furthermore, it can be shown that the entire set of generators given by Eq.@eq:314 is closed under the Lie bracket, as summarized in TABLE 1. To the best of our knowledge, this is the first time that a closed Lie algebra of invariant transformations for the relativistic Boltzmann equation has been presented.
| \(\hat{L}^{(0)}_{ij}\) | \(\hat{L}^{(1)}_i\) | \(\hat{L}^{(2)}_i\) | \(\hat{L}^{(3)}\) | \(\hat{L}^{(4)}\) | \(\hat{L}^{(5)}\) | |
|---|---|---|---|---|---|---|
| \(\hat{L}^{(0)}_{mn}\) | \(g_{jm} \hat{L}^{(0)}_{in} + g_{in} \hat{L}^{(0)}_{jm} + g_{jn} \hat{L}^{(0)}_{mi} + g_{im} \hat{L}^{(0)}_{nj}\) | \(\backslash\) | \(\backslash\) | \(\backslash\) | \(\backslash\) | \(\backslash\) |
| \(\hat{L}^{(1)}_m\) | \(g_{jm} \hat{L}^{(1)}_i - g_{im} \hat{L}^{(1)}_j\) | \(0\) | \(\backslash\) | \(\backslash\) | \(\backslash\) | \(\backslash\) |
| \(\hat{L}^{(2)}_m\) | \(g_{jm} \hat{L}^{(2)}_i - g_{im} \hat{L}^{(2)}_j\) | \(-g_{im} \hat{L}^{(3)}\) | \(\hat{L}^{(0)}_{mi}\) | \(\backslash\) | \(\backslash\) | \(\backslash\) |
| \(\hat{L}^{(3)}\) | \(0\) | \(0\) | \(- \hat{L}^{(1)}_i\) | \(0\) | \(\backslash\) | \(\backslash\) |
| \(\hat{L}^{(4)}\) | \(0\) | \(\hat{L}^{(1)}_i\) | \(0\) | \(\hat{L}^{(3)}\) | \(0\) | \(\backslash\) |
| \(\hat{L}^{(5)}\) | \(0\) | \(v \hat{L}^{(1)}_i\) | \(0\) | \(v \hat{L}^{(3)}\) | \(0\) | \(0\) |
The given Lie algebra formed by \(L^{(0)}\) to \(L^{(5)}\) has a one-to-one correspondence with its non-relativistic counterpart [29], [30] except for the following element \(\hat{L}^{(6)}\). Specifically, the non-relativistic Boltzmann equation admits an additional generator, \(\hat{L}^{(6)}\), associated with a nonlinear transformation, \[\begin{align} \begin{aligned} & \hat{L}^{(6)} = t^2 \frac{\partial}{\partial t} + t \boldsymbol{x} \frac{\partial}{\partial \boldsymbol{x}} + (\boldsymbol{x} - \boldsymbol{v} t) \frac{\partial}{\partial \boldsymbol{v}}, \\ & f^{(6)}_{\theta}(t, \boldsymbol{x}, \boldsymbol{v}) = f(\frac{t}{1 - \theta t}, \frac{\boldsymbol{x}}{1 - \theta t}, \boldsymbol{v} + \theta (\boldsymbol{x} - \boldsymbol{v} t)), \label{315} \end{aligned} \end{align}\tag{44}\] where \(\boldsymbol{v}\) is the particle velocity. According to Ref.[29], [30], the above generator is valid only when the system is under the power law potential \(U(r) \propto r^{-m}\) with \(m = 2\), see [29], [30] for more details. However, due to the stringent requirements of relativistic covariance and the on-shell condition for particles, no generator analogous to \(\hat{L}^{(6)}\) can be constructed in the relativistic regime.
In this work, we provide a novel and efficient approach to constructing the relativistic BKW solution of the Boltzmann equation. We first introduce a class of ansatz functions for the single-particle distribution function, and then demonstrate that within this specific functional space, only two forms yield exact solutions: the equilibrium distribution and the BKW-type solution. This analysis not only gives a streamlined derivation compared to conventional moment methods but also offers insight into the uniqueness of these solutions.
Furthermore, guided by physical intuition and the structure of known analytical solutions, particularly the NTFP, we have systematically derived the closed Lie algebra of invariant transformations admitted by the relativistic Boltzmann equation. To the best of our knowledge, this is the first time such a closed and comprehensive symmetry algebra has been constructed in the context of relativistic kinetic theory. The resulting algebra, spanned by generators \(L^{(0)}\) to \(L^{(5)}\), encompasses Poincaré symmetry (\(\hat{L}^{(0)}\) – \(\hat{L}^{(3)}\)), spacetime scaling (\(L^{(4)}\)), and a more general momentum-rescaling transformation (\(L^{(5)}\)) tied to the specific interactions of scaling cross-section.
There are several avenues for future work that can be pursued based on the current findings. First, we can explore whether there is a better approach to provide a first-principle determination of the Lie algebra, which might help us exhaustively enumerate all physically relevant invariant Lie algebras. Secondly, the properties of the linearized collision operator around the BKW solution are also worth investigating. How does it differ from the linearized collision operator around the equilibrium state? What distinct behaviors might emerge in the linear response around the BKW solution? As an exact analytical solution of the Boltzmann equation, the BKW solution holds the potential to offer insights into how to probe the linear response behavior of non-equilibrium states within kinetic theory. Additionally, as another self-similar solution of the relativistic Boltzmann equation, the connection between the BKW solution and the NTFP has yet to be explored. This presents an important direction for future research that could yield significant new insights.
: We express our gratitude for the valuable discussions with Qiuze Sun, Tianzhe Zhou, Xiaotian Ma, Haiyang Shao, Shuzhe Shi, Qi Chen, Wan Wu, Luyao Li, Baiting Tian, Puyuan Bai, Weiyao Ke, Yi Yin, and others. This work was financially supported by the National Natural Science Foundation of China under Grant Nos. 12035006 and 12505149.
We use \(f_n(\tau, p^0)\) to denote the single-particle distribution function given in Eq.@eq:222 , whose highest power of \(p^0\) is \(n\). Then we have \[\begin{align} \begin{aligned} & f_n(\tau, p_f^0) f_n(\tau, k_f^0) - f_n(\tau, p^0) f_n(\tau, k^0) \\ = & f_{n - 1}(\tau, p_f^0) f_{n - 1}(\tau, k_f^0) - f_{n - 1}(\tau, p^0) f_{n - 1}(\tau, k^0) \\ &+ A_n(\tau) \left[(p_f^0)^n f_{n - 1}(\tau, k_f^0) + (k_f^0)^n f_{n - 1}(\tau, p_f^0) - (p_f^0)^n f_{n - 1}(\tau, k_f^0) - (p_f^0)^n f_{n - 1}(\tau, k_f^0)\right] \\ & + A_n^2(\tau) \left[(p_f^0)^n (k_f^0)^n - (p^0)^n (k^0)^n\right] \\ = & f_{n - 1}(\tau, p^0 - \omega) f_{n - 1}(\tau, k^0 + \omega) - f_{n - 1}(\tau, p^0) f_{n - 1}(\tau, k^0) \\ & + A_n(\tau) \left[(p^0 - \omega)^n f_{n - 1}(\tau, k^0 + \omega) + (k^0 + \omega)^n f_{n - 1}(\tau, p^0 - \omega) \right.\\ &\left. - (p^0)^n f_{n - 1}(\tau, k^0) - (k^0)^n f_{n - 1}(\tau, p^0)\right] + A_n^2(\tau) \left[(p^0 - \omega)^n (k^0 + \omega)^n - (p^0)^n (k^0)^n\right], \label{A1} \end{aligned} \end{align}\tag{45}\] where the above expression will be utilized to calculate the collision kernel. When \(n = 2\), the trial function takes the following form \[\begin{align} \begin{aligned} f_2(\tau ,p^0) = e^{- \frac{p^0}{\alpha(\tau)}} \left(A_0(\tau) + A_1(\tau) p^0 + A_2(\tau) (p^0)^2\right), \label{A2} \end{aligned} \end{align}\tag{46}\]
and the collision term in Eq.@eq:221 becomes
\[\begin{align} \begin{aligned} C[f_2] = & C[f_1] + \frac{1}{(2 \pi)^2} e^{- \frac{|\boldsymbol{p}|}{\alpha(\tau)}} \\ & \times \left[ \int_{\frac{|\boldsymbol{q}| - \omega}{2}}^{+ \infty}\mathrm{d}|\boldsymbol{k}|\int_{-|\boldsymbol{q}|}^{+|\boldsymbol{q}|} \mathrm{d} \omega \int_0^{|\boldsymbol{p}|}\mathrm{d}|\boldsymbol{q}| e^{- \frac{|\boldsymbol{k}|}{\alpha(\tau)}} |\boldsymbol{k}| T(\tau, \boldsymbol{p}, \boldsymbol{k}, \omega) \right.\left. \left(1 - (\frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|}) (\frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|})\right) \right.\\ &\left. + \int_{\frac{|\boldsymbol{q}| - \omega}{2}}^{+ \infty}\mathrm{d}|\boldsymbol{k}| \int_{-|\boldsymbol{q}|}^{2|\boldsymbol{p}| - |\boldsymbol{q}|}\mathrm{d} \omega \int_{|\boldsymbol{p}|}^{+ \infty}\mathrm{d}|\boldsymbol{q}| e^{- \frac{|\boldsymbol{k}|}{\alpha(\tau)}} |\boldsymbol{k}| T(\tau, \boldsymbol{p}, \boldsymbol{k}, \omega) \right.\left. \left(1 - (\frac{\omega}{|\boldsymbol{q}|} - \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{p}||\boldsymbol{q}|}) (\frac{\omega}{|\boldsymbol{q}|} + \frac{\omega^2 - |\boldsymbol{q}|^2}{2|\boldsymbol{k}||\boldsymbol{q}|})\right)\right] \\ = & e^{- \frac{p^0}{\alpha(\tau)}} \left[R_1(\tau) p^0 + R_2(\tau) (p^0)^2 + R_3(\tau) (p^0)^3 + R_4(\tau) (p^0)^4 + R_5(\tau) (p^0)^5\right], \label{A3} \end{aligned} \end{align}\tag{47}\]
where \[\begin{align} \begin{aligned} T(\tau, \boldsymbol{p}, \boldsymbol{k}, \omega) \equiv & A_2(\tau) \left[(|\boldsymbol{p}| - \omega)^2 f_1(\tau, |\boldsymbol{k}| + \omega) + (|\boldsymbol{k}| + \omega)^2 f_1(\tau, |\boldsymbol{p}| - \omega) \right.\left. - |\boldsymbol{p}|^2 f_1(\tau, |\boldsymbol{k}|) - |\boldsymbol{k}|^2 f_1(\tau, |\boldsymbol{p}|)\right] \\ & + A_2^2(\tau) \left[(|\boldsymbol{p}| - \omega)^2 (|\boldsymbol{k}| + \omega)^2 - |\boldsymbol{p}|^2 |\boldsymbol{k}|^2\right], \label{A4} \end{aligned} \end{align}\tag{48}\]
and the following shorthand notations are introduced: \[\begin{align} \begin{aligned} & R_1(\tau) \equiv \frac{2 \alpha^5(\tau)}{\pi^2} \left[A_1^2(\tau) - 2 A_0(\tau) A_2(\tau) + 5 A_1(\tau) A_2(\tau) \alpha(\tau) + 6 A_2^2(\tau) \alpha^2(\tau)\right], \\ & R_2(\tau) \equiv \frac{2 \alpha^4(\tau)}{3 \pi^2} \left[-2 A_1^2(\tau) + 4 A_0(\tau) A_2(\tau) - 5 A_1(\tau) A_2(\tau) \alpha(\tau) + 16 A_2^2(\tau) \alpha^2(\tau)\right], \\ & R_3(\tau) \equiv \frac{\alpha^3(\tau)}{6 \pi^2} \left[A_1^2(\tau) - 2 A_0(\tau) A_2(\tau) - 5 A_1(\tau) A_2(\tau) \alpha(\tau) - 44 A_2^2(\tau) \alpha^2(\tau)\right], \\ & R_4(\tau) \equiv \frac{A_2(\tau) \alpha^3(\tau)}{30 \pi^2} \left[5 A_1(\tau) + 16 A_2(\tau) \alpha(\tau)\right], \\ & R_5(\tau) \equiv \frac{A_2^2(\tau) \alpha^2(\tau)}{30 \pi^2}. \label{A5} \end{aligned} \end{align}\tag{49}\] By substituting Eq.@eq:A2 into the left-hand side of Eq.@eq:220 , combining the result with Eq.@eq:A3 , and matching the coefficients of identical powers of \(p^0\) on both sides, we obtain a new set of ordinary differential equations that closely resembles Eq.@eq:228 : \[\begin{align} \begin{aligned} A_0'(\tau) = \frac{2 A_1^2(\tau) \alpha^5(\tau)}{\pi^2} - \frac{4 A_0(\tau) A_2(\tau) \alpha^5(\tau)}{\pi^2} + \frac{10 A_1(\tau) A_2(\tau) \alpha^6(\tau)}{\pi^2} + \frac{12 A_2^2(\tau) \alpha^7(\tau)}{\pi^2}, \label{A6} \end{aligned} \end{align}\tag{50}\]
\[\begin{align} \begin{aligned} A_1'(\tau) + \frac{A_0(\tau) \alpha'(\tau)}{\alpha^2(\tau)} & = - \frac{4 A_1^2(\tau) \alpha^4(\tau)}{3 \pi^2} + \frac{8 A_0(\tau) A_2(\tau) \alpha^4(\tau)}{3 \pi^2} - \frac{10 A_1(\tau) A_2(\tau) \alpha^5(\tau)}{3 \pi^2} + \frac{32 A_2^2(\tau) \alpha^6(\tau)}{3 \pi^2}, \label{A7} \end{aligned} \end{align}\tag{51}\]
\[\begin{align} \begin{aligned} A_2'(\tau) + \frac{A_1(\tau) \alpha'(\tau)}{\alpha^2(\tau)} & = \frac{A_1^2(\tau) \alpha^3(\tau)}{6 \pi^2} - \frac{A_0(\tau) A_2(\tau) \alpha^3(\tau)}{3 \pi^2} - \frac{5 A_1(\tau) A_2(\tau) \alpha^4(\tau)}{6 \pi^2} - \frac{22 A_2^2(\tau) \alpha^5(\tau)}{3 \pi^2}, \label{A8} \end{aligned} \end{align}\tag{52}\]
\[\begin{align} \begin{aligned} \frac{A_2(\tau) \alpha'(\tau)}{\alpha^2(\tau)} = \frac{A_1(\tau) A_2(\tau) \alpha^3(\tau)}{6 \pi^2} + \frac{8 A_2^2(\tau) \alpha^4(\tau)}{15 \pi^2}, \label{A9} \end{aligned} \end{align}\tag{53}\]
\[\begin{align} \begin{aligned} \frac{A_2^2(\tau) \alpha^3(\tau)}{30 \pi^2} = 0. \label{A10} \end{aligned} \end{align}\tag{54}\] From Eq.@eq:A10 , we directly obtain a non-trivial constraint \[\begin{align} \begin{aligned} A_2(\tau) = 0, \label{A11} \end{aligned} \end{align}\tag{55}\] and find that Eq.@eq:A2 reduces to Eq.@eq:223 , Eqs.@eq:A6 – 52 reduce to Eq.@eq:228 , Eqs.@eq:A9 and 54 vanish identically. In other words, \(f_n(\tau, p^0)\) for \(n = 2\) is not a solution of the Boltzmann equation, unless it degenerates into the form corresponding to \(n = 1\).
When \(n = 3\), repeating the above steps yields another set of ordinary differential equations \[\begin{align} \begin{aligned} A_0'(\tau) = & \frac{2 A_1^2(\tau) \alpha^5(\tau)}{\pi^2} - \frac{4 A_0(\tau) A_2(\tau) \alpha^5(\tau)}{\pi^2} + \frac{10 A_1(\tau) A_2(\tau) \alpha^6(\tau)}{\pi^2} \\ & + \frac{12 A_2^2(\tau) \alpha^7(\tau)}{\pi^2} - \frac{30 A_0(\tau) A_3(\tau) \alpha^6(\tau)}{\pi^2} + \frac{36 A_1(\tau) A_3(\tau) \alpha^7(\tau)}{\pi^2} \\ & + \frac{84 A_2(\tau) A_3(\tau) \alpha^8(\tau)}{\pi^2} + \frac{144 A_3^2(\tau) \alpha^9(\tau)}{\pi^2}, \label{A12} \end{aligned} \end{align}\tag{56}\]
\[\begin{align} \begin{aligned} A_1'(\tau) + \frac{A_0(\tau) \alpha'(\tau)}{\alpha^2(\tau)} = & - \frac{4 A_1^2(\tau) \alpha^4(\tau)}{3 \pi^2} + \frac{8 A_0(\tau) A_2(\tau) \alpha^4(\tau)}{3 \pi^2} - \frac{10 A_1(\tau) A_2(\tau) \alpha^5(\tau)}{3 \pi^2} \\ & + \frac{32 A_2^2(\tau) \alpha^6(\tau)}{3 \pi^2} + \frac{10 A_0(\tau) A_3(\tau) \alpha^5(\tau)}{\pi^2} - \frac{28 A_1(\tau) A_3(\tau) \alpha^6(\tau)}{\pi^2} \\ & + \frac{76 A_2(\tau) A_3(\tau) \alpha^7(\tau)}{\pi^2} + \frac{132 A_3^2(\tau) \alpha^8(\tau)}{\pi^2}, \label{A13} \end{aligned} \end{align}\tag{57}\]
\[\begin{align} \begin{aligned} A_2'(\tau) + \frac{A_1(\tau) \alpha'(\tau)}{\alpha^2(\tau)} = & \frac{A_1^2(\tau) \alpha^3(\tau)}{6 \pi^2} - \frac{A_0(\tau) A_2(\tau) \alpha^3(\tau)}{3 \pi^2} - \frac{5 A_1(\tau) A_2(\tau) \alpha^4(\tau)}{6 \pi^2} \\ & - \frac{22 A_2^2(\tau) \alpha^5(\tau)}{3 \pi^2} + \frac{5 A_0(\tau) A_3(\tau) \alpha^4(\tau)}{2 \pi^2} + \frac{8 A_0(\tau) A_3(\tau) \alpha^5(\tau)}{\pi^2} \\ & - \frac{26 A_2(\tau) A_3(\tau) \alpha^6(\tau)}{\pi^2} + \frac{60 A_3^2(\tau) \alpha^7(\tau)}{\pi^2}, \label{A14} \end{aligned} \end{align}\tag{58}\]
\[\begin{align} \begin{aligned} A_3'(\tau) + \frac{A_2(\tau) \alpha'(\tau)}{\alpha^2(\tau)} = & \frac{A_1(\tau) A_2(\tau) \alpha^3(\tau)}{6 \pi^2} + \frac{8 A_2^2(\tau) \alpha^4(\tau)}{15 \pi^2} - \frac{A_0(\tau) A_3(\tau) \alpha^3(\tau)}{2 \pi^2} \\ & - \frac{7 A_1(\tau) A_3(\tau) \alpha^4(\tau)}{5 \pi^2} - \frac{26 A_2(\tau) A_3(\tau) \alpha^5(\tau)}{5 \pi^2} - \frac{42 A_3^2(\tau) \alpha^6(\tau)}{\pi^2}, \label{A15} \end{aligned} \end{align}\tag{59}\]
\[\begin{align} \begin{aligned} \frac{A_3(\tau) \alpha'(\tau)}{\alpha^2(\tau)} = \frac{A_2^2(\tau) \alpha^3(\tau)}{30 \pi^2} + \frac{A_1(\tau) A_3(\tau) \alpha^3(\tau)}{10 \pi^2} + \frac{19 A_2(\tau) A_3(\tau) \alpha^4(\tau)}{30 \pi^2} + \frac{2 A_3^2(\tau) \alpha^5(\tau)}{\pi^2}, \label{A16} \end{aligned} \end{align}\tag{60}\]
\[\begin{align} \begin{aligned} \frac{A_2(\tau) A_3(\tau) \alpha^3(\tau)}{30 \pi^2} + \frac{11 A_3^2(\tau) \alpha^4(\tau)}{70 \pi^2} = 0, \label{A17} \end{aligned} \end{align}\tag{61}\]
\[\begin{align} \begin{aligned} \frac{A_3^2(\tau) \alpha^3(\tau)}{140 \pi^2} = 0. \label{A18} \end{aligned} \end{align}\tag{62}\] Once again from Eq.@eq:A18 , we arrive at \[\begin{align} \begin{aligned} A_3(\tau) = 0. \label{A19} \end{aligned} \end{align}\tag{63}\] Therefore \(f_n(\tau,p^0)\) with \(n=3\) still degenerates into the form of \(n=1\) similarly.
Finally, we note that for the case of \(n \geq 2\), the term \((p^0 - \omega)^n \omega^n\) coming from \((p^0 - \omega)^n (k^0 + \omega)^n\) in Eq.@eq:A1 gives the highest power of \(p^0\) after integration, whose coefficients have the form of \(A_n^2(\tau) \alpha^3(\tau) /(c_n \pi^2)\) where \(c_n\) is a constant. However, on the left-hand side of Eq.@eq:21 , this highest-power term does not appear, resulting in \[\begin{align} \begin{aligned} \frac{A_n^2(\tau) \alpha^3(\tau)}{c_n \pi^2} = 0, \label{A20} \end{aligned} \end{align}\tag{64}\] and \[\begin{align} \begin{aligned} A_n(\tau) = 0. \label{A21} \end{aligned} \end{align}\tag{65}\] Consequently, both \(f_n(\tau, p^0)\) and its associated system of ordinary differential equations reduce to the same form as previously analyzed, ultimately degenerating into Eq.@eq:223 and Eq.@eq:228 , respectively.