Cosmology Statistics Kernels

/-- **THEOREM (Bose–Einstein distribution).**  The ensemble-mean
occupation of a bosonic mode is `⟨n⟩_B = 1/(eᵗ−1)`: the weighted
geometric series `Σ n·xⁿ = x/(1−x)²` divided by `Z_B = (1−x)⁻¹`. -/
theorem boseOccupation_eq {t : ℝ} (ht : 0 < t) :
    boseOccupation t = 1 / (Real.exp t - 1) := by
  unfold boseOccupation
  rw [bosePartition_eq ht]
  simp only [boltzmannWeight_pow]
  have hx0 : (0 : ℝ) ≤ Real.exp (-t) := le_of_lt (Real.exp_pos _)
  have hx1 : Real.exp (-t) < 1 := exp_neg_lt_one ht
  rw [tsum_coe_mul_geometric_of_norm_lt_one
    (by rw [Real.norm_of_nonneg hx0]; exact hx1)]
  have hy : 1 < Real.exp t := one_lt_exp ht
  have hy0 : Real.exp t ≠ 0 := ne_of_gt (Real.exp_pos t)
  have hy1 : Real.exp t - 1 ≠ 0 := by linarith
  have h1x : 1 - Real.exp (-t) ≠ 0 := by linarith
  rw [Real.exp_neg]
  rw [Real.exp_neg] at h1x
  field_simp

/-- **THEOREM (Fermi–Dirac distribution).**  The ensemble-mean occupation
of a fermionic mode is `⟨n⟩_F = 1/(eᵗ+1)`, for every `t` (the two-state
sum needs no convergence condition). -/
theorem fermiOccupation_eq (t : ℝ) :
    fermiOccupation t = 1 / (Real.exp t + 1) := by
  unfold fermiOccupation boltzmannWeight
  rw [fermiPartition_eq]
  simp only [Finset.sum_range_succ, Finset.sum_range_zero,
    Nat.cast_zero, Nat.cast_one, zero_mul, one_mul, neg_zero,
    Real.exp_zero, zero_add, add_zero]
  have hy0 : Real.exp t ≠ 0 := ne_of_gt (Real.exp_pos t)
  have hpos : (0 : ℝ) < 1 + Real.exp (-t) := by positivity
  rw [Real.exp_neg]
  rw [Real.exp_neg] at hpos
  field_simp

/-- **THEOREM: the Bose pressure kernel is the log partition function.**
`−ln(1−e^{−t}) = ln Z_B(t)`.  The kernel that was a definition in
`PhaseSpaceReduction` is the grand-canonical `ln Z` of one mode. -/
theorem boseLogKernel_eq_log_partition {t : ℝ} (ht : 0 < t) :
    boseLogKernel t = Real.log (bosePartition t) := by
  rw [bosePartition_eq ht, Real.log_inv]
  rfl

/-- **THEOREM: the Fermi pressure kernel is the log partition function.**
`ln(1+e^{−t}) = ln Z_F(t)`. -/
theorem fermiLogKernel_eq_log_partition (t : ℝ) :
    fermiLogKernel t = Real.log (fermiPartition t) := by
  rw [fermiPartition_eq]
  rfl

/-- **THEOREM (grand-canonical consistency, Bose).**  The mean occupation
is minus the derivative of the log partition function with respect to the
dimensionless energy: `d/dt[−ln(1−e^{−t})] = −⟨n⟩_B(t)`.  The pressure
and energy kernels are therefore *one* input, not two. -/
theorem boseLogKernel_hasDerivAt {t : ℝ} (ht : 0 < t) :
    HasDerivAt boseLogKernel (-(boseOccupation t)) t := by
  have h1 : HasDerivAt (fun s : ℝ => -s) (-1) t := (hasDerivAt_id t).neg
  have h2 : HasDerivAt (fun s : ℝ => Real.exp (-s))
      (Real.exp (-t) * (-1)) t := (Real.hasDerivAt_exp (-t)).comp t h1
  have h3 : HasDerivAt (fun s : ℝ => 1 - Real.exp (-s))
      (0 - Real.exp (-t) * (-1)) t := (hasDerivAt_const t 1).sub h2
  have hlt : Real.exp (-t) < 1 := exp_neg_lt_one ht
  have hne : 1 - Real.exp (-t) ≠ 0 := by linarith
  have h4 := (h3.log hne).neg
  have heq : -((0 - Real.exp (-t) * (-1)) / (1 - Real.exp (-t)))
      = -(boseOccupation t) := by
    rw [boseOccupation_eq ht]
    have hy0 : Real.exp t ≠ 0 := ne_of_gt (Real.exp_pos t)
    have hy1 : Real.exp t - 1 ≠ 0 := by
      have := one_lt_exp ht; linarith
    rw [Real.exp_neg]
    rw [Real.exp_neg] at hne
    field_simp
    ring
  rw [← heq]
  exact h4

