/-
Copyright (c) 2026 Michael R. Douglas. All rights reserved.
Released under Apache 2.0 license as described in the file LICENSE.

# Multivariate Wick Polynomials in the Eigenbasis

Identifies multivariate Wick polynomials in the lattice GFF with
multivariate Hermite polynomials in the eigenbasis-orthogonalized
coordinates `(ξ_k)` from `GaussianField/StandardGaussianBridge.lean`.

The setup: the lattice GFF measure pushes forward via
`gffOrthonormalProj` to the standard multivariate Gaussian
`Π_k gaussianReal 0 1`. Under this change of variables, multivariate
Wick polynomials become multivariate Hermite polynomials, since:

* Each `ξ_k` is standard `N(0,1)` (variance 1).
* The 1D Wick monomial of variance 1 equals the probabilist's Hermite
  polynomial: `:x^n:_1 = He_n(x)` (consequence of `wick_eq_hermiteR`
  in `SchwartzNuclear/HermiteWick.lean`).
* Independence of distinct `ξ_k` makes
  `:ξ^α:_1 := ∏_k :ξ_k^{α_k}:_1` a multivariate Wick monomial that is
  also a multivariate Hermite polynomial.

The site Wick monomial `:φ(x)^d:_{c_a}` then admits an explicit
expansion in eigenbasis multi-indices, giving the chaos decomposition
of the GFF interaction.

## Main definitions

- `gffMultiWickMonomial` — multivariate Wick monomial indexed by a
  multi-index over eigenbasis sites.
- `gffMultiHermiteValue` — multivariate Hermite polynomial in the
  orthogonalized variables (purely as a function of `ω`).

## Main theorems

- `gffMultiWickMonomial_eq_hermiteMulti` — the two coincide.
- `siteWickMonomial_eigenbasis_expansion` — site Wick monomial as an
  explicit linear combination of `gffMultiWickMonomial`s.
- `interaction_centered_in_chaosLE` — the GFF interaction
  $V_a - \mathbb E V_a$ is a finite linear combination of Wick
  monomials of total degree $\le \deg P$, hence (after the change of
  variables) lies in `wienerChaosLE _ (deg P)`.
- `gff_wickPower_two_site_inner` — two-**site** Wick power formula:
  $\int :\varphi(x)^n: \, :\varphi(y)^m: \, d\mu_{\mathrm{GFF}}
  = \delta_{nm}\, n!\, C(x,y)^n$ (Janson–Hilbert orthogonality).
- `gff_wickPower_two_smeared_inner` — its **smeared** generalization to
  arbitrary test fields $f, g$:
  $\int :\varphi(f)^n: \, :\varphi(g)^m: \, d\mu_{\mathrm{GFF}}
  = \delta_{nm}\, n!\, \langle\varphi(f)\varphi(g)\rangle^n$
  (the Wick/Mehler kernel; $n=m=4$ gives $4!\,\langle\varphi(f)\varphi(g)\rangle^4$).
  Proved by $\omega$-linearity from the multinomial Wick expansion
  `wickMonomial_pow_sum_expansion_of_totalDegree`, reusing the site proof
  verbatim with $\gamma_j(x)\mapsto\Gamma_j(f)=\sum_x f(x)\gamma_j(x)$.

## References

- S. Janson, *Gaussian Hilbert Spaces*, Cambridge (1997), §3.1
  (Hermite polynomials) and Theorem 3.21 (orthogonality).
- Glimm and Jaffe, *Quantum Physics*, §6.1 (Wick ordering).
- `SchwartzNuclear/HermiteWick.lean` (this repo) — `wick_eq_hermiteR`,
  the 1D building block.

## Status

The identification + orthogonality theorems are **proved** (no axioms);
the two-site and smeared two-field Wick power formulas
(`gff_wickPower_two_site_inner`, `gff_wickPower_two_smeared_inner`,
2026-06-05) are sorry-free and axiom-clean. QFT-specific
specializations (e.g. interaction polynomial $V_a$ in the chaos)
live downstream in pphi2's `PolynomialChaosBridge.lean`; this file
provides only the generic Gaussian-Hilbert-space identifications.
-/

import GaussianField.StandardGaussianBridge
import GaussianField.Wick
import SchwartzNuclear.HermiteWick
import SchwartzNuclear.HermiteFunctions

noncomputable section

namespace GaussianField

open MeasureTheory ProbabilityTheory

variable (d N : ℕ) [NeZero N]

/-- Total degree of a multi-index over `FinLatticeSites d N`. -/
def MultiIndexLattice.totalDegree {d N : ℕ} [NeZero N]
    (α : FinLatticeSites d N → ℕ) : ℕ :=
  haveI : Fintype (ZMod N) := ZMod.fintype N
  haveI : Fintype (FinLatticeSites d N) := Pi.instFintype
  ∑ k, α k

/-- All multi-indices on `FinLatticeSites d N` with each component
bounded by `D + 1` (so total degree is bounded by `D · |Λ|`).
Constrained per-coordinate because `FinLatticeSites d N → ℕ` is not a
`Fintype`; the constraint gives a finite Finset enumeration via
`Fintype.piFinset`. -/
noncomputable def multiIndicesUpToDegree (d N D : ℕ) [NeZero N] :
    Finset (FinLatticeSites d N → ℕ) :=
  haveI : Fintype (ZMod N) := ZMod.fintype N
  haveI : Fintype (FinLatticeSites d N) := Pi.instFintype
  Fintype.piFinset (fun _ : FinLatticeSites d N => Finset.range (D + 1))

/-- Multi-indices on `FinLatticeSites d N` with total degree exactly `k`
(and each coord in `Finset.range (k + 1)`). -/
noncomputable def multiIndicesOfDegree (d N k : ℕ) [NeZero N] :
    Finset (FinLatticeSites d N → ℕ) :=
  (multiIndicesUpToDegree d N k).filter
    (fun α => MultiIndexLattice.totalDegree α = k)

/-- Multivariate Wick monomial in the orthogonalized GFF coordinates
$(\xi_k)$. For a multi-index $\alpha : \mathrm{FinLatticeSites}\ d\ N
\to \mathbb N$, this is the product of 1D unit-variance Wick monomials:

  `:ξ^α:_1 := ∏_k :ξ_k^{α_k}:_1`.

(All variances are 1 because each `ξ_k` is standard Gaussian under the
lattice GFF pushforward.) -/
noncomputable def gffMultiWickMonomial
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (α : FinLatticeSites d N → ℕ)
    (ω : Configuration (FinLatticeField d N)) : ℝ :=
  ∏ k : FinLatticeSites d N,
    wickMonomial (α k) 1 (gffOrthonormalCoord d N a mass ha hmass k ω)

/-- **GFF multivariate Wick monomial = multivariate Hermite in
orthogonalized variables.**

Specifically, since each `ξ_k` is standard `N(0,1)` (variance 1) under
the GFF pushforward, the 1D Wick monomial `:ξ_k^{α_k}:_1` equals the
probabilist's Hermite polynomial $\mathrm{He}_{\alpha_k}(\xi_k)$ by
`wick_eq_hermiteR` (with $c = 1$, $\sqrt c = 1$). Taking the product
over $k$:

  `gffMultiWickMonomial α ω = ∏_k He_{α_k}(ξ_k(ω))`.

This is exactly the multivariate Hermite polynomial of multi-index
$\alpha$ evaluated at the orthogonalized vector. -/
theorem gffMultiWickMonomial_eq_hermite_product
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (α : FinLatticeSites d N → ℕ)
    (ω : Configuration (FinLatticeField d N)) :
    gffMultiWickMonomial d N a mass ha hmass α ω =
      ∏ k : FinLatticeSites d N,
        ((Polynomial.hermite (α k)).map (Int.castRingHom ℝ)).eval
          (gffOrthonormalCoord d N a mass ha hmass k ω) := by
  unfold gffMultiWickMonomial
  refine Finset.prod_congr rfl ?_
  intro k _
  rw [wick_eq_hermiteR (α k) 1 (by norm_num : (0:ℝ) < 1)]
  unfold scaledHermite
  rw [Real.sqrt_one, one_pow, one_mul, div_one]

/-- 1D **Wick monomial orthogonality** under the standard Gaussian
`gaussianReal 0 1`:

  `∫ :x^m:_1 · :x^n:_1 ∂γ = δ_{m,n} · m!`.

Reduced to `J_eq` (Hermite L²-inner product against the unnormalized
Gaussian weight) via `wick_eq_hermiteR` (with `c = 1`, `√c = 1`,
`x/√c = x`), then converted from Lebesgue to `gaussianReal 0 1` via
the explicit PDF `(√(2π))⁻¹ · exp(-x²/2)`. -/
theorem wickMonomial_inner_gaussianReal_one (m n : ℕ) :
    ∫ x : ℝ, wickMonomial m 1 x * wickMonomial n 1 x ∂(gaussianReal 0 1) =
    if m = n then (m.factorial : ℝ) else 0 := by
  -- Step 1: wickMonomial k 1 x = (hermiteR k).eval x.
  have h_wick : ∀ (k : ℕ) (x : ℝ),
      wickMonomial k 1 x = (hermiteR k).eval x := by
    intro k x
    rw [wick_eq_hermiteR k 1 (by norm_num : (0:ℝ) < 1)]
    show Real.sqrt 1 ^ k * (hermiteR k).eval (x / Real.sqrt 1) = _
    rw [Real.sqrt_one, one_pow, one_mul, div_one]
  simp_rw [h_wick]
  -- Step 2: Convert gaussianReal 0 1 integral to Lebesgue with the PDF.
  rw [integral_gaussianReal_eq_integral_smul (one_ne_zero)]
  simp_rw [smul_eq_mul]
  -- Step 3: Identify gaussianPDFReal 0 1 x = (√(2π))⁻¹ * gaussian x.
  have h_pdf : ∀ x : ℝ,
      gaussianPDFReal 0 1 x = (Real.sqrt (2 * Real.pi))⁻¹ * gaussian x := by
    intro x
    rw [gaussianPDFReal_def]
    simp only [NNReal.coe_one, mul_one, sub_zero]
    show (Real.sqrt (2 * Real.pi))⁻¹ * Real.exp (-x^2 / 2) =
        (Real.sqrt (2 * Real.pi))⁻¹ * Real.exp (-(x^2 / 2))
    congr 2
    ring
  simp_rw [h_pdf]
  -- Step 4: Pull out the (√(2π))⁻¹ constant; what's left is `J m n`.
  rw [show (fun x : ℝ =>
        (Real.sqrt (2 * Real.pi))⁻¹ * gaussian x *
          ((hermiteR m).eval x * (hermiteR n).eval x)) =
      (fun x : ℝ => (Real.sqrt (2 * Real.pi))⁻¹ *
        ((hermiteR m).eval x * ((hermiteR n).eval x * gaussian x))) from by
    funext x; ring]
  rw [integral_const_mul]
  show (Real.sqrt (2 * Real.pi))⁻¹ * J m n = _
  rw [J_eq n m]
  -- Goal: (√(2π))⁻¹ * (if m = n then m! * √(2π) else 0) = if m = n then m! else 0
  split_ifs with h
  · have h2pi_pos : (0 : ℝ) < Real.sqrt (2 * Real.pi) :=
      Real.sqrt_pos.mpr (by positivity)
    field_simp
  · ring

/-- **Multivariate Wick orthogonality on the standard product Gaussian.**

For any finite index set `ι` and multi-indices `α, β : ι → ℕ`,
  `∫ ∏_j :ξ_j^{α_j}:_1 · ∏_j :ξ_j^{β_j}:_1 ∂(Π_j gaussianReal 0 1)
     = δ_{α, β} · ∏_j α_j!`.

The product-of-Wicks form is the building block for both
`gffMultiWickMonomial_orthogonality` (after the GFF pushforward) and
the abstract Janson 2-site Wick power formula
(`janson_two_site_wick_power_inner` below) used in the canonical
setup. Direct combination of `Fubini` (`integral_fintype_prod_eq_prod`)
with the 1D orthogonality `wickMonomial_inner_gaussianReal_one`. -/
theorem multiWickMonomial_pi_gaussianReal_inner
    {ι : Type*} [Fintype ι] [DecidableEq ι]
    (α β : ι → ℕ) :
    ∫ ξ : ι → ℝ,
      (∏ j, wickMonomial (α j) 1 (ξ j)) *
      (∏ j, wickMonomial (β j) 1 (ξ j))
      ∂(MeasureTheory.Measure.pi (fun _ : ι => gaussianReal 0 1)) =
    if α = β then ((∏ j, ((α j).factorial : ℕ) : ℕ) : ℝ) else 0 := by
  -- Combine the two products into one product of pairs
  have h_eq : (fun ξ : ι → ℝ =>
      (∏ j, wickMonomial (α j) 1 (ξ j)) * (∏ j, wickMonomial (β j) 1 (ξ j))) =
      (fun ξ : ι → ℝ =>
        ∏ j, wickMonomial (α j) 1 (ξ j) * wickMonomial (β j) 1 (ξ j)) := by
    funext ξ
    rw [← Finset.prod_mul_distrib]
  rw [h_eq]
  -- Apply Fubini for product measures
  rw [integral_fintype_prod_eq_prod
    (f := fun j (x : ℝ) => wickMonomial (α j) 1 x * wickMonomial (β j) 1 x)]
  -- Apply 1D orthogonality to each factor
  simp_rw [wickMonomial_inner_gaussianReal_one]
  -- Collapse per-coordinate indicators into the multi-index indicator
  by_cases hαβ : α = β
  · rw [if_pos hαβ]
    rw [show (∏ j : ι,
        if α j = β j then (((α j).factorial : ℕ) : ℝ) else 0) =
        ∏ j : ι, ((α j).factorial : ℝ) from by
      refine Finset.prod_congr rfl ?_
      intro j _
      rw [if_pos (by rw [hαβ] : α j = β j)]]
    push_cast
    rfl
  · rw [if_neg hαβ]
    obtain ⟨j, hj⟩ : ∃ j, α j ≠ β j := by
      by_contra h
      push Not at h
      exact hαβ (funext h)
    apply Finset.prod_eq_zero (Finset.mem_univ j)
    rw [if_neg hj]

/-- **Orthogonality of GFF multivariate Wick monomials under the lattice
GFF measure.**

  `∫ :ξ^α: · :ξ^β: dμ_GFF = δ_{αβ} · ∏_k α_k!`

This is the multivariate Wick orthogonality: distinct multi-indices
give zero, and the diagonal is a product of factorials.

**Proof:** Push forward to `Π_k gaussianReal 0 1` via
`gffOrthonormalProj_pushforward_eq_stdGaussian`, decompose into a
product of 1D integrals via `integral_fintype_prod_eq_prod`, and apply
the 1D orthogonality `wickMonomial_inner_gaussianReal_one` to each
factor. -/
theorem gffMultiWickMonomial_orthogonality
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (α β : FinLatticeSites d N → ℕ) :
    ∫ ω, gffMultiWickMonomial d N a mass ha hmass α ω *
        gffMultiWickMonomial d N a mass ha hmass β ω
        ∂(latticeGaussianMeasure d N a mass ha hmass) =
    (if α = β then ((∏ k, (α k).factorial : ℕ) : ℝ) else 0) := by
  have h_meas := gffOrthonormalProj_measurable d N a mass ha hmass
  -- Helper: each `wickMonomial k 1` is continuous (it's `(hermiteR k).eval`).
  have h_wick_cont : ∀ k : ℕ, Continuous (wickMonomial k 1) := by
    intro k
    have h_eq : (fun x : ℝ => wickMonomial k 1 x) = (hermiteR k).eval := by
      funext x
      rw [wick_eq_hermiteR k 1 (by norm_num : (0:ℝ) < 1)]
      show Real.sqrt 1 ^ k * (hermiteR k).eval (x / Real.sqrt 1) = _
      rw [Real.sqrt_one, one_pow, one_mul, div_one]
    rw [show wickMonomial k 1 = fun x : ℝ => wickMonomial k 1 x from rfl, h_eq]
    exact (hermiteR k).continuous
  -- Step 1: combine the two products into one product of pairs (using the
  -- definitional equality `gffOrthonormalProj ω k = gffOrthonormalCoord k ω`).
  have h_eq : ∀ ω : Configuration (FinLatticeField d N),
      gffMultiWickMonomial d N a mass ha hmass α ω *
        gffMultiWickMonomial d N a mass ha hmass β ω =
      ∏ k, wickMonomial (α k) 1 (gffOrthonormalProj d N a mass ha hmass ω k) *
            wickMonomial (β k) 1 (gffOrthonormalProj d N a mass ha hmass ω k) := by
    intro ω
    unfold gffMultiWickMonomial gffOrthonormalProj
    rw [← Finset.prod_mul_distrib]
  simp_rw [h_eq]
  -- Step 2: Use integral_map and the pushforward equation.
  have h_strong_meas : AEStronglyMeasurable
      (fun y : FinLatticeSites d N → ℝ =>
        ∏ k, wickMonomial (α k) 1 (y k) * wickMonomial (β k) 1 (y k))
      ((latticeGaussianMeasure d N a mass ha hmass).map
        (gffOrthonormalProj d N a mass ha hmass)) := by
    apply Continuous.aestronglyMeasurable
    apply continuous_finsetProd
    intro k _
    exact ((h_wick_cont _).comp (continuous_apply k)).mul
      ((h_wick_cont _).comp (continuous_apply k))
  rw [← integral_map h_meas.aemeasurable
    (f := fun y : FinLatticeSites d N → ℝ =>
      ∏ k, wickMonomial (α k) 1 (y k) * wickMonomial (β k) 1 (y k))
    h_strong_meas]
  rw [gffOrthonormalProj_pushforward_eq_stdGaussian d N a mass ha hmass]
  -- Step 3: apply Fubini to the product over k.
  rw [integral_fintype_prod_eq_prod
    (f := fun k (x : ℝ) => wickMonomial (α k) 1 x * wickMonomial (β k) 1 x)]
  -- Step 4: each factor is the 1D orthogonality.
  simp_rw [wickMonomial_inner_gaussianReal_one]
  -- Step 4: combine the per-coordinate indicators into a single multi-index indicator.
  by_cases hαβ : α = β
  · rw [if_pos hαβ]
    rw [show (∏ k : FinLatticeSites d N,
        if α k = β k then (((α k).factorial : ℕ) : ℝ) else 0) =
        ∏ k : FinLatticeSites d N, ((α k).factorial : ℝ) from by
      refine Finset.prod_congr rfl ?_
      intro k _
      rw [if_pos (by rw [hαβ] : α k = β k)]]
    push_cast
    rfl
  · rw [if_neg hαβ]
    -- Some k has α k ≠ β k, so that factor is 0.
    obtain ⟨k, hk⟩ : ∃ k, α k ≠ β k := by
      by_contra h
      push Not at h
      exact hαβ (funext h)
    apply Finset.prod_eq_zero (Finset.mem_univ k)
    rw [if_neg hk]

/-- The GFF site (auto-)variance at site `x`:
`c_a(x) = ⟨φ(x), φ(x)⟩ = (a^d)⁻¹ · Σ_k λ_k⁻¹ · e_k(x)²`.

This is the diagonal of the lattice GFF covariance, i.e. the variance
that the local Wick subtraction must subtract for `:φ(x)^k:_{c_a}` to
be the QFT-correct Wick monomial. The site-Wick-expansion axiom below
is true only for this specific value of `c`; for any other `c`, the
expansion picks up lower-degree terms. -/
noncomputable def gffSiteVariance
    (a mass : ℝ) (_ha : 0 < a) (_hmass : 0 < mass)
    (x : FinLatticeSites d N) : ℝ :=
  haveI : Fintype (ZMod N) := ZMod.fintype N
  haveI : Fintype (FinLatticeSites d N) := Pi.instFintype
  (a ^ d)⁻¹ *
    ∑ k : FinLatticeSites d N,
      (massEigenvalues d N a mass k)⁻¹ *
        ((massEigenvectorBasis d N a mass k : EuclideanSpace ℝ _) x) ^ 2

/-! ## Spectral expansion helper for the site-Wick expansion

The orthogonalized eigenvector coefficient `γ_j(x) = e_j(x) / √(a^d λ_j)`
re-expresses `gffSiteVariance` as `∑ γ_j(x)²` (Parseval-style). This
identity decouples `gffSiteVariance` from the explicit `(a^d)⁻¹` and
`λ_j⁻¹` factors. -/

/-- The orthogonalized eigenvector coefficient: `γ_j(x) = e_j(x) / √(a^d λ_j)`. -/
private noncomputable def gffEigenCoeff
    (a mass : ℝ) (j x : FinLatticeSites d N) : ℝ :=
  (massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) x /
    Real.sqrt (a^d * massEigenvalues d N a mass j)

/-- Pointwise completeness of the eigenbasis:
`∑ j, e_j(x) · e_j(y) = δ_{xy}`. The (x,y) entry of `M Mᵀ = I` for
the orthonormal basis matrix `M`. -/
lemma eigenbasis_completeness
    (a mass : ℝ) (x y : FinLatticeSites d N) :
    ∑ j, (massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) x *
        (massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) y =
    if y = x then (1 : ℝ) else 0 := by
  -- sum_repr applied to (EuclideanSpace.single x 1):
  -- ∑ j, b.repr (single x 1) j • b j = single x 1.
  have h_repr := (massEigenvectorBasis d N a mass).sum_repr
    (EuclideanSpace.single x (1 : ℝ))
  -- Replace the abstract repr coefficient with e_j(x).
  have h_repr_eq : ∀ j, (massEigenvectorBasis d N a mass).repr
      (EuclideanSpace.single x (1 : ℝ)) j =
      (massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) x := by
    intro j
    rw [OrthonormalBasis.repr_apply_apply]
    -- ⟨b j, single x 1⟩ via the dotProduct convention used elsewhere:
    -- inner a b = b.ofLp ⬝ᵥ star a.ofLp.
    change ((EuclideanSpace.single x (1 : ℝ)).ofLp ⬝ᵥ
      star (massEigenvectorBasis d N a mass j).ofLp) = _
    simp [dotProduct, star_trivial,
      ite_mul, one_mul, zero_mul, Finset.sum_ite_eq', Finset.mem_univ]
  simp_rw [h_repr_eq] at h_repr
  -- h_repr : ∑ j, e_j(x) • b_j = single x 1 (in EuclideanSpace).
  -- Cast both sides to Pi via ofLp.
  have h_ofLp : ((∑ j, (massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) x •
        (massEigenvectorBasis d N a mass j) : EuclideanSpace ℝ _) :
        FinLatticeSites d N → ℝ) =
      ((EuclideanSpace.single x (1 : ℝ)) : FinLatticeSites d N → ℝ) :=
    congrArg WithLp.ofLp h_repr
  -- Apply at y.
  have h_y := congrFun h_ofLp y
  -- Compute LHS: (∑_j c_j • b_j).ofLp y = ∑_j c_j * b_j.ofLp y.
  -- Convert (∑_j c_j • b_j).ofLp y → ∑_j c_j * (b_j).ofLp y.
  simp only [WithLp.ofLp_sum, WithLp.ofLp_smul] at h_y
  -- Now h_y is `(∑_j c_j • (b_j).ofLp) y = (single x 1).ofLp y`.
  rw [show (∑ j, ((massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) x •
        ((massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) :
          FinLatticeSites d N → ℝ))) y =
      ∑ j, (massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) x *
        (massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) y from by
    rw [Finset.sum_apply]
    refine Finset.sum_congr rfl ?_
    intro j _
    rfl] at h_y
  -- Compute RHS: (single x 1).ofLp y = if y = x then 1 else 0.
  rw [show ((EuclideanSpace.single x (1 : ℝ)) : FinLatticeSites d N → ℝ) y =
      if y = x then (1 : ℝ) else 0 from by rw [PiLp.single_apply]] at h_y
  exact h_y

/-- Spectral expansion: `ω(δ_x) = ∑_j γ_j(x) · ξ_j(ω)`. -/
private lemma omega_eval_delta_eq_sum_gamma_xi
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (x : FinLatticeSites d N)
    (ω : Configuration (FinLatticeField d N)) :
    ω (Pi.single x (1 : ℝ)) =
    ∑ j, gffEigenCoeff d N a mass j x *
        gffOrthonormalCoord d N a mass ha hmass j ω := by
  -- Step 1: Pi.single x 1 = ∑ j, e_j(x) • (e_j : Pi) (pointwise via eigenbasis_completeness).
  have h_pi_eq : (Pi.single x (1 : ℝ) : FinLatticeField d N) =
      ∑ j, (massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) x •
        ((massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) :
          FinLatticeSites d N → ℝ) := by
    funext y
    rw [show (Pi.single x (1 : ℝ) : FinLatticeField d N) y =
        (if y = x then (1 : ℝ) else 0) from Pi.single_apply _ _ _]
    rw [← eigenbasis_completeness d N a mass x y]
    rw [Finset.sum_apply]
    refine Finset.sum_congr rfl ?_
    intro j _
    rfl
  -- Step 2: apply ω-linearity.
  rw [h_pi_eq, map_sum]
  refine Finset.sum_congr rfl ?_
  intro j _
  rw [map_smul]
  show (massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) x •
      ω (fun y => (massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) y) =
    gffEigenCoeff d N a mass j x *
      gffOrthonormalCoord d N a mass ha hmass j ω
  unfold gffEigenCoeff gffOrthonormalCoord
  simp only [smul_eq_mul]
  have h_pos : (0 : ℝ) < a^d * massEigenvalues d N a mass j :=
    mul_pos (pow_pos ha d)
      (massOperatorMatrix_eigenvalues_pos d N a mass ha hmass j)
  have h_sqrt_ne : Real.sqrt (a^d * massEigenvalues d N a mass j) ≠ 0 :=
    (Real.sqrt_pos.mpr h_pos).ne'
  field_simp

/-- The site variance equals `∑_j γ_j(x)²`. -/
private lemma gffSiteVariance_eq_sum_gamma_sq
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass) (x : FinLatticeSites d N) :
    gffSiteVariance d N a mass ha hmass x =
    ∑ j, (gffEigenCoeff d N a mass j x)^2 := by
  unfold gffSiteVariance gffEigenCoeff
  rw [Finset.mul_sum]
  refine Finset.sum_congr rfl ?_
  intro j _
  have h_pos : (0 : ℝ) < a^d * massEigenvalues d N a mass j :=
    mul_pos (pow_pos ha d)
      (massOperatorMatrix_eigenvalues_pos d N a mass ha hmass j)
  have h_sqrt_sq : Real.sqrt (a^d * massEigenvalues d N a mass j) ^ 2 =
      a^d * massEigenvalues d N a mass j :=
    Real.sq_sqrt h_pos.le
  rw [show ((massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) x /
        Real.sqrt (a^d * massEigenvalues d N a mass j))^2 =
      ((massEigenvectorBasis d N a mass j : EuclideanSpace ℝ _) x)^2 /
        (Real.sqrt (a^d * massEigenvalues d N a mass j))^2 from by
    rw [div_pow]]
  rw [h_sqrt_sq, mul_comm (a^d) (massEigenvalues d N a mass j)]
  field_simp

/-- The explicit coefficient appearing in the multinomial Wick expansion:
`(k!/∏α_j!) · ∏ γ_j(x)^{α_j}` for multi-indices of total degree `k`,
zero otherwise. Used in `siteWickMonomial_eigenbasis_expansion`. -/
private noncomputable def siteWickExpansionCoeff
    (a mass : ℝ) (k : ℕ) (x : FinLatticeSites d N)
    (α : FinLatticeSites d N → ℕ) : ℝ :=
  if MultiIndexLattice.totalDegree α = k then
    ((k.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
      ∏ j, gffEigenCoeff d N a mass j x ^ (α j)
  else 0

/-- The explicit-coefficient function vanishes outside the right degree. -/
private lemma siteWickExpansionCoeff_eq_zero
    (a mass : ℝ) (k : ℕ) (x : FinLatticeSites d N)
    (α : FinLatticeSites d N → ℕ)
    (hα : MultiIndexLattice.totalDegree α ≠ k) :
    siteWickExpansionCoeff d N a mass k x α = 0 := by
  unfold siteWickExpansionCoeff
  rw [if_neg hα]

end GaussianField

/-! ## Multivariate Wick multinomial expansion (generic, not GFF-specific)

The statement below is a textbook fact about the 1D `wickMonomial`
recursion over an arbitrary `Fintype` index set — no lattice or
Gaussian-field content is involved. It naturally lives alongside
`wickMonomial` itself and is a clean candidate for upstreaming to
`Mathlib/Analysis/SpecialFunctions/.../Hermite/`. -/

/-- Multi-indices on a finite index set `ι` of total degree exactly `k`.
Each component is automatically bounded by `k`, giving a finite
enumeration via `Fintype.piFinset`. -/
noncomputable def multiIndicesOfTotalDegree
    (ι : Type*) [Fintype ι] [DecidableEq ι] (k : ℕ) :
    Finset (ι → ℕ) :=
  (Fintype.piFinset (fun _ : ι => Finset.range (k + 1))).filter
    (fun α => ∑ j, α j = k)

/-- **Multivariate Wick multinomial expansion** (textbook lemma, taken
as axiom).

For any finite index set `ι`, any `γ ξ : ι → ℝ` and any `k : ℕ`,

  `wickMonomial k (∑ γ²) (∑ γ_j ξ_j) =
     ∑_{|α|=k} (k! / ∏ α_j!) · (∏ γ_j^{α_j}) · ∏_j wickMonomial α_j 1 ξ_j`.

This is a polynomial identity in the `γ_j`'s and `ξ_j`'s — no
probabilistic or analytic content. Equivalent to matching
`t^k`-coefficients in the formal generating-function identity

  `exp(t · ∑ γ_j ξ_j − (∑ γ_j²) · t²/2)
     = ∏_j exp(t · γ_j ξ_j − γ_j² · t² / 2)`

after expanding each side via the 1D Wick exponential
`exp(t · x − v · t²/2) = ∑_n (t^n/n!) · wickMonomial n v x`.

**References:** Janson, *Gaussian Hilbert Spaces* (CUP 1997), §3.1
(univariate Hermite) and §3.4 (multivariate generalisation); Glimm &
Jaffe, *Quantum Physics*, §6.1 (Wick ordering); Nourdin & Peccati,
*Normal Approximations with Malliavin Calculus*, §2.7.

**Proof:** First specialise the 1D Wick recursion to a binomial-type
addition formula
`W_n(c₁+c₂, x+y) = ∑_k C(n,k) W_k(c₁,x) W_{n-k}(c₂,y)`
(`wickMonomial_add_add` in `SchwartzNuclear/HermiteWick.lean`),
proved by two-step induction on `n` plus Pascal's rule
(`Nat.choose_succ_succ`) and the choose-absorption identities
(`Nat.add_one_mul_choose_eq`).  Then induct on the support Finset
`s : Finset ι`, peeling off one index `j₀` at a time via
`wickMonomial_add_add` and `wickMonomial_homogeneity`, and invoke
the Function.update bijection
`(m, α) ↔ Function.update α j₀ m`
to pass between the multi-indices `α'` of total degree `k` and the
nested decomposition `(α' j₀, α' restricted to s)`.  The coefficient
identity at the final step is `C(k,m) · (k-m)! · m! = k!`.

The proof itself lives in `SchwartzNuclear/HermiteWick.lean`; this is
just a wrapper that uses the locally-defined `multiIndicesOfTotalDegree`. -/
theorem wickMonomial_pow_sum_expansion_of_totalDegree
    {ι : Type*} [Fintype ι] [DecidableEq ι]
    (γ ξ : ι → ℝ) (k : ℕ) :
    wickMonomial k (∑ j, (γ j) ^ 2) (∑ j, γ j * ξ j) =
    ∑ α ∈ multiIndicesOfTotalDegree ι k,
      ((k.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
      (∏ j, γ j ^ (α j)) *
      (∏ j, wickMonomial (α j) 1 (ξ j)) :=
  wickMonomial_pow_sum_expansion γ ξ k

namespace GaussianField

variable (d N : ℕ) [NeZero N]

/-- The lattice-side `multiIndicesOfDegree d N k` equals the generic
`multiIndicesOfTotalDegree (FinLatticeSites d N) k`. The two differ
only by the `MultiIndexLattice.totalDegree` wrapper unfolding to
`∑ j, α j`. -/
private lemma multiIndicesOfDegree_eq_generic (k : ℕ) :
    multiIndicesOfDegree d N k =
      multiIndicesOfTotalDegree (FinLatticeSites d N) k := by
  classical
  ext α
  unfold multiIndicesOfDegree multiIndicesOfTotalDegree
  rw [Finset.mem_filter, Finset.mem_filter]
  constructor
  · rintro ⟨hα, hdegree⟩
    refine ⟨?_, ?_⟩
    · rw [Fintype.mem_piFinset]
      intro a
      rw [multiIndicesUpToDegree, Fintype.mem_piFinset] at hα
      exact hα a
    · simpa [MultiIndexLattice.totalDegree] using hdegree
  · rintro ⟨hα, hdegree⟩
    refine ⟨?_, ?_⟩
    · rw [multiIndicesUpToDegree, Fintype.mem_piFinset]
      exact Fintype.mem_piFinset.mp hα
    · simpa [MultiIndexLattice.totalDegree] using hdegree

/-- **Site Wick monomial expansion in the eigenbasis.**

The single-site Wick monomial `:φ(x)^k:_{c_a(x)}`
(with the local Wick constant `c_a(x) = gffSiteVariance d N a mass ha hmass x
= (a^d)⁻¹ Σ_k λ_k⁻¹ e_k(x)²`) admits an explicit linear expansion in
eigenbasis multi-indices:

  `:φ(x)^k:_{c_a(x)} = ∑_{|α| = k} coeff(α, x, k) · :ξ^α:_1`,

where the coefficients are explicit polynomial combinations of the
eigenvector values `e_k(x)` and eigenvalue powers `λ_k^{-1/2}`. The
`totalDegree α = k` constraint is the key content: only multi-indices
of *exact* total degree `k` appear, because the local Wick subtraction
with the matched site variance `c_a(x)` cancels exactly the lower-degree
contractions.

**Proof:** Combine (1) the spectral expansion
`omega_eval_delta_eq_sum_gamma_xi`: `ω(δ_x) = ∑_j γ_j(x) · ξ_j(ω)`,
(2) the variance identity
`gffSiteVariance_eq_sum_gamma_sq`: `c_a(x) = ∑_j γ_j(x)²`,
and (3) the multivariate Wick multinomial
expansion `wickMonomial_pow_sum_expansion`. -/
theorem siteWickMonomial_eigenbasis_expansion
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (k : ℕ) (x : FinLatticeSites d N) :
    ∃ (coeff : (FinLatticeSites d N → ℕ) → ℝ),
      (∀ α, coeff α ≠ 0 → MultiIndexLattice.totalDegree α = k) ∧
      ∀ ω : Configuration (FinLatticeField d N),
        wickMonomial k (gffSiteVariance d N a mass ha hmass x)
            (ω (Pi.single x 1)) =
        ∑ α ∈ multiIndicesOfDegree d N k, coeff α *
          gffMultiWickMonomial d N a mass ha hmass α ω := by
  refine ⟨siteWickExpansionCoeff d N a mass k x, ?_, ?_⟩
  · -- The coefficient vanishes off the right-degree slice.
    intro α hα
    by_contra hcontra
    exact hα (siteWickExpansionCoeff_eq_zero d N a mass k x α hcontra)
  · -- The expansion equation.
    intro ω
    -- Substitute the spectral expansion of ω(δ_x) and gffSiteVariance.
    rw [gffSiteVariance_eq_sum_gamma_sq d N a mass ha hmass x,
        omega_eval_delta_eq_sum_gamma_xi d N a mass ha hmass x ω]
    -- LHS now: wickMonomial k (∑ γ²) (∑ γ_j · ξ_j(ω)). Apply the (generic)
    -- multinomial axiom — it produces a sum over `multiIndicesOfTotalDegree`,
    -- which we then identify with `multiIndicesOfDegree d N k`.
    rw [wickMonomial_pow_sum_expansion_of_totalDegree
          (γ := fun j => gffEigenCoeff d N a mass j x)
          (ξ := fun j => gffOrthonormalCoord d N a mass ha hmass j ω) (k := k)]
    rw [multiIndicesOfDegree_eq_generic d N k]
    -- Match summand by summand.
    refine Finset.sum_congr rfl ?_
    intro α hα
    have h_deg : MultiIndexLattice.totalDegree α = k := by
      unfold multiIndicesOfTotalDegree at hα
      simp only [Finset.mem_filter] at hα
      exact hα.2
    unfold siteWickExpansionCoeff
    rw [if_pos h_deg]
    -- Both sides: (k!/∏α!) · (∏ γ_j(x)^{α_j}) · (∏_j wickMonomial α_j 1 (ξ_j ω)),
    -- where on the RHS the wickMonomial product is wrapped as `gffMultiWickMonomial α ω`.
    rfl

/-! ## 2-site Wick power formula on the lattice GFF

The Janson-Hilbert two-site formula specialised to single-site Wick
powers under the lattice GFF measure:

  `∫ :φ(x)^n:_{c_a(x)} · :φ(y)^m:_{c_a(y)} dμ_GFF
     = δ_{n,m} · n! · gffPositionCovariance(x, y)^n`

where `gffPositionCovariance(x, y) = Σ_j γ_j(x) · γ_j(y) =
(a^d)⁻¹ · Σ_j (e_j(x) e_j(y) / λ_j)` is the position-space covariance
kernel (the (x, y)-entry of the inverse mass operator, scaled).

Used by pphi2 to discharge the cross-term orthogonality
(`canonicalCrossTerm_inner_eq_zero`) and the diagonal piece needed by
the per-cross-term L² bound (`canonicalCrossTerm_l2_sq_le`).

**Proof sketch** (deferred):
1. Substitute `gffSiteVariance_eq_sum_gamma_sq` and
   `omega_eval_delta_eq_sum_gamma_xi` to rewrite each Wick monomial
   as `wickMonomial k (Σ γ_j(x)²) (Σ γ_j(x) ξ_j(ω))`.
2. Apply `wickMonomial_pow_sum_expansion_of_totalDegree` to expand
   each side into a multi-index sum over `gffMultiWickMonomial`.
3. Multiply, integrate (linearity), apply
   `gffMultiWickMonomial_orthogonality` to vanish off-diagonal pairs.
4. The diagonal `α = β` only contributes when `|α| = n = m`. The
   resulting sum reduces via the multinomial theorem
     `Σ_{|α|=n} (n! / ∏ α_j!) · ∏ a_j^{α_j} = (Σ_j a_j)^n`
   applied to `a_j = γ_j(x) γ_j(y)` to give `n! · (Σ γ_j(x) γ_j(y))^n`.
-/

/-- Position-space covariance kernel of the lattice GFF:
`C(x, y) = Σ_j γ_j(x) · γ_j(y) = (a^d)⁻¹ · Σ_j (e_j(x) e_j(y)) / λ_j`.

The (x, y)-entry of the inverse mass operator (scaled by `(a^d)⁻¹`),
expressed via the orthogonalised eigenvector coefficients. The
diagonal `gffPositionCovariance x x = gffSiteVariance x` is the
on-site variance (the auto-covariance Wick subtracts). -/
noncomputable def gffPositionCovariance
    (a mass : ℝ) (x y : FinLatticeSites d N) : ℝ :=
  ∑ j, gffEigenCoeff d N a mass j x * gffEigenCoeff d N a mass j y

/-- The position covariance is symmetric. -/
lemma gffPositionCovariance_symm
    (a mass : ℝ) (x y : FinLatticeSites d N) :
    gffPositionCovariance d N a mass x y =
    gffPositionCovariance d N a mass y x := by
  unfold gffPositionCovariance
  refine Finset.sum_congr rfl fun j _ => ?_
  ring

/-- The position covariance on the diagonal equals the site variance. -/
lemma gffPositionCovariance_self
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass) (x : FinLatticeSites d N) :
    gffPositionCovariance d N a mass x x =
    gffSiteVariance d N a mass ha hmass x := by
  unfold gffPositionCovariance
  rw [gffSiteVariance_eq_sum_gamma_sq d N a mass ha hmass x]
  refine Finset.sum_congr rfl fun j _ => ?_
  ring

/-- **Eigenbasis ↔ operator covariance bridge.** The eigenbasis position covariance
`gffPositionCovariance x y = Σ_j γ_j(x)γ_j(y)` equals the Glimm–Jaffe operator covariance
`⟨T_GJ δ_x, T_GJ δ_y⟩ = covariance (latticeCovarianceGJ) δ_x δ_y`. Both are
`(a^d)⁻¹ Σ_k λ_k⁻¹ e_k(x) e_k(y)`. This connects the smeared Wick kernel (stated in the eigenbasis)
to the operator-form lattice covariance used downstream (e.g. pphi2's `wickConstant`). -/
lemma gffPositionCovariance_eq_covarianceGJ
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass) (x y : FinLatticeSites d N) :
    gffPositionCovariance d N a mass x y =
    GaussianField.covariance (latticeCovarianceGJ d N a mass ha hmass)
      (Pi.single x (1 : ℝ)) (Pi.single y (1 : ℝ)) := by
  rw [lattice_covariance_GJ_eq_spectral d N a mass ha hmass]
  unfold gffPositionCovariance gffEigenCoeff
  rw [Finset.mul_sum]
  refine Finset.sum_congr rfl fun k _ => ?_
  simp only [Pi.single_apply, mul_ite, mul_one, mul_zero, Finset.sum_ite_eq',
    Finset.mem_univ, if_true]
  have hpos : (0 : ℝ) < a ^ d * massEigenvalues d N a mass k :=
    mul_pos (pow_pos ha d) (massOperatorMatrix_eigenvalues_pos d N a mass ha hmass k)
  have hsq : Real.sqrt (a ^ d * massEigenvalues d N a mass k) *
      Real.sqrt (a ^ d * massEigenvalues d N a mass k) = a ^ d * massEigenvalues d N a mass k :=
    Real.mul_self_sqrt hpos.le
  have ha_ne : (a ^ d : ℝ) ≠ 0 := (pow_pos ha d).ne'
  have hlam_ne : massEigenvalues d N a mass k ≠ 0 :=
    (massOperatorMatrix_eigenvalues_pos d N a mass ha hmass k).ne'
  rw [div_mul_div_comm, hsq]
  field_simp

/-- Diagonal form of the eigenbasis ↔ operator covariance bridge:
`gffSiteVariance x = ⟨T_GJ δ_x, T_GJ δ_x⟩`. -/
lemma gffSiteVariance_eq_covarianceGJ
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass) (x : FinLatticeSites d N) :
    gffSiteVariance d N a mass ha hmass x =
    GaussianField.covariance (latticeCovarianceGJ d N a mass ha hmass)
      (Pi.single x (1 : ℝ)) (Pi.single x (1 : ℝ)) := by
  rw [← gffPositionCovariance_self d N a mass ha hmass x,
    gffPositionCovariance_eq_covarianceGJ d N a mass ha hmass x x]

/-- Eigenbasis expansion of a single-site Wick monomial as an explicit
sum over multi-indices of total degree `n`. Refines
`siteWickMonomial_eigenbasis_expansion` by exposing the explicit
coefficient `(n! / ∏ α_j!) · ∏ γ_j(x)^{α_j}` rather than wrapping it
behind an existential. -/
private lemma wickMonomial_at_site_eq_eigen_sum
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (n : ℕ) (x : FinLatticeSites d N)
    (ω : Configuration (FinLatticeField d N)) :
    wickMonomial n (gffSiteVariance d N a mass ha hmass x)
        (ω (Pi.single x 1)) =
    ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
      ((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
        (∏ j, gffEigenCoeff d N a mass j x ^ (α j)) *
        gffMultiWickMonomial d N a mass ha hmass α ω := by
  rw [gffSiteVariance_eq_sum_gamma_sq d N a mass ha hmass x,
      omega_eval_delta_eq_sum_gamma_xi d N a mass ha hmass x ω,
      wickMonomial_pow_sum_expansion_of_totalDegree
        (fun j => gffEigenCoeff d N a mass j x)
        (fun j => gffOrthonormalCoord d N a mass ha hmass j ω) n]
  rfl

private lemma integrable_pow_gaussianReal_one (k : ℕ) :
    MeasureTheory.Integrable (fun x : ℝ => x ^ k)
      (ProbabilityTheory.gaussianReal (0 : ℝ) (1 : NNReal)) := by
  cases k with
  | zero =>
      simp only [pow_zero]
      exact MeasureTheory.integrable_const (1 : ℝ)
  | succ k =>
      have h_mem0 :=
        ProbabilityTheory.memLp_id_gaussianReal
          (μ := (0 : ℝ)) (v := (1 : NNReal)) (p := (Nat.succ k : NNReal))
      have h_mem : MeasureTheory.MemLp id (Nat.succ k)
          (ProbabilityTheory.gaussianReal (0 : ℝ) (1 : NNReal)) := h_mem0
      have h_norm : MeasureTheory.Integrable (fun x : ℝ => ‖x‖ ^ Nat.succ k)
          (ProbabilityTheory.gaussianReal (0 : ℝ) (1 : NNReal)) := by
        simpa using h_mem.integrable_norm_pow'
      have h_meas :
          MeasureTheory.AEStronglyMeasurable (fun x : ℝ => x ^ Nat.succ k)
            (ProbabilityTheory.gaussianReal (0 : ℝ) (1 : NNReal)) :=
        (continuous_id.pow _).aestronglyMeasurable
      rw [← MeasureTheory.integrable_norm_iff h_meas]
      simpa [norm_pow] using h_norm

private lemma integrable_polynomial_gaussianReal_one (p : Polynomial ℝ) :
    MeasureTheory.Integrable (fun x : ℝ => p.eval x)
      (ProbabilityTheory.gaussianReal (0 : ℝ) (1 : NNReal)) := by
  have h_eval :
      (fun x : ℝ => p.eval x) =
        (fun x => ∑ i ∈ Finset.range (p.natDegree + 1), p.coeff i * x ^ i) := by
    ext x
    rw [Polynomial.eval_eq_sum_range]
  rw [h_eval]
  exact MeasureTheory.integrable_finsetSum _ fun i hi =>
    (integrable_pow_gaussianReal_one i).const_mul (p.coeff i)

private lemma integrable_wickMonomial_mul_gaussianReal_one (m n : ℕ) :
    MeasureTheory.Integrable (fun x : ℝ => wickMonomial m 1 x * wickMonomial n 1 x)
      (ProbabilityTheory.gaussianReal (0 : ℝ) (1 : NNReal)) := by
  have h_eq :
      (fun x : ℝ => wickMonomial m 1 x * wickMonomial n 1 x) =
        fun x => ((hermiteR m) * (hermiteR n)).eval x := by
    funext x
    rw [Polynomial.eval_mul, wick_eq_hermiteR m 1 (by norm_num : (0 : ℝ) < 1),
      wick_eq_hermiteR n 1 (by norm_num : (0 : ℝ) < 1)]
    simp [scaledHermite]
  rw [h_eq]
  exact integrable_polynomial_gaussianReal_one ((hermiteR m) * (hermiteR n))

private lemma integrable_gffMultiWickMonomial_mul
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (α β : FinLatticeSites d N → ℕ) :
    MeasureTheory.Integrable
      (fun ω => gffMultiWickMonomial d N a mass ha hmass α ω *
        gffMultiWickMonomial d N a mass ha hmass β ω)
      (latticeGaussianMeasure d N a mass ha hmass) := by
  let f : (FinLatticeSites d N → ℝ) → ℝ := fun y =>
    ∏ k, wickMonomial (α k) 1 (y k) * wickMonomial (β k) 1 (y k)
  have hf :
      MeasureTheory.Integrable f
        (MeasureTheory.Measure.pi fun _ : FinLatticeSites d N =>
          ProbabilityTheory.gaussianReal (0 : ℝ) (1 : NNReal)) := by
    simpa [f] using
      (MeasureTheory.Integrable.fintype_prod
        (μ := fun _ : FinLatticeSites d N =>
          ProbabilityTheory.gaussianReal (0 : ℝ) (1 : NNReal))
        (f := fun k x => wickMonomial (α k) 1 x * wickMonomial (β k) 1 x)
        (fun k => integrable_wickMonomial_mul_gaussianReal_one (α k) (β k)))
  have h_push :
      MeasureTheory.Integrable f
        ((latticeGaussianMeasure d N a mass ha hmass).map
          (gffOrthonormalProj d N a mass ha hmass)) := by
    rw [gffOrthonormalProj_pushforward_eq_stdGaussian d N a mass ha hmass]
    exact hf
  refine (h_push.comp_measurable
      (gffOrthonormalProj_measurable d N a mass ha hmass)).congr ?_
  filter_upwards with ω
  simp [f, gffOrthonormalProj, gffMultiWickMonomial, ← Finset.prod_mul_distrib]

private lemma multiIndicesOfTotalDegree_eq_piAntidiag (n : ℕ) :
    multiIndicesOfTotalDegree (FinLatticeSites d N) n =
      Finset.piAntidiag Finset.univ n := by
  ext α
  constructor
  · intro hα
    unfold multiIndicesOfTotalDegree at hα
    simp only [Finset.mem_filter] at hα
    exact (Finset.mem_piAntidiag).2 ⟨hα.2, fun i hi => Finset.mem_univ i⟩
  · intro hα
    have hsum : ∑ j, α j = n := (Finset.mem_piAntidiag.mp hα).1
    simp only [multiIndicesOfTotalDegree, Finset.mem_filter, Fintype.mem_piFinset,
      Finset.mem_range]
    refine ⟨?_, hsum⟩
    intro i
    have h_le : α i ≤ ∑ j, α j := by
      simpa using
        (Finset.single_le_sum (fun j _ => Nat.zero_le (α j)) (Finset.mem_univ i))
    exact Nat.lt_succ_of_le (hsum ▸ h_le)

private lemma sum_pow_eq_sum_multiIndices (f : FinLatticeSites d N → ℝ) (n : ℕ) :
    (∑ j, f j) ^ n =
      ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
        ((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
          ∏ j, f j ^ α j := by
  rw [multiIndicesOfTotalDegree_eq_piAntidiag d N n, Finset.sum_pow_eq_sum_piAntidiag]
  refine Finset.sum_congr rfl ?_
  intro α hα
  have hsum : ∑ j, α j = n := (Finset.mem_piAntidiag.mp hα).1
  have hprod_pos : (0 : ℝ) < ∏ j, ((α j).factorial : ℝ) := by
    apply Finset.prod_pos
    intro j hj
    exact_mod_cast Nat.factorial_pos (α j)
  have hprod_ne : (∏ j, ((α j).factorial : ℝ)) ≠ 0 := ne_of_gt hprod_pos
  have hmult :
      (Nat.multinomial Finset.univ α : ℝ) =
        ((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) := by
    apply (eq_div_iff hprod_ne).2
    calc
      (Nat.multinomial Finset.univ α : ℝ) * ∏ j, ((α j).factorial : ℝ)
          = (∏ j, ((α j).factorial : ℝ)) *
              (Nat.multinomial Finset.univ α : ℝ) := by ring
      _ = (n.factorial : ℝ) := by
        have hspec :
            (∏ j ∈ (Finset.univ : Finset (FinLatticeSites d N)), (α j).factorial) *
              Nat.multinomial Finset.univ α =
            (∑ j ∈ (Finset.univ : Finset (FinLatticeSites d N)), α j).factorial :=
          Nat.multinomial_spec Finset.univ α
        exact_mod_cast (hspec.trans (by simp [hsum]))
  simp [hmult]

/-- **2-site Wick power formula on the lattice GFF.** For sites `x, y`
and Wick powers `n, m`, the integral of the product of single-site
Wick monomials under the lattice GFF measure equals `n! · C(x, y)^n`
when `n = m`, zero otherwise. The Janson–Hilbert orthogonality applied
in the eigenbasis.

Used downstream to discharge cross-term orthogonality and the
diagonal piece of per-cross-term L² bounds. -/
theorem gff_wickPower_two_site_inner
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (n m : ℕ) (x y : FinLatticeSites d N) :
    ∫ ω, wickMonomial n (gffSiteVariance d N a mass ha hmass x)
            (ω (Pi.single x 1)) *
          wickMonomial m (gffSiteVariance d N a mass ha hmass y)
            (ω (Pi.single y 1))
        ∂(latticeGaussianMeasure d N a mass ha hmass) =
    if n = m then
      (n.factorial : ℝ) * (gffPositionCovariance d N a mass x y) ^ n
    else 0 := by
  classical
  let coeff :=
    fun (k : ℕ) (z : FinLatticeSites d N) (α : FinLatticeSites d N → ℕ) =>
      ((k.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
        ∏ j, gffEigenCoeff d N a mass j z ^ α j
  let diagFac := fun (α : FinLatticeSites d N → ℕ) => ∏ j, ((α j).factorial : ℝ)
  have h_summand_int :
      ∀ α β, MeasureTheory.Integrable
        (fun ω =>
          (coeff n x α * gffMultiWickMonomial d N a mass ha hmass α ω) *
            (coeff m y β * gffMultiWickMonomial d N a mass ha hmass β ω))
        (latticeGaussianMeasure d N a mass ha hmass) := by
    intro α β
    have h_base := (integrable_gffMultiWickMonomial_mul d N a mass ha hmass α β).const_mul
      (coeff n x α * coeff m y β)
    refine h_base.congr ?_
    filter_upwards with ω
    ring
  have h_expand :
      ∀ ω : Configuration (FinLatticeField d N),
        wickMonomial n (gffSiteVariance d N a mass ha hmass x) (ω (Pi.single x 1)) *
            wickMonomial m (gffSiteVariance d N a mass ha hmass y) (ω (Pi.single y 1)) =
          (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
            coeff n x α * gffMultiWickMonomial d N a mass ha hmass α ω) *
          (∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) m,
            coeff m y β * gffMultiWickMonomial d N a mass ha hmass β ω) := by
    intro ω
    rw [wickMonomial_at_site_eq_eigen_sum d N a mass ha hmass n x ω,
      wickMonomial_at_site_eq_eigen_sum d N a mass ha hmass m y ω]
  simp_rw [h_expand]
  have h_distrib :
      ∀ ω : Configuration (FinLatticeField d N),
        (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
            coeff n x α * gffMultiWickMonomial d N a mass ha hmass α ω) *
          (∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) m,
            coeff m y β * gffMultiWickMonomial d N a mass ha hmass β ω) =
        ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          ∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) m,
            (coeff n x α * gffMultiWickMonomial d N a mass ha hmass α ω) *
              (coeff m y β * gffMultiWickMonomial d N a mass ha hmass β ω) := by
    intro ω
    rw [Finset.sum_mul_sum]
  simp_rw [h_distrib]
  rw [MeasureTheory.integral_finsetSum _ (fun α hα => by
    apply MeasureTheory.integrable_finsetSum
    intro β hβ
    exact h_summand_int α β)]
  have h_outer :
      (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          ∫ ω, ∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) m,
            (coeff n x α * gffMultiWickMonomial d N a mass ha hmass α ω) *
              (coeff m y β * gffMultiWickMonomial d N a mass ha hmass β ω)
          ∂(latticeGaussianMeasure d N a mass ha hmass)) =
        ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          ∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) m,
            (coeff n x α * coeff m y β) *
              (if α = β then diagFac α else 0) := by
    refine Finset.sum_congr rfl ?_
    intro α hα
    rw [MeasureTheory.integral_finsetSum _ (fun β hβ => h_summand_int α β)]
    refine Finset.sum_congr rfl ?_
    intro β hβ
    rw [show
        (fun ω =>
          (coeff n x α * gffMultiWickMonomial d N a mass ha hmass α ω) *
            (coeff m y β * gffMultiWickMonomial d N a mass ha hmass β ω)) =
          (fun ω =>
            (coeff n x α * coeff m y β) *
              (gffMultiWickMonomial d N a mass ha hmass α ω *
                gffMultiWickMonomial d N a mass ha hmass β ω)) from by
        ext ω
        ring]
    rw [MeasureTheory.integral_const_mul,
      gffMultiWickMonomial_orthogonality d N a mass ha hmass]
    simp [diagFac]
  rw [h_outer]
  by_cases hnm : n = m
  · subst hnm
    rw [if_pos rfl]
    have h_collapse :
        (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          ∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
            (coeff n x α * coeff n y β) * (if α = β then diagFac α else 0)) =
        ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          (coeff n x α * coeff n y α) * diagFac α := by
      refine Finset.sum_congr rfl ?_
      intro α hα
      rw [Finset.sum_eq_single_of_mem α hα]
      · simp [diagFac]
      · intro β hβ hβα
        have hneq : α ≠ β := by
          intro h
          exact hβα h.symm
        simp [hneq, diagFac]
    have h_diag_term :
        ∀ α, α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n →
          (coeff n x α * coeff n y α) * diagFac α =
            (n.factorial : ℝ) *
              (((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
                ∏ j, (gffEigenCoeff d N a mass j x *
                  gffEigenCoeff d N a mass j y) ^ α j) := by
      intro α hα
      have hprod_mul :
          (∏ j, gffEigenCoeff d N a mass j x ^ α j) *
              (∏ j, gffEigenCoeff d N a mass j y ^ α j) =
            ∏ j, (gffEigenCoeff d N a mass j x *
              gffEigenCoeff d N a mass j y) ^ α j := by
        rw [← Finset.prod_mul_distrib]
        refine Finset.prod_congr rfl ?_
        intro j hj
        rw [mul_pow]
      have hprod_pos : (0 : ℝ) < ∏ j, ((α j).factorial : ℝ) := by
        apply Finset.prod_pos
        intro j hj
        exact_mod_cast Nat.factorial_pos (α j)
      have hprod_ne : (∏ j, ((α j).factorial : ℝ)) ≠ 0 := ne_of_gt hprod_pos
      calc
        (coeff n x α * coeff n y α) * diagFac α
            = (n.factorial : ℝ) *
                (((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
                  ((∏ j, gffEigenCoeff d N a mass j x ^ α j) *
                    (∏ j, gffEigenCoeff d N a mass j y ^ α j))) := by
              unfold coeff diagFac
              field_simp [hprod_ne]
        _ = (n.factorial : ℝ) *
              (((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
                ∏ j, (gffEigenCoeff d N a mass j x *
                  gffEigenCoeff d N a mass j y) ^ α j) := by
              simp [hprod_mul]
    calc
      (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          ∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
            (coeff n x α * coeff n y β) * (if α = β then diagFac α else 0))
          =
        ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          (coeff n x α * coeff n y α) * diagFac α := h_collapse
      _ =
        ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          (n.factorial : ℝ) *
            (((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
              ∏ j, (gffEigenCoeff d N a mass j x *
                gffEigenCoeff d N a mass j y) ^ α j) := by
            refine Finset.sum_congr rfl ?_
            intro α hα
            exact h_diag_term α hα
      _ = (n.factorial : ℝ) *
          (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
            ((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
              ∏ j, (gffEigenCoeff d N a mass j x *
                gffEigenCoeff d N a mass j y) ^ α j) := by
          rw [← Finset.mul_sum]
      _ = (n.factorial : ℝ) *
          (∑ j, gffEigenCoeff d N a mass j x *
            gffEigenCoeff d N a mass j y) ^ n := by
          rw [sum_pow_eq_sum_multiIndices d N
            (fun j => gffEigenCoeff d N a mass j x * gffEigenCoeff d N a mass j y) n]
      _ = (n.factorial : ℝ) * (gffPositionCovariance d N a mass x y) ^ n := by
          rw [gffPositionCovariance]
  · rw [if_neg hnm]
    apply Finset.sum_eq_zero
    intro α hα
    apply Finset.sum_eq_zero
    intro β hβ
    by_cases hαβ : α = β
    · subst hαβ
      unfold multiIndicesOfTotalDegree at hα hβ
      simp only [Finset.mem_filter] at hα hβ
      exact absurd (hα.2.symm.trans hβ.2) hnm
    · simp [hαβ, diagFac]

/-! ## Smeared (test-function) Wick power inner product

Generalises `gff_wickPower_two_site_inner` from single sites `δ_x` to arbitrary
smeared lattice fields `f, g : FinLatticeField d N`, by `ω`-linearity. The
single-site case is recovered with `f = Pi.single x 1`
(then `gffSmearedCoeff … f j = gffEigenCoeff … j x`). Used by pphi2 to evaluate
the connected four-point `u₄` of the interacting lattice measure on a general
(smeared) test function, where the `gff_wickPower_two_site_inner` form does not
apply because the embedded test function is supported on many sites. -/

/-- Smeared eigen-coefficient `Γ_j(f) = ∑_x f(x) · γ_j(x)`. Public API for the smeared Wick kernel. -/
noncomputable def gffSmearedCoeff
    (a mass : ℝ) (f : FinLatticeField d N) (j : FinLatticeSites d N) : ℝ :=
  ∑ x, f x * gffEigenCoeff d N a mass j x

/-- Smeared position covariance `∑_j Γ_j(f) Γ_j(g) = ∑_{x,y} f(x) g(y) C(x,y) = ⟨φ(f)φ(g)⟩`.
Public API: the `n=m` value of `gff_wickPower_two_smeared_inner`; rewrite via the position-space
covariance with `gffSmearedCovariance_eq_sum_position`. -/
noncomputable def gffSmearedCovariance
    (a mass : ℝ) (f g : FinLatticeField d N) : ℝ :=
  ∑ j, gffSmearedCoeff d N a mass f j * gffSmearedCoeff d N a mass g j

/-- `gffSmearedCovariance f g = ∑_{x,y} f(x) g(y) C(x,y)` — the smeared covariance via the
position-space covariance kernel `gffPositionCovariance`. Connects the smeared Wick kernel to the
position-space `(C_a f)(z)` used downstream. -/
lemma gffSmearedCovariance_eq_sum_position
    (a mass : ℝ) (f g : FinLatticeField d N) :
    gffSmearedCovariance d N a mass f g =
      ∑ x, ∑ y, f x * g y * gffPositionCovariance d N a mass x y := by
  unfold gffSmearedCovariance gffSmearedCoeff gffPositionCovariance
  simp_rw [Finset.sum_mul_sum]
  rw [Finset.sum_comm]
  refine Finset.sum_congr rfl fun x _ => ?_
  rw [Finset.sum_comm]
  refine Finset.sum_congr rfl fun y _ => ?_
  rw [Finset.mul_sum]
  refine Finset.sum_congr rfl fun j _ => ?_
  ring

/-- `gffSmearedCovariance f (δ_z) = ∑_x f(x) C(x,z) = (C_a f)(z)` — the single-site specialization
giving the position-space covariance applied to `f` at the vertex site `z`. -/
lemma gffSmearedCovariance_single_right
    (a mass : ℝ) (f : FinLatticeField d N) (z : FinLatticeSites d N) :
    gffSmearedCovariance d N a mass f (Pi.single z (1 : ℝ)) =
      ∑ x, f x * gffPositionCovariance d N a mass x z := by
  rw [gffSmearedCovariance_eq_sum_position]
  refine Finset.sum_congr rfl fun x _ => ?_
  simp only [Pi.single_apply, mul_ite, mul_one, mul_zero, ite_mul, zero_mul,
    Finset.sum_ite_eq', Finset.mem_univ, if_true]

/-- `gffSmearedCovariance f f = ∑_j Γ_j(f)²` — the smeared field variance `Var ω(f)`, i.e. the Wick
constant in `gff_wickPower_two_smeared_inner`. -/
lemma gffSmearedCovariance_self
    (a mass : ℝ) (f : FinLatticeField d N) :
    gffSmearedCovariance d N a mass f f = ∑ j, (gffSmearedCoeff d N a mass f j) ^ 2 := by
  unfold gffSmearedCovariance
  refine Finset.sum_congr rfl fun j _ => ?_
  rw [sq]

/-- `ω(f) = ∑_j Γ_j(f) · ξ_j(ω)` — smeared spectral expansion, by `ω`-linearity. -/
private lemma omega_eval_smeared_eq_sum_gamma_xi
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (f : FinLatticeField d N) (ω : Configuration (FinLatticeField d N)) :
    ω f = ∑ j, gffSmearedCoeff d N a mass f j *
        gffOrthonormalCoord d N a mass ha hmass j ω := by
  have hf : f = ∑ x, f x • (Pi.single x (1 : ℝ) : FinLatticeField d N) := by
    funext y
    simp only [Finset.sum_apply, Pi.smul_apply, Pi.single_apply, smul_eq_mul,
      mul_ite, mul_one, mul_zero, Finset.sum_ite_eq, Finset.mem_univ, if_true]
  conv_lhs => rw [hf, map_sum]
  simp_rw [map_smul, smul_eq_mul,
    omega_eval_delta_eq_sum_gamma_xi d N a mass ha hmass, Finset.mul_sum]
  rw [Finset.sum_comm]
  refine Finset.sum_congr rfl fun j _ => ?_
  rw [gffSmearedCoeff, Finset.sum_mul]
  refine Finset.sum_congr rfl fun x _ => ?_
  ring

/-- Smeared analogue of `wickMonomial_at_site_eq_eigen_sum`: the Wick monomial of
the smeared field, with matched Wick constant `∑_j Γ_j(f)² = Var ω(f)`, expands
over eigenbasis multi-indices of total degree `n`. -/
private lemma wickMonomial_at_smeared_eq_eigen_sum
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (n : ℕ) (f : FinLatticeField d N)
    (ω : Configuration (FinLatticeField d N)) :
    wickMonomial n (∑ j, (gffSmearedCoeff d N a mass f j) ^ 2) (ω f) =
      ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
        (((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
          ∏ j, gffSmearedCoeff d N a mass f j ^ (α j)) *
          gffMultiWickMonomial d N a mass ha hmass α ω := by
  rw [omega_eval_smeared_eq_sum_gamma_xi d N a mass ha hmass f ω,
    wickMonomial_pow_sum_expansion_of_totalDegree
      (γ := fun j => gffSmearedCoeff d N a mass f j)
      (ξ := fun j => gffOrthonormalCoord d N a mass ha hmass j ω) (k := n)]
  refine Finset.sum_congr rfl fun α _ => ?_
  rfl

/-- **Smeared two-field Wick power formula on the lattice GFF.** For smeared
lattice fields `f, g` and Wick powers `n, m` (each with the matched Wick constant
`∑_j Γ_j(·)² = Var ω(·)`), the integral of the product under the lattice GFF
measure equals `n! · (∑_j Γ_j(f) Γ_j(g))^n = n! · ⟨φ(f) φ(g)⟩^n` when `n = m`,
zero otherwise. Generalises `gff_wickPower_two_site_inner` from sites to test
functions; the single-site case is `f = Pi.single x 1`, `g = Pi.single y 1`.

This is the connected-correlator input for the perturbative four-point: with
`n = m = 4` it gives `∫ :φ(f)⁴: :φ(g)⁴: dμ_GFF = 4! · ⟨φ(f)φ(g)⟩⁴`, the
Wick/Mehler kernel. -/
theorem gff_wickPower_two_smeared_inner
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (n m : ℕ) (f g : FinLatticeField d N) :
    ∫ ω, wickMonomial n (∑ j, (gffSmearedCoeff d N a mass f j) ^ 2) (ω f) *
          wickMonomial m (∑ j, (gffSmearedCoeff d N a mass g j) ^ 2) (ω g)
        ∂(latticeGaussianMeasure d N a mass ha hmass) =
    if n = m then
      (n.factorial : ℝ) * (gffSmearedCovariance d N a mass f g) ^ n
    else 0 := by
  classical
  let coeff :=
    fun (k : ℕ) (h : FinLatticeField d N) (α : FinLatticeSites d N → ℕ) =>
      ((k.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
        ∏ j, gffSmearedCoeff d N a mass h j ^ α j
  let diagFac := fun (α : FinLatticeSites d N → ℕ) => ∏ j, ((α j).factorial : ℝ)
  have h_summand_int :
      ∀ α β, MeasureTheory.Integrable
        (fun ω =>
          (coeff n f α * gffMultiWickMonomial d N a mass ha hmass α ω) *
            (coeff m g β * gffMultiWickMonomial d N a mass ha hmass β ω))
        (latticeGaussianMeasure d N a mass ha hmass) := by
    intro α β
    have h_base := (integrable_gffMultiWickMonomial_mul d N a mass ha hmass α β).const_mul
      (coeff n f α * coeff m g β)
    refine h_base.congr ?_
    filter_upwards with ω
    ring
  have h_expand :
      ∀ ω : Configuration (FinLatticeField d N),
        wickMonomial n (∑ j, (gffSmearedCoeff d N a mass f j) ^ 2) (ω f) *
            wickMonomial m (∑ j, (gffSmearedCoeff d N a mass g j) ^ 2) (ω g) =
          (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
            coeff n f α * gffMultiWickMonomial d N a mass ha hmass α ω) *
          (∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) m,
            coeff m g β * gffMultiWickMonomial d N a mass ha hmass β ω) := by
    intro ω
    rw [wickMonomial_at_smeared_eq_eigen_sum d N a mass ha hmass n f ω,
      wickMonomial_at_smeared_eq_eigen_sum d N a mass ha hmass m g ω]
  simp_rw [h_expand]
  have h_distrib :
      ∀ ω : Configuration (FinLatticeField d N),
        (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
            coeff n f α * gffMultiWickMonomial d N a mass ha hmass α ω) *
          (∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) m,
            coeff m g β * gffMultiWickMonomial d N a mass ha hmass β ω) =
        ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          ∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) m,
            (coeff n f α * gffMultiWickMonomial d N a mass ha hmass α ω) *
              (coeff m g β * gffMultiWickMonomial d N a mass ha hmass β ω) := by
    intro ω
    rw [Finset.sum_mul_sum]
  simp_rw [h_distrib]
  rw [MeasureTheory.integral_finsetSum _ (fun α hα => by
    apply MeasureTheory.integrable_finsetSum
    intro β hβ
    exact h_summand_int α β)]
  have h_outer :
      (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          ∫ ω, ∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) m,
            (coeff n f α * gffMultiWickMonomial d N a mass ha hmass α ω) *
              (coeff m g β * gffMultiWickMonomial d N a mass ha hmass β ω)
          ∂(latticeGaussianMeasure d N a mass ha hmass)) =
        ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          ∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) m,
            (coeff n f α * coeff m g β) *
              (if α = β then diagFac α else 0) := by
    refine Finset.sum_congr rfl ?_
    intro α hα
    rw [MeasureTheory.integral_finsetSum _ (fun β hβ => h_summand_int α β)]
    refine Finset.sum_congr rfl ?_
    intro β hβ
    rw [show
        (fun ω =>
          (coeff n f α * gffMultiWickMonomial d N a mass ha hmass α ω) *
            (coeff m g β * gffMultiWickMonomial d N a mass ha hmass β ω)) =
          (fun ω =>
            (coeff n f α * coeff m g β) *
              (gffMultiWickMonomial d N a mass ha hmass α ω *
                gffMultiWickMonomial d N a mass ha hmass β ω)) from by
        ext ω
        ring]
    rw [MeasureTheory.integral_const_mul,
      gffMultiWickMonomial_orthogonality d N a mass ha hmass]
    simp [diagFac]
  rw [h_outer]
  by_cases hnm : n = m
  · subst hnm
    rw [if_pos rfl]
    have h_collapse :
        (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          ∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
            (coeff n f α * coeff n g β) * (if α = β then diagFac α else 0)) =
        ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          (coeff n f α * coeff n g α) * diagFac α := by
      refine Finset.sum_congr rfl ?_
      intro α hα
      rw [Finset.sum_eq_single_of_mem α hα]
      · simp [diagFac]
      · intro β hβ hβα
        have hneq : α ≠ β := by
          intro h
          exact hβα h.symm
        simp [hneq, diagFac]
    have h_diag_term :
        ∀ α, α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n →
          (coeff n f α * coeff n g α) * diagFac α =
            (n.factorial : ℝ) *
              (((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
                ∏ j, (gffSmearedCoeff d N a mass f j *
                  gffSmearedCoeff d N a mass g j) ^ α j) := by
      intro α hα
      have hprod_mul :
          (∏ j, gffSmearedCoeff d N a mass f j ^ α j) *
              (∏ j, gffSmearedCoeff d N a mass g j ^ α j) =
            ∏ j, (gffSmearedCoeff d N a mass f j *
              gffSmearedCoeff d N a mass g j) ^ α j := by
        rw [← Finset.prod_mul_distrib]
        refine Finset.prod_congr rfl ?_
        intro j hj
        rw [mul_pow]
      have hprod_pos : (0 : ℝ) < ∏ j, ((α j).factorial : ℝ) := by
        apply Finset.prod_pos
        intro j hj
        exact_mod_cast Nat.factorial_pos (α j)
      have hprod_ne : (∏ j, ((α j).factorial : ℝ)) ≠ 0 := ne_of_gt hprod_pos
      calc
        (coeff n f α * coeff n g α) * diagFac α
            = (n.factorial : ℝ) *
                (((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
                  ((∏ j, gffSmearedCoeff d N a mass f j ^ α j) *
                    (∏ j, gffSmearedCoeff d N a mass g j ^ α j))) := by
              unfold coeff diagFac
              field_simp [hprod_ne]
        _ = (n.factorial : ℝ) *
              (((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
                ∏ j, (gffSmearedCoeff d N a mass f j *
                  gffSmearedCoeff d N a mass g j) ^ α j) := by
              rw [hprod_mul]
    calc
      (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          ∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
            (coeff n f α * coeff n g β) * (if α = β then diagFac α else 0))
          =
        ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          (coeff n f α * coeff n g α) * diagFac α := h_collapse
      _ =
        ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
          (n.factorial : ℝ) *
            (((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
              ∏ j, (gffSmearedCoeff d N a mass f j *
                gffSmearedCoeff d N a mass g j) ^ α j) := by
          refine Finset.sum_congr rfl ?_
          intro α hα
          exact h_diag_term α hα
      _ = (n.factorial : ℝ) *
          (∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
            ((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
              ∏ j, (gffSmearedCoeff d N a mass f j *
                gffSmearedCoeff d N a mass g j) ^ α j) := by
          rw [← Finset.mul_sum]
      _ = (n.factorial : ℝ) *
          (∑ j, gffSmearedCoeff d N a mass f j *
            gffSmearedCoeff d N a mass g j) ^ n := by
          rw [sum_pow_eq_sum_multiIndices d N
            (fun j => gffSmearedCoeff d N a mass f j *
              gffSmearedCoeff d N a mass g j) n]
      _ = (n.factorial : ℝ) * (gffSmearedCovariance d N a mass f g) ^ n := by
          rw [gffSmearedCovariance]
  · rw [if_neg hnm]
    apply Finset.sum_eq_zero
    intro α hα
    apply Finset.sum_eq_zero
    intro β hβ
    by_cases hαβ : α = β
    · subst hαβ
      unfold multiIndicesOfTotalDegree at hα hβ
      simp only [Finset.mem_filter] at hα hβ
      exact absurd (hα.2.symm.trans hβ.2) hnm
    · simp [hαβ, diagFac]

/-- **Integrability of a product of two smeared Wick monomials** under the lattice GFF. The linearity
prerequisite for evaluating `∫ :φ(f)⁴: · wickPolynomial P c (φ_z)` term by term (pulling the integral
through the interaction polynomial's finite sum of Wick monomials). Proved by expanding both Wick
monomials into eigenbasis multi-index sums (`wickMonomial_at_smeared_eq_eigen_sum`) and using the
per-pair integrability `integrable_gffMultiWickMonomial_mul`. -/
theorem integrable_wickMonomial_smeared_mul
    (a mass : ℝ) (ha : 0 < a) (hmass : 0 < mass)
    (n m : ℕ) (f g : FinLatticeField d N) :
    MeasureTheory.Integrable
      (fun ω => wickMonomial n (∑ j, (gffSmearedCoeff d N a mass f j) ^ 2) (ω f) *
        wickMonomial m (∑ j, (gffSmearedCoeff d N a mass g j) ^ 2) (ω g))
      (latticeGaussianMeasure d N a mass ha hmass) := by
  have hcongr : ∀ ω, wickMonomial n (∑ j, (gffSmearedCoeff d N a mass f j) ^ 2) (ω f) *
        wickMonomial m (∑ j, (gffSmearedCoeff d N a mass g j) ^ 2) (ω g) =
      ∑ α ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) n,
        ∑ β ∈ multiIndicesOfTotalDegree (FinLatticeSites d N) m,
          ((((n.factorial : ℝ) / ∏ j, ((α j).factorial : ℝ)) *
              ∏ j, gffSmearedCoeff d N a mass f j ^ α j) *
            (((m.factorial : ℝ) / ∏ j, ((β j).factorial : ℝ)) *
              ∏ j, gffSmearedCoeff d N a mass g j ^ β j)) *
            (gffMultiWickMonomial d N a mass ha hmass α ω *
              gffMultiWickMonomial d N a mass ha hmass β ω) := by
    intro ω
    rw [wickMonomial_at_smeared_eq_eigen_sum d N a mass ha hmass n f ω,
        wickMonomial_at_smeared_eq_eigen_sum d N a mass ha hmass m g ω,
        Finset.sum_mul_sum]
    refine Finset.sum_congr rfl fun α _ => Finset.sum_congr rfl fun β _ => ?_
    ring
  rw [MeasureTheory.integrable_congr (Filter.Eventually.of_forall hcongr)]
  apply MeasureTheory.integrable_finsetSum
  intro α _
  apply MeasureTheory.integrable_finsetSum
  intro β _
  exact (integrable_gffMultiWickMonomial_mul d N a mass ha hmass α β).const_mul _

end GaussianField