Mini projects

Stream: PhysLean

[Joseph Tooby-Smith]

For people interested in contributing to HepLean, but don't know where to start or what to do, I am going to use this thread to post mini project ideas. This is in addition to what HepLean currently has in the form of:

• A TODO list
• A graph of informal lemmas

which can also be used as a source of ideas of how to contribute.

Also if you have ideas for mini projects suitable for HepLean but you do not have the time or capacity to do them yourself please post them here.

[Joseph Tooby-Smith]

The two-higgs doublet model potential

Aim: Prove Theorem 1, 2 & 4 of this paper: https://arxiv.org/pdf/hep-ph/0605184, or variations of.

Physics context: The two-Higgs doublet model is a model of the universe with a proposed extra Higgs boson - giving two Higgs fields in total. Associated with these two Higgs fields is a potential whose minimum and boundness is very important to physics (it is what determines mass in the universe). The theorems in the above paper prove results related to these minima and boundness properties.

Math context: The actual maths involved is simple algebraic manipulation with potentially a small bit of calculus and group theory.

Current progress: Some of the set-up is already in HepLean but not much substantial progress has been made to proving these results.

Note: I don't believe you need to know much/any physics to do this project.

[Joseph Tooby-Smith]

Properties of grand unified theory groups

Aim: Prove group theoretic properties of the grand unified theories: SU(5), Spin(10) and Pati-Salam. Examples of such properties can be found in https://math.ucr.edu/home/baez/guts.pdf. They also appear as informal results in HepLean. See e.g. here.

Physics context: Grand unified theories are popular models beyond the standard model which unify all the fermions into the universe into the representation of a semi-simple or simple Lie group.

Math context: Group theory of several low-dimensional Lie groups (don't actually need Lie group structures) and maps between them.

Note: Like the project above, I don't believe you need to know much/any physics for this.

[Joseph Tooby-Smith]

Properties of Pauli matrices

Aim: Prove the properties of the Pauli matrices in Eq. 2.23 - 2.30 of this using index notation. An example of a similar calculation already in HepLean is this calculation.

[Joseph Tooby-Smith]

A start on optics

Aim: To convert the relevant parts of this into Lean.

[Lode Vermeulen]

Hi! I've been watching this project for a while now and I thought I'd have a try at one of the TODOs (semiformal_result (6YBYP): force). To begin with, I defined the force like this:

noncomputable def force (S : HarmonicOscillator) (x : Time → ℝ) : Time → ℝ := fun t =>
  - (deriv (potentialEnergy S x) (x t))

This gives me issues when trying to prove

lemma force_is_linear (x : Time → ℝ) : force S x = - S.k • x := by
  funext t
  unfold force potentialEnergy
  rw [deriv_const_mul_field, deriv_pow'']
  ring_nf

Namely, the proof state at the end here is -(S.k * x (x t) * deriv x (x t)) = (-S.k • x) t, which does not make sense to me (x (x t) seems like gibberish, other than that the two sides seem equal). It seems to me that a way around this issue is to refactor force and potentialEnergy to not be functions taking t, but I am not sure about this. Is there maybe some way to get around this issue of t being in the way in the proof itself, such that no refactoring is necessary?

[Joseph Tooby-Smith]

Hi Lode! Thanks for working on this :).

I personally can't think of a way of doing this without refactoring the potential energy and force to not depend on t. I think physically this might actually be the strictly correct thing to do in this case.

In particular I think the potential energy in general should be a function of current position of the particle (and possibly its current velocity - although not in this case), but it shouldn't be a function of the trajectory of the particle, so it probably needs refactoring in any case.

[ZhiKai Pong]

correct me if I'm wrong, but I think it might be a good idea to distinguish between the position x from the position function "x as a function of t" (not sure what's the best name, something like pos_fun might be too confusing). then you have a def saying x = pos_fun t.

(x should technically be of type Space 1 but I realized I can't operate directly on it e.g. taking the square, and I have to take the norm first which complicates things so maybe for the 1D case it's fine to just use ℝ?)

[Joseph Tooby-Smith]

@ZhiKai Pong you mean notationally wise? Maybe xₜ? But I also think sticking with x is ok too - it's what you see sometimes in the literature, albeit it's confusing.