/-- **THEOREM (grand-canonical consistency, Fermi).**
`d/dt[ln(1+e^{−t})] = −⟨n⟩_F(t)`, for every `t`. -/
theorem fermiLogKernel_hasDerivAt (t : ℝ) :
    HasDerivAt fermiLogKernel (-(fermiOccupation t)) t := by
  have h1 : HasDerivAt (fun s : ℝ => -s) (-1) t := (hasDerivAt_id t).neg
  have h2 : HasDerivAt (fun s : ℝ => Real.exp (-s))
      (Real.exp (-t) * (-1)) t := (Real.hasDerivAt_exp (-t)).comp t h1
  have h3 : HasDerivAt (fun s : ℝ => 1 + Real.exp (-s))
      (0 + Real.exp (-t) * (-1)) t := (hasDerivAt_const t 1).add h2
  have hpos : (0 : ℝ) < 1 + Real.exp (-t) := by positivity
  have h4 := h3.log (ne_of_gt hpos)
  have heq : (0 + Real.exp (-t) * (-1)) / (1 + Real.exp (-t))
      = -(fermiOccupation t) := by
    rw [fermiOccupation_eq]
    have hy0 : Real.exp t ≠ 0 := ne_of_gt (Real.exp_pos t)
    rw [Real.exp_neg]
    rw [Real.exp_neg] at hpos
    field_simp
    ring
  rw [← heq]
  exact h4

/-- **Pauli bound.**  A fermionic mode is never more than singly
occupied on average: `⟨n⟩_F < 1`.  This is the statistical shadow of the
occupancy restriction `n ∈ {0,1}`. -/
theorem fermiOccupation_lt_one (t : ℝ) : fermiOccupation t < 1 := by
  rw [fermiOccupation_eq]
  have h := Real.exp_pos t
  rw [div_lt_one (by linarith)]
  linarith

/-- **CAPSTONE (pressure).**  The plasma pressure of the η_B chain equals
the phase-space integral of `T·ln Z_mode(E/T)` — the grand-canonical
pressure `P = (T/V)·ln Z` — with `Z_B = Σ_{n∈ℕ} e^{−nE/T}` and
`Z_F = Σ_{n∈{0,1}} e^{−nE/T}`.  The log kernels are gone as inputs; only
the Gibbs weight and the occupancy sets remain. -/
theorem plasmaPressure_from_partitionFunction (gB gF : ℝ) {T : ℝ}
    (hT : 0 < T) :
    phaseSpaceDensity 3 gB T (fun t => Real.log (bosePartition t))
      + phaseSpaceDensity 3 gF T (fun t => Real.log (fermiPartition t))
      = GrandPotential.plasmaPressure gB gF T := by
  have hB := phaseSpaceDensity_congr_pos gB T hT
    (fun t => Real.log (bosePartition t)) boseLogKernel
    (fun t ht => (boseLogKernel_eq_log_partition ht).symm)
  have hF := phaseSpaceDensity_congr_pos gF T hT
    (fun t => Real.log (fermiPartition t)) fermiLogKernel
    (fun t _ => (fermiLogKernel_eq_log_partition t).symm)
  rw [hB, hF]
  exact PhaseSpaceReduction.plasmaPressure_from_phaseSpace gB gF hT

/-- **CAPSTONE (energy).**  The plasma energy equals the phase-space
integral of `E·⟨n⟩(E/T)`: mean occupation times mode energy, with `⟨n⟩`
derived from the same partition functions. -/
theorem plasmaEnergy_from_occupation (gB gF : ℝ) {T : ℝ} (hT : 0 < T) :
    phaseSpaceDensity 3 gB T (fun t => t * boseOccupation t)
      + phaseSpaceDensity 3 gF T (fun t => t * fermiOccupation t)
      = GrandPotential.plasmaEnergy gB gF T := by
  have hB := phaseSpaceDensity_congr_pos gB T hT
    (fun t => t * boseOccupation t) boseEnergyKernel
    (fun t ht => (boseEnergyKernel_eq_occupation ht).symm)
  have hF := phaseSpaceDensity_congr_pos gF T hT
    (fun t => t * fermiOccupation t) fermiEnergyKernel
    (fun t _ => (fermiEnergyKernel_eq_occupation t).symm)
  rw [hB, hF]
  exact PhaseSpaceReduction.plasmaEnergy_from_phaseSpace gB gF hT

/-- Gibbs weight of the `n`-quanta microstate of a mode at dimensionless
energy `t = E/T` (equivalently `β·E`), chemical potential zero:
`w_n = e^{−n·t}`. -/
noncomputable def boltzmannWeight (t : ℝ) (n : ℕ) : ℝ := Real.exp (-(n * t))