refactor: Move statmech

Physlib.lean

@@ -382,6 +382,9 @@ public import Physlib.StatisticalMechanics.CanonicalEnsemble.Basic
 public import Physlib.StatisticalMechanics.CanonicalEnsemble.Finite
 public import Physlib.StatisticalMechanics.CanonicalEnsemble.Lemmas
 public import Physlib.StatisticalMechanics.CanonicalEnsemble.TwoState
+public import Physlib.StatisticalMechanics.MicroCanonicalEnsemble.Basic
+public import Physlib.StatisticalMechanics.MicroCanonicalEnsemble.IdealGas
+public import Physlib.StatisticalMechanics.MicroCanonicalEnsemble.ThermoQuantities
 public import Physlib.StringTheory.Basic
 public import Physlib.StringTheory.FTheory.SU5.Basic
 public import Physlib.StringTheory.FTheory.SU5.Charges.AnomalyFree

Physlib/StatisticalMechanics/MicroCanonicalEnsemble/Basic.lean

@@ -0,0 +1,50 @@
+/-
+Copyright (c) 2025 Alex Meiburg. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Alex Meiburg
+-/
+module
+
+public import Mathlib.Data.Fintype.Defs
+public import Mathlib.Data.Real.Basic
+public import Physlib.Meta.Sorry
+/-!
+
+## The Microcanonical Ensemble
+
+-/
+
+@[expose] public section
+
+section ThermodynamicEnsemble
+
+/-- The Hamiltonian for a microcanonical ensemble, parameterized by any extrinsic parameters of data
+    D. Since it's an arbitrary type, 'T' would be nice, but we use 'D' for data to avoid confusion
+    with temperature. We define Hamiltonians as an energy function on a real manifold, whose
+    dimension n possibly depends on the data D. Energy is a WithTop ℝ, ⊤ is used to mean "excluded
+    as a possibility in phase space". This formalization does exclude the possibility of discrete
+    degrees of freedom; it is not _completely_ general. A better formalization could have an
+    arbitrary measurable space.
+-/
+structure MicroHamiltonian (D : Type) where
+  /-- For extrinsic parameters D (e.g. the number of particles of different chemical species,
+    the shape of the box), how many continuous degrees of freedom are there? -/
+  dim : D → Type
+  /-- We require that the number of degrees of freedom is finite. -/
+  [dim_fin : ∀ d, Fintype (dim d)]
+  /-- Given the configuration, what is its energy? -/
+  H : {d : D} → (dim d → ℝ) → WithTop ℝ
+
+/-- We add the dim_fin to the instance cache so that things like the measure can be synthesized -/
+instance microHamiltonianFintype {D} (H : MicroHamiltonian D) (d : D) : Fintype (H.dim d) :=
+  H.dim_fin d
+
+/-- The standard microcanonical ensemble Hamiltonian, where the data is the particle number N and
+  the volume V. -/
+abbrev NVEHamiltonian := MicroHamiltonian (ℕ × ℝ)
+
+/-- Helper to get the number in an N-V Hamiltonian -/
+def NVEHamiltonian.N : (ℕ × ℝ) → ℕ := Prod.fst
+
+/-- Helper to get the volume in an N-V Hamiltonian -/
+def NVEHamiltonian.V : (ℕ × ℝ) → ℝ := Prod.snd

QuantumInfo/StatMech/IdealGas.lean → Physlib/StatisticalMechanics/MicroCanonicalEnsemble/IdealGas.lean

@@ -5,21 +5,28 @@ Authors: Alex Meiburg
 -/
 module
 
-public import QuantumInfo.StatMech.ThermoQuantities
+public import Physlib.StatisticalMechanics.MicroCanonicalEnsemble.ThermoQuantities
 public import Mathlib.Analysis.SpecialFunctions.Gaussian.FourierTransform
+/-!
 
+## Ideal gas as a Micro Canonical Ensemble
+
+In this module we give the
+-/
 @[expose] public section
 
 noncomputable section
 
 --! Specializing to an ideal gas of distinguishable particles.
 
-/-- The Hamiltonian for an ideal gas: particles live in a cube of volume V^(1/3), and each contributes an energy p^2/2.
-The per-particle mass is normalized to 1. -/
+/-- The Hamiltonian for an ideal gas: particles live in a cube of volume V^(1/3), and each
+  contributes an energy p^2/2. The per-particle mass is normalized to 1. -/
 def IdealGas : NVEHamiltonian where
