Gradient Expansion Formalism (GEF)

• Gradient Expansion Formalism (GEF) is a systematic method that reorganizes complex dynamics by expanding in small spacetime gradients rather than using a full mode-by-mode analysis.
• Its applications span cosmology, hydrodynamics, spintronics, and quantum transport, where it simplifies infinite hierarchies into manageable, order-by-order truncations.
• Practical implementations of GEF, such as in inflationary magnetogenesis and leptogenesis, demonstrate controlled approximations with error margins reduced to a few percent.

Gradient Expansion Formalism (GEF) denotes a family of techniques in which a dynamical problem is reorganized by the number of spacetime or spatial gradients, rather than solved directly in its full mode-by-mode form. In the literature, the name appears in several technically distinct settings: inflationary gauge-field production and magnetogenesis (Sobol et al., 2020), axion inflation with backreaction and Schwinger damping (Gorbar et al., 2021), superhorizon cosmological perturbation theory beyond the separate-universe limit (Naruko et al., 2012), closed-time-path kinetic theory for leptogenesis (Garny et al., 2010), causal relativistic hydrodynamics (Lahiri, 2019), gauge-covariant spin-transport theory (Shitade, 2017), and other gradient-based expansions in density-functional theory and gradient flow (Benites et al., 30 Jan 2026). This suggests that GEF is best understood not as a single universal formalism, but as a recurring methodological pattern: identify a small-gradient regime, construct an ordered hierarchy of operators or correlators, and truncate that hierarchy in a controlled way.

1. Common organizing principle

Across applications, the expansion variable is always tied to a separation of scales. In cosmology on superhorizon scales one expands in a long-wavelength parameter such as $\varepsilon \equiv 1/(HL)\ll1$ or $\sigma \equiv k/(aH)\ll1$; in relativistic hydrodynamics the small parameter is the Knudsen number $Kn\sim \ell_m/L\ll1$; in closed-time-path kinetic theory gradients enter through the Moyal-type operator $\diamondsuit$ and correspond to powers of $\hbar$; in inflationary gauge-field production the hierarchy is often indexed by the number of spatial curls acting on bilinear correlators (Naruko et al., 2012).

Domain	Ordering variable	Reformulation
Inflationary gauge fields	number of curls or gradients; superhorizon amplification	infinite chain of bilinear correlators
Superhorizon perturbations	$\varepsilon=1/(HL)\ll1$ or $\sigma=k/(aH)\ll1$	ADM system expanded through $O(\varepsilon^2)$
Leptogenesis	gradients $\partial_X$ generated by $\diamondsuit$	Wigner-space KB equations
Relativistic hydrodynamics	$\sigma \equiv k/(aH)\ll1$0	constitutive relations by gradient order
Spin torques	spacetime gradients and field strengths	gauge-covariant Wigner/Dyson expansion

The practical benefit is similar in all cases. The original problem usually involves either a large set of Fourier modes, nonlocal two-point functions, or an enormous basis of tensor structures. GEF replaces that description by a hierarchy of local or spatially homogeneous quantities with a clear order-by-order meaning. The main technical challenge is then closure: the exact hierarchy is typically infinite, and one must specify an asymptotic relation, a relaxation ansatz, or a finite-order truncation.

2. Inflationary magnetogenesis in the kinetic-coupling model

In the kinetic-coupling model of inflationary magnetogenesis, the basic action is

$\sigma \equiv k/(aH)\ll1$1

with a spatially flat FLRW background, homogeneous inflaton $\sigma \equiv k/(aH)\ll1$2, and time-dependent kinetic coupling $\sigma \equiv k/(aH)\ll1$3. The physical electric and magnetic fields are

$\sigma \equiv k/(aH)\ll1$4

The motivation for GEF in this setting is that only long-wavelength superhorizon modes undergo amplification, while the standard mode-by-mode treatment breaks down once backreaction or Schwinger pair production becomes important (Sobol et al., 2020).

The formalism introduces three infinite families of bilinear vacuum expectation values,

