pith. machine review for the scientific record. sign in

IndisputableMonolith.Gravity.CausalKernelChain

IndisputableMonolith/Gravity/CausalKernelChain.lean · 306 lines · 10 declarations

show as:
view math explainer →

open module explainer GitHub source

Explainer status: pending

   1import Mathlib
   2import IndisputableMonolith.Gravity.CaldeiraLeggett
   3
   4/-!
   5# Causal-kernel chain: time-domain exponential kernel → transfer function → limits
   6
   7This module formalizes the **single-timescale exponential** memory kernel
   8(Debye/single-pole response) and its frequency-domain properties.
   9
  10Goal
  11----
  12Formalize, end-to-end:
  13
  14- the time-domain exponential kernel (causal, \(t \ge 0\)),
  15- its frequency response / transfer function \(H(i\omega)\),
  16- the steady-state (ω→0) and Newtonian (ω→∞) limits for the response \(C(\omega)=\Re H(i\omega)\).
  17
  18Scope / limitations
  19-------------------
  20This is the **Debye** (single-pole) realization only. The paper allows broader spectral
  21densities; those require additional Fourier/Laplace machinery and are not yet formalized.
  22-/
  23
  24namespace IndisputableMonolith
  25namespace Gravity
  26namespace CausalKernelChain
  27
  28open scoped Real
  29open Complex
  30
  31noncomputable section
  32
  33/-! ## Core analytic lemmas: ∫ exp(-a t) dt and exp(-a B) → 0 -/
  34
  35/-- Finite-interval complex exponential integral:
  36\[
  37\int_0^B e^{-a t}\,dt = \frac{1 - e^{-aB}}{a}.
  38\]
  39
  40We state it in the convenient `t • (-a)` form, where `t : ℝ` and `a : ℂ`. -/
  41theorem integral_exp_smul_neg (a : ℂ) (ha : a ≠ 0) (B : ℝ) :
  42    ∫ t in (0 : ℝ)..B, Complex.exp (t • (-a))
  43      = (Complex.exp (B • (-a)) - 1) * (-a)⁻¹ := by
  44  -- Antiderivative F(t) = exp(t•(-a)) * (-a)⁻¹
  45  have hderiv : ∀ x ∈ Set.uIcc (0 : ℝ) B,
  46      HasDerivAt (fun t : ℝ => Complex.exp (t • (-a)) * (-a)⁻¹)
  47        (Complex.exp (x • (-a))) x := by
  48    intro x hx
  49    have h_id : HasDerivAt (fun t : ℝ => t) (1 : ℝ) x := by
  50      simpa using (hasDerivAt_id x)
  51    have h_inner :
  52        HasDerivAt (fun t : ℝ => t • (-a)) ((1 : ℝ) • (-a)) x :=
  53      (HasDerivAt.smul_const (𝕜 := ℝ) (𝕜' := ℝ) (F := ℂ) (x := x) h_id (-a))
  54    have h_inner' : HasDerivAt (fun t : ℝ => t • (-a)) (-a) x := by
  55      simpa using h_inner
  56    have h_exp :
  57        HasDerivAt (fun t : ℝ => Complex.exp (t • (-a)))
  58          (Complex.exp (x • (-a)) * (-a)) x := by
  59      simpa [Function.comp] using ((Complex.hasDerivAt_exp (x • (-a))).comp x h_inner')
  60    have hmul :
  61        HasDerivAt (fun t : ℝ => (fun t => Complex.exp (t • (-a))) t * (-a)⁻¹)
  62          ((Complex.exp (x • (-a)) * (-a)) * (-a)⁻¹) x :=
  63      (HasDerivAt.mul_const (𝕜 := ℝ) (𝔸 := ℂ) (x := x) h_exp ((-a)⁻¹))
  64    -- cancel `(-a) * (-a)⁻¹` using `ha`
  65    simpa [ha] using hmul
  66
  67  have hcont : Continuous (fun t : ℝ => Complex.exp (t • (-a))) := by
  68    fun_prop
  69  have hint :
  70      IntervalIntegrable (fun t : ℝ => Complex.exp (t • (-a)))
  71        MeasureTheory.volume (0 : ℝ) B :=
  72    hcont.intervalIntegrable 0 B
  73
  74  have hFTC :=
  75    intervalIntegral.integral_eq_sub_of_hasDerivAt
  76      (a := (0 : ℝ)) (b := B)
  77      (f := fun t : ℝ => Complex.exp (t • (-a)) * (-a)⁻¹)
  78      (f' := fun t : ℝ => Complex.exp (t • (-a)))
  79      hderiv hint
  80
  81  -- Rewrite F(B) - F(0) into the desired algebraic form.
  82  calc
  83    ∫ t in (0 : ℝ)..B, Complex.exp (t • (-a))
  84        = Complex.exp (B • (-a)) * (-a)⁻¹ - Complex.exp ((0 : ℝ) • (-a)) * (-a)⁻¹ := by
  85            simpa using hFTC
  86    _ = Complex.exp (B • (-a)) * (-a)⁻¹ - (1 : ℂ) * (-a)⁻¹ := by simp
  87    _ = (Complex.exp (B • (-a)) - 1) * (-a)⁻¹ := by ring
  88
  89
  90/-- Helper: norm of `exp (-(↑B * a))` is `exp (-(B * Re a))`. -/
  91lemma norm_exp_neg_mul_ofReal (a : ℂ) (B : ℝ) :
  92    ‖Complex.exp (-( (B : ℂ) * a))‖ = Real.exp (-(B * a.re)) := by
  93  have hre : (-( (B : ℂ) * a)).re = -(B * a.re) := by
  94    simp [Complex.neg_re, Complex.mul_re, Complex.ofReal_re, Complex.ofReal_im, mul_comm, mul_left_comm]
  95  simpa [Complex.norm_exp, hre]
  96
  97
  98/-- If `Re(a) > 0`, then `exp (-(↑B * a)) → 0` as `B → +∞`. -/
  99theorem tendsto_exp_neg_mul_ofReal_atTop (a : ℂ) (ha : 0 < a.re) :
 100    Filter.Tendsto (fun B : ℝ => Complex.exp (-( (B : ℂ) * a))) Filter.atTop (nhds (0 : ℂ)) := by
 101  -- Reduce to the norm tending to 0.
 102  refine (tendsto_zero_iff_norm_tendsto_zero).2 ?_
 103
 104  have hmul : Filter.Tendsto (fun B : ℝ => B * a.re) Filter.atTop Filter.atTop := by
 105    simpa using ((Filter.tendsto_id).atTop_mul_const ha)
 106  have hneg : Filter.Tendsto (fun B : ℝ => -(B * a.re)) Filter.atTop Filter.atBot :=
 107    (Filter.tendsto_neg_atTop_atBot.comp hmul)
 108  have hexp : Filter.Tendsto (fun B : ℝ => Real.exp (-(B * a.re))) Filter.atTop (nhds 0) :=
 109    (Real.tendsto_exp_atBot.comp hneg)
 110
 111  have hnorm_eq :
 112      (fun B : ℝ => ‖Complex.exp (-( (B : ℂ) * a))‖) = fun B : ℝ => Real.exp (-(B * a.re)) := by
 113    funext B
 114    simpa using (norm_exp_neg_mul_ofReal a B)
 115
 116  simpa [hnorm_eq] using hexp
 117
 118
 119/-! ## Debye kernel → transfer function -/
 120
 121/-- The (complex) transfer function for a single-pole (Debye) response:
 122\[
 123H(i\omega) = 1 + \frac{\Delta}{1 + i\omega\tau}.
 124\]
 125This matches the structure in `Gravity.CaldeiraLeggett.TransferFunction`, but is a complex-valued
 126frequency response rather than its extracted real part. -/
 127def transfer_function_complex (H : CaldeiraLeggett.TransferFunction) (ω : ℝ) : ℂ :=
 128  (1 : ℂ) + (H.Δ : ℂ) / ((1 : ℂ) + Complex.I * (ω : ℂ) * (H.τ : ℂ))
 129
 130
 131/-- The Debye exponential kernel for a single-timescale response:
 132\[
 133\Gamma(t) = \frac{\Delta}{\tau} e^{-t/\tau},\quad t \ge 0.
 134\]
 135We treat it as a function on `ℝ` and integrate it on `[0,B]` (then take `B → ∞`). -/
 136def debye_kernel (H : CaldeiraLeggett.TransferFunction) (t : ℝ) : ℝ :=
 137  (H.Δ / H.τ) * Real.exp (-t / H.τ)
 138
 139
 140/-- Truncated (finite-horizon) frequency response contribution from the Debye kernel:
 141\[
 142K_B(\omega)=\int_0^B \Gamma(t)\,e^{-i\omega t}\,dt.
 143\]
 144The full transfer function is `1 + K_∞(ω)`. -/
 145def kernel_response_trunc (H : CaldeiraLeggett.TransferFunction) (ω B : ℝ) : ℂ :=
 146  ∫ t in (0 : ℝ)..B,
 147    ((debye_kernel H t : ℝ) : ℂ) * Complex.exp (-(Complex.I * (ω : ℂ) * (t : ℂ)))
 148
 149
 150/-!
 151### Bridge lemma (frequency-domain closed form)
 152
 153For τ>0, define `a = (1/τ) + i ω`. Then
 154
 155  exp(-t/τ) * exp(-i ω t) = exp(-(a * t)).
 156
 157The truncated integral can be evaluated in closed form using `integral_exp_smul_neg`,
 158then the `B → ∞` limit is obtained from `tendsto_exp_neg_mul_ofReal_atTop`.
 159-/
 160
 161/-! ### Laplace transform limit and transfer-function bridge -/
 162
 163/-- The complex exponent `a = (1/τ) + i ω` appearing in the Debye kernel transform. -/
 164def laplaceExponent (H : CaldeiraLeggett.TransferFunction) (ω : ℝ) : ℂ :=
 165  ((1 / H.τ : ℝ) : ℂ) + Complex.I * (ω : ℂ)
 166
 167
 168/-- Truncated Debye-kernel response tends to its closed form as `B → ∞`. -/
 169theorem kernel_response_limit (H : CaldeiraLeggett.TransferFunction) (ω : ℝ) :
 170    Filter.Tendsto (fun B => kernel_response_trunc H ω B) Filter.atTop
 171      (nhds ((H.Δ : ℂ) / ((1 : ℂ) + Complex.I * (ω : ℂ) * (H.τ : ℂ)))) := by
 172  -- Abbreviate `a = (1/τ) + iω`.
 173  set a : ℂ := laplaceExponent H ω with ha_def
 174
 175  have ha_re : 0 < a.re := by
 176    have h : 0 < (1 / H.τ) := one_div_pos.2 H.τ_pos
 177    -- `a.re = 1/τ` since `Re(iω)=0`.
 178    simpa [ha_def, laplaceExponent] using h
 179
 180  have ha : a ≠ 0 := by
 181    have hne : a.re ≠ 0 := ne_of_gt ha_re
 182    intro h0
 183    apply hne
 184    simpa [h0]
 185
 186  -- Closed form for each truncation bound `B`.
 187  have hclosed (B : ℝ) :
 188      kernel_response_trunc H ω B =
 189        (H.Δ / H.τ : ℂ) * ((Complex.exp (B • (-a)) - 1) * (-a)⁻¹) := by
 190    -- Rewrite the integrand into the `exp (t • (-a))` shape and apply `integral_exp_smul_neg`.
 191    unfold kernel_response_trunc
 192    -- Push the real kernel into `ℂ`, and turn `Real.exp` into `Complex.exp`.
 193    simp_rw [debye_kernel]
 194    simp_rw [Complex.ofReal_mul, Complex.ofReal_div, Complex.ofReal_exp]
 195    -- Combine the two exponentials.
 196    simp_rw [← Complex.exp_add, ← mul_assoc, ← mul_left_comm, ← mul_comm]
 197    -- Pull out the constant factor and evaluate the remaining exponential integral.
 198    -- (At this point the integrand is exactly `exp (t • (-a))` after simp.)
 199    simp [ha_def, laplaceExponent, Complex.real_smul, intervalIntegral.integral_const_mul,
 200      integral_exp_smul_neg a ha B, mul_assoc, mul_left_comm, mul_comm, add_assoc, add_left_comm,
 201      add_comm]
 202
 203  -- Reduce to the closed form and take `B → ∞`.
 204  have hExp :
 205      Filter.Tendsto (fun B : ℝ => Complex.exp (B • (-a))) Filter.atTop (nhds (0 : ℂ)) := by
 206    -- `B • (-a) = -((B:ℂ) * a)`, and `Re(a) > 0`.
 207    have :=
 208      tendsto_exp_neg_mul_ofReal_atTop a ha_re
 209    simpa [Complex.real_smul, mul_assoc, mul_left_comm, mul_comm] using this
 210
 211  have hExpSub :
 212      Filter.Tendsto (fun B : ℝ => Complex.exp (B • (-a)) - (1 : ℂ)) Filter.atTop (nhds (-1 : ℂ)) :=
 213    (hExp.sub tendsto_const_nhds)
 214
 215  have hExpMul :
 216      Filter.Tendsto (fun B : ℝ => (Complex.exp (B • (-a)) - (1 : ℂ)) * (-a)⁻¹) Filter.atTop
 217        (nhds ((-1 : ℂ) * (-a)⁻¹)) :=
 218    (Filter.Tendsto.mul_const (-a)⁻¹ hExpSub)
 219
 220  have hMain :
 221      Filter.Tendsto (fun B : ℝ => (H.Δ / H.τ : ℂ) * ((Complex.exp (B • (-a)) - 1) * (-a)⁻¹))
 222        Filter.atTop (nhds ((H.Δ / H.τ : ℂ) * ((-1 : ℂ) * (-a)⁻¹))) :=
 223    (Filter.Tendsto.const_mul (H.Δ / H.τ : ℂ) hExpMul)
 224
 225  -- Rewrite the limit into the desired `Δ / (1 + iωτ)` form.
 226  have hlim_simplify :
 227      (H.Δ / H.τ : ℂ) * ((-1 : ℂ) * (-a)⁻¹) =
 228        (H.Δ : ℂ) / ((1 : ℂ) + Complex.I * (ω : ℂ) * (H.τ : ℂ)) := by
 229    -- `(-1) * (-a)⁻¹ = a⁻¹`, then regroup as a division by `τ*a = 1 + iωτ`.
 230    have : ((-1 : ℂ) * (-a)⁻¹) = a⁻¹ := by
 231      -- `(-a)⁻¹ = -(a⁻¹)`
 232      simp
 233    -- Rewrite `(Δ/τ) * a⁻¹` as `Δ / (τ * a)`, and compute `τ * a`.
 234    -- Use `ha_def` to expand `a = (1/τ) + iω`.
 235    -- `τ * ((1/τ) + iω) = 1 + iωτ`.
 236    -- Finally, cast `Δ/τ` into `ℂ` consistently.
 237    simp [this, ha_def, laplaceExponent, div_div, div_eq_mul_inv, mul_add, add_mul,
 238      mul_assoc, mul_left_comm, mul_comm]
 239
 240  -- Assemble the tendsto statement.
 241  have hMain' :
 242      Filter.Tendsto (fun B : ℝ => kernel_response_trunc H ω B) Filter.atTop
 243        (nhds ((H.Δ : ℂ) / ((1 : ℂ) + Complex.I * (ω : ℂ) * (H.τ : ℂ)))) := by
 244    -- Rewrite by the pointwise closed form, then apply `hMain`.
 245    have hcongr :
 246        (fun B : ℝ => kernel_response_trunc H ω B) =
 247          fun B : ℝ => (H.Δ / H.τ : ℂ) * ((Complex.exp (B • (-a)) - 1) * (-a)⁻¹) := by
 248      funext B
 249      simpa [hclosed B]
 250    -- Transfer the limit.
 251    simpa [hcongr, hlim_simplify] using hMain
 252
 253  exact hMain'
 254
 255
 256/-- `transfer_function_complex` is exactly the Debye single-pole transfer function form. -/
 257theorem transfer_function_eq_one_plus_kernel (H : CaldeiraLeggett.TransferFunction) (ω : ℝ) :
 258    transfer_function_complex H ω =
 259      (1 : ℂ) + (H.Δ : ℂ) / ((1 : ℂ) + Complex.I * (ω : ℂ) * (H.τ : ℂ)) := by
 260  simp [transfer_function_complex]
 261
 262
 263/-- The paper’s real-valued response function is the real part of the complex transfer function. -/
 264theorem response_function_is_real_part (H : CaldeiraLeggett.TransferFunction) (ω : ℝ) :
 265    CaldeiraLeggett.response_function H ω = (transfer_function_complex H ω).re := by
 266  -- Unfold both sides.
 267  unfold CaldeiraLeggett.response_function transfer_function_complex
 268  -- Let the complex denominator be `w = 1 + i ω τ`.
 269  set w : ℂ := (1 : ℂ) + Complex.I * (ω : ℂ) * (H.τ : ℂ) with hw
 270  have wre : w.re = 1 := by
 271    -- `Re(i * real) = 0`.
 272    simp [hw, mul_assoc, mul_left_comm, mul_comm]
 273  have wnormSq : Complex.normSq w = 1 + (ω * H.τ) ^ (2 : ℕ) := by
 274    -- `normSq w = w.re^2 + w.im^2`, and here `w.re = 1`, `w.im = ωτ`.
 275    -- We compute directly using `normSq_apply`.
 276    have : w.im = ω * H.τ := by
 277      simp [hw, mul_assoc, mul_left_comm, mul_comm]
 278    -- Expand `normSq` and simplify.
 279    -- Use `pow_two` (as `x^2 = x*x`) in the reverse direction.
 280    calc
 281      Complex.normSq w = w.re * w.re + w.im * w.im := by
 282        simpa [Complex.normSq_apply]
 283      _ = (1 : ℝ) * 1 + (ω * H.τ) * (ω * H.τ) := by
 284        simp [wre, this]
 285      _ = 1 + (ω * H.τ) ^ (2 : ℕ) := by
 286        simp [pow_two, mul_assoc]
 287
 288  -- Reduce the real part using `Complex.div_re`.
 289  -- Since `H.Δ` is real, its imaginary part is 0.
 290  have hdiv :
 291      ((H.Δ : ℂ) / w).re = H.Δ / (1 + (ω * H.τ) ^ (2 : ℕ)) := by
 292    -- Apply `div_re` and simplify with `wre` and `wnormSq`.
 293    simp [Complex.div_re, wre, wnormSq, hw, mul_assoc, mul_left_comm, mul_comm, add_assoc,
 294      add_left_comm, add_comm]
 295
 296  -- Finish: real part of `1 + z` is `1 + Re(z)`.
 297  -- Note: `simp` can unfold `.re` of addition definitionaly once `hdiv` is in place.
 298  -- We rewrite the divisor `w` back into the original denominator.
 299  simp [hw, hdiv, Complex.add_re]
 300
 301end
 302
 303end CausalKernelChain
 304end Gravity
 305end IndisputableMonolith
 306

source mirrored from github.com/jonwashburn/shape-of-logic