Cosmology Grand Potential

/-- **Euler relation (derived).** `T·s = ρ + P` is an algebraic identity of
the Legendre-transform structure — not an independent equilibrium postulate. -/
theorem potential_euler (P s : ℝ → ℝ) (x : ℝ) :
    x * s x = energyOf P s x + P x := by
  simp only [energyOf]
  ring

/-- **Gibbs–Duhem relation (derived).** Along any temperature trajectory
`T(t)`, the pressure obeys `p′ = s(T)·T′`: this is the chain rule applied to
`s = dP/dT`, not an independent postulate. -/
theorem potential_gibbs_duhem
    {P s T : ℝ → ℝ} {T' t : ℝ}
    (hP : HasDerivAt P (s (T t)) (T t))
    (hT : HasDerivAt T T' t) :
    HasDerivAt (fun u => P (T u)) (s (T t) * T') t := by
  simpa [Function.comp] using hP.comp t hT

/-- **THEOREM (plasma pressure closed form).** The log-kernel integrals
(`π⁴/45`, `7π⁴/360`, both derived via Mellin transforms) collapse the
pressure to `P = (π²/90)·(g_B + (7/8)·g_F)·T⁴`.  The `7/8` is the same
fermionic weight that appears in entropy and energy — here it comes out of
the pressure channel independently. -/
theorem plasmaPressure_eq (gB gF T : ℝ) :
    plasmaPressure gB gF T = π ^ 2 / 90 * (gB + 7 / 8 * gF) * T ^ 4 := by
  unfold plasmaPressure
  rw [RadiationEntropyRelation.boseLog_integral_value,
    RadiationEntropyRelation.fermiLog_integral_value]
  have hπ : (π : ℝ) ≠ 0 := Real.pi_ne_zero
  field_simp
  ring

/-- **THEOREM (the plasma is a potential fluid).** The temperature derivative
of the statistical-mechanical pressure is *exactly* the `radiationEntropy` of
`NeutrinoDilution` — the object previously built from the independent entropy
integrals `∫σ_B`, `∫σ_F`.  So `s = dP/dT` holds between two independently
derived statistical-mechanical quantities: the potential structure of the
plasma is a theorem, not a definition. -/
theorem plasmaPressure_potential (gB gF x : ℝ) :
    HasDerivAt (fun T => plasmaPressure gB gF T)
      (NeutrinoDilution.radiationEntropy gB gF x) x := by
  rw [NeutrinoDilution.radiationEntropy_eq]
  have hfun : (fun T => plasmaPressure gB gF T)
      = fun T => π ^ 2 / 90 * (gB + 7 / 8 * gF) * T ^ 4 :=
    funext fun T => plasmaPressure_eq gB gF T
  rw [hfun]
  have hpow : HasDerivAt (fun T : ℝ => T ^ 4) (4 * x ^ 3) x := by
    simpa using hasDerivAt_pow 4 x
  have h := hpow.const_mul (π ^ 2 / 90 * (gB + 7 / 8 * gF))
  have hval : π ^ 2 / 90 * (gB + 7 / 8 * gF) * (4 * x ^ 3)
      = 2 * π ^ 2 / 45 * (gB + 7 / 8 * gF) * x ^ 3 := by ring
  rw [hval] at h
  exact h

/-- **CAPSTONE (dilution from the potential).** `(T_ν/T_γ)³ = 4/11` from:
the coupled sector has a pressure potential with `s = dP/dT` (local
equilibrium, the *only* thermodynamic input), both sectors satisfy their FRW
continuity equations, and the boundary data (plasma dof `2+4 → 2` across e±
annihilation, shared temperature at decoupling).  Compared with
`EntropyConservationFRW.dilution_from_frw`, the Euler and Gibbs–Duhem
hypotheses are gone — they are theorems of the potential structure. -/
theorem dilution_from_potential
    {P s : ℝ → ℝ} {T a Tν : ℝ → ℝ} {sT T' a' Tν' ρν' : ℝ → ℝ}
    {t₁ t₂ : ℝ} {T₁ Tγ αν : ℝ}
    (hTt : ∀ t, T t ≠ 0) (hαν : αν ≠ 0) (hTνt : ∀ t, Tν t ≠ 0)
    (ha₂ : a t₂ ≠ 0) (hTγ : Tγ ≠ 0)
    (hP : ∀ t, HasDerivAt P (s (T t)) (T t))
    (hs : ∀ t, HasDerivAt s (sT t) (T t))
    (hT : ∀ t, HasDerivAt T (T' t) t) (ha : ∀ t, HasDerivAt a (a' t) t)
    (hTν : ∀ t, HasDerivAt Tν (Tν' t) t)
    (hρν : ∀ t, HasDerivAt (fun u => αν * Tν u ^ 4) (ρν' t) t)
    (hcont : ∀ t, a t * (T t * sT t * T' t)
        = -3 * a' t * (energyOf P s (T t) + P (T t)))
    (hcontν : ∀ t, a t * ρν' t
        = -3 * a' t * (αν * Tν t ^ 4 + αν * Tν t ^ 4 / 3))
    (hbefore : s (T t₁) = NeutrinoDilution.radiationEntropy 2 4 T₁)
    (hafter : s (T t₂) = NeutrinoDilution.radiationEntropy 2 0 Tγ)
    (hshare : Tν t₁ = T₁) :
    (Tν t₂ / Tγ) ^ 3 = 4 / 11 := by
  -- Adiabaticity of the coupled sector: derived from the potential structure.
  have hcons : NeutrinoDilution.radiationEntropy 2 4 T₁ * a t₁ ^ 3
      = NeutrinoDilution.radiationEntropy 2 0 Tγ * a t₂ ^ 3 := by
    have h := potential_entropy_constant hTt hP hs hT ha hcont t₁ t₂
    rw [hbefore, hafter] at h
    exact h
  -- Free streaming of the neutrino sector (derived in EntropyConservationFRW).
  have hfree : a t₂ * Tν t₂ = a t₁ * T₁ := by
    have h := EntropyConservationFRW.radiation_aT_constant hαν hTνt hTν ha
      hρν hcontν t₂ t₁
    rw [hshare] at h
    exact h
  exact NeutrinoDilution.dilution_from_entropy_conservation ha₂ hTγ hcons hfree

/-- **THEOREM (radiation equation of state derived).** `p = ρ/3` is forced by
the grand-canonical integrals — it is not an input anywhere in the chain. -/
theorem plasma_eos (gB gF T : ℝ) :
    plasmaPressure gB gF T = plasmaEnergy gB gF T / 3 := by
  rw [plasmaEnergy_eq, plasmaPressure_eq]
  ring