$\sigma \equiv k/(aH)\ll1$5

$\sigma \equiv k/(aH)\ll1$6

Here $\sigma \equiv k/(aH)\ll1$7 and $\sigma \equiv k/(aH)\ll1$8 generalize electric and magnetic energy-density bilinears by inserting $\sigma \equiv k/(aH)\ll1$9 curls, while $Kn\sim \ell_m/L\ll1$0 is an $Kn\sim \ell_m/L\ll1$1–$Kn\sim \ell_m/L\ll1$2 mixed correlator. At $Kn\sim \ell_m/L\ll1$3, $Kn\sim \ell_m/L\ll1$4 is twice the total electromagnetic energy density, so that

$Kn\sim \ell_m/L\ll1$5

From Maxwell’s equations one obtains an infinite chain of coupled ODEs,

$Kn\sim \ell_m/L\ll1$6

$Kn\sim \ell_m/L\ll1$7

$Kn\sim \ell_m/L\ll1$8

The $Kn\sim \ell_m/L\ll1$9 terms describe mode enhancement from the kinetic coupling, and the appearance of $\diamondsuit$0 moments in the $\diamondsuit$1th equation makes the system infinite. Physically, each $\diamondsuit$2 receives contributions from all superhorizon modes, while new modes continuously enter at the horizon-crossing scale $\diamondsuit$3.

Closure is achieved by exploiting the large-$\diamondsuit$4 asymptotics

$\diamondsuit$5

up to $\diamondsuit$6 corrections. In practice one increases the maximal order $\diamondsuit$7 until observables converge at the percent level. For Starobinsky inflation with both Ratra power-law and nonmonotonic cosh-type couplings, the reported benchmarks are that $\diamondsuit$8–8 yields electric and magnetic energy densities within $\diamondsuit$9 of exact Fourier-mode integration, while $\hbar$0 reduces residual errors to a few percent across increasing, decreasing, and nonmonotonic couplings, both with and without strong backreaction. Convergence is faster for the dominant component: electric for decreasing $\hbar$1, magnetic for increasing $\hbar$2 (Sobol et al., 2020).

The same work also states the principal limitations. The truncation relies on approximate mode functions near horizon crossing and on de Sitter or small-slow-roll approximations; errors of a few percent remain. The method cannot be straightforwardly extended beyond inflation, where $\hbar$3 may freeze and spatial inhomogeneity of $\hbar$4 becomes important. Strong Schwinger pair production would require additional damping terms not yet included.

3. Axion inflation, Schwinger damping, gravitational waves, and software implementations

A closely related but distinct GEF was developed for gauge-field production during axion inflation. Instead of the kinetic coupling $\hbar$5, the sourced Maxwell system contains an axial term and, in fermionic settings, an Ohmic current $\hbar$6. The bilinear hierarchy is again written in terms of

$\hbar$7

with $\hbar$8, $\hbar$9, and $\varepsilon=1/(HL)\ll1$0 as the helicity density (Gorbar et al., 2021).

In this version, backreaction and Schwinger pair production enter directly into the hierarchy. The conductivity produces friction terms proportional to $\varepsilon=1/(HL)\ll1$1 in the $\varepsilon=1/(HL)\ll1$2 and $\varepsilon=1/(HL)\ll1$3 equations, and also appears as a sink term $\varepsilon=1/(HL)\ll1$4 in the energy balance. A new nonlocal-in-time quantity,

$\varepsilon=1/(HL)\ll1$5

multiplies the short-scale Bunch–Davies vacuum. Because $\varepsilon=1/(HL)\ll1$6 depends on the full history of $\varepsilon=1/(HL)\ll1$7, the description becomes inherently nonlocal in time. The same framework also yields Whittaker-function solutions for the mode functions, which are then used to construct boundary terms for the bilinear chain (Gorbar et al., 2021).

For pure axion inflation, a later implementation normalizes the correlators by a time-dependent cutoff momentum

$\varepsilon=1/(HL)\ll1$8

