Cosmology Statistics Kernels Plasma Pressure From Partition Function

/-- **Bose partition function** (geometric series): for `t > 0`,
`Z_B(t) = (1 − e^{−t})⁻¹`. -/
theorem bosePartition_eq {t : ℝ} (ht : 0 < t) :
    bosePartition t = (1 - Real.exp (-t))⁻¹ := by
  unfold bosePartition
  simp only [boltzmannWeight_pow]
  exact tsum_geometric_of_lt_one (le_of_lt (Real.exp_pos _))
    (exp_neg_lt_one ht)

/-- **Fermi partition function** (two occupancy states):
`Z_F(t) = 1 + e^{−t}`, for every `t`. -/
theorem fermiPartition_eq (t : ℝ) :
    fermiPartition t = 1 + Real.exp (-t) := by
  unfold fermiPartition boltzmannWeight
  simp [Finset.sum_range_succ]

/-- **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 (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 (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

/-- **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