/-- The affine coordinate has the normalization required at both annular endpoints. -/
theorem annulusParameter_mul (q : ℝ) (hq : q≠1) : annulusParameter ⟨q, hq⟩*(1-q)=1 := by
  unfold annulusParameter
  exact inv_mul_cancel₀ (sub_ne_zero.mpr (Ne.symm hq))

/-- A nondegenerate annulus puts the third singularity strictly beyond one. -/
theorem one_lt_annulusParameter {q : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1) :
    1 < annulusParameter ⟨q, ne_of_lt hq.2⟩ := by
  unfold annulusParameter
  rw [one_lt_inv₀ (by linarith [hq.2])]
  linarith [hq.1]

/-- The entire unit parameter interval has a nonnegative squared radius. -/
theorem zernikeAnnulusCoordinates_radical_nonneg
    (q : {q : ℝ // 0 ≤ q ∧ q ≠ 1}) (x : Icc (0 : ℝ) 1) :
    0 ≤ 1-(x : ℝ)/annulusParameter ⟨q,q.property.2⟩ := by
  have he : 1-(x : ℝ)/annulusParameter ⟨q,q.property.2⟩ =
      (1-(x : ℝ))+(q : ℝ)*(x : ℝ) := by
    simp only [annulusParameter, div_inv_eq_mul]
    ring
  rw [he]
  exact add_nonneg (sub_nonneg.mpr x.property.2) (mul_nonneg q.property.1 x.property.1)

/-- All polynomial moments exist for the auxiliary product measure. -/
@[aesop (rule_sets := [PolynomialIntegrability]) unsafe 80% apply]
theorem hasFiniteMoments_zernikeAnnulusProductMeasure {q a b : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1)
    (ha : -1<a) (hb : -1<b) : HasFiniteMoments (zernikeAnnulusProductMeasure ⟨q,hq.1.le,ne_of_lt hq.2⟩ a b) :=
  fun p ↦ integrable_zernikeAnnulusProductMeasure hq ha hb (p.continuous_eval.pow 2)

/-- The geometric scalar agrees with the Cartesian radial Laplacian. -/
theorem semiclassicalJacobiDirichletMap_eq (t γ : ℝ) (ht : t≠0) (p : ℝ[X]) :
    semiclassicalJacobiDirichletMap ⟨t,ht⟩ γ p =
      semiclassicalJacobiAnnularLaplacianMap t γ (annularDirichletWeight ⟨t,ht⟩*p) := by
  simp only [semiclassicalJacobiDirichletMap, semiclassicalJacobiAnnularLaplacianMap,
    LinearMap.smul_apply, LinearMap.comp_apply]
  have he : annularDirichletWeight ⟨t,ht⟩*p = t⁻¹^2 • (X*(1-X)*p) := by
    simp [annularDirichletWeight, smul_eq_C_mul, mul_assoc]
  rw [he, map_smul, map_smul, smul_smul]
  congr 1
  field_simp

/-- Squared Dirichlet energy in one radial sector. -/
noncomputable def semiclassicalJacobiDirichletEnergy (t : Set.Ioi (1 : ℝ)) (γ : ℝ)
    (p : ℝ[X]) : ℝ :=
  -integralInner (semiclassicalJacobiMeasure t 1 1 γ)
    (IntegrablePair.of (semiclassicalJacobiMeasure t 1 1 γ)
      (semiclassicalJacobiDirichletMap ⟨t,ne_of_gt (lt_trans zero_lt_one t.property)⟩ γ p) p
      (h := (hasFiniteMoments_semiclassicalJacobiMeasure_one_one
        (γ := γ) t.property).integrable_mul _ _))

/-- Energy orthogonality gives the exact Pythagorean error identity. -/
theorem semiclassicalJacobiDirichletEnergy_pythagoras {t : ℝ} (ht : 1<t) (γ : ℝ)
    (p q : ℝ[X])
    (h : integralInner (semiclassicalJacobiMeasure ⟨t, by assumption⟩ 1 1 γ) (IntegrablePair.of (semiclassicalJacobiMeasure ⟨t, by assumption⟩ 1 1 γ) (semiclassicalJacobiDirichletMap ⟨t,ne_of_gt (lt_trans zero_lt_one ht)⟩ γ p) q (h := (hasFiniteMoments_semiclassicalJacobiMeasure_one_one (γ := γ) ht).integrable_mul _ _))=0) :
    semiclassicalJacobiDirichletEnergy ⟨t, ht⟩ γ (p+q) =
      semiclassicalJacobiDirichletEnergy ⟨t, ht⟩ γ p+semiclassicalJacobiDirichletEnergy ⟨t, ht⟩ γ q := by
  have hμ := hasFiniteMoments_semiclassicalJacobiMeasure (γ := γ) ht
    (by norm_num : (-1 : ℝ)<1) (by norm_num : (-1 : ℝ)<1)
  have hr (u v w : ℝ[X]) : integralInner (semiclassicalJacobiMeasure ⟨t, by assumption⟩ 1 1 γ) (IntegrablePair.of (semiclassicalJacobiMeasure ⟨t, by assumption⟩ 1 1 γ) u (v+w) (h := (hasFiniteMoments_semiclassicalJacobiMeasure_one_one (γ := γ) ht).integrable_mul _ _)) =
      integralInner (semiclassicalJacobiMeasure ⟨t, by assumption⟩ 1 1 γ) (IntegrablePair.of (semiclassicalJacobiMeasure ⟨t, by assumption⟩ 1 1 γ) u v (h := (hasFiniteMoments_semiclassicalJacobiMeasure_one_one (γ := γ) ht).integrable_mul _ _))+
        integralInner (semiclassicalJacobiMeasure ⟨t, by assumption⟩ 1 1 γ) (IntegrablePair.of (semiclassicalJacobiMeasure ⟨t, by assumption⟩ 1 1 γ) u w (h := (hasFiniteMoments_semiclassicalJacobiMeasure_one_one (γ := γ) ht).integrable_mul _ _)) := by
    rw [integralInner_comm, integralInner_add_left hμ]
    rw [integralInner_comm _ v u, integralInner_comm _ w u]
  unfold semiclassicalJacobiDirichletEnergy
  simp only [map_add]
  rw [integralInner_add_left hμ, hr, hr,
    integralInner_semiclassicalJacobi_dirichlet_symm ht γ q p, h]
  ring

/-- Galerkin orthogonality makes the discrete polynomial solution energy-optimal. -/
theorem semiclassicalJacobiDirichletGalerkin_bestApproximation {t : ℝ} (ht : 1<t)
    (γ : ℝ) (N : ℕ) (u : ℝ[X]) (v : Fin N → ℝ)
    (hv : semiclassicalJacobiDirichletSection ht γ N *ᵥ v =
      semiclassicalJacobiDirichletLoad ht γ N (semiclassicalJacobiDirichletMap ⟨t,ne_of_gt (lt_trans zero_lt_one ht)⟩ γ u))
    (w : Fin N → ℝ) :
    semiclassicalJacobiDirichletEnergy ⟨t, ht⟩ γ (u-semiclassicalJacobiDirichletTrial ht γ N v) ≤
      semiclassicalJacobiDirichletEnergy ⟨t, ht⟩ γ (u-semiclassicalJacobiDirichletTrial ht γ N w) := by
  let vP := semiclassicalJacobiDirichletTrial ht γ N v
  let wP := semiclassicalJacobiDirichletTrial ht γ N w
  have h := (semiclassicalJacobiDirichletSection_eq_iff_galerkin ht γ N _ v).mp hv (v-w)
  have he : semiclassicalJacobiDirichletTrial ht γ N (v-w)=vP-wP := by
    simp [semiclassicalJacobiDirichletTrial, vP, wP, sub_smul, Finset.sum_sub_distrib]
  simp only [he, ← map_sub] at h
  have ho : integralInner (semiclassicalJacobiMeasure ⟨t, by assumption⟩ 1 1 γ) (IntegrablePair.of (semiclassicalJacobiMeasure ⟨t, by assumption⟩ 1 1 γ) (semiclassicalJacobiDirichletMap ⟨t,ne_of_gt (lt_trans zero_lt_one ht)⟩ γ (u-vP)) (vP-wP) (h := (hasFiniteMoments_semiclassicalJacobiMeasure_one_one (γ := γ) ht).integrable_mul _ _))=0 := by
    have hm : u-vP=-(vP-u) := by abel
    simp only [hm, map_neg, ← neg_one_smul ℝ
      (semiclassicalJacobiDirichletMap ⟨t,ne_of_gt (lt_trans zero_lt_one ht)⟩ γ (vP-u))]
    rw [integralInner_smul_left (h := (hasFiniteMoments_semiclassicalJacobiMeasure_one_one (γ := γ) ht).integrable_mul _ _), h, mul_zero]
  have hp := semiclassicalJacobiDirichletEnergy_pythagoras ht γ (u-vP) (vP-wP) ho
  have hsum : (u-vP)+(vP-wP)=u-wP := by abel
  simp only [hsum] at hp
  have hn := semiclassicalJacobiDirichletEnergy_nonneg ht γ (vP-wP)
  change semiclassicalJacobiDirichletEnergy ⟨t, ht⟩ γ (u-vP) ≤
    semiclassicalJacobiDirichletEnergy ⟨t, ht⟩ γ (u-wP)
  linarith

/-- The completed inner product extends the polynomial Dirichlet form exactly. -/
theorem inner_embed (t : Set.Ioi (1 : ℝ)) (γ : ℝ) (p q : PolynomialSpace t γ) :
    inner ℝ (embed t γ p) (embed t γ q) =
      -integralInner (semiclassicalJacobiMeasure t 1 1 γ) (IntegrablePair.of (semiclassicalJacobiMeasure t 1 1 γ) (semiclassicalJacobiDirichletMap ⟨t,ne_of_gt (lt_trans zero_lt_one t.property)⟩ γ (polynomialEquiv t γ p)) (polynomialEquiv t γ q) (h := (hasFiniteMoments_semiclassicalJacobiMeasure_one_one (γ := γ) (gt_one t)).integrable_mul _ _)) := by
  change inner ℝ (p : Space t γ) (q : Space t γ) = _
  rw [UniformSpace.Completion.inner_coe]
  rfl

/-- Radial polynomials are dense in the chosen energy completion. -/
theorem denseRange_embed (t : Set.Ioi (1 : ℝ)) (γ : ℝ) : DenseRange (embed t γ) :=
  UniformSpace.Completion.denseRange_coe

/-- Distinct angular sectors are orthogonal for the planar negative-Laplacian pairing. -/
theorem integralInner_annulusDirichletLaplacian_eq_zero {q : ℝ}
    (hq : q ∈ Ioo (0 : ℝ) 1) {u v : Bool} {m n : ℕ}
    (h : u ≠ v ∨ m ≠ n) (p z : Polynomial ℝ) :
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (diskLaplacian (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ u m p)) (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ v n z) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) = 0 := by
  dsimp only [integralInner, IntegrablePair.of]
  rw [diskLaplacian_annulusDirichletLift hq, ← annulusLift_dirichletWeight hq]
  exact integralInner_annulusRadialLift_eq_zero_of_angular (a := 0) (b := 0) hq (by norm_num) (by norm_num) _ _ h

/-- The radial energy completion agrees with the polynomial planar pairing after scaling. -/
theorem neg_integralInner_annulusDirichletLaplacian {q : ℝ}
    (hq : q ∈ Ioo (0 : ℝ) 1) (u : Bool) (m : ℕ) (p z : Polynomial ℝ) :
    -integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (diskLaplacian (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ u m p)) (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ u m z) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) =
      annulusEnergyScale ⟨q,ne_of_lt hq.2⟩ u m * inner ℝ
        (SemiclassicalJacobiEnergy.polynomialEmbed ⟨annulusParameter ⟨q,ne_of_lt hq.2⟩, one_lt_annulusParameter hq⟩ m p)
        (SemiclassicalJacobiEnergy.polynomialEmbed ⟨annulusParameter ⟨q,ne_of_lt hq.2⟩, one_lt_annulusParameter hq⟩ m z) := by
  simp only [diskLaplacian_annulusDirichletLift hq]
  rw [integralInner_annulusLift_dirichletLift hq,
    SemiclassicalJacobiEnergy.inner_polynomialEmbed]
  ring

/-- Different valid angular coordinates are orthogonal. -/
theorem inner_polynomialEmbed_eq_zero (q : Set.Ioo (0 : ℝ) 1) {s r : Sector}
    (h : s ≠ r) (p z : Polynomial ℝ) :
    inner ℝ (polynomialEmbed q s p) (polynomialEmbed q r z) = 0 := by
  simp [polynomialEmbed, sectorEmbed, real_inner_smul_left, inner_smul_right,
    lp.inner_single_left, lp.single_apply, h]

/-- The combined angular and radial polynomial trial spaces have dense union. -/
theorem trialSpace_closure_iSup (q : Set.Ioo (0 : ℝ) 1) :
    (⨆ c, trialSpace q c).topologicalClosure = ⊤ :=
  lp.HilbertSumGalerkin.trialSpace_closure_iSup _
    (fun s ↦ SemiclassicalJacobiEnergy.trialSpace_mono (parameter q) s.val.2)
    (fun s ↦ SemiclassicalJacobiEnergy.trialSpace_closure_iSup (parameter q) s.val.2)

/-- The Cartesian Laplacian product rule, with its cross-gradient term explicit. -/
theorem diskLaplacian_mul (p z : MvPolynomial (Fin 2) ℝ) :
    diskLaplacian (p*z) = diskLaplacian p*z + p*diskLaplacian z + 2*diskGradientProduct p z := by
  simp only [diskLaplacian, diskGradientProduct, pderiv_mul, map_add]
  ring