chosen to enclose all tachyonically amplified modes. The resulting ODE tower is solved together with the homogeneous inflaton and Friedmann equations, under the explicit assumption of vanishing axion gradients, $\varepsilon=1/(HL)\ll1$9. Truncation is imposed at $\sigma=k/(aH)\ll1$0, with a closure parameter $\sigma=k/(aH)\ll1$1–100 chosen so that the hierarchy converges at the percent level. The practical accuracy criterion is

$\sigma=k/(aH)\ll1$2

with a target $\sigma=k/(aH)\ll1$3 and reported performance of $\sigma=k/(aH)\ll1$4 in practice; if needed, the system is reinitialized from the mode-by-mode solution (Eckardstein et al., 1 Aug 2025).

This formulation was then used for gravitational-wave benchmarks. In pure axion inflation, a comprehensive parameter scan near the onset of strong backreaction found that GW signals within the reach of future GW interferometers can only be realized in parameter regions that also lead to strong backreaction and that are in conflict with the upper limit on $\sigma=k/(aH)\ll1$5 (Eckardstein et al., 1 Aug 2025). In fermionic axion inflation, by contrast, Schwinger pair creation of charged fermions damps gauge-field production, attenuates the GW signal, and can bring the signal into the sensitivity reach of LISA and ET without violating the upper limit on $\sigma=k/(aH)\ll1$6. The same study identifies a new backreaction regime in which the axion velocity and energy-density components exhibit damped oscillations around the slow-roll trajectory, moderated by fermion production (Eckardstein et al., 29 Sep 2025).

The numerical infrastructure for this line of work is the GEFF package, “The Gradient Expansion Formalism Factory,” a Python package that provides ready-to-use model files for pure axion inflation and fermionic axion inflation, methods to solve GEF equations, an integrated error estimator, a self-correction algorithm, and tools to study gauge-field-induced primordial gravitational waves (Eckardstein, 14 Oct 2025).

4. Superhorizon cosmology and extensions beyond the standard $\sigma=k/(aH)\ll1$7 limit

In nonlinear cosmological perturbation theory, GEF is the long-wavelength expansion of the Einstein–matter system in powers of

$\sigma=k/(aH)\ll1$8

with each spatial derivative counting as $\sigma=k/(aH)\ll1$9. In ADM variables one expands all fields as

$O(\varepsilon^2)$0

so that the separate-universe approximation is the leading term, while $O(\varepsilon^2)$1 provides the first nontrivial superhorizon correction. For multi-component scalar fields with general kinetic terms and potentials, the formalism yields the Hamiltonian and momentum constraints, dynamical equations for $O(\varepsilon^2)$2 and $O(\varepsilon^2)$3, and fully nonlinear gauge-transformation rules valid through $O(\varepsilon^2)$4 (Naruko et al., 2012).

A central structural statement is that the usual $O(\varepsilon^2)$5 formalism is precisely the $O(\varepsilon^2)$6 part of the gradient expansion. At leading order, each Hubble patch evolves like an FLRW universe with local proper time, and the curvature perturbation is identified with a local e-fold difference. At next order, spatial gradients of the leading solution source $O(\varepsilon^2)$7, $O(\varepsilon^2)$8, and matter corrections. In the multi-field construction, the practical procedure is: solve the leading FLRW-like system in uniform-$O(\varepsilon^2)$9 gauge, solve the next-to-leading equations with spatial gradients as sources, transform nonlinearly to uniform-$\partial_X$0 gauge, and extract the physical curvature perturbation $\partial_X$1 (Naruko et al., 2012).

The same general framework has been extended to single-field $\partial_X$2-inflation, multi-field models, Galileon or G-inflation, and scalar-tensor theories. In the review formulation, the nonlinear curvature perturbation obeys a generalized Mukhanov–Sasaki-type equation of the schematic form

$\partial_X$3

where $\partial_X$4 is the nonlinear spatial-curvature term. The formalism is explicitly designed to treat both scalar and tensor superhorizon modes and to go beyond the pure separate-universe result (Takamizu, 2018).