-  --The dimension of the manifold is 6 times the number of particles: three for position, three for momentum.
+  --The dimension of the manifold is 6 times the number of particles: three for position, three for
+  --momentum.
   dim := fun (n,_) ↦ Fin n × (Fin 3 ⊕ Fin 3)
-  --The energy is ∞ if any positions are outside the cube, otherwise it's the sum of the momenta squared over 2.
+  --The energy is ∞ if any positions are outside the cube, otherwise it's the sum of the momenta
+  --squared over 2.
   H := fun {d} config ↦
     let (n,V) := d
     let R := V^(1/3:ℝ) / 2 --half-sidelength of a cubical box
@@ -37,28 +44,30 @@ variable (n : ℕ) {V β T : ℝ}
 set_option backward.isDefEq.respectTransparency false in
 open MeasureTheory in
 /-- The partition function Z for an ideal gas. -/
-theorem PartitionZ_eq (hV : 0 < V) (hβ : 0 < β) :
-    IdealGas.PartitionZ (n,V) β = V^n * (2 * Real.pi / β)^(3 * n / 2 : ℝ) := by
-  rw [PartitionZ, IdealGas]
+lemma partitionZ_eq (hV : 0 < V) (hβ : 0 < β) :
+    IdealGas.partitionZ (n,V) β = V ^ n * (2 * Real.pi / β) ^ (3 * n / 2 : ℝ) := by
+  rw [partitionZ, IdealGas]
   simp only [Finset.univ_product_univ, one_div, ite_eq_right_iff, WithTop.sum_eq_top,
     Finset.mem_univ, WithTop.coe_ne_top, and_false, exists_false, imp_false, not_forall, not_le,
     neg_mul]
   have h₀ : ∀ (config:Fin n × (Fin 3 ⊕ Fin 3) → ℝ) proof,
       ((if ∀ (i : Fin n) (ax : Fin 3), |config (i, Sum.inl ax)| ≤ V ^ (3 : ℝ)⁻¹ / 2 then
                   ∑ x : Fin n × Fin 3, config (x.1, Sum.inr x.2) ^ 2 / (2 :ℝ)
-                else ⊤) : WithTop ℝ).untop proof = ∑ x : Fin n × Fin 3, config (x.1, Sum.inr x.2) ^ 2 / (2 :ℝ) := by
+                else ⊤) : WithTop ℝ).untop proof =
+                ∑ x : Fin n × Fin 3, config (x.1, Sum.inr x.2) ^ 2 / (2 :ℝ) := by
     intro config proof
     rw [WithTop.untop_eq_iff]
     split_ifs with h
     · simp
     · simp [h] at proof
   simp only [h₀, dite_eq_ite]; clear h₀
-
-  let eq_pm : MeasurableEquiv ((Fin n × Fin 3 → ℝ) × (Fin n × Fin 3 → ℝ)) (Fin n × (Fin 3 ⊕ Fin 3) → ℝ) :=
-    let e1 := (MeasurableEquiv.sumPiEquivProdPi (α := fun (_ : (Fin n × Fin 3) ⊕ (Fin n × Fin 3)) ↦ ℝ))
-    let e2 := (MeasurableEquiv.piCongrLeft _ (MeasurableEquiv.prodSumDistrib (Fin n) (Fin 3) (Fin 3))).symm
+  let eq_pm : MeasurableEquiv ((Fin n × Fin 3 → ℝ) ×
+      (Fin n × Fin 3 → ℝ)) (Fin n × (Fin 3 ⊕ Fin 3) → ℝ) :=
+    let e1 := (MeasurableEquiv.sumPiEquivProdPi
+      (α := fun (_ : (Fin n × Fin 3) ⊕ (Fin n × Fin 3)) ↦ ℝ))
+    let e2 := (MeasurableEquiv.piCongrLeft _
+      (MeasurableEquiv.prodSumDistrib (Fin n) (Fin 3) (Fin 3))).symm
     e1.symm.trans e2
-
   have h_preserve : MeasurePreserving eq_pm := by
     unfold eq_pm
     -- fun_prop --this *should* be a fun_prop!
@@ -69,16 +78,13 @@ theorem PartitionZ_eq (hV : 0 < V) (hβ : 0 < β) :
     · apply MeasurePreserving.symm
       apply measurePreserving_sumPiEquivProdPi
   rw [← MeasurePreserving.integral_comp h_preserve eq_pm.measurableEmbedding]; clear h_preserve
-
   rw [show volume = Measure.prod volume volume from rfl]
-
   have h_eval_eq_pm : ∀ (x y i p_i), eq_pm (x, y) (i, Sum.inl p_i) = x (i, p_i) := by
     intros; rfl
   have h_eval_eq_pm' : ∀ (x y i m_i), eq_pm (x, y) (i, Sum.inr m_i) = y (i, m_i) := by
     intros; rfl
   simp_rw [h_eval_eq_pm, h_eval_eq_pm']
   clear h_eval_eq_pm h_eval_eq_pm'
-
   have h_measurable_box : Measurable fun (a : (Fin n × Fin 3 → ℝ))
       => ∃ x_1 x_2, V ^ (3⁻¹:ℝ) / 2 < |a (x_1, x_2)| := by
     simp_rw [← Classical.not_forall_not, not_not, not_lt, abs_le]
@@ -88,7 +94,6 @@ theorem PartitionZ_eq (hV : 0 < V) (hβ : 0 < β) :
     apply Measurable.forall
     intro j
     refine Measurable.comp (measurableSet_setOf.mp measurableSet_Icc) (measurable_pi_apply (i, j))
-
   have h_measurability : Measurable fun x : (Fin n × Fin 3 → ℝ) × (Fin n × Fin 3 → ℝ) =>
       if ∃ x_1 x_2, V ^ (3⁻¹:ℝ) / 2 < |x.1 (x_1, x_2)| then 0
       else Real.exp (-(β * ∑ x_1 : Fin n × Fin 3, x.2 (x_1.1, x_1.2) ^ 2 / 2)) := by
@@ -97,7 +102,6 @@ theorem PartitionZ_eq (hV : 0 < V) (hβ : 0 < β) :
       convert Measurable.comp h_measurable_box measurable_fst
     · fun_prop
     · fun_prop
-
   rw [MeasureTheory.integral_eq_lintegral_of_nonneg_ae]
   rotate_left
   · apply Filter.Eventually.of_forall
@@ -108,7 +112,6 @@ theorem PartitionZ_eq (hV : 0 < V) (hβ : 0 < β) :
   rw [MeasureTheory.lintegral_prod]; swap
   · exact (Measurable.comp (g := ENNReal.ofReal)
       ENNReal.measurable_ofReal h_measurability).aemeasurable
-
   conv =>
     enter [1, 1, 2, x, 2, y]
     rw [← ite_not _ _ (0:ℝ), ← boole_mul _ (Real.exp _)]
@@ -118,8 +121,6 @@ theorem PartitionZ_eq (hV : 0 < V) (hβ : 0 < β) :
     enter [1, 1, 2, x]
     simp only [not_exists, not_lt, Prod.mk.eta]
     rw [MeasureTheory.lintegral_const_mul' _ _ (ENNReal.ofReal_ne_top)]
-
-
   rw [MeasureTheory.lintegral_mul_const, ENNReal.toReal_mul]
   rw [← MeasureTheory.integral_eq_lintegral_of_nonneg_ae]
   rw [← MeasureTheory.integral_eq_lintegral_of_nonneg_ae]
@@ -144,7 +145,6 @@ theorem PartitionZ_eq (hV : 0 < V) (hβ : 0 < β) :
       fun_prop
     · fun_prop
     · fun_prop
-
   congr 1
   · --Volume of the box
     have h_integrand_prod : ∀ (a : Fin n × Fin 3 → ℝ),
@@ -166,65 +166,71 @@ theorem PartitionZ_eq (hV : 0 < V) (hβ : 0 < β) :
         (measurableSet_Icc (a := -(V ^ (3⁻¹ : ℝ) / 2)) (b := (V ^ (3⁻¹ : ℝ) / 2)))
       simp_rw [Set.indicator] at h_indicator
       simp_rw [abs_le, ← Set.mem_Icc, h_indicator]
-      simp
+      simp only [integral_const, MeasurableSet.univ, measureReal_restrict_apply, Set.univ_inter,
+        Real.volume_real_Icc, sub_neg_eq_add, add_halves, smul_eq_mul, mul_one, sup_eq_left,
+        ge_iff_le]
       positivity
     rw [h_integral_1d]; clear h_integral_1d
     rw [← Real.rpow_mul_natCast hV.le]
     field_simp
     simp
   · --Gaussian integral
     have h_gaussian :=
-      GaussianFourier.integral_rexp_neg_mul_sq_norm (V := PiLp 2 (fun (_ : Fin n × Fin 3) ↦ ℝ)) (half_pos hβ)
+      GaussianFourier.integral_rexp_neg_mul_sq_norm
+        (V := PiLp 2 (fun (_ : Fin n × Fin 3) ↦ ℝ)) (half_pos hβ)
     apply (Eq.trans ?_ h_gaussian).trans ?_
     · have := EuclideanSpace.volume_preserving_symm_measurableEquiv_toLp (Fin n × Fin 3)
       rw [← this.integral_comp (MeasurableEquiv.measurableEmbedding _)]
       congr! 3 with x
-      simp_rw [div_eq_inv_mul, ← Finset.mul_sum, ← mul_assoc, neg_mul, mul_comm, PiLp.norm_sq_eq_of_L2]
+      simp_rw [div_eq_inv_mul, ← Finset.mul_sum, ← mul_assoc, neg_mul, mul_comm,
+        PiLp.norm_sq_eq_of_L2]
       congr! 3
       simp only [Real.norm_eq_abs, sq_abs]
       congr
     · field_simp
       congr
-      simp
+      simp only [finrank_euclideanSpace, Fintype.card_prod, Fintype.card_fin, Nat.cast_mul,
+        Nat.cast_ofNat]
       ring_nf
 
 /-- The Helmholtz Free Energy A for an ideal gas. -/
-theorem HelmholtzA_eq (hV : 0 < V) (hT : 0 < T) : IdealGas.HelmholtzA (n,V) T =
+lemma helmholtzA_eq (hV : 0 < V) (hT : 0 < T) : IdealGas.helmholtzA (n,V) T =
     -n * T * (Real.log V + (3/2) * Real.log (2 * Real.pi * T)) := by
-  rw [HelmholtzA, PartitionZT, PartitionZ_eq n hV (one_div_pos.mpr hT), Real.log_mul,
+  rw [helmholtzA, partitionZT, partitionZ_eq n hV (one_div_pos.mpr hT), Real.log_mul,
     Real.log_pow, Real.log_rpow, one_div, div_inv_eq_mul]
   ring_nf
   all_goals positivity
 
-theorem ZIntegrable (hV : 0 < V) (hβ : 0 < β) : IdealGas.ZIntegrable (n,V) β := by
-  have hZpos : 0 < PartitionZ IdealGas (n, V) β := by
-    rw [PartitionZ_eq n hV hβ]
+lemma ZIntegrable (hV : 0 < V) (hβ : 0 < β) : IdealGas.ZIntegrable (n,V) β := by
+  have hZpos : 0 < partitionZ IdealGas (n, V) β := by
+    rw [partitionZ_eq n hV hβ]
     positivity
   constructor
   · apply MeasureTheory.Integrable.of_integral_ne_zero
-    rw [← PartitionZ]
+    rw [← partitionZ]
     exact hZpos.ne'
   · exact hZpos.ne'
 
 /-- The ideal gas law: PV = nRT. In our unitsless system, R = 1.-/
-theorem IdealGasLaw (hV : 0 < V) (hT : 0 < T) :
-    let P := IdealGas.Pressure (n,V) T;
+theorem ideal_gas_law (hV : 0 < V) (hT : 0 < T) :
+    let P := IdealGas.pressure (n,V) T;
     let R := 1;
     P * V = n * R * T := by
-  dsimp [Pressure]
+  dsimp [pressure]
   rw [← derivWithin_of_isOpen (s := Set.Ioi 0) isOpen_Ioi hV]
-  rw [derivWithin_congr (f := fun V' ↦ -n * T * (Real.log V' + (3/2) * Real.log (2 * Real.pi * T))) ?_ ?_]
+  rw [derivWithin_congr (f := fun V' ↦ -n * T *
+    (Real.log V' + (3/2) * Real.log (2 * Real.pi * T))) ?_ ?_]
   rw [derivWithin_of_isOpen (s := Set.Ioi 0) isOpen_Ioi hV]
   erw [deriv_mul (by fun_prop) (by fun_prop (disch := exact hV.ne'))]
   simp [field]
   ring_nf
   · intro _ hV'
     dsimp
-    rw [HelmholtzA_eq n hV' hT]
+    rw [helmholtzA_eq n hV' hT]
     ring_nf
-  · exact HelmholtzA_eq n hV hT
+  · exact helmholtzA_eq n hV hT
 
--- Now proving e.g. Boyle's Law ("for an ideal gas with a fixed particle number, P and V are inversely proportional")
--- is a trivial consequence of the ideal gas law.
+-- Now proving e.g. Boyle's Law ("for an ideal gas with a fixed particle number, P and V are
+-- inversely proportional") is a trivial consequence of the ideal gas law.
 
 end IdealGas