/-- The integral of a radial-sector Laplacian consists only of the constant-mode boundary flux. -/
theorem integral_diskLaplacian_diskRadialLift {q : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1)
    (u : Bool) (m : ℕ) (p : Polynomial ℝ) :
    ∫ v, eval v (diskLaplacian (diskRadialLift u p m)) ∂zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0 =
      if u = false ∧ m = 0 then
        4*Real.pi*(p.derivative.eval 1-q*p.derivative.eval q) else 0 := by
  have hone : diskRadialLift false (1 : Polynomial ℝ) 0 = 1 := by
    simp [diskRadialLift, diskHarmonics]
  have he (f : MvPolynomial (Fin 2) ℝ) :
      (∫ v, eval v f ∂zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) =
        integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) f (diskRadialLift false 1 0) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) := by
    simp [hone, integralInner]
  rw [he]
  simp only [diskLaplacian_diskRadialLift]
  rw [integralInner_annulusRadialLift hq (by norm_num) (by norm_num)]
  by_cases h : u = false ∧ m = 0
  · obtain ⟨rfl, rfl⟩ := h
    simp only [and_self, ↓reduceIte, Real.rpow_zero, Nat.cast_zero, zero_add, pow_zero,
      Polynomial.eval_one, mul_one, one_mul, zernikeAngularNormSq, Bool.false_eq_true]
    have hd := intervalIntegral.integral_eq_sub_of_hasDerivAt
      (a := q) (b := 1) (fun x _ ↦ (Polynomial.X*p.derivative).hasDerivAt x)
      ((Polynomial.X*p.derivative).derivative.continuous.intervalIntegrable q 1)
    have hi : (∫ s in q..1,
        (4*(Polynomial.X*p.derivative.derivative+p.derivative)).eval s) =
        4*(p.derivative.eval 1-q*p.derivative.eval q) := by
      calc
        _ = (4 : ℝ) * ∫ s in q..1, (Polynomial.X*p.derivative).derivative.eval s := by
          rw [← intervalIntegral.integral_const_mul]
          apply intervalIntegral.integral_congr
          intro s hs
          simp [Polynomial.derivative_mul]
          ring
        _ = _ := by rw [hd]; simp
    rw [intervalIntegral.integral_of_le hq.2.le, ← restrict_Ioo_eq_restrict_Ioc,
      ← integral_subtype_comap measurableSet_Ioo] at hi
    have hh := congrArg (fun a : ℝ ↦ (1/2)*a*(2*Real.pi)) hi
    convert hh using 1
    ring
  · simp [h]

/-- Any Cartesian polynomial with a squared boundary factor has zero integrated Laplacian. -/
theorem integral_diskLaplacian_boundary_sq {q : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1)
    (p : MvPolynomial (Fin 2) ℝ) :
    ∫ v, eval v (diskLaplacian (annulusBoundaryWeight q ^ 2 * p))
      ∂zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0 = 0 := by
  let μ := zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0
  let hμ := hasFiniteMoments_zernikeAnnulusMeasure hq
    (by norm_num : (-1 : ℝ) < 0) (by norm_num : (-1 : ℝ) < 0)
  let L : MvPolynomial (Fin 2) ℝ →ₗ[ℝ] ℝ :=
    ((integralBilinForm μ hμ).flip 1).comp
      (diskLaplacianLinearMap.comp (LinearMap.mulLeft ℝ (annulusBoundaryWeight q ^ 2)))
  have he (z : MvPolynomial (Fin 2) ℝ) : L z =
      ∫ v, eval v (diskLaplacian (annulusBoundaryWeight q ^ 2*z)) ∂μ := by
    simp [L, integralBilinForm, integralInner, diskLaplacianLinearMap]
  have hk : LinearMap.ker L = ⊤ := eq_top_of_annulusLift_mem hq _ (by
    intro u m z
    rw [LinearMap.mem_ker, he]
    exact integral_diskLaplacian_boundary_sq_annulusLift hq u m z)
  have hp : p ∈ LinearMap.ker L := by rw [hk]; trivial
  rw [LinearMap.mem_ker, he] at hp
  exact hp

/-- The polynomial negative-Laplacian pairing is symmetric on the entire boundary-weighted space. -/
theorem integralInner_annulus_boundary_laplacian_symm {q : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1)
    (p z : MvPolynomial (Fin 2) ℝ) :
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (diskLaplacian (annulusBoundaryWeight q*p)) (annulusBoundaryWeight q*z) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) =
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (diskLaplacian (annulusBoundaryWeight q*z)) (annulusBoundaryWeight q*p) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) := by
  let μ := zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0
  let hμ := hasFiniteMoments_zernikeAnnulusMeasure hq
    (by norm_num : (-1 : ℝ) < 0) (by norm_num : (-1 : ℝ) < 0)
  let B := LinearMap.mulLeft ℝ (annulusBoundaryWeight q)
  let F := (integralBilinForm μ hμ).compl₁₂ (diskLaplacianLinearMap.comp B) B
  have hm (u v : Bool) (m n : ℕ) (a b : Polynomial ℝ) :
      (F-F.flip) (annulusLift ⟨q,ne_of_lt hq.2⟩ u m a) (annulusLift ⟨q,ne_of_lt hq.2⟩ v n b) = 0 := by
    change integralInner μ (IntegrablePair.of μ (diskLaplacian (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ u m a)) (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ v n b) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) -
      integralInner μ (IntegrablePair.of μ (diskLaplacian (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ v n b)) (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ u m a) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) = 0
    by_cases h : u=v ∧ m=n
    · obtain ⟨rfl, rfl⟩ := h
      dsimp [μ]
      simp only [diskLaplacian_annulusDirichletLift hq]
      rw [integralInner_annulusLift_dirichletLift hq, integralInner_annulusLift_dirichletLift hq,
        integralInner_semiclassicalJacobi_dirichlet_symm (one_lt_annulusParameter hq)]
      exact sub_self _
    · have hne : u≠v ∨ m≠n := by tauto
      rw [integralInner_annulusDirichletLaplacian_eq_zero hq hne,
        integralInner_annulusDirichletLaplacian_eq_zero hq (by tauto)]
      exact sub_self _
  have hk : LinearMap.ker (F-F.flip) = ⊤ := eq_top_of_annulusLift_mem hq _ (by
    intro u m a
    rw [LinearMap.mem_ker]
    have hs : LinearMap.ker ((F-F.flip) (annulusLift ⟨q,ne_of_lt hq.2⟩ u m a)) = ⊤ :=
      eq_top_of_annulusLift_mem hq _ (by intro v n b; exact hm u v m n a b)
    apply LinearMap.ext
    intro y
    have hy : y ∈ LinearMap.ker ((F-F.flip) (annulusLift ⟨q,ne_of_lt hq.2⟩ u m a)) := by rw [hs]; trivial
    exact hy)
  have hp : p ∈ LinearMap.ker (F-F.flip) := by rw [hk]; trivial
  have he := LinearMap.congr_fun (LinearMap.mem_ker.mp hp) z
  exact sub_eq_zero.mp he

/-- Green's first identity for all polynomial trials vanishing by the annular boundary factor. -/
theorem integralInner_annulus_gradient_eq_neg_laplacian {q : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1)
    (p z : MvPolynomial (Fin 2) ℝ) :
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (pderiv 0 (annulusBoundaryWeight q*p)) (pderiv 0 (annulusBoundaryWeight q*z)) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) +
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (pderiv 1 (annulusBoundaryWeight q*p)) (pderiv 1 (annulusBoundaryWeight q*z)) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) =
    -integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (diskLaplacian (annulusBoundaryWeight q*p)) (annulusBoundaryWeight q*z) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) := by
  let μ := zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0
  let hμ := hasFiniteMoments_zernikeAnnulusMeasure hq
    (by norm_num : (-1 : ℝ) < 0) (by norm_num : (-1 : ℝ) < 0)
  let J : MvPolynomial (Fin 2) ℝ →ₗ[ℝ] ℝ := (integralBilinForm μ hμ).flip 1
  have hJ (f g : MvPolynomial (Fin 2) ℝ) : J (f*g) = integralInner μ (IntegrablePair.of μ f g (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) := by
    simp [J, integralBilinForm, integralInner]
  let f := annulusBoundaryWeight q*p
  let g := annulusBoundaryWeight q*z
  have hz : J (diskLaplacian (f*g)) = 0 := by
    have he : f*g = annulusBoundaryWeight q^2*(p*z) := by dsimp [f,g]; ring
    rw [he]
    simpa [J, integralBilinForm, integralInner] using integral_diskLaplacian_boundary_sq hq (p*z)
  rw [diskLaplacian_mul, map_add, map_add, hJ, hJ] at hz
  have ht : J (2*diskGradientProduct f g) = 2*J (diskGradientProduct f g) := by
    rw [two_mul, map_add, two_mul]
  rw [ht, diskGradientProduct, map_add, hJ, hJ, integralInner_comm μ f,
    integralInner_annulus_boundary_laplacian_symm hq z p] at hz
  change integralInner μ (IntegrablePair.of μ (pderiv 0 f) (pderiv 0 g) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) +
    integralInner μ (IntegrablePair.of μ (pderiv 1 f) (pderiv 1 g) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) = -integralInner μ (IntegrablePair.of μ (diskLaplacian f) g (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _))
  linarith

/-- Distinct angular sectors are orthogonal for Cartesian gradient energy. -/
theorem integralInner_annulusDirichlet_gradient_eq_zero {q : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1)
    {u v : Bool} {m n : ℕ} (h : u ≠ v ∨ m ≠ n) (p z : Polynomial ℝ) :
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (pderiv 0 (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ u m p)) (pderiv 0 (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ v n z)) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) +
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (pderiv 1 (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ u m p)) (pderiv 1 (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ v n z)) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) = 0 := by
  calc
    _ = -integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (diskLaplacian (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ u m p)) (annulusDirichletLift ⟨q,ne_of_lt hq.2⟩ v n z) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) :=
      integralInner_annulus_gradient_eq_neg_laplacian hq _ _
    _ = 0 := by rw [integralInner_annulusDirichletLaplacian_eq_zero hq h p z, neg_zero]

/-- Weighted L² space of endpoint-weighted radial derivatives. -/
abbrev DerivativeSpace (t : Set.Ioi (1 : ℝ)) (γ : ℝ) :=
  Lp ℝ 2 (semiclassicalJacobiMeasure t 0 0 (γ + 1))

/-- The scaled endpoint-weighted derivative as an L² class. -/
def polynomialDerivative (t : Set.Ioi (1 : ℝ)) (γ : ℝ) :
    PolynomialSpace t γ →ₗ[ℝ] DerivativeSpace t γ :=
  Real.sqrt (4 / (t : ℝ)) •
    (polynomialToL2 (derivativeMeasure_moments t γ)).comp
      (semiclassicalJacobiABWeightedDerivativeMap.comp (polynomialEquiv t γ).toLinearMap)

/-- The completed energy space embeds isometrically into concrete weighted L² derivatives. -/
def derivativeRealization (t : Set.Ioi (1 : ℝ)) (γ : ℝ) :
    Space t γ →ₗᵢ[ℝ] DerivativeSpace t γ := (polynomialDerivativeIsometry t γ).fromCompletion

/-- The completed realization has the prescribed value on polynomial vectors. -/
theorem derivativeRealization_embed (t : Set.Ioi (1 : ℝ)) (γ : ℝ)
    (p : PolynomialSpace t γ) :
    derivativeRealization t γ (embed t γ p) = polynomialDerivative t γ p :=
  (polynomialDerivativeIsometry t γ).fromCompletion_apply_coe p

/-- The realization is faithful: no two completed energy vectors have the same derivative class. -/
theorem derivativeRealization_injective (t : Set.Ioi (1 : ℝ)) (γ : ℝ) :
    Function.Injective (derivativeRealization t γ) := (derivativeRealization t γ).injective

/-- The realized space is exactly the L² closure of the scaled polynomial derivatives. -/
theorem range_derivativeRealization (t : Set.Ioi (1 : ℝ)) (γ : ℝ) :
    Set.range (derivativeRealization t γ) = closure (Set.range (polynomialDerivative t γ)) := by
  have he : derivativeRealization t γ '' Set.range (embed t γ) =
      Set.range (polynomialDerivative t γ) := by
    rw [← Set.range_comp]
    congr 1
    funext p
    exact derivativeRealization_embed t γ p
  apply Set.Subset.antisymm
  · have h := (derivativeRealization t γ).continuous.range_subset_closure_image_dense
      (denseRange_embed t γ)
    rwa [he] at h
  · apply closure_minimal _ (isClosed_range_derivativeRealization t γ)
    rintro y ⟨p, rfl⟩
    exact ⟨embed t γ p, derivativeRealization_embed t γ p⟩

/-- The realization preserves the completed energy inner product. -/
theorem inner_derivativeRealization (t : Set.Ioi (1 : ℝ)) (γ : ℝ) (x y : Space t γ) :
    inner ℝ (derivativeRealization t γ x) (derivativeRealization t γ y) = inner ℝ x y :=
  (derivativeRealization t γ).inner_map_map x y

/-- The derivative-data load is the concrete weighted L² pairing. -/
theorem derivativeLoad_apply (t : Set.Ioi (1 : ℝ)) (γ : ℝ)
    (g : DerivativeSpace t γ) (x : Space t γ) :
    derivativeLoad t γ g x = inner ℝ g (derivativeRealization t γ x) := rfl

/-- The derivative-data load is bounded by the L² norm of its data. -/
theorem norm_derivativeLoad_le (t : Set.Ioi (1 : ℝ)) (γ : ℝ) (g : DerivativeSpace t γ) :
    ‖derivativeLoad t γ g‖ ≤ ‖g‖ := by
  apply ContinuousLinearMap.opNorm_le_bound _ (norm_nonneg g)
  intro x
  rw [derivativeLoad_apply]
  exact (norm_inner_le_norm g _).trans_eq (by rw [(derivativeRealization t γ).norm_map])

/-- The finite solves are uniformly bounded for every square-integrable derivative datum. -/
theorem norm_coefficients_derivativeLoad_le (t : Set.Ioi (1 : ℝ)) (γ : ℝ) (N : ℕ)
    (g : DerivativeSpace t γ) :
    ‖trialEmbed t γ N (coefficients t γ N (derivativeLoad t γ g))‖ ≤ ‖g‖ :=
  (norm_trialEmbed_coefficients_le t γ N _).trans (norm_derivativeLoad_le t γ g)

/-- The finite coefficient solves converge as concrete weighted L² derivative classes. -/
theorem tendsto_derivativeRealization_coefficients (t : Set.Ioi (1 : ℝ)) (γ : ℝ)
    (ℓ : Space t γ →L[ℝ] ℝ) :
    Tendsto (fun N ↦ derivativeRealization t γ
      (trialEmbed t γ N (coefficients t γ N ℓ))) atTop
        (𝓝 (derivativeRealization t γ (solution t γ ℓ))) :=
  (derivativeRealization t γ).continuous.tendsto _ |>.comp (tendsto_trialEmbed_coefficients t γ ℓ)

/-- Weighted L² derivative classes in all valid real angular sectors. -/
abbrev DerivativeSpace (q : Set.Ioo (0 : ℝ) 1) :=
  lp (fun s : Sector ↦ SemiclassicalJacobiEnergy.DerivativeSpace (parameter q) s.val.2) 2

/-- Realize all energy coordinates simultaneously as concrete weighted L² derivatives. -/
def derivativeRealization (q : Set.Ioo (0 : ℝ) 1) : Space q →ₗᵢ[ℝ] DerivativeSpace q where
  toFun x := ⟨fun s ↦ SemiclassicalJacobiEnergy.derivativeRealization (parameter q) s.val.2 (x s),
    memℓp_norm_iff.mp (by simpa only [LinearIsometry.norm_map] using (lp.memℓp x).norm)⟩
  map_add' x y := by
    apply lp.ext
    funext s
    exact map_add (SemiclassicalJacobiEnergy.derivativeRealization (parameter q) s.val.2) _ _
  map_smul' a x := by
    apply lp.ext
    funext s
    exact map_smul (SemiclassicalJacobiEnergy.derivativeRealization (parameter q) s.val.2) a _
  norm_map' x := by
    rw [lp.norm_eq_tsum_rpow (by norm_num), lp.norm_eq_tsum_rpow (by norm_num)]
    congr 1
    apply tsum_congr
    intro s
    change ‖SemiclassicalJacobiEnergy.derivativeRealization (parameter q) s.val.2 (x s)‖ ^
      (2 : ENNReal).toReal = ‖x s‖ ^ (2 : ENNReal).toReal
    rw [LinearIsometry.norm_map]

/-- Completed annular energy has a faithful realization in the Hilbert sum of derivative classes. -/
theorem derivativeRealization_injective (q : Set.Ioo (0 : ℝ) 1) :
    Function.Injective (derivativeRealization q) := (derivativeRealization q).injective

/-- The realized annular derivative space is closed. -/
theorem isClosed_range_derivativeRealization (q : Set.Ioo (0 : ℝ) 1) :
    IsClosed (Set.range (derivativeRealization q)) :=
  (derivativeRealization q).isometry.isClosedEmbedding.isClosed_range

/-- The finite solver's energy error equals its concrete derivative-space error. -/
theorem norm_sub_derivativeRealization_finiteSolve (q : Set.Ioo (0 : ℝ) 1)
    (c : Finset Sector × ℕ) (ℓ : Space q →L[ℝ] ℝ) :
    ‖derivativeRealization q (solution q ℓ) - derivativeRealization q (finiteSolve q c ℓ)‖ =
      ‖solution q ℓ - finiteSolve q c ℓ‖ := by
  rw [← map_sub, LinearIsometry.norm_map]

/-- Joint angular and radial convergence also holds in the concrete derivative Hilbert sum. -/
theorem tendsto_derivativeRealization_finiteSolve (q : Set.Ioo (0 : ℝ) 1)
    (ℓ : Space q →L[ℝ] ℝ) :
    Tendsto (fun N ↦ derivativeRealization q (finiteSolve q (frequencyCutoff N, N) ℓ))
      atTop (𝓝 (derivativeRealization q (solution q ℓ))) :=
  (derivativeRealization q).continuous.tendsto _ |>.comp
    (tendsto_finiteSolve_frequencyCutoff q ℓ)

/-- Physical L² classes on the annulus, with its unweighted area measure. -/
abbrev PhysicalL2 (q : Set.Ioo (0 : ℝ) 1) := Lp ℝ 2 (zernikeAnnulusMeasure ⟨q,q.property.1.le,q.property.2⟩ 0 0)

/-- The physical L² realization of one completed radial sector. -/
def radialL2 (q : Set.Ioo (0 : ℝ) 1) (s : Sector) : RadialSpace q s →L[ℝ] PhysicalL2 q :=
  (radialPolynomialL2 q s).extendOfNorm
    (SemiclassicalJacobiEnergy.polynomialEmbed (parameter q) s.val.2)

/-- The extension retains the prescribed boundary-weighted polynomial values. -/
theorem radialL2_polynomialEmbed (q : Set.Ioo (0 : ℝ) 1) (s : Sector) (p : Polynomial ℝ) :
    radialL2 q s (SemiclassicalJacobiEnergy.polynomialEmbed (parameter q) s.val.2 p) =
      radialPolynomialL2 q s p := by
  apply LinearMap.extendOfNorm_eq
  · exact SemiclassicalJacobiEnergy.denseRange_polynomialEmbed _ _
  · exact ⟨_, norm_radialPolynomialL2_le q s⟩

/-- Poincaré persists on the completed radial energy space. -/
theorem norm_radialL2_le (q : Set.Ioo (0 : ℝ) 1) (s : Sector) (x : RadialSpace q s) :
    ‖radialL2 q s x‖ ≤ Real.sqrt (annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2) * ‖x‖ :=
  LinearMap.norm_extendOfNorm_apply_le
    (SemiclassicalJacobiEnergy.denseRange_polynomialEmbed _ _) _ (norm_radialPolynomialL2_le q s) x

/-- Ordinary physical L² body forcing defines a continuous load on a completed radial sector. -/
def radialBodyLoad (q : Set.Ioo (0 : ℝ) 1) (s : Sector) (g : PhysicalL2 q) :
    RadialSpace q s →L[ℝ] ℝ :=
  (InnerProductSpace.toDual ℝ (PhysicalL2 q) g).comp (radialL2 q s)

/-- On polynomial trials, the body load is the usual physical integral. -/
theorem radialBodyLoad_polynomialEmbed (q : Set.Ioo (0 : ℝ) 1) (s : Sector)
    (g : PhysicalL2 q) (p : Polynomial ℝ) :
    radialBodyLoad q s g (SemiclassicalJacobiEnergy.polynomialEmbed (parameter q) s.val.2 p) =
      ∫ v, g v * eval v (annulusDirichletLift ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2 p)
        ∂zernikeAnnulusMeasure ⟨q,q.property.1.le,q.property.2⟩ 0 0 := by
  rw [radialBodyLoad_apply, radialL2_polynomialEmbed, L2.inner_def]
  apply integral_congr_ae
  filter_upwards [MvPolynomial.coe_polynomialToL2 (physicalMeasure_moments q)
    (annulusDirichletLift ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2 p)] with v hv
  change inner ℝ (g v) ((MvPolynomial.polynomialToL2 (physicalMeasure_moments q)
    (annulusDirichletLift ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2 p)) v) = _
  simp [hv, mul_comm]

/-- The body-load bound includes exactly the physical angular scale. -/
theorem norm_radialBodyLoad_le (q : Set.Ioo (0 : ℝ) 1) (s : Sector) (g : PhysicalL2 q) :
    ‖radialBodyLoad q s g‖ ≤ ‖g‖ * Real.sqrt (annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2) := by
  apply ContinuousLinearMap.opNorm_le_bound _ (by positivity)
  intro x
  rw [radialBodyLoad_apply]
  exact (norm_inner_le_norm g _).trans ((mul_le_mul_of_nonneg_left
    (norm_radialL2_le q s x) (norm_nonneg g)).trans_eq (mul_assoc _ _ _).symm)

/-- The physical Poisson load expressed in the solver's unscaled radial energy coordinates. -/
def radialPoissonLoad (q : Set.Ioo (0 : ℝ) 1) (s : Sector) (g : PhysicalL2 q) :
    RadialSpace q s →L[ℝ] ℝ :=
  (annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2)⁻¹ • radialBodyLoad q s g

/-- The radial solution satisfies the correctly scaled physical weak equation in its sector. -/
theorem radialPoisson_solution_weak (q : Set.Ioo (0 : ℝ) 1) (s : Sector)
    (g : PhysicalL2 q) (z : RadialSpace q s) :
    annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2 * inner ℝ
      (SemiclassicalJacobiEnergy.solution (parameter q) s.val.2 (radialPoissonLoad q s g)) z =
      inner ℝ g (radialL2 q s z) := by
  rw [SemiclassicalJacobiEnergy.solution_weak]
  change annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2 *
    ((annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2)⁻¹ * inner ℝ g (radialL2 q s z)) = _
  rw [← mul_assoc, mul_inv_cancel₀
    (annulusEnergyScale_pos q.property s.val.1 s.val.2 s.property).ne', one_mul]

/-- The existing radial coefficient solves converge in physical annular L². -/
theorem tendsto_radialL2_coefficients (q : Set.Ioo (0 : ℝ) 1) (s : Sector)
    (ℓ : RadialSpace q s →L[ℝ] ℝ) :
    Tendsto (fun N ↦ radialL2 q s
      (SemiclassicalJacobiEnergy.trialEmbed (parameter q) s.val.2 N
        (SemiclassicalJacobiEnergy.coefficients (parameter q) s.val.2 N ℓ))) atTop
      (𝓝 (radialL2 q s (SemiclassicalJacobiEnergy.solution (parameter q) s.val.2 ℓ))) :=
  (radialL2 q s).continuous.tendsto _ |>.comp
    (SemiclassicalJacobiEnergy.tendsto_trialEmbed_coefficients (parameter q) s.val.2 ℓ)

/-- A radial energy vector whose physical L² class vanishes is zero. -/
theorem radialL2_injective (q : Set.Ioo (0 : ℝ) 1) (s : Sector) :
    Function.Injective (radialL2 q s) := by
  apply (LinearMap.ker_eq_bot (f := (radialL2 q s).toLinearMap)).mp
  rw [LinearMap.ker_eq_bot']
  intro x hx
  apply ext_inner_right ℝ
  intro y
  simp only [inner_zero_left]
  refine (SemiclassicalJacobiEnergy.denseRange_polynomialEmbed (parameter q) s.val.2).induction
    ?_ (isClosed_eq (by fun_prop) continuous_const) y
  rintro _ ⟨p, rfl⟩
  have he := inner_radialL2_polynomialLaplacian q s x p
  change radialL2 q s x = 0 at hx
  rw [hx, inner_zero_left, neg_zero] at he
  exact (mul_eq_zero.mp he).resolve_left
    (annulusEnergyScale_pos q.property s.val.1 s.val.2 s.property).ne'

/-- Physical L² realization of a normalized angular energy coordinate. -/
def sectorL2 (q : Set.Ioo (0 : ℝ) 1) (s : Sector) : RadialSpace q s →L[ℝ] PhysicalL2 q :=
  (Real.sqrt (annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2))⁻¹ • radialL2 q s

/-- Every normalized sector realization is a contraction. -/
theorem norm_sectorL2_le (q : Set.Ioo (0 : ℝ) 1) (s : Sector) (x : RadialSpace q s) :
    ‖sectorL2 q s x‖ ≤ ‖x‖ := by
  have hs : 0 < Real.sqrt (annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2) :=
    Real.sqrt_pos.mpr (annulusEnergyScale_pos q.property s.val.1 s.val.2 s.property)
  simp only [sectorL2, smul_apply, norm_smul, Real.norm_eq_abs,
    abs_of_nonneg (inv_nonneg.mpr hs.le)]
  calc
    _ ≤ (Real.sqrt (annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2))⁻¹ *
        (Real.sqrt (annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2) * ‖x‖) :=
      mul_le_mul_of_nonneg_left (norm_radialL2_le q s x) (inv_nonneg.mpr hs.le)
    _ = ‖x‖ := by rw [← mul_assoc, inv_mul_cancel₀ hs.ne', one_mul]

/-- The physical L² realization of the full completed annular energy space. -/
def physicalL2 (q : Set.Ioo (0 : ℝ) 1) : Space q →L[ℝ] PhysicalL2 q :=
  lp.orthogonalSynthesis (sectorL2 q) (norm_sectorL2_le q) (inner_sectorL2_eq_zero q)

/-- The full energy completion has a faithful physical L² realization. -/
theorem physicalL2_injective (q : Set.Ioo (0 : ℝ) 1) :
    Function.Injective (physicalL2 q) :=
  lp.orthogonalSynthesis_injective _ _ _ (sectorL2_injective q)

/-- The normalized physical sector series converges unconditionally in L². -/
theorem hasSum_physicalL2 (q : Set.Ioo (0 : ℝ) 1) (x : Space q) :
    HasSum (fun s ↦ sectorL2 q s (x s)) (physicalL2 q x) :=
  lp.hasSum_orthogonalSynthesis _ _ _ x

/-- The global Poincaré bound persists on the full angular energy completion. -/
theorem norm_physicalL2_le (q : Set.Ioo (0 : ℝ) 1) (x : Space q) :
    ‖physicalL2 q x‖ ≤ ‖x‖ := lp.norm_orthogonalSynthesis_le _ _ _ x

/-- Physical polynomial embedding cancels the coordinate normalization exactly. -/
theorem physicalL2_polynomialEmbed (q : Set.Ioo (0 : ℝ) 1) (s : Sector) (p : Polynomial ℝ) :
    physicalL2 q (polynomialEmbed q s p) = radialPolynomialL2 q s p := by
  have hs : Real.sqrt (annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2) ≠ 0 :=
    (Real.sqrt_pos.mpr (annulusEnergyScale_pos q.property s.val.1 s.val.2 s.property)).ne'
  simp only [polynomialEmbed, LinearMap.smul_apply, LinearMap.comp_apply,
    ContinuousLinearMap.coe_coe, map_smul, physicalL2_sectorEmbed, sectorL2,
    smul_apply, smul_smul, mul_inv_cancel₀ hs, one_smul]
  exact radialL2_polynomialEmbed q s p

/-- A physical L² body force acts by pairing with the realized energy vector. -/
def bodyLoad (q : Set.Ioo (0 : ℝ) 1) (g : PhysicalL2 q) : Space q →L[ℝ] ℝ :=
  (InnerProductSpace.toDual ℝ (PhysicalL2 q) g).comp (physicalL2 q)

/-- Against a physical polynomial trial, the load is its ordinary area integral. -/
theorem bodyLoad_polynomialEmbed (q : Set.Ioo (0 : ℝ) 1) (g : PhysicalL2 q)
    (s : Sector) (p : Polynomial ℝ) :
    bodyLoad q g (polynomialEmbed q s p) =
      ∫ v, g v * eval v (annulusDirichletLift ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2 p)
        ∂zernikeAnnulusMeasure ⟨q,q.property.1.le,q.property.2⟩ 0 0 := by
  rw [bodyLoad_apply, physicalL2_polynomialEmbed]
  rw [← radialBodyLoad_polynomialEmbed q s g p, radialBodyLoad_apply, radialL2_polynomialEmbed]

/-- The radial load passed to the assembled solver has the inverse square-root normalization. -/
theorem sectorLoad_bodyLoad (q : Set.Ioo (0 : ℝ) 1) (g : PhysicalL2 q) (s : Sector) :
    sectorLoad q (bodyLoad q g) s =
      (Real.sqrt (annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2))⁻¹ • radialBodyLoad q s g := by
  ext x
  change inner ℝ g (physicalL2 q (sectorEmbed q s x)) = _
  rw [physicalL2_sectorEmbed]
  simp [sectorL2, inner_smul_right, radialBodyLoad_apply]

/-- The physical body-force problem has a unique solution in the chosen energy completion. -/
theorem existsUnique_body_solution (q : Set.Ioo (0 : ℝ) 1) (g : PhysicalL2 q) :
    ∃! x : Space q, ∀ z : Space q, inner ℝ x z = inner ℝ g (physicalL2 q z) :=
  existsUnique_solution q (bodyLoad q g)

/-- The full body-force solution is bounded by the physical L² norm of its datum. -/
theorem norm_solution_bodyLoad_le (q : Set.Ioo (0 : ℝ) 1) (g : PhysicalL2 q) :
    ‖solution q (bodyLoad q g)‖ ≤ ‖g‖ :=
  ((InnerProductSpace.toDual ℝ (Space q)).symm.norm_map _).le.trans (norm_bodyLoad_le q g)

/-- Finite body-force solutions are uniformly bounded in the full energy norm. -/
theorem norm_finiteSolve_bodyLoad_le (q : Set.Ioo (0 : ℝ) 1)
    (c : Finset Sector × ℕ) (g : PhysicalL2 q) :
    ‖finiteSolve q c (bodyLoad q g)‖ ≤ ‖g‖ :=
  (norm_finiteSolve_le q c _).trans (norm_bodyLoad_le q g)

/-- Physical L² stability is uniform in both truncation parameters. -/
theorem norm_physicalL2_finiteSolve_bodyLoad_le (q : Set.Ioo (0 : ℝ) 1)
    (c : Finset Sector × ℕ) (g : PhysicalL2 q) :
    ‖physicalL2 q (finiteSolve q c (bodyLoad q g))‖ ≤ ‖g‖ :=
  (norm_physicalL2_le q _).trans (norm_finiteSolve_bodyLoad_le q c g)

/-- The realized finite solution is precisely the sum of its computed radial reconstructions. -/
theorem physicalL2_finiteSolve (q : Set.Ioo (0 : ℝ) 1) (c : Finset Sector × ℕ)
    (ℓ : Space q →L[ℝ] ℝ) :
    physicalL2 q (finiteSolve q c ℓ) =
      ∑ s ∈ c.1, sectorL2 q s
        (SemiclassicalJacobiEnergy.trialEmbed (parameter q) s.val.2 c.2
          (SemiclassicalJacobiEnergy.coefficients (parameter q) s.val.2 c.2 (sectorLoad q ℓ s))) := by
  simp only [finiteSolve, map_sum, physicalL2_sectorEmbed]

/-- Energy best approximation controls the physical L² error of the actual finite solve. -/
theorem physicalL2_finiteSolve_bestApproximation (q : Set.Ioo (0 : ℝ) 1)
    (c : Finset Sector × ℕ) (ℓ : Space q →L[ℝ] ℝ) (y : Space q) (hy : y ∈ trialSpace q c) :
    ‖physicalL2 q (solution q ℓ) - physicalL2 q (finiteSolve q c ℓ)‖ ≤ ‖solution q ℓ-y‖ :=
  (norm_physicalL2_error_le q c ℓ).trans (finiteSolve_bestApproximation q c ℓ y hy)

/-- Joint angular and radial coefficient convergence holds in physical L². -/
theorem tendsto_physicalL2_finiteSolve (q : Set.Ioo (0 : ℝ) 1) (ℓ : Space q →L[ℝ] ℝ) :
    Tendsto (fun c : Finset Sector × ℕ ↦ physicalL2 q (finiteSolve q c ℓ)) atTop
      (𝓝 (physicalL2 q (solution q ℓ))) :=
  (physicalL2 q).continuous.tendsto _ |>.comp (tendsto_finiteSolve q ℓ)

/-- The common frequency and radial-degree cutoff also converges in physical L². -/
theorem tendsto_physicalL2_finiteSolve_frequencyCutoff (q : Set.Ioo (0 : ℝ) 1)
    (ℓ : Space q →L[ℝ] ℝ) :
    Tendsto (fun N ↦ physicalL2 q (finiteSolve q (frequencyCutoff N, N) ℓ)) atTop
      (𝓝 (physicalL2 q (solution q ℓ))) :=
  (tendsto_physicalL2_finiteSolve q ℓ).comp tendsto_frequencyCutoff

/-- Two physical L² components, with the Euclidean sum norm for Cartesian gradients. -/
abbrev GradientL2 (q : Set.Ioo (0 : ℝ) 1) := lp (fun _ : Fin 2 ↦ PhysicalL2 q) 2

/-- A polynomial's two Cartesian derivatives as physical L² classes. -/
def polynomialGradient (q : Set.Ioo (0 : ℝ) 1) :
    MvPolynomial (Fin 2) ℝ →ₗ[ℝ] GradientL2 q :=
  (lp.singleContinuousLinearMap ℝ (fun _ : Fin 2 ↦ PhysicalL2 q) 2 0).toLinearMap.comp
      ((MvPolynomial.polynomialToL2 (physicalMeasure_moments q)).comp (pderiv 0).toLinearMap) +
    (lp.singleContinuousLinearMap ℝ (fun _ : Fin 2 ↦ PhysicalL2 q) 2 1).toLinearMap.comp
      ((MvPolynomial.polynomialToL2 (physicalMeasure_moments q)).comp (pderiv 1).toLinearMap)

/-- The vector-field inner product is the integral of the Cartesian gradient dot product. -/
theorem inner_polynomialGradient (q : Set.Ioo (0 : ℝ) 1) (p z : MvPolynomial (Fin 2) ℝ) :
    inner ℝ (polynomialGradient q p) (polynomialGradient q z) =
      integralInner (zernikeAnnulusMeasure ⟨q,q.property.1.le,q.property.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,q.property.1.le,q.property.2⟩ 0 0) (pderiv 0 p) (pderiv 0 z) (h := (hasFiniteMoments_zernikeAnnulusMeasure q.property (by norm_num) (by norm_num)).integrable_mul _ _)) +
      integralInner (zernikeAnnulusMeasure ⟨q,q.property.1.le,q.property.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,q.property.1.le,q.property.2⟩ 0 0) (pderiv 1 p) (pderiv 1 z) (h := (hasFiniteMoments_zernikeAnnulusMeasure q.property (by norm_num) (by norm_num)).integrable_mul _ _)) := by
  rw [lp.inner_eq_tsum, tsum_fintype, Fin.sum_univ_two]
  simp only [polynomialGradient_apply, MvPolynomial.inner_polynomialToL2, integralInner]

/-- Cartesian gradient inner products agree with scaled radial energy. -/
theorem inner_radialPolynomialGradient (q : Set.Ioo (0 : ℝ) 1) (s : Sector)
    (p z : Polynomial ℝ) :
    inner ℝ (radialPolynomialGradient q s p) (radialPolynomialGradient q s z) =
      annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2 * inner ℝ
        (SemiclassicalJacobiEnergy.polynomialEmbed (parameter q) s.val.2 p)
        (SemiclassicalJacobiEnergy.polynomialEmbed (parameter q) s.val.2 z) := by
  change inner ℝ (polynomialGradient q _) (polynomialGradient q _) = _
  rw [inner_polynomialGradient]
  exact integralInner_annulusDirichlet_gradient q.property s.val.1 s.val.2 p z

/-- The normalization removes the physical sector factor from gradient energy. -/
theorem inner_normalizedPolynomialGradient (q : Set.Ioo (0 : ℝ) 1) (s : Sector)
    (p z : Polynomial ℝ) :
    inner ℝ (normalizedPolynomialGradient q s p) (normalizedPolynomialGradient q s z) =
      inner ℝ (SemiclassicalJacobiEnergy.polynomialEmbed (parameter q) s.val.2 p)
        (SemiclassicalJacobiEnergy.polynomialEmbed (parameter q) s.val.2 z) := by
  have hs := Real.sq_sqrt (annulusEnergyScale_pos q.property s.val.1 s.val.2 s.property).le
  have hn : Real.sqrt (annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2) ≠ 0 :=
    (Real.sqrt_pos.mpr (annulusEnergyScale_pos q.property s.val.1 s.val.2 s.property)).ne'
  simp only [normalizedPolynomialGradient, LinearMap.smul_apply, real_inner_smul_left,
    inner_smul_right, inner_radialPolynomialGradient]
  field_simp
  rw [hs]

/-- Completed normalized sector energy is realized isometrically as Cartesian gradient classes. -/
def sectorGradient (q : Set.Ioo (0 : ℝ) 1) (s : Sector) :
    RadialSpace q s →ₗᵢ[ℝ] GradientL2 q :=
  (normalizedPolynomialGradient q s).extendOfIsometry
    (SemiclassicalJacobiEnergy.denseRange_polynomialEmbed _ _)
    (norm_normalizedPolynomialGradient q s)

/-- The completed sector gradient agrees with its polynomial Cartesian derivatives. -/
theorem sectorGradient_polynomialEmbed (q : Set.Ioo (0 : ℝ) 1) (s : Sector) (p : Polynomial ℝ) :
    sectorGradient q s (SemiclassicalJacobiEnergy.polynomialEmbed (parameter q) s.val.2 p) =
      normalizedPolynomialGradient q s p :=
  LinearMap.extendOfIsometry_eq _ _ _ p

/-- The full annular energy norm is the norm of its physical Cartesian gradient classes. -/
def gradientL2 (q : Set.Ioo (0 : ℝ) 1) : Space q →ₗᵢ[ℝ] GradientL2 q :=
  (orthogonalFamily_sectorGradient q).linearIsometry

/-- Cartesian gradients are the unconditionally convergent sum of their normalized sectors. -/
theorem hasSum_gradientL2 (q : Set.Ioo (0 : ℝ) 1) (x : Space q) :
    HasSum (fun s ↦ sectorGradient q s (x s)) (gradientL2 q x) :=
  (orthogonalFamily_sectorGradient q).hasSum_linearIsometry x

/-- The gradient realization preserves the full energy inner product. -/
theorem inner_gradientL2 (q : Set.Ioo (0 : ℝ) 1) (x y : Space q) :
    inner ℝ (gradientL2 q x) (gradientL2 q y) = inner ℝ x y :=
  (gradientL2 q).inner_map_map x y

/-- No completed energy vector is lost in the Cartesian gradient realization. -/
theorem gradientL2_injective (q : Set.Ioo (0 : ℝ) 1) :
    Function.Injective (gradientL2 q) := (gradientL2 q).injective

/-- The completed Cartesian gradient space is closed in the physical vector-field L² space. -/
theorem isClosed_range_gradientL2 (q : Set.Ioo (0 : ℝ) 1) :
    IsClosed (Set.range (gradientL2 q)) :=
  (gradientL2 q).isometry.isClosedEmbedding.isClosed_range

/-- Physical polynomial energy vectors realize precisely their Cartesian polynomial gradients. -/
theorem gradientL2_polynomialEmbed (q : Set.Ioo (0 : ℝ) 1) (s : Sector) (p : Polynomial ℝ) :
    gradientL2 q (polynomialEmbed q s p) =
      polynomialGradient q (annulusDirichletLift ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2 p) := by
  have hs : Real.sqrt (annulusEnergyScale ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2) ≠ 0 :=
    (Real.sqrt_pos.mpr (annulusEnergyScale_pos q.property s.val.1 s.val.2 s.property)).ne'
  simp only [polynomialEmbed, LinearMap.smul_apply, LinearMap.comp_apply,
    ContinuousLinearMap.coe_coe, map_smul, gradientL2_sectorEmbed, sectorGradient_polynomialEmbed,
    normalizedPolynomialGradient, smul_smul, mul_inv_cancel₀ hs, one_smul]
  rfl

/-- The simultaneous function and gradient realization is an isometry for the product max norm. -/
def functionGradient (q : Set.Ioo (0 : ℝ) 1) :
    Space q →ₗᵢ[ℝ] PhysicalL2 q × GradientL2 q where
  __ := (physicalL2 q).toLinearMap.prod (gradientL2 q).toLinearMap
  norm_map' x := by
    change max ‖physicalL2 q x‖ ‖gradientL2 q x‖ = ‖x‖
    rw [LinearIsometry.norm_map, max_eq_right (norm_physicalL2_le q x)]

/-- The simultaneous physical function/gradient graph is closed. -/
theorem isClosed_range_functionGradient (q : Set.Ioo (0 : ℝ) 1) :
    IsClosed (Set.range (functionGradient q)) :=
  (functionGradient q).isometry.isClosedEmbedding.isClosed_range

/-- The completed physical body-force equation is a pairing of physical gradient classes. -/
theorem gradient_body_solution_weak (q : Set.Ioo (0 : ℝ) 1) (g : PhysicalL2 q) (z : Space q) :
    inner ℝ (gradientL2 q (solution q (bodyLoad q g))) (gradientL2 q z) =
      inner ℝ g (physicalL2 q z) := by
  rw [inner_gradientL2]
  exact body_solution_weak q g z

/-- The actual finite reconstructions have precisely the energy error in Cartesian gradient L². -/
theorem norm_gradientL2_error (q : Set.Ioo (0 : ℝ) 1) (c : Finset Sector × ℕ)
    (ℓ : Space q →L[ℝ] ℝ) :
    ‖gradientL2 q (solution q ℓ) - gradientL2 q (finiteSolve q c ℓ)‖ =
      ‖solution q ℓ - finiteSolve q c ℓ‖ := by
  rw [← map_sub, LinearIsometry.norm_map]

/-- Common angular and radial cutoffs converge simultaneously as functions and Cartesian gradients. -/
theorem tendsto_functionGradient_finiteSolve (q : Set.Ioo (0 : ℝ) 1)
    (ℓ : Space q →L[ℝ] ℝ) :
    Tendsto (fun N ↦ functionGradient q (finiteSolve q (frequencyCutoff N, N) ℓ)) atTop
      (𝓝 (functionGradient q (solution q ℓ))) :=
  (functionGradient q).continuous.tendsto _ |>.comp (tendsto_finiteSolve_frequencyCutoff q ℓ)

/-- Pairing a polynomial with a physical test gives the ordinary Lebesgue integral. -/
theorem inner_polynomial_physicalTestL2 (q : Ioo (0 : ℝ) 1)
    (p : MvPolynomial (Fin 2) ℝ) (φ : TestFunction (domain q) ℝ ⊤) :
    inner ℝ (polynomialToL2 (physicalMeasure_moments q) p) (physicalTestL2 q φ) =
      ∫ x, eval x p * φ x := by
  rw [L2.inner_def]
  calc
    _ = ∫ x, eval x p*φ x ∂zernikeAnnulusMeasure ⟨q,q.property.1.le,q.property.2⟩ 0 0 := by
      apply integral_congr_ae
      filter_upwards [coe_polynomialToL2 (physicalMeasure_moments q) p,
        coe_physicalTestL2 q φ] with x hp hφ
      simp [hp, hφ, mul_comm]
    _ = _ := by
      rw [zernikeAnnulusMeasure_zero_zero_eq_volume q.property]
      apply setIntegral_eq_integral_of_forall_compl_eq_zero
      intro x hx
      have hφ : φ x = 0 := by
        by_contra hn
        exact hx (φ.tsupport_subset (subset_tsupport φ (Function.mem_support.mpr hn)))
      rw [hφ, mul_zero]

/-- Polynomial Cartesian gradients satisfy the weak-derivative identity on the annulus. -/
theorem polynomial_weakDerivative (q : Ioo (0 : ℝ) 1) (p : MvPolynomial (Fin 2) ℝ)
    (φ : TestFunction (domain q) ℝ ⊤) (i : Fin 2) :
    inner ℝ (polynomialToL2 (physicalMeasure_moments q) p)
      (physicalTestL2 q (DirichletSpace.testDerivative (domain q) i φ)) =
      -inner ℝ (polynomialGradient q p i) (physicalTestL2 q φ) := by
  rw [polynomialGradient_apply, inner_polynomial_physicalTestL2,
    inner_polynomial_physicalTestL2]
  exact integral_eval_mul_testDerivative (domain q) p φ i

/-- Every completed annular energy vector has the prescribed Cartesian weak derivatives. -/
theorem gradientL2_weakDerivative (q : Ioo (0 : ℝ) 1) (x : Space q)
    (φ : TestFunction (domain q) ℝ ⊤) (i : Fin 2) :
    inner ℝ (physicalL2 q x)
      (physicalTestL2 q (DirichletSpace.testDerivative (domain q) i φ)) =
      -inner ℝ (gradientL2 q x i) (physicalTestL2 q φ) := by
  let L : Space q →L[ℝ] ℝ :=
    (innerSLFlip ℝ (physicalTestL2 q (DirichletSpace.testDerivative (domain q) i φ))).comp
      (physicalL2 q)
  let R : Space q →L[ℝ] ℝ :=
    -((innerSLFlip ℝ (physicalTestL2 q φ)).comp
      ((lp.evalCLM ℝ (fun _ : Fin 2 ↦ PhysicalL2 q) 2 i).comp
        (gradientL2 q).toContinuousLinearMap))
  have hp (s : Sector) (p : Polynomial ℝ) :
      L (polynomialEmbed q s p) = R (polynomialEmbed q s p) := by
    have h := polynomial_weakDerivative q (annulusDirichletLift ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2 p) φ i
    simpa only [L, R, ContinuousLinearMap.comp_apply, innerSLFlip_apply_apply,
      neg_apply, lp.evalCLM, lp.evalₗ_apply, LinearMap.mkContinuous_apply,
      LinearIsometry.coe_toContinuousLinearMap, physicalL2_polynomialEmbed,
      gradientL2_polynomialEmbed, radialPolynomialL2, LinearMap.comp_apply] using h
  have hLR : L = R := by
    apply lp.ext_continuousLinearMap (by norm_num : (2 : ENNReal) ≠ ⊤)
    intro s
    apply ContinuousLinearMap.ext
    intro y
    change L (sectorEmbed q s y) = R (sectorEmbed q s y)
    refine (SemiclassicalJacobiEnergy.denseRange_polynomialEmbed (parameter q) s.val.2).induction
      ?_ (isClosed_eq (by fun_prop) (by fun_prop)) y
    rintro _ ⟨p, rfl⟩
    have h := hp s p
    simp only [polynomialEmbed, LinearMap.smul_apply, LinearMap.comp_apply,
      ContinuousLinearMap.coe_coe, map_smul, smul_eq_mul] at h
    exact mul_left_cancel₀
      (Real.sqrt_pos.mpr (annulusEnergyScale_pos q.property s.val.1 s.val.2 s.property)).ne' h
  have h := congrArg (fun T : Space q →L[ℝ] ℝ ↦ T x) hLR
  simpa only [L, R, ContinuousLinearMap.comp_apply, innerSLFlip_apply_apply,
    neg_apply, lp.evalCLM, lp.evalₗ_apply, LinearMap.mkContinuous_apply,
    LinearIsometry.coe_toContinuousLinearMap] using h

/-- Positive cutoff scales give compact support strictly inside the annulus. -/
theorem cutoffTrial_support (q : Ioo (0 : ℝ) 1) {a : ℝ} (ha : 0 < a)
    (p : MvPolynomial (Fin 2) ℝ) :
    HasCompactSupport (cutoffTrial q a p) ∧ tsupport (cutoffTrial q a p) ⊆ domain q := by
  let K : Set E := {x | (q : ℝ)+1/a ≤ x 0^2+x 1^2 ∧ x 0^2+x 1^2 ≤ 1-1/a}
  have hc : IsClosed K := by
    have hr : Continuous (fun x : E ↦ x 0^2+x 1^2) := by fun_prop
    exact (isClosed_le continuous_const hr).inter (isClosed_le hr continuous_const)
  have hs : tsupport (cutoffTrial q a p) ⊆ K := closure_minimal (cutoff_support ha p) hc
  have hk : K ⊆ box := by
    intro x hx
    apply radius_mem_box
    have h := hx.2
    have : 0 < 1/a := by positivity
    linarith
  refine ⟨(isCompact_Icc.of_isClosed_subset (isClosed_tsupport _) (hs.trans hk)), ?_⟩
  intro x hx
  have h := hs hx
  change (q : ℝ) < x 0^2+x 1^2 ∧ x 0^2+x 1^2 < 1
  have : 0 < 1/a := by positivity
  constructor <;> linarith [h.1,h.2]

/-- At every interior point the approximants eventually agree locally with the trial. -/
theorem cutoffTest_eventuallyEq (q : Ioo (0 : ℝ) 1) (p : MvPolynomial (Fin 2) ℝ)
    {x : E} (hx : x ∈ domain q) :
    ∀ᶠ n : ℕ in atTop,
      (cutoffTest q p n : E → ℝ) =ᶠ[𝓝 x] (fun y ↦ eval y (annulusBoundaryWeight q*p)) := by
  have h0 : 0 < x 0^2+x 1^2-(q : ℝ) := sub_pos.mpr hx.1
  have h1 : 0 < 1-(x 0^2+x 1^2) := sub_pos.mpr hx.2
  filter_upwards [eventually_scaled_gt h0,eventually_scaled_gt h1] with n hn0 hn1
  have hc0 : Continuous (fun y : E ↦ ((n : ℝ)+1)*(y 0^2+y 1^2-(q : ℝ))) := by fun_prop
  have hc1 : Continuous (fun y : E ↦ ((n : ℝ)+1)*(1-(y 0^2+y 1^2))) := by fun_prop
  filter_upwards [hc0.continuousAt.eventually (eventually_gt_nhds hn0),
    hc1.continuousAt.eventually (eventually_gt_nhds hn1)] with y hy0 hy1
  change cutoffTrial q ((n : ℝ)+1) p y = _
  simp only [cutoffTrial,boundaryFactor_eq hy0.le,boundaryFactor_eq hy1.le,
    annulusBoundaryWeight,eval_mul,eval_sub,map_one,eval_add,eval_pow,eval_X,eval_C]
  ring

/-- Each Cartesian first derivative converges in L² to the polynomial derivative. -/
theorem tendsto_cutoffTest_derivative_L2 (q : Ioo (0 : ℝ) 1)
    (p : MvPolynomial (Fin 2) ℝ) (i : Fin 2) :
    Tendsto (fun n ↦ physicalTestL2 q (DirichletSpace.testDerivative (domain q) i
      (cutoffTest q p n))) atTop
      (𝓝 (polynomialToL2 (physicalMeasure_moments q) (pderiv i (annulusBoundaryWeight q*p)))) := by
  obtain ⟨C,_,hC⟩ := cutoffTrial_uniform_bound q p
  apply tendsto_testL2_of_bound q
    (fun n ↦ DirichletSpace.testDerivative (domain q) i (cutoffTest q p n))
    (pderiv i (annulusBoundaryWeight q*p)) (C := C)
  · intro n x hx
    rw [DirichletSpace.testDerivative_apply]
    have h := (fderiv ℝ (cutoffTest q p n) x).le_opNorm (Pi.single i 1)
    simp only [Pi.norm_single, norm_one, mul_one] at h
    exact h.trans (hC n x (radius_mem_box hx.2.le)).2
  · intro x hx
    apply (tendsto_congr' ?_).mpr tendsto_const_nhds
    filter_upwards [cutoffTest_eventuallyEq q p hx] with n hn
    rw [DirichletSpace.testDerivative_apply, hn.fderiv_eq, fderiv_eval_single]

/-- The equality of polar and area measures identifies their function/gradient spaces. -/
def areaGraphEquiv (q : Ioo (0 : ℝ) 1) :
    PhysicalL2 q × GradientL2 q ≃ₗᵢ[ℝ]
      DirichletSpace.FunctionL2 (domain q) × DirichletSpace.GradientL2 (domain q) :=
  L2.graphMeasureCongr 2 (zernikeAnnulusMeasure_zero_zero_eq_volume q.property)

/-- The Cartesian gradient of a smooth test in the physical L² spaces. -/
def physicalTestGradient (q : Ioo (0 : ℝ) 1) :
    TestFunction (domain q) ℝ ⊤ →ₗ[ℝ] GradientL2 q :=
  ∑ i : Fin 2, (lp.singleContinuousLinearMap ℝ (fun _ : Fin 2 ↦ PhysicalL2 q) 2 i).toLinearMap.comp
    ((physicalTestL2 q).comp (DirichletSpace.testDerivative (domain q) i).toLinearMap)

/-- The full Cartesian gradient of the smooth cutoffs converges in the Hilbert sum. -/
theorem tendsto_cutoffTest_gradient_L2 (q : Ioo (0 : ℝ) 1)
    (p : MvPolynomial (Fin 2) ℝ) :
    Tendsto (fun n ↦ physicalTestGradient q (cutoffTest q p n)) atTop
      (𝓝 (polynomialGradient q (annulusBoundaryWeight q*p))) := by
  have h : Tendsto (fun n ↦ ∑ i : Fin 2, lp.single 2 i
      (physicalTestL2 q (DirichletSpace.testDerivative (domain q) i (cutoffTest q p n)))) atTop
      (𝓝 (∑ i : Fin 2, lp.single 2 i
        (polynomialToL2 (physicalMeasure_moments q) (pderiv i (annulusBoundaryWeight q*p))))) := by
    apply tendsto_finsetSum
    intro i _
    exact (lp.singleContinuousLinearMap ℝ (fun _ : Fin 2 ↦ PhysicalL2 q) 2 i).continuous.tendsto _
      |>.comp (tendsto_cutoffTest_derivative_L2 q p i)
  simpa [physicalTestGradient, polynomialGradient, Fin.sum_univ_two] using h

/-- Every boundary-vanishing Cartesian polynomial belongs to the zero-boundary graph. -/
theorem polynomialGraph_mem_dirichlet (q : Ioo (0 : ℝ) 1)
    (p : MvPolynomial (Fin 2) ℝ) :
    areaGraphEquiv q (polynomialToL2 (physicalMeasure_moments q) (annulusBoundaryWeight q*p),
      polynomialGradient q (annulusBoundaryWeight q*p)) ∈ DirichletSpace.graph (domain q) := by
  have ht := (areaGraphEquiv q).continuous.tendsto _ |>.comp
    ((tendsto_cutoffTest_L2 q p).prodMk_nhds (tendsto_cutoffTest_gradient_L2 q p))
  apply (DirichletSpace.graph_isClosed (domain q)).mem_of_tendsto ht
  exact Filter.Eventually.of_forall (fun n ↦ by
    dsimp only [Function.comp_apply]
    rw [areaGraphEquiv_testGraph]
    exact DirichletSpace.testGraph_mem (domain q) _)

/-- Every completed annular energy vector belongs to the zero-boundary Sobolev graph. -/
theorem functionGradient_mem_dirichlet (q : Ioo (0 : ℝ) 1) (x : Space q) :
    areaGraphEquiv q (functionGradient q x) ∈ DirichletSpace.graph (domain q) := by
  let S : Submodule ℝ (Space q) := (DirichletSpace.graph (domain q)).comap
    ((areaGraphEquiv q).toLinearMap.comp (functionGradient q).toLinearMap)
  have hS : IsClosed (S : Set (Space q)) :=
    (DirichletSpace.graph_isClosed (domain q)).preimage
      ((areaGraphEquiv q).continuous.comp (functionGradient q).continuous)
  have hp (s : Sector) (p : Polynomial ℝ) : polynomialEmbed q s p ∈ S := by
    have h := polynomialGraph_mem_dirichlet q (annulusLift ⟨q,ne_of_lt q.property.2⟩ s.val.1 s.val.2 p)
    change areaGraphEquiv q (physicalL2 q (polynomialEmbed q s p),
      gradientL2 q (polynomialEmbed q s p)) ∈ DirichletSpace.graph (domain q)
    rw [physicalL2_polynomialEmbed, gradientL2_polynomialEmbed]
    exact h
  have hs (s : Sector) (y : RadialSpace q s) : sectorEmbed q s y ∈ S := by
    refine (Polynomial.SemiclassicalJacobiEnergy.denseRange_polynomialEmbed
      (parameter q) s.val.2).induction ?_ (hS.preimage (sectorEmbed q s).continuous) y
    rintro _ ⟨p, rfl⟩
    have h := hp s p
    simp only [polynomialEmbed, LinearMap.smul_apply, LinearMap.comp_apply,
      ContinuousLinearMap.coe_coe] at h
    exact (S.smul_mem_iff
      (Real.sqrt_pos.mpr (annulusEnergyScale_pos q.property s.val.1 s.val.2 s.property)).ne').mp h
  apply hS.mem_of_tendsto (lp.hasSum_single (by norm_num : (2 : ENNReal) ≠ ⊤) x)
  exact Filter.Eventually.of_forall (fun t ↦ S.sum_mem (fun i _ ↦ hs i (x i)))

/-- Boundary-vanishing polynomials approximate every smooth annular test and its gradient in L². -/
theorem exists_boundary_polynomial_approximation (q : Ioo (0 : ℝ) 1)
    (φ : TestFunction (domain q) ℝ ⊤) :
    ∃ p : ℕ → MvPolynomial (Fin 2) ℝ,
      Tendsto (fun n ↦ (polynomialToL2 (physicalMeasure_moments q) (annulusBoundaryWeight q*p n),
        polynomialGradient q (annulusBoundaryWeight q*p n))) atTop
        (𝓝 (physicalTestL2 q φ, physicalTestGradient q φ)) := by
  obtain ⟨g,hg,he⟩ := exists_smooth_boundary_quotient q φ
  choose p hp using (fun n : ℕ ↦ exists_polynomial_C1_near_on_square hg
    (show (0 : ℝ) < 1/((n : ℝ)+1) by positivity))
  obtain ⟨M,hM,hB⟩ := boundary_bound q
  have hb (n : ℕ) (x : E) (hx : x ∈ box) := boundary_mul_error q φ hg he (p n) hM
    (by positivity : (0 : ℝ) ≤ 1/((n : ℝ)+1)) hB
    (fun x hx ↦ ⟨(hp n x hx).1.le,fun i ↦ ((hp n x hx).2 i).le⟩) x hx
  have ht0 : Tendsto (fun n ↦ polynomialToL2 (physicalMeasure_moments q) (annulusBoundaryWeight q*p n))
      atTop (𝓝 (physicalTestL2 q φ)) :=
    tendsto_polynomial_L2_of_error q (test_memLp q φ) _ (by positivity : 0 ≤ 2*M)
      (fun n x hx ↦ (hb n x hx).1)
  have ht (i : Fin 2) : Tendsto
      (fun n ↦ polynomialToL2 (physicalMeasure_moments q) (pderiv i (annulusBoundaryWeight q*p n)))
      atTop (𝓝 (physicalTestL2 q (DirichletSpace.testDerivative (domain q) i φ))) :=
    tendsto_polynomial_L2_of_error q (test_memLp q _) _ (by positivity : 0 ≤ 2*M)
      (fun n x hx ↦ (hb n x hx).2 i)
  have htG : Tendsto (fun n ↦ polynomialGradient q (annulusBoundaryWeight q*p n)) atTop
      (𝓝 (physicalTestGradient q φ)) := by
    have h : Tendsto (fun n ↦ ∑ i : Fin 2, lp.single 2 i
        (polynomialToL2 (physicalMeasure_moments q) (pderiv i (annulusBoundaryWeight q*p n)))) atTop
        (𝓝 (∑ i : Fin 2, lp.single 2 i (physicalTestL2 q (DirichletSpace.testDerivative (domain q) i φ)))) := by
      apply tendsto_finsetSum
      intro i _
      exact (lp.singleContinuousLinearMap ℝ (fun _ : Fin 2 ↦ PhysicalL2 q) 2 i).continuous.tendsto _
        |>.comp (ht i)
    simpa [physicalTestGradient,polynomialGradient,Fin.sum_univ_two] using h
  exact ⟨p,ht0.prodMk_nhds htG⟩

/-- Every smooth annular test graph belongs to the completed polynomial energy image. -/
theorem physicalTestGraph_mem_range (q : Ioo (0 : ℝ) 1)
    (φ : TestFunction (domain q) ℝ ⊤) :
    (physicalTestL2 q φ,physicalTestGradient q φ) ∈ Set.range (functionGradient q) := by
  obtain ⟨p,hp⟩ := exists_boundary_polynomial_approximation q φ
  apply (isClosed_range_functionGradient q).mem_of_tendsto hp
  exact Filter.Eventually.of_forall (fun n ↦ exists_functionGradient_polynomial q (p n))

/-- A smooth cutoff of a linearly vanishing factor at scale a. -/
def boundaryFactor (a r : ℝ) : ℝ := r * smoothTransition (a*r-1)

/-- The cutoff boundary factor is smooth. -/
theorem boundaryFactor_contDiff (a : ℝ) : ContDiff ℝ ∞ (boundaryFactor a) := by
  unfold boundaryFactor
  exact contDiff_id.mul (smoothTransition.contDiff.comp ((contDiff_const.mul contDiff_id).sub contDiff_const))

/-- The derivative of the cutoff boundary factor. -/
theorem hasDerivAt_boundaryFactor (a r : ℝ) : HasDerivAt (boundaryFactor a)
    (smoothTransition (a*r-1) + (a*r)*deriv smoothTransition (a*r-1)) r := by
  have ht : HasDerivAt smoothTransition (deriv smoothTransition (a*r-1)) (a*r-1) := ((smoothTransition.contDiff (n := 1)).differentiable (by norm_num)).differentiableAt.hasDerivAt
  unfold boundaryFactor
  convert (hasDerivAt_id r).mul
    (ht.comp r (((hasDerivAt_id r).const_mul a).sub_const 1)) using 1
  · rfl
  · simp only [one_mul, mul_one, Function.comp_apply, id_eq]
    ring

/-- Cutting off does not increase the absolute value of the boundary factor. -/
theorem abs_boundaryFactor_le (a r : ℝ) : |boundaryFactor a r| ≤ |r| := by
  rw [boundaryFactor, abs_mul, abs_of_nonneg (smoothTransition.nonneg _)]
  exact mul_le_of_le_one_right (abs_nonneg _) (smoothTransition.le_one _)

/-- The derivative bound is independent of the cutoff scale. -/
theorem abs_deriv_boundaryFactor_le {M : ℝ}
    (hM : ∀ x, |x * deriv smoothTransition (x-1)| ≤ 2*M) (a r : ℝ) :
    |deriv (boundaryFactor a) r| ≤ 1+2*M := by
  rw [(hasDerivAt_boundaryFactor a r).deriv]
  exact (abs_add_le _ _).trans (add_le_add
    (by rw [abs_of_nonneg (smoothTransition.nonneg _)]; exact smoothTransition.le_one _) (hM (a*r)))

/-- The cutoff factor vanishes in the boundary layer. -/
theorem boundaryFactor_zero {a r : ℝ} (h : a*r ≤ 1) : boundaryFactor a r = 0 := by
  rw [boundaryFactor, smoothTransition.zero_of_nonpos (by linarith), mul_zero]

/-- Outside the transition layer the cutoff equals the original factor. -/
theorem boundaryFactor_eq {a r : ℝ} (h : 2 ≤ a*r) : boundaryFactor a r = r := by
  rw [boundaryFactor, smoothTransition.one_of_one_le (by linarith), mul_one]

/-- Outside the transition layer the cutoff derivative equals one. -/
theorem deriv_boundaryFactor_eq {a r : ℝ} (h : 2 < a*r) :
    deriv (boundaryFactor a) r = 1 := by
  rw [(hasDerivAt_boundaryFactor a r).deriv,
    smoothTransition.one_of_one_le (by linarith),
    transition_deriv_zero (Or.inr (by linarith)), mul_zero, add_zero]

/-- Equal measures give the same L² function/gradient graph space. -/
def graphMeasureCongr (n : ℕ) (h : μ = ν) :
    (Lp ℝ 2 μ × lp (fun _ : Fin n ↦ Lp ℝ 2 μ) 2) ≃ₗᵢ[ℝ]
    (Lp ℝ 2 ν × lp (fun _ : Fin n ↦ Lp ℝ 2 ν) 2) := by
  subst ν
  exact LinearIsometryEquiv.refl ℝ _

/-- Measure transport preserves the scalar L² inner product separately. -/
theorem graphMeasureCongr_inner_fst (n : ℕ) (h : μ = ν)
    (u v : Lp ℝ 2 μ × lp (fun _ : Fin n ↦ Lp ℝ 2 μ) 2) :
    inner ℝ (graphMeasureCongr n h u).1 (graphMeasureCongr n h v).1 = inner ℝ u.1 v.1 := by
  subst ν
  rfl

/-- Measure transport preserves the gradient Hilbert-sum inner product separately. -/
theorem graphMeasureCongr_inner_snd (n : ℕ) (h : μ = ν)
    (u v : Lp ℝ 2 μ × lp (fun _ : Fin n ↦ Lp ℝ 2 μ) 2) :
    inner ℝ (graphMeasureCongr n h u).2 (graphMeasureCongr n h v).2 = inner ℝ u.2 v.2 := by
  subst ν
  rfl

/-- The algebraic partial derivative equals the Cartesian Fréchet derivative. -/
theorem fderiv_eval_single (p : MvPolynomial (Fin n) ℝ) (i : Fin n) (x : Fin n → ℝ) :
    fderiv ℝ (fun y : Fin n → ℝ ↦ eval y p) x (Pi.single i 1) = eval x (pderiv i p) := by
  induction p using MvPolynomial.induction_on with
  | C c => simp
  | add p q hp hq =>
    simp only [eval_add]
    rw [fderiv_fun_add (p.contDiff_eval.differentiable (by simp)).differentiableAt
      (q.contDiff_eval.differentiable (by simp)).differentiableAt]
    simp [hp, hq]
  | mul_X p j hp =>
    simp only [eval_mul, eval_X]
    rw [fderiv_fun_mul (p.contDiff_eval.differentiable (by simp)).differentiableAt
      (differentiableAt_apply j x)]
    rw [(hasFDerivAt_apply j x).fderiv]
    by_cases h : j = i <;> simp [hp, pderiv_X, h]

/-- Polynomial evaluation satisfies integration by parts against compactly supported smooth tests. -/
theorem integral_eval_mul_testDerivative (p : MvPolynomial (Fin n) ℝ)
    (φ : TestFunction Ω ℝ ⊤) (i : Fin n) :
    (∫ x, eval x p * testDerivative Ω i φ x) =
      -(∫ x, eval x (pderiv i p) * φ x) := by
  have hi (p : MvPolynomial (Fin n) ℝ) (φ : TestFunction Ω ℝ ⊤) :
      Integrable (fun x ↦ eval x p * φ x) :=
    (p.continuous_eval.mul φ.continuous).integrable_of_hasCompactSupport
      φ.hasCompactSupport.mul_left
  have h := integral_mul_fderiv_eq_neg_fderiv_mul_of_integrable
    (μ := volume) (f := fun x ↦ eval x p) (g := fun x ↦ φ x) (v := Pi.single i 1)
    (by simpa only [fderiv_eval_single] using hi (pderiv i p) φ)
    (by simpa only [testDerivative_apply] using hi p (testDerivative Ω i φ))
    (hi p φ)
    (fun x _ ↦ (p.contDiff_eval.differentiable (by simp)).differentiableAt)
    (fun x _ ↦ (φ.contDiff.differentiable (by simp)).differentiableAt)
  simpa only [testDerivative_apply, fderiv_eval_single] using h

/-- Smooth compactly supported function/gradient pairs belong to the closed graph. -/
theorem testGraph_mem (f : TestFunction Ω ℝ ⊤) : testGraph Ω f ∈ graph Ω :=
  Submodule.le_topologicalClosure _ ⟨f, rfl⟩

/-- The zero-boundary Sobolev graph is closed in the product L² topology. -/
theorem graph_isClosed : IsClosed (graph Ω : Set (FunctionL2 Ω × GradientL2 Ω)) :=
  Submodule.isClosed_topologicalClosure _

/-- The L² representative agrees almost everywhere with the test function. -/
theorem coe_testL2 (f : TestFunction Ω ℝ ⊤) :
    (testL2 Ω f : (Fin n → ℝ) → ℝ) =ᵐ[volume.restrict (Ω : Set (Fin n → ℝ))] f :=
  (f.continuous.memLp_of_hasCompactSupport f.hasCompactSupport).coeFn_toLp

/-- Test-function differentiation is Cartesian Fréchet differentiation. -/
theorem testDerivative_apply (f : TestFunction Ω ℝ ⊤) (i : Fin n) (x : Fin n → ℝ) :
    testDerivative Ω i f x = fderiv ℝ f x (Pi.single i 1) := by
  rw [testDerivative, TestFunction.lineDerivCLM_apply_of_le (by simp)]
  exact (f.contDiff.differentiable (by simp)).differentiableAt.lineDeriv_eq_fderiv

/-- The L² pairing of tests is their ordinary Lebesgue integral on coordinate space. -/
theorem inner_testL2 (f g : TestFunction Ω ℝ ⊤) :
    inner ℝ (testL2 Ω f) (testL2 Ω g) = ∫ x, f x*g x := by
  rw [L2.inner_def]
  calc
    _ = ∫ x in (Ω : Set (Fin n → ℝ)), f x*g x := by
      apply integral_congr_ae
      filter_upwards [coe_testL2 Ω f, coe_testL2 Ω g] with x hf hg
      simp [hf, hg, mul_comm]
    _ = _ := setIntegral_eq_integral_of_forall_compl_eq_zero (fun x hx ↦ by
      have hf : f x = 0 := by
        by_contra hn
        exact hx (f.tsupport_subset (subset_tsupport f (Function.mem_support.mpr hn)))
      rw [hf, zero_mul])

/-- Compactly supported smooth tests satisfy integration by parts. -/
theorem inner_testDerivative (f g : TestFunction Ω ℝ ⊤) (i : Fin n) :
    inner ℝ (testL2 Ω f) (testL2 Ω (testDerivative Ω i g)) =
      -inner ℝ (testL2 Ω (testDerivative Ω i f)) (testL2 Ω g) := by
  rw [inner_testL2, inner_testL2]
  have hi (a b : TestFunction Ω ℝ ⊤) : Integrable (fun x ↦ a x*b x) :=
    (a.continuous.mul b.continuous).integrable_of_hasCompactSupport
      a.hasCompactSupport.mul_right
  simp_rw [testDerivative_apply]
  apply integral_mul_fderiv_eq_neg_fderiv_mul_of_integrable
  · simpa only [testDerivative_apply] using hi (testDerivative Ω i f) g
  · simpa only [testDerivative_apply] using hi f (testDerivative Ω i g)
  · exact hi f g
  · exact fun x _ ↦ (f.contDiff.differentiable (by simp)).differentiableAt
  · exact fun x _ ↦ (g.contDiff.differentiable (by simp)).differentiableAt

/-- Every pair in the closed Sobolev graph satisfies the weak-derivative identities. -/
theorem graph_weakDerivative {u : FunctionL2 Ω × GradientL2 Ω} (hu : u ∈ graph Ω)
    (φ : TestFunction Ω ℝ ⊤) (i : Fin n) :
    inner ℝ u.1 (testL2 Ω (testDerivative Ω i φ)) =
      -inner ℝ (u.2 i) (testL2 Ω φ) := by
  have hc : IsClosed {u : FunctionL2 Ω × GradientL2 Ω |
      inner ℝ u.1 (testL2 Ω (testDerivative Ω i φ)) =
        -inner ℝ (u.2 i) (testL2 Ω φ)} := by
    apply isClosed_eq
    · fun_prop
    · exact (((lp.evalCLM ℝ (fun _ : Fin n ↦ FunctionL2 Ω) 2 i).continuous.comp
        continuous_snd).inner continuous_const).neg
  apply closure_minimal (s := Set.range (testGraph Ω)) ?_ hc hu
  rintro _ ⟨f, rfl⟩
  change inner ℝ (testL2 Ω f) (testL2 Ω (testDerivative Ω i φ)) =
    -inner ℝ (testGradient Ω f i) (testL2 Ω φ)
  rw [testGradient_apply]
  exact inner_testDerivative Ω f φ i

/-- Euler differentiation preserves radial angular sectors. -/
theorem diskEuler_diskRadialLift (u : Bool) (p : Polynomial ℝ) (m : ℕ) :
    diskEuler (diskRadialLift u p m) =
      diskRadialLift u (2*Polynomial.X*p.derivative+(m : Polynomial ℝ)*p) m := by
  have h (z : MvPolynomial (Fin 2) ℝ)
      (hz : X 0*pderiv 0 z+X 1*pderiv 1 z = (m : MvPolynomial (Fin 2) ℝ)*z) :
      diskEuler (p.eval₂ C (X 0^2+X 1^2)*z) =
        (2*Polynomial.X*p.derivative+(m : Polynomial ℝ)*p).eval₂ C (X 0^2+X 1^2)*z := by
    change X 0*pderiv 0 (p.eval₂ C (X 0^2+X 1^2)*z) +
      X 1*pderiv 1 (p.eval₂ C (X 0^2+X 1^2)*z) = _
    simp only [pderiv_mul, pderiv_zero_radial, pderiv_one_radial]
    simp only [Polynomial.eval₂_add, Polynomial.eval₂_mul, Polynomial.eval₂_X,
      Polynomial.eval₂_ofNat, Polynomial.eval₂_natCast]
    linear_combination p.eval₂ C (X 0^2+X 1^2)*hz
  cases u <;> simp only [diskRadialLift, Bool.false_eq_true, ↓reduceIte]
  · exact h _ (euler_diskHarmonics m).1
  · exact h _ (euler_diskHarmonics m).2

/-- Only the constant angular sector contributes to the integrated radial divergence. -/
theorem integral_diskRadialDivergence_diskRadialLift {q : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1)
    (u : Bool) (m : ℕ) (p : Polynomial ℝ) :
    ∫ v, eval v (diskRadialDivergence (diskRadialLift u p m)) ∂zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0 =
      if u = false ∧ m = 0 then 2*Real.pi*(p.eval 1-q*p.eval q) else 0 := by
  have hone : diskRadialLift false (1 : Polynomial ℝ) 0 = 1 := by
    simp [diskRadialLift, diskHarmonics]
  have he (f : MvPolynomial (Fin 2) ℝ) :
      (∫ v, eval v f ∂zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) =
        integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) f (diskRadialLift false 1 0) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) := by
    simp [hone, integralInner]
  rw [he]
  simp only [diskRadialDivergence_diskRadialLift]
  rw [integralInner_annulusRadialLift hq (by norm_num) (by norm_num)]
  by_cases h : u = false ∧ m = 0
  · obtain ⟨rfl, rfl⟩ := h
    simp only [and_self, ↓reduceIte, Real.rpow_zero, Nat.cast_zero, pow_zero,
      Polynomial.eval_one, mul_one, one_mul, zernikeAngularNormSq, Bool.false_eq_true,
      zero_mul, add_zero]
    have hd := intervalIntegral.integral_eq_sub_of_hasDerivAt
      (a := q) (b := 1) (fun x _ ↦ (Polynomial.X*p).hasDerivAt x)
      ((Polynomial.X*p).derivative.continuous.intervalIntegrable q 1)
    have hi : (∫ s in q..1, (2*(Polynomial.X*p.derivative+p)).eval s) =
        2*(p.eval 1-q*p.eval q) := by
      calc
        _ = (2 : ℝ) * ∫ s in q..1, (Polynomial.X*p).derivative.eval s := by
          rw [← intervalIntegral.integral_const_mul]
          apply intervalIntegral.integral_congr
          intro s hs
          simp [Polynomial.derivative_mul]
          ring
        _ = _ := by rw [hd]; simp
    rw [intervalIntegral.integral_of_le hq.2.le, ← restrict_Ioo_eq_restrict_Ioc,
      ← integral_subtype_comap measurableSet_Ioo] at hi
    have hh := congrArg (fun a : ℝ ↦ (1/2)*a*(2*Real.pi)) hi
    convert hh using 1
    ring
  · simp [h]

/-- A squared boundary factor annihilates radial divergence flux. -/
theorem integral_diskRadialDivergence_boundary_sq {q : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1)
    (p : MvPolynomial (Fin 2) ℝ) :
    ∫ v, eval v (diskRadialDivergence (annulusBoundaryWeight q ^ 2*p))
      ∂zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0 = 0 := by
  let μ := zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0
  let hμ := hasFiniteMoments_zernikeAnnulusMeasure hq
    (by norm_num : (-1 : ℝ) < 0) (by norm_num : (-1 : ℝ) < 0)
  let L : MvPolynomial (Fin 2) ℝ →ₗ[ℝ] ℝ :=
    ((integralBilinForm μ hμ).flip 1).comp
      (diskRadialDivergence.comp (LinearMap.mulLeft ℝ (annulusBoundaryWeight q ^ 2)))
  have he (z : MvPolynomial (Fin 2) ℝ) : L z =
      ∫ v, eval v (diskRadialDivergence (annulusBoundaryWeight q ^ 2*z)) ∂μ := by
    simp [L, integralBilinForm, integralInner]
  have hk : LinearMap.ker L = ⊤ := eq_top_of_annulusLift_mem hq _ (by
    intro u m z
    rw [LinearMap.mem_ker, he]
    change (∫ v, eval v (diskRadialDivergence
      (annulusBoundaryWeight q ^ 2 * diskRadialLift u _ m)) ∂zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) = 0
    rw [annulusBoundaryWeight_sq_mul_radialLift, integral_diskRadialDivergence_diskRadialLift hq]
    split_ifs <;> simp)
  have hp : p ∈ LinearMap.ker L := by rw [hk]; trivial
  rw [LinearMap.mem_ker, he] at hp
  exact hp

/-- The radial divergence identity gives the polynomial dilation integration-by-parts formula. -/
theorem integralInner_annulus_boundary_euler {q : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1)
    (p : MvPolynomial (Fin 2) ℝ) :
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (annulusBoundaryWeight q*p) (diskEuler (annulusBoundaryWeight q*p)) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) =
      -integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (annulusBoundaryWeight q*p) (annulusBoundaryWeight q*p) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) := by
  let μ := zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0
  let hμ := hasFiniteMoments_zernikeAnnulusMeasure hq
    (by norm_num : (-1 : ℝ) < 0) (by norm_num : (-1 : ℝ) < 0)
  let J : MvPolynomial (Fin 2) ℝ →ₗ[ℝ] ℝ := (integralBilinForm μ hμ).flip 1
  have hJ (f g : MvPolynomial (Fin 2) ℝ) : J (f*g) = integralInner μ (IntegrablePair.of μ f g (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) := by
    simp [J, integralBilinForm, integralInner]
  let f := annulusBoundaryWeight q*p
  have hz : J (diskRadialDivergence (f*f)) = 0 := by
    have he : f*f = annulusBoundaryWeight q^2*(p*p) := by dsimp [f]; ring
    rw [he]
    simpa [J, integralBilinForm, integralInner] using
      integral_diskRadialDivergence_boundary_sq hq (p*p)
  have he : diskRadialDivergence (f*f) = (f*f+f*f)+(f*diskEuler f+f*diskEuler f) := by
    change (2 : ℝ) • (f*f)+diskEuler (f*f) = _
    rw [diskEuler_mul, show (2 : ℝ) = 1+1 by norm_num, add_smul, one_smul]
    ring
  rw [he, map_add, map_add, map_add, hJ, hJ] at hz
  change integralInner μ (IntegrablePair.of μ f (diskEuler f) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) = -integralInner μ (IntegrablePair.of μ f f (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _))
  linarith

/-- A degree- and mode-independent Poincaré estimate for boundary-weighted annular polynomials. -/
theorem integralInner_annulus_boundary_le_gradient {q : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1)
    (p : MvPolynomial (Fin 2) ℝ) :
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (annulusBoundaryWeight q*p) (annulusBoundaryWeight q*p) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) ≤
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (pderiv 0 (annulusBoundaryWeight q*p)) (pderiv 0 (annulusBoundaryWeight q*p)) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) +
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (pderiv 1 (annulusBoundaryWeight q*p)) (pderiv 1 (annulusBoundaryWeight q*p)) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) := by
  let μ := zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0
  let hμ := hasFiniteMoments_zernikeAnnulusMeasure hq
    (by norm_num : (-1 : ℝ) < 0) (by norm_num : (-1 : ℝ) < 0)
  let J : MvPolynomial (Fin 2) ℝ →ₗ[ℝ] ℝ := (integralBilinForm μ hμ).flip 1
  have hJ (f g : MvPolynomial (Fin 2) ℝ) : J (f*g) = integralInner μ (IntegrablePair.of μ f g (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) := by
    simp [J, integralBilinForm, integralInner]
  let f := annulusBoundaryWeight q*p
  let P := diskGradientProduct f f+f*f+(f*diskEuler f+f*diskEuler f)
  have hn : 0 ≤ J P := by
    change 0 ≤ ∫ v, eval v P * eval v 1 ∂μ
    simp only [map_one, mul_one]
    apply integral_nonneg_zernikeAnnulusMeasure hq 0 0 P
      (_hi := by simpa only [map_one, mul_one] using hμ.integrable_mul P 1)
    intro v hv
    dsimp [P]
    simp only [diskGradientProduct, map_add, map_mul]
    change 0 ≤ eval v (pderiv 0 f)*eval v (pderiv 0 f) +
      eval v (pderiv 1 f)*eval v (pderiv 1 f) + eval v f*eval v f +
      (eval v f*eval v (diskEuler f)+eval v f*eval v (diskEuler f))
    have he : eval v (diskEuler f) = v 0*eval v (pderiv 0 f)+v 1*eval v (pderiv 1 f) := by
      change eval v (X 0*pderiv 0 f+X 1*pderiv 1 f) = _
      simp
    rw [he]
    have hid (a b c x y : ℝ) :
        b*b+c*c+a*a+(a*(x*b+y*c)+a*(x*b+y*c)) =
        (a+x*b+y*c)^2+(y*b-x*c)^2+(1-(x^2+y^2))*(b^2+c^2) := by ring
    rw [hid]
    exact add_nonneg (add_nonneg (sq_nonneg _) (sq_nonneg _))
      (mul_nonneg (sub_nonneg.mpr hv) (add_nonneg (sq_nonneg _) (sq_nonneg _)))
  dsimp [P] at hn
  rw [map_add, map_add, map_add, diskGradientProduct, map_add, hJ, hJ, hJ, hJ] at hn
  have he := integralInner_annulus_boundary_euler hq p
  change integralInner μ (IntegrablePair.of μ f (diskEuler f) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) = -integralInner μ (IntegrablePair.of μ f f (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) at he
  change integralInner μ (IntegrablePair.of μ f f (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) ≤
    integralInner μ (IntegrablePair.of μ (pderiv 0 f) (pderiv 0 f) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _))+integralInner μ (IntegrablePair.of μ (pderiv 1 f) (pderiv 1 f) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _))
  linarith

/-- Equivalently, physical polynomial mass is bounded by the negative-Laplacian energy. -/
theorem integralInner_annulus_boundary_le_neg_laplacian {q : ℝ} (hq : q ∈ Ioo (0 : ℝ) 1)
    (p : MvPolynomial (Fin 2) ℝ) :
    integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (annulusBoundaryWeight q*p) (annulusBoundaryWeight q*p) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) ≤
    -integralInner (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (IntegrablePair.of (zernikeAnnulusMeasure ⟨q,hq.1.le,hq.2⟩ 0 0) (diskLaplacian (annulusBoundaryWeight q*p)) (annulusBoundaryWeight q*p) (h := (hasFiniteMoments_zernikeAnnulusMeasure hq (by norm_num) (by norm_num)).integrable_mul _ _)) := by
  rw [← integralInner_annulus_gradient_eq_neg_laplacian hq]
  exact integralInner_annulus_boundary_le_gradient hq p

/-- Polynomial evaluation as a linear map to weighted L² (Dunkl–Xu, Chapter 3). -/
noncomputable def polynomialToL2 : MvPolynomial σ ℝ →ₗ[ℝ] Lp ℝ 2 μ where
  toFun p := (memLp_eval_of_hasFiniteMoments hμ p).toLp (fun x ↦ eval x p)
  map_add' p q := by
    apply Lp.ext
    filter_upwards [(memLp_eval_of_hasFiniteMoments hμ (p + q)).coeFn_toLp,
      (memLp_eval_of_hasFiniteMoments hμ p).coeFn_toLp, (memLp_eval_of_hasFiniteMoments hμ q).coeFn_toLp,
      Lp.coeFn_add ((memLp_eval_of_hasFiniteMoments hμ p).toLp _) ((memLp_eval_of_hasFiniteMoments hμ q).toLp _)]
      with x hpq hp hq hadd
    rw [hpq, hadd]
    simp only [Pi.add_apply, hp, hq, eval_add]
  map_smul' c p := by
    change (memLp_eval_of_hasFiniteMoments hμ (c • p)).toLp _ = c • (memLp_eval_of_hasFiniteMoments hμ p).toLp _
    apply Lp.ext
    filter_upwards [(memLp_eval_of_hasFiniteMoments hμ (c • p)).coeFn_toLp,
      (memLp_eval_of_hasFiniteMoments hμ p).coeFn_toLp,
      Lp.coeFn_smul c ((memLp_eval_of_hasFiniteMoments hμ p).toLp _)] with x hcp hp hsmul
    rw [hcp, hsmul]
    simp only [Pi.smul_apply, hp, smul_eq_C_mul, map_mul, eval_C, smul_eq_mul]

/-- An almost-everywhere representative of polynomial evaluation (Dunkl–Xu, Chapter 3). -/
theorem coe_polynomialToL2 (p : MvPolynomial σ ℝ) : polynomialToL2 hμ p =ᵐ[μ] fun x ↦ eval x p :=
  (memLp_eval_of_hasFiniteMoments hμ p).coeFn_toLp

/-- The polynomial L² inner product is the weighted product integral (Dunkl–Xu, Chapter 3). -/
theorem inner_polynomialToL2 (p q : MvPolynomial σ ℝ) :
    ⟪polynomialToL2 hμ p, polynomialToL2 hμ q⟫_ℝ = ∫ x, eval x p * eval x q ∂μ := by
  rw [L2.inner_def]
  apply integral_congr_ae
  filter_upwards [coe_polynomialToL2 hμ p, coe_polynomialToL2 hμ q] with x hp hq
  simp [hp, hq, mul_comm]

/-- Apply the contractions coordinatewise, retaining their range subspaces. -/
def rangeCoordinates : lp G 2 →L[ℝ] lp (fun i ↦ LinearMap.range (T i).toLinearMap) 2 :=
  LinearMap.mkContinuous
    { toFun := fun x ↦ ⟨fun i ↦ ⟨T i (x i), LinearMap.mem_range_self _ _⟩,
        (lp.memℓp x).mono' (fun i ↦ hT i (x i))⟩
      map_add' := by
        intro x y
        apply lp.ext
        funext i
        apply Subtype.ext
        exact map_add (T i) _ _
      map_smul' := by
        intro c x
        apply lp.ext
        funext i
        apply Subtype.ext
        exact map_smul (T i) c _ }
    1 (by
      intro x
      simpa only [one_mul] using
        (lp.norm_mono (by norm_num : (2 : ENNReal) ≠ 0) (fun i ↦ hT i (x i))))

/-- Orthogonal contractions synthesize to a continuous map from the full Hilbert sum. -/
def orthogonalSynthesis : lp G 2 →L[ℝ] E :=
  (orthogonalFamily_range T horth).linearIsometry.toContinuousLinearMap.comp (rangeCoordinates T hT)

/-- Synthesis is the unconditionally convergent sum of the coordinate images. -/
theorem hasSum_orthogonalSynthesis (x : lp G 2) :
    HasSum (fun i ↦ T i (x i)) (orthogonalSynthesis T hT horth x) :=
  (orthogonalFamily_range T horth).hasSum_linearIsometry (rangeCoordinates T hT x)

/-- Orthogonal synthesis does not increase the Hilbert-sum norm. -/
theorem norm_orthogonalSynthesis_le (x : lp G 2) :
    ‖orthogonalSynthesis T hT horth x‖ ≤ ‖x‖ := by
  change ‖(orthogonalFamily_range T horth).linearIsometry (rangeCoordinates T hT x)‖ ≤ _
  rw [LinearIsometry.norm_map]
  exact lp.norm_mono (by norm_num) (fun i ↦ hT i (x i))

/-- A single nonzero coordinate synthesizes to precisely its operator image. -/
theorem orthogonalSynthesis_single [DecidableEq ι] (i : ι) (x : G i) :
    orthogonalSynthesis T hT horth (lp.single 2 i x) = T i x := by
  rw [orthogonalSynthesis_apply, tsum_eq_single i]
  · simp [lp.single_apply]
  · intro j hj
    simp [lp.single_apply, hj]

/-- Synthesis is faithful when each coordinate operator is faithful. -/
theorem orthogonalSynthesis_injective (hi : ∀ i, Function.Injective (T i)) :
    Function.Injective (orthogonalSynthesis T hT horth) := by
  intro x y h
  have he := (orthogonalFamily_range T horth).linearIsometry.injective h
  apply lp.ext
  funext i
  apply hi i
  exact congrArg (fun z : lp (fun i ↦ LinearMap.range (T i).toLinearMap) 2 ↦ (z i : E)) he

/-- Trial vectors have finitely many sectors and a finite trial vector in each sector. -/
def trialSpace (c : Finset ι × ℕ) : Submodule ℝ (lp E 2) :=
  ⨆ i : c.1, (W i c.2).map (lp.lsingle (𝕜 := ℝ) (E := E) 2 (i : ι))

/-- Independent finite sector projections reconstruct the combined orthogonal projection. -/
theorem sum_projection_eq (c : Finset ι × ℕ) (x : lp E 2) :
    ∑ i ∈ c.1, lp.single 2 i ((W i c.2).starProjection (x i)) =
      (trialSpace W c).starProjection x := by
  symm
  apply Submodule.eq_starProjection_of_mem_of_inner_eq_zero
  · exact Submodule.sum_mem _ fun i hi ↦
      single_mem_trialSpace W c hi ((W i c.2).starProjection_apply_mem _)
  · intro z hz
    unfold trialSpace at hz
    refine Submodule.iSup_induction _ (motive := fun z ↦
      inner ℝ (x - ∑ j ∈ c.1, lp.single 2 j ((W j c.2).starProjection (x j))) z = 0)
      hz ?_ ?_ ?_
    · intro i z hz
      obtain ⟨y, hy, rfl⟩ := hz
      change inner ℝ (x - ∑ j ∈ c.1, lp.single 2 j ((W j c.2).starProjection (x j)))
        (lp.single 2 i y) = 0
      rw [lp.inner_single_right]
      simp only [lp.coeFn_sub, Pi.sub_apply, lp.coeFn_sum, Finset.sum_apply,
        lp.coeFn_single, Finset.sum_pi_single, ite_eq_left i.property]
      exact (W i c.2).starProjection_inner_eq_zero _ _ hy
    · simp
    · intro a b ha hb
      simp [inner_add_right, ha, hb]

/-- Dense radial trial spaces give a dense union of combined finite trial spaces. -/
theorem trialSpace_closure_iSup (hW : ∀ i, Monotone (W i))
    (hd : ∀ i, (⨆ N, W i N).topologicalClosure = ⊤) :
    (⨆ c, trialSpace W c).topologicalClosure = ⊤ := by
  let K := (⨆ c, trialSpace W c).topologicalClosure
  have hK : IsClosed (K : Set (lp E 2)) := Submodule.isClosed_topologicalClosure _
  have hs (i : ι) (y : E i) : lp.single 2 i y ∈ K := by
    have ht := (lp.singleContinuousLinearMap ℝ E 2 i).continuous.tendsto y |>.comp
      (Submodule.starProjection_tendsto_self (W i) (hW i) y (by rw [hd i]))
    apply hK.mem_of_tendsto ht
    apply Eventually.of_forall
    intro N
    apply Submodule.le_topologicalClosure
    apply (le_iSup (trialSpace W) ({i}, N))
    exact single_mem_trialSpace W ({i}, N) (Finset.mem_singleton_self i)
      ((W i N).starProjection_apply_mem _)
  apply le_antisymm le_top
  intro x hx
  apply hK.mem_of_tendsto (lp.hasSum_single (by norm_num) x)
  exact Eventually.of_forall fun S ↦ Submodule.sum_mem _ fun i hi ↦ hs i (x i)

/-- Actual sector projections converge jointly as the finite sector set and radial degree grow. -/
theorem tendsto_sum_projection (hW : ∀ i, Monotone (W i))
    (hd : ∀ i, (⨆ N, W i N).topologicalClosure = ⊤) (x : lp E 2) :
    Tendsto (fun c : Finset ι × ℕ ↦
      ∑ i ∈ c.1, lp.single 2 i ((W i c.2).starProjection (x i))) atTop (𝓝 x) := by
  simp only [sum_projection_eq]
  exact Submodule.starProjection_tendsto_self (trialSpace W) (trialSpace_mono W hW) x
    (by rw [trialSpace_closure_iSup W hW hd])

/-- A positive symmetric polynomial operator has positive definite Galerkin matrices. -/
theorem posDef_integralInner_operator (hμ : HasFiniteMoments μ)
    (L : ℝ[X] →ₗ[ℝ] ℝ[X]) (f : ι → ℝ[X]) (hf : LinearIndependent ℝ f)
    (hsym : ∀ p q, integralInner μ (IntegrablePair.of μ (L p) q) = integralInner μ (IntegrablePair.of μ (L q) p))
    (hpos : ∀ p, p≠0 → 0 < integralInner μ (IntegrablePair.of μ (L p) p)) :
    Matrix.PosDef (Matrix.of (fun i j ↦ integralInner μ (IntegrablePair.of μ (L (f j)) (f i)))) := by
  rw [Matrix.posDef_iff_dotProduct_mulVec]
  constructor
  · rw [Matrix.isHermitian_iff_isSymm, Matrix.IsSymm.ext_iff]
    intro i j
    exact hsym _ _
  · intro v hv
    have hp : (∑ i, v i • f i)≠0 := by
      intro h
      apply hv
      funext i
      exact (Fintype.linearIndependent_iff.mp hf v h) i
    simpa only [star_trivial, ← integralInner_operator_sum hμ] using hpos _ hp

/-- Finite polynomial synthesis turns the matrix rows into tested residuals. -/
theorem integralInner_operator_sum_apply (hμ : HasFiniteMoments μ)
    (L : ℝ[X] →ₗ[ℝ] ℝ[X]) (f : ι → ℝ[X]) (v : ι → ℝ) (i : ι) :
    integralInner μ (IntegrablePair.of μ (L (∑ j, v j • f j)) (f i)) =
      ((Matrix.of (fun i j ↦ integralInner μ (IntegrablePair.of μ (L (f j)) (f i)))) *ᵥ v) i := by
  have hsmul (p q : ℝ[X]) (a : ℝ) :=
    integralInner_smul_left μ p q a (h := hμ.integrable_mul p q)
  simp only [map_sum, map_smul, integralInner_sum_left hμ, hsmul,
    mulVec, dotProduct, Matrix.of_apply]
  apply Finset.sum_congr rfl
  intro j hj
  ring