For shift-symmetric Galileon-type actions, the gradient expansion was worked out in synchronous gauge up to second order without imposing extra conditions at first order. In that setting, the first-order curvature perturbation

$\partial_X$5

is exactly conserved in time at $\partial_X$6, while the second-order solution contains explicit decaying modes proportional to $\partial_X$7 and $\partial_X$8 during late quasi-de Sitter evolution (Frusciante et al., 2013).

A recent development reformulates finite-$\partial_X$9 corrections to $\diamondsuit$0 by adding an effective source to the background Klein–Gordon equation,

$\diamondsuit$1

so that

$\diamondsuit$2

This construction is presented as a way to recover the full Mukhanov–Sasaki equation once one linearizes in $\diamondsuit$3, to include $\diamondsuit$4, $\diamondsuit$5, and higher corrections systematically, and to compute quantities such as the equilateral non-Gaussianity parameter

$\diamondsuit$6

The paper states that this GEF–$\diamondsuit$7 method reproduces the full Mukhanov–Sasaki solution for both the power spectrum and $\diamondsuit$8 even when matching is performed at horizon crossing, where standard $\diamondsuit$9 fails (Ahmadi et al., 31 Jan 2026).

5. Nonequilibrium kinetics, hydrodynamics, and transport theory

In leptogenesis, the relevant GEF is the gradient expansion of the Kadanoff–Baym equations in Wigner space. With center and relative coordinates

$\sigma \equiv k/(aH)\ll1$00

the Wigner transform turns fast relative-coordinate oscillations into momentum dependence and slow center-coordinate variation into spacetime gradients. The Moyal-type operator

$\sigma \equiv k/(aH)\ll1$01

organizes the expansion. At zeroth order, $\sigma \equiv k/(aH)\ll1$02, and one recovers the standard Boltzmann gain-minus-loss terms. At first order, one obtains additional Poisson-bracket structures and a new CP-violating source $\sigma \equiv k/(aH)\ll1$03 proportional to the Hubble rate $\sigma \equiv k/(aH)\ll1$04, nonzero even when all species are in local thermal equilibrium, but vanishing in static equilibrium in accordance with the Sakharov conditions. The paper concludes that this source is typically tiny in standard thermal leptogenesis and becomes relevant only in the limit of ultra-strong washout (Garny et al., 2010).

In relativistic hydrodynamics, GEF is a systematic derivative expansion of the conserved currents around local thermodynamic equilibrium. Zeroth order gives perfect-fluid hydrodynamics, first order gives relativistic Navier–Stokes, and second order adds all possible scalar, vector, and tensor structures consistent with the frame and symmetries. In the Eckart frame, the dissipative fluxes are expanded as

$\sigma \equiv k/(aH)\ll1$05

$\sigma \equiv k/(aH)\ll1$06

where the second-order basis contains curvature invariants such as $\sigma \equiv k/(aH)\ll1$07, $\sigma \equiv k/(aH)\ll1$08, $\sigma \equiv k/(aH)\ll1$09, and $\sigma \equiv k/(aH)\ll1$10. To restore hyperbolicity and causality, one upgrades these constitutive relations to relaxation-type equations with times $\sigma \equiv k/(aH)\ll1$11, $\sigma \equiv k/(aH)\ll1$12, and $\sigma \equiv k/(aH)\ll1$13, in the spirit of Müller–Israel–Stewart theory (Lahiri, 2019).

A related structural refinement proves two theorem-level simplifications for the gradient expansion of relativistic hydrodynamics. For non-conformal fluids, the ordering of transverse derivatives is irrelevant up to curvature tensors and lower-order gradients already present in the basis. For conformal fluids, the longitudinal projection of the Weyl-covariant derivative acting on hydrodynamic variables can be eliminated, although this statement does not apply to curvature tensors. These theorems sharply reduce the bookkeeping of independent transport structures at higher order (Diles, 2020).

