---
title: "Lefschetz Thimbles, MCMC, Resummation, and Spinfoam Numerics"
subtitle: "A half-day numerical workshop"
author: "Hongguang Liu"
format:
  revealjs:
    theme: [default]
    css: workshop-slides.css
    slide-number: false
    chalkboard: true
    code-fold: false
    code-overflow: wrap
    toc: false
execute:
  enabled: true
  eval: true
  echo: true
  freeze: auto
  warning: false
  message: false
---

## Why This Workshop Needs Numerics

Spinfoam amplitudes are not just formal objects.  In practice we need to evaluate high-dimensional sums and integrals with phases, constraints, gauge redundancies, and saddle-point structure.

The workshop goal is to build a usable numerical mental model:

- How oscillatory integrals become sampling problems.
- How complex critical points enter actual estimates.
- How MCMC, resummation, and summation acceleration fit into the same workflow.

## Structure

1. Probability measures, Monte Carlo, and MCMC.
2. Spinfoam numerics and the sign problem.
3. Complex Gaussian warm-up: contour deformation without mystery.
4. Picard-Lefschetz theory at the level needed for computation.
5. Finite-flow thimble algorithms.
6. Airy integral: complete executable example.
7. Diagnostics on flowed manifolds.
8. Advanced methods: Worldvolume TLTM and intersection numbers.
9. Resummation, resurgence, and spinfoam-oriented examples.

## Probability Measure

Monte Carlo needs a probability measure:

$$
d\mu(x)=p(x)\,dx,\qquad
p(x)\ge0,\qquad
\int p(x)\,dx=1.
$$

Expectations are then

$$
\mathbb E_\mu[O]
=
\int O(x)\,d\mu(x)
=
\int O(x)p(x)\,dx .
$$

With independent samples $x_k\sim \mu$,

$$
\widehat O_N=\frac1N\sum_{k=1}^{N}O(x_k),
\qquad
\operatorname{Var}(\widehat O_N)=\frac{\operatorname{Var}_\mu(O)}{N}.
$$

This is the basic reason we insist on a positive normalized measure.

## Why MCMC Is Needed

In physics we often know an unnormalized density

$$
\pi(x)=\frac{1}{Z}e^{-S(x)},
\qquad
Z=\int e^{-S(x)}dx,
$$

but not the normalization $Z$.

MCMC constructs correlated samples from $\pi$ using only ratios such as

$$
\frac{\pi(y)}{\pi(x)}
=
\exp[-S(y)+S(x)].
$$

The normalization cancels.  That is why MCMC is useful before we ever
complexify the variables.

## Markov Chains and Stationarity

![](assets/metropolis-transition.svg){fig-alt="Metropolis proposal, accept, and reject transition schematic." width=72%}

A Markov chain is defined by a transition kernel

$$
P(x,dy)=\Pr(X_{n+1}\in dy\mid X_n=x).
$$

The target distribution $\pi$ is stationary if

$$
\pi(dy)=\int \pi(dx)\,P(x,dy).
$$

For a finite state space this is

$$
\pi_j=\sum_i\pi_iP_{ij},
\qquad
\pi P=\pi.
$$

If the chain is irreducible and aperiodic, samples converge to the target
distribution rather than remembering the starting point.

## Runnable Check: A Discrete Chain

```{julia}
using LinearAlgebra

P = [0.88 0.12 0.00;
     0.08 0.84 0.08;
     0.00 0.18 0.82]

eig = eigen(P')
k = argmin(abs.(eig.values .- 1))
pi_stat = real(eig.vectors[:, k])
pi_stat ./= sum(pi_stat)

(stationary = pi_stat,
 residual = pi_stat' * P - pi_stat')
```

This cell turns $\pi P=\pi$ into a number before introducing
continuous-state Metropolis sampling.

## Detailed Balance

A sufficient condition for stationarity is

$$
\pi(x)dx P(x,dy)=\pi(y)dy P(y,dx).
$$

For a proposal density $q(y\mid x)$, we can introduce anacceptance probability
$\alpha(x,y)$, s.t.,

$$
P(x,dy)=q(y\mid x)\alpha(x,y)dy
+r(x)\delta_x(dy),
$$
where
$$
r(x)=1-\int q(y\mid x)\alpha(x,y)dy
$$
is the mathematical version of reject and stay.

The detailed balance condition becomes
$$
\pi(x)q(y\mid x)\alpha(x,y) = \pi(y)q(x\mid y)\alpha(y,x).
$$
## Metropolis-Hastings Rule

Choose

$$
\alpha(x,y)
=
\min\left[
1,
\frac{\pi(y)q(x\mid y)}{\pi(x)q(y\mid x)}
\right].
$$

Then detailed balance $\pi(x)q(y\mid x)\alpha(x,y) = \pi(y)q(x\mid y)\alpha(y,x)$ holds.  

If $q(y\mid x)=q(x\mid y)$,

$$
\alpha(x,y)
=
\min\left[1,\frac{\pi(y)}{\pi(x)}\right].
$$

For $\pi(x)\propto e^{-S(x)}$,

$$
\alpha(x,y)
=
\min\left[1,\exp[-S(y)+S(x)]\right].
$$

## Runnable Metropolis Core

```{julia}
using Random, Statistics

function metropolis(logpi; nsteps=6_000, step=0.8, x0=0.0)
    x = x0
    logpx = logpi(x)
    xs = Vector{Float64}(undef, nsteps)
    accepts = 0
    for n in 1:nsteps
        y = x + step * randn()
        logpy = logpi(y)
        if log(rand()) < logpy - logpx
            x = y
            logpx = logpy
            accepts += 1
        end
        xs[n] = x
    end
    xs, accepts / nsteps
end

logpi_gaussian(x) = -0.5 * x^2
chain, acc = metropolis(logpi_gaussian; step=1.0)
(acceptance = acc, mean = mean(chain), variance = var(chain))
```

Only log-density differences enter.  This is the feature that later makes
Metropolis useful for flowed thimble targets.

## Complex Weights Are Not Measures

Now consider an oscillatory or complex weight

$$
w(x)=f(x)e^{-S(x)},\qquad S=S_R+iS_I.
$$

The signed or complex object

$$
w(x)\,dx
$$

is not a probability measure.  It is not positive, and it cannot be used
directly as a sampling distribution.

We must choose a positive sampling measure and put the complex phase into
the estimator.  This is the origin of the sign problem.

## The Basic Numerical Object

We repeatedly meet integrals of the form

$$
A = \int_{\mathcal C} d^n x\, f(x) e^{-S(x)},
\qquad S(x)=S_R(x)+iS_I(x).
$$

If $S_I$ varies rapidly, then $e^{-S}$ is not a probability density.

Direct Monte Carlo samples from some positive density $p(x)$ and estimates

$$
A = \left\langle \frac{f(x)e^{-S(x)}}{p(x)} \right\rangle_p .
$$

The estimator fails when phase cancellations dominate.

## Sign Problem as Signal-to-Noise Problem

![](assets/sign-problem-phases.svg){fig-alt="Many phase vectors cancel to a small average phase." width=82%}

A common reweighting form is 

$$
\langle O\rangle
=\frac{\langle O e^{-i\theta}\rangle_{|w|}}
{\langle e^{-i\theta}\rangle_{|w|}}.
$$

where we take $|w| = e^(-S_R)$ and $\theta = S_I$.

The denominator is the average phase.
If

$$
|\langle e^{-i\theta}\rangle_{|w|}| \sim e^{-cV},
$$

then the number of samples must grow exponentially with volume or dimension.

More explicitly, for phase samples of unit modulus,

$$
\operatorname{Var}\!\left(\overline{e^{i\theta}}\right)
\simeq \frac{1-|\langle e^{i\theta}\rangle|^2}{N}.
$$

To keep a fixed relative error in the denominator, one needs roughly

$$
N\gtrsim |\langle e^{i\theta}\rangle|^{-2}
\sim e^{2cV}.
$$


## Why Toy Models Are Necessary

Spinfoam examples are too expensive for every method to be introduced there first.

We use a ladder of toy models:

| Model | What it teaches |
|---|---|
| Complex Gaussian | exact contour shift and phase removal |
| Airy integral | saddles, Stokes sectors, finite-flow thimbles |
| Two-saddle quartic | multi-thimble weights, TLTM ergodicity, Borel singularities |
| 0D phi-four | asymptotic series and Borel-Pade |
| Double well QM | instanton/nonperturbative saddle |
| Lattice scalar | field-theory average phase decay |
| Summation boost | cutoff extrapolation and sequence acceleration |
| Multiple shooting toy | unstable upward flow and intersection-number conditioning |

## Complex Gaussian Warm-Up

Integral:

$$
I(b)=\int_{-\infty}^{\infty}dx\,\exp\left[-\frac{x^2}{2}+ibx\right]\\
= e^{-b^2/2} \int_{-\infty}^{\infty}dx\,\exp\left[-\frac{1}{2}(x - ib)^2\right] = e^{-b^2/2} \int_{-\infty - ib}^{\infty - ib}dz\,\exp\left[-\frac{z^2}{2}\right]
=\sqrt{2\pi}e^{-b^2/2}.
$$

This is the simplest exact sign-problem example and the cleanest contour-shift calculation in the deck.

## Gaussian: Runnable Core

```{julia}
using QuadGK

gaussian_exact(b) = sqrt(2 * pi) * exp(-b^2 / 2)
gaussian_reweighted(b) = quadgk(x -> exp(-x^2 / 2 + im*b*x), -Inf, Inf)[1]
gaussian_shifted(b) = quadgk(x -> exp(-((x + im*b)^2)/2 + im*b*(x + im*b)), -Inf, Inf)[1]

b = 3.0
(exact = gaussian_exact(b),
 real_axis = gaussian_reweighted(b),
 shifted = gaussian_shifted(b))
```

The shifted contour is the same integral but no longer
has the violent $e^{ibx}$ phase.

## Gaussian: Phase Cancellation Plot

```{julia}
using Plots

xs = range(-5, 5; length=500)
b = 3.0
phase = exp.(im .* b .* xs)
amp = exp.(-xs.^2 ./ 2)

plot(xs, amp .* real.(phase); label="real part", lw=3,
     xlabel="x", ylabel="weighted integrand",
     title="Oscillatory real-axis integrand")
plot!(xs, amp .* imag.(phase); label="imaginary part", lw=3)
plot!(xs, amp; label="Gaussian envelope", lw=2, ls=:dash)
```

This plot is where the sign problem stops being a slogan: the envelope is
large, but the signed area is small.

## Complex Gaussian Saddle

![](assets/gaussian-contour-shift.svg){fig-alt="Original real Gaussian contour shifted to pass through the complex saddle." width=82%}

Write

$$
S(z)=\frac{z^2}{2}-ibz.
$$

The saddle is

$$
\partial_zS=z-ib=0,
\qquad z_*=ib.
$$

Shift the contour to

$$
z=x+ib.
$$

Then

$$
S(x+ib)=\frac{x^2}{2}+\frac{b^2}{2}.
$$

The oscillation disappears exactly.

The contour equality being used is

$$
\int_{\mathbb R} dz\,e^{-S(z)}
=\int_{\mathbb R+ib} dz\,e^{-S(z)},
$$

because the integrand is entire and the vertical closing segments vanish.

## Lesson from the Gaussian

This model shows the cleanest version of the thimble idea:

- The original real contour is oscillatory.
- The deformed contour passes through the saddle.
- The action phase becomes constant.
- The integral is unchanged by contour deformation.
- The Jacobian is trivial here, but nontrivial in general.

Airy is the next step because it has multiple saddles and Stokes geometry.

## Where This Appears in Spinfoams

In EPRL-type amplitudes one sees:

- Representation sums over spins and intermediate labels.
- Group/spinor integrals with complex phases.
- Boundary data that may or may not admit real Regge-like critical points.
- Large-spin asymptotics controlled by saddle points.

The numerical methods must handle both discrete sums and complex continuous integrals.

## Three Numerical Regimes

| Regime | Typical Method | Main Difficulty |
|---|---|---|
| Small spins | exact or truncated sums | combinatorics and memory |
| Intermediate | integral representation + MCMC | sign problem and high dimension |
| Large spins | saddle/asymptotics | critical points and Stokes data |

The workshop focuses on the middle column while keeping the other two connected.

## Minimal Spinfoam Link

The spinfoam vertex amplitude in coherent-state integral form schematically resembles

$$
A_v(j_f,\xi_{ef})
=\int \prod_e dg_e \prod_{ef} dz_{ef}\; \mu(g,z)\,e^{S[j,\xi,g,z]}.
$$

For large spins $j_f=\lambda k_f$,

$$
S[j,\xi,g,z]=\lambda S[k,\xi,g,z].
$$

The large parameter makes saddle methods natural.

## Critical Points in the Large-Spin Limit

Critical points solve

$$
\partial S = 0.
$$

They encode:

- Closure conditions.
- Parallel transport conditions.
- Geometric or non-geometric boundary data.
- Real or complex discrete geometries.

The same equations also determine the local directions used in thimble sampling.

## Why Complex Critical Points Matter

Complex saddles can be relevant when:

- No real saddle exists for the chosen boundary data.
- A subdominant branch controls tunneling-like behavior.
- A Stokes phenomenon changes which saddles contribute.
- A perturbative expansion around one saddle is ambiguous and needs another saddle sector.

They are not optional decorations; they can control the answer.

## Complexification

Start from real variables $x^i\in\mathbb R$ and extend them to complex variables

$$
z^i=x^i+i y^i.
$$

Assume $S(z)$ and $f(z)$ are holomorphic in the relevant domain and singularities are controlled.

Then one can deform integration cycles without changing the integral, as long as no singularity is crossed and boundary behavior remains valid.

## Picard-Lefschetz Decomposition

![](assets/pl-flow-cycles.svg){fig-alt="Downward thimbles and upward dual cycles attached to complex saddles." width=82%}

The original cycle decomposes as

$$
\mathcal C \sim \sum_\sigma n_\sigma \mathcal J_\sigma.
$$

The integral becomes

$$
\int_{\mathcal C} dz\,e^{-S(z)}O(z)
=
\sum_\sigma n_\sigma
\int_{\mathcal J_\sigma} dz\,e^{-S(z)}O(z).
$$

Here:

- $z_\sigma$ is a critical point.
- $\mathcal J_\sigma$ is the downward thimble.
- $\mathcal K_\sigma$ is the upward dual cycle.
- $n_\sigma=\langle \mathcal C,\mathcal K_\sigma\rangle$ is an intersection number.

## Downward and Upward Flows

With holomorphic $S$, use the antiholomorphic flow

$$
\frac{dz^i}{dt}=\overline{\frac{\partial S}{\partial z^i}}.
$$

Then

$$
\frac{dS}{dt}=\sum_i \left|\frac{\partial S}{\partial z^i}\right|^2.
$$

Therefore

$$
\frac{dS_R}{dt}\ge 0,
\qquad
\frac{dS_I}{dt}=0.
$$

This is steepest ascent for $S_R$.

## Thimble Geometry

![](assets/ppt-thimble-geometry.png){fig-alt="Cropped thimble geometry figure showing steepest-flow paths approaching a critical point." width=52%}

The same equation can be written in real and imaginary components,
$z^i=z_R^i+i z_I^i$:

$$
\frac{dz_R^i}{dt}
=\frac{\partial S_R}{\partial z_R^i}
=\frac{\partial S_I}{\partial z_I^i},
\qquad
\frac{dz_I^i}{dt}
=\frac{\partial S_R}{\partial z_I^i}
=-\frac{\partial S_I}{\partial z_R^i}.
$$

So the flow is simultaneously gradient flow for $S_R$ and Hamiltonian
flow for $S_I$.  This is the geometric reason the phase is conserved on
an exact thimble.

## Which Flow Defines the Thimble?

The thimble attached to $z_\sigma$ is the set of points that flow back to $z_\sigma$ under steepest ascent, equivalently points reached from the saddle by finite positive flow from the tangent space.

The dual cycle is generated by the opposite flow direction.

For computation, it is often easier to start near the saddle and flow outward.

## Finite-Time Thimble Approximation

![](assets/ppt-thimble-finite-time.png){fig-alt="Cropped thimble figure showing a tangent-space neighborhood flowed toward a thimble." width=56%}

The exact thimble is reached only in an infinite-flow limit.  Numerically
we replace it by a finite-flow surface:

$$
V_\sigma \xrightarrow{\;\mathrm{SA\;flow\;for\;}T\;} \widetilde{\mathcal J}_{\sigma,T},
\qquad
\widetilde{\mathcal J}_{\sigma,T}\to\mathcal J_\sigma
\quad (T\to\infty).
$$

Here $V_\sigma$ is a neighborhood in the tangent space at the critical
point.  The practical approximation is controlled by both the size of
$V_\sigma$ and the flow time $T$.

## Why the Phase Is Better on a Thimble

On an exact thimble,

$$
S_I(z)=S_I(z_\sigma).
$$

Thus the original oscillating factor becomes

$$
e^{-S(z)}=e^{-S_R(z)}e^{-iS_I(z_\sigma)}.
$$

The main phase oscillation is removed.  A residual phase remains from the embedding measure.

## Local Expansion Near a Saddle

Let $z=z_\sigma+\eta$.  Then

$$
S(z)=S(z_\sigma)+\frac12 H_{ij}\eta^i\eta^j+O(\eta^3),
$$

where

$$
H_{ij}=\frac{\partial^2 S}{\partial z^i\partial z^j}\bigg|_{z_\sigma}.
$$

The tangent directions are determined by the complex symmetric Hessian.

## Linearized Flow Near a Critical Point

![](assets/ppt-linearized-directions.png){fig-alt="Cropped figure showing repulsive and attractive directions near a complex critical point." width=70%}

In a small neighborhood $V_\sigma$, write $z=z_\sigma+\omega$:

$$
S(z)=S(z_\sigma)+\frac12 H_{ij}\omega^i\omega^j+O(\omega^3).
$$

The linearized steepest-ascent equation is

$$
\frac{d\omega}{dt}=\overline{H\omega}.
$$

After splitting $\omega=\omega_R+i\omega_I$ and
$H=H_R+iH_I$, this becomes a real system

$$
\frac{d}{dt}
\begin{pmatrix}\omega_R\\ \omega_I\end{pmatrix}
=
\begin{pmatrix}
H_R & -H_I\\
-H_I & -H_R
\end{pmatrix}
\begin{pmatrix}\omega_R\\ \omega_I\end{pmatrix}.
$$

Its eigenvalues appear in pairs $(\lambda,-\lambda)$; the positive
directions are the outgoing thimble directions for steepest ascent.

## Takagi Factorization

For a complex symmetric Hessian,

$$
H\rho_a=\lambda_a \overline{\rho_a},
\qquad \lambda_a\ge0.
$$

The positive Takagi directions span the local thimble tangent space.

Near the saddle,

$$
\eta=\sum_a \rho_a x^a,
\qquad x^a\in\mathbb R.
$$

## One-Dimensional Shortcut

If there is one complex variable and Hessian $H=|H|e^{i\phi}$, a Takagi direction is

$$
\rho=e^{-i\phi/2}.
$$

Then

$$
H\rho^2=|H|>0.
$$

This is used in the Airy example below.

## Airy Model Used Below

The one-dimensional test action is

$$
S(z;a)=-i\left(\frac{z^3}{3}+az\right),
\qquad a>0.
$$

The exact real-axis integral is

$$
\int_{-\infty}^{\infty}dz\,e^{-S(z;a)}
=2\pi\operatorname{Ai}(a).
$$

The critical equation is

$$
\partial_zS=-i(z^2+a)=0,
\qquad
z_\pm=\pm i\sqrt a.
$$

This gives an analytic benchmark for thimble geometry, finite flow, MCMC,
and residual phase.

## Takagi Direction: Runnable Check

For the Airy saddle $z_\sigma=i\sqrt a$,

$$
S''(z_\sigma)=-2iz_\sigma=2\sqrt a.
$$

So the Hessian is real and positive, and the local tangent direction is
the real axis in the saddle-centered coordinate.

```{julia}
a = 1.0
zc = im * sqrt(a)
H = -2im * zc
phi = angle(H)
rho = exp(-im * phi / 2)

(zc = zc,
 H = H,
 takagi_direction = rho,
 quadratic_coefficient = H * rho^2)
```

For a general one-dimensional complex Hessian, the last number should be
real and positive.  That is the local Gaussian direction we flow from.

## Finite-Flow Map

![](assets/finite-flow-map.svg){fig-alt="Finite flow maps tangent coordinates to a curved flowed contour with Jacobian." width=82%}

Define

$$
z_0(x)=z_\sigma+\sum_a\rho_a x^a,
$$

and flow for time $T$:

$$
z_T(x)=\mathcal C_T(z_0(x)).
$$

This gives a finite-flow approximation to the thimble.

As $T$ grows, sign oscillations improve, but ODE stiffness and cost increase.

## From Tangent Space to Flowed Surface

![](assets/ppt-thimble-tangent-flow.png){fig-alt="Cropped figure showing finite flow from a tangent-space patch to a flowed thimble approximation." width=58%}

The finite-flow construction has two independent cutoffs:

$$
x\in V_\sigma\subset\mathbb R^n,
\qquad
t\in[0,T].
$$

The approximation improves when $V_\sigma$ is large enough to cover the
relevant Gaussian support and $T$ is large enough to align the phase:

$$
\int_{\mathcal J_\sigma} dz\,e^{-S(z)}O(z)
\approx
\int_{V_\sigma} dx\,
\det J_T(x)\,e^{-S(z_T(x))}O(z_T(x)).
$$

The numerical danger is that increasing $T$ also makes the target sharper
and the Jacobian flow less stable.

## Finite-Flow Replacement

The finite-flow replacement is

$$
\int_{\mathcal C} d^n z\,e^{-S(z)}O(z)
\approx
\int_{\mathbb R^n} d^n x\;
\det J_T(x)\,e^{-S(z_T(x))}O(z_T(x)).
$$

This is the formula the live Airy cells implement.

## Jacobian Flow

The Jacobian matrix is

$$
J^i{}_a(t)=\frac{\partial z^i_t}{\partial x^a}.
$$

It obeys

$$
\frac{dJ^i{}_a}{dt}
=\overline{H_{ij}(z_t)J^j{}_a},
\qquad
J^i{}_a(0)=\rho^i_a.
$$

This equation is as important as the coordinate flow.

## Jacobian Flow: Geometric Picture

![](assets/ppt-jacobian-flow-cubes.png){fig-alt="Cropped figure showing local volume elements deforming along flow lines." width=56%}

The Jacobian is not optional bookkeeping.  It is the change of variables
from the flat tangent coordinate $x$ to the curved surface:

$$
J^i{}_a(t,x)=\frac{\partial z^i_t(x)}{\partial x^a}.
$$

Differentiating the coordinate flow gives

$$
\frac{dJ^i{}_a}{dt}
=
\overline{
\frac{\partial^2 S}{\partial z^i\partial z^j}(z_t)
J^j{}_a
}.
$$

Thus one proposal in tangent space carries two flows: $z_t(x)$ and
$J_t(x)$.

## Flowed Integral

The original thimble integral becomes

$$
\int_{\mathcal J_\sigma} d^n z\,f(z)e^{-S(z)}
\approx
\int_{\mathbb R^n} d^n x\,
\det J(x) f(z_T(x))e^{-S(z_T(x))}.
$$

For multiple thimbles,

$$
A\approx \sum_\sigma n_\sigma
\int d^n x\,\det J_\sigma(x) f(z_{\sigma,T}(x))e^{-S(z_{\sigma,T}(x))}.
$$

## Effective Action

Define

$$
S_{\rm eff}(x)=S_R(z_T(x))-\log|\det J(x)|.
$$

Then

$$
p_{\rm eff}(x)\propto e^{-S_{\rm eff}(x)}
$$

is a positive target distribution.

This is the object sampled by MCMC.

## Residual Phase

The remaining phase is

$$
e^{i\theta_{\rm res}(x)}
=\frac{\det J(x)}{|\det J(x)|}
 e^{-iS_I(z_T(x))}.
$$

For an exact thimble, $S_I$ is constant, but $\arg\det J$ can still vary.

The observable estimator is

$$
\langle O\rangle
\approx
\frac{\langle O(z_T(x))e^{i\theta_{\rm res}(x)}\rangle_{\rm eff}}
{\langle e^{i\theta_{\rm res}(x)}\rangle_{\rm eff}}.
$$

## Algorithm: Single-Thimble Finite-Flow MCMC

1. Find a critical point $z_\sigma$.
2. Compute Hessian and Takagi directions $\rho_a$.
3. Choose flow time $T$ and tangent coordinate scale.
4. For each proposed $x$, solve flow and Jacobian ODEs.
5. Compute $S_{\rm eff}(x)$ and residual phase.
6. Run MCMC with target $e^{-S_{\rm eff}}$.
7. Estimate observables by residual-phase reweighting.

## Pseudocode: One MCMC Proposal

```text
input: current tangent coordinate x, current S_eff
propose x_new = x + proposal_noise
z0 = z_sigma + rho * x_new
solve dz/dt = conjugate(dS/dz)
solve dJ/dt = conjugate(H(z) J)
compute S_eff_new = Re S(z_T) - log |det J_T|
accept with min(1, exp[-S_eff_new + S_eff])
if accepted: store z_T, det J_T, residual phase
```

This is simple conceptually but expensive because every proposal contains ODE solves.

For expensive models cache or reuse:

- Critical point data.
- Hessian and Takagi basis at the saddle.
- Boundary data and gauge-fixing substitutions.
- ODE solver configuration and tolerances.
- Failed flow regions, if proposals often hit singularities.
- Diagnostic metadata: flow time, proposal scale, seed, chain id.



## Cost

Each MCMC proposal costs:

- One nonlinear ODE solve for $z_T(x)$.
- One Jacobian ODE solve, dimension $n\times n$.
- One action and Hessian evaluation along the path.

For spinfoam-like dimensions this is expensive.  This motivates adaptive samplers, approximations, and precomputed benchmarks.


For $n$ complex variables:

- The coordinate flow has $2n$ real degrees of freedom.
- The Jacobian is an $n\times n$ complex matrix.
- A naive Jacobian flow therefore adds $2n^2$ real ODE variables.

Thus a single proposal may require solving an ODE system with

$$
2n+2n^2
$$

real variables, before counting action/Hessian evaluation cost.

## Why High Dimension Is Hard

The challenge is not only $n$.

Important costs are:

- Hessian evaluation along every flow trajectory.
- Stiffness at large flow time.
- Determinant computation and phase tracking.
- Poor local proposals when the flowed distribution is anisotropic.
- Multiple separated thimbles or modes.

This is why a clean toy model pipeline is necessary before spinfoam-scale calculations.

## Practical Approximation Ladder

In increasing cost:

1. Saddle-point approximation using Hessian determinant.
2. Perturbative expansion around one saddle.
3. Finite-flow quadrature in one or two dimensions.
4. Single-thimble MCMC with exact Jacobian flow.
5. Approximate Jacobian or tangent-space proposals.
6. Tempered or worldvolume thimble sampling.
7. Multi-thimble sampling with intersection data.

The lecture examples are chosen to climb this ladder gradually.

## Choosing Flow Time

Small $T$:

- Easier ODE solve.
- Less deformation.
- Larger sign problem.

Large $T$:

- Better phase alignment.
- Sharper target distribution.
- More stiffness and Jacobian instability.

There is an optimal numerical window, not an infinite-flow ideal.

## Choosing Proposal Scale

For random-walk Metropolis,

$$
x'=x+\epsilon \xi,
\qquad \xi\sim\mathcal N(0,I).
$$

Acceptance probability:

$$
\alpha=\min\left(1,e^{-S_{\rm eff}(x')+S_{\rm eff}(x)}\right).
$$

The proposal scale must be tuned using acceptance and autocorrelation, not acceptance alone.

## Diagnostics Checklist

Report:

- Acceptance rate.
- Burn-in behavior.
- Trace plots.
- Autocorrelation time.
- Effective sample size.
- Average residual phase.
- Flow failure rate.
- Sensitivity to $T$ and proposal scale.

## Airy Integral

We now return to the Airy action introduced in the thimble section:

$$
S(z;a)=-i\left(\frac{z^3}{3}+az\right).
$$

It has exact benchmark

$$
\int_{-\infty}^{\infty}dz\,e^{-S(z;a)}=2\pi\operatorname{Ai}(a).
$$

## Airy: Runnable Critical Points

```{julia}
using SpecialFunctions

S(z, a) = -im * (z^3 / 3 + a * z)
dS(z, a) = -im * (z^2 + a)

critical_points(a) = (im * sqrt(a), -im * sqrt(a))

a = 1.0
zplus, zminus = critical_points(a)
(zplus = zplus,
 S_zplus = S(zplus, a),
 zminus = zminus,
 S_zminus = S(zminus, a),
 exact = 2 * pi * airyai(a))
```

This cell fixes the saddle data and the exact number used in the following
finite-flow checks.

## Airy Critical Points

The saddle equation is

$$
\partial_zS=-i(z^2+a)=0.
$$

Thus

$$
z_\pm=\pm i\sqrt a
\qquad (a>0).
$$

For the standard Airy contour and $a>0$, the contributing saddle used here
is $z_+=i\sqrt a$.

## Airy Flow Equations

Write $z=x+iy$.  Since

$$
\dot z=\overline{-i(z^2+a)},
$$

one obtains an explicit two-dimensional real ODE.

The code uses the complex form directly and stores $(x,y)$ for the ODE solver.

## Airy Code Skeleton

```{julia}
S(z, a) = -im * (z^3 / 3 + a * z)
dS(z, a) = -im * (z^2 + a)
d2S(z, a) = -2im * z

function airy_flow_rhs!(du, u, p, t)
    z = complex(u[1], u[2])
    dz = conj(dS(z, p.a))
    du[1] = real(dz)
    du[2] = imag(dz)
end
```

## Airy: Runnable Finite Flow Cell

```{julia}
using DifferentialEquations
using LinearAlgebra

dS(z, a) = -im * (z^2 + a)
d2S(z, a) = -2im * z

function flow_rhs!(du, u, p, t)
    z = complex(u[1], u[2])
    J = complex(u[3], u[4])
    dz = conj(dS(z, p.a))
    dJ = conj(d2S(z, p.a) * J)
    du[1] = real(dz)
    du[2] = imag(dz)
    du[3] = real(dJ)
    du[4] = imag(dJ)
end

function flow_point(x; a=1.0, T=0.35)
    zc = im * sqrt(a)
    rho = 1.0 + 0im
    T == 0 && return zc + rho*x, rho
    u0 = [real(zc + rho*x), imag(zc + rho*x), real(rho), imag(rho)]
    sol = solve(ODEProblem(flow_rhs!, u0, (0.0, T), (; a)), Tsit5();
                reltol=1e-8, abstol=1e-10, save_everystep=false)
    u = sol.u[end]
    return complex(u[1], u[2]), complex(u[3], u[4])
end

flow_point(0.2; a=1.0, T=0.35)
```

The same flow map now determines the geometry, the finite-flow quadrature,
the MCMC target, and the residual phase.

## Airy: Plot the Flowed Contour

```{julia}
xs_flow = range(-2.0, 2.0; length=80)
points_T0 = [im * sqrt(1.0) + x for x in xs_flow]
points_T = first.(flow_point.(xs_flow; a=1.0, T=0.35))

xs_grid = range(-2.4, 2.4; length=160)
ys_grid = range(-1.4, 1.8; length=140)
reS = [real(S(x + im*y, 1.0)) for y in ys_grid, x in xs_grid]
imS = [imag(S(x + im*y, 1.0)) for y in ys_grid, x in xs_grid]

p_re = contour(xs_grid, ys_grid, reS; levels=18, fill=false,
               xlabel="Re z", ylabel="Im z", title="Re S and flowed contour",
               color=:viridis, linewidth=1.2, aspect_ratio=:equal,
               colorbar=false, label=false)
plot!(p_re, real.(points_T0), imag.(points_T0); label="tangent T=0", lw=3, color=:black, ls=:dash)
plot!(p_re, real.(points_T), imag.(points_T); label="flowed T=0.35", lw=4, color=:red)
scatter!(p_re, [0, 0], [1, -1]; label="critical points", ms=5, color=:white, markerstrokecolor=:black)

p_im = contour(xs_grid, ys_grid, imS; levels=18, fill=false,
               xlabel="Re z", ylabel="Im z", title="Im S contours",
               color=:plasma, linewidth=1.2, aspect_ratio=:equal,
               colorbar=false, label=false)
plot!(p_im, real.(points_T), imag.(points_T); label="flowed T=0.35", lw=4, color=:red)
scatter!(p_im, [0, 0], [1, -1]; label="critical points", ms=5, color=:white, markerstrokecolor=:black)

plot(p_re, p_im; layout=(1, 2), size=(1250, 520))
```

The red curve is the finite-flow approximation $\Sigma_T$.  The left panel
shows how it follows the steepest-growth geometry of $S_R$ near the saddle;
the right panel checks that $S_I$ varies much less along the flowed contour
than along a generic complex curve.

## Airy: Runnable Finite-Flow Quadrature

```{julia}
using QuadGK

function airy_weight(x; a=1.0, T=0.35)
    z, J = flow_point(x; a, T)
    J * exp(-S(z, a))
end

function airy_finite_flow_integral(; a=1.0, T=0.35, xmax=1.2)
    quadgk(x -> airy_weight(x; a, T), -xmax, xmax; rtol=1e-6)[1]
end

Ts = [0.0, 0.1, 0.2, 0.35]
exact_airy = 2 * pi * airyai(1.0)
[(T = T, estimate = airy_finite_flow_integral(T=T), error = abs(airy_finite_flow_integral(T=T) - exact_airy))
 for T in Ts]
```

The comparison separates three errors: finite $T$, finite `xmax`, and
ODE/quadrature tolerance.

## Airy Saddle Data

For $a>0$,

$$
z_\pm=\pm i\sqrt a.
$$

At the contributing saddle $z_+=i\sqrt a$,

$$
S(z_+)=\frac{2}{3}a^{3/2}.
$$

The other saddle has the opposite real action.  This immediately explains why saddle relevance is contour-dependent.

## Airy: Flow-Time Scan

For each flow time,

$$
I_T=\int_{-x_{\max}}^{x_{\max}}dx\,J_T(x)e^{-S(z_T(x))}.
$$

Compare:

$$
I_T \quad \text{vs.}\quad 2\pi Ai(a).
$$

This demonstrates finite-flow bias, truncation in tangent coordinate, and numerical ODE error.

The same scan also gives a stability check: a useful $T$ window should not
change the estimate beyond the expected truncation and solver tolerances.

## From Flowed Integral to MCMC Target

After the finite-flow map, the one-thimble observable has the form

$$
\langle O\rangle
\approx
\frac{
\int dx\,O(z_T(x))\,e^{i\theta_{\rm res}(x)}e^{-S_{\rm eff}(x)}
}{
\int dx\,e^{i\theta_{\rm res}(x)}e^{-S_{\rm eff}(x)}
}.
$$

The positive target is

$$
\pi(x)=\frac{1}{Z_{\rm eff}}e^{-S_{\rm eff}(x)},
\qquad
Z_{\rm eff}=\int dx\,e^{-S_{\rm eff}(x)}.
$$

MCMC is useful precisely because the normalization $Z_{\rm eff}$ is not
needed to generate the chain.

## Airy: Metropolis on the Flowed Target

Now reuse the Airy finite-flow map.  For one tangent coordinate,

$$
S_{\rm eff}(x)=\operatorname{Re}S(z_T(x))-\log|J_T(x)|.
$$

Thus the log target is

$$
\log\pi(x)=-S_{\rm eff}(x)
=-\operatorname{Re}S(z_T(x))+\log|J_T(x)|.
$$

```{julia}
function logpi_airy_flowed(x; a=1.0, T=0.35)
    z, J = flow_point(x; a, T)
    -real(S(z, a)) + log(abs(J))
end

airy_chain, airy_acc = metropolis(x -> logpi_airy_flowed(x; a=1.0, T=0.35);
                                  nsteps=3_000, step=0.45, x0=0.0)

(acceptance = airy_acc,
 mean_tangent = mean(airy_chain),
 std_tangent = std(airy_chain))
```

This is the smallest live version of a thimble MCMC run: every proposed
$x$ is mapped to $z_T(x)$ before the accept/reject step.

## Airy: Residual-Phase Estimator

For an unnormalized integral, plain MCMC still needs the normalization
of $e^{-S_{\rm eff}}$.  For normalized observables, however, that
normalization cancels:

$$
\langle O\rangle
\approx
\frac{\langle O(z_T(x))e^{i\theta_{\rm res}(x)}\rangle_{\pi}}
{\langle e^{i\theta_{\rm res}(x)}\rangle_{\pi}},
$$

where

$$
e^{i\theta_{\rm res}(x)}
=\frac{J_T(x)}{|J_T(x)|}e^{-i\,\operatorname{Im}S(z_T(x))}.
$$

```{julia}
function residual_phase_airy(x; a=1.0, T=0.35)
    z, J = flow_point(x; a, T)
    J / abs(J) * exp(-im * imag(S(z, a)))
end

phases = residual_phase_airy.(airy_chain; a=1.0, T=0.35)
mean_phase = mean(phases)
(average_phase = mean_phase, abs_average_phase = abs(mean_phase))
```

This is the diagnostic that tells us whether the flowed contour has
actually reduced the sign problem.

## MCMC Is Not Just a Sampling Button

A Markov chain can look plausible and still be wrong.

Common failure modes:

- The chain stays in one mode.
- Burn-in is too short.
- Proposal scale hides long autocorrelation.
- Residual phase denominator is too small.
- The sampled thimble is not the full contributing cycle.

## Autocorrelation

For samples $x_k$, define

$$
\rho(t)=\frac{\langle (x_k-\bar x)(x_{k+t}-\bar x)\rangle}{\langle (x_k-\bar x)^2\rangle}.
$$

The integrated autocorrelation time is approximately

$$
\tau_{\rm int}=1+2\sum_{t=1}^{t_*}\rho(t).
$$

Effective sample size:

$$
N_{\rm eff}\approx \frac{N}{\tau_{\rm int}}.
$$

## Reading MCMC Diagnostics

It uses a simple multimodal target as a stand-in for a difficult flowed distribution.

The point is not physics realism; the point is to teach how to read diagnostics.

## MCMC Diagnostics: Runnable Core

```{julia}
using Statistics
using StatsBase

function autocorr_1d(x, maxlag)
    y = x .- mean(x)
    denom = sum(abs2, y)
    [sum(y[1:end-l] .* y[1+l:end]) / denom for l in 0:maxlag]
end

function tau_int(ac)
    cutoff = findfirst(<(0), ac)
    m = isnothing(cutoff) ? length(ac) : max(1, cutoff - 1)
    1 + 2 * sum(ac[2:m])
end

samples = randn(5_000)
ac = autocorr_1d(samples, 50)
(mean = mean(samples), std = std(samples), tau = tau_int(ac), ess = length(samples) / tau_int(ac))
```

Replace `samples` by the tangent-coordinate chain from the Airy run and
the same diagnostics apply.

## MCMC Diagnostics: Plot the Chain

![](assets/mcmc-diagnostics-pipeline.svg){fig-alt="MCMC diagnostics pipeline from chain to trace, histogram, autocorrelation, ESS, and decision." width=82%}

```{julia}
using Random

function rw_metropolis(logp; nsteps=4_000, step=0.65, x0=0.0)
    x = x0
    xs = Vector{Float64}(undef, nsteps)
    accepts = 0
    for n in 1:nsteps
        proposal = x + step * randn()
        if log(rand()) < logp(proposal) - logp(x)
            x = proposal
            accepts += 1
        end
        xs[n] = x
    end
    xs, accepts / nsteps
end

logp_mixture(x) = log(0.55 * exp(-(x + 2.0)^2 / 0.7) +
                      0.45 * exp(-(x - 2.2)^2 / 0.9))
chain, acc = rw_metropolis(logp_mixture)
(acceptance = acc, mean = mean(chain), std = std(chain))
```

## MCMC Diagnostics: Trace, Histogram, ACF

```{julia}
ac = autocorr_1d(chain, 80)
p1 = plot(chain; label=false, xlabel="step", ylabel="x", title="trace")
p2 = histogram(chain; bins=45, normalize=true, label=false,
               xlabel="x", ylabel="density", title="histogram")
p3 = plot(0:80, ac; marker=:circle, label=false,
          xlabel="lag", ylabel="acf", title="autocorrelation")
plot(p1, p2, p3; layout=(1, 3), size=(1100, 330))
```

Now the diagnostic vocabulary has pictures attached: a sticky trace,
an undersampled mode, and a slow ACF are different failure modes.

## Adaptive Metropolis

A practical adaptive rule adjusts the proposal covariance during warmup:

$$
\Sigma_{\rm prop}=s_d\,\widehat{\operatorname{Cov}}(x)+\epsilon I,
\qquad s_d\approx \frac{2.38^2}{d}.
$$

For thimble simulations, adaptation must be handled carefully because each proposal is expensive.

## Differential Evolution MCMC

For multiple chains $x_i$, propose

$$
x_i'=x_i+\gamma(x_j-x_k)+\epsilon.
$$

This is useful when the target has complicated covariance or multiple scales.

The original spinfoam propagator code used a Differential Evolution Adaptive Metropolis strategy.

## Worldvolume TLTM Motivation

The finite-flow method above fixes one flow time $T$ and samples one
flowed surface $\Sigma_T$.

That choice creates a real algorithmic tension:

- Small $T$: easier mixing but worse sign problem.
- Large $T$: better sign problem but harder mixing.
- Very large $T$: ODEs and Jacobians become expensive or unstable.

Worldvolume tempering uses flow time itself as a Monte Carlo coordinate.

The sampler can move through easy-mixing slices and still measure on
better sign-problem slices.

There is a second issue.  A single-thimble construction starts near one
saddle and then asks how much that saddle contributes.  TLTM instead flows
the original real integration manifold:

$$
x\in\mathbb R^n
\longmapsto
z_t(x).
$$

The relative weights of several contributing saddle regions are inherited
from the original integration cycle.  They are not assigned by hand from
separate single-saddle simulations.

## Worldvolume TLTM

Fukuma and Matsumoto propose sampling the worldvolume

$$
\mathcal R=\bigcup_{t\in[T_0,T_1]}\Sigma_t.
$$

The flow time $t$ becomes a dynamical tempering coordinate.

![](assets/worldvolume-tltm-schematic.svg){fig-alt="Schematic stack of flowed manifolds forming a worldvolume." width=88%}

Geometry schematic redrawn from the construction in Fukuma-Matsumoto,
[arXiv:2012.08468](https://arxiv.org/abs/2012.08468).

## Real-Manifold Flow and Saddle Weights

For each fixed $t$, the flowed surface is a deformation of the original
cycle:

$$
\Sigma_t=\{z_t(x)\mid x\in\mathbb R^n\},
\qquad
\Sigma_0=\mathbb R^n.
$$

If no singularity is crossed,

$$
\int_{\mathbb R^n} dx\,e^{-S(x)}
=
\int_{\Sigma_t} dz\,e^{-S(z)}.
$$

When the original cycle decomposes as

$$
\mathbb R^n\sim \sum_\sigma n_\sigma\mathcal J_\sigma,
$$

the flowed real manifold carries all sectors with $n_\sigma\neq0$ together.
This is different from choosing one saddle and sampling only its local
tangent space.

## Multi-Thimble Ergodicity

At large flow time, different contributing regions may become separated by
large barriers on one fixed $\Sigma_T$.

Then a Markov chain can have good local acceptance but still fail to move
between saddle regions.

Worldvolume sampling uses the extra coordinate $t$ to connect two regimes:

$$
t\approx 0:
\text{regions connected but phase problem large},
$$

$$
t\gg0:
\text{phase improved but regions separated}.
$$

This is why TLTM is also a multi-thimble ergodicity method, not only a
method for choosing flow time.

## Two-Saddle Model for the Next Questions

Consider a one-dimensional complex integral with two real wells:

$$
Z(g,\eta)=\int_{\mathbb R}dx\,e^{-S(x;g,\eta)},
\qquad
S(z;g,\eta)=\frac{(z^2-1)^2}{4g}+i\eta z .
$$

The saddle equation is

$$
\partial_zS=\frac{z(z^2-1)}{g}+i\eta=0 .
$$

For small $\eta$, two critical points are near $z=\pm1$.  A
single-thimble simulation begun near only one saddle does not determine
the relative weight of the other.  Flowing the original real axis deforms
both regions together:

$$
\dot z=\overline{\partial_zS(z)},\qquad
\Sigma_t=\{z_t(x)\mid x\in\mathbb R\}.
$$

## Two-Saddle Model: Critical Data

```{julia}
S_two(z; g=1.0, eta=0.8) = (z^2 - 1)^2 / (4g) + im * eta * z
dS_two(z; g=1.0, eta=0.8) = z * (z^2 - 1) / g + im * eta
d2S_two(z; g=1.0, eta=0.8) = (3z^2 - 1) / g

function newton_saddle(z0; g=1.0, eta=0.8, niter=30)
    z = complex(z0)
    for _ in 1:niter
        z -= dS_two(z; g, eta) / d2S_two(z; g, eta)
    end
    z
end

two_saddles = sort([newton_saddle(z0) for z0 in (-1.2, 0.0, 1.2)]; by=real)
labels = ["left", "middle", "right"]

[(label = label,
  Re_z = round(real(z), digits=4),
  Im_z = round(imag(z), digits=4),
  Re_S = round(real(S_two(z)), digits=4),
  Im_S = round(imag(S_two(z)), digits=4))
 for (label, z) in zip(labels, two_saddles)]
```

The two outer saddles have equal $\operatorname{Re}S$ and opposite
$\operatorname{Im}S$.  The middle saddle is different saddle data; whether
it contributes is a question about the integration cycle, not about
Newton's method.

## Two-Saddle Model: Flowed Real Axis

```{julia}
#| fig-width: 12
#| fig-height: 5.2
#| out-width: "100%"

using DifferentialEquations

function two_flow_rhs!(du, u, p, t)
    z = complex(u[1], u[2])
    dz = conj(dS_two(z; p.g, p.eta))
    du[1] = real(dz)
    du[2] = imag(dz)
end

function two_flow_point(x; g=1.0, eta=0.8, T=0.08)
    T == 0 && return complex(x, 0.0)
    prob = ODEProblem(two_flow_rhs!, [float(x), 0.0], (0.0, T), (; g, eta))
    sol = solve(prob, Tsit5(); reltol=1e-8, abstol=1e-10,
                save_everystep=false)
    u = sol.u[end]
    complex(u[1], u[2])
end

xs_two = range(-1.45, 1.45; length=120)
flowed_two = [two_flow_point(x; T=0.08) for x in xs_two]

xs_grid_two = range(-1.8, 1.8; length=160)
ys_grid_two = range(-0.75, 0.75; length=140)
reS_two = [real(S_two(x + im*y)) for y in ys_grid_two, x in xs_grid_two]

p_geom = contour(xs_grid_two, ys_grid_two, reS_two; levels=24,
                 xlabel="Re z", ylabel="Im z",
                 title="Re S contours and flowed real axis",
                 color=:viridis, colorbar=false, aspect_ratio=:equal,
                 linewidth=1.1, label=false)
plot!(p_geom, xs_two, fill(0.0, length(xs_two));
      label="original real axis", lw=3, color=:black, ls=:dash)
plot!(p_geom, real.(flowed_two), imag.(flowed_two);
      label="flowed Σ_T", lw=4, color=:red)
scatter!(p_geom, real.(two_saddles), imag.(two_saddles);
         label="critical points", ms=5, color=:white,
         markerstrokecolor=:black)

proxy_V(x, T) = real(S_two(two_flow_point(x; T)))
p_target = plot(; xlabel="x on original real axis",
                ylabel="Re S(z_t(x)) - minimum",
                title="fixed-slice target proxy")
for t in (0.0, 0.035, 0.08)
    vals = [proxy_V(x, t) for x in xs_two]
    plot!(p_target, xs_two, vals .- minimum(vals);
          lw=3, label="t=$(round(t, digits=3))")
end

plot(p_geom, p_target; layout=(1, 2), size=(1250, 520))
```

The right panel is only the geometric part of the target.  The full
positive worldvolume weight also contains $|\det J_t|$, $\alpha_t$, and
$W(t)$.  The plot isolates the important point: one fixed large-$t$ slice
can become effectively bimodal.

## Two-Saddle Model: What TLTM Changes

At a fixed flow time,

$$
p_T(x)\propto
|\det J_T(x)|\,e^{-\operatorname{Re}S(z_T(x))}
$$

may have two separated high-probability regions.  A local Markov chain can
then equilibrate inside one region while rarely tunneling to the other.

Worldvolume sampling replaces this by

$$
p(x,t)\propto
\alpha_t(x)|\det J_t(x)|
\exp[-\operatorname{Re}S(z_t(x))-W(t)] .
$$

A chain can move toward smaller $t$, cross through a more connected slice,
and return to larger $t$.  Since $\Sigma_t$ is the flow of the original
real cycle, the relative saddle regions are transported with the contour
rather than inserted as separate hand-tuned weights.

## Fixed Slice or Flow-Time Volume?

![](assets/fixed-vs-worldvolume.svg){fig-alt="Comparison between fixed-flow sampling on one slice and worldvolume sampling over flow time." width=88%}

The fixed-flow algorithm makes one hard choice before sampling:

$$
T=T_{\rm chosen}.
$$

Worldvolume tempering turns that choice into a coordinate:

$$
(x,t)\sim p(x,t),
\qquad t\in[T_0,T_1].
$$

This is the conceptual move that connects thimble geometry to tempering.

## Worldvolume: Paper Figure

![](assets/fukuma-matsumoto-fig1.png){fig-alt="Fig. 1 from Fukuma and Matsumoto showing the worldvolume R between flowed surfaces." width=76%}

Source: Fukuma and Matsumoto, Fig. 1, arXiv:2012.08468.

## Worldvolume Integral

Let $x\in\mathbb R^n$ parametrize the original real variables and
$z_t(x)$ be the antiholomorphic gradient-flow map.

On one flowed surface,

$$
Z_t=\int dx\;\det J_t(x)\,e^{-S(z_t(x))}.
$$

Worldvolume TLTM inserts a tempering weight $e^{-W(t)}$ and samples

$$
Z_{\mathcal R}^{\rm pq}=
\int_{T_0}^{T_1}dt\,e^{-W(t)}
\int dx\;\alpha_t(x)\,|\det J_t(x)|\,e^{-\operatorname{Re}S(z_t(x))}.
$$

The factor $\alpha_t(x)$ is the local normal factor converting
$dx\,dt$ into the induced volume element on $\mathcal R$.

## Worldvolume Geometry: ADM Data

![](assets/fukuma-matsumoto-fig2-3.png){fig-alt="Fukuma-Matsumoto Figs. 2 and 3 showing tangent basis and ADM geometry on the worldvolume." width=78%}

Coordinates on $\mathcal R$ are

$$
\xi^\mu=(t,x^a),\qquad z=z_t(x).
$$

The tangent basis is

$$
E_0=\frac{\partial z}{\partial t},\qquad
E_a=\frac{\partial z}{\partial x^a}.
$$

The induced metric is

$$
g_{\mu\nu}=E_\mu\cdot E_\nu.
$$

## Worldvolume Geometry: Lapse and Volume

In ADM form,

$$
ds^2=\alpha^2dt^2+
\gamma_{ab}\left(dx^a+\beta^a dt\right)
\left(dx^b+\beta^b dt\right).
$$

Here $\gamma_{ab}=E_a\cdot E_b$ is the metric on $\Sigma_t$,
$\beta^a$ is the shift, and $\alpha$ is the lapse.

The volume element is

$$
Dz=\sqrt{\det g}\,dt\,dx
=\alpha\,|\det J_t(x)|\,dt\,dx.
$$

This is why $\alpha$ appears in worldvolume reweighting.

## Worldvolume Reweighting

The sampled distribution is positive:

$$
p(x,t)\propto
e^{-V(x,t)},\qquad
V(x,t)=\operatorname{Re}S(z_t(x))+W(t)-\log|\det J_t(x)|-\log\alpha_t(x).
$$

The observable is recovered with a residual phase and tempering factor.

Schematic estimator:

$$
\langle O\rangle
=
\frac{\langle O(z_t)\,A(z_t)\,e^{i\theta(x,t)}\rangle_{\mathcal R}}
{\langle A(z_t)\,e^{i\theta(x,t)}\rangle_{\mathcal R}}.
$$

The exact notation varies by implementation, but the message is stable:
sampling and measurement live on different weighted objects.

For the paper's notation, the reweighting factor has the structure

$$
A(z)=\alpha^{-1}(z)
\exp\!\left[i\varphi(z)-i\,\operatorname{Im}S(z)\right],
\qquad
e^{i\varphi(z)}=\frac{\det J_t(x)}{|\det J_t(x)|}.
$$

So the determinant phase is not gone; it has been isolated into a
measurable residual factor.

## Worldvolume: HMC Constraint Picture

![](assets/fukuma-matsumoto-fig5.png){fig-alt="Fukuma-Matsumoto Fig. 5 showing how constrained molecular dynamics finds the next point on the worldvolume." width=78%}

One molecular-dynamics step proposes a point and then solves constraints
so that the result remains on $\mathcal R$.

Schematic RATTLE form:

$$
\pi_{1/2}=\pi-\frac{\Delta s}{2}\partial V(z)-\lambda^a F_a(z),
$$

$$
z'=z+\Delta s\,\pi_{1/2},\qquad z'\in\mathcal R,
$$

$$
\pi'=\pi_{1/2}-\frac{\Delta s}{2}\partial V(z')-\lambda'^a F_a(z'),
\qquad \pi'\in T_{z'}\mathcal R.
$$

The multipliers $\lambda^a,\lambda'^a$ enforce the constraints.

## What Worldvolume Sampling Buys

It addresses two high-dimensional problems:

- Multimodal or sharply curved thimble distributions.
- Poor tunneling between regions at a single fixed flow time.
- Repeated expensive determinant evaluation during configuration generation.

In practice one tunes $W(t)$ so that the histogram of $t$ is not trapped
near either endpoint.

The fixed-flow Airy example already contains the pieces that a
worldvolume implementation generalizes: flow map, Jacobian, effective
action, residual phase, and MCMC diagnostics.

## Fixed-Flow Versus Worldvolume Sampling

At fixed flow time, the simulation lives on one slice:

$$
x\in\mathbb R^n
\longmapsto
z_T(x)\in\Sigma_T .
$$

Worldvolume sampling instead uses

$$
(x,t)\in\mathbb R^n\times[T_0,T_1]
\longmapsto
z_t(x)\in\mathcal R .
$$

This adds one coordinate but removes a rigid algorithmic choice. The chain
can move at smaller $t$, where proposals are gentler, and still measure
configurations at larger $t$, where the phase is better controlled.

## Worldvolume Probability Target

The positive sampling weight has the schematic form

$$
p(x,t)\propto
\alpha_t(x)\,|\det J_t(x)|\,
\exp[-\operatorname{Re}S(z_t(x))-W(t)] .
$$

Equivalently,

$$
V(x,t)=
\operatorname{Re}S(z_t(x))
+W(t)
-\log|\det J_t(x)|
-\log\alpha_t(x).
$$

The role of $W(t)$ is not physics; it is tempering.  It is tuned so the
chain visits the useful range of flow times instead of sticking at one end.

## What Is Actually Measured?

Sampling uses $p(x,t)$, but the original integral still contains the complex
orientation and action phase.

The residual complex factor is schematically

$$
e^{i\theta_{\mathcal R}(x,t)}
=
\frac{\det J_t(x)}{|\det J_t(x)|}
\exp[-i\,\operatorname{Im}S(z_t(x))].
$$

Because the worldvolume measure contains the normal factor $\alpha_t$, the
measurement also carries an $\alpha_t^{-1}$ factor.  The important lesson is:

$$
\text{worldvolume sampling improves exploration, but it does not erase
the residual phase.}
$$

## Worldvolume Algorithm Skeleton

For an implementation, the clean separation is:

1. Define the flow map $z_t(x)$ and tangent Jacobian $J_t(x)$.
2. Add a tempering potential $W(t)$.
3. Sample $(x,t)$ on $\mathcal R$ with HMC or Metropolis proposals.
4. Include the normal factor and residual phase in measurement.
5. Check the $t$ histogram, average phase, ESS, and flow-time stability.

For large spinfoam variables, steps 1 and 4 are the expensive ones.

This is why a worldvolume implementation is usually prepared as a research
calculation rather than improvised during the first pass through the method.

## Worldvolume: Runnable Tempering Toy

```{julia}
#| fig-width: 8.8
#| fig-height: 4.8
#| out-width: "86%"

flow_times = range(0.0, 4.0; length=100)
mixing_proxy = exp.(-0.9 .* flow_times)
phase_proxy = 1 .- exp.(-1.25 .* flow_times)
overlap_proxy = mixing_proxy .* phase_proxy

plot(flow_times, mixing_proxy; label="mixing proxy", lw=3,
     xlabel="flow time T", ylabel="relative scale",
     title="Why sample a flow-time worldvolume?")
plot!(flow_times, phase_proxy; label="phase proxy", lw=3)
plot!(flow_times, overlap_proxy; label="useful overlap", lw=3)
```

This minimal proxy separates the two pressures that TLTM balances:
mixing prefers small flow time, while the residual phase usually improves
at larger flow time.

## From Worldvolume to Intersection Numbers

Worldvolume sampling reduces the need to simulate each thimble separately:
it samples a deformation of the original cycle.

It does not remove the conceptual question

$$
\mathbb R^n\sim\sum_\sigma n_\sigma\mathcal J_\sigma .
$$

For interpretation, asymptotics, and parameter continuation, we still need
to know which upward cycles intersect the original contour:

$$
n_\sigma=\langle \mathcal C,\mathcal K_\sigma\rangle .
$$

The two-saddle model makes the issue concrete: both outer saddles are
visible locally, but the global cycle decides which saddle sectors belong
to the original integral and how their contributions jump across Stokes
walls.

## Which Saddles Contribute?

A saddle exists if

$$
\partial S(z_\sigma)=0.
$$

But it contributes only if

$$
n_\sigma=\langle\mathcal C,\mathcal K_\sigma\rangle\neq0.
$$

This is one of the central practical difficulties for complex critical points.

![](assets/intersection-thimbles-schematic.svg){fig-alt="Schematic original contour crossing upward cycles, selecting contributing thimbles." width=88%}

## Why Critical Points Are Not Enough

Solving

$$
\partial S(z_\sigma)=0
$$

only gives candidates.  The original integration cycle decides which
candidates are present in the decomposition:

$$
\mathcal C \sim
\sum_{\sigma} n_\sigma\mathcal J_\sigma .
$$

Two saddles can have similar local actions and Hessians, but one may have
$n_\sigma=0$.  Including it would change the problem.

The local task is saddle finding.  The global task is cycle decomposition.

## Intersection Numbers in One Sentence

Picard-Lefschetz decomposition is not just a list of saddles:

$$
\mathcal C=\sum_\sigma n_\sigma\mathcal J_\sigma.
$$

The integer $n_\sigma$ is an oriented intersection number between the
original cycle $\mathcal C$ and the upward cycle $\mathcal K_\sigma$.

Equivalently,

$$
n_\sigma=\#\left(\mathcal C\cap\mathcal K_\sigma\right)_{\rm oriented}.
$$

Therefore:

- Newton may find many complex critical points.
- Only saddles with $n_\sigma\neq0$ enter the original integral.
- Stokes phenomena occur when these integers jump as parameters vary.

For spinfoams this matters because complex critical points can exist even
when real Regge-like critical points are absent or suppressed.

## Stokes Phenomena

![](assets/stokes-jump-schematic.svg){fig-alt="Schematic Stokes jump in which an additional upward cycle intersects the original contour." width=88%}

Intersection numbers can change when parameters cross a Stokes wall.  A
typical condition is that two saddles align in phase:

$$
\operatorname{Im}S(z_\sigma)
=
\operatorname{Im}S(z_\tau).
$$

Then a flow line can connect the two critical points, and the thimble
decomposition jumps:

$$
n_\sigma \longrightarrow n_\sigma \pm n_\tau .
$$

This is why parameter continuation must monitor the global flow geometry,
not only the saddle equations.

## Upward Flow Instability

To compute $n_\sigma$, one studies upward cycles.

Direct shooting is unstable:

- Small initial errors grow exponentially.
- Long flow times amplify numerical noise.
- Boundary conditions are global.

This is why boundary-value methods are useful.

![](assets/multiple-shooting-schematic.svg){fig-alt="A long unstable flow split into short multiple-shooting segments." width=88%}

If the upward flow is written as

$$
\dot z=F(z)=\overline{\partial S(z)},
$$

then direct shooting tries to tune one initial tangent vector $v$ so that

$$
z(T;z_\sigma+\epsilon v)\in\mathcal C.
$$

The sensitivity is governed by the variational equation

$$
\dot{\delta z}=DF(z(t))\,\delta z,
$$

which is the source of exponential error growth.

## Multiple Shooting for Intersection Numbers

Shoji and Trailovic propose a multiple-shooting strategy:

1. Split a long flow into short segments.
2. Treat segment endpoints as unknowns.
3. Enforce flow matching between segments.
4. Enforce saddle and original-cycle boundary conditions.
5. Solve the nonlinear system.

Reference: [arXiv:2510.06334](https://arxiv.org/abs/2510.06334).

## Multiple Shooting: Residual Equations

Let $\Phi_h$ be the short-time flow map over one segment.

Unknown segment endpoints:

$$
Z=(z_0,z_1,\ldots,z_K).
$$

Matching residuals:

$$
R_k(Z)=z_{k+1}-\Phi_h(z_k)=0,
\qquad k=0,\ldots,K-1.
$$

Boundary conditions:

$$
z_0=z_\sigma+\epsilon v,\qquad
z_K\in\mathcal C.
$$

Solve all residuals together by Newton or trust-region methods.

## Multiple Shooting: Unknowns and Constraints

A useful way to explain the implementation is:

| Object | Numerical role |
|---|---|
| Segment endpoints | variables in a nonlinear solve |
| Short-time flow maps | local matching constraints |
| Saddle boundary data | starts on $\mathcal K_\sigma$ |
| Original-cycle condition | endpoint lies on $\mathcal C$ |
| Orientation/Jacobian sign | determines the sign of $n_\sigma$ |

This turns one unstable long initial-value problem into a structured
boundary-value problem.

The price is a larger nonlinear system, but it is better conditioned and
more parallelizable.

## Multiple Shooting: Why Mention It?

It answers a question that the basic lecture otherwise leaves unresolved:

> Which complex saddles should be included?

This is essential for real-time path integrals, Stokes phenomena, and complex spinfoam critical points.

## Executable Multiple-Shooting Toy

The toy model uses an unstable scalar flow

$$
\dot y=\lambda y.
$$

Single shooting is exponentially sensitive to initial errors, while multiple shooting imposes short-time matching constraints.  This is the numerical idea behind more serious intersection-number computations.

For the toy flow, a segment of length $h=T/K$ gives

$$
y_{k+1}=e^{\lambda h}y_k,
$$

so the multiple-shooting equations are

$$
y_{k+1}-e^{\lambda h}y_k=0,\qquad y_K=1.
$$

## Multiple Shooting: Runnable Core

```{julia}
using LinearAlgebra

function multiple_shooting_matrix(lambda, T, K)
    h = T / K
    short = exp(lambda * h)
    A = zeros(K + 1, K + 1)
    b = zeros(K + 1)
    for k in 1:K
        A[k, k] = -short
        A[k, k+1] = 1.0
    end
    A[K+1, K+1] = 1.0
    b[K+1] = 1.0
    A, b
end

lambda, T = 2.0, 8.0
[(K = K, cond = cond(multiple_shooting_matrix(lambda, T, K)[1]))
 for K in (1, 2, 4, 8, 16)]
```

The output shows the numerical moral: replacing one long unstable map by
short matching constraints changes the conditioning problem.

## Multiple Shooting: Conditioning Plot

```{julia}
Ks = [1, 2, 4, 8, 16, 32]
conds = [cond(multiple_shooting_matrix(lambda, T, K)[1]) for K in Ks]

plot(Ks, conds; marker=:circle, lw=3, xscale=:log2, yscale=:log10,
     xlabel="number of segments K", ylabel="cond(A)",
     title="Multiple-shooting linear system conditioning", label=false)
```

For nonlinear thimble flows the matrix is larger, but the numerical message
is the same: local matching constraints are less fragile than one long
shooting map.

## Resummation: Why It Belongs Here

Saddle expansions are often divergent:

$$
F(g)\sim\sum_{n=0}^\infty a_n g^n,
\qquad a_n\sim n!.
$$

This does not mean the expansion is useless.

It means the expansion must be interpreted asymptotically or resummed.

## Why Perturbation Series Diverge

Suppose an expansion around one saddle has large-order behavior

$$
a_n\sim C\,\frac{\Gamma(n+\beta)}{A^{n+\beta}}.
$$

Then the $n$th term behaves roughly as

$$
a_n g^n
\sim
C\,\Gamma(n+\beta)\left(\frac{g}{A}\right)^n .
$$

For any fixed nonzero $g$, the factorial eventually wins over the power.
So the radius of convergence is usually zero.

The best truncation occurs near the smallest term:

$$
n_*\sim \frac{|A|}{|g|},
\qquad
\text{remainder}\sim e^{-A/g}.
$$

That exponential scale is already pointing to another saddle.

## Divergence as Missing Saddle Data

For a multi-saddle integral with small parameter $g$,

$$
Z(g)=\int_{\mathcal C} dz\,e^{-S(z)/g},
$$

a transseries has the schematic form

$$
Z(g)\sim
\sum_\sigma n_\sigma
e^{-S_\sigma/g}
\Phi_\sigma(g),
\qquad
\Phi_\sigma(g)=\sum_{n\ge0}a_{\sigma,n}g^n .
$$

Large-order growth of the expansion around saddle $\sigma$ is controlled by
action differences

$$
A_{\tau\sigma}=S_\tau-S_\sigma .
$$

This is the resummation version of Picard-Lefschetz theory: perturbation
theory around one saddle carries information about the other saddles.

## Same Two-Saddle Data in the Borel Plane

Return to the two-saddle action

$$
S(z;g,\eta)=\frac{(z^2-1)^2}{4g}+i\eta z .
$$

Let $\sigma$ be the left saddle and $\tau$ the right saddle.  The action
difference

$$
A_{\tau\sigma}=S_\tau-S_\sigma
$$

has already appeared three times:

- in thimble geometry, it controls relative exponential weights;
- in Stokes phenomena, its phase controls when saddle sectors mix;
- in resurgence, it locates Borel singularities near $\xi=A_{\tau\sigma}$.

Thus the same complex-plane data helps decide how to sample, which
intersection numbers matter, and how to resum a divergent expansion.

## Two-Saddle Borel Data: Runnable Check

```{julia}
#| fig-width: 6.8
#| fig-height: 5.4
#| out-width: "68%"

outer = (two_saddles[1], two_saddles[3])
A_lr = S_two(outer[2]) - S_two(outer[1])
A_rl = -A_lr
stokes_angle = angle(A_lr)
display((A_right_left = A_lr,
         stokes_angle_radians = stokes_angle,
         stokes_angle_over_pi = stokes_angle / pi))

pA = plot(; xlabel="Re ξ", ylabel="Im ξ",
          title="Borel singularity directions from saddle actions",
          aspect_ratio=:equal, legend=:outerright)
scatter!(pA, [real(A_lr)], [imag(A_lr)];
         label="A_right,left", ms=8, color=:red)
scatter!(pA, [real(A_rl)], [imag(A_rl)];
         label="A_left,right", ms=8, color=:blue)
plot!(pA, [0, real(A_lr)], [0, imag(A_lr)];
      lw=3, color=:red, label="Stokes ray")
plot!(pA, [0, real(A_rl)], [0, imag(A_rl)];
      lw=3, color=:blue, label=false)
hline!(pA, [0.0]; color=:black, lw=1, ls=:dash, label=false)
vline!(pA, [0.0]; color=:black, lw=1, ls=:dash, label=false)
pA
```

For real positive coupling, the Laplace ray is the positive real
$\xi$-axis.  In this parameter choice the outer-saddle singularity is
off that ray.  Rotating a parameter or the coupling can place it on the
ray, where lateral Borel sums and thimble coefficients jump together.

## Borel Transform

![](assets/borel-plane-saddles.svg){fig-alt="Borel plane singularities representing other saddles and possible Stokes phenomena." width=82%}

Define

$$
\mathcal B F(\xi)=\sum_{n=0}^\infty \frac{a_n}{n!}\xi^n.
$$

Then formally

$$
\mathcal S F(g)=\int_0^\infty d\xi\,e^{-\xi/g}\mathcal B F(\xi).
$$

Singularities in the Borel plane encode other saddles.

If the nearest singularity sits at $\xi=A$, then large order often has
the schematic form

$$
a_n\sim \frac{\Gamma(n+\beta)}{A^{n+\beta}}
\left(c_0+\frac{c_1}{n}+\cdots\right).
$$

That is the resurgence bridge to thimble data: the perturbative series
around one saddle knows about actions of other saddles.

## Why Borel Singularities Matter

![](assets/borel-stokes-line.svg){fig-alt="Borel plane pole on the Laplace ray with lateral contours above and below it." width=86%}

If

$$
a_n\sim \frac{n!}{A^{n+1}},
$$

then the Borel transform behaves like

$$
\mathcal B F(\xi)
\sim
\sum_{n=0}^{\infty}\frac{\xi^n}{A^{n+1}}
=
\frac{1}{A-\xi}.
$$

The factorial divergence has become a pole in the Borel plane.

For $g>0$, the Laplace integral

$$
\int_0^\infty d\xi\,e^{-\xi/g}\frac{1}{A-\xi}
$$

is ambiguous if $A$ lies on the positive real axis. Going above or below
the pole gives different answers:

$$
\mathcal S_{0^+}F(g)-\mathcal S_{0^-}F(g)
=
2\pi i\,e^{-A/g}
\,\operatorname{Res}_{\xi=A}\mathcal BF(\xi).
$$

The ambiguity has the same exponential size as the missing nonperturbative
saddle sector.

## Stokes Lines in the Borel Plane

For complex coupling

$$
g=|g|e^{i\theta},
$$

the Borel-Laplace sum is taken along the ray

$$
\xi=e^{i\theta}s,\qquad s\in[0,\infty).
$$

A Stokes line occurs when a Borel singularity sits on this integration ray:

$$
\arg g=\arg A_{\tau\sigma}.
$$

Equivalently, two saddle exponentials have aligned phases:

$$
\operatorname{Im}\frac{S_\tau-S_\sigma}{g}=0.
$$

Across such a line, the lateral sums jump and the transseries coefficient
of another saddle must jump with them. This is the resurgent counterpart
of a Picard-Lefschetz Stokes jump.

## What Complex-Plane Structure Buys Us

Finding singularities in the complex Borel plane does not make the original
series convergent. It does something more useful:

- It identifies the exponential scales $e^{-A/g}$ missing from perturbation theory.
- It tells us where lateral resummation is necessary.
- It predicts which saddle sectors can cancel Borel ambiguities.
- It gives numerical diagnostics for Pade poles and integration contours.

This is why resummation is not separate from thimbles. Both are ways of
organizing the same complex saddle structure.

## Quantum Mechanics and Field Theory

In the double-well quantum-mechanics path integral, perturbation theory
around one minimum is not the whole answer. Instantons contribute

$$
\Delta E\sim e^{-S_{\rm inst}/\hbar},
$$

and instanton-anti-instanton sectors cancel ambiguities of perturbation
theory around the vacuum.

In quantum field theory the same logic appears with more structure:

- instantons and complex saddles give semiclassical Borel singularities
- renormalons encode factorial growth from momentum regions
- Stokes phenomena control how nonperturbative sectors enter observables

The practical lesson for numerics is modest but important: inspect the
complex structure before trusting a truncated perturbative expansion.

## Borel-Pade in Practice

Given finitely many coefficients:

1. Build truncated Borel series.
2. Replace it by a Pade approximant.
3. Perform the Laplace integral numerically.
4. Compare to direct quadrature if available.

The next example keeps the arithmetic small enough to run during the
lecture, while still showing why divergent series can be useful.

## 0D Phi-Four Integral

Integral:

$$
Z(g)=\int_{-\infty}^{\infty}dx\,
\exp\left[-\frac{x^2}{2}-gx^4\right].
$$

It compares exact quadrature, optimal truncation, and Borel-Pade.

## 0D Phi-Four: Runnable Core

```{julia}
using QuadGK

Z_exact(g) = quadgk(x -> exp(-x^2 / 2 - g*x^4), -Inf, Inf)[1]

function double_factorial_odd(m)
    m <= 0 && return big(1)
    prod(big(k) for k in 1:2:m)
end

coeff(n) = Float64((-1)^n * double_factorial_odd(4 * n - 1) / factorial(big(n)))

function perturbative_Z(g, N)
    sqrt(2 * pi) * sum(coeff(n) * g^n for n in 0:N)
end

g = 0.03
[(N = N, estimate = perturbative_Z(g, N)) for N in 0:8],
Z_exact(g)
```

The table is deliberately simple: as $N$ grows, the estimates first improve
and then eventually drift because the perturbative series is asymptotic.

## 0D Phi-Four: Coefficients

Expanding around the Gaussian saddle gives

$$
Z(g)=\sqrt{2\pi}\sum_{n=0}^\infty
(-g)^n\frac{(4n-1)!!}{n!}.
$$

The coefficients grow factorially because

$$
(4n-1)!! \sim \frac{(4n)!}{2^{2n}(2n)!}.
$$

This is the simplest place to see divergent perturbation theory numerically.

## 0D Phi-Four: Three Estimates

For each $g$, compare:

1. Exact quadrature.
2. Truncated perturbation theory.
3. Borel-Pade resummation.

The best truncation order decreases as $g$ grows, while Borel-Pade remains
useful in a wider window.

## 0D Phi-Four: Minimal Borel-Pade Cell

For

$$
F(g)=\sum_{n=0}^{\infty}a_ng^n,
$$

use

$$
\mathcal B_NF(\xi)=\sum_{n=0}^N\frac{a_n}{n!}\xi^n,
\qquad
\mathcal S_{\rm BP}F(g)=
\int_0^\infty du\,e^{-u}\,[L/M]_{\mathcal B}(gu).
$$

```{julia}
using LinearAlgebra

function pade_from_series(a, L, M)
    @assert length(a) >= L + M + 1
    A = [a[L + row - col + 1] for row in 1:M, col in 1:M]
    b = [-a[L + row + 1] for row in 1:M]
    qtail = A \ b
    q = vcat(1.0, qtail)
    p = [sum(q[j+1] * a[k-j+1] for j in 0:min(k, M)) for k in 0:L]
    p, q
end

polyval(c, x) = sum(c[k+1] * x^k for k in 0:(length(c)-1))
padeval(p, q, x) = polyval(p, x) / polyval(q, x)

function borel_pade_Z(g; L=4, M=4)
    c = [coeff(n) / factorial(n) for n in 0:(L+M)]
    p, q = pade_from_series(c, L, M)
    sqrt(2*pi) * quadgk(u -> exp(-u) * padeval(p, q, g*u), 0, Inf; rtol=1e-7)[1]
end

g = 0.03
(exact = Z_exact(g),
 truncation_N8 = perturbative_Z(g, 8),
 borel_pade = borel_pade_Z(g))
```

This is a compact classroom version, not a production resummation routine.
In real work one scans Pade orders and checks for poles near the integration
ray.

## Resummation Caveat

Borel-Pade is not magic:

- Poles may appear near the integration contour.
- Lateral resummation may be needed.
- Ambiguities signal missing nonperturbative sectors.
- Numerical agreement must be checked against an independent benchmark when possible.

## Resurgence Link

Perturbation theory around one saddle carries information about other saddles through large-order growth.

Thimble decompositions and resurgence are two views of the same analytic structure:

- Thimbles organize integration cycles.
- Resurgence organizes asymptotic expansions.
- Stokes phenomena change how sectors combine.

## Quantum Mechanics Example

A useful next step beyond 0D is Euclidean quantum mechanics.

For a double-well potential,

$$
S_E[x]=\int d\tau\left[\frac12\dot x^2+\lambda(x^2-a^2)^2\right].
$$

Instantons are nonperturbative saddles connecting the two minima.

## Double-Well Instanton Numerics

It demonstrates:

- Discretized Euclidean action.
- Relaxation/gradient descent to an instanton-like path.
- Action comparison with the analytic continuum estimate.
- Relation to nonperturbative saddle sectors.

## Double Well: Runnable Core

```{julia}
function double_well_action(x; beta=8.0, a=1.0, lambda=1.0)
    N = length(x)
    h = beta / (N - 1)
    kinetic = sum((x[i+1] - x[i])^2 / (2h) for i in 1:(N-1))
    potential = h * sum(lambda * (x[i]^2 - a^2)^2 for i in 2:(N-1))
    kinetic + potential
end

N = 101
beta = 8.0
tau = range(-beta/2, beta/2; length=N)
x_trial = tanh.(tau)
(action = double_well_action(collect(x_trial); beta), endpoints = (x_trial[1], x_trial[end]))
```

The trial path is already close to the continuum instanton shape.  A short
relaxation step would improve the discretized action while keeping the
boundary values fixed.

## Double Well: One Relaxation Step

The finite-dimensional gradient is taken only over interior points:

$$
x_1=-a,\qquad x_N=a,\qquad
x_i\leftarrow x_i-\eta\,\frac{\partial S_E}{\partial x_i}.
$$

For the discretized action,

$$
\frac{\partial S_E}{\partial x_i}
=
\frac{2x_i-x_{i-1}-x_{i+1}}{h}
+4h\lambda x_i(x_i^2-a^2),
\qquad 2\le i\le N-1.
$$

```{julia}
function relax_double_well(x; beta=8.0, a=1.0, lambda=1.0, eta=0.02, sweeps=200)
    y = copy(x)
    N = length(y)
    h = beta / (N - 1)
    for _ in 1:sweeps
        grad = zeros(N)
        for i in 2:(N-1)
            grad[i] = (2*y[i] - y[i-1] - y[i+1]) / h +
                      4*h*lambda*y[i]*(y[i]^2 - a^2)
        end
        y[2:end-1] .-= eta .* grad[2:end-1]
        y[1], y[end] = -a, a
    end
    y
end

x0 = collect(range(-1.0, 1.0; length=N))
x_relaxed = relax_double_well(x0; beta)
(initial = double_well_action(x0; beta),
 relaxed = double_well_action(x_relaxed; beta))
```

## Double Well: Discretized Action

With lattice spacing $h$, use

$$
S_E[x]\approx\sum_{i=1}^{N-1}\frac{(x_{i+1}-x_i)^2}{2h}
+ h\sum_{i=1}^N \lambda(x_i^2-a^2)^2.
$$

Fix boundary conditions:

$$
x(-\beta/2)=-a,
\qquad
x(\beta/2)=a.
$$

Relaxing the interior points finds an instanton-like path.

## Double Well: Why This Supports the Lecture

This example connects:

- Saddles beyond the perturbative vacuum.
- Nonperturbative exponential weights.
- Large-order/resurgence intuition.
- Path-integral discretization as a finite-dimensional numerical problem.

It is the bridge from 0D examples to field-theory-style integrals.

## Field Theory Toy Example

For a lattice scalar field with complex source,

$$
S[\phi]=\sum_i\left[\frac12 m^2\phi_i^2+\lambda\phi_i^4+i h\phi_i\right]
+\frac{\kappa}{2}\sum_i(\phi_{i+1}-\phi_i)^2.
$$

Direct sampling of $e^{-S}$ has a phase problem.

A thimble method would complexify all $\phi_i$.

## Lattice Scalar Reweighting

It demonstrates:

- Real-field Metropolis sampling of the phase-quenched theory.
- Average phase decay with source strength and lattice size.
- A concrete sign-problem diagnostic before using thimbles.

## Lattice Scalar: Runnable Core

```{julia}
using Statistics

function average_phase_gaussian(N, h, m2; nsamples=20_000)
    sigma = 1 / sqrt(m2)
    phases = [exp(-im * h * sum(sigma .* randn(N))) for _ in 1:nsamples]
    mean(phases)
end

[(N = N, phase = average_phase_gaussian(N, 0.4, 1.0; nsamples=2_000))
 for N in (4, 8, 16, 32)]
```

Even this Gaussian toy version shows the volume trend that makes
phase-quenched reweighting fragile.

## Lattice Scalar: Phase-Quenched Estimator

Write

$$
S[\phi]=S_0[\phi]+ih\sum_i\phi_i.
$$

Sampling $e^{-S_0}$ gives

$$
\langle O\rangle
=\frac{\langle Oe^{-ih\sum_i\phi_i}\rangle_0}
{\langle e^{-ih\sum_i\phi_i}\rangle_0}.
$$

For independent Gaussian fields this denominator is analytic:

$$
\left\langle e^{-ih\sum_i\phi_i}\right\rangle_0
=
\exp\left[-\frac{Nh^2}{2m^2}\right].
$$

That formula is the simplest quantitative version of the exponential
average-phase problem.

## Lattice Scalar: Why It Is Not Yet a Thimble Demo

This example intentionally stops before contour deformation.

Its role is to show the problem that thimbles are supposed to solve:

- The phase-quenched chain may sample well.
- The reweighted observable may still be unusable.
- The average phase is the key diagnostic.

A full lattice thimble implementation would complexify all field variables and solve high-dimensional flow equations.

## Summation Problems in Spinfoams

Not every problem is a continuous integral.

Spinfoam numerics also includes sums like

$$
A(\Lambda)=\sum_{j\le\Lambda} a_j.
$$

The cutoff sequence may converge slowly or oscillate.

We need extrapolation and sequence acceleration.

## Summation Boost Methods

Useful tools:

- Richardson extrapolation for power-law cutoff errors.
- Shanks transformation for nearly geometric tails.
- Wynn epsilon algorithm for nonlinear acceleration.
- Pade approximants for rational reconstruction.

These methods need validation against known limits or stability windows.

## Summation Boost Example

It includes:

- Alternating Leibniz series.
- Power-law cutoff sequence.
- Shanks and Richardson examples.
- Guidance for spin-sum extrapolation diagnostics.

## Summation Boost: Runnable Core

```{julia}
function shanks(seq)
    [seq[k+1] - (seq[k+1] - seq[k])^2 /
     (seq[k+1] - 2 * seq[k] + seq[k-1]) for k in 2:(length(seq)-1)]
end

target = pi / 4
seq = [sum((-1.0)^k / (2 * k + 1) for k in 0:N) for N in 2:30]
boosted = shanks(seq)

(raw_error = abs(seq[end] - target),
 shanks_error = abs(boosted[end] - target))
```

For spin sums, the same idea is useful only after checking that the
cutoff error model is plausible.

## Spinfoam Data Bridge

![](assets/spinfoam-pipeline.svg){fig-alt="Pipeline from spinfoam boundary data through action, saddles, cycle data, sampling, diagnostics, and prediction." width=88%}

The bridge is not to run the full propagator calculation live.  The bridge
is to identify where each numerical object in the thimble lecture appears
in the spinfoam workflow.

Useful local material includes boundary data, action/Hessian construction,
and booster or vertex-amplitude data.  Those files are treated as source
material for prepared examples, not as the main classroom interface.

## Spinfoam: Data Inspection Cell

```{julia}
using JSON3

theta_phi = JSON3.read(read("../theta_phi_one_data.json", String))
areas = JSON3.read(read("../data_school/areadata.json", String))
kappa = JSON3.read(read("../data_school/kappadata.json", String))

(tetrahedra = length(theta_phi),
 angle_table_shape = (length(theta_phi), length(theta_phi[1]), length(theta_phi[1][1])),
 area_blocks = length(areas),
 first_angle_pair = theta_phi[1][1],
 kappa_blocks = length(kappa))
```

This cell is deliberately modest: it verifies that the boundary-data
objects have the expected nested structure before any expensive symbolic or
Monte Carlo step is attempted.

## Spinfoam: Minimal Computational Pipeline

For a real project, organize the computation as:

1. Generate or load boundary data.
2. Build action variables and gauge fixing.
3. Lambdify action, gradient, and Hessian.
4. Solve real critical equations as a baseline.
5. Search complex critical points by continuation or Newton methods.
6. Classify Hessian/Takagi data.
7. Decide contributing saddles.
8. Only then run finite-flow MCMC.

## Spinfoam: What Needs Precomputation

The following should usually be precomputed before a lecture:

- Boundary data and closure checks.
- Critical point candidates.
- Hessian spectra.
- Representative flow-time tests.
- Short pilot MCMC chains.

Live coding should focus on interpreting these objects, not waiting for heavy symbolic algebra.

## What We Do Not Run Live

A full spinfoam propagator thimble-MCMC run is not appropriate for a half-day first workshop.

Reasons:

- Dimension is high.
- Each MCMC proposal is expensive.
- Gauge fixing and data preparation are nontrivial.
- Diagnostics require long chains.

Instead, we build the complete method on smaller models and show exactly where spinfoam enters.

## Practical Workflow for a New Model

1. Write the holomorphic action and measure.
2. Identify singularities and allowed contour deformations.
3. Solve critical equations.
4. Classify saddles and Hessian directions.
5. Determine or estimate contributing thimbles.
6. Choose finite-flow or worldvolume sampling.
7. Validate against exact limits, perturbation theory, or known asymptotics.
8. Only then scale to high dimension.

## Validation Hierarchy

Use increasingly hard checks:

- Exact analytic answer if available.
- Direct quadrature in low dimension.
- Perturbative expansion and optimal truncation.
- Borel-Pade or other resummation.
- Known saddle asymptotics.
- Independent MCMC chains.
- Flow-time and cutoff stability.

## Failure Modes and What They Mean

| Symptom | Likely Cause |
|---|---|
| Low acceptance | proposal too large or target too sharp |
| High acceptance but bad ESS | proposal too small |
| Average phase near zero | residual sign problem |
| Strong flow-time drift | finite-flow bias or missing thimbles |
| ODE failures | too large flow time or singular region |
| Wrong exact benchmark | wrong saddle or contour |

## Method Pipeline

The examples build one pipeline:

$$
\text{positive measure}
\to
\text{MCMC}
\to
\text{sign problem}
\to
\text{complex saddles}
\to
\text{flowed geometry}
\to
\text{residual phase}
\to
\text{diagnostics}.
$$

Resummation and resurgence provide an independent view of the same complex
saddle structure through large-order growth and Borel singularities.

## Reading Map

For this workshop:

- Original spinfoam numerics slides in this repository.
- Witten, analytic continuation and Picard-Lefschetz theory.
- Fukuma and Matsumoto, Worldvolume TLTM, arXiv:2012.08468.
- Shoji and Trailovic, multiple shooting for intersection numbers, arXiv:2510.06334.
- Standard references on Borel resummation and resurgence.

## Takeaway

The numerical strategy is not a single algorithm.

It is a pipeline:

$$
\text{critical points}
\to
\text{cycle information}
\to
\text{flowed geometry}
\to
\text{MCMC}
\to
\text{residual phase and diagnostics}
\to
\text{asymptotic/resummed checks}.
$$

Spinfoam numerics needs all parts of this pipeline.

