diff --git a/PhysLean.lean b/PhysLean.lean
index 9a8b6ca79..a67007d9d 100644
--- a/PhysLean.lean
+++ b/PhysLean.lean
@@ -36,6 +36,7 @@ import PhysLean.Electromagnetism.LorentzAction
 import PhysLean.Electromagnetism.MaxwellEquations
 import PhysLean.Electromagnetism.Potential
 import PhysLean.Electromagnetism.Wave
+import PhysLean.Mathematics.Analysis.ContDiff
 import PhysLean.Mathematics.Calculus.AdjFDeriv
 import PhysLean.Mathematics.Calculus.Divergence
 import PhysLean.Mathematics.DataStructures.FourTree.Basic
@@ -49,8 +50,15 @@ import PhysLean.Mathematics.Distribution.PowMul
 import PhysLean.Mathematics.FDerivCurry
 import PhysLean.Mathematics.Fin
 import PhysLean.Mathematics.Fin.Involutions
+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.Manifold.PartialHomeomorph
+import PhysLean.Mathematics.Geometry.Metric.PseudoRiemannian.Chart
 import PhysLean.Mathematics.Geometry.Metric.PseudoRiemannian.Defs
 import PhysLean.Mathematics.Geometry.Metric.Riemannian.Defs
+import PhysLean.Mathematics.LinearAlgebra.BilinearForm
 import PhysLean.Mathematics.InnerProductSpace.Adjoint
 import PhysLean.Mathematics.InnerProductSpace.Basic
 import PhysLean.Mathematics.InnerProductSpace.Calculus
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/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