In spintronics, GEF appears in a gauge-covariant Keldysh–Wigner formalism. One defines a Wilson-line-dressed Green’s function and expands the Moyal-star product in covariant derivatives $\sigma \equiv k/(aH)\ll1$14 and field strengths $\sigma \equiv k/(aH)\ll1$15. The Dyson equation is then solved order by order in $\sigma \equiv k/(aH)\ll1$16. The spin density

$\sigma \equiv k/(aH)\ll1$17

yields the torque

$\sigma \equiv k/(aH)\ll1$18

In the reported ordering, spin renormalization $\sigma \equiv k/(aH)\ll1$19 and Gilbert damping $\sigma \equiv k/(aH)\ll1$20 arise from the first-order term $\sigma \equiv k/(aH)\ll1$21, while the spin-transfer torque and the nonadiabatic $\sigma \equiv k/(aH)\ll1$22-term arise from the mixed second-order term $\sigma \equiv k/(aH)\ll1$23 (Shitade, 2017).

6. Other realizations, controversies, and scope

In density-functional theory for the inhomogeneous electron gas, the gradient-expansion approximation derives the $\sigma \equiv k/(aH)\ll1$24 correction to the exchange-correlation functional from the $\sigma \equiv k/(aH)\ll1$25 term in the proper-polarization function. A key conclusion is that the separate exchange and correlation contributions to the coefficient $\sigma \equiv k/(aH)\ll1$26 are regulator-scheme dependent even when the regulator is set to zero at the end, whereas the total

$\sigma \equiv k/(aH)\ll1$27

is regulator-scheme independent and unique. The paper therefore concludes that it is incorrect to separate exchange and correlation when constructing a generalized-gradient-approximation contribution to the density functional (Benites et al., 30 Jan 2026). This is a concrete controversy internal to one branch of gradient-expansion literature: the expansion itself is not disputed, but the decomposition of its coefficient is.

A different recent use of the term appears in the study of gradient flow in large learning problems. There, the loss under gradient flow is expanded as a formal power series in time, with coefficients encoded by diagrams akin to Feynman diagrams. For CP tensor decomposition, the large-size limit exposes distinct extreme lazy and rich regimes, including free evolution, NTK, and under- and over-parameterized mean-field phases; the resulting generating functions can, in many cases, be reduced to first-order PDEs solvable by the method of characteristics (Yarotsky et al., 4 Feb 2026). Although technically far removed from cosmology or hydrodynamics, this usage preserves the same core logic: a complicated nonlinear evolution is reorganized into a hierarchy indexed by a controlled expansion.

The inflationary gauge-field literature also illustrates the scope for extension. The magnetogenesis hierarchy has been proposed for axial coupling $\sigma \equiv k/(aH)\ll1$28, non-Abelian gauge fields, higher-spin systems, and gravitational-wave generation by the same hierarchy of curls (Sobol et al., 2020). A later extension to a massive vector field introduces additional bilinears for the longitudinal sector, couples them to transverse correlators and the inflaton, and finds that pure mass coupling enhances only longitudinal modes, while adding kinetic coupling changes whether transverse or longitudinal polarizations dominate the energy density. For increasing coupling functions, the longitudinal component can rapidly propel the system into the strong backreaction regime (Lysenko et al., 29 Sep 2025).

The main limitations are correspondingly context-dependent. Inflationary GEFs often assume quasi-de Sitter evolution, adiabatically varying couplings, and a homogeneous inflaton; they generally cease to be straightforward beyond inflation, in preheating, or once inflaton gradients become large (Sobol et al., 2020). In axion-inflation benchmarks, the homogeneous-axion approximation captures the onset of strong backreaction but may deviate after large $\sigma \equiv k/(aH)\ll1$29 develops (Eckardstein et al., 1 Aug 2025). This suggests a general rule that holds across the literature: GEF is most powerful when there is a clean hierarchy between slow macroscopic evolution and fast microscopic structure, and when the truncation can be calibrated against either exact mode integration, lattice simulation, or an independent microscopic theory.