Stating Schrodinger's equation

Stream: mathlib4

[Joseph Tooby-Smith]

In the project PhysLean, I am interested in making steps in the direction of quantum mechanics. One thing that would need to be done in this direction is state the Schrodinger's equation in Euclidean space:

For a twice-differentiable function $\psi : \mathbb{R}^3 \to  \mathbb{C}$, and a  smooth function $V :  \mathbb{R}^3 \to  \mathbb{R}$ (possibly defined everywhere except a finite set of points e.g. $1/||x||$) we say that $\psi$ satisfies the Schrodinger equation with potential $V$ if there exists an $E \in \mathbb{R}$ such that 

$(- \frac{\hbar^2}{2m} \nabla^2 + V (x)) \psi= E \psi$
Here $\hbar$ and $m$ are constants in $\mathbb{R}$, and $\nabla^2$ is the Laplacian operator on Euclidean space, such that in the standard coordinates $(x_1, x_2, x_3)$ of $\mathbb{R}^3$

$\nabla^2 = \frac{\partial^2}{\partial x_1^2}  + \frac{\partial^2}{\partial x_2^2}  + \frac{\partial^2}{\partial x_3^2}$

As a first goal I am interested in just writing down this problem in a way that can be later used to e.g. write down or find solutions to the Schrodinger equation.

Are there any suggestions of what should be (and shouldn't be) used from Mathlib to do this?

(This question follows on from this discussion #HepLean > Hydrogen atom)

[Joseph Tooby-Smith]

(There is a modification of this problem to the explicitly time-dependent version of Schrodinger's equation which shouldn't be too hard once the above is done)

[Tomas Skrivan]

Here is my version of time dependent Schrodinger equation

def SchrodingerEquation (m h : ℝ) (V : ℝ³ → ℝ) (ψ : ℝ → ℝ³ → ℂ) :=
  ∀ t x, deriv (ψ · x) t = (-h^2/(2*m)) • Δ (ψ t) x + V x • ψ t x

It should generalize to arbitrary vector valued wave function defined over inner product space.

code

import Mathlib.Analysis.Calculus.FDeriv.Basic
import Mathlib.Analysis.Calculus.Gradient.Basic
import Mathlib.Analysis.InnerProductSpace.PiL2
import Mathlib.LinearAlgebra.Trace
import Mathlib.Analysis.InnerProductSpace.Adjoint

notation "ℝ³" => EuclideanSpace ℝ (Fin 3)

-- simplest version of laplace
noncomputable
def laplace (f : ℝ³ → ℝ) (x : ℝ³) : ℝ :=
  LinearMap.trace ℝ ℝ³ (fderiv ℝ (gradient f) x)

-- component wise for complex valued function
noncomputable
def laplace' (f : ℝ³ → ℂ) (x : ℝ³) : ℂ :=
  ⟨laplace (fun x => (f x).re) x, laplace (fun x => (f x).im) x⟩

variable {𝕜 : Type*} [RCLike 𝕜]
  {E} [NormedAddCommGroup E] [InnerProductSpace 𝕜 E] [CompleteSpace E]
  {F} [NormedAddCommGroup F] [InnerProductSpace 𝕜 F] [CompleteSpace F]

-- continuous verion of `LinearMap.trace`
def ContinuousLinearMap.trace : (E →L[𝕜] E) →L[𝕜] 𝕜 := sorry

def ContinuousLinearMap.adjoint' : (E →L[𝕜] F) ≃L[𝕜] (F →L[𝕜] E) := sorry

open InnerProductSpace in
/-- Compute trace of hessian of `F` valued function -/
noncomputable
def hessianTrace (H : E →L[𝕜] E →L[𝕜] F) : F :=
  (toDual 𝕜 F).symm <| ContinuousLinearMap.trace ∘L (((ContinuousLinearMap.adjoint' (𝕜:=𝕜) (E:=E) (F:=F)) ∘L H).flip)

/-- laplace for general vector valued function - version where we take trace of hessian -/
noncomputable
def laplace'' (f : E → F) (x : E) : F :=  hessianTrace (fderiv 𝕜 (fderiv 𝕜 f) x)


variable (𝕜)
/-- vector valued gradient -/
noncomputable
def vgrad (f : E → F) (x : E) : F →L[𝕜] E := (fderiv 𝕜 f x).adjoint
variable {𝕜}

open InnerProductSpace in
/-- vector valued divergence -/
noncomputable
def vdiv (f : E → F →L[𝕜] E) (x : E) : F :=
  (toDual 𝕜 F).symm <| ContinuousLinearMap.trace ∘L (fderiv 𝕜 f x).flip

variable (𝕜)
/-- laplace for general vector valued function - laplace = divergence of gradient  -/
noncomputable
def laplace''' (f : E → F) (x : E) : F := vdiv (vgrad 𝕜 f) x
variable {𝕜}

prefix:max "Δ" => laplace'' (𝕜:=ℝ)

def SchrodingerEquation (m h : ℝ) (V : ℝ³ → ℝ) (ψ : ℝ → ℝ³ → ℂ) :=
  ∀ t x, deriv (ψ · x) t = (-h^2/(2*m)) • Δ (ψ t) x + V x • ψ t x

In particular, I'm missing continuous version of LinearMap.trace. Also having continuous equivalence in ContinuousLinearMap.adjoint would be useful, it can be done without it but the code looks a bit uglier.

[Tomas Skrivan]

I gave two possible definitions of laplace for functions between inner product spaces.

1. laplace is trace of hessian

   • you compute fderiv twice and then you have to compute component wise trace of bilinear map
   • pros: computing derivatives via fderiv is probably the easiest as most of the API is designed for it right now
   • cons: not sure how easy/difficult is to compute the component wise trace

2. laplace is divergence of gradient

   • you have to compute vector valued gradient vgrad and vector valued divergence vdiv
   • pros: this will generalize to other equations like heat equation with variable coefficients where you multiply the gradient by some position dependent matrix before you pass it to divergence
   • cons: you have to develop whole new API for vgrad and vdiv

[Joseph Tooby-Smith]

Thanks for this @Tomas Skrivan ! I'm going to have a go trying to adapt what you have to a specific case, and see where I can take it. I'll report back if I get to somewhere interesting.

[Tomas Skrivan]

I'm curious, what is the specific case?

[Tomas Skrivan]

The first stress test would be to show that planar wave is a solution for zero potential.

I would just copy API of LinearMap.trace for ContinuousLinearMap.trace and see if we can show the planar wave is a solution.

It would be nice if some mathlib expert could comment why there is no ContinuousLinearMap.trace or explain why it is not necessary. Or is it the usual explanation, no one needed it before?

[Joseph Tooby-Smith]

I was thinking of doing Landau levels i.e. a particle in a plane with a uniform magnetic field perpendicular to that plane. Just because it doesn't have the singularity issues in its potential.

The annoying thing is that you need to modify the Laplacian slightly for this case (the magnetic field gives single derivative terms in the Hamiltonian). The solutions are also quite nice - related to Hermite polynomials.

But yes - your suggestion of the case when $V = 0$ is probably a better first case! :)

[Joseph Tooby-Smith]

(Actually even easier then Landau levels - the harmonic oscillator would probably be the next best thing to try - then wouldn't need to change the Laplacian)

[Tomas Skrivan]

Here is my attempt at Schrodinger equation in electro-magnetic field

open InnerProductSpace

/-- Turn vector potential `A` into vector potential operator -/
noncomputable
def vecPotOp (A : ℝ³ → ℝ³) : (ℝ³ → ℂ) → (ℝ³ → ℝ³ →L[ℝ] ℂ) :=
  fun ψ x => ⟨⟨⟨fun v => ⟪A x,v⟫_ℝ * ψ x,sorry⟩,sorry⟩, sorry⟩

open Complex in
/-- `vecPotOp' A` is adjoint of `vecPotOp A` -/
noncomputable
def vecPotOp' (A : ℝ³ → ℝ³) : (ℝ³ → ℝ³ →L[ℝ] ℂ) → (ℝ³ → ℂ) :=
  fun v x => v x (A x)

/-- vector valued divergence - adjoint of `fderiv` -/
noncomputable
def vdiv' (f : E → E →L[𝕜] F) (x : E) : F := hessianTrace (fderiv 𝕜 f x)

open Complex in
def SchroedinterEquationInMagneticField (m h q : ℝ) (φ : ℝ³ → ℝ) (A : ℝ³ → ℝ³) (ψ : ℝ → ℝ³ → ℂ) :=
  let D  := (-h*I) • fderiv ℝ - q • vecPotOp A
  let D' := (-h*I) • vdiv' - q • vecPotOp' A
  ∀ t x,
     deriv (ψ · x) t = (-h^2/(2*m)) • D' (D (ψ t ·)) x + φ x • ψ t x

[Tomas Skrivan]

everything together for easy Lean 4 playground access

code

import Mathlib.Analysis.Calculus.FDeriv.Basic
import Mathlib.Analysis.Calculus.Gradient.Basic
import Mathlib.Analysis.InnerProductSpace.PiL2
import Mathlib.LinearAlgebra.Trace
import Mathlib.Analysis.InnerProductSpace.Adjoint

notation "ℝ³" => EuclideanSpace ℝ (Fin 3)

-- simplest version of laplace
noncomputable
def laplace (f : ℝ³ → ℝ) (x : ℝ³) : ℝ :=
  LinearMap.trace ℝ ℝ³ (fderiv ℝ (gradient f) x)

-- component wise for complex valued function
noncomputable
def laplace' (f : ℝ³ → ℂ) (x : ℝ³) : ℂ :=
  ⟨laplace (fun x => (f x).re) x, laplace (fun x => (f x).im) x⟩

variable {𝕜 : Type*} [RCLike 𝕜]
  {E} [NormedAddCommGroup E] [InnerProductSpace 𝕜 E] [CompleteSpace E]
  {F} [NormedAddCommGroup F] [InnerProductSpace 𝕜 F] [CompleteSpace F]

-- continuous verion of `LinearMap.trace`
def ContinuousLinearMap.trace : (E →L[𝕜] E) →L[𝕜] 𝕜 := sorry

def ContinuousLinearMap.adjoint' : (E →L[𝕜] F) ≃L[𝕜] (F →L[𝕜] E) := sorry

open InnerProductSpace in
/-- Compute trace of hessian of `F` valued function -/
noncomputable
def hessianTrace (H : E →L[𝕜] E →L[𝕜] F) : F :=
  (toDual 𝕜 F).symm <| ContinuousLinearMap.trace ∘L (((ContinuousLinearMap.adjoint' (𝕜:=𝕜) (E:=E) (F:=F)) ∘L H).flip)

-- laplace for general vector valued function
noncomputable
def laplace'' (f : E → F) (x : E) : F :=  hessianTrace (fderiv 𝕜 (fderiv 𝕜 f) x)

prefix:max "Δ" => laplace'' (𝕜:=ℝ)

def SchroedinterEquation (m h : ℝ) (V : ℝ³ → ℝ) (ψ : ℝ → ℝ³ → ℂ) :=
  ∀ t x, deriv (ψ · x) t = (-h^2/(2*m)) • Δ (ψ t) x + V x • ψ t x


variable (𝕜)
/-- vector valued gradient -/
noncomputable
def vgrad (f : E → F) (x : E) : F →L[𝕜] E := (fderiv 𝕜 f x).adjoint
variable {𝕜}

open InnerProductSpace in
/-- vector valued divergence -/
noncomputable
def vdiv (f : E → F →L[𝕜] E) (x : E) : F :=
  (toDual 𝕜 F).symm <| ContinuousLinearMap.trace ∘L (fderiv 𝕜 f x).flip


open InnerProductSpace

noncomputable
def vecPotOp (A : ℝ³ → ℝ³) : (ℝ³ → ℂ) → (ℝ³ → ℝ³ →L[ℝ] ℂ) :=
  fun ψ x => ⟨⟨⟨fun v => ⟪A x,v⟫_ℝ * ψ x,sorry⟩,sorry⟩, sorry⟩

open Complex in
/-- `vecPotOp' A` is adjoint of `vecPotOp A` -/
noncomputable
def vecPotOp' (A : ℝ³ → ℝ³) : (ℝ³ → ℝ³ →L[ℝ] ℂ) → (ℝ³ → ℂ) :=
  fun v x => v x (A x)

open InnerProductSpace in
/-- vector valued divergence -/
noncomputable
def vdiv' (f : E → E →L[𝕜] F) (x : E) : F := hessianTrace (fderiv 𝕜 f x)

open Complex in
def SchroedinterEquationInMagneticField (m h q : ℝ) (φ : ℝ³ → ℝ) (A : ℝ³ → ℝ³) (ψ : ℝ → ℝ³ → ℂ) :=
  let D  := (-h*I) • fderiv ℝ - q • vecPotOp A
  let D' := (-h*I) • vdiv' - q • vecPotOp' A
  ∀ t x,
     deriv (ψ · x) t = (-h^2/(2*m)) • D' (D (ψ t ·)) x + φ x • ψ t x

[Tomas Skrivan]

(deleted)

[Tomas Skrivan]

Also I would like to point out that

fderiv : (E → F) → (E → E →L[𝕜] F)
vdiv' : (E → E →L[𝕜] F) → (E → F)

are transpose of each other in a particular sense.

We can interpret fderiv as function between test functions and vdiv' as distributional derivative

vgrad : 𝒟(E,F) → 𝒟(E, E →L[𝕜] F)
vdiv' : 𝒟'(E, E →L[𝕜] F)  → 𝒟'(E,F)

(here I'm using vector valued test functions and distributions, reference)
For function between test functions we can talk about transpose and in this sense -vdiv' is transpose of fderiv

[Tomas Skrivan]

In this light D' is really transpose of D and it makes sense why people write $( -ih\nabla - q A)^2$ which usually means $( -ih\nabla - q A)^\dagger ( -ih\nabla - q A)$ .

[Tomas Skrivan]

(I did bunch of edits because in my first attempt I defined D using vgrad but then realized that I can't multiply it by complex unit I so I redid everything using fderiv instead of vgrad.)

[Tomas Skrivan]

I'm trying to show that plane wave is a solution

def planeWave (ω : ℝ) (k : ℝ³) (t : ℝ) (x : ℝ³) : ℂ := Complex.exp (- (⟪k,x⟫_ℝ - ω*t) * I)

and I'm trying to show that

theorem planeWave_is_solution (m k : ℝ) (ω : ℝ) (k : ℝ³) :
    SchroedinterEquation m (fun _ => 0) (planeWave ω k) := ...

(one should pick appropriate ω based on m and k, I just havn't done that yet)

After unfolding all the definitions you get this equality

(fderiv ℝ (fun x_1 => cexp (-(↑⟪k, x⟫_ℝ - ↑ω * ↑x_1) * I)) t) 1 =
    (-ℏ ^ 2 / (2 * m)) • hessianTrace (fderiv ℝ (fderiv ℝ fun x => cexp (-(↑⟪k, x⟫_ℝ - ↑ω * ↑t) * I)) x) +
      (fun x => 0) x • cexp (-(↑⟪k, x⟫_ℝ - ↑ω * ↑t) * I)

I have managed to eliminate derivative on lhs and inner derivative on rhs.
Now I'm hitting a problem that there is no ContinuousLinearMap for fun x => ⟪y, x⟫_ℝ and that it is differentiable as a function in y + what is its derivative.

[Tomas Skrivan]

Here is my progress, there are bunch of missing theorems or they are not in the correct form.

code

import Mathlib.Analysis.Calculus.FDeriv.Basic
import Mathlib.Analysis.Calculus.Gradient.Basic
import Mathlib.Analysis.InnerProductSpace.PiL2
import Mathlib.LinearAlgebra.Trace
import Mathlib.Analysis.InnerProductSpace.Adjoint

import Mathlib.Analysis.Calculus.FDeriv.Mul
import Mathlib.Analysis.Calculus.FDeriv.Add
import Mathlib.Analysis.SpecialFunctions.ExpDeriv
import Mathlib.Analysis.InnerProductSpace.Calculus

notation "ℝ³" => EuclideanSpace ℝ (Fin 3)

-- simplest version of laplace
noncomputable
def laplace (f : ℝ³ → ℝ) (x : ℝ³) : ℝ :=
  LinearMap.trace ℝ ℝ³ (fderiv ℝ (gradient f) x)

-- component wise for complex valued function
noncomputable
def laplace' (f : ℝ³ → ℂ) (x : ℝ³) : ℂ :=
  ⟨laplace (fun x => (f x).re) x, laplace (fun x => (f x).im) x⟩

variable {𝕜 : Type*} [RCLike 𝕜]
  {E} [NormedAddCommGroup E] [InnerProductSpace 𝕜 E] [CompleteSpace E]
  {F} [NormedAddCommGroup F] [InnerProductSpace 𝕜 F] [CompleteSpace F]

-- continuous verion of `LinearMap.trace`
def ContinuousLinearMap.trace : (E →L[𝕜] E) →L[𝕜] 𝕜 := sorry

def ContinuousLinearMap.adjoint' : (E →L[𝕜] F) ≃L[𝕜] (F →L[𝕜] E) := sorry

open InnerProductSpace in
/-- Compute trace of hessian of `F` valued function -/
noncomputable
def hessianTrace (H : E →L[𝕜] E →L[𝕜] F) : F :=
  (toDual 𝕜 F).symm <| ContinuousLinearMap.trace ∘L (((ContinuousLinearMap.adjoint' (𝕜:=𝕜) (E:=E) (F:=F)) ∘L H).flip)

-- laplace for general vector valued function
noncomputable
def laplace'' (f : E → F) (x : E) : F :=  hessianTrace (fderiv 𝕜 (fderiv 𝕜 f) x)

prefix:max "Δ" => laplace'' (𝕜:=ℝ)

/-- Reduced Planc's constant -/
opaque ℏ : ℝ

def SchroedinterEquation (m : ℝ) (V : ℝ³ → ℝ) (ψ : ℝ → ℝ³ → ℂ) :=
  ∀ t x, deriv (ψ · x) t = (-ℏ^2/(2*m)) • Δ (ψ t) x + V x • ψ t x

open InnerProductSpace

-- there is `fderiv_exp` and `deriv_exp` but no `fderiv_exp`
@[simp]
theorem fderiv_cexp
    {𝕜 : Type*} [NontriviallyNormedField 𝕜] [NormedAlgebra 𝕜 ℂ]
    {E} [NormedAddCommGroup E] [NormedSpace 𝕜 E]
    (f : E → ℂ) (hc : Differentiable 𝕜 f) :
    fderiv 𝕜 (fun x => Complex.exp (f x)) = fun x => Complex.exp (f x) • fderiv 𝕜 f x :=
  sorry

@[fun_prop]
theorem differentiable_complex_ofReal_over_real : Differentiable ℝ Complex.ofReal := sorry

@[fun_prop]
theorem differentiableAt_complex_ofReal_over_real (x : ℝ) : DifferentiableAt ℝ Complex.ofReal x := sorry

@[simp]
theorem fderiv_complex_ofReal
    {E} [NormedAddCommGroup E] [NormedSpace ℝ E]
    (f : E → ℝ) (hf : Differentiable ℝ f) :
    fderiv ℝ (fun x => Complex.ofReal (f x))
    =
    fun x => ⟨⟨⟨fun dx => Complex.ofReal (fderiv ℝ f x dx),sorry⟩,sorry⟩,sorry⟩  := sorry

-- simp does not pick up this with `fderiv_complex_ofReal`
@[simp]
theorem fderiv_complex_ofReal_simple :
    fderiv ℝ (Complex.ofReal : ℝ → ℂ)
    =
    fun x => ⟨⟨⟨fun dx => Complex.ofReal dx,sorry⟩,sorry⟩,sorry⟩ := by
  rw[fderiv_complex_ofReal]; simp; fun_prop

attribute [fun_prop] differentiable_inner

-- non applied version of `fderiv_inner_apply`
@[simp]
theorem fderiv_inner
    {E} [NormedAddCommGroup E] [InnerProductSpace ℝ E] [CompleteSpace E]
    {F} [NormedAddCommGroup F] [InnerProductSpace ℝ F] [CompleteSpace F]
    (f : E → F) (hf : Differentiable ℝ f)
    (g : E → F) (hg : Differentiable ℝ g) :
    fderiv ℝ (fun x => ⟪f x, g x⟫_ℝ)
    =
    fun x => ⟨⟨⟨fun dx => ⟪f x, fderiv ℝ g x dx⟫_ℝ + ⟪fderiv ℝ f x dx, g x⟫_ℝ,sorry⟩,sorry⟩,sorry⟩ := by
  funext x; ext dx
  simp (disch:=fun_prop) [fderiv_inner_apply]

-- simp misses `fderiv_inner` because of eta expansion issues in DiscTree
@[simp]
theorem fderiv_inner_left_apply
    {E} [NormedAddCommGroup E] [InnerProductSpace ℝ E] [CompleteSpace E] (y : E) :
    fderiv ℝ (fun x => ⟪y, x⟫_ℝ)
    =
    fun x => ⟨⟨⟨fun dx => ⟪y,dx⟫_ℝ,sorry⟩,sorry⟩,sorry⟩ := by
  rw[fderiv_inner]
  simp
  all_goals fun_prop

-- `fun_prop` seems to have trouble using `Differentiable.comp` because of `Function.comp`
@[fun_prop]
theorem Differentiable.comp'
    {E} [NormedAddCommGroup E] [NormedSpace ℝ E]
    {F} [NormedAddCommGroup F] [NormedSpace ℝ F]
    {G} [NormedAddCommGroup G] [NormedSpace ℝ G]
    (f : E → F) (hf : Differentiable ℝ f)
    (g : F → G) (hg : Differentiable ℝ g) :
    Differentiable ℝ (fun x => g (f x)) := by fun_prop

-- `fun_prop` seems to have trouble using `DifferentiableAt.comp` because of `Function.comp`
@[fun_prop]
theorem DifferentiableAt.comp'
    {E} [NormedAddCommGroup E] [NormedSpace ℝ E]
    {F} [NormedAddCommGroup F] [NormedSpace ℝ F]
    {G} [NormedAddCommGroup G] [NormedSpace ℝ G]
    (x : E)
    (f : E → F) (hf : DifferentiableAt ℝ f x)
    (g : F → G) (hg : DifferentiableAt ℝ g (f x)) :
    DifferentiableAt ℝ (fun x => g (f x)) x := by fun_prop

-- `fun_prop` seems to have trouble using `Differentiable.id` because of `id`
@[fun_prop]
theorem Differentiable.id'
    {E} [NormedAddCommGroup E] [NormedSpace ℝ E] :
    Differentiable ℝ (fun x : E => x) := by fun_prop

-- `Differentiable.mul` is not sufficiently general
@[fun_prop]
theorem Differentiable.mul'
    {𝕜 : Type*} [NontriviallyNormedField 𝕜]
    {𝕜' : Type*} [NontriviallyNormedField 𝕜'] [NormedAlgebra 𝕜 𝕜']
    {E} [NormedAddCommGroup E] [NormedSpace 𝕜 E] [CompleteSpace E]
    (f g : E → 𝕜') :
    Differentiable 𝕜 (fun x : E => f x * g x) := by sorry

-- `fderiv_mul` is not sufficiently general
@[simp]
theorem fderiv_mul''
    {𝕜 : Type*} [NontriviallyNormedField 𝕜]
    {𝕜' : Type*} [NontriviallyNormedField 𝕜'] [NormedAlgebra 𝕜 𝕜']
    {E} [NormedAddCommGroup E] [NormedSpace 𝕜 E] [CompleteSpace E]
    (f g : E → 𝕜') (hf : Differentiable 𝕜 f) (hg : Differentiable 𝕜 g) :
    fderiv 𝕜 (fun x : E => f x * g x)
    =
    fun x => f x • fderiv 𝕜 g x + g x • fderiv 𝕜 f x := by sorry

@[simp]
theorem smul_zero_clm (c : ℂ) : c • (0 : ℝ³ →L[ℝ] ℂ) = 0 := by ext dx; simp

open Complex
noncomputable
def planeWave (ω : ℝ) (k : ℝ³) (t : ℝ) (x : ℝ³) : ℂ := Complex.exp (- (⟪k,x⟫_ℝ - ω*t) * I)

theorem planeWave_is_solution (m k : ℝ) (ω : ℝ) (k : ℝ³) :
    SchroedinterEquation m (fun _ => 0) (planeWave ω k) := by
  unfold SchroedinterEquation laplace'' planeWave
  intro t x; unfold deriv
  simp (disch:=fun_prop) [-PiLp.inner_apply,fderiv_sub,-fderiv_eq_smul_deriv]

[Joseph Tooby-Smith]

Hi @Tomas Skrivan ! All this looks great - thanks! I'm going to go through it later today and try and fill in some of the sorries.

I've been working on a simpler case corresponding to the 1d Quantum Harmonic Oscillator. This is easier since you can just use deriv instead of fderiv. I've got as far as stating the Schrodinger equation, the eigenfunctions and eigenvalues and proving that the eigenfunctions and eigenvalues are indeed eigenfunctions and eigenvalues. I haven't proved any completeness relations or orthogonality conditions etc.

It exceeds the character length of Zulip (didn't even know this existed until now :)), so I can't paste it here - but it is in this PR.

[Joseph Tooby-Smith]

Naturally any suggestions for improvements either here or on the Github would be much appreciated :).

[Tomas Skrivan]

Looks good, few comments:

1. Not sure if defining hermitePolynomial with three cases has any advantage over two cases. If you do it in two cases then you have to use the fact that 0-1=0 which I guess is that reason for three cases.
2. You compute lots the derivatives by hand unnecessarily, simp can do most of the work for you, see code bellow.
3. Formulate the derivative of hermitePolynomial in "compositional" form (do the same for Complex.ofReal

lemma deriv_hermitePolynomial' (x : ℝ)
    (f : ℝ → ℝ) (hf : DifferentiableAt ℝ f x) (n : ℕ)  :
    deriv (fun x => hermitePolynomial n (f x)) x
    =
    (2 * n * hermitePolynomial (n - 1) (f x)) * deriv f x := ...

with this computing derivatives can be done with simp, for example

lemma deriv_hermitePolynomial_succ (m ℏ ω : ℝ) (n : ℕ) :
    deriv (fun x => (hermitePolynomial (n + 1) (Real.sqrt (m * ω / ℏ) * x))) =
    fun x =>
    (Real.sqrt (m * ω / ℏ)) * 2 * (n + 1) *
    hermitePolynomial n (Real.sqrt (m * ω / ℏ) * x) := by

  funext x
  simp (disch:=fun_prop)
  rw[deriv_const_mul (hd:=by fun_prop)]
  simp;
  ring

the proof is morally just funext x; simp (disch:=fun_prop); ring but simp fails to apply deriv_const_mul because of indexing issues.

---

Here is my version of deriv_hermitePolynomial and derived compositional forms of that lemma.

code

import Mathlib.Analysis.Calculus.Deriv.Mul
import Mathlib.Analysis.Calculus.Deriv.Add
import Mathlib.Analysis.Calculus.Deriv.Pow


/-- The Hermite polynomials defined via their recursive relation. -/
def hermitePolynomial : (n : ℕ) → (x : ℝ) → ℝ
  | 0, _ => 1
  | 1, x => 2 * x
  | (n + 2), x => 2 * x * hermitePolynomial (n + 1) x - 2 * (n + 1) * hermitePolynomial n x

@[fun_prop]
lemma hermitePolynomial_differentiableAt : (n : ℕ) → (x : ℝ) →
    DifferentiableAt ℝ (hermitePolynomial n) x
  | 0, x => by
    simp [hermitePolynomial]
  | 1, x => by
    simp [hermitePolynomial]
    fun_prop
  | n + 2, x => by
    simp [hermitePolynomial]
    have ih := hermitePolynomial_differentiableAt (n + 1) x
    have ih' := hermitePolynomial_differentiableAt n x
    fun_prop

lemma hermitePolynomial_add_succ_succ (n : ℕ) :
    hermitePolynomial (n + 2) = 2 • (fun x => x) * hermitePolynomial (n + 1) -
    (2 * (n + 1) : ℝ) • hermitePolynomial n := by
  simp [hermitePolynomial]
  rfl

lemma hermitePolynomial_succ (n : ℕ) (x : ℝ) :
    hermitePolynomial (n + 1) x = 2 * x * hermitePolynomial n x - 2 * n * hermitePolynomial (n - 1) x := by
  induction n <;> simp[hermitePolynomial]

-- simp runs `fun_prop` with some odd reducibility setting so `differentiableAt_id` fails to unify
@[fun_prop]
theorem differentiableAt_id''
    {𝕜 : Type*} [NontriviallyNormedField 𝕜]
    {E : Type*} [NormedAddCommGroup E] [NormedSpace 𝕜 E] (x : E) :
    DifferentiableAt 𝕜 (fun y => y) x := by fun_prop

lemma deriv_hermitePolynomial (n : ℕ) :
    deriv (hermitePolynomial n) = 2 * n * hermitePolynomial (n - 1) := by
  match n with
  | 0 =>
    funext x; simp [hermitePolynomial]
  | 1 =>
    funext x; simp[hermitePolynomial]
    rw[deriv_const_mul (hd:=by fun_prop)] -- ugh this does not fire because of DiscrTree eta epansion :(
    simp (disch:=fun_prop)
  | n+2 =>
    have hn := deriv_hermitePolynomial n
    have hn' := deriv_hermitePolynomial (n+1)
    funext x; simp[hermitePolynomial]
    simp (disch:=fun_prop) [deriv_sub]
    rw[deriv_const_mul (hd:=by fun_prop)] -- ugh again hitting DiscrTree issue
    simp (disch:=fun_prop) [hn,hn',hermitePolynomial_succ]
    ring

lemma fderiv_hermitePolynomial
    {E : Type*} [NormedAddCommGroup E] [NormedSpace ℝ E] (x : E)
    (f : E → ℝ) (hf : DifferentiableAt ℝ f x) (n : ℕ)  :
    fderiv ℝ (fun x => hermitePolynomial n (f x)) x
    =
    (2 * n * hermitePolynomial (n - 1) (f x)) • fderiv ℝ f x := by
  have h := fderiv_comp (𝕜:=ℝ) (g:=hermitePolynomial n) (f:=f) (hg:=by fun_prop) (hf:=by fun_prop)
  simp +unfoldPartialApp [Function.comp] at h
  ext dx
  simp[deriv_hermitePolynomial,h]
  ring

@[simp]
lemma deriv_hermitePolynomial' (x : ℝ)
    (f : ℝ → ℝ) (hf : DifferentiableAt ℝ f x) (n : ℕ)  :
    deriv (fun x => hermitePolynomial n (f x)) x
    =
    (2 * n * hermitePolynomial (n - 1) (f x)) * deriv f x := by
  unfold deriv
  rw[fderiv_hermitePolynomial (hf:=by fun_prop)]
  rfl


-- we can now compute derivatives of `hermitePolynomial` with simp
lemma deriv_hermitePolynomial_succ (m ℏ ω : ℝ) (n : ℕ) :
    deriv (fun x => (hermitePolynomial (n + 1) (Real.sqrt (m * ω / ℏ) * x))) =
    fun x =>
    (Real.sqrt (m * ω / ℏ)) * 2 * (n + 1) *
    hermitePolynomial n (Real.sqrt (m * ω / ℏ) * x) := by

  funext x
  simp (disch:=fun_prop)
  rw[deriv_const_mul (hd:=by fun_prop)] -- ugh again hitting DiscrTree is
  simp;
  ring

[Joseph Tooby-Smith]

@Tomas Skrivan Great! I also just found that Mathlib has the probabilists version of Hermite polynomials here and I think inspiration can also be taken from there.

It's still probably worth defining the "physicists" version explicitly though.

[Tomas Skrivan]

Few tips for computing derivatives with simp. I usually run simp it with

set_option trace.Meta.Tactic.simp.discharge true in

to see which theorem fails because simp can't discharging differentiability. There seems to be some reducibility issues with id and Function.comp while running fun_prop as discharger. For example this works

example : DifferentiableAt ℝ (fun y => 2 * y * hermitePolynomial (n + 1) y) x := by fun_prop

but the same goal fails when simp tries to discharge it with fun_prop.

---

It is also very useful to run

set_option trace.Meta.Tactic.simp.unify true in

when you know a certain theorem should be applied but it is not. So turn that option when calling for example simp only [deriv_sub]. Without only you will get lots of noise.

---

One thing to keep in mind is that simp has issues with eta expansions when indexing theorems. It won't apply deriv_const_mul or deriv_mul to deriv (HMul.hMull 2). The issue is not with unification so it won't come up trace.Meta.Tactic.simp.unify. So either use rw[deriv_const_mul] to simplify deriv (HMul.hMull 2) as in the above code or you can state a new theorem deriv (HMul.hMull c) = (HMul.hMull c) and mark it with simp attribute.

[Tomas Skrivan]

Joseph Tooby-Smith said:

Tomas Skrivan
It's still probably worth defining the "physicists" version explicitly though.

I think it is a good idea to have probabilistic and physicist version such that we can state formulas as they are in textbooks. Not sure though if we should pick one canonical version for proofs though.

[Kevin Buzzard]

Tomas Skrivan said:

It would be nice if some mathlib expert could comment why there is no ContinuousLinearMap.trace or explain why it is not necessary. Or is it the usual explanation, no one needed it before?

For general vector spaces there is no concept of trace, one needs finite-dimensionality (or some restriction on the class of continuous linear maps you're considering). And in the finite-dimensional case continuity is automatic (all linear maps are continuous).

[Tomas Skrivan]

Kevin Buzzard said:

Tomas Skrivan said:

It would be nice if some mathlib expert could comment why there is no ContinuousLinearMap.trace or explain why it is not necessary. Or is it the usual explanation, no one needed it before?

For general vector spaces there is no concept of trace, one needs finite-dimensionality (or some restriction on the class of continuous linear maps you're considering). And in the finite-dimensional case continuity is automatic (all linear maps are continuous).

That was my suspicion but then I do no know how the API is supposed to work. How can I easily turn f : E →ₗ[ℝ] E →ₗ[ℝ] → F into E →L[ℝ] E →L[ℝ] → F for E and F finite dimensional? The linear algebra concepts are done with E →ₗ[ℝ] → F but derivatives are with E →L[ℝ] → F and I just do not know how to make them interact nicely.

(In SciLean I just work with unbundled maps everywhere and use fun_prop to prove any property on the fly when necessary but mathlib is not designed this way so I'm a bit lost)

[Kevin Buzzard]

I would imagine the analysts have had to deal with this problem at some point, so hopefully they will chime in. On my part I can mention docs#IsModuleTopology.continuous_of_linearMap , which says that any linear map between vector spaces (or more generally modules) with topologies is automatically continuous if source and target have what I call the "module topology", although this a theorem rather than an upgrading of a linear map to a continuous linear map.

[Tomas Skrivan]

Ohh some time ago I used docs#LinearMap.continuous_of_finiteDimensional what is the difference? Looks like that the IsModuleTopology version is lot less restrictive on the ring. Anyway, I'm really interested in the practicalities of actually working with the (continuous) linear maps, mathematically it all make sense.

[Kevin Buzzard]

The theorem you quote says that if the ground field is complete then continuity is automatic. My theorem is that if the ground ring is an arbitrary topological ring and the vector space is also arbitrary but but the topology on the vector space satisfies some kind of natural condition then continuity is automatic. If your ground field is complete then you should probably be using the one you suggested as then you won't have to check that the the topology on the module is the natural one. The point is that given any module over any topological ring there is some kind of natural topology one can put on it, but this seems like overkill here.

[Joseph Tooby-Smith]

Here is my refactoring of the physicists Hermite polynomials based on @Tomas Skrivan's comments above, and on Mathlibs Hermite polynomials:

code

import Mathlib.Analysis.Calculus.Deriv.Polynomial
import Mathlib.Analysis.InnerProductSpace.Basic
/-!
# Physicists Hermite Polynomial
This file may eventually be upstreamed to Mathlib.
-/

open Polynomial
namespace PhysLean

/-- The Physicists Hermite  polynomial. -/
noncomputable def physHermite : ℕ → Polynomial ℤ
  | 0 => 1
  | n + 1 => 2 • X * physHermite n - derivative (physHermite n)

lemma physHermite_succ (n : ℕ) :
    physHermite (n + 1) = 2 • X * physHermite n - derivative (physHermite n) := by
  simp [physHermite]

lemma physHermite_eq_iterate (n : ℕ) :
    physHermite n = (fun p => 2 * X * p - derivative p)^[n] 1 := by
  induction n with
  | zero => rfl
  | succ n ih => simp [Function.iterate_succ_apply', ← ih, physHermite_succ]

@[simp]
lemma physHermite_zero : physHermite 0 = C 1 := rfl

lemma physHermite_one : physHermite 1 = 2 * X := by
  rw [physHermite_succ, physHermite_zero]
  simp

lemma derivative_physHermite_succ : (n : ℕ) →
    derivative (physHermite (n + 1)) = 2 * (n + 1) • physHermite n
  | 0 => by
    simp [physHermite_one]
  | n + 1 => by
    rw [physHermite_succ]
    simp only [nsmul_eq_mul, Nat.cast_ofNat, derivative_sub, derivative_mul, derivative_ofNat,
      zero_mul, derivative_X, mul_one, zero_add, Nat.cast_add, Nat.cast_one]
    rw [derivative_physHermite_succ]
    simp only [physHermite_succ, nsmul_eq_mul, Nat.cast_ofNat, Nat.cast_add, Nat.cast_one,
      derivative_mul, derivative_ofNat, zero_mul, derivative_add, derivative_natCast,
      derivative_one, add_zero, zero_add]
    ring

lemma derivative_physHermite : (n : ℕ) →
    derivative (physHermite n) = 2 * n • physHermite (n - 1)
  | 0 => by
    simp
  | n + 1 => by
    rw [derivative_physHermite_succ]
    simp

lemma physHermite_succ' (n : ℕ) :
    physHermite (n + 1) = 2 • X * physHermite n - 2 * n • physHermite (n - 1) := by
  rw [physHermite_succ, derivative_physHermite]

noncomputable def physHermiteFun (n : ℕ) : ℝ → ℝ := fun x => aeval x (physHermite n)

lemma physHermiteFun_eq_aeval_physHermite (n : ℕ) :
    physHermiteFun n = fun x => aeval x (physHermite n) := rfl

lemma physHermiteFun_succ_fun (n : ℕ) :
    physHermiteFun (n + 1) = 2 • (fun x => x) *
    physHermiteFun n - (2 * n : ℝ) • physHermiteFun (n - 1) := by
  ext x
  simp [physHermite_succ', aeval, physHermiteFun, mul_assoc]

lemma physHermiteFun_succ (n : ℕ) :
    physHermiteFun (n + 1) = fun x =>  2 • x *
    physHermiteFun n x -
    (2 * n : ℝ) • physHermiteFun (n - 1) x := by
  ext x
  simp [physHermite_succ', aeval, mul_assoc, physHermiteFun]

@[fun_prop]
lemma physHermiteFun_differentiableAt : (n : ℕ) → (x : ℝ) →
    DifferentiableAt ℝ (physHermiteFun n) x
  | 0, x => by
    simp [physHermite, physHermiteFun_eq_aeval_physHermite]
  | 1, x => by
    simp [physHermite, physHermiteFun_eq_aeval_physHermite, aeval]
    fun_prop
  | n + 2, x => by
    rw [physHermiteFun_succ]
    have ih := physHermiteFun_differentiableAt n x
    have ih' := physHermiteFun_differentiableAt (n + 1) x
    fun_prop

lemma deriv_physHermiteFun (n : ℕ) :
    deriv (physHermiteFun n) = 2 * n * (physHermiteFun (n - 1)) := by
  ext x
  rw [physHermiteFun_eq_aeval_physHermite]
  rw [Polynomial.deriv_aeval (physHermite n)]
  rw [derivative_physHermite]
  simp [aeval, mul_assoc, physHermiteFun_eq_aeval_physHermite]


lemma fderiv_physHermiteFun
    {E : Type*} [NormedAddCommGroup E] [NormedSpace ℝ E] (x : E)
    (f : E → ℝ) (hf : DifferentiableAt ℝ f x) (n : ℕ)  :
    fderiv ℝ (fun x => physHermiteFun n (f x)) x
    = (2 * n * physHermiteFun  (n - 1) (f x)) • fderiv ℝ f x := by
  have h := fderiv_comp (𝕜:=ℝ) (g:=physHermiteFun n) (f:=f) (hg:=by fun_prop) (hf:=by fun_prop)
  simp +unfoldPartialApp [Function.comp] at h
  ext dx
  simp [deriv_physHermiteFun,h]
  ring

@[simp]
lemma deriv_physHermiteFun' (x : ℝ)
    (f : ℝ → ℝ) (hf : DifferentiableAt ℝ f x) (n : ℕ)  :
    deriv (fun x => physHermiteFun n (f x)) x
    = (2 * n * physHermiteFun  (n - 1) (f x)) * deriv f x := by
  unfold deriv
  rw[fderiv_physHermiteFun (hf:=by fun_prop)]
  rfl

end PhysLean

[Tomas Skrivan]

I find it strange that it is defined as polynomial over integers. I would expect to define it as Polynomial ℝ and having coercion CoeFun (Polynomial ℝ) (fun _ => (ℝ → ℝ)). I guess there is a reason for doing it this way. I see that it is more flexible but I do not like that we have physHermite and physHermiteFun

[Joseph Tooby-Smith]

I did it this way to match Polynomial.hermite. I agree having both physHermite and physHermiteFun does not seem ideal. Though I don't think writing physHermiteFun as fun x => aeval x (physHermite n) everywhere is also great.

[Joseph Tooby-Smith]

@Tomas Skrivan Going to push this just to PhysLean today. I'll add your name to the 'authors...' bit of the PhysicistsHermite file, unless you don't want me too.

I haven't managed to show completeness of the eigenfunctions yet - this seems pretty tricky to do in general (so much so that most resources for the Harmonic oscillator skip any such proof). Ideally it would be down to showing that $\{x^p e^{- c * x^2} \mid p \in \mathbb{N}\}$ for fixed real $0<c$ form a basis of the 1d Hilbert space of square integrable (equivalence classes of) functions. Or rather that if $f$ is square integrable such that if

$\int^\infty_{-\infty} f(x) x^p e^{- c * x^2} = 0$

for all $p$ then $f$ is almost-everywhere zero.

[Tomas Skrivan]

I just skimmed through the code but it looks good to me.

Of course, ideally we would have a general theory of orthogonal polynomials in mathlib and then all proofs about eigenfunctions would be one liners. But don't let the perfect be the enemy of good ;) so I think it will be a great reference for someone else who might want to formalize some physics and don't know where to start.

[Joseph Tooby-Smith]

Thanks :).

Agreed. I've been using the 'Mathematics' file of PhysLean for things like this - where we know one day it will be in Mathlib but it would be nice to have certain specific results now to do something with .

[Joseph Tooby-Smith]

Ok, I managed to prove completeness of the eigenvectors for the harmonic oscillator - well more properly I've reduced it to a single theorem in analysis called Plancherel's theorem which is not yet in Mathlib.

The PR is here, it is not ready to push yet - as it needs a major tidy & refactor.