Actually I think really one should use the inner-product here, since this is what would generalize to higher dimensions, namely:

/-- The potential energy of the harmonic oscillator is `1/2 k x ^ 2` -/
noncomputable def potentialEnergy (x : Space 1) :  ℝ := 1/2 * S.k * ⟪x, x⟫_ℝ

This will involve refactoring a lot of the equations in that file to replace ℝ with Space 1 though, but I think in the end it will be worth it.

So in short, @Lode Vermeulen , to do this TODO item I think you can either refactor potential energy to what I wrote above (which would be more work), or for now refactor it to:

/-- The potential energy of the harmonic oscillator is `1/2 k x ^ 2` -/
noncomputable def potentialEnergy (x : ℝ) :  ℝ := 1/2 * S.k * x ^ 2

and define force and force_is_linear from there. Either would be a very nice contribution :).

[ZhiKai Pong]

as in I think(x : ℝ) and (xₜ : Time → ℝ)are different things but both should be used in the definitions, where it is possible to write deriv xₜ tfor velocity etc., but allowing potentialEnergy to only depend on x

[Lode Vermeulen]

Thanks for your responses! I have now tried to refactor some of the definitions to use Space 1 (refactoring potentialEnergy to just ℝ was a success everywhere, but I would like to get the more general Space to work). It seems to mostly go ok, but I am having some issues again in the definition of force. I now have

noncomputable def potentialEnergy (x : Space 1) : ℝ :=
  1/2 * S.k * ⟪x, x⟫_ℝ

noncomputable def force {d : Fin 1} (S : HarmonicOscillator) (x : Space 1) : Space 1 :=
  - ∂[d] (potentialEnergy S x) x

It intuitively seems right that potentialEnergy has type ℝ and force has type Space 1 (should be a vector in the same space as the position). I am not entirely sure how to correctly state the force, as in every attempt there seems to have been an issue caused by the fact that potentialEnergy has type ℝ. I think physically it is fine for it to have that, as the derivative of a scalar function with respect to a vector (in this case the 1D x) should still return a vector (which force expects). Note that I have tried many different ways of defining the derivative (deriv, fderiv, Space.deriv, gradient and ∂[d]).

[ZhiKai Pong]

@Lode Vermeulen good point, I think the correct definition here should actually be using Space.grad (which gives you a vector) instead of deriv.

something like this:

open InnerProductSpace
open Space
/-- The potential energy of the harmonic oscillator is `1/2 k x ^ 2` -/
noncomputable def potentialEnergy (x : Space 1) : ℝ :=
  1/2 * S.k * ⟪x, x⟫_ℝ

noncomputable def force (S : HarmonicOscillator) (x : Space 1) : EuclideanSpace ℝ (Fin 1) :=
  - ∇ (potentialEnergy S) x

[ZhiKai Pong]

Space should really only be for spatial coordinates - although Space 1 and EuclideanSpace ℝ (Fin 1) is exactly the same in terms of mathematical behaviour, I believe the rationale is that having this distinction makes things clearer, so force here being a vector should have type EuclideanSpace ℝ (Fin 1). It is also the reason we defined our custom Space.grad etc. for working with functions in physical Space.

[Lode Vermeulen]

Thanks, it all seems to work well now. I have been stuck at trying to prove force_is_linear for some time now. I have been able to get the goal to -∇ (1 / 2 • S.k • fun x => ⟪x, x⟫_ℝ) x = -S.k • x. My plan was to apply a new helper lemma that I added in VectorIdentities.lean:

@[simp]
lemma grad_const_smul (f : Space d → ℝ) (hf : Differentiable ℝ f) :
    ∇ (c • f) = c • ∇ f := by
  unfold grad Space.deriv
  ext x i
  rw [fderiv_const_smul']
  simp only [ContinuousLinearMap.coe_smul', Pi.smul_apply, smul_eq_mul]
  exact hf.differentiableAt

But I just don't seem to be able to apply this lemma, for some reason. For context, I defined have hd : Differentiable ℝ (fun (x : Space 1) => ⟪x, x⟫_ℝ) (which I have just sorry'd in the proof for now), which is necessary as hf in grad_const_mul. Could someone give me some pointers on how to resolve this?

[Joseph Tooby-Smith]

(deleted: sorry I missed the 1 / 2 - length of my window was just the wrong length)

[Joseph Tooby-Smith]

Using erw [grad_const_smul] works for me, but it isn't ideal

[Joseph Tooby-Smith]

Ah, so 1 / 2 • S.k • fun x => ⟪x, x⟫_ℝ is not (1/2) • S.k • fun x => ⟪x, x⟫_ℝ it's 1/ (2 • S.k • fun x => ⟪x, x⟫_ℝ). I'm guessing you actually want the former.

[Joseph Tooby-Smith]

Actually one needs to be even more explicit: (1 / (2 : ℝ)) • S.k • fun x => ⟪x, x⟫_ℝ to get things to work, otherwise it takes (1/2) to be division in the natural numbers :(.

[Joseph Tooby-Smith]

Here are a couple of results which could be added to VectorIdentities to help:

Code

@[simp]
lemma deriv_component_same (μ : Fin d) (x : Space d) :
    (deriv μ (fun x => x μ) x) = 1 := by
  simp [deriv, fderiv_coord_eq_proj_comp, ContinuousLinearMap.proj]
  erw [LinearMap.proj_apply]
  simp

lemma deriv_component_diff (μ ν : Fin d) (x : Space d) (h : μ ≠ ν) :
    (deriv μ (fun x => x ν) x) = 0 := by
  simp [deriv, fderiv_coord_eq_proj_comp, ContinuousLinearMap.proj]
  erw [LinearMap.proj_apply]
  simp [h]
  omega

open InnerProductSpace

lemma inner_differentiable {d : ℕ} :
    Differentiable ℝ (fun y : Space d => ⟪y, y⟫_ℝ) := by
  simp
  fun_prop

lemma grad_inner {d : ℕ} :
    ∇ (fun y : Space d => ⟪y, y⟫_ℝ) = fun z => 2 • z := by
  ext z i
  simp [Space.grad]
  rw [deriv]
  rw [fderiv_sum]
  · simp
    rw [Finset.sum_eq_single i]
    · trans (fderiv ℝ (fun y => y i ^ 2) z) (EuclideanSpace.single i 1)
      · congr
        funext y
        ring
      trans deriv i ((fun x => x^ 2) ∘ fun y => y i ) z
      · rfl
      rw [deriv, fderiv_comp]
      · simp
        rw [← deriv_eq]
        simp
      · fun_prop
      · fun_prop
    · intro b _ hb
      trans (fderiv ℝ (fun y => y b ^ 2) z) (EuclideanSpace.single i 1)
      · congr
        funext y
        ring
      trans deriv i ((fun x => x^ 2) ∘ fun y => y b) z
      · rfl
      rw [deriv, fderiv_comp]
      simp
      · left
        rw [← deriv_eq]
        rw [deriv_component_diff]
        omega
      · fun_prop
      · fun_prop
    · simp
  · intro i
    fun_prop

(Need to import import Mathlib.Analysis.SpecialFunctions.Pow.Deriv)

[ZhiKai Pong]

for some reason there's no fderiv_pow, should be simpler if we just define that instead

[ZhiKai Pong]

I imagine it can probably be defined in the exact same way as e.g. docs#fderiv_cos

[ZhiKai Pong]

might still be a good idea to have grad_inner as you defined here, I'm just saying the proof can probably be much simpler. maybe I can look into PRing fderiv_pow to mathlib when I have time.

[Joseph Tooby-Smith]

One should actually be able to use docs#fderiv_inner_apply in grad_inner to simplify it aswell.

[ZhiKai Pong]

here's how I'll write this for now:

Code

lemma grad_const_smul (c : ℝ) (f : Space d → ℝ) (hf : Differentiable ℝ f) :
    ∇ (c • f) = c • ∇ f := by
  unfold grad
  ext x i
  rw [deriv_smul]
  simp only [Pi.smul_apply, smul_eq_mul]
  fun_prop

lemma force_is_linear (x : Space 1) :
    -∇ ((1 / (2:ℝ)) • S.k • fun x => ⟪x, x⟫_ℝ) x = -S.k • x := by
  rw [grad_const_smul, grad_const_smul]
  unfold grad Space.deriv
  ext i
  simp [-PiLp.inner_apply]
  rw [fderiv_inner_apply]
  fin_cases i
  simp
  norm_num
  ring
  repeat fun_prop
  · apply Differentiable.inner
    repeat fun_prop
  · change Differentiable ℝ (fun x => S.k • ⟪x, x⟫_ℝ)
    apply Differentiable.const_smul
    apply Differentiable.inner
    repeat fun_prop

but feel free to implement some of the suggestions above which helps with the API and leads to less unfolding in the future :)

[Joseph Tooby-Smith]

btw I think grad_const_smul already exists as grad_smul in VectorIdentities

[ZhiKai Pong]

(! I completely forgot about that lol)

[Lode Vermeulen]

Thanks for these!

Those conflicts between ℕ and ℝ are really sneaky.

Joseph Tooby-Smith said:

btw I think grad_const_smul already exists as grad_smul in VectorIdentities

Completely read over this as well.

[Lode Vermeulen]

Could I get push permission? I'd like to create a PR with this code.

[Lode Vermeulen]

These changes have caused some issues in/HarmonicOscillator/Solutions.lean, as Space 1 is now used instead of ℝ. Mostly these are pretty fixable. I am having trouble with the following:

open Complex in
lemma amplitude_eq_norm (IC : InitialConditions) :
    S.amplitude IC = ‖(↑IC.x₀ + -↑IC.v₀ / ↑S.ω * Complex.I)‖ := by
  rw [amplitude_eq]
  trans √(IC.x₀ ^ 2 + (- IC.v₀ / S.ω) ^ 2)
    · ring_nf
    · simp [← Complex.norm_add_mul_I]

Which I have translated to

structure InitialConditions where
  /-- The initial position of the harmonic oscillator. -/
  x₀ :  Space 1
  /-- The initial velocity of the harmonic oscillator. -/
  v₀ : Space 1

open Complex in
lemma amplitude_eq_norm (IC : InitialConditions) :
    S.amplitude IC = ‖(((↑IC.x₀) : (EuclideanSpace ℂ (Fin 1))) - ((1:ℂ) / ↑S.ω) • Complex.I • (↑IC.v₀ : EuclideanSpace ℂ (Fin 1)))‖ := by

Using the coercion defined in /Space/Basic.lean:

instance : Coe (EuclideanSpace ℝ (Fin d)) (EuclideanSpace ℂ (Fin d)) where
  coe v := fun i => (v i : ℂ)

Which I got using some help from Gemini. This does allow me to do the coercion, i.e. the statement of amplitude_eq_norm does not give me any errors. Now I still need to do the proof, which now includes the coercion (so specifically, coe v := fun i => (v i : ℂ)). Before trying to prove it, I first want to figure out whether the statements are even equal. I don't understand coercions and this instance statement well enough, so I hoped someone could let me know if this coercion is even correct in the first place.

EDIT: The proof state (after rw [amplitude_eq]) is the following:

S : HarmonicOscillator
IC : InitialConditions
⊢ √(⟪IC.x₀, IC.x₀⟫_ℝ + ⟪(1 / S.ω) • IC.v₀, (1 / S.ω) • IC.v₀⟫_ℝ) =
  ‖(fun i => ↑(IC.x₀ i)) - (1 / ↑S.ω) • I • fun i => ↑(IC.v₀ i)‖

[Joseph Tooby-Smith]

On a train, so only a partial answer. But here is a suggestion that might make things easier. Instead of the coercion you have defined, note the following: IC.x₀ is a vector with one element, and IC.x₀ 0 : ℝ is that element. So if you change the theorem to

lemma amplitude_eq_norm (IC : InitialConditions) :
    S.amplitude IC = ‖(((IC.x₀ 0)  - ((1:ℂ) / ↑S.ω) • Complex.I • (IC.v₀ 0))‖ := by

the coercion between ℝ and ℂ should take care of things. This might require changing instances of IC.x₀ appearing earlier in the file to IC.x₀ 0, etc.

[Lode Vermeulen]

Thanks, that seems to work.

[Joseph Tooby-Smith]

Instead of this thread for mini project ideas, I've updated the website to include this page:

https://physlean.com/ProjectIdeas

If anyone has any suggestions or ideas of what they would like to see there, let me know.

[Lode Vermeulen]

I am running into some more problems now.

1. With the above definition of amplitude_eq_norm, I am unsure how to prove identities like ‖IC.x₀‖ = IC.x₀ 0.
2. I am unsure how to correctly state the following lemma:

lemma sol_eq_amplitude_mul_cos_phase (IC : InitialConditions) :
    S.sol IC = fun t => S.amplitude IC * cos (S.ω * t + S.phase IC)

As the LHS has type Time → Space 1 and the RHS Time → ℝ . The way to fix this is to transform the type of amplitude to Space 1, but this is not intuitive (it should be ℝ) nor clear how to do so.

3. potentialEnergy is now of type ℝ, but sol_potentialEnergy is of type Time → ℝ (which also does not make sense to just remove as S.sol IC also takes time as an input). Especially considering the following, I guess I am a bit stuck with this
   Joseph Tooby-Smith said:

I personally can't think of a way of doing this without refactoring the potential energy and force to not depend on t.

[Joseph Tooby-Smith]

Regarding your points:

1. I don't think the identity ‖IC.x₀‖ = IC.x₀ 0 should be true, what I think is true is ‖IC.x₀‖ = |IC.x₀ 0|. Is there a context in which you need ‖IC.x₀‖ = IC.x₀ 0?
2. I would rewrite this as something like:

lemma sol_eq_amplitude_mul_cos_phase (IC : InitialConditions) :
    S.sol IC = fun t => S.amplitude IC • (fun _ =>  cos (S.ω * t + S.phase IC))

That is, trivially make the RHS a vector.

3. I think you need to replace sol_potentialEnergy with:

lemma sol_potentialEnergy (IC : InitialConditions) (t : Time) :  S.potentialEnergy (S.sol IC t) = 1/2 * (S.k * ‖IC.x₀‖ ^ 2 + S.m * ‖IC.v₀‖ ^2) * cos (S.ω * t + S.phase IC) ^ 2 := by ...

An aside: Ultimately I think we should repalce the lagrangian to be a function of x : Space 1 and a velocityv : ℝ (or a type we make for velocity), so that it will work with the Euler-Lagrange equations and variational calculus we have. It partially for this reason we want the potential to be a function of just the position not the trajectory.

[Lode Vermeulen]

Thanks! For 1., you're right, I actually need ‖IC.x₀‖ ^ 2=IC.x₀ 0 ^ 2 (which should be true) instead in

lemma amplitude_eq_norm (IC : InitialConditions) :
    S.amplitude IC = ‖((IC.x₀ 0)  - ((1:ℂ) / ↑S.ω) • (IC.v₀ 0) • Complex.I)‖ := by
  rw [amplitude_eq]
  trans √(‖IC.x₀‖^2 + (‖IC.v₀‖/S.ω)^2)
  · ring
  · simp only [Complex.norm_add_mul_I, norm_eq_sqrt_sq_add_sq]
    simp only [Fin.isValue, one_div, real_smul, smul_eq_mul, sub_re, ofReal_re, mul_re, inv_re,
      normSq_ofReal, div_self_mul_self', I_re, mul_zero, ofReal_im, I_im, mul_one, sub_self, inv_im,
      neg_zero, zero_div, mul_im, add_zero, zero_mul, sub_zero, sub_im, zero_sub, even_two,
      Even.neg_pow]

[Lode Vermeulen]

I still don't know how to prove that fact though :sweat_smile:

[Joseph Tooby-Smith]

Do you mean ‖IC.x₀‖ ^ 2=IC.x₀ 0 ^ 2? if so the below is an equivalent lemma.

import Mathlib

open InnerProductSpace
example (x : EuclideanSpace ℝ (Fin 1)) :
    ‖x‖ ^ 2 = (x 0) ^2 := by
  rw [@PiLp.norm_sq_eq_of_L2] -- found using rw?
  simp

[ZhiKai Pong]

Great work for finishing this! I have a few minor comments/questions, but the PR has been merged already and maybe @Joseph Tooby-Smith can answer those from a design philosophy perspective so I thought maybe I'll just ask them here.

1. Mathlib.Algebra.Lie.OfAssociative, Mathlib.Analysis.CStarAlgebra.Classes, Mathlib.Data.Real.StarOrdered can be removed without breaking anything (I'm also curious why they're there in the first place). I thought after PR#570 the linter checks for redundant imports, maybe I'm misunderstanding how it works. (#min_imports only catches the first one I think, is it because the file contains sorry?)

2. polarCoord is used to define amplitude and phase which feels a bit weird to me since these quantities have nothing to do with coordinates, personally I'd rather define it via the norm and Complex.Arg directly without invoking polarCoord.

3. I feel like the overall exposition right now is a bit strange, it looks like "given a system with these parameters, define arbitrary quantities and relations, then show everything works out" but doesn't really feel like deriving properties of a harmonic oscillator. Maybe this is a byproduct of the Euler-Lagrange equations not ready yet and the overall presentation can be refactored afterwards (but I also think another starting point would be Newton's second law i.e. F = ma = -kx which isn't really used here). To give a concrete example, let's consider the amplitude. The way I would think about the amplitude is the maximal displacement, which is the extrema of the position. I would expect some notion of "maximum" in the definition, then show that it equals the norm of initial condition (probably involving differentiation). Not saying any of the current results are wrong, but the current exposition just doesn't seem that "physics" to me. Maybe I'm completely overthinking this but would like to know what you think :)

[Joseph Tooby-Smith]

Hey @ZhiKai Pong , let me address your comments in turn, all of which I agree with.

1. The redundant imports linter only checks for transitive imports (i.e. when an import is included already via another file), it is not quite the same as #min_imports, which is harder to implement a linter for because of the number of false-positives. I'm guessing these imports were added via a copy-and-paste (probably by me), but should be removed. Really every so often we should run:

lake exe shake --cfg ./scripts/physleanNoShake.json PhysLean

and clean up the results.

2. I'm happy for these to be changed, really if we can define things without needing complex numbers at all, I think that would be better.

3. Indeed, I think the most important TODO items here are still left: that is showing that the solutions actually solve the equation of motion and that they are unique given their boundary conditions. I agree with your comment about the amplitude, and I think it should be defined as the maximum displacement, and more generally that the overall structure of the file can be made better. This is my bad.

[Joseph Tooby-Smith]

Surely the phase also has a physical definition too, based on the time t = 0 and the how long it takes to get from t=0 to the amplitude.

[Lode Vermeulen]

I'd love to help on any tasks that come out of this!

I think these are some good points. Indeed implementing the Euler-Lagrange equations would probably help make everything more physics-y. Furthermore, I was wondering about extending these results to n dimensions to make it as general as possible. At first glance it seems like this would not be too hard if we just use the norm in calculations everywhere. It might also be interesting to introduce dampening/driving coefficients. I feel like those things would be nice if this would be used as an API, right?

[Joseph Tooby-Smith]

Yeah, I agree, any of these seem useful. For the Euler-Lagrange equations that starting point should likely be eulerLagrangeOp, the statement that this is an extrema of the action is just the lemma euler_lagrange_varGradient. (see this file).

I think physically what makes the most sense is the following:

1. Obtain the Euler-Lagrange equations from eulerLagrangeOp.
2. Show that for a given set of initial conditions, the Euler lagrange equations have a unique solution. (not sure how difficult this is - but once we know how to do this in this example, it will be very useful elsewhere.)
3. Using choice, define the sol to be said unique solution.
4. Then show that sol has the form we have.

This kind of follows ZhiKai's suggested approach with the amplitude.

[Joseph Tooby-Smith]

Uniqueness can probably be proved similar to docs#ODE_solution_unique_univ

[Joseph Tooby-Smith]

I added another project to: https://physlean.com/ProjectIdeas, which is actually similar in nature to the one discussed here regarding the harmonic oscillator. That is: Formalize the classical mechanics of a pendulum. If anyone wants to have a go at this, that would be great - should be a nice 'relatively' simple project (although I note even simple looking projects can turn out difficult :)).