diff --git a/PhysLean.lean b/PhysLean.lean
index 71682c2db..621f78fac 100644
--- a/PhysLean.lean
+++ b/PhysLean.lean
@@ -18,9 +18,19 @@ import PhysLean.Electromagnetism.Homogeneous
 import PhysLean.Electromagnetism.LorentzAction
 import PhysLean.Electromagnetism.MaxwellEquations
 import PhysLean.Electromagnetism.Wave
+import PhysLean.Mathematics.LinearAlgebra.BilinearForm
+import PhysLean.Mathematics.Analysis.ContDiff
 import PhysLean.Mathematics.FDerivCurry
 import PhysLean.Mathematics.Fin
 import PhysLean.Mathematics.Fin.Involutions
+import PhysLean.Mathematics.Geometry.Manifold.PartialHomeomorph
+import PhysLean.Mathematics.Geometry.Manifold.Chart.Utilities
+import PhysLean.Mathematics.Geometry.Manifold.Chart.BilinearSmoothness
+import PhysLean.Mathematics.Geometry.Manifold.Chart.Smoothness
+import PhysLean.Mathematics.Geometry.Manifold.Chart.CoordinateTransformations
+import PhysLean.Mathematics.Geometry.Metric.PseudoRiemannian.Defs
+import PhysLean.Mathematics.Geometry.Metric.PseudoRiemannian.Chart
+import PhysLean.Mathematics.Geometry.Metric.Riemannian.Defs
 import PhysLean.Mathematics.LinearMaps
 import PhysLean.Mathematics.List
 import PhysLean.Mathematics.List.InsertIdx
diff --git a/PhysLean/Mathematics/Analysis/ContDiff.lean b/PhysLean/Mathematics/Analysis/ContDiff.lean
new file mode 100644
index 000000000..f4438c13f
--- /dev/null
+++ b/PhysLean/Mathematics/Analysis/ContDiff.lean
@@ -0,0 +1,230 @@
+/-
+Copyright (c) 2025 Matteo Cipollina. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Matteo Cipollina
+-/
+
+import Mathlib.Analysis.Calculus.ContDiff.Operations
+import Mathlib.LinearAlgebra.Dual.Lemmas
+import Mathlib.LinearAlgebra.FreeModule.PID
+import Mathlib.RingTheory.Henselian
+import PhysLean.Mathematics.LinearAlgebra.BilinearForm
+
+/-!
+# Smoothness (ContDiff) Utilities
+
+This file provides utility lemmas and constructions for working with smooth
+functions (`ContDiff`) and continuity in the context of normed and finite-dimensional
+vector spaces over a nontrivially normed field.
+
+-/
+
+namespace ContDiff
+
+variable {𝕜 X V : Type*} [NontriviallyNormedField 𝕜]
+variable [NormedAddCommGroup X] [NormedSpace 𝕜 X]
+variable [NormedAddCommGroup V] [NormedSpace 𝕜 V]
+variable {f : X → V} {s : Set X} {x₀ : X} {n : WithTop ℕ∞}
+
+-- First direction: if f is C^n, then φ ∘ f is C^n for any continuous linear functional φ
+lemma comp_continuous_linear_apply_right
+    (hf : ContDiffWithinAt 𝕜 n f s x₀) (φ : V →L[𝕜] 𝕜) :
+    ContDiffWithinAt 𝕜 n (φ ∘ f) s x₀ :=
+  ContDiffWithinAt.comp x₀ φ.contDiff.contDiffWithinAt hf (Set.mapsTo_univ _ _)
+
+-- Second direction: in finite dimensions, if all projections are C^n, then f is C^n
+lemma of_forall_coord [FiniteDimensional 𝕜 V] [CompleteSpace 𝕜]
+    (h : ∀ φ : V →L[𝕜] 𝕜, ContDiffWithinAt 𝕜 n (φ ∘ f) s x₀) :
+    ContDiffWithinAt 𝕜 n f s x₀ := by
+  let b := Module.finBasis 𝕜 V
+  let equiv := b.equivFunL
+  suffices ContDiffWithinAt 𝕜 n (equiv ∘ f) s x₀ by
+    have hequiv_symm_smooth : ContDiff 𝕜 ⊤ equiv.symm := ContinuousLinearEquiv.contDiff equiv.symm
+    have hequiv_symm_smooth_n : ContDiff 𝕜 n equiv.symm :=
+      ContDiff.of_le hequiv_symm_smooth (le_top : n ≤ ⊤)
+    have h_eq : f = equiv.symm ∘ (equiv ∘ f) := by
+      ext x; simp only [Function.comp_apply, ContinuousLinearEquiv.symm_apply_apply];
+    rw [h_eq]
+    apply ContDiffWithinAt.comp x₀ hequiv_symm_smooth_n.contDiffWithinAt this (Set.mapsTo_univ _ _)
+  apply contDiffWithinAt_pi.mpr
+  intro i
+  let coord_i : V →L[𝕜] 𝕜 := LinearMap.toContinuousLinearMap (b.coord i)
+  exact h coord_i
+
+-- Full bidirectional lemma
+lemma iff_forall_coord [FiniteDimensional 𝕜 V] [CompleteSpace 𝕜] :
+    ContDiffWithinAt 𝕜 n f s x₀ ↔
+    ∀ φ : V →L[𝕜] 𝕜, ContDiffWithinAt 𝕜 n (φ ∘ f) s x₀ := by
+  constructor
+  · exact comp_continuous_linear_apply_right
+  · exact of_forall_coord
+
+end ContDiff
+
+section ContinuityBounds
+
+variable {𝕜 E F FHom : Type*} [NormedField 𝕜]
+
+@[simp]
+lemma AddMonoidHomClass.coe_fn_to_addMonoidHom
+    [FunLike FHom E F] [AddZeroClass E] [AddZeroClass F]
+    [AddMonoidHomClass FHom E F] (φ : FHom) :
+    ⇑(AddMonoidHomClass.toAddMonoidHom φ) = ⇑φ := by
+  rfl
+
+variable [NormedAddCommGroup E] [NormedSpace 𝕜 E]
+variable [NormedAddCommGroup F] [NormedSpace 𝕜 F]
+
+/-- A bounded additive map is continuous at zero. -/
+lemma AddMonoidHom.continuousAt_zero_of_bound
+    (φ : AddMonoidHom E F) {C : ℝ} (h : ∀ x, ‖φ x‖ ≤ C * ‖x‖) :
+    ContinuousAt φ 0 := by
+  rw [Metric.continuousAt_iff]
+  intro ε εpos
+  simp only [map_zero φ, dist_zero_right]
+  by_cases hE : Subsingleton E
+  · use 1
+    refine ⟨zero_lt_one, fun y _hy_norm_lt_one => ?_⟩
+    rw [@Subsingleton.elim E hE y 0, map_zero φ, norm_zero]
+    exact εpos
+  · have C_nonneg : 0 ≤ C := by
+      obtain ⟨x_ne, y_ne, h_x_ne_y⟩ : ∃ x y : E, x ≠ y := by
+        contrapose! hE; exact { allEq := hE }
+      let z := x_ne - y_ne
+      have hz_ne_zero : z ≠ 0 := sub_ne_zero_of_ne h_x_ne_y
+      have hz_norm_pos : 0 < ‖z‖ := norm_pos_iff.mpr hz_ne_zero
+      by_contra hC_is_neg
+      push_neg at hC_is_neg
+      have h_phi_z_bound := h z
+      have H1 : 0 ≤ C * ‖z‖ := le_trans (norm_nonneg (φ z)) h_phi_z_bound
+      have H2 : C * ‖z‖ < 0 := mul_neg_of_neg_of_pos hC_is_neg hz_norm_pos
+      linarith [H1, H2]
+    by_cases hC_eq_zero : C = 0
+    · have phi_is_zero : φ = 0 := by
+        ext x_val
+        have h_phi_x_val_bound := h x_val
+        rw [hC_eq_zero, zero_mul] at h_phi_x_val_bound
+        exact norm_le_zero_iff.mp h_phi_x_val_bound
+      use 1
+      refine ⟨zero_lt_one, fun y _hy_norm_lt_one => ?_⟩
+      rw [phi_is_zero, AddMonoidHom.zero_apply, norm_zero]
+      exact εpos
+    · have C_pos : 0 < C := lt_of_le_of_ne C_nonneg fun a => hC_eq_zero (_root_.id (Eq.symm a))
+      use ε / C
+      refine ⟨div_pos εpos C_pos, fun y hy_norm_lt_delta => ?_⟩
+      calc
+        ‖φ y‖ ≤ C * ‖y‖ := h y
+        _ < C * (ε / C) := mul_lt_mul_of_pos_left hy_norm_lt_delta C_pos
+        _ = ε := by rw [mul_div_cancel₀ ε hC_eq_zero]
+
+omit [NormedSpace 𝕜 F] in
+/-- A semi-linear map that is linearly bounded by the norm of its input is continuous. -/
+lemma SemilinearMapClass.continuous_of_bound {𝕜₂ : Type*} [NormedField 𝕜₂] [NormedSpace 𝕜₂ F]
+    [FunLike FHom E F] {σ : 𝕜 →+* 𝕜₂} [SemilinearMapClass FHom σ E F]
+    {φ : FHom} {C : ℝ} (h : ∀ x, ‖φ x‖ ≤ C * ‖x‖) : Continuous φ := by
+  haveI : AddMonoidHomClass FHom E F := inferInstance
+  let φ_add_hom : AddMonoidHom E F := AddMonoidHomClass.toAddMonoidHom φ
+  exact continuous_of_continuousAt_zero φ_add_hom
+    (AddMonoidHom.continuousAt_zero_of_bound φ_add_hom h)
+
+/-- A function that is linearly bounded by the norm of its input is continuous. -/
+lemma AddMonoidHomClass.continuous_of_bound' [FunLike FHom E F] [AddMonoidHomClass FHom E F]
+    {φ : FHom} {C : ℝ} (h : ∀ x, ‖φ x‖ ≤ C * ‖x‖) : Continuous φ := by
+  let φ_add_hom : AddMonoidHom E F := AddMonoidHomClass.toAddMonoidHom φ
+  exact continuous_of_continuousAt_zero φ_add_hom
+    (AddMonoidHom.continuousAt_zero_of_bound φ_add_hom h)
+
+end ContinuityBounds
+
+namespace ContinuousLinearMap
+
+variable {X₁ E₁ F₁ G₁ E₁' F₁' : Type*} [NontriviallyNormedField 𝕜₁]
+  [NormedAddCommGroup X₁] [NormedSpace 𝕜₁ X₁]
+  [NormedAddCommGroup E₁] [NormedSpace 𝕜₁ E₁]
+  [NormedAddCommGroup F₁] [NormedSpace 𝕜₁ F₁]
+  [NormedAddCommGroup G₁] [NormedSpace 𝕜₁ G₁]
+  [NormedAddCommGroup E₁'] [NormedSpace 𝕜₁ E₁']
+  [NormedAddCommGroup F₁'] [NormedSpace 𝕜₁ F₁']
+  {n₁ : WithTop ℕ∞}
+
+/-- The `ContinuousLinearMap.bilinearComp` operation is smooth.
+    Given smooth functions `f : X₁ → (E₁ →L[𝕜₁] F₁ →L[𝕜₁] G₁)`, `g : X₁ → (E₁' →L[𝕜₁] E₁)`,
+    and `h : X₁ → (F₁' →L[𝕜₁] F₁)`, the composition `x ↦ (f x).bilinearComp (g x) (h x)`
+    is smooth. -/
+lemma contDiff_bilinearComp
+    {f : X₁ → E₁ →L[𝕜₁] F₁ →L[𝕜₁] G₁} {g : X₁ → E₁' →L[𝕜₁] E₁} {h : X₁ → F₁' →L[𝕜₁] F₁}
+    (hf : ContDiff 𝕜₁ n₁ f) (hg : ContDiff 𝕜₁ n₁ g) (hh : ContDiff 𝕜₁ n₁ h) :
+    ContDiff 𝕜₁ n₁ fun x => (f x).bilinearComp (g x) (h x) := by
+  have h1 : ContDiff 𝕜₁ n₁ (fun x ↦ (f x).comp (g x)) := ContDiff.clm_comp hf hg
+  let L_flip1 := ContinuousLinearMap.flipₗᵢ 𝕜₁ E₁' F₁ G₁
+  have eq_flip : ∀ x, L_flip1 ((f x).comp (g x)) = ((f x).comp (g x)).flip := by
+    intro x
+    rfl
+  have h2 : ContDiff 𝕜₁ n₁ (fun x => ((f x).comp (g x)).flip) := by
+    have hL₁ : ContDiff 𝕜₁ n₁ L_flip1 :=
+      (ContinuousLinearMap.contDiff (𝕜 := 𝕜₁) (E := E₁' →L[𝕜₁] F₁ →L[𝕜₁] G₁)
+        (F := F₁ →L[𝕜₁] E₁' →L[𝕜₁] G₁) L_flip1).of_le le_top
+    have h2' : ContDiff 𝕜₁ n₁ (fun x => L_flip1 ((f x).comp (g x))) :=
+      ContDiff.comp hL₁ h1
+    exact (funext eq_flip).symm ▸ h2'
+  have h3 : ContDiff 𝕜₁ n₁ (fun x => (((f x).comp (g x)).flip).comp (h x)) :=
+    ContDiff.clm_comp h2 hh
+  let L_flip2 := ContinuousLinearMap.flipₗᵢ 𝕜₁ F₁' E₁' G₁
+  have eq_flip2 : ∀ x, L_flip2 ((((f x).comp (g x)).flip).comp (h x)) =
+      ((((f x).comp (g x)).flip).comp (h x)).flip := by
+    intro x
+    rfl
+  have h4 : ContDiff 𝕜₁ n₁ (fun x => ((((f x).comp (g x)).flip).comp (h x)).flip) := by
+    have hL₂ : ContDiff 𝕜₁ n₁ L_flip2 :=
+      (ContinuousLinearMap.contDiff (𝕜 := 𝕜₁) (E := F₁' →L[𝕜₁] E₁' →L[𝕜₁] G₁)
+        (F := E₁' →L[𝕜₁] F₁' →L[𝕜₁] G₁) L_flip2).of_le le_top
+    have h4' := ContDiff.comp hL₂ h3
+    exact (funext eq_flip2).symm ▸ h4'
+  exact h4
+
+variable {X₁ E₁ F₁ G₁ E₁' F₁' : Type*} [NontriviallyNormedField 𝕜₁]
+  [NormedAddCommGroup X₁] [NormedSpace 𝕜₁ X₁]
+  [NormedAddCommGroup E₁] [NormedSpace 𝕜₁ E₁] [FiniteDimensional 𝕜₁ E₁]
+  [NormedAddCommGroup F₁] [NormedSpace 𝕜₁ F₁] [FiniteDimensional 𝕜₁ F₁]
+  [NormedAddCommGroup G₁] [NormedSpace 𝕜₁ G₁] [FiniteDimensional 𝕜₁ G₁]
+  [NormedAddCommGroup E₁'] [NormedSpace 𝕜₁ E₁'] [FiniteDimensional 𝕜₁ E₁']
+  [NormedAddCommGroup F₁'] [NormedSpace 𝕜₁ F₁'] [FiniteDimensional 𝕜₁ F₁']
+  {n₁ : WithTop ℕ∞}
+
+/-- The "flip" operation on continuous bilinear maps is smooth. -/
+lemma flip_contDiff {F₁ F₂ R : Type*}
+    [NormedAddCommGroup F₁] [NormedSpace ℝ F₁]
+    [NormedAddCommGroup F₂] [NormedSpace ℝ F₂]
+    [NormedAddCommGroup R] [NormedSpace ℝ R] :
+    ContDiff ℝ ⊤ (fun f : F₁ →L[ℝ] F₂ →L[ℝ] R => ContinuousLinearMap.flip f) := by
+  let flip_clm :=
+    (ContinuousLinearMap.flipₗᵢ ℝ F₁ F₂ R).toContinuousLinearEquiv.toContinuousLinearMap
+  exact
+    @ContinuousLinearMap.contDiff ℝ _
+      (F₁ →L[ℝ] F₂ →L[ℝ] R) _ _ (F₂ →L[ℝ] F₁ →L[ℝ] R) _ _ _ flip_clm
+
+/-- Composition of a bilinear map with a linear map in the first argument is smooth. -/
+lemma comp_first_contDiff {F₁ F₂ F₃ R : Type*}
+    [NormedAddCommGroup F₁] [NormedSpace ℝ F₁]
+    [NormedAddCommGroup F₂] [NormedSpace ℝ F₂]
+    [NormedAddCommGroup F₃] [NormedSpace ℝ F₃]
+    [NormedAddCommGroup R] [NormedSpace ℝ R] :
+    ContDiff ℝ ⊤ (fun p : (F₂ →L[ℝ] F₃ →L[ℝ] R) × (F₁ →L[ℝ] F₂) =>
+      ContinuousLinearMap.comp p.1 p.2) := by
+  exact ContDiff.clm_comp contDiff_fst contDiff_snd
+
+variable {E₁_₂ : Type*} {E₂_₂ : Type*} {R₂ : Type*}
+variable [NormedAddCommGroup E₁_₂] [NormedSpace ℝ E₁_₂]
+variable [NormedAddCommGroup E₂_₂] [NormedSpace ℝ E₂_₂]
+variable [NormedAddCommGroup R₂] [NormedSpace ℝ R₂]
+
+/-- The pullback of a bilinear map by a linear map is smooth with respect to both arguments. -/
+theorem contDiff_pullbackBilinear_op :
+    ContDiff ℝ ⊤ (fun p : (E₂_₂ →L[ℝ] E₂_₂ →L[ℝ] R₂) × (E₁_₂ →L[ℝ] E₂_₂) =>
+      BilinearForm.pullback p.1 p.2) := by
+  apply contDiff_bilinearComp
+  · exact (contDiff_fst (E := (E₂_₂ →L[ℝ] E₂_₂ →L[ℝ] R₂)) (F := (E₁_₂ →L[ℝ] E₂_₂))).of_le le_top
+  · exact (contDiff_snd (E := (E₂_₂ →L[ℝ] E₂_₂ →L[ℝ] R₂)) (F := (E₁_₂ →L[ℝ] E₂_₂))).of_le le_top
+  · exact (contDiff_snd (E := (E₂_₂ →L[ℝ] E₂_₂ →L[ℝ] R₂)) (F := (E₁_₂ →L[ℝ] E₂_₂))).of_le le_top
+
+end ContinuousLinearMap
diff --git a/PhysLean/Mathematics/Geometry/Manifold/Chart/BilinearSmoothness.lean b/PhysLean/Mathematics/Geometry/Manifold/Chart/BilinearSmoothness.lean
new file mode 100644
index 000000000..a5fd7fda9
--- /dev/null
+++ b/PhysLean/Mathematics/Geometry/Manifold/Chart/BilinearSmoothness.lean
@@ -0,0 +1,135 @@
+/-
+Copyright (c) 2025 Matteo Cipollina. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Matteo Cipollina
+-/
+
+import Mathlib.Algebra.Lie.OfAssociative
+import Mathlib.Analysis.RCLike.Lemmas
+import Mathlib.Geometry.Manifold.IsManifold.Basic
+import PhysLean.Mathematics.Analysis.ContDiff
+
+/-!
+# Smoothness of Bilinear Forms in Chart Coordinates
+
+This file contains lemmas about the smoothness of bilinear forms in chart coordinates.
+-/
+noncomputable section
+open BilinearForm
+open Filter
+
+variable {E : Type v} {M : Type v} {n : WithTop ℕ∞}
+variable [NormedAddCommGroup E] [NormedSpace ℝ E]
+variable [TopologicalSpace M] [ChartedSpace E M] [T2Space M]
+variable {I : ModelWithCorners ℝ E E} [ModelWithCorners.Boundaryless I]
+variable [inst_mani_smooth : IsManifold I (n + 1) M] -- For C^{n+1} manifold for C^n metric
+variable [inst_tangent_findim : ∀ (x : M), FiniteDimensional ℝ (TangentSpace I x)]
+
+noncomputable instance (x : M) : NormedAddCommGroup (TangentSpace I x) :=
+  show NormedAddCommGroup E from inferInstance
+
+noncomputable instance (x : M) : NormedSpace ℝ (TangentSpace I x) :=
+  show NormedSpace ℝ E from inferInstance
+namespace Manifold.ChartSmoothness
+
+open BilinearForm ContDiff ContinuousLinearMap
+
+variable {X : Type*} [NormedAddCommGroup X] [NormedSpace ℝ X]
+variable {f_bilin : X → E →L[ℝ] E →L[ℝ] ℝ} {s_set : Set X} {x₀_pt : X}
+-- n is from the outer scope (smoothness of the metric)
+
+lemma contDiffWithinAt_eval_bilinear_apply (hf : ContDiffWithinAt ℝ n f_bilin s_set x₀_pt)
+    (v w : E) :
+    ContDiffWithinAt ℝ n (fun x => f_bilin x v w) s_set x₀_pt := by
+  let eval_vw : (E →L[ℝ] E →L[ℝ] ℝ) →L[ℝ] ℝ := BilinearForm.eval_vectors_continuousLinear v w
+  exact (eval_vw.contDiff.of_le le_top).contDiffWithinAt.comp x₀_pt hf (Set.mapsTo_univ _ _)
+
+variable[FiniteDimensional ℝ E]
+
+lemma contDiffWithinAt_bilinear_apply_iff_forall_coord :
+    (∀ v w, ContDiffWithinAt ℝ n (fun x => f_bilin x v w) s_set x₀_pt) →
+    ContDiffWithinAt ℝ n f_bilin s_set x₀_pt := by
+  intro h_comp
+  rw [ContDiff.iff_forall_coord (V := E →L[ℝ] E →L[ℝ] ℝ) (𝕜 := ℝ)]
+  intro φ
+  let b := Module.finBasis ℝ E
+  obtain ⟨e_forms, h_e_forms_def⟩ := BilinearForm.elementary_bilinear_forms_def b
+  have h_f_decomp : ∀ (x : X), f_bilin x =
+      ∑ i ∈ Finset.univ, ∑ j ∈ Finset.univ, ((f_bilin x) (b i) (b j)) • e_forms i j := by
+    intro x
+    obtain ⟨e, h_e_prop, h_decomp⟩ := BilinearForm.decomposition b (f_bilin x)
+    have e_eq_e_forms : e = e_forms := by
+      ext i j v w
+      rw [h_e_prop i j v w, h_e_forms_def i j v w]
+    rw [e_eq_e_forms] at h_decomp
+    exact h_decomp
+  have h_phi_f : ∀ x, φ (f_bilin x) =
+      ∑ i ∈ Finset.univ, ∑ j ∈ Finset.univ, ((f_bilin x) (b i) (b j)) * φ (e_forms i j) := by
+    intro x
+    rw [h_f_decomp x]
+    have h_expand : φ (∑ i ∈ Finset.univ,
+      ∑ j ∈ Finset.univ, ((f_bilin x) (b i) (b j)) • e_forms i j) =
+        ∑ i ∈ Finset.univ, φ (∑ j ∈ Finset.univ, ((f_bilin x) (b i) (b j)) • e_forms i j) :=
+        ContinuousLinearMap.map_finset_sum φ Finset.univ (fun i => ∑ j ∈ Finset.univ,
+        ((f_bilin x) (b i) (b j)) • e_forms i j)
+    rw [h_expand]
+    apply Finset.sum_congr rfl
+    intro i _
+    have h_expand_inner : φ (∑ j ∈ Finset.univ, ((f_bilin x) (b i) (b j)) • e_forms i j) =
+                          ∑ j ∈ Finset.univ, φ (((f_bilin x) (b i) (b j)) • e_forms i j) :=
+      ContinuousLinearMap.map_finset_sum φ Finset.univ (fun j =>
+        ((f_bilin x) (b i) (b j)) • e_forms i j)
+    rw [h_expand_inner]
+    apply Finset.sum_congr rfl
+    intro j _
+    rw [ContinuousLinearMap.map_smul]
+    rw [smul_eq_mul]
+    rw [← h_f_decomp]
+  have h_goal : ContDiffWithinAt ℝ n (fun x => φ (f_bilin x)) s_set x₀_pt := by
+    simp only [h_phi_f]
+    apply ContDiffWithinAt.sum
+    intro i _
+    apply ContDiffWithinAt.sum
+    intro j _
+    have h_term_smooth : ContDiffWithinAt ℝ n (fun x => f_bilin x (b i) (b j)) s_set x₀_pt :=
+      h_comp (b i) (b j)
+    have h_const : ContDiffWithinAt ℝ n (fun x => φ (e_forms i j)) s_set x₀_pt :=
+      contDiffWithinAt_const
+    exact ContDiffWithinAt.mul h_term_smooth h_const
+  exact h_goal
+
+/--
+A bilinear form `f_bilin : X → E → E → F` is continuously differentiable of order `n_level`
+at a point `x₀_pt` within a set `s_set` if and only if for all vectors `v w : E`, the function
+`x ↦ f_bilin x v w` is continuously differentiable of order `n_level` at `x₀_pt` within `s_set`.
+
+This provides an equivalence between the continuous differentiability of a bilinear map
+and the continuous differentiability of all its partial evaluations.
+-/
+theorem contDiffWithinAt_bilinear_iff {n_level : WithTop ℕ∞} :
+    (∀ (v w : E), ContDiffWithinAt ℝ n_level (fun x => f_bilin x v w) s_set x₀_pt) ↔
+    (ContDiffWithinAt ℝ n_level f_bilin s_set x₀_pt) := by
+  constructor
+  · intro h; exact contDiffWithinAt_bilinear_apply_iff_forall_coord h
+  · intro h; exact contDiffWithinAt_eval_bilinear_apply h
+
+/--
+A bilinear form `f_bilin : X → E → E → F` is continuously differentiable of order `n_level`
+on a set `s` if and only if for all vectors `v w : E`, the function `x ↦ f_bilin x v w`
+is continuously differentiable of order `n_level` on `s`.
+
+This extends `contDiffWithinAt_bilinear_iff` from a single point to an entire set.
+-/
+theorem contDiffOn_bilinear_iff {n_level : WithTop ℕ∞} (s : Set X) :
+    (∀ (v w : E), ContDiffOn ℝ n_level (fun x => f_bilin x v w) s) ↔
+    (ContDiffOn ℝ n_level f_bilin s) := by
+  simp only [ContDiffOn, contDiffWithinAt_bilinear_iff (n_level := n_level)]
+  constructor
+  · intro h x hx
+    apply contDiffWithinAt_bilinear_apply_iff_forall_coord
+    intro v w
+    exact h v w x hx
+  · intro h v w x hx
+    exact contDiffWithinAt_eval_bilinear_apply (h x hx) v w
+
+end Manifold.ChartSmoothness
diff --git a/PhysLean/Mathematics/Geometry/Manifold/Chart/CoordinateTransformations.lean b/PhysLean/Mathematics/Geometry/Manifold/Chart/CoordinateTransformations.lean
new file mode 100644
index 000000000..1aadd2ae5
--- /dev/null
+++ b/PhysLean/Mathematics/Geometry/Manifold/Chart/CoordinateTransformations.lean
@@ -0,0 +1,160 @@
+/-
+Copyright (c) 2025 Matteo Cipollina. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Matteo Cipollina
+-/
+
+import PhysLean.Mathematics.Geometry.Manifold.Chart.Smoothness
+import PhysLean.Mathematics.Geometry.Manifold.Chart.Utilities
+import PhysLean.Mathematics.Geometry.Metric.PseudoRiemannian.Defs
+
+/-!
+# Coordinate Transformations and Transition Maps
+
+This file contains lemmas about coordinate transformations and transition maps between charts.
+-/
+
+namespace Manifold.Chart
+open ContDiff ContinuousLinearMap
+open PartialHomeomorph
+open Filter Manifold
+
+variable {E : Type v} {M : Type v} {n : WithTop ℕ∞}
+variable [NormedAddCommGroup E]
+variable [TopologicalSpace M]
+
+/-- The domain where the chart transition map `φ = e' ∘ e.symm` is well-defined forms
+a neighborhood of `y`. -/
+lemma chartMetric_transition_domain_in_nhds (e e' : PartialHomeomorph M E) (y : E)
+    (hy_target : y ∈ e.target) (hz_source' : e.symm y ∈ e'.source) :
+    let _ := e' ∘ e.symm
+    y ∈ e.target ∩ e.symm ⁻¹' e'.source := by
+  intro φ
+  simp only [Set.mem_inter_iff, Set.mem_preimage]
+  exact ⟨hy_target, hz_source'⟩
+
+/-- When applying the transition map `φ = e' ∘ e.symm` to a point `z` in the domain,
+the result is in `e'.target`. -/
+lemma chartMetric_transition_map_target (e e' : PartialHomeomorph M E) {z : E}
+    (hz_source' : e.symm z ∈ e'.source) :
+    let φ := e' ∘ e.symm
+    φ z ∈ e'.target := by
+  intro φ
+  exact e'.mapsTo hz_source'
+
+variable [NormedSpace ℝ E]
+variable [ChartedSpace E M]
+variable {I : ModelWithCorners ℝ E E}
+variable [inst_mani_smooth : IsManifold I (n + 1) M] -- For C^{n+1} manifold for C^n metric
+variable [inst_tangent_findim : ∀ (x : M), FiniteDimensional ℝ (TangentSpace I x)]
+variable {X : Type*} [NormedAddCommGroup X] [NormedSpace ℝ X]
+variable {f_bilin : X → E →L[ℝ] E →L[ℝ] ℝ} {s_set : Set X} {x₀_pt : X}
+
+/-- Core derivative chain rule for manifold derivatives that combines correctly with the metric. -/
+theorem manifold_derivative_chain_rule (g : PseudoRiemannianMetric E E M n I)
+    (e e' : PartialHomeomorph M E) {z : E}
+    (hz_target : z ∈ e.target) (hz_source' : e.symm z ∈ e'.source)
+    (he : e ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (he' : e' ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (hn1 : 1 ≤ n) :
+    let φ := e' ∘ e.symm
+    let ψ := e.symm.trans e'
+    let x_z := e.symm z
+    ∀ v w : E,
+      g.val x_z (mfderiv I I e.symm z v) (mfderiv I I e.symm z w) =
+        g.val x_z (mfderiv I I e'.symm (φ z) ((mfderiv I I ψ z) v))
+                 (mfderiv I I e'.symm (φ z) ((mfderiv I I ψ z) w)) := by
+  intro φ ψ x_z v w
+  have hmani : IsManifold I n M := by
+    letI := inst_mani_smooth
+    have h_le : n ≤ n + 1 := by simp [self_le_add_right]
+    exact IsManifold.of_le h_le
+  have h_chain_rule_v : mfderiv I I e'.symm (φ z) ((mfderiv I I ψ z) v) =
+                        mfderiv I I e.symm z v := by
+    let x := e.symm z
+    have h_z_eq : z = e x := Eq.symm (PartialHomeomorph.right_inv e hz_target)
+    have hx_e : x ∈ e.source := e.map_target hz_target
+    have h_fun_eq : e.symm =ᶠ[nhds z] e'.symm ∘ φ := by
+      have h_z_eq_ex : z = e x := by
+        exact h_z_eq
+      have h_lemma : e.symm =ᶠ[nhds (e x)] e'.symm ∘ ↑(e.symm.trans e') :=
+        chart_inverse_composition e e' x hx_e hz_source'
+      have h_φ_eq : φ = ↑(e.symm.trans e') := by
+        ext y
+        simp [φ, trans_apply, Function.comp_apply]
+      rw [← h_z_eq_ex, ← h_φ_eq] at h_lemma
+      exact h_lemma
+    have h_deriv_eq : mfderiv I I e.symm z = mfderiv I I (e'.symm ∘ φ) z :=
+      Filter.EventuallyEq.mfderiv_eq h_fun_eq
+    have h_chain : mfderiv I I (e'.symm ∘ φ) z v =
+                   mfderiv I I e'.symm (φ z) (mfderiv I I φ z v) := by
+      have h_comp : mfderiv I I (e'.symm ∘ φ) z =
+                    (mfderiv I I e'.symm (φ z)).comp (mfderiv I I φ z) := by
+        have hφz_target : φ z ∈ e'.target :=
+          chartMetric_phi_maps_to_target e e' hz_source'
+        have hφ_diff : MDifferentiableAt I I φ z := by
+          have h_φ_eq_ψ : φ = ψ := rfl
+          rw [h_φ_eq_ψ]
+          let x := e.symm z
+          have hx_e := e.map_target hz_target
+          have h_smooth := chart_transition_smooth e e' x hx_e hz_source' n he he'
+          simpa [h_z_eq] using h_smooth.mdifferentiableAt hn1
+        have he'_symm_diff : MDifferentiableAt I I e'.symm (φ z) := by
+          have hmani : IsManifold I n M := by
+            letI := inst_mani_smooth
+            have h : n ≤ n + 1 := by
+              simp only [self_le_add_right, φ]
+            exact IsManifold.of_le h
+          have hn1 : 1 ≤ n := by exact hn1
+          have h_symm_smooth : ContMDiffOn I I n e'.symm e'.target :=
+            contMDiffOn_symm_of_mem_maximalAtlas he'
+          exact
+            (h_symm_smooth.contMDiffAt (e'.open_target.mem_nhds hφz_target)).mdifferentiableAt hn1
+        exact mfderiv_comp_of_eq he'_symm_diff hφ_diff rfl
+      have h_apply : mfderiv I I (e'.symm ∘ φ) z v =
+                     (mfderiv I I e'.symm (φ z)).comp (mfderiv I I φ z) v := by
+        rw [h_comp]
+        exact rfl
+      rw [h_apply]
+      rfl
+    have h_φ_eq_ψ : φ = ψ := rfl
+    rw [← h_deriv_eq, h_φ_eq_ψ] at h_chain
+    exact h_chain.symm
+  have h_chain_rule_w : mfderiv I I e'.symm (φ z) ((mfderiv I I ψ z) w) =
+                      mfderiv I I e.symm z w := by
+    let x := e.symm z
+    have h_z_eq : z = e x := Eq.symm (PartialHomeomorph.right_inv e hz_target)
+    have hx_e : x ∈ e.source := e.map_target hz_target
+    have h_fun_eq : e.symm =ᶠ[nhds z] e'.symm ∘ φ := by
+      have h_z_eq_ex : z = e x := h_z_eq
+      have h_lemma := chart_inverse_composition e e' x hx_e hz_source'
+      have h_φ_eq : φ = (e.symm.trans e') := rfl
+      simpa [h_z_eq_ex, h_φ_eq] using h_lemma
+    have h_deriv_eq : mfderiv I I e.symm z = mfderiv I I (e'.symm ∘ φ) z :=
+      Filter.EventuallyEq.mfderiv_eq h_fun_eq
+    let h_comp := mfderiv_chart_transition_helper e e' x hx_e hz_source' hn1 hmani he he'
+    have h_comp_w : mfderiv I I (e'.symm ∘ φ) z w =
+                   (mfderiv I I e'.symm (φ z)).comp (mfderiv I I φ z) w := by
+      have h_comp_map : mfderiv I I (e'.symm ∘ φ) z =
+                        (mfderiv I I e'.symm (φ z)).comp (mfderiv I I φ z) := by
+        have hφz_target : φ z ∈ e'.target :=
+          chartMetric_phi_maps_to_target e e' hz_source'
+        have hφ_diff : MDifferentiableAt I I φ z := by
+          have h_φ_eq_ψ : φ = ψ := rfl
+          rw [h_φ_eq_ψ]
+          have h_smooth := chart_transition_smooth e e' x hx_e hz_source' n he he'
+          simpa [h_z_eq] using h_smooth.mdifferentiableAt hn1
+        have he'_symm_diff : MDifferentiableAt I I e'.symm (φ z) := by
+          have h_symm_smooth : ContMDiffOn I I n e'.symm e'.target :=
+            contMDiffOn_symm_of_mem_maximalAtlas he'
+          exact
+            (h_symm_smooth.contMDiffAt (e'.open_target.mem_nhds hφz_target)).mdifferentiableAt hn1
+        exact mfderiv_comp_of_eq he'_symm_diff hφ_diff rfl
+      rw [h_comp_map]
+      rfl
+    have h_φ_eq_ψ : φ = ψ := rfl
+    rw [← h_deriv_eq, h_φ_eq_ψ] at h_comp_w
+    exact h_comp_w.symm
+  rw [h_chain_rule_v, h_chain_rule_w]
+
+end Manifold.Chart
diff --git a/PhysLean/Mathematics/Geometry/Manifold/Chart/Smoothness.lean b/PhysLean/Mathematics/Geometry/Manifold/Chart/Smoothness.lean
new file mode 100644
index 000000000..e0d5f960c
--- /dev/null
+++ b/PhysLean/Mathematics/Geometry/Manifold/Chart/Smoothness.lean
@@ -0,0 +1,135 @@
+/-
+Copyright (c) 2025 Matteo Cipollina. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Matteo Cipollina
+-/
+
+import Mathlib.Geometry.Manifold.MFDeriv.Basic
+import PhysLean.Mathematics.Geometry.Manifold.PartialHomeomorph
+
+/-!
+# Smoothness of Charts and Transition Maps
+
+This file contains lemmas about the smoothness and differentiability of charts and transition maps.
+-/
+
+namespace Manifold.Chart
+
+variable {E : Type v} {M : Type v} {n : WithTop ℕ∞}
+variable [NormedAddCommGroup E]
+variable [TopologicalSpace M]
+
+open PartialHomeomorph ContinuousLinearMap PartialEquiv
+
+/-- Coordinate transformation identity: e' x = φ(e x) where φ is the transition map. -/
+lemma chart_coordinate_transform (e e' : PartialHomeomorph M E) (x : M)
+    (hx_e : x ∈ e.source) :
+    let φ := e.symm.trans e';
+    let y := e x;
+    e' x = φ y := (apply_transition_map_eq_chart_apply e e' x hx_e).symm
+
+/-- Chart inverse composition identity: e.symm = e'.symm ∘ φ in a neighborhood of e x,
+    where φ is the transition map. -/
+lemma chart_inverse_composition (e e' : PartialHomeomorph M E) (x : M)
+    (hx_e : x ∈ e.source) (hx_e' : x ∈ e'.source) :
+    let φ := e.symm.trans e';
+    let y := e x;
+    e.symm =ᶠ[nhds y] e'.symm ∘ φ := by
+  intro φ y
+  have hy_dom_φ : y ∈ φ.source :=
+    transition_map_source_contains_image_point e e' x hx_e hx_e'
+  filter_upwards [φ.open_source.mem_nhds hy_dom_φ] with z hz
+  exact apply_symm_eq_apply_symm_comp_transition e e' z hz
+
+variable [NormedSpace ℝ E]
+variable [ChartedSpace E M]
+variable {I : ModelWithCorners ℝ E E}
+
+open PartialHomeomorph ContinuousLinearMap PartialEquiv
+
+/-- Transition maps between charts in the maximal atlas are smooth. -/
+lemma chart_transition_smooth (e e' : PartialHomeomorph M E) (x : M)
+    (hx_e : x ∈ e.source) (hx_e' : x ∈ e'.source)
+    (n : WithTop ℕ∞)
+    (he : e ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (he' : e' ∈ (contDiffGroupoid n I).maximalAtlas M) :
+    let φ := e.symm.trans e';
+    let y := e x;
+    ContMDiffAt I I n φ y := by
+  intro φ y
+  have hy_dom_φ : y ∈ φ.source :=
+    transition_map_source_contains_image_point e e' x hx_e hx_e'
+  have h_φ_groupoid : φ ∈ contDiffGroupoid n I :=
+    (contDiffGroupoid n I).compatible_of_mem_maximalAtlas he he'
+  have h_φ_smooth : ContMDiffOn I I n φ φ.source :=
+    contMDiffOn_of_mem_contDiffGroupoid h_φ_groupoid
+  exact ContMDiffOn.contMDiffAt h_φ_smooth (φ.open_source.mem_nhds hy_dom_φ)
+
+variable {E : Type v} {M : Type v} {n : WithTop ℕ∞}
+variable [NormedAddCommGroup E] [NormedSpace ℝ E] --[FiniteDimensional ℝ E]
+variable [TopologicalSpace M] [ChartedSpace E M] --[T2Space M]
+variable {I : ModelWithCorners ℝ E E} --[ModelWithCorners.Boundaryless I]
+
+/-- Charts in the maximal atlas are differentiable on their source. -/
+lemma chart_differentiable_on_source
+    (e : PartialHomeomorph M E) (hmani : IsManifold I n M)
+    (he : e ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (hn1 : 1 ≤ n) :
+    MDifferentiableOn I I e e.source := by
+  intro x hx
+  have h_smooth : ContMDiffAt I I n e x :=
+    (contMDiffOn_of_mem_maximalAtlas he).contMDiffAt (e.open_source.mem_nhds hx)
+  exact (h_smooth.mdifferentiableWithinAt hn1).mono (by exact fun ⦃a⦄ a => trivial)
+
+lemma chart_symm_differentiable_on_target
+    (e : PartialHomeomorph M E) (hmani : IsManifold I n M)
+    (he : e ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (hn1 : 1 ≤ n) :
+    MDifferentiableOn I I e.symm e.target := by
+  intro y hy
+  have h_smooth : ContMDiffAt I I n e.symm y :=
+    (contMDiffOn_symm_of_mem_maximalAtlas he).contMDiffAt (e.open_target.mem_nhds hy)
+  let h_md := h_smooth.mdifferentiableWithinAt hn1
+  exact h_md.mono (Set.subset_univ _)
+
+/-- The manifold derivative of a chart transition map can be computed from charts. -/
+lemma mfderiv_chart_transition_helper (e e' : PartialHomeomorph M E) (x : M)
+    (hx_e : x ∈ e.source) (hx_e' : x ∈ e'.source)
+    (hn1 : 1 ≤ n) (hmani : IsManifold I n M)
+    (he : e ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (he' : e' ∈ (contDiffGroupoid n I).maximalAtlas M) :
+    let φ := e.symm.trans e';
+    let y := e x;
+    mfderiv I I (e'.symm ∘ φ) y = (mfderiv I I e'.symm (φ y)).comp (mfderiv I I φ y) := by
+  intro φ y
+  have h_φ_y_eq : φ y = e' x := (chart_coordinate_transform e e' x hx_e).symm
+  have hφ_diff : ContMDiffAt I I n φ y :=
+    chart_transition_smooth e e' x hx_e hx_e' n he he'
+  have hφ_mdiff : MDifferentiableAt I I φ y := hφ_diff.mdifferentiableAt hn1
+  have he'_symm_diff : MDifferentiableAt I I e'.symm (e' x) :=
+    ((chart_symm_differentiable_on_target e' hmani he' hn1) (e' x)
+      (e'.mapsTo hx_e')).mdifferentiableAt
+      (e'.open_target.mem_nhds (e'.mapsTo hx_e'))
+  exact mfderiv_comp_of_eq he'_symm_diff hφ_mdiff h_φ_y_eq
+
+/-- Chain rule for composition of chart inverses with transition maps. -/
+theorem mfderiv_chart_transition (e e' : PartialHomeomorph M E) (x : M)
+    (hx_e : x ∈ e.source) (hx_e' : x ∈ e'.source)
+    (hn1 : 1 ≤ n) (hmani : IsManifold I n M)
+    (he : e ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (he' : e' ∈ (contDiffGroupoid n I).maximalAtlas M) :
+    let φ := e.symm.trans e';
+    let y := e x;
+    mfderiv I I e.symm y = (mfderiv I I e'.symm (e' x)).comp (mfderiv I I φ y) := by
+  intro φ y
+  have h_fun_eq : e.symm =ᶠ[nhds y] e'.symm ∘ φ :=
+    chart_inverse_composition e e' x hx_e hx_e'
+  have h_φ_y_eq : φ y = e' x :=
+    (chart_coordinate_transform e e' x hx_e).symm
+  have h_deriv_eq : mfderiv I I e.symm y = mfderiv I I (e'.symm ∘ φ) y :=
+    Filter.EventuallyEq.mfderiv_eq h_fun_eq
+  rw [h_deriv_eq]
+  have h_comp := mfderiv_chart_transition_helper e e' x hx_e hx_e' hn1 hmani he he'
+  rw [h_comp, h_φ_y_eq]
+
+end Manifold.Chart
diff --git a/PhysLean/Mathematics/Geometry/Manifold/Chart/Utilities.lean b/PhysLean/Mathematics/Geometry/Manifold/Chart/Utilities.lean
new file mode 100644
index 000000000..fc1600faa
--- /dev/null
+++ b/PhysLean/Mathematics/Geometry/Manifold/Chart/Utilities.lean
@@ -0,0 +1,59 @@
+/-
+Copyright (c) 2025 Matteo Cipollina. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Matteo Cipollina
+-/
+
+import Mathlib.Analysis.Normed.Group.Basic
+import Mathlib.Topology.PartialHomeomorph
+/-!
+# Utilities for Charts and Transition Maps
+
+This file contains general utility lemmas for charts and transition maps.
+-/
+
+namespace Manifold.Chart
+variable {E : Type v} {M : Type v} {n : WithTop ℕ∞}
+variable [NormedAddCommGroup E]
+variable [TopologicalSpace M]
+
+/-- The domain where the chart transition is well-defined forms a neighborhood of y. -/
+lemma chartMetric_domain_in_nhds (e e' : PartialHomeomorph M E) (y : E)
+    (hy : y ∈ e.target) (hx_source' : e.symm y ∈ e'.source) :
+    let s := e.target ∩ e.symm ⁻¹' e'.source
+    s ∈ nhds y := by
+  intro s
+  apply Filter.inter_mem
+  · exact e.open_target.mem_nhds hy
+  · exact (e.continuousAt_symm hy).preimage_mem_nhds (e'.open_source.mem_nhds hx_source')
+
+/-- When applying φ = e' ∘ e.symm to a point z in the domain s, the result is in e'.target -/
+lemma chartMetric_phi_in_target (e e' : PartialHomeomorph M E) {z : E}
+    (hz_source' : e.symm z ∈ e'.source) :
+    let φ := e'.toPartialEquiv ∘ e.symm
+    φ z ∈ e'.target := by
+  intro φ
+  unfold φ
+  exact e'.mapsTo hz_source'
+
+/-- When applying φ = e' ∘ e.symm to a point z, the result is in e'.target -/
+lemma chartMetric_phi_maps_to_target (e e' : PartialHomeomorph M E) {z : E}
+    (hz_source' : e.symm z ∈ e'.source) :
+    let φ := e' ∘ e.symm
+    φ z ∈ e'.target := by
+  intro φ
+  unfold φ
+  exact e'.mapsTo hz_source'
+
+/-- The composition of chart maps preserves points: e'.symm (e' (e.symm z)) = e.symm z -/
+lemma chart_map_composition_identity (e e' : PartialHomeomorph M E) {z : E}
+    (hz_source' : e.symm z ∈ e'.source) :
+    let φ := e' ∘ e.symm
+    let x_z := e.symm z
+    e'.symm (φ z) = x_z := by
+  intro φ x_z
+  unfold φ
+  simp only [Function.comp_apply, PartialHomeomorph.coe_coe]
+  rw [e'.left_inv hz_source']
+
+end Manifold.Chart
diff --git a/PhysLean/Mathematics/Geometry/Manifold/PartialHomeomorph.lean b/PhysLean/Mathematics/Geometry/Manifold/PartialHomeomorph.lean
new file mode 100644
index 000000000..413363084
--- /dev/null
+++ b/PhysLean/Mathematics/Geometry/Manifold/PartialHomeomorph.lean
@@ -0,0 +1,209 @@
+/-
+Copyright (c) 2025 Matteo Cipollina. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Matteo Cipollina
+-/
+
+import Mathlib.Geometry.Manifold.ContMDiff.Atlas
+import Mathlib.Geometry.Manifold.ContMDiff.NormedSpace
+
+/-!
+# PartialHomeomorph Utilities for Manifolds
+
+This file contains lemmas for `PartialHomeomorph` specifically relevant to manifolds.
+-/
+
+namespace PartialHomeomorph
+
+open Set ModelWithCorners PartialEquiv
+
+variable {𝕜 : Type*} [NontriviallyNormedField 𝕜]
+variable {E : Type*} [NormedAddCommGroup E] [NormedSpace 𝕜 E]
+variable {H : Type*} [NormedAddCommGroup H] [NormedSpace 𝕜 H] [TopologicalSpace H]
+variable {M : Type*} [TopologicalSpace M] [ChartedSpace H M]
+variable {I : ModelWithCorners 𝕜 E H} {n : WithTop ℕ∞}
+
+/-- The toPartialEquiv of the symm of a partial homeomorphism equals the symm of its
+toPartialEquiv. -/
+@[simp, mfld_simps]
+theorem toPartialEquiv_symm {X Y : Type*} [TopologicalSpace X] [TopologicalSpace Y]
+    (e : PartialHomeomorph X Y) :
+    e.symm.toPartialEquiv = e.toPartialEquiv.symm :=
+  rfl
+
+/-- The function component of the symm of a partial homeomorphism equals its invFun. -/
+theorem coe_symm_eq_invFun {X Y : Type*} [TopologicalSpace X] [TopologicalSpace Y]
+    (e : PartialHomeomorph X Y) : ⇑e.symm = e.invFun :=
+  rfl
+
+omit [NormedAddCommGroup H] [NormedSpace 𝕜 H] in
+/-- If a partial homeomorphism belongs to the `contDiffGroupoid n I`, then the forward map is
+`C^n` on its source. This extracts the smoothness property directly from groupoid membership. -/
+theorem contMDiffOn_of_mem_groupoid {e : PartialHomeomorph H H}
+    (he : e ∈ contDiffGroupoid n I) : ContMDiffOn I I n e e.source :=
+  contMDiffOn_of_mem_contDiffGroupoid he
+
+omit [NormedAddCommGroup H] [NormedSpace 𝕜 H] in
+/-- The coordinate change map `φ = e' ∘ e.symm` between two charts `e` and `e'` is `C^n` smooth
+    at a point `y = e x` if `x` is in the intersection of the sources of `e` and `e'`. -/
+lemma contMDiffAt_coord_change
+    {e e' : PartialHomeomorph M H} {x : M}
+    (hx_e_source : x ∈ e.source) (hx_e'_source : x ∈ e'.source)
+    (he : e ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (he' : e' ∈ (contDiffGroupoid n I).maximalAtlas M) :
+    ContMDiffAt I I n (e' ∘ e.symm) (e x) := by
+  let y := e x
+  have hy_target : y ∈ e.target := e.mapsTo hx_e_source
+  let transition_domain := e.target ∩ e.symm ⁻¹' e'.source
+  have hy_domain : y ∈ transition_domain := ⟨hy_target, by
+    rw [@Set.mem_preimage];
+    rw [e.left_inv hx_e_source];
+    exact hx_e'_source⟩
+  have open_domain : IsOpen transition_domain := e.isOpen_inter_preimage_symm e'.open_source
+  have h_compatible := StructureGroupoid.compatible_of_mem_maximalAtlas he he'
+  have h_transition_smooth : ContMDiffOn I I n (e' ∘ e.symm) transition_domain :=
+    contMDiffOn_of_mem_groupoid h_compatible
+  exact h_transition_smooth.contMDiffAt (open_domain.mem_nhds hy_domain)
+
+omit [NormedAddCommGroup H] [NormedSpace 𝕜 H] in
+/-- The coordinate change map between atlas charts is smooth in the manifold sense. -/
+lemma contMDiffAt_atlas_coord_change [IsManifold I n M]
+    {e e' : PartialHomeomorph M H}
+    (he : e ∈ atlas H M) (he' : e' ∈ atlas H M) {x : M}
+    (hx_e_source : x ∈ e.source) (hx_e'_source : x ∈ e'.source) :
+    ContMDiffAt I I n (e' ∘ e.symm) (e x) := by
+  let transition_domain := e.target ∩ e.symm ⁻¹' e'.source
+  have open_domain : IsOpen transition_domain := e.isOpen_inter_preimage_symm e'.open_source
+  have hy_target : e x ∈ e.target := e.mapsTo hx_e_source
+  have hy_domain : e x ∈ transition_domain := ⟨hy_target, by
+    show e.symm (e x) ∈ e'.source
+    rw [e.left_inv hx_e_source]
+    exact hx_e'_source⟩
+  have h_compatible : e.symm.trans e' ∈ contDiffGroupoid n I := HasGroupoid.compatible he he'
+  have h_transition_smooth : ContMDiffOn I I n (e' ∘ e.symm) transition_domain :=
+    contMDiffOn_of_mem_groupoid h_compatible
+  exact h_transition_smooth.contMDiffAt (open_domain.mem_nhds hy_domain)
+
+/-- The coordinate change map `φ = e' ∘ e.symm` between two charts `e` and `e'` taken from the
+    atlas is `C^n` smooth at a point `y = e x` if `x` is in the intersection of their sources. -/
+lemma contDiffAt_atlas_coord_change
+    {𝕜 : Type*} [NontriviallyNormedField 𝕜] [NormedSpace 𝕜 E]
+    {n : WithTop ℕ∞} [ChartedSpace E M]
+    [IsManifold (modelWithCornersSelf 𝕜 E) n M]
+    {e e' : PartialHomeomorph M E}
+    (he : e ∈ atlas E M) (he' : e' ∈ atlas E M) {x : M}
+    (hx_e_source : x ∈ e.source) (hx_e'_source : x ∈ e'.source) :
+    let y := e x
+    ContDiffAt 𝕜 n (e' ∘ e.symm) y := by
+  intro y
+  have hy_target : y ∈ e.target := e.mapsTo hx_e_source
+  let transition_domain := e.target ∩ e.symm ⁻¹' e'.source
+  have hy_domain : y ∈ transition_domain := by
+    refine ⟨hy_target, ?_⟩
+    simp only [Set.mem_preimage]
+    rw [e.left_inv hx_e_source]
+    exact hx_e'_source
+  have h_compatible : e.symm.trans e' ∈ contDiffGroupoid n (modelWithCornersSelf 𝕜 E) :=
+    HasGroupoid.compatible he he'
+  have h_transition_mfd_smooth : ContMDiffOn (modelWithCornersSelf 𝕜 E)
+      (modelWithCornersSelf 𝕜 E) n (e' ∘ e.symm) transition_domain :=
+    contMDiffOn_of_mem_contDiffGroupoid h_compatible
+  have h_transition_smooth : ContDiffOn 𝕜 n (e' ∘ e.symm) transition_domain :=
+    (contMDiffOn_iff_contDiffOn).mp h_transition_mfd_smooth
+  have open_domain : IsOpen transition_domain :=
+    ContinuousOn.isOpen_inter_preimage e.continuousOn_invFun e.open_target e'.open_source
+  exact h_transition_smooth.contDiffAt (IsOpen.mem_nhds open_domain hy_domain)
+
+lemma extChartAt_eq_extend_chartAt {E' : Type*} [NormedAddCommGroup E'] [NormedSpace 𝕜 E']
+    (I' : ModelWithCorners 𝕜 E' E') (x_m : M) [ChartedSpace E' M] :
+    extChartAt I' x_m = (chartAt E' x_m).extend I' :=
+  rfl
+
+variable {F : Type*} [NormedAddCommGroup F] [NormedSpace 𝕜 F]
+variable [K_bl : I.Boundaryless] -- Changed K to I
+
+omit [NormedSpace 𝕜 E] in
+/-- Domain of the transition map `e.symm.trans e'` contains `e x`. -/
+lemma transition_map_source_contains_image_point {E_chart : Type*} [TopologicalSpace E_chart]
+    (e e' : PartialHomeomorph M E_chart) (x : M)
+    (hx_e_source : x ∈ e.source) (hx_e'_source : x ∈ e'.source) :
+    (e x) ∈ (e.symm.toPartialEquiv.trans e'.toPartialEquiv).source := by
+  rw [@PartialEquiv.trans_source'']
+  simp only [symm_toPartialEquiv, PartialEquiv.symm_symm, toFun_eq_coe, PartialEquiv.symm_target,
+    Set.mem_image, Set.mem_inter_iff]
+  rw [← @symm_image_target_eq_source]
+  rw [@symm_image_target_eq_source]
+  apply Exists.intro
+  · apply And.intro
+    on_goal 2 => {rfl
+    }
+    · simp_all only [and_self]
+
+variable {f_map : H → H} {y_pt : H} -- Renamed f, y to avoid clash
+
+omit [NormedAddCommGroup H] [NormedSpace 𝕜 H] K_bl in
+lemma contDiffAt_chart_expression_of_contMDiffAt [I.Boundaryless]
+    (h_contMDiff : ContMDiffAt I I n f_map y_pt) :
+    ContDiffAt 𝕜 n (I ∘ f_map ∘ I.symm) (I y_pt) := by
+  rw [@contMDiffAt_iff] at h_contMDiff
+  obtain ⟨_, h_diff⟩ := h_contMDiff
+  have h_fun_eq : (extChartAt I (f_map y_pt)) ∘ f_map ∘ (extChartAt I y_pt).symm =
+      I ∘ f_map ∘ I.symm := by
+    simp only [extChartAt, extend, refl_partialEquiv, PartialEquiv.refl_source,
+      singletonChartedSpace_chartAt_eq, PartialEquiv.refl_trans, toPartialEquiv_coe,
+      toPartialEquiv_coe_symm]
+  have h_point_eq : (extChartAt I y_pt) y_pt = I y_pt := by
+    simp only [extChartAt, extend, refl_partialEquiv, PartialEquiv.refl_source,
+      singletonChartedSpace_chartAt_eq, PartialEquiv.refl_trans, toPartialEquiv_coe]
+  have h_diff' : ContDiffWithinAt 𝕜 n (I ∘ f_map ∘ I.symm) (Set.range I) (I y_pt) := by
+    rw [← h_fun_eq, ← h_point_eq]
+    exact h_diff
+  apply ContDiffWithinAt.contDiffAt h_diff'
+  have : I y_pt ∈ interior (Set.range I) := by
+    have h_range : Set.range I = Set.univ := I.range_eq_univ
+    simp [h_range, interior_univ]
+  exact mem_interior_iff_mem_nhds.mp this
+
+variable {E : Type*} [NormedAddCommGroup E] [NormedSpace ℝ E]
+variable {I : ModelWithCorners ℝ E E} [I.Boundaryless]
+variable {n : WithTop ℕ∞}
+
+lemma contDiff_of_boundaryless_model_on_self :
+    ContDiff ℝ n (modelWithCornersSelf ℝ E).toFun := by
+  let I := modelWithCornersSelf ℝ E
+  have h_mem_groupoid : PartialHomeomorph.refl E ∈ contDiffGroupoid n I :=
+    StructureGroupoid.id_mem (contDiffGroupoid n I)
+  have h_contMDiffOn : ContMDiffOn I I n (PartialHomeomorph.refl E) Set.univ :=
+    contMDiffOn_of_mem_contDiffGroupoid h_mem_groupoid
+  rw [contMDiffOn_iff_contDiffOn] at h_contMDiffOn
+  simpa [contDiffOn_univ] using h_contMDiffOn
+
+lemma contDiff_symm_of_boundaryless_model_on_self :
+    ContDiff ℝ n (modelWithCornersSelf ℝ E).symm.toFun := by
+  let I := modelWithCornersSelf ℝ E
+  have h_mem_groupoid : PartialHomeomorph.refl E ∈ contDiffGroupoid n I :=
+    StructureGroupoid.id_mem (contDiffGroupoid n I)
+  have h_mdf := contMDiffOn_of_mem_contDiffGroupoid h_mem_groupoid
+  exact contDiffOn_univ.1 (contMDiffOn_iff_contDiffOn.1 h_mdf)
+
+omit [NormedSpace ℝ E]  in
+/-- Functional equality e.symm = e'.symm ∘ φ-/
+lemma apply_symm_eq_apply_symm_comp_transition (e e' : PartialHomeomorph M E) :
+    let φ_pe := e.symm.toPartialEquiv.trans e'.toPartialEquiv;
+    ∀ z ∈ φ_pe.source, e.symm z = (e'.symm ∘ ⇑φ_pe) z := by
+  intro φ_pe z hz_in_φ_source
+  rw [PartialEquiv.trans_source, Set.mem_inter_iff, Set.mem_preimage] at hz_in_φ_source
+  rcases hz_in_φ_source with ⟨_, hz_esymm_z_in_eprime_source⟩
+  simp only [Function.comp_apply, PartialEquiv.trans_apply,
+             PartialHomeomorph.coe_coe]
+  exact Eq.symm (e'.left_inv' hz_esymm_z_in_eprime_source)
+
+omit [NormedSpace ℝ E]  in
+/-- Image of y under φ-/
+lemma apply_transition_map_eq_chart_apply (e e' : PartialHomeomorph M E) (x : M)
+    (hx_e_source : x ∈ e.source) :
+    (e.symm.toPartialEquiv.trans e'.toPartialEquiv) (e x) = e' x := by
+  simp only [PartialEquiv.trans_apply, PartialHomeomorph.coe_coe,
+             e.left_inv hx_e_source]
+
+end PartialHomeomorph
diff --git a/PhysLean/Mathematics/Geometry/Metric/PseudoRiemannian/Chart.lean b/PhysLean/Mathematics/Geometry/Metric/PseudoRiemannian/Chart.lean
new file mode 100644
index 000000000..3bd811af9
--- /dev/null
+++ b/PhysLean/Mathematics/Geometry/Metric/PseudoRiemannian/Chart.lean
@@ -0,0 +1,157 @@
+/-
+Copyright (c) 2025 Matteo Cipollina. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Matteo Cipollina
+-/
+
+import Mathlib.Algebra.Lie.OfAssociative
+import PhysLean.Mathematics.LinearAlgebra.BilinearForm
+import PhysLean.Mathematics.Geometry.Manifold.Chart.Smoothness
+import PhysLean.Mathematics.Geometry.Metric.PseudoRiemannian.Defs
+
+/-!
+# Pseudo-Riemannian Metric in Chart Coordinates
+
+This file defines `chartMetric`, which expresses a pseudo-Riemannian metric `g`
+on a manifold `M` in local chart coordinates `e`. It establishes the
+tensor transformation law for `chartMetric` under a change of chart,
+`chartMetric_coord_change`. Auxiliary lemmas concern the smoothness of
+bilinear forms when viewed as maps into function spaces, and properties
+of chart transitions necessary for proving the main transformation result.
+-/
+
+namespace PseudoRiemannianMetric
+
+noncomputable section
+open BilinearForm PartialHomeomorph ContinuousLinearMap PartialEquiv ContDiff Manifold.Chart
+
+open Filter
+
+variable {E : Type v} {M : Type v} {n : WithTop ℕ∞}
+variable [NormedAddCommGroup E] [NormedSpace ℝ E]
+variable [TopologicalSpace M] [ChartedSpace E M]
+variable {I : ModelWithCorners ℝ E E}
+variable [inst_mani_smooth : IsManifold I (n + 1) M]
+variable [inst_tangent_findim : ∀ (x : M), FiniteDimensional ℝ (TangentSpace I x)]
+
+noncomputable instance (x : M) : NormedAddCommGroup (TangentSpace I x) :=
+  show NormedAddCommGroup E from inferInstance
+
+noncomputable instance (x : M) : NormedSpace ℝ (TangentSpace I x) :=
+  show NormedSpace ℝ E from inferInstance
+
+/-- Given a pseudo-Riemannian metric `g` on a manifold `M` and a chart `e : PartialHomeomorph M E`,
+`chartMetric g e y` computes the metric tensor in the local coordinates provided by `e`.
+Let `x := e.symm y` be the point on `M` corresponding to chart coordinates `y : E`.
+Let `Df_y := mfderiv I I e.symm y : E →L[ℝ] TangentSpace I x` be the differential of `e.symm` at
+`y`. Then `chartMetric g e y` is the bilinear form on `E` defined as the pullback of `g.val x` by
+`Df_y`. That is, for `v, w : E`, `(chartMetric g e y) v w = g.val x (Df_y v) (Df_y w)`.
+If `(b_i)` is a basis for `E`, then the values `(chartMetric g e y) (b_i) (b_j)` are the
+components `g_{ij}(y)` of the metric `g` at `x` with respect to the basis for `T_xM`
+induced by `(b_i)` via `Df_y`.
+If `y ∉ e.target`, the result is defined as the zero bilinear form. -/
+def chartMetric (g : PseudoRiemannianMetric E E M n I) (e : PartialHomeomorph M E) (y : E) :
+    E →L[ℝ] E →L[ℝ] ℝ :=
+  letI := Classical.propDecidable
+  dite (y ∈ e.target)
+    (fun _ : y ∈ e.target =>
+      let x := e.symm y
+      letI : FiniteDimensional ℝ (TangentSpace I x) := inferInstance
+      let Df : E →L[ℝ] TangentSpace I x := mfderiv I I e.symm y
+      BilinearForm.pullback (g.val x) Df)
+    (fun _ => (0 : E →L[ℝ] E →L[ℝ] ℝ))
+
+/-- Evaluation of `chartMetric` at `y ∈ e.target`.
+For `v, w : E`, `chartMetric g e y v w = g_x(D(e⁻¹)_y v, D(e⁻¹)_y w)`, where `x = e⁻¹y`. -/
+@[simp]
+lemma chartMetric_apply_of_mem (g : PseudoRiemannianMetric E E M n I) (e : PartialHomeomorph M E)
+    {y : E} (hy : y ∈ e.target) (v w : E) :
+    chartMetric g e y v w = g.val (e.symm y) (mfderiv I I e.symm y v) (mfderiv I I e.symm y w) := by
+  letI := Classical.propDecidable
+  rw [chartMetric, dite_eq_ite, if_pos hy]
+  simp only [BilinearForm.pullback, ContinuousLinearMap.bilinearComp_apply]
+  exact rfl
+
+/-- If `y ∉ e.target`, `chartMetric g e y` is the zero bilinear form. -/
+lemma chartMetric_apply_of_not_mem
+    (g : PseudoRiemannianMetric E E M n I)
+    (e : PartialHomeomorph M E)
+    {y : E} (hy : y ∉ e.target) (v w : E) :
+    chartMetric g e y v w = 0 := by
+  simp only [chartMetric, hy, ↓reduceDIte, ContinuousLinearMap.zero_apply]
+
+/-- The metric `g_x(D(e⁻¹)_y v, D(e⁻¹)_y w)` can be expressed using a second chart `e'`
+via the chain rule for `D(e⁻¹) = D(e'⁻¹) ∘ D(φ)`.
+`x` is a point on the manifold, `y = ex`, `y' = e'x`.
+`φ = e' ∘ e⁻¹` is the transition map. `1 ≤ n`. -/
+lemma chartMetric_at_corresponding_points (g : PseudoRiemannianMetric E E M n I)
+    (e e' : PartialHomeomorph M E) (x : M)
+    (hx_e : x ∈ e.source) (hx_e' : x ∈ e'.source)
+    (hmani : IsManifold I n M)
+    (hn1 : 1 ≤ n)
+    (he : e ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (he' : e' ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (v w : E) :
+    let y := e x
+    let y' := e' x
+    g.val x (mfderiv I I e.symm y v) (mfderiv I I e.symm y w) =
+    g.val x (mfderiv I I e'.symm y' (mfderiv I I (e.symm.trans e') y v))
+            (mfderiv I I e'.symm y' (mfderiv I I (e.symm.trans e') y w)) := by
+  intro y y'
+  have h_chain := mfderiv_chart_transition e e' x hx_e hx_e'
+                    hn1 hmani he he'
+  rw [h_chain]
+  congr
+
+/-- Transformation law for the metric components under chart change.
+    `g_e(y)(v,w) = g_{e'}(y')(D(φ)_y v, D(φ)_y w)`, where `y = ex`, `y' = e'x`, `φ = e' ∘ e⁻¹`.
+    Assumes `e, e'` are in a `C^n` maximal atlas (`1 ≤ n`). -/
+lemma chartMetric_bilinear_pullback (g : PseudoRiemannianMetric E E M n I)
+    (e e' : PartialHomeomorph M E) (x : M)
+    (hx_e : x ∈ e.source) (hx_e' : x ∈ e'.source)
+    (hmani : IsManifold I n M)
+    (hn1 : 1 ≤ n)
+    (he : e ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (he' : e' ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (v w : E) :
+    let y := e x
+    let y' := e' x
+    let φ := e.symm.trans e'
+    chartMetric g e y v w =
+    (BilinearForm.pullback (chartMetric g e' y') (mfderiv I I φ y)) v w := by
+  intro y y' φ
+  have hy_target : y ∈ e.target := e.mapsTo hx_e
+  have hy'_target : y' ∈ e'.target := e'.mapsTo hx_e'
+  have h_x_e : e.symm y = x := e.left_inv hx_e
+  have h_x_e' : e'.symm y' = x := e'.left_inv hx_e'
+  simp only [chartMetric_apply_of_mem g e hy_target,
+    chartMetric_apply_of_mem g e' hy'_target,
+    BilinearForm.pullback, ContinuousLinearMap.bilinearComp_apply]
+  rw [h_x_e]
+  rw [h_x_e']
+  have h_phi_def : φ = e.symm.trans e' := rfl
+  rw [h_phi_def]
+  exact chartMetric_at_corresponding_points g e e' x hx_e hx_e' hmani hn1 he he' v w
+
+/-- Transformation law for `chartMetric`: `(g_e)_y = ((e'∘e⁻¹)^* (g_{e'})_{e'y})_y`.
+Given a metric `g`, and two charts `e, e'`, the representation of `g` in chart `e` at `y = ex`
+is the pullback by the transition map `φ = e' ∘ e⁻¹` (evaluated at `y`)
+of the representation of `g` in chart `e'` (evaluated at `y' = e'x = φy`).
+Assumes charts are `C^n` compatible, `1 ≤ n`. -/
+theorem chartMetric_coord_change (g : PseudoRiemannianMetric E E M n I)
+    (e e' : PartialHomeomorph M E)
+    (x_pt : M)
+    (hx_e_source : x_pt ∈ e.source) (hx_e'_source : x_pt ∈ e'.source)
+    (hn1 : 1 ≤ n)
+    (hmani_chart : IsManifold I n M)
+    (he_atlas : e ∈ (contDiffGroupoid n I).maximalAtlas M)
+    (he'_atlas : e' ∈ (contDiffGroupoid n I).maximalAtlas M) :
+    chartMetric g e (e x_pt) =
+      BilinearForm.pullback (chartMetric g e' (e' x_pt))
+        (mfderiv I I (e.symm.trans e') (e x_pt)) := by
+  let φ := e.symm.trans e'
+  let y := e x_pt
+  let y' := e' x_pt
+  ext v w
+  exact chartMetric_bilinear_pullback g e e' x_pt hx_e_source hx_e'_source
+                                     hmani_chart hn1 he_atlas he'_atlas v w
diff --git a/PhysLean/Mathematics/Geometry/Metric/PseudoRiemannian/Defs.lean b/PhysLean/Mathematics/Geometry/Metric/PseudoRiemannian/Defs.lean
new file mode 100644
index 000000000..f5661003f
--- /dev/null
+++ b/PhysLean/Mathematics/Geometry/Metric/PseudoRiemannian/Defs.lean
@@ -0,0 +1,567 @@
+/-
+Copyright (c) 2025 Matteo Cipollina. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Matteo Cipollina
+-/
+
+import Mathlib.Analysis.InnerProductSpace.Basic
+import Mathlib.Analysis.RCLike.Lemmas
+import Mathlib.Geometry.Manifold.MFDeriv.Defs
+import Mathlib.LinearAlgebra.BilinearForm.Properties
+import Mathlib.LinearAlgebra.QuadraticForm.Real
+import Mathlib.Topology.LocallyConstant.Basic
+
+/-!
+# Pseudo-Riemannian Metrics on Smooth Manifolds
+
+This file formalizes pseudo-Riemannian metrics on smooth manifolds and establishes their basic
+properties.
+
+A pseudo-Riemannian metric equips a manifold with a smoothly varying, non-degenerate, symmetric
+bilinear form of *constant index* on each tangent space, generalizing the concept of an inner
+product space to curved spaces. The index here refers to `QuadraticForm.negDim`, the dimension
+of a maximal negative definite subspace.
+
+## Main Definitions
+
+* `PseudoRiemannianMetric E H M n I`: A structure representing a `C^n` pseudo-Riemannian metric
+  on a manifold `M` modelled on `(E, H)` with model with corners `I`. It consists of a family
+  of non-degenerate, symmetric, continuous bilinear forms `gₓ` on each tangent space `TₓM`,
+  varying `C^n`-smoothly with `x` and having a locally constant negative dimension (`negDim`).
+  The model space `E` must be finite-dimensional, and the manifold `M` must be `C^{n+1}` smooth.
+
+* `PseudoRiemannianMetric.flatEquiv g x`: The "musical isomorphism" from the tangent space at `x`
+  to its dual space, representing the canonical isomorphism induced by the metric.
+
+* `PseudoRiemannianMetric.sharpEquiv g x`: The inverse of the flat isomorphism, mapping from
+  the dual space back to the tangent space.
+
+* `PseudoRiemannianMetric.toQuadraticForm g x`: The quadratic form `v ↦ gₓ(v, v)` associated
+  with the metric at point `x`.
+-/
+
+section PseudoRiemannianMetric
+
+noncomputable section
+
+universe u v w
+
+open Bundle Set Finset Function Filter Module Topology ContinuousLinearMap
+open scoped Manifold Bundle LinearMap Dual
+
+namespace QuadraticForm
+
+variable {K : Type*} [Field K]
+
+/-- The negative dimension (or index) of a quadratic form is the dimension
+    of a maximal negative definite subspace. -/
+noncomputable def negDim {E : Type*} [AddCommGroup E]
+    [Module ℝ E] [FiniteDimensional ℝ E]
+    (q : QuadraticForm ℝ E) : ℕ := by classical
+  let P : (Fin (finrank ℝ E) → SignType) → Prop := fun w =>
+      QuadraticMap.Equivalent q (QuadraticMap.weightedSumSquares ℝ fun i => (w i : ℝ))
+  let h_exists : ∃ w, P w := QuadraticForm.equivalent_signType_weighted_sum_squared q
+  let w := Classical.choose h_exists
+  exact Finset.card (Finset.filter (fun i => w i = SignType.neg) Finset.univ)
+
+/-- For a standard basis vector in a weighted sum of squares, only one term in the sum
+    is nonzero. -/
+lemma QuadraticMap.weightedSumSquares_basis_vector {E : Type*} [AddCommGroup E]
+    [Module ℝ E] {weights : Fin (finrank ℝ E) → ℝ}
+    {i : Fin (finrank ℝ E)} (v : Fin (finrank ℝ E) → ℝ)
+    (hv : ∀ j, v j = if j = i then 1 else 0) :
+    QuadraticMap.weightedSumSquares ℝ weights v = weights i := by
+  simp only [QuadraticMap.weightedSumSquares_apply]
+  rw [Finset.sum_eq_single i]
+  · simp only [hv i, ↓reduceIte, mul_one, smul_eq_mul]
+  · intro j _ hj
+    simp only [hv j, if_neg hj, mul_zero, smul_eq_mul]
+  · simp only [Finset.mem_univ, not_true_eq_false, smul_eq_mul, mul_eq_zero, or_self,
+    IsEmpty.forall_iff]
+
+/-- When a quadratic form is equivalent to a weighted sum of squares,
+    negative weights correspond to vectors where the form takes negative values. -/
+lemma neg_weight_implies_neg_value {E : Type*} [AddCommGroup E] [Module ℝ E]
+    {q : QuadraticForm ℝ E} {w : Fin (finrank ℝ E) → SignType}
+    (h_equiv : QuadraticMap.Equivalent q (QuadraticMap.weightedSumSquares ℝ fun i => (w i : ℝ)))
+    {i : Fin (finrank ℝ E)} (hi : w i = SignType.neg) :
+    ∃ v : E, v ≠ 0 ∧ q v < 0 := by
+  let f := Classical.choice h_equiv
+  let v_std : Fin (finrank ℝ E) → ℝ := fun j => if j = i then 1 else 0
+  let v := f.symm v_std
+  have hv_ne_zero : v ≠ 0 := by
+    intro h
+    have : f v = f 0 := by rw [h]
+    have : f (f.symm v_std) = f 0 := by rw [← this]
+    have : v_std = 0 := by
+      rw [← f.apply_symm_apply v_std]
+      exact Eq.trans this (map_zero f)
+    have : v_std i = 0 := by rw [this]; rfl
+    simp only [↓reduceIte, one_ne_zero, v_std] at this
+  have hq_neg : q v < 0 := by
+    have heq : q v = QuadraticMap.weightedSumSquares ℝ (fun j => (w j : ℝ)) v_std :=
+      QuadraticMap.IsometryEquiv.map_app f.symm v_std
+    have hw : QuadraticMap.weightedSumSquares ℝ (fun j => (w j : ℝ)) v_std = (w i : ℝ) := by
+      apply QuadraticMap.weightedSumSquares_basis_vector v_std
+      intro j; simp only [v_std]
+    rw [heq, hw]
+    have : (w i : ℝ) = -1 := by simp only [hi, SignType.neg_eq_neg_one, SignType.coe_neg,
+      SignType.coe_one, v_std]
+    rw [this]
+    exact neg_one_lt_zero
+  exact ⟨v, hv_ne_zero, hq_neg⟩
+
+/-- A positive definite quadratic form cannot have any negative weights
+    in its diagonal representation. -/
+lemma posDef_no_neg_weights {E : Type*} [AddCommGroup E] [Module ℝ E]
+    {q : QuadraticForm ℝ E} (hq : q.PosDef)
+    {w : Fin (finrank ℝ E) → SignType}
+    (h_equiv : QuadraticMap.Equivalent q (QuadraticMap.weightedSumSquares ℝ fun i => (w i : ℝ))) :
+    ∀ i, w i ≠ SignType.neg := by
+  intro i hi
+  obtain ⟨v, hv_ne_zero, hq_neg⟩ := QuadraticForm.neg_weight_implies_neg_value h_equiv hi
+  have hq_pos : 0 < q v := hq v hv_ne_zero
+  exact lt_asymm hq_neg hq_pos
+
+/-- For a positive definite quadratic form, the negative dimension (index) is zero. -/
+theorem rankNeg_eq_zero {E : Type*} [AddCommGroup E]
+    [Module ℝ E] [FiniteDimensional ℝ E] {q : QuadraticForm ℝ E} (hq : q.PosDef) :
+    q.negDim = 0 := by
+  haveI : Invertible (2 : ℝ) := inferInstance
+  unfold QuadraticForm.negDim
+  have h_exists := equivalent_signType_weighted_sum_squared q
+  let w := Classical.choose h_exists
+  have h_equiv : QuadraticMap.Equivalent q
+      (QuadraticMap.weightedSumSquares ℝ fun i => (w i : ℝ)) :=
+    Classical.choose_spec h_exists
+  have h_no_neg : ∀ i, w i ≠ SignType.neg :=
+    QuadraticForm.posDef_no_neg_weights hq h_equiv
+  simp [Finset.card_eq_zero, Finset.filter_eq_empty_iff, h_no_neg]
+  exact fun ⦃x⦄ => h_no_neg x
+
+end QuadraticForm
+
+/-- Helper function to convert the metric tensor at `x` to a quadratic form. -/
+private def pseudoRiemannianMetricValToQuadraticForm
+    {E : Type v} [NormedAddCommGroup E] [NormedSpace ℝ E]
+    {H : Type w} [TopologicalSpace H]
+    {M : Type w} [TopologicalSpace M] [ChartedSpace H M]
+    {I : ModelWithCorners ℝ E H}
+    (val : ∀ (x : M), TangentSpace I x →L[ℝ] (TangentSpace I x →L[ℝ] ℝ))
+    (symm : ∀ (x : M) (v w : TangentSpace I x), (val x v) w = (val x w) v)
+    (x : M) : QuadraticForm ℝ (TangentSpace I x) where
+  toFun v := val x v v
+  toFun_smul a v := by
+    simp only [ContinuousLinearMap.map_smul, ContinuousLinearMap.smul_apply, smul_smul, pow_two]
+  exists_companion' :=
+      ⟨ LinearMap.mk₂ ℝ (fun v y => val x v y + val x y v)
+        (fun v₁ v₂ y => by simp only [map_add, add_apply]; ring)
+        (fun a v y => by simp only [map_smul, smul_apply]; ring_nf; exact Eq.symm (smul_add ..))
+        (fun v y₁ y₂ => by simp only [map_add, add_apply]; ring)
+        (fun a v y => by simp only [map_smul, smul_apply]; ring_nf; exact Eq.symm (smul_add ..)),
+            by
+              intro v y
+              simp only [LinearMap.mk₂_apply, ContinuousLinearMap.map_add,
+                ContinuousLinearMap.add_apply, symm x]
+              ring⟩
+
+/-- A pseudo-Riemannian metric of smoothness class `C^n` on a manifold `M` modelled on `(E, H)`
+with model `I`. This structure defines a smoothly varying, non-degenerate, symmetric,
+continuous bilinear form `gₓ` of constant negative dimension on the tangent space `TₓM`
+at each point `x`. Requires `M` to be `C^{n+1}` smooth.-/
+@[ext]
+structure PseudoRiemannianMetric
+    (E : Type v) (H : Type w) (M : Type w) (n : WithTop ℕ∞)
+    [inst_norm_grp_E : NormedAddCommGroup E]
+    [inst_norm_sp_E : NormedSpace ℝ E]
+    [inst_top_H : TopologicalSpace H]
+    [inst_top_M : TopologicalSpace M]
+    [inst_chart_M : ChartedSpace H M]
+    [inst_chart_E : ChartedSpace H E]
+    (I : ModelWithCorners ℝ E H)
+    [inst_mani : IsManifold I (n + 1) M]
+    [inst_tangent_findim : ∀ (x : M), FiniteDimensional ℝ (TangentSpace I x)] :
+      Type (max u v w) where
+  /-- The metric tensor at each point `x : M`, represented as a continuous linear map
+      `TₓM →L[ℝ] (TₓM →L[ℝ] ℝ)`. Applying it twice, `(val x v) w`, yields `gₓ(v, w)`. -/
+  val : ∀ (x : M), TangentSpace I x →L[ℝ] (TangentSpace I x →L[ℝ] ℝ)
+  /-- The metric is symmetric: `gₓ(v, w) = gₓ(w, v)`. -/
+  symm : ∀ (x : M) (v w : TangentSpace I x), (val x v) w = (val x w) v
+  /-- The metric is non-degenerate: if `gₓ(v, w) = 0` for all `w`, then `v = 0`. -/
+  nondegenerate : ∀ (x : M) (v : TangentSpace I x), (∀ w : TangentSpace I x,
+    (val x v) w = 0) → v = 0
+  /-- The metric varies smoothly: Expressed in local coordinates via any chart `e`, the function
+      `y ↦ g_{e.symm y}(mfderiv I I e.symm y v, mfderiv I I e.symm y w)` is `C^n` smooth on the
+      chart's target `e.target` for any constant vectors `v, w` in the model space `E`. -/
+  smooth_in_charts' : ∀ (x₀ : M) (v w : E),
+    let e := extChartAt I x₀
+    ContDiffWithinAt ℝ n
+      (fun y => val (e.symm y) (mfderiv I I e.symm y v) (mfderiv I I e.symm y w))
+      (e.target) (e x₀)
+  /-- The negative dimension (`QuadraticForm.negDim`) of the metric's quadratic form is
+      locally constant. On a connected manifold, this implies it is constant globally. -/
+  negDim_isLocallyConstant :
+    IsLocallyConstant (fun x : M =>
+      have : FiniteDimensional ℝ (TangentSpace I x) := inferInstance
+      (pseudoRiemannianMetricValToQuadraticForm val symm x).negDim)
+
+namespace PseudoRiemannianMetric
+
+variable {E : Type v} {H : Type w} {M : Type w} {n : WithTop ℕ∞}
+variable [NormedAddCommGroup E] [NormedSpace ℝ E]
+variable [TopologicalSpace H] [TopologicalSpace M] [ChartedSpace H M] [ChartedSpace H E]
+variable {I : ModelWithCorners ℝ E H}
+variable [IsManifold I (n + 1) M]
+variable [inst_tangent_findim : ∀ (x : M), FiniteDimensional ℝ (TangentSpace I x)]
+variable {g : PseudoRiemannianMetric E H M n I}
+
+/-- Convert the metric's continuous linear map representation `val x` to the algebraic
+    `LinearMap.BilinForm`. -/
+def toBilinForm (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    LinearMap.BilinForm ℝ (TangentSpace I x) where
+  toFun := λ v => { toFun := λ w => g.val x v w,
+                    map_add' := λ w₁ w₂ => by
+                      simp only [ContinuousLinearMap.map_add, ContinuousLinearMap.add_apply],
+                    map_smul' := λ c w => by
+                      simp only [map_smul, smul_eq_mul, RingHom.id_apply] }
+  map_add' := λ v₁ v₂ => by
+    ext w
+    simp only [map_add, add_apply, LinearMap.coe_mk, AddHom.coe_mk, LinearMap.add_apply]
+  map_smul' := λ c v => by
+    ext w
+    simp only [map_smul, coe_smul', Pi.smul_apply, smul_eq_mul, LinearMap.coe_mk, AddHom.coe_mk,
+      RingHom.id_apply, LinearMap.smul_apply]
+
+/-- Convert a pseudo-Riemannian metric at a point `x` to a quadratic form `v ↦ gₓ(v, v)`. -/
+def toQuadraticForm (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    QuadraticForm ℝ (TangentSpace I x) :=
+  pseudoRiemannianMetricValToQuadraticForm g.val g.symm x
+
+-- Coercion from PseudoRiemannianMetric to its function representation.
+instance coeFunInst : CoeFun (PseudoRiemannianMetric E H M n I)
+        (fun _ => ∀ x : M, TangentSpace I x →L[ℝ] (TangentSpace I x →L[ℝ] ℝ)) where
+   coe g := g.val
+
+@[simp]
+lemma toBilinForm_apply (g : PseudoRiemannianMetric E H M n I) (x : M)
+    (v w : TangentSpace I x) :
+  toBilinForm g x v w = g.val x v w := rfl
+
+@[simp]
+lemma toQuadraticForm_apply (g : PseudoRiemannianMetric E H M n I) (x : M)
+    (v : TangentSpace I x) :
+  toQuadraticForm g x v = g.val x v v := rfl
+
+@[simp]
+lemma toBilinForm_isSymm (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    (toBilinForm g x).IsSymm := by
+  intro v w; simp only [toBilinForm_apply]; exact g.symm x v w
+
+@[simp]
+lemma toBilinForm_nondegenerate (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    (toBilinForm g x).Nondegenerate := by
+  intro v hv; simp_rw [toBilinForm_apply] at hv; exact g.nondegenerate x v hv
+
+/-- The "musical" isomorphism (index lowering) from the tangent space to its dual,
+    induced by a pseudo-Riemannian metric. -/
+def flat (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    TangentSpace I x →ₗ[ℝ] (TangentSpace I x →L[ℝ] ℝ) :=
+  { toFun := λ v => g.val x v,
+    map_add' := λ v w => by simp only [ContinuousLinearMap.map_add],
+    map_smul' := λ a v => by simp only [ContinuousLinearMap.map_smul]; rfl }
+
+@[simp]
+lemma flat_apply (g : PseudoRiemannianMetric E H M n I) (x : M) (v w : TangentSpace I x) :
+    (flat g x v) w = g.val x v w := by rfl
+
+/-- The musical isomorphism as a continuous linear map. -/
+def flatL (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    TangentSpace I x →L[ℝ] (TangentSpace I x →L[ℝ] ℝ) where
+  toFun := λ v => g.val x v
+  map_add' := λ v w => by simp only [ContinuousLinearMap.map_add]
+  map_smul' := λ a v => by simp only [ContinuousLinearMap.map_smul]; rfl
+  cont := ContinuousLinearMap.continuous (g.val x)
+
+@[simp]
+lemma flatL_apply (g : PseudoRiemannianMetric E H M n I) (x : M) (v w : TangentSpace I x) :
+    (flatL g x v) w = g.val x v w := rfl
+
+@[simp]
+lemma flat_inj (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    Function.Injective (flat g x) := by
+  rw [← LinearMap.ker_eq_bot]; apply LinearMap.ker_eq_bot'.mpr
+  intro v hv; apply g.nondegenerate x v; intro w; exact DFunLike.congr_fun hv w
+
+@[simp]
+lemma flatL_inj (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    Function.Injective (flatL g x) :=
+  flat_inj g x
+
+/-- In a finite-dimensional normed space, the continuous dual is linearly equivalent
+    to the algebraic dual. -/
+def ContinuousLinearMap.equivModuleDual (𝕜 E : Type*) [NontriviallyNormedField 𝕜] [CompleteSpace 𝕜]
+    [NormedAddCommGroup E] [NormedSpace 𝕜 E] [FiniteDimensional 𝕜 E] :
+    (E →L[𝕜] 𝕜) ≃ₗ[𝕜] Module.Dual 𝕜 E where
+  toFun f := f.toLinearMap
+  invFun φ :=
+    { toFun := φ
+      map_add' := LinearMap.map_add φ
+      map_smul' := LinearMap.map_smul φ
+      cont := φ.continuous_of_finiteDimensional }
+  map_add' f g := rfl
+  map_smul' c f := rfl
+  left_inv f := by ext; rfl
+  right_inv φ := rfl
+
+/-- For a finite-dimensional normed space, the dimension of the continuous dual
+equals the dimension of the original space. -/
+lemma finrank_continuousDual_eq_finrank {𝕜 E : Type*} [NontriviallyNormedField 𝕜]
+    [NormedAddCommGroup E] [NormedSpace 𝕜 E] [FiniteDimensional 𝕜 E] [CompleteSpace 𝕜] :
+    finrank 𝕜 (E →L[𝕜] 𝕜) = finrank 𝕜 E := by
+  have h1 : (E →L[𝕜] 𝕜) ≃ₗ[𝕜] Module.Dual 𝕜 E := by
+    exact ContinuousLinearMap.equivModuleDual 𝕜 E
+  have h2 : finrank 𝕜 (Module.Dual 𝕜 E) = finrank 𝕜 E := by
+    exact finrank_linearMap_self 𝕜 𝕜 E
+  rw [LinearEquiv.finrank_eq h1, h2]
+
+@[simp]
+lemma flatL_surj
+    (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    Function.Surjective (g.flatL x) := by
+  haveI : FiniteDimensional ℝ (TangentSpace I x) := inst_tangent_findim x
+  have h_finrank_eq : finrank ℝ (TangentSpace I x) = finrank ℝ (TangentSpace I x →L[ℝ] ℝ) := by
+    have h_dual_eq : finrank ℝ (TangentSpace I x →L[ℝ] ℝ) = finrank ℝ (Module.Dual ℝ
+    (TangentSpace I x)) := by
+      let to_dual : (TangentSpace I x →L[ℝ] ℝ) → Module.Dual ℝ (TangentSpace I x) :=
+        fun f => f.toLinearMap
+      let from_dual : Module.Dual ℝ (TangentSpace I x) → (TangentSpace I x →L[ℝ] ℝ) := fun f =>
+        ContinuousLinearMap.mk f (by
+          apply LinearMap.continuous_of_finiteDimensional)
+      let equiv : (TangentSpace I x →L[ℝ] ℝ) ≃ₗ[ℝ] Module.Dual ℝ (TangentSpace I x) :=
+      { toFun := to_dual,
+        invFun := from_dual,
+        map_add' := fun f g => by
+          ext v; unfold to_dual; simp only [LinearMap.add_apply]; rfl,
+        map_smul' := fun c f => by
+          ext v; unfold to_dual; simp only [LinearMap.smul_apply]; rfl,
+        left_inv := fun f => by
+          ext v; unfold to_dual from_dual; simp,
+        right_inv := fun f => by
+          ext v; unfold to_dual from_dual; simp }
+      exact LinearEquiv.finrank_eq equiv
+    rw [h_dual_eq, ← Subspace.dual_finrank_eq]
+  exact (LinearMap.injective_iff_surjective_of_finrank_eq_finrank h_finrank_eq).mp (flatL_inj g x)
+
+/-- The "musical" isomorphism (index lowering) from the tangent space to its dual,
+    as a continuous linear equivalence. -/
+def flatEquiv
+    (g : PseudoRiemannianMetric E H M n I)
+    (x : M) :
+    TangentSpace I x ≃L[ℝ] (TangentSpace I x →L[ℝ] ℝ) :=
+  LinearEquiv.toContinuousLinearEquiv
+    (LinearEquiv.ofBijective
+      ((g.flatL x).toLinearMap)
+      ⟨g.flatL_inj x, g.flatL_surj x⟩)
+
+lemma coe_flatEquiv
+    (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    (g.flatEquiv x : TangentSpace I x →ₗ[ℝ] (TangentSpace I x →L[ℝ] ℝ)) = g.flatL x := rfl
+
+@[simp]
+lemma flatEquiv_apply
+    (g : PseudoRiemannianMetric E H M n I) (x : M) (v w : TangentSpace I x) :
+    (g.flatEquiv x v) w = g.val x v w := rfl
+
+/-- The "musical" isomorphism (index raising) from the dual of the tangent space to the
+    tangent space, induced by a pseudo-Riemannian metric. This is the inverse of `flatEquiv`. -/
+def sharpEquiv
+    (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    (TangentSpace I x →L[ℝ] ℝ) ≃L[ℝ] TangentSpace I x :=
+  (g.flatEquiv x).symm
+
+/-- The index raising map `sharp` as a continuous linear map. -/
+def sharpL
+    (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    (TangentSpace I x →L[ℝ] ℝ) →L[ℝ] TangentSpace I x := (g.sharpEquiv x).toContinuousLinearMap
+
+lemma sharpL_eq_toContinuousLinearMap
+    (g : PseudoRiemannianMetric E H M n I) (x : M) :
+  g.sharpL x = (g.sharpEquiv x).toContinuousLinearMap := rfl
+
+lemma coe_sharpEquiv
+    (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    (g.sharpEquiv x : (TangentSpace I x →L[ℝ] ℝ) →L[ℝ] TangentSpace I x) = g.sharpL x := rfl
+
+/-- The index raising map `sharp` as a linear map. -/
+noncomputable def sharp
+    (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    (TangentSpace I x →L[ℝ] ℝ) →ₗ[ℝ] TangentSpace I x := (g.sharpEquiv x).toLinearEquiv.toLinearMap
+
+@[simp]
+lemma sharpL_apply_flatL
+    (g : PseudoRiemannianMetric E H M n I) (x : M) (v : TangentSpace I x) :
+    g.sharpL x (g.flatL x v) = v :=
+  (g.flatEquiv x).left_inv v
+
+@[simp]
+lemma flatL_apply_sharpL
+    (g : PseudoRiemannianMetric E H M n I) (x : M) (ω : TangentSpace I x →L[ℝ] ℝ) :
+    g.flatL x (g.sharpL x ω) = ω := (g.flatEquiv x).right_inv ω
+
+/-- Applying `sharp` then `flat` recovers the original covector. -/
+@[simp]
+lemma flat_sharp_apply
+    (g : PseudoRiemannianMetric E H M n I) (x : M) (ω : TangentSpace I x →L[ℝ] ℝ) :
+    g.flat x (g.sharp x ω) = ω := by
+  have := flatL_apply_sharpL g x ω
+  simp only [sharp, sharpL, flat, flatL, coe_flatEquiv]; simp only [coe_sharpEquiv,
+             ContinuousLinearEquiv.coe_coe, LinearEquiv.coe_coe] at this ⊢
+  exact this
+
+@[simp]
+lemma sharp_flat_apply
+    (g : PseudoRiemannianMetric E H M n I) (x : M) (v : TangentSpace I x) :
+    g.sharp x (g.flat x v) = v := by
+  have := sharpL_apply_flatL g x v
+  simp only [sharp, sharpL, flat, flatL]; simp only [coe_flatEquiv, coe_sharpEquiv,
+             ContinuousLinearEquiv.coe_coe, LinearEquiv.coe_coe] at this ⊢
+  exact this
+
+/-- The metric evaluated at `sharp ω₁` and `sharp ω₂`. -/
+@[simp]
+lemma apply_sharp_sharp
+    (g : PseudoRiemannianMetric E H M n I) (x : M) (ω₁ ω₂ : TangentSpace I x →L[ℝ] ℝ) :
+    g.val x (g.sharpL x ω₁) (g.sharpL x ω₂) = ω₁ (g.sharpL x ω₂) := by
+  rw [← flatL_apply g x (g.sharpL x ω₁)]
+  rw [flatL_apply_sharpL g x ω₁]
+
+/-- The metric evaluated at `v` and `sharp ω`. -/
+lemma apply_vec_sharp
+    (g : PseudoRiemannianMetric E H M n I) (x : M) (v : TangentSpace I x)
+    (ω : TangentSpace I x →L[ℝ] ℝ) :
+    g.val x v (g.sharpL x ω) = ω v := by
+  rw [g.symm x v (g.sharpL x ω)]
+  rw [← flatL_apply g x (g.sharpL x ω)]
+  rw [flatL_apply_sharpL g x ω]
+
+section Cotangent
+
+variable {E : Type v} {H : Type w} {M : Type w} {n : WithTop ℕ∞}
+variable [NormedAddCommGroup E] [NormedSpace ℝ E]
+variable [TopologicalSpace H] [TopologicalSpace M] [ChartedSpace H M] [ChartedSpace H E]
+variable {I : ModelWithCorners ℝ E H}
+variable [IsManifold I (n + 1) M]
+variable [inst_tangent_findim : ∀ (x : M), FiniteDimensional ℝ (TangentSpace I x)]
+
+/-- The value of the induced metric on the cotangent space at point `x`. -/
+noncomputable def cotangentMetricVal (g : PseudoRiemannianMetric E H M n I) (x : M)
+    (ω₁ ω₂ : TangentSpace I x →L[ℝ] ℝ) : ℝ :=
+  g.val x (g.sharpL x ω₁) (g.sharpL x ω₂)
+
+@[simp]
+lemma cotangentMetricVal_eq_apply_sharp (g : PseudoRiemannianMetric E H M n I) (x : M)
+    (ω₁ ω₂ : TangentSpace I x →L[ℝ] ℝ) :
+  cotangentMetricVal g x ω₁ ω₂ = ω₁ (g.sharpL x ω₂) := by
+  rw [cotangentMetricVal, apply_sharp_sharp]
+
+lemma cotangentMetricVal_symm (g : PseudoRiemannianMetric E H M n I) (x : M)
+    (ω₁ ω₂ : TangentSpace I x →L[ℝ] ℝ) :
+  cotangentMetricVal g x ω₁ ω₂ = cotangentMetricVal g x ω₂ ω₁ := by
+  unfold cotangentMetricVal
+  rw [g.symm x (g.sharpL x ω₁) (g.sharpL x ω₂)]
+
+/-- The induced metric on the cotangent space at point `x` as a bilinear form. -/
+noncomputable def cotangentToBilinForm (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    LinearMap.BilinForm ℝ (TangentSpace I x →L[ℝ] ℝ) where
+  toFun ω₁ := { toFun := λ ω₂ => cotangentMetricVal g x ω₁ ω₂,
+                map_add' := λ ω₂ ω₃ => by
+                  simp only [cotangentMetricVal,
+                    ContinuousLinearMap.map_add,
+                    ContinuousLinearMap.add_apply],
+                map_smul' := λ c ω₂ => by
+                  simp only [cotangentMetricVal,
+                    map_smul, smul_eq_mul, RingHom.id_apply] }
+  map_add' := λ ω₁ ω₂ => by
+    ext ω₃
+    simp only [cotangentMetricVal,
+      ContinuousLinearMap.map_add,
+      ContinuousLinearMap.add_apply,
+      LinearMap.coe_mk, AddHom.coe_mk, LinearMap.add_apply]
+  map_smul' := λ c ω₁ => by
+    ext ω₂
+    simp only [cotangentMetricVal,
+      ContinuousLinearMap.map_smul,
+      ContinuousLinearMap.smul_apply,
+      LinearMap.coe_mk, AddHom.coe_mk,
+      RingHom.id_apply, LinearMap.smul_apply]
+
+/-- The induced metric on the cotangent space at point `x` as a quadratic form. -/
+noncomputable def cotangentToQuadraticForm (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    QuadraticForm ℝ (TangentSpace I x →L[ℝ] ℝ) where
+  toFun ω := cotangentMetricVal g x ω ω
+  toFun_smul a ω := by
+    simp only [cotangentMetricVal,
+      ContinuousLinearMap.map_smul,
+      ContinuousLinearMap.smul_apply,
+      smul_smul, pow_two]
+  exists_companion' :=
+      ⟨ LinearMap.mk₂ ℝ (fun ω₁ ω₂ =>
+          cotangentMetricVal g x ω₁ ω₂ + cotangentMetricVal g x ω₂ ω₁)
+        (fun ω₁ ω₂ ω₃ => by simp only [cotangentMetricVal, map_add, add_apply]; ring)
+        (fun a ω₁ ω₂ => by
+          simp only [cotangentMetricVal, map_smul, smul_apply];
+          ring_nf; exact Eq.symm (smul_add ..))
+        (fun ω₁ ω₂ ω₃ => by
+          simp only [cotangentMetricVal, map_add, add_apply]; ring)
+        (fun a ω₁ ω₂ => by
+          simp only [cotangentMetricVal, map_smul, smul_apply]; ring_nf;
+          exact Eq.symm (smul_add ..)),
+          by
+            intro ω₁ ω₂
+            simp only [LinearMap.mk₂_apply, cotangentMetricVal]
+            simp only [ContinuousLinearMap.map_add, ContinuousLinearMap.add_apply]
+            ring⟩
+
+@[simp]
+lemma cotangentToBilinForm_apply (g : PseudoRiemannianMetric E H M n I) (x : M)
+    (ω₁ ω₂ : TangentSpace I x →L[ℝ] ℝ) :
+  cotangentToBilinForm g x ω₁ ω₂ = cotangentMetricVal g x ω₁ ω₂ := rfl
+
+@[simp]
+lemma cotangentToQuadraticForm_apply (g : PseudoRiemannianMetric E H M n I) (x : M)
+    (ω : TangentSpace I x →L[ℝ] ℝ) :
+  cotangentToQuadraticForm g x ω = cotangentMetricVal g x ω ω := rfl
+
+@[simp]
+lemma cotangentToBilinForm_isSymm (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    (cotangentToBilinForm g x).IsSymm := by
+  intro ω₁ ω₂; simp only [cotangentToBilinForm_apply]; exact cotangentMetricVal_symm g x ω₁ ω₂
+
+/-- The cotangent metric is non-degenerate: if `cotangentMetricVal g x ω v = 0` for all `v`,
+    then `ω = 0`. -/
+lemma cotangentMetricVal_nondegenerate (g : PseudoRiemannianMetric E H M n I) (x : M)
+    (ω : TangentSpace I x →L[ℝ] ℝ) (h : ∀ v : TangentSpace I x →L[ℝ] ℝ,
+      cotangentMetricVal g x ω v = 0) :
+    ω = 0 := by
+  apply ContinuousLinearMap.ext
+  intro v
+  have h_forall : ∀ w : TangentSpace I x, ω w = 0 := by
+    intro w
+    let ω' : TangentSpace I x →L[ℝ] ℝ := g.flatL x w
+    have this : g.sharpL x ω' = w := by
+      simp only [ω', sharpL_apply_flatL]
+    have h_apply : cotangentMetricVal g x ω ω' = 0 := h ω'
+    simp only [cotangentMetricVal_eq_apply_sharp] at h_apply
+    rw [this] at h_apply
+    exact h_apply
+  exact h_forall v
+
+@[simp]
+lemma cotangentToBilinForm_nondegenerate (g : PseudoRiemannianMetric E H M n I) (x : M) :
+    (cotangentToBilinForm g x).Nondegenerate := by
+  intro ω hω
+  apply cotangentMetricVal_nondegenerate g x ω
+  intro v
+  exact hω v
+
+end Cotangent
diff --git a/PhysLean/Mathematics/Geometry/Metric/Riemannian/Defs.lean b/PhysLean/Mathematics/Geometry/Metric/Riemannian/Defs.lean
new file mode 100644
index 000000000..2974fb6af
--- /dev/null
+++ b/PhysLean/Mathematics/Geometry/Metric/Riemannian/Defs.lean
@@ -0,0 +1,204 @@
+/-
+Copyright (c) 2025 Matteo Cipollina. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Matteo Cipollina
+-/
+
+import Mathlib.Algebra.Lie.OfAssociative
+import Mathlib.MeasureTheory.Integral.IntervalIntegral.Basic
+import Mathlib.MeasureTheory.Integral.IntervalIntegral.FundThmCalculus
+import Mathlib.MeasureTheory.Integral.IntervalIntegral.IntegrationByParts
+import PhysLean.Mathematics.Geometry.Metric.PseudoRiemannian.Chart
+
+/-!
+# Riemannian Metric Definitions
+
+This module defines the Riemannian metric, building on pseudo-Riemannian metrics.
+-/
+namespace PseudoRiemannianMetric
+section RiemannianMetric
+
+open Bundle Set Finset Function Filter Module Topology ContinuousLinearMap
+open scoped Manifold Bundle LinearMap Dual
+open PseudoRiemannianMetric InnerProductSpace
+
+noncomputable section
+
+variable {E : Type v} [NormedAddCommGroup E] [NormedSpace ℝ E]
+variable {H : Type w} [TopologicalSpace H]
+variable {M : Type w} [TopologicalSpace M] [ChartedSpace H M] [ChartedSpace H E]
+variable {I : ModelWithCorners ℝ E H} {n : ℕ∞}
+
+/-- A Riemannian metric of smoothness class `n` on a manifold `M` over `ℝ`.
+    This extends a pseudo-Riemannian metric with the positive definiteness condition. -/
+@[ext]
+structure RiemannianMetric
+  (I : ModelWithCorners ℝ E H) (n : ℕ∞) (M : Type w)
+  [TopologicalSpace M] [ChartedSpace H M] [IsManifold I (n + 1) M]
+  [inst_tangent_findim : ∀ (x : M), FiniteDimensional ℝ (TangentSpace I x)] : Type _ where
+  /-- The underlying pseudo-Riemannian metric. -/
+  toPseudoRiemannianMetric : PseudoRiemannianMetric E H M (n) I
+  /-- `gₓ(v, v) > 0` for all nonzero `v`. -/
+  pos_def' : ∀ x v, v ≠ 0 → toPseudoRiemannianMetric.val x v v > 0
+
+namespace RiemannianMetric
+
+variable {I : ModelWithCorners ℝ E H} {n : ℕ∞} {M : Type w}
+variable [TopologicalSpace M] [ChartedSpace H M] [IsManifold I (n + 1) M]
+variable [inst_tangent_findim : ∀ (x : M), FiniteDimensional ℝ (TangentSpace I x)]
+
+/-- Coercion from RiemannianMetric to its underlying PseudoRiemannianMetric. -/
+instance : Coe (RiemannianMetric I n M) (PseudoRiemannianMetric E H M (n) I) where
+  coe g := g.toPseudoRiemannianMetric
+
+@[simp]
+lemma pos_def (g : RiemannianMetric I n M) (x : M) (v : TangentSpace I x)
+    (hv : v ≠ 0) :
+  (g.toPseudoRiemannianMetric.val x v) v > 0 := g.pos_def' x v hv
+
+/-- The inverse of the musical isomorphism (index raising), which is an isomorphism
+    in the Riemannian case. This is well-defined because the metric is positive definite. -/
+def sharp
+    (g : RiemannianMetric I n M) (x : M) :
+    (TangentSpace I x →L[ℝ] ℝ) ≃L[ℝ] TangentSpace I x :=
+  g.toPseudoRiemannianMetric.sharpEquiv x
+
+/-- Express a Riemannian metric at a point as a quadratic form. -/
+abbrev toQuadraticForm (g : RiemannianMetric I n M) (x : M) :
+    QuadraticForm ℝ (TangentSpace I x) :=
+  g.toPseudoRiemannianMetric.toQuadraticForm x
+
+/-- The quadratic form associated with a Riemannian metric is positive definite. -/
+@[simp]
+lemma toQuadraticForm_posDef (g : RiemannianMetric I n M) (x : M) :
+    (g.toQuadraticForm x).PosDef :=
+  λ v hv => g.pos_def x v hv
+
+/-- The application of a Riemannian metric's quadratic form to a vector. -/
+@[simp]
+lemma toQuadraticForm_apply (g : RiemannianMetric I n M) (x : M)
+    (v : TangentSpace I x) :
+    g.toQuadraticForm x v = g.toPseudoRiemannianMetric.val x v v := by
+  simp only [toQuadraticForm, PseudoRiemannianMetric.toQuadraticForm_apply]
+
+lemma riemannian_metric_negDim_zero (g : RiemannianMetric I n M) (x : M) :
+    (g.toQuadraticForm x).negDim = 0 := by
+  apply QuadraticForm.rankNeg_eq_zero
+  exact g.toQuadraticForm_posDef x
+
+/-- The inner product on the tangent space at point `x` induced by the Riemannian metric `g`. -/
+def inner (g : RiemannianMetric I n M) (x : M) (v w : TangentSpace I x) : ℝ :=
+  g.toPseudoRiemannianMetric.val x v w
+
+@[simp]
+lemma inner_apply (g : RiemannianMetric I n M) (x : M) (v w : TangentSpace I x) :
+    inner g x v w = g.toPseudoRiemannianMetric.val x v w := rfl
+
+variable (g : RiemannianMetric I n M) (x : M)
+
+/-- Pointwise symmetry of the inner product. -/
+lemma inner_symm (v w : TangentSpace I x) :
+    g.inner x v w = g.inner x w v := by
+  simp only [inner_apply]
+  exact g.toPseudoRiemannianMetric.symm x v w
+
+/-- The inner product space core for the tangent space at a point, derived from the
+Riemannian metric. -/
+noncomputable def tangentInnerCore (g : RiemannianMetric I n M) (x : M) :
+    InnerProductSpace.Core ℝ (TangentSpace I x) where
+  inner := λ v w => g.inner x v w
+  conj_inner_symm := λ v w => by
+    simp only [inner_apply, conj_trivial]
+    exact g.toPseudoRiemannianMetric.symm x w v
+  re_inner_nonneg := λ v => by
+    simp only [inner_apply, RCLike.re_to_real]
+    by_cases hv : v = 0
+    · simp [hv, inner_apply, map_zero]
+    · exact le_of_lt (g.pos_def x v hv)
+  add_left := λ u v w => by
+    simp only [inner_apply, map_add, ContinuousLinearMap.add_apply]
+  smul_left := λ r u v => by
+    simp only [inner_apply, map_smul, conj_trivial]
+    rfl
+  definite := fun v (h_inner_zero : g.inner x v v = 0) => by
+    by_contra h_v_ne_zero
+    have h_pos : g.inner x v v > 0 := g.pos_def x v h_v_ne_zero
+    linarith [h_inner_zero, h_pos]
+
+/- We Define NormedAddCommGroup and InnerProductSpace structures locally.
+ We DO NOT mark them as instance.-/
+
+/-- Creates a `NormedAddCommGroup` structure on the tangent space at a point,
+    derived from a Riemannian metric. -/
+noncomputable def TangentSpace.metricNormedAddCommGroup (g : RiemannianMetric I n M) (x : M) :
+    NormedAddCommGroup (TangentSpace I x) :=
+  @InnerProductSpace.Core.toNormedAddCommGroup ℝ (TangentSpace I x) _ _ _ (tangentInnerCore g x)
+
+/-- Creates an `InnerProductSpace` structure on the tangent space at a point,
+    derived from a Riemannian metric. Alternative implementation using `letI`. -/
+noncomputable def TangentSpace.metricInnerProductSpace' (g : RiemannianMetric I n M) (x : M) :
+    letI := TangentSpace.metricNormedAddCommGroup g x
+    InnerProductSpace ℝ (TangentSpace I x) :=
+  InnerProductSpace.ofCore (tangentInnerCore g x)
+
+/-- Creates an `InnerProductSpace` structure on the tangent space at a point,
+    derived from a Riemannian metric. -/
+noncomputable def TangentSpace.metricInnerProductSpace (g : RiemannianMetric I n M) (x : M) :
+    let _ := TangentSpace.metricNormedAddCommGroup g x
+    InnerProductSpace ℝ (TangentSpace I x) :=
+  let inner_core := tangentInnerCore g x
+  let _ := TangentSpace.metricNormedAddCommGroup g x
+  InnerProductSpace.ofCore inner_core
+
+/-- The norm on a tangent space induced by a Riemannian metric, defined as the square root
+    of the inner product of a vector with itself. -/
+noncomputable def norm (g : RiemannianMetric I n M) (x : M) (v : TangentSpace I x) : ℝ :=
+  Real.sqrt (g.inner x v v)
+
+-- Example using the norm function
+example (g : RiemannianMetric I n M) (x : M) (v : TangentSpace I x) :
+    norm g x v ≥ 0 := Real.sqrt_nonneg _
+
+-- Example showing how to use the metric inner product space
+example (g : RiemannianMetric I n M) (x : M) (v w : TangentSpace I x) :
+    (TangentSpace.metricInnerProductSpace g x).inner v w = g.inner x v w := by
+  letI := TangentSpace.metricInnerProductSpace g x
+  rfl
+
+/-- Helper function to compute the norm on a tangent space from a Riemannian metric,
+    using the underlying `NormedAddCommGroup` structure. -/
+noncomputable def norm' (g : RiemannianMetric I n M) (x : M) (v : TangentSpace I x) : ℝ :=
+  let normed_group := TangentSpace.metricNormedAddCommGroup g x
+  @Norm.norm (TangentSpace I x) (@NormedAddCommGroup.toNorm (TangentSpace I x) normed_group) v
+
+-- Example: Using a custom norm function instead of the notation
+example (g : RiemannianMetric I n M) (x : M) (v : TangentSpace I x) :
+    norm g x v ≥ 0 := by
+  unfold norm
+  apply Real.sqrt_nonneg
+
+example (g : RiemannianMetric I n M) (x : M) (v : TangentSpace I x) : ℝ :=
+  letI := TangentSpace.metricNormedAddCommGroup g x
+  ‖v‖
+
+example (g : RiemannianMetric I n M) (x : M) (v : TangentSpace I x) : ℝ :=
+  let normed_group := TangentSpace.metricNormedAddCommGroup g x
+  @Norm.norm (TangentSpace I x) (@NormedAddCommGroup.toNorm (TangentSpace I x) normed_group) v
+
+lemma norm_eq_norm_of_metricNormedAddCommGroup (g : RiemannianMetric I n M) (x : M)
+    (v : TangentSpace I x) : norm g x v = @Norm.norm (TangentSpace I x)
+    (@NormedAddCommGroup.toNorm _ (TangentSpace.metricNormedAddCommGroup g x)) v := by
+  unfold norm
+  let normed_group := TangentSpace.metricNormedAddCommGroup g x
+  unfold TangentSpace.metricNormedAddCommGroup
+  simp only [inner_apply]
+  rfl
+
+/-- Calculates the length of a curve `γ` between parameters `t₀` and `t₁`
+    using the Riemannian metric `g`. This is defined as the integral of the norm of
+    the tangent vector along the curve. -/
+def curveLength (g : RiemannianMetric I n M) (γ : ℝ → M) (t₀ t₁ : ℝ)
+    {IDE : ModelWithCorners ℝ ℝ ℝ}[ChartedSpace ℝ ℝ] : ℝ :=
+  ∫ t in t₀..t₁, norm g (γ t) ((mfderiv IDE I γ t) ((1 : ℝ) : TangentSpace IDE t))
+
+end RiemannianMetric
diff --git a/PhysLean/Mathematics/LinearAlgebra/BilinearForm.lean b/PhysLean/Mathematics/LinearAlgebra/BilinearForm.lean
new file mode 100644
index 000000000..06e7bc63b
--- /dev/null
+++ b/PhysLean/Mathematics/LinearAlgebra/BilinearForm.lean
@@ -0,0 +1,185 @@
+/-
+Copyright (c) 2025 Matteo Cipollina. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Matteo Cipollina
+-/
+
+import Mathlib.Analysis.InnerProductSpace.Basic
+import Mathlib.Analysis.Normed.Module.FiniteDimension
+import Mathlib.Analysis.Normed.Operator.BoundedLinearMaps
+import Mathlib.Topology.Algebra.Module.ModuleTopology
+import Mathlib.Topology.Metrizable.CompletelyMetrizable
+
+/-!
+# Bilinear Form Utilities
+
+This file contains general lemmas and definitions related to bilinear forms,
+particularly in the context of finite-dimensional normed spaces.
+-/
+
+namespace BilinearForm
+
+noncomputable section
+
+variable {E : Type*} [NormedAddCommGroup E] [NormedSpace ℝ E]
+
+/-- Given a bilinear map `B` and a linear map `A`, constructs the pullback of `B` by `A`,
+which is the bilinear map `(v, w) ↦ B (A v) (A w)`. -/
+abbrev pullback {F₁ F₂ R : Type*} [NormedAddCommGroup F₁] [NormedSpace ℝ F₁]
+    [NormedAddCommGroup F₂] [NormedSpace ℝ F₂] [NormedAddCommGroup R] [NormedSpace ℝ R]
+    (B : F₂ →L[ℝ] F₂ →L[ℝ] R) (A : F₁ →L[ℝ] F₂) : F₁ →L[ℝ] F₁ →L[ℝ] R :=
+  ContinuousLinearMap.bilinearComp B A A
+
+/-- For a finite-dimensional space E with basis b, constructs elementary bilinear forms e_ij
+    such that e_ij(v, w) = (b.coord i)(v) * (b.coord j)(w). -/
+lemma elementary_bilinear_forms_def [FiniteDimensional ℝ E]
+    (b : Basis (Fin (Module.finrank ℝ E)) ℝ E) :
+    ∃ (e : (Fin (Module.finrank ℝ E)) → (Fin (Module.finrank ℝ E)) → (E →L[ℝ] E →L[ℝ] ℝ)),
+    ∀ (i j : Fin (Module.finrank ℝ E)) (v w : E),
+    e i j v w = (b.coord i) v * (b.coord j) w := by
+  let n := Module.finrank ℝ E
+  let idx := Fin n
+  let b_dual (i : idx) := b.coord i
+  let e (i j : idx) : E →L[ℝ] E →L[ℝ] ℝ :=
+    let fi : E →L[ℝ] ℝ := LinearMap.toContinuousLinearMap (b_dual i)
+    let fj : E →L[ℝ] ℝ := LinearMap.toContinuousLinearMap (b_dual j)
+    let mul_map : ℝ →L[ℝ] ℝ →L[ℝ] ℝ := ContinuousLinearMap.mul ℝ ℝ
+    ContinuousLinearMap.bilinearComp mul_map fi fj
+  have h_prop : ∀ (i j : idx) (v w : E), e i j v w = (b.coord i) v * (b.coord j) w := by
+    intro i j v w
+    simp only [Basis.coord_apply]
+    rfl
+  exact ⟨e, h_prop⟩
+
+/-- Every vector in a finite-dimensional space can be written as a sum
+    of basis vectors with coordinates given by the dual basis. -/
+lemma vector_basis_expansion
+    (b : Basis (Fin (Module.finrank ℝ E)) ℝ E) (v : E) :
+    v = ∑ i ∈ Finset.univ, (b.coord i) v • b i := by
+  rw [← b.sum_repr v]
+  congr
+  ext i
+  simp only [b.coord_apply, Finsupp.single_apply]
+  rw [@Basis.repr_sum_self]
+
+lemma sum_smul_left (L : E →L[ℝ] E →L[ℝ] ℝ) {ι : Type*}
+    (s : Finset ι) (c : ι → ℝ) (v : ι → E) (w : E) :
+    L (∑ i ∈ s, c i • v i) w = ∑ i ∈ s, c i • L (v i) w := by
+  simp_rw [map_sum, ContinuousLinearMap.map_smul]
+  exact ContinuousLinearMap.sum_apply s (fun d => c d • L (v d)) w
+
+lemma sum_mul_left (L : E →L[ℝ] E →L[ℝ] ℝ) {ι : Type*}
+    (s : Finset ι) (c : ι → ℝ) (v : ι → E) (w : E) :
+    L (∑ i ∈ s, c i • v i) w = ∑ i ∈ s, c i * L (v i) w := by
+  rw [sum_smul_left L s c v w]
+  simp_rw [smul_eq_mul]
+
+lemma sum_smul_right (L : E →L[ℝ] E →L[ℝ] ℝ) (v : E) {ι : Type*}
+    (t : Finset ι) (d : ι → ℝ) (w : ι → E) :
+    L v (∑ j ∈ t, d j • w j) = ∑ j ∈ t, d j • L v (w j) := by
+  simp_rw [map_sum, ContinuousLinearMap.map_smul]
+
+lemma sum_mul_right (L : E →L[ℝ] E →L[ℝ] ℝ) (v : E) {ι : Type*}
+    (t : Finset ι) (d : ι → ℝ) (w : ι → E) :
+    L v (∑ j ∈ t, d j • w j) = ∑ j ∈ t, d j * L v (w j) := by
+  rw [sum_smul_right L v t d w]
+  simp_rw [smul_eq_mul]
+
+lemma sum_smul_left_right (L : E →L[ℝ] E →L[ℝ] ℝ) {ι₁ ι₂ : Type*}
+    (s : Finset ι₁) (t : Finset ι₂)
+    (c : ι₁ → ℝ) (v : ι₁ → E) (d : ι₂ → ℝ) (w : ι₂ → E) :
+    L (∑ i ∈ s, c i • v i) (∑ j ∈ t, d j • w j) =
+      ∑ i ∈ s, ∑ j ∈ t, c i • d j • L (v i) (w j) := by
+  simp_rw [sum_smul_left, sum_smul_right, Finset.smul_sum]
+
+lemma sum_mul_left_right (L : E →L[ℝ] E →L[ℝ] ℝ) {ι₁ ι₂ : Type*}
+    (s : Finset ι₁) (t : Finset ι₂)
+    (c : ι₁ → ℝ) (v : ι₁ → E) (d : ι₂ → ℝ) (w : ι₂ → E) :
+    L (∑ i ∈ s, c i • v i) (∑ j ∈ t, d j • w j) =
+      ∑ i ∈ s, ∑ j ∈ t, c i * d j * L (v i) (w j) := by
+  rw [sum_smul_left_right L s t c v d w]
+  apply Finset.sum_congr rfl
+  intro i _
+  apply Finset.sum_congr rfl
+  intro j _
+  simp_rw [smul_eq_mul, mul_assoc]
+
+lemma on_basis_expansions
+    (b : Basis (Fin (Module.finrank ℝ E)) ℝ E)
+    (L : E →L[ℝ] E →L[ℝ] ℝ) (v w : E) :
+    L v w = ∑ i ∈ Finset.univ, ∑ j ∈ Finset.univ,
+      (b.coord i v) * (b.coord j w) * L (b i) (b j) := by
+  rw [vector_basis_expansion b v, vector_basis_expansion b w]
+  rw [sum_mul_left_right L Finset.univ Finset.univ
+      (fun i => b.coord i v) b (fun j => b.coord j w) b]
+  simp only [← vector_basis_expansion b v, ← vector_basis_expansion b w]
+
+/-- The sum of elementary bilinear forms weighted by coefficients,
+    when applied to vectors, equals a weighted sum of the products of
+    coordinate functions. -/
+lemma elementary_bilinear_forms_sum_apply
+    (b : Basis (Fin (Module.finrank ℝ E)) ℝ E)
+    (e : (Fin (Module.finrank ℝ E)) → (Fin (Module.finrank ℝ E)) → (E →L[ℝ] E →L[ℝ] ℝ))
+    (h_e : ∀ i j v w, e i j v w = (b.coord i) v * (b.coord j) w)
+    (c : (Fin (Module.finrank ℝ E)) → (Fin (Module.finrank ℝ E)) → ℝ)
+    (v w : E) :
+    ((∑ i ∈ Finset.univ, ∑ j ∈ Finset.univ, c i j • e i j) v) w =
+    ∑ i ∈ Finset.univ, ∑ j ∈ Finset.univ, c i j * ((b.coord i) v * (b.coord j) w) := by
+  simp only [ContinuousLinearMap.sum_apply, ContinuousLinearMap.smul_apply]
+  apply Finset.sum_congr rfl
+  intro i _
+  apply Finset.sum_congr rfl
+  intro j _
+  congr
+  exact h_e i j v w
+
+/-- A simple algebraic identity for rearranging terms in double sums with multiple scalars. -/
+lemma Finset.rearrange_double_sum_coefficients {α β R : Type*} [CommSemiring R]
+    {s : Finset α} {t : Finset β} {f : α → β → R} {g : α → β → R} {h : α → β → R} :
+    (∑ i ∈ s, ∑ j ∈ t, f i j * g i j * h i j) =
+    (∑ i ∈ s, ∑ j ∈ t, h i j * (f i j * g i j)) := by
+  apply Finset.sum_congr rfl; intro i _
+  apply Finset.sum_congr rfl; intro j _
+  ring
+
+/-- Any bilinear form can be decomposed as a sum of elementary bilinear forms
+    weighted by the form's values on basis elements. -/
+theorem decomposition [FiniteDimensional ℝ E] (b : Basis (Fin (Module.finrank ℝ E)) ℝ E) :
+    ∀ (L : E →L[ℝ] E →L[ℝ] ℝ),
+    ∃ (e : (Fin (Module.finrank ℝ E)) → (Fin (Module.finrank ℝ E)) → (E →L[ℝ] E →L[ℝ] ℝ)),
+    (∀ i j v w, e i j v w = (b.coord i) v * (b.coord j) w) ∧
+    L = ∑ i ∈ Finset.univ, ∑ j ∈ Finset.univ, L (b i) (b j) • e i j := by
+  intro L
+  obtain ⟨e, h_e⟩ := elementary_bilinear_forms_def b
+  have h_decomp : L = ∑ i ∈ Finset.univ, ∑ j ∈ Finset.univ, L (b i) (b j) • e i j := by
+    ext v w
+    have expansion := on_basis_expansions b L v w
+    have right_side := elementary_bilinear_forms_sum_apply b e h_e (λ i j => L (b i) (b j)) v w
+    rw [expansion, Finset.rearrange_double_sum_coefficients]
+    exact right_side.symm
+  exact ⟨e, ⟨h_e, h_decomp⟩⟩
+
+/-- Given two vectors `v` and `w`, `eval_vectors_continuousLinear v w` is the continuous linear map
+that evaluates a bilinear form `B` at `(v, w)`. -/
+def eval_vectors_continuousLinear (v w : E) :
+    (E →L[ℝ] E →L[ℝ] ℝ) →L[ℝ] ℝ :=
+  { toFun := fun B => B v w,
+    map_add' := fun f g => by simp [ContinuousLinearMap.add_apply],
+    map_smul' := fun c f => by simp [ContinuousLinearMap.smul_apply] }
+
+namespace ContinuousLinearMap
+
+lemma map_sum_clm {R S M M₂ : Type*}
+    [Semiring R] [Semiring S] [AddCommMonoid M] [AddCommMonoid M₂]
+    {σ : R →+* S} [Module R M] [Module S M₂] (f : M →ₛₗ[σ] M₂) {ι : Type*} (t : Finset ι)
+    (g : ι → M) : f (∑ i ∈ t, g i) = ∑ i ∈ t, f (g i) :=
+  _root_.map_sum f g t
+
+lemma map_finset_sum {𝕜 : Type*} [NontriviallyNormedField 𝕜]
+    {E F : Type*} [NormedAddCommGroup E] [NormedSpace 𝕜 E]
+    [NormedAddCommGroup F] [NormedSpace 𝕜 F]
+    (f : E →L[𝕜] F) {ι : Type*} (s : Finset ι) (g : ι → E) :
+    f (∑ i ∈ s, g i) = ∑ i ∈ s, f (g i) :=
+  _root_.map_sum f g s
+
+end ContinuousLinearMap