Tensor Network Weingarten Calculus

Haar Measure & Weingarten Calculus, numerical solution.

Open In Colab

Tensor Network Weingarten Calculus

This notebook builds fluency with the tensor diagrammatic rules for computing Haar integrals over the unitary group. Every analytical step is derived symbolically in SymPy and independently verified by Monte Carlo sampling of Haar-random unitaries via NumPy/SciPy.

Background

A Haar average over the unitary group is the integral of a polynomial built from the entries \(U_{ij}\) and their complex conjugates, taken uniformly over all unitaries of a fixed dimension. Such integrals look forbidding because the group is continuous and high-dimensional, yet they collapse to something finite and combinatorial. The Weingarten calculus is the tool that carries out this collapse.

The rule expresses an average with \(k\) copies of \(U\) and \(k\) of \(U^\dagger\) as a sum over pairs of permutations of \(k\) objects. One permutation records how the output indices of the \(U\) factors are matched to those of the \(U^\dagger\) factors, the other does the same for the input indices, and each pairing carries a rational weight given by the Weingarten function. The dimension \(D\) enters only through powers that count the closed loops formed when the index matchings are drawn as a tensor diagram.

This mirrors Wick’s theorem for Gaussian integrals, where expectations reduce to sums over ways of pairing variables. Here the pairings are permutations rather than simple contractions, and the Weingarten weights replace the Gaussian propagator. The identity permutation, the uncrossed pairing, carries the largest weight and dominates at large \(D\), while crossings are suppressed by powers of \(1/D\), which makes the leading behavior of many quantities easy to read off.

The notebook works the lowest moments by hand within this scheme, climbing from the first-order trace product through the second-order purity to \(\mathbb{E}[|\mathrm{Tr}\,U|^6]=6\), and checks each step against symbolic algebra and Monte Carlo sampling.

Mathematical background

Weingarten averages have a diagrammatic form that makes the bookkeeping systematic. Each \(U\) and \(U^\dagger\) is a vertex, the average pairs every \(U\) with a \(U^\dagger\) through a permutation \(\sigma\) of the row indices and \(\tau\) of the column indices, and each closed loop of contracted indices contributes a factor \(D\). The general moment reads

\[ \int U^{\otimes t}\, X\, U^{\dagger\otimes t}\, d\mu = \sum_{\sigma,\tau \in S_t} \mathrm{Wg}(\sigma\tau^{-1}, D)\;\; \sigma\, X\, \tau^{-1}, \]

with \(\mathrm{Wg}\) the inverse of the \(t!\times t!\) Gram matrix of permutation overlaps. At large \(D\) the Weingarten function behaves as \(\mathrm{Wg}(\sigma, D) \sim D^{-t-|\sigma|}\), where \(|\sigma|\) is the minimal number of transpositions in \(\sigma\), so non-identity pairings are suppressed by powers of \(1/D\). The leading contributions are the planar (non-crossing) contractions, and each crossing costs a factor \(1/D\). This is the genus expansion that organizes matrix models and their link to two-dimensional surfaces. The third trace moment evaluated in the notebook is the smallest case in which a non-planar diagram appears, so it is a clean place to watch the \(1/D\) hierarchy at work.

1. Physical Motivation

Many quantities in quantum information, purities, frame potentials, OTOCs, channel fidelities, reduce to polynomial averages over Haar-random unitaries. The Weingarten calculus is a systematic tool for computing these averages. It is the unitary-group analogue of Wick’s theorem for Gaussian integrals: just as Gaussian expectations reduce to sums over pair contractions, Haar expectations reduce to sums over permutation pairings, weighted by the Weingarten function.

The tensor network (Penrose) representation makes these sums geometrically transparent. The goal of this notebook is to master the three-step recipe: 1. Draw the unevaluated tensor network 2. Replace \(U/U^*\) pairs with permutation pairings 3. Count closed loops (each loop = factor of \(D\))

2. Mathematical Preliminaries (Self-Contained)

2.1 The symmetric group \(S_k\)

The symmetric group \(S_k\) is the group of all permutations of \(k\) objects \(\{1, 2, \ldots, k\}\). It has \(k!\) elements.

  • \(S_1 = \{\mathrm{id}\}\): one element (the identity).
  • \(S_2 = \{\mathrm{id}, (12)\}\): two elements, the identity and the swap of objects 1 and 2.
  • \(S_3\): six elements, \(\mathrm{id}\), three transpositions \((12), (13), (23)\), and two 3-cycles \((123), (132)\).

A permutation \(\pi \in S_k\) maps each index \(a \in \{1, \ldots, k\}\) to \(\pi(a)\). The cycle structure of \(\pi\) decomposes it into disjoint cycles. For example, \((123)\) means \(1 \to 2 \to 3 \to 1\), a single cycle of length 3. The number of cycles determines how many independent index sums survive after applying Kronecker deltas.

2.2 The Weingarten function

The Weingarten function \(\mathrm{Wg}(\pi, D)\) is a real-valued function on \(S_k\) (depending on Hilbert space dimension \(D\)) that serves as the weight for each permutation pairing in the Haar integral. It is defined implicitly as the inverse of the Gram matrix: \[\sum_{\sigma \in S_k} \mathrm{Wg}(\sigma, D) \cdot D^{\#\mathrm{cycles}(\pi \sigma^{-1})} = \delta_{\pi, \mathrm{id}}.\]

Key properties: - \(\mathrm{Wg}(\pi, D)\) depends only on the cycle type (conjugacy class) of \(\pi\) - The identity permutation has the largest weight, it is the “planar” (uncrossed) pairing - Non-identity permutations (“crossings”) are suppressed by powers of \(1/D\)

2.3 The Weingarten formula

The \(k\)-th moment of the Haar measure over \(U(D)\) is: \[\boxed{\mathbb{E}_{\mathrm{Haar}}[U_{i_1 j_1} \cdots U_{i_k j_k} U^*_{\ell_1 m_1} \cdots U^*_{\ell_k m_k}] = \sum_{\pi,\sigma \in S_k} \mathrm{Wg}(\pi^{-1}\sigma, D) \prod_{a=1}^k \delta_{i_a \ell_{\pi(a)}} \delta_{j_a m_{\sigma(a)}}}\]

Interpretation: The permutation \(\pi\) controls how the output (left) indices of \(U\) are paired with those of \(U^*\). The permutation \(\sigma\) does the same for input (right) indices. Each pairing is weighted by \(\mathrm{Wg}(\pi^{-1}\sigma, D)\).

2.4 Explicit Weingarten values

\(k\) Permutation \(\tau\) Cycle type \(\mathrm{Wg}(\tau, D)\)
1 \(\mathrm{id}\) \((1)\) \(1/D\)
2 \(\mathrm{id}\) \((1)(2)\) \(1/(D^2-1)\)
2 \((12)\) \((12)\) \(-1/[D(D^2-1)]\)
3 \(\mathrm{id}\) \((1)(2)(3)\) \((D^2-2)/[D(D^2-1)(D^2-4)]\)
3 any transposition \((2,1)\) \(-1/[(D^2-1)(D^2-4)]\)
3 any 3-cycle \((3)\) \(2/[D(D^2-1)(D^2-4)]\)

2.5 Tensor network rules

  • Box = matrix/tensor; wire = index; connected wire = contracted (summed) index
  • Closed loop (wire with no open ends) = contributes a factor of \(D\)
  • Cup (semicircular arc connecting two wires) = Kronecker delta \(\delta_{ij}\)
Code 
# Imports used throughout
import numpy as np
from scipy.stats import unitary_group
import sympy as sp
from sympy import symbols, Matrix, trace as sp_trace, simplify, expand, Symbol
from itertools import permutations
from collections import defaultdict, Counter

# === Utility functions for permutation algebra ===

def compose(p1, p2):
    """Compose two permutations (as tuples): (p1 ∘ p2)(i) = p1(p2(i))."""
    return tuple(p1[p2[i]] for i in range(len(p1)))

def inverse(p):
    """Inverse of a permutation (as tuple)."""
    r = [0]*len(p)
    for i, v in enumerate(p): r[v] = i
    return tuple(r)

def cycle_type(p):
    """Return the cycle type of permutation p as a sorted tuple."""
    visited = set(); lengths = []
    for s in range(len(p)):
        if s not in visited:
            l = 0; k = s
            while k not in visited:
                visited.add(k); k = p[k]; l += 1
            lengths.append(l)
    return tuple(sorted(lengths, reverse=True))

def num_cycles(p):
    """Count the number of cycles in permutation p."""
    visited = set(); c = 0
    for s in range(len(p)):
        if s not in visited:
            c += 1; k = s
            while k not in visited:
                visited.add(k); k = p[k]
    return c

# Demonstrate S_2 and S_3
print('=== The symmetric group S_k ===')
print()
S2 = list(permutations(range(2)))
print(f'S_2 has {len(S2)} elements:')
for p in S2:
    print(f'  {p} - cycle type {cycle_type(p)}, {num_cycles(p)} cycle(s)')

print()
S3 = list(permutations(range(3)))
print(f'S_3 has {len(S3)} elements:')
for p in S3:
    print(f'  {p} - cycle type {cycle_type(p)}, {num_cycles(p)} cycle(s)')
=== The symmetric group S_k ===

S_2 has 2 elements:
  (0, 1) - cycle type (1, 1), 2 cycle(s)
  (1, 0) - cycle type (2,), 1 cycle(s)

S_3 has 6 elements:
  (0, 1, 2) - cycle type (1, 1, 1), 3 cycle(s)
  (0, 2, 1) - cycle type (2, 1), 2 cycle(s)
  (1, 0, 2) - cycle type (2, 1), 2 cycle(s)
  (1, 2, 0) - cycle type (3,), 1 cycle(s)
  (2, 0, 1) - cycle type (3,), 1 cycle(s)
  (2, 1, 0) - cycle type (2, 1), 2 cycle(s)

3. Part (a): First Moment, \(\mathbb{E}[\mathrm{Tr}(UAU^\dagger B)]\)

Analytical derivation

Goal: Compute \(\mathbb{E}_U[\mathrm{Tr}(\hat{U}\hat{A}\hat{U}^\dagger\hat{B})]\) for fixed matrices \(\hat{A}, \hat{B}\).

Step 1, Index expansion: \[\mathrm{Tr}(UAU^\dagger B) = \sum_{i,j,k,l} U_{ij}\, A_{jk}\, U^*_{lk}\, B_{li}\]

Step 2, Apply \(k=1\) Weingarten (\(|S_1| = 1\), only \(\pi = \sigma = \mathrm{id}\)): \[\mathbb{E}[U_{ij}U^*_{lk}] = \frac{1}{D}\,\delta_{il}\,\delta_{jk}\]

Step 3, Contract deltas: - \(\delta_{jk}\) closes a loop through \(\hat{A}\): \(\sum_j A_{jj} = \mathrm{Tr}(\hat{A})\) - \(\delta_{il}\) closes a loop through \(\hat{B}\): \(\sum_i B_{ii} = \mathrm{Tr}(\hat{B})\)

Result: \(\mathbb{E}[\mathrm{Tr}(UAU^\dagger B)] = \frac{\mathrm{Tr}(\hat{A})\,\mathrm{Tr}(\hat{B})}{D}\)

Code 
print('=== Part (a): SymPy symbolic verification ===')
print()

# Step 1: Build symbolic matrices (D=3 for concreteness)
D_val = 3
A = Matrix([[symbols(f'a{i}{j}') for j in range(D_val)] for i in range(D_val)])
B = Matrix([[symbols(f'b{i}{j}') for j in range(D_val)] for i in range(D_val)])

# Step 2: Apply Weingarten formula by brute-force index contraction
print('Step 2: Applying E[U_{ij}U*_{lk}] = δ_{il}δ_{jk}/D index by index...')
result_indices = sp.Integer(0)
for i in range(D_val):
    for j in range(D_val):
        for k in range(D_val):
            for l in range(D_val):
                if i == l and j == k:  # δ_{il} * δ_{jk}
                    result_indices += A[j,k] * B[l,i] / D_val

# Step 3: Compare with closed-form
result_formula = sp_trace(A) * sp_trace(B) / D_val

print(f'  Index contraction: {simplify(result_indices)}')
print(f'  Formula Tr(A)Tr(B)/D: {simplify(result_formula)}')
assert expand(result_indices - result_formula) == 0
print('  MATCH')
=== Part (a): SymPy symbolic verification ===

Step 2: Applying E[U_{ij}U*_{lk}] = δ_{il}δ_{jk}/D index by index...
  Index contraction: a00*b00/3 + a00*b11/3 + a00*b22/3 + a11*b00/3 + a11*b11/3 + a11*b22/3 + a22*b00/3 + a22*b11/3 + a22*b22/3
  Formula Tr(A)Tr(B)/D: (a00 + a11 + a22)*(b00 + b11 + b22)/3
  MATCH
Code 
print('=== Part (a): NumPy Monte Carlo cross-check ===')
np.random.seed(42)

for D_val in [2, 4, 8]:
    # Random complex matrices
    A_n = np.random.randn(D_val, D_val) + 1j*np.random.randn(D_val, D_val)
    B_n = np.random.randn(D_val, D_val) + 1j*np.random.randn(D_val, D_val)

    # Monte Carlo: sample Haar-random U, compute Tr(UAU†B)
    mc_vals = [np.trace(U @ A_n @ U.conj().T @ B_n)
               for U in [unitary_group.rvs(D_val) for _ in range(15000)]]
    mc_mean = np.mean(mc_vals)

    # Analytical
    exact = np.trace(A_n) * np.trace(B_n) / D_val

    print(f'  D={D_val}: MC = {mc_mean:.4f}, analytic = {exact:.4f}, '
          f'|err| = {abs(mc_mean - exact):.4f}')
    assert abs(mc_mean - exact) < 0.3 * max(abs(exact), 1.0)

print('  All first-moment checks PASS')
=== Part (a): NumPy Monte Carlo cross-check ===
  D=2: MC = -0.8851-0.5541j, analytic = -0.8590-0.5729j, |err| = 0.0321
  D=4: MC = -0.3010+0.1774j, analytic = -0.3099+0.1877j, |err| = 0.0136
  D=8: MC = 0.2752-0.0297j, analytic = 0.2306+0.0000j, |err| = 0.0536
  All first-moment checks PASS

4. Part (b): Second Moment, Haar-Averaged Purity

Analytical derivation

Goal: Derive \(\mathbb{E}[\mathrm{Tr}(\hat{\rho}_A^2)]\) for a Haar-random pure state \(|\psi\rangle = U|0\rangle\) on \(\mathcal{H}_A \otimes \mathcal{H}_B\) with \(\dim\mathcal{H}_A = D_A\), \(\dim\mathcal{H}_B = D_B\), \(D = D_A D_B\).

Step 1, Setup: Using the SWAP trick, \(\mathrm{Tr}(\hat{\rho}_A^2) = \mathrm{Tr}(\mathbb{F}_A \cdot \hat{\rho}^{\otimes 2})\). This requires \(k=2\) copies of \(U\) and \(U^*\).

Step 2, Index structure: All input indices are fixed to the reference state \(|0\rangle\), meaning \(j_1 = j_2 = m_1 = m_2 = 0\). This makes \(\delta_{j_a, m_{\sigma(a)}} = \delta_{0,0} = 1\) for every \(\sigma \in S_2\). Therefore the \(\sigma\) sum is free: every permutation \(\sigma\) contributes equally.

Step 3, Sigma sum: \[\sum_{\sigma \in S_2} \mathrm{Wg}(\pi^{-1}\sigma, D) = \sum_{\tau \in S_2} \mathrm{Wg}(\tau, D) = \frac{1}{D^2-1} - \frac{1}{D(D^2-1)} = \frac{1}{D(D+1)}\] This is the same for every \(\pi\) (since we’re summing \(\tau = \pi^{-1}\sigma\) over all of \(S_2\)).

Step 4, Loop counting: The output wiring has SWAP on \(A\) and identity on \(B\). The effective permutation for output pairing \(\pi\) is:

\(\pi\) Eff. on \(A\): \(\pi \circ \mathrm{SWAP}\) Eff. on \(B\): \(\pi \circ \mathrm{id}\) Loop factor
\(\mathrm{id}\) \((12)\): 1 cycle → \(D_A^1\) \(\mathrm{id}\): 2 cycles → \(D_B^2\) \(D_A D_B^2\)
\((12)\) \(\mathrm{id}\): 2 cycles → \(D_A^2\) \((12)\): 1 cycle → \(D_B^1\) \(D_A^2 D_B\)

Step 5, Assembly: \[\mathbb{E}[\mathrm{Tr}(\hat{\rho}_A^2)] = \frac{1}{D(D+1)}(D_A D_B^2 + D_A^2 D_B) = \frac{D_A D_B(D_A + D_B)}{D_A D_B(D_A D_B + 1)} = \frac{D_A + D_B}{D_A D_B + 1}\]

Code 
print('=== Part (b): Step-by-step purity derivation ===')
print()

# Step 1: Define permutation composition
perm_id = {0:0, 1:1}
perm_12 = {0:1, 1:0}

def compose_d(p1, p2): return {k: p1[p2[k]] for k in p2}
def ncyc_d(p):
    visited, c = set(), 0
    for s in p:
        if s not in visited:
            c += 1; k = s
            while k not in visited: visited.add(k); k = p[k]
    return c

# Step 2: Show the sigma sum
print('Step 2: Sigma sum (free because reference state fixes all inputs)')
D_sym = Symbol('D', positive=True)
wg_sum = 1/(D_sym**2 - 1) + (-1)/(D_sym*(D_sym**2 - 1))
print(f'  Wg(id) + Wg((12)) = {simplify(wg_sum)} = 1/[D(D+1)]')
assert simplify(wg_sum - 1/(D_sym*(D_sym+1))) == 0
print('  Confirmed')
print()

# Step 3: Loop counting table
ext_A = perm_12  # SWAP
ext_B = perm_id  # identity

print('Step 3: Cycle counting table')
print('     π | Eff on A (π∘SWAP) | Eff on B (π∘id) | Loop factor')
print('-' * 60)
total_loops = 0
for name, pi in [('id', perm_id), ('(12)', perm_12)]:
    eA = compose_d(pi, ext_A)
    eB = compose_d(pi, ext_B)
    cA, cB = ncyc_d(eA), ncyc_d(eB)
    print(f'{name:>6} | {cA} cycle(s) → D_A^{cA}    | {cB} cycle(s) → D_B^{cB}  | D_A^{cA}·D_B^{cB}')
print()

# Step 4: Numerical cross-check
print('Step 4: Numerical verification')
for dA, dB in [(2,2), (2,4), (3,3), (4,4), (2,8)]:
    D = dA * dB
    wg_s = 1.0/(D*(D+1))  # sum_tau Wg(tau)
    loops = dA*dB**2 + dA**2*dB
    purity = wg_s * loops
    exact = (dA+dB) / (dA*dB + 1.0)
    print(f'  DA={dA},DB={dB}: analytic = {purity:.8f}, formula = {exact:.8f}, '
          f'match = {abs(purity-exact) < 1e-12}')
    assert abs(purity - exact) < 1e-12
=== Part (b): Step-by-step purity derivation ===

Step 2: Sigma sum (free because reference state fixes all inputs)
  Wg(id) + Wg((12)) = 1/(D*(D + 1)) = 1/[D(D+1)]
  Confirmed

Step 3: Cycle counting table
     π | Eff on A (π∘SWAP) | Eff on B (π∘id) | Loop factor
------------------------------------------------------------
    id | 1 cycle(s) → D_A^1    | 2 cycle(s) → D_B^2  | D_A^1·D_B^2
  (12) | 2 cycle(s) → D_A^2    | 1 cycle(s) → D_B^1  | D_A^2·D_B^1

Step 4: Numerical verification
  DA=2,DB=2: analytic = 0.80000000, formula = 0.80000000, match = True
  DA=2,DB=4: analytic = 0.66666667, formula = 0.66666667, match = True
  DA=3,DB=3: analytic = 0.60000000, formula = 0.60000000, match = True
  DA=4,DB=4: analytic = 0.47058824, formula = 0.47058824, match = True
  DA=2,DB=8: analytic = 0.58823529, formula = 0.58823529, match = True
Code 
print('=== Part (b): NumPy Monte Carlo cross-check ===')
np.random.seed(42)

for dA, dB in [(2,2), (2,4), (3,3)]:
    D = dA * dB
    purities = []
    for _ in range(10000):
        psi = unitary_group.rvs(D)[:, 0]  # Haar-random |ψ⟩ = U|0⟩
        rho_A = psi.reshape(dA, dB) @ psi.reshape(dA, dB).conj().T
        purities.append(np.real(np.trace(rho_A @ rho_A)))
    mc = np.mean(purities)
    mc_err = np.std(purities) / np.sqrt(len(purities))
    exact = (dA + dB) / (dA*dB + 1.0)
    print(f'  DA={dA},DB={dB}: MC = {mc:.4f} ± {mc_err:.4f}, '
          f'analytic = {exact:.4f}, |err| = {abs(mc - exact):.4f}')
    assert abs(mc - exact) < 3 * mc_err + 0.01  # within ~3σ

print('  All purity checks PASS')
=== Part (b): NumPy Monte Carlo cross-check ===
  DA=2,DB=2: MC = 0.8004 ± 0.0013, analytic = 0.8000, |err| = 0.0004
  DA=2,DB=4: MC = 0.6664 ± 0.0010, analytic = 0.6667, |err| = 0.0003
  DA=3,DB=3: MC = 0.6004 ± 0.0010, analytic = 0.6000, |err| = 0.0004
  All purity checks PASS

5. Part (c): Third Moment, \(\mathbb{E}[|\mathrm{Tr}(U)|^6]\) via Full \(S_3\) Calculus

Analytical derivation

Goal: Compute \(\mathbb{E}[|\mathrm{Tr}(U)|^6]\) for Haar-random \(U \in U(D)\).

Step 1, Index expansion: \(|\mathrm{Tr}(U)|^6 = \sum_{i_1,i_2,i_3,\ell_1,\ell_2,\ell_3} U_{i_1 i_1} U_{i_2 i_2} U_{i_3 i_3} U^*_{\ell_1 \ell_1} U^*_{\ell_2 \ell_2} U^*_{\ell_3 \ell_3}\).

This sets \(j_a = i_a\) and \(m_a = \ell_a\) (the trace constraint: all indices are “diagonal”).

Step 2, Weingarten sum: For each \((\pi, \sigma) \in S_3 \times S_3\), the delta product \(\prod_a \delta_{i_a, \ell_{\pi(a)}} \delta_{i_a, \ell_{\sigma(a)}}\) forces \(\ell_{\pi(a)} = \ell_{\sigma(a)}\) for all \(a\). The number of free summation indices equals \(D^{\#\text{cycles}(\pi^{-1}\sigma)}\).

Step 3, Group by cycle type of \(\pi^{-1}\sigma\): Since \(\mathrm{Wg}\) depends only on cycle type, we count how many of the \(36\) pairs \((\pi, \sigma)\) give each cycle type of \(\pi^{-1}\sigma\).

Code 
print('=== Part (c): Full S₃ Weingarten calculus ===')
print()

# Step 1: Enumerate S_3 and show cycle types
S3 = list(permutations(range(3)))
print(f'Step 1: S₃ has {len(S3)} elements, giving {len(S3)**2} pairs (π,σ)')
print()
ct_count = Counter(cycle_type(p) for p in S3)
for ct, n in sorted(ct_count.items()):
    print(f'  Cycle type {ct}: {n} permutations')

# Step 2: Count (π,σ) pairs by cycle type of π^{-1}σ
print()
print('Step 2: Distribution of π⁻¹σ cycle types across all 36 pairs')

pair_types = defaultdict(lambda: {'count': 0, 'ncyc': 0})
for pi in S3:
    for sig in S3:
        prod = compose(inverse(pi), sig)
        ct = cycle_type(prod)
        pair_types[ct]['count'] += 1
        pair_types[ct]['ncyc'] = num_cycles(prod)

print('  cycle type | # pairs | # cycles | loop factor')
print('-' * 50)
for ct in sorted(pair_types.keys()):
    d = pair_types[ct]
    cnt = d["count"]
    ncyc = d["ncyc"]
    print(f'{str(ct):>12} | {cnt:7d} | {ncyc:8d} | D^{ncyc}')
=== Part (c): Full S₃ Weingarten calculus ===

Step 1: S₃ has 6 elements, giving 36 pairs (π,σ)

  Cycle type (1, 1, 1): 1 permutations
  Cycle type (2, 1): 3 permutations
  Cycle type (3,): 2 permutations

Step 2: Distribution of π⁻¹σ cycle types across all 36 pairs
  cycle type | # pairs | # cycles | loop factor
--------------------------------------------------
   (1, 1, 1) |       6 |        3 | D^3
      (2, 1) |      18 |        2 | D^2
        (3,) |      12 |        1 | D^1
Code 
# Step 3: Compute E[|Tr(U)|^6] for several D
print('Step 3: Explicit evaluation of the Weingarten sum')
print()

def wg3_values(D):
    """Return k=3 Weingarten function values by cycle type."""
    dn = D * (D**2 - 1) * (D**2 - 4)
    return {
        (1,1,1): (D**2 - 2) / dn,
        (2,1): -1.0 / ((D**2 - 1) * (D**2 - 4)),
        (3,): 2.0 / dn
    }

print('   D | ', end='')
for ct in sorted(pair_types.keys()):
    n = pair_types[ct]['count']
    nc = pair_types[ct]['ncyc']
    ct_str = str(ct)
    print(f'{n} x Wg({ct_str}) x D^{nc}', end=' | ')
print('Total')
print('-' * 80)

for D in [3, 4, 5, 10, 50, 200]:
    W = wg3_values(D)
    terms = []
    total = 0.0
    for ct in sorted(pair_types.keys()):
        count = pair_types[ct]['count']
        nc = pair_types[ct]['ncyc']
        contrib = count * W[ct] * D**nc
        terms.append(f'{contrib:+10.4f}')
        total += contrib
    print(f'{D:4d} | ' + ' | '.join(terms) + f' | {total:6.1f}')

print()
print('Result: E[|Tr(U)|⁶] = 6 = 3! for all D ≥ 3')
print('This is the CUE moment theorem: E[|Tr(U)|^{2k}] = k!')
Step 3: Explicit evaluation of the Weingarten sum

   D | 6 x Wg((1, 1, 1)) x D^3 | 18 x Wg((2, 1)) x D^2 | 12 x Wg((3,)) x D^1 | Total
--------------------------------------------------------------------------------
   3 |    +9.4500 |    -4.0500 |    +0.6000 |    6.0
   4 |    +7.4667 |    -1.6000 |    +0.1333 |    6.0
   5 |    +6.8452 |    -0.8929 |    +0.0476 |    6.0
  10 |    +6.1869 |    -0.1894 |    +0.0025 |    6.0
  50 |    +6.0072 |    -0.0072 |    +0.0000 |    6.0
 200 |    +6.0005 |    -0.0005 |    +0.0000 |    6.0

Result: E[|Tr(U)|⁶] = 6 = 3! for all D ≥ 3
This is the CUE moment theorem: E[|Tr(U)|^{2k}] = k!
Code 
# Step 4: NumPy Monte Carlo cross-check
print('Step 4: NumPy Monte Carlo cross-check')
np.random.seed(42)

for D in [3, 5, 10]:
    vals = [abs(np.trace(unitary_group.rvs(D)))**6 for _ in range(20000)]
    mc = np.mean(vals)
    mc_err = np.std(vals) / np.sqrt(len(vals))
    print(f'  D={D:3d}: MC = {mc:.3f} ± {mc_err:.3f}, analytic = 6.000, '
          f'|err| = {abs(mc - 6):.3f}')
    assert abs(mc - 6.0) < 3*mc_err + 0.2

print('  All sixth-moment checks PASS')
Step 4: NumPy Monte Carlo cross-check
  D=  3: MC = 6.155 ± 0.159, analytic = 6.000, |err| = 0.155
  D=  5: MC = 6.331 ± 0.209, analytic = 6.000, |err| = 0.331
  D= 10: MC = 6.089 ± 0.178, analytic = 6.000, |err| = 0.089
  All sixth-moment checks PASS

6. Physical Interpretation

The Gaussian limit of the trace moments

This result connects Haar-random unitaries to complex Gaussian random variables. For a single complex Gaussian \(z \sim \mathcal{N}_\mathbb{C}(0,1)\), the moments are \(\mathbb{E}[|z|^{2k}] = k!\) (by Wick’s theorem). In the large-\(D\) limit, \(\mathrm{Tr}(U)\) for a Haar-random \(U \in U(D)\) behaves like a standard complex Gaussian, a result of Diaconis and Shahshahani that holds despite the strong correlations between the eigenphases \(e^{i\theta_j}\) of \(U\).

The remarkable fact is that the Weingarten calculus gives \(k!\) exactly for all \(D \geq k\), not just asymptotically. This is because the Weingarten function is the exact inverse of the representation-theoretic Gram matrix, and the exact cancellations between planar and non-planar contributions produce dimension-independent trace moments.

The topological structure

  1. Planar pairings (identity permutation, uncrossed wires) produce the most loops and dominate at large \(D\)
  2. Non-planar pairings (crossings) produce fewer loops and are suppressed by \(1/D\) per crossing
  3. The \(1/D\) suppression hierarchy mirrors the genus expansion of matrix models, a deep connection between random matrix theory and two-dimensional quantum gravity

Takeaway

Each Haar average collapses to a weighted sum over permutations, with loop counting fixing the dimension factors and the Weingarten function fixing the weights, reproducing the closed forms through \(\mathbb{E}[|\mathrm{Tr}\,U|^6] = 6\). Symbolic algebra and Monte Carlo sampling agree at every moment.

References

  • Collins and Nechita, Random quantum channels I: graphical calculus and the Bell state phenomenon, Communications in Mathematical Physics (2010).
  • Collins et al., The Weingarten calculus, Notices of the American Mathematical Society (2022).
  • Fukuda et al., RTNI, A symbolic integrator for Haar-random tensor networks, Journal of Physics A: Mathematical and Theoretical (2019). arXiv:1902.08539
  • Puchała and Miszczak, Symbolic integration with respect to the Haar measure on the unitary group, Bulletin of the Polish Academy of Sciences: Technical Sciences (2017). arXiv:1109.4244
  • Płodzień, Information Scrambling, Quantum Chaos, and Haar-Random States, lecture notes (2025). arXiv:2511.14397