6. Reheating

Reading time: ~55 minutes | Difficulty: Advanced

Prerequisites: Inflaton dynamics, slow-roll inflation, quantum field theory basics

Reheating is the critical bridge between the inflationary epoch and the hot Big Bang. After inflation ends, the universe is cold, empty, and dominated by a coherently oscillating inflaton field. Reheating converts the inflaton's energy into a thermal bath of Standard Model particles, setting the initial conditions for Big Bang nucleosynthesis, baryogenesis, and all subsequent cosmological evolution.

The Reheating Problem

During inflation, the inflaton field $\phi$ slowly rolls down its potential, driving exponential expansion that dilutes all pre-existing matter and radiation to negligible densities. When slow-roll conditions are violated ($\epsilon \sim 1$), inflation ends, but the universe is far from the hot, radiation-dominated state required by Big Bang nucleosynthesis. The inflaton must decay into lighter particles, which then thermalize to produce the hot Big Bang plasma. This process—reheating—determines the highest temperature the universe reaches after inflation, with profound consequences for particle physics and cosmology.

The physics of reheating proceeds in three stages: (1) the inflaton oscillates coherently around its potential minimum, behaving as a fluid with an effective equation of state; (2) the inflaton transfers energy to other fields via perturbative decay or non-perturbative parametric resonance (preheating); and (3) the produced particles scatter and thermalize, establishing thermal equilibrium at the reheating temperature $T_{\text{rh}}$.

Key question: The reheating temperature $T_{\text{rh}}$ must satisfy $T_{\text{rh}} \gtrsim 4$ MeV for successful nucleosynthesis, but must be low enough to avoid overproduction of gravitinos and other relics. This tension constrains the inflaton's couplings to the Standard Model.

Derivation 1: Inflaton Oscillations

After slow-roll ends, the inflaton $\phi$ oscillates around the minimum of its potential. We start from the Klein-Gordon equation in an expanding background:

$$\ddot{\phi} + 3H\dot{\phi} + V'(\phi) = 0$$

where $H = \dot{a}/a$ is the Hubble parameter and $V'(\phi) = dV/d\phi$. The Friedmann equation gives $3M_P^2 H^2 = \frac{1}{2}\dot{\phi}^2 + V(\phi)$.

Quadratic Potential: $V = \frac{1}{2}m^2\phi^2$

For a massive inflaton with $V(\phi) = \frac{1}{2}m^2\phi^2$, the equation of motion becomes:

$$\ddot{\phi} + 3H\dot{\phi} + m^2\phi = 0$$

When $m \gg H$ (the oscillation frequency far exceeds the expansion rate), we employ the WKB ansatz:

$$\phi(t) = \Phi(t)\sin(mt)$$

where $\Phi(t)$ is a slowly-varying amplitude with $\dot{\Phi}/\Phi \ll m$. Substituting into the Klein-Gordon equation and averaging over one oscillation period $\Delta t = 2\pi/m$:

$$\ddot{\phi} = \ddot{\Phi}\sin(mt) + 2m\dot{\Phi}\cos(mt) - m^2\Phi\sin(mt)$$

$$3H\dot{\phi} = 3H\dot{\Phi}\sin(mt) + 3Hm\Phi\cos(mt)$$

The $m^2\phi$ term cancels the $-m^2\Phi\sin(mt)$ piece. Collecting the slowly-varying terms and averaging (the $\sin\cos$ cross-terms vanish over a full period):

$$\dot{\Phi} + \frac{3}{2}H\Phi = 0$$

This is the equation for the amplitude decay. Since the energy density averages to $\langle\rho\rangle = \frac{1}{2}m^2\Phi^2$, we find:

$$\dot{\Phi} + \frac{3}{2}H\Phi = 0 \quad \Longrightarrow \quad \Phi \propto a^{-3/2}$$

Therefore $\langle\rho\rangle = \frac{1}{2}m^2\Phi^2 \propto a^{-3}$, which is precisely the scaling of pressureless matter. This makes physical sense: a coherently oscillating scalar field with a quadratic potential is equivalent to a collection of zero-momentum particles (a condensate of non-relativistic inflaton quanta).

Effective Equation of State

The time-averaged pressure is $\langle p \rangle = \langle\frac{1}{2}\dot{\phi}^2 - V\rangle$. For the quadratic potential, the virial theorem gives $\langle\dot{\phi}^2/2\rangle = \langle V\rangle$, so $\langle p\rangle = 0$ and $\langle w\rangle = 0$.

Quartic Potential: $V = \frac{1}{4}\lambda\phi^4$

For a self-interacting inflaton with $V(\phi) = \frac{1}{4}\lambda\phi^4$, the oscillation frequency depends on the amplitude, making the analysis more subtle. The virial theorem for a $\phi^4$ potential gives:

$$\langle\dot{\phi}^2/2\rangle = 2\langle V\rangle \quad \Longrightarrow \quad \langle\rho\rangle = \langle\dot{\phi}^2/2\rangle + \langle V\rangle = 3\langle V\rangle$$

The time-averaged pressure is:

$$\langle p\rangle = \langle\dot{\phi}^2/2 - V\rangle = 2\langle V\rangle - \langle V\rangle = \langle V\rangle = \frac{1}{3}\langle\rho\rangle$$

So $\langle w\rangle = 1/3$, which means radiation-like behavior: $\langle\rho\rangle \propto a^{-4}$. A coherently oscillating $\phi^4$ condensate mimics a relativistic gas.

General Monomial: $V \propto \phi^n$

For a general power-law potential $V(\phi) = \lambda_n |\phi|^n / n$, the virial theorem gives the time-averaged kinetic and potential energies:

$$\left\langle\frac{1}{2}\dot{\phi}^2\right\rangle = \frac{n}{2}\langle V\rangle$$

The effective equation of state follows:

$$\langle w\rangle = \frac{\langle p\rangle}{\langle\rho\rangle} = \frac{\langle\dot{\phi}^2/2\rangle - \langle V\rangle}{\langle\dot{\phi}^2/2\rangle + \langle V\rangle} = \frac{\frac{n}{2} - 1}{\frac{n}{2} + 1} = \frac{n - 2}{n + 2}$$

General Result

For $V \propto \phi^n$ near the minimum, the time-averaged equation of state is:

$$\boxed{\langle w \rangle = \frac{n-2}{n+2}}$$

Check: $n=2$ gives $w=0$ (matter), $n=4$ gives $w=1/3$ (radiation), $n \to \infty$ gives $w \to 1$ (stiff fluid). The energy density scales as $\langle\rho\rangle \propto a^{-6n/(n+2)}$.

Derivation 2: Perturbative Decay

The simplest model of reheating assumes the inflaton couples to a lighter scalar field $\chi$ via the interaction Lagrangian:

$$\mathcal{L}_{\text{int}} = -\frac{1}{2}g^2\phi\chi^2$$

This trilinear coupling allows the process $\phi \to \chi\chi$: one inflaton quantum decays into two $\chi$ particles.

Decay Rate from Fermi's Golden Rule

The decay rate is computed using Fermi's golden rule. The matrix element for $\phi \to \chi\chi$ from the vertex $g^2\phi\chi^2$ is simply $|\mathcal{M}|^2 = g^4$. The two-body decay rate in the rest frame of $\phi$ is:

$$\Gamma_\phi = \frac{|\mathcal{M}|^2}{8\pi m_\phi} = \frac{g^4}{8\pi m}$$

where we assumed $m_\chi \ll m_\phi \equiv m$ so the final-state particles are effectively massless. For fermionic decay channels $\phi \to \bar{\psi}\psi$ with Yukawa coupling $h\phi\bar{\psi}\psi$, one obtains $\Gamma_\phi = h^2 m/(8\pi)$.

Modified Klein-Gordon Equation

The backreaction of particle production on the inflaton condensate can be modeled by adding a friction term to the equation of motion:

$$\ddot{\phi} + (3H + \Gamma_\phi)\dot{\phi} + V'(\phi) = 0$$

The $\Gamma_\phi\dot{\phi}$ term represents the energy drain from the coherent oscillations into produced particles. The energy densities of the inflaton and radiation evolve as:

$$\dot{\rho}_\phi + 3H(1 + w)\rho_\phi = -\Gamma_\phi\rho_\phi$$

$$\dot{\rho}_r + 4H\rho_r = +\Gamma_\phi\rho_\phi$$

where $w = 0$ for a quadratic potential. The inflaton energy is drained at rate $\Gamma_\phi$ and transferred to radiation.

Reheating Completes When $\Gamma_\phi \sim H$

Initially $H \gg \Gamma_\phi$, so the decay is negligible compared to Hubble friction. As the universe expands, $H$ decreases (for matter domination, $H = 2/(3t)$). Reheating completes when the decay rate catches up to the expansion rate:

$$\Gamma_\phi \sim H_{\text{rh}} \quad \Longrightarrow \quad \Gamma_\phi^2 \sim \frac{\rho_{\text{rh}}}{3M_P^2}$$

Deriving the Reheating Temperature

At the moment of reheating, the energy density has been converted to a thermal bath of relativistic particles with:

$$\rho_{\text{rh}} = \frac{\pi^2}{30}g_* T_{\text{rh}}^4$$

Setting $\Gamma_\phi = H_{\text{rh}}$ and using $H^2 = \rho/(3M_P^2)$:

$$\Gamma_\phi^2 = \frac{\rho_{\text{rh}}}{3M_P^2} = \frac{\pi^2 g_*}{90} \frac{T_{\text{rh}}^4}{M_P^2}$$

$$T_{\text{rh}}^4 = \frac{90}{\pi^2 g_*}\Gamma_\phi^2 M_P^2$$

Reheating Temperature

$$\boxed{T_{\text{rh}} \sim \left(\frac{90}{\pi^2 g_*}\right)^{1/4}\sqrt{\Gamma_\phi M_P} \sim \sqrt{\Gamma_\phi M_P}}$$

Substituting $\Gamma_\phi = g^4/(8\pi m)$ for the scalar coupling:

$$T_{\text{rh}} \sim g^2\sqrt{\frac{m \cdot M_P}{8\pi}} \sim g^2\sqrt{m \cdot M_P}$$

For typical values $m \sim 10^{13}$ GeV and $g^2 \sim 10^{-3}$, this gives $T_{\text{rh}} \sim 10^{9}$ GeV. The reheating temperature is highly sensitive to the coupling constant, which is essentially a free parameter in most inflationary models.

Derivation 3: Parametric Resonance (Preheating)

The perturbative treatment above misses an important effect: when the inflaton oscillation amplitude is large, its coupling to other fields can cause explosive, non-perturbative particle production through parametric resonance. This phase, called preheating, can transfer energy far more efficiently than perturbative decay.

Setup: $g^2\phi^2\chi^2$ Coupling

Consider the interaction $\mathcal{L}_{\text{int}} = -\frac{1}{2}g^2\phi^2\chi^2$. The $\chi$ field satisfies:

$$\ddot{\chi} + 3H\dot{\chi} - \frac{\nabla^2\chi}{a^2} + (m_\chi^2 + g^2\phi^2)\chi = 0$$

Decomposing into Fourier modes $\chi_k$ and substituting $\phi(t) = \Phi(t)\sin(mt)$, with the rescaling $X_k = a^{3/2}\chi_k$ to absorb Hubble friction:

$$\ddot{X}_k + \left[\frac{k^2}{a^2} + m_\chi^2 + g^2\Phi^2\sin^2(mt) - \frac{9}{4}H^2 - \frac{3}{2}\dot{H}\right]X_k = 0$$

Neglecting the Hubble correction terms (valid when $m \gg H$) and setting $m_\chi = 0$, we use the identity $\sin^2(mt) = \frac{1}{2}(1 - \cos(2mt))$ and change variables to $z = mt$:

The Mathieu Equation

$$\boxed{\frac{d^2 X_k}{dz^2} + \left[A_k - 2q\cos(2z)\right]X_k = 0}$$

This is the Mathieu equation, with parameters:

$$A_k = \frac{k^2}{a^2 m^2} + 2q$$

$$q = \frac{g^2\Phi^2}{4m^2}$$

The Mathieu equation has a rich structure of instability bands: for certain ranges of $A_k$ and $q$, the solutions grow exponentially. In cosmology, we distinguish two regimes:

Narrow Resonance ($q \ll 1$)

When $q \ll 1$, the instability bands are narrow, centered near $A_k \approx \ell^2$ for integer $\ell$. The dominant instability is the first band near $A_k \approx 1$, corresponding to momenta $k^2/a^2 \approx m^2 - 2qm^2$. The growth rate is:

$$\mu_k \approx \frac{q}{2} \quad \Longrightarrow \quad X_k \propto e^{qmt/2}$$

This is essentially perturbative decay dressed up: the particle number grows linearly (since $n_k \propto |X_k|^2 \propto e^{qmt}$ with small $q$).

Broad Resonance ($q \gg 1$)

When $q \gg 1$, the instability bands overlap and the resonance becomes stochastic. This is the physically interesting regime for preheating. The key features are:

The effective mass of $\chi$ varies as $m_{\text{eff}}^2(t) = g^2\Phi^2\sin^2(mt)$, passing through zero twice per oscillation.
Each time $\phi$ passes through zero, $m_{\text{eff}} \to 0$ and there is a burst of non-adiabatic particle production, analogous to the Schwinger effect.
The particle occupation number grows exponentially:

$$n_k \propto e^{2\mu_k m t}$$

where the Floquet exponent $\mu_k$ can be $\mathcal{O}(0.1)$ for optimal momenta, giving enormous particle production in just $\sim 10$ oscillation periods.

The number density of produced $\chi$ particles after $N$ oscillations can be estimated as:

$$n_\chi \sim \frac{1}{(2\pi)^3}\int_0^{k_*} 4\pi k^2 \, e^{2\mu_k m t} \, dk$$

where $k_* \sim \sqrt{g\Phi m}$ is the maximum resonant momentum. The exponential growth terminates when backreaction effects become important: the produced particles modify the effective mass of $\chi$ and shift the resonance bands, or rescattering redistributes particles to non-resonant momenta.

Physical picture: Broad resonance is like a parametric amplifier. Each zero-crossing of $\phi$ acts as a “pump” that coherently amplifies the $\chi$ field. The amplification is Bose-enhanced: the probability of producing particles in a mode already occupied by $n_k$ quanta is proportional to $(1 + n_k)$, leading to exponential growth.

Derivation 4: Thermalization and the Reheating Temperature

After preheating (or perturbative decay) produces particles, the resulting state is highly non-thermal: the particles are concentrated in specific momentum bands determined by the resonance structure. Full thermalization requires redistribution of energy across all momenta through scattering processes.

Thermalization Timescale

The thermalization rate is governed by 2-to-2 scattering. For particles with number density $n$, scattering cross-section $\sigma$, and typical velocity $v$:

$$\Gamma_{\text{th}} = n\sigma v$$

Thermalization occurs when $\Gamma_{\text{th}} > H$. For gauge-mediated scattering with cross-section $\sigma \sim \alpha^2/E^2$ and relativistic particles ($v \sim 1$):

$$t_{\text{th}} \sim \frac{1}{n\sigma v} \sim \frac{E^2}{n\alpha^2}$$

Since the particles produced by preheating have typical energies $E \sim \sqrt{gm\Phi}$ and number density $n \sim n_\chi$, thermalization can occur relatively quickly if the coupling constants are not too small.

The BBN Constraint

Big Bang nucleosynthesis (BBN) requires the universe to be in thermal equilibrium by $T \sim 1$ MeV, with a radiation-dominated equation of state. The Friedmann equation at BBN gives:

$$H_{\text{BBN}}^2 = \frac{\pi^2 g_*}{90}\frac{T_{\text{BBN}}^4}{M_P^2}$$

$$H_{\text{BBN}} \sim \frac{T_{\text{BBN}}^2}{M_P} \sim \frac{(1\text{ MeV})^2}{2.4 \times 10^{18}\text{ GeV}} \sim 10^{-25}\text{ GeV}$$

For reheating to complete before BBN, we need $\Gamma_\phi > H_{\text{BBN}}$, which translates to:

BBN Lower Bound on Reheating Temperature

$$\boxed{T_{\text{rh}} \gtrsim 4 \text{ MeV}}$$

This is the absolute minimum reheating temperature consistent with the observed light element abundances. More precise analyses including neutrino decoupling strengthen this to $T_{\text{rh}} \gtrsim 4$ MeV.

The Gravitino Problem

In supersymmetric (SUSY) extensions of the Standard Model, the gravitino (the spin-3/2 superpartner of the graviton) is produced thermally at a rate proportional to $T$. The gravitino abundance is:

$$\frac{n_{3/2}}{s} \sim \frac{T_{\text{rh}}}{M_P}$$

If the gravitino is unstable with mass $m_{3/2} \sim 100$ GeV–1 TeV, it decays during or after BBN, potentially destroying the successful light element predictions. Requiring the gravitino abundance not to spoil BBN gives:

Gravitino Upper Bound

$$\boxed{T_{\text{rh}} \lesssim 10^9 \text{ GeV}}$$

This creates tension with thermal leptogenesis, which requires $T_{\text{rh}} \gtrsim 10^9$ GeV for generating the matter-antimatter asymmetry. Resolving this tension is an active area of research in SUSY cosmology.

The allowed window for the reheating temperature in supersymmetric models is therefore:

$$4 \text{ MeV} \lesssim T_{\text{rh}} \lesssim 10^9 \text{ GeV}$$

spanning some 15 orders of magnitude—a vast range that highlights our ignorance about this critical epoch.

Applications

Baryogenesis via Leptogenesis

The matter-antimatter asymmetry of the universe, quantified by $\eta_B = n_B/s \sim 6 \times 10^{-10}$, requires new physics beyond the Standard Model. One of the most attractive mechanisms is thermal leptogenesis (Fukugita & Yanagida, 1986): heavy right-handed neutrinos $N_i$ with masses $M_i$ are produced in the thermal bath and decay asymmetrically into leptons and Higgs bosons.

For the lightest right-handed neutrino mass $M_1 \gtrsim 10^9$ GeV (the Davidson-Ibarra bound), thermal production requires:

$$T_{\text{rh}} \gtrsim M_1 \gtrsim 10^9 \text{ GeV}$$

This is precisely the scale where the gravitino bound becomes relevant, creating the famous leptogenesis-gravitino tension in supersymmetric theories.

Dark Matter Production

Reheating also determines the abundance of dark matter. For Weakly Interacting Massive Particles (WIMPs), the standard thermal freeze-out calculation requires $T_{\text{rh}} \gg T_{\text{f.o.}} \sim m_{\text{DM}}/20$. For a 100 GeV WIMP, this means $T_{\text{rh}} \gg 5$ GeV.

However, dark matter can also be produced non-thermally during preheating itself, or through the direct decay of the inflaton: $\phi \to \text{DM} + \text{DM}$. In the latter case, the dark matter abundance is:

$$\Omega_{\text{DM}} \propto \frac{m_{\text{DM}} \, \text{Br}(\phi \to \text{DM})}{m_\phi} \cdot \frac{T_{\text{rh}}}{m_\phi}$$

This provides a direct connection between the inflaton's couplings, the reheating temperature, and the observed dark matter density.

Observational Signatures in the CMB

The duration of reheating affects how many e-folds of inflation correspond to observable scales. The pivot scale $k_* = 0.05$ Mpc$^{-1}$ left the horizon $N_*$ e-folds before the end of inflation, where:

$$N_* \approx 55 - \frac{1}{3}\ln\left(\frac{T_{\text{rh}}}{10^{16}\text{ GeV}}\right) - \frac{1-3w_{\text{rh}}}{12(1+w_{\text{rh}})}\ln\left(\frac{\rho_{\text{end}}}{\rho_{\text{rh}}}\right)$$

A lower reheating temperature means more e-folds of reheating, which shifts $N_*$ and consequently the predicted spectral index $n_s$ and tensor-to-scalar ratio $r$. Current Planck data constrains the duration of reheating for specific inflationary models.

Historical Context

1982

Dolgov & Linde; Abbott, Ellis & Nanopoulos

First systematic studies of inflaton decay and reheating. Introduced the perturbative approach with the $\Gamma_\phi$ friction term and estimated the reheating temperature. Recognized that the inflaton must couple to Standard Model fields for inflation to connect to the hot Big Bang.

1982

Albrecht, Steinhardt, Turner & Wilczek

Studied reheating in “new inflation” models and established the basic framework relating $T_{\text{rh}}$ to the inflaton mass and coupling constants.

1990

Traschen & Brandenberger

Discovered that non-perturbative effects can dramatically enhance particle production during inflaton oscillations. This was an early precursor to the preheating mechanism, showing that parametric amplification could be important.

1994

Kofman, Linde & Starobinsky (First Preheating Paper)

Introduced the concept of “preheating”—a non-perturbative phase of explosive particle production via parametric resonance, preceding the perturbative decay. Showed that this can be far more efficient than the old perturbative calculation.

1997

Kofman, Linde & Starobinsky (Comprehensive Analysis)

Detailed treatment of broad and narrow parametric resonance, the stochastic resonance regime, backreaction effects, and the Mathieu equation formalism. This landmark paper established the modern understanding of preheating.

2008+

Lattice Simulations Era

Groups using codes like LATTICEEASY (Felder & Tkachev) and DEFROST (Frolov) performed full 3D lattice simulations of preheating, confirming the basic picture while revealing new phenomena: turbulent thermalization, formation of oscillons, and non-equilibrium phase transitions.

Python Simulations

Below we numerically solve the coupled inflaton-radiation system and map the Mathieu instability chart. All simulations use only numpy for numerical computation.

Inflaton Oscillation and Perturbative Decay

Python

script.py127 lines

import numpy as np

# ============================================================
# Inflaton oscillation + perturbative decay into radiation
# Solve: d(rho_phi)/dt = -3H(1+w)*rho_phi - Gamma*rho_phi
#        d(rho_r)/dt   = -4H*rho_r + Gamma*rho_phi
#        H^2 = (rho_phi + rho_r) / (3*M_P^2)
# ============================================================

# Natural units: M_P = 1, m = 1e-6 M_P
M_P = 1.0
m = 1e-6          # inflaton mass
g2 = 1e-4         # coupling g^2
Gamma = g2**2 * m / (8 * np.pi)  # perturbative decay rate

print("Reheating via Perturbative Inflaton Decay")
print("=" * 55)
print(f"Inflaton mass m     = {m:.2e} M_P")
print(f"Coupling g^2        = {g2:.2e}")
print(f"Decay rate Gamma    = {Gamma:.4e} M_P")
print(f"T_rh estimate       ~ sqrt(Gamma * M_P) = {np.sqrt(Gamma * M_P):.4e} M_P")
print()

# Initial conditions: inflaton-dominated, rho_phi ~ m^2 Phi^2 / 2
phi0 = 0.1 * M_P
rho_phi0 = 0.5 * m**2 * phi0**2
rho_r0 = 1e-30  # tiny seed
H0 = np.sqrt((rho_phi0 + rho_r0) / (3.0 * M_P**2))

# Solve the field equation directly:
# phi'' + (3H + Gamma) phi' + m^2 phi = 0
# plus Friedmann: H^2 = (0.5 phi'^2 + 0.5 m^2 phi^2 + rho_r)/(3 M_P^2)
# d(rho_r)/dt = Gamma * phi'^2  (energy transfer)

# State: [phi, phi_dot, rho_r, a]
def derivs(state, t):
    phi_val, phi_dot, rho_r, a = state
    rho_phi_val = 0.5 * phi_dot**2 + 0.5 * m**2 * phi_val**2
    rho_total = rho_phi_val + rho_r
    if rho_total <= 0:
        rho_total = 1e-50
    H = np.sqrt(rho_total / (3.0 * M_P**2))

dphi = phi_dot
    dphi_dot = -(3.0 * H + Gamma) * phi_dot - m**2 * phi_val
    drho_r = -4.0 * H * rho_r + Gamma * phi_dot**2
    da = H * a
    return np.array([dphi, dphi_dot, drho_r, da])

# RK4 integrator
def rk4_step(state, t, dt, f):
    k1 = f(state, t)
    k2 = f(state + 0.5 * dt * k1, t + 0.5 * dt)
    k3 = f(state + 0.5 * dt * k2, t + 0.5 * dt)
    k4 = f(state + dt * k3, t + dt)
    return state + (dt / 6.0) * (k1 + 2*k2 + 2*k3 + k4)

# Time integration
N_steps = 80000
t_end = 4.0 / Gamma  # integrate well past reheating
dt = t_end / N_steps

state = np.array([phi0, 0.0, rho_r0, 1.0])
t = 0.0

# Storage for output
times = []
phi_vals = []
rho_phi_vals = []
rho_r_vals = []
T_vals = []

g_star = 106.75  # relativistic degrees of freedom

sample_every = 200
for i in range(N_steps):
    if i % sample_every == 0:
        phi_val, phi_dot, rho_r, a = state
        rho_phi_val = 0.5 * phi_dot**2 + 0.5 * m**2 * phi_val**2
        # Temperature from rho_r = (pi^2/30) g_* T^4
        if rho_r > 0:
            T_val = (30.0 * rho_r / (np.pi**2 * g_star))**0.25
        else:
            T_val = 0.0
        times.append(t * Gamma)  # in units of 1/Gamma
        phi_vals.append(phi_val / phi0)
        rho_phi_vals.append(rho_phi_val)
        rho_r_vals.append(rho_r)
        T_vals.append(T_val)

state = rk4_step(state, t, dt, derivs)
    # Prevent negative rho_r
    if state[2] < 0:
        state[2] = 0.0
    t += dt

# Print results table
print(f"{'t*Gamma':>10s}  {'phi/phi0':>10s}  {'rho_phi':>12s}  {'rho_r':>12s}  {'T/M_P':>12s}")
print("-" * 62)

# Sample ~20 rows from the data
n_data = len(times)
indices = np.linspace(0, n_data - 1, 25, dtype=int)
for idx in indices:
    print(f"{times[idx]:10.4f}  {phi_vals[idx]:10.4f}  {rho_phi_vals[idx]:12.4e}  {rho_r_vals[idx]:12.4e}  {T_vals[idx]:12.4e}")

# Find reheating time (when rho_r first exceeds rho_phi)
rh_time = None
for i in range(len(times)):
    if rho_r_vals[i] > rho_phi_vals[i] and rho_phi_vals[i] > 0:
        rh_time = times[i]
        T_rh = T_vals[i]
        break

print()
if rh_time is not None:
    print(f"Reheating occurs at t*Gamma = {rh_time:.4f}")
    print(f"Reheating temperature T_rh  = {T_rh:.4e} M_P")
    print(f"Analytic estimate:   T_rh  ~ {np.sqrt(Gamma * M_P):.4e} M_P")
else:
    print("Reheating not yet complete in simulation window")

# Final state
print(f"\nFinal rho_phi / rho_r = {rho_phi_vals[-1] / max(rho_r_vals[-1], 1e-100):.4e}")
print(f"Final temperature     = {T_vals[-1]:.4e} M_P")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Mathieu Instability Chart — Parametric Resonance Bands

Python

script.py130 lines

import numpy as np

# ============================================================
# Mathieu Equation Stability Chart
# d^2 X/dz^2 + [A - 2q cos(2z)] X = 0
#
# We compute the Floquet exponent mu(A,q) by integrating
# over one period and extracting the growth rate.
# ============================================================

print("Mathieu Equation Instability Chart")
print("Parametric Resonance Bands for Preheating")
print("=" * 55)

def compute_floquet(A, q, N_steps=2000):
    """Integrate Mathieu eq over one period [0, pi] and extract Floquet exponent."""
    dz = np.pi / N_steps

# Two linearly independent solutions
    # Solution 1: X=1, dX/dz=0 at z=0
    # Solution 2: X=0, dX/dz=1 at z=0
    X1, dX1 = 1.0, 0.0
    X2, dX2 = 0.0, 1.0

for i in range(N_steps):
        z = i * dz

omega2 = A - 2.0 * q * np.cos(2.0 * z)

# RK4 for X'' + omega2 X = 0
        # State: (X, V) where V = dX/dz
        def step_rk4(X, V):
            def f(X_val, V_val, z_val):
                return V_val, -(A - 2.0 * q * np.cos(2.0 * z_val)) * X_val

k1x, k1v = f(X, V, z)
            k2x, k2v = f(X + 0.5*dz*k1x, V + 0.5*dz*k1v, z + 0.5*dz)
            k3x, k3v = f(X + 0.5*dz*k2x, V + 0.5*dz*k2v, z + 0.5*dz)
            k4x, k4v = f(X + dz*k3x, V + dz*k3v, z + dz)

X_new = X + (dz/6.0)*(k1x + 2*k2x + 2*k3x + k4x)
            V_new = V + (dz/6.0)*(k1v + 2*k2v + 2*k3v + k4v)
            return X_new, V_new

X1, dX1 = step_rk4(X1, dX1)
        X2, dX2 = step_rk4(X2, dX2)

# Monodromy matrix M = [[X1, X2], [dX1, dX2]]
    trace = X1 + dX2

# Floquet exponent: eigenvalues of M are exp(+/- mu*pi)
    # |trace| > 2 means instability
    if abs(trace) > 2.0:
        mu = np.arccosh(abs(trace) / 2.0) / np.pi
    else:
        mu = 0.0

return mu

# Scan the (A, q) plane
N_A = 80
N_q = 60
A_max = 12.0
q_max = 8.0

A_vals = np.linspace(0.01, A_max, N_A)
q_vals = np.linspace(0.01, q_max, N_q)

print(f"\nScanning {N_A} x {N_q} = {N_A*N_q} points in (A, q) plane...")
print(f"A range: [0, {A_max}],  q range: [0, {q_max}]")

mu_grid = np.zeros((N_q, N_A))
for j in range(N_q):
    for i in range(N_A):
        mu_grid[j, i] = compute_floquet(A_vals[i], q_vals[j])

# Print the instability chart as ASCII
print(f"\nInstability Chart (Floquet exponent mu > 0 = unstable)")
print(f"Rows: q (bottom=0, top={q_max:.0f}), Cols: A (left=0, right={A_max:.0f})")
print()

# Print column header
header = "q\A  "
for i in range(0, N_A, 8):
    header += f"{A_vals[i]:5.1f}   "
print(header)
print("-" * (6 + N_A))

for j in range(N_q - 1, -1, -3):
    row = f"{q_vals[j]:4.1f} "
    for i in range(N_A):
        mu = mu_grid[j, i]
        if mu > 0.15:
            row += "#"
        elif mu > 0.05:
            row += "+"
        elif mu > 0.01:
            row += "."
        else:
            row += " "
    print(row)

print()
print("Legend: # = strong instability (mu > 0.15)")
print("        + = moderate instability (0.05 < mu < 0.15)")
print("        . = weak instability (0.01 < mu < 0.05)")
print("        (space) = stable region (mu ~ 0)")

# Print maximum Floquet exponents for selected q values
print(f"\n{'q':>6s}  {'max mu':>8s}  {'A at max mu':>12s}  {'Regime':>15s}")
print("-" * 48)
for j in range(0, N_q, 5):
    q = q_vals[j]
    mu_row = mu_grid[j, :]
    max_mu = np.max(mu_row)
    A_at_max = A_vals[np.argmax(mu_row)]
    regime = "narrow" if q < 1 else ("broad" if q > 3 else "intermediate")
    print(f"{q:6.2f}  {max_mu:8.4f}  {A_at_max:12.2f}  {regime:>15s}")

# Physical interpretation
print("\n" + "=" * 55)
print("Physical Interpretation for Preheating:")
print(f"  q = g^2 Phi^2 / (4 m^2)")
print(f"  For g=0.01, Phi=0.1 M_P, m=1e-6 M_P:")
print(f"  q = {0.01**2 * (0.1)**2 / (4 * (1e-6)**2):.2e}  (BROAD resonance)")
print(f"  Exponential growth: n_k ~ exp(2 mu m t)")
print(f"  For mu ~ 0.15, after 10 oscillations (mt ~ 60):")
print(f"  n_k ~ exp({2*0.15*60:.0f}) ~ 10^{2*0.15*60/np.log(10):.0f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Reheating Temperature and Equation of State Evolution

Python

script.py133 lines

import numpy as np

# ============================================================
# Equation of state during inflaton oscillations
# For V = lambda |phi|^n / n, <w> = (n-2)/(n+2)
# Also: full numerical solution of oscillating inflaton
# ============================================================

print("Equation of State for Oscillating Inflaton")
print("=" * 55)

# Part 1: Verify w = (n-2)/(n+2) numerically for various n
print(f"\n{'n':>4s}  {'w_analytic':>12s}  {'w_numeric':>12s}  {'rho scaling':>14s}")
print("-" * 48)

for n in [2, 3, 4, 6, 8]:
    w_analytic = (n - 2.0) / (n + 2.0)
    rho_exp = 6.0 * n / (n + 2.0)

# Numerically solve phi'' + V'(phi) = 0 (no expansion, just verify virial)
    # V(phi) = |phi|^n / n,  V'(phi) = |phi|^(n-1) * sign(phi)
    phi_val = 1.0
    dphi = 0.0
    dt_sim = 0.001
    N_sim = 50000
    KE_sum = 0.0
    PE_sum = 0.0

for step in range(N_sim):
        PE = abs(phi_val)**n / n
        KE = 0.5 * dphi**2

if step > 5000:  # skip transient
            KE_sum += KE
            PE_sum += PE

# Leapfrog / Verlet
        if abs(phi_val) > 1e-15:
            force = -abs(phi_val)**(n-1) * np.sign(phi_val)
        else:
            force = 0.0
        dphi += force * dt_sim
        phi_val += dphi * dt_sim

avg_KE = KE_sum / (N_sim - 5000)
    avg_PE = PE_sum / (N_sim - 5000)
    w_numeric = (avg_KE - avg_PE) / (avg_KE + avg_PE)

print(f"{n:4d}  {w_analytic:12.6f}  {w_numeric:12.6f}  rho ~ a^(-{rho_exp:.2f})")

# Part 2: Full inflaton decay with scale factor evolution
print("\n" + "=" * 55)
print("Full Reheating Evolution: phi(t), rho_phi(t), rho_r(t), T(t)")
print("=" * 55)

M_P = 1.0
m = 1e-6
g2_values = [1e-3, 1e-4, 1e-5]

for g2 in g2_values:
    Gamma = g2**2 * m / (8.0 * np.pi)
    T_rh_est = (90.0 / (np.pi**2 * 106.75))**0.25 * np.sqrt(Gamma * M_P)

print(f"\n--- g^2 = {g2:.0e} ---")
    print(f"Gamma = {Gamma:.4e},  T_rh(est) = {T_rh_est:.4e} M_P")

# Solve coupled system
    phi_val = 0.1 * M_P
    dphi = 0.0
    rho_r = 1e-40
    a = 1.0

dt_sim = 0.5 / m  # small compared to oscillation period
    N_sim = min(int(6.0 / Gamma / dt_sim), 500000)

t = 0.0
    print(f"  {'t*Gamma':>10s}  {'|phi|/phi0':>10s}  {'rho_phi':>12s}  {'rho_r':>12s}  {'T':>12s}")
    print("  " + "-" * 62)

rh_found = False
    print_every = max(N_sim // 20, 1)

for step in range(N_sim):
        rho_phi = 0.5 * dphi**2 + 0.5 * m**2 * phi_val**2
        rho_total = rho_phi + rho_r
        if rho_total < 1e-60:
            break
        H = np.sqrt(rho_total / (3.0 * M_P**2))

if step % print_every == 0:
            T_val = (30.0 * max(rho_r, 0) / (np.pi**2 * 106.75))**0.25
            print(f"  {t*Gamma:10.4f}  {abs(phi_val)/(0.1*M_P):10.4f}  {rho_phi:12.4e}  {max(rho_r,0):12.4e}  {T_val:12.4e}")

if not rh_found and rho_r > rho_phi and rho_phi > 0:
            T_rh_num = (30.0 * rho_r / (np.pi**2 * 106.75))**0.25
            rh_found = True

# RK4 step
        def rhs(st):
            p, v, rr, aa = st
            rp = 0.5*v**2 + 0.5*m**2*p**2
            rt = rp + max(rr, 0)
            if rt < 1e-60:
                return np.array([0,0,0,0])
            HH = np.sqrt(rt / (3.0*M_P**2))
            dp = v
            dv = -(3.0*HH + Gamma)*v - m**2*p
            drr = -4.0*HH*max(rr,0) + Gamma*v**2
            da = HH*aa
            return np.array([dp, dv, drr, da])

st = np.array([phi_val, dphi, rho_r, a])
        k1 = rhs(st)
        k2 = rhs(st + 0.5*dt_sim*k1)
        k3 = rhs(st + 0.5*dt_sim*k2)
        k4 = rhs(st + dt_sim*k3)
        st = st + (dt_sim/6.0)*(k1 + 2*k2 + 2*k3 + k4)

phi_val, dphi, rho_r, a = st
        if rho_r < 0:
            rho_r = 0.0
        t += dt_sim

if rh_found:
        print(f"  => Reheating T = {T_rh_num:.4e} M_P (analytic: {T_rh_est:.4e})")
    else:
        print(f"  => Reheating not complete in window")

print("\n" + "=" * 55)
print("Summary: T_rh ~ g^2 * sqrt(m * M_P)")
print("Higher coupling => faster decay => higher T_rh")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Inflationary Models Observational Signatures →