Lefschetz Thimbles, MCMC, Resummation, and Spinfoam Numerics

A half-day numerical workshop

Hongguang Liu

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

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

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')

(stationary = [0.31578947368420995, 0.4736842105263162, 0.21052631578947384], residual = adjoint([1.1102230246251565e-16, -1.1102230246251565e-16, 0.0]))

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

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))

Precompiling packages...
    399.7 ms  ✓ Statistics
  1 dependency successfully precompiled in 0 seconds

(acceptance = 0.713, mean = 0.020448485912465595, variance = 1.0178967550807554)

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

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

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))

Precompiling packages...
   1098.8 ms  ✓ OrderedCollections
   1120.0 ms  ✓ DataStructures
   1234.7 ms  ✓ QuadGK
  3 dependencies successfully precompiled in 3 seconds

(exact = 0.027846124825536066, real_axis = 0.027846124825536056 + 1.3072501110915319e-17im, shifted = 0.0278461248256268 + 0.0im)

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

Gaussian: Phase Cancellation Plot

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)

Precompiling packages...
    316.2 ms  ✓ LaTeXStrings
    300.3 ms  ✓ Reexport
    347.1 ms  ✓ TensorCore
    324.5 ms  ✓ OpenLibm_jll
    351.9 ms  ✓ StatsAPI
    348.9 ms  ✓ Contour
    347.9 ms  ✓ Measures
    348.1 ms  ✓ StableRNGs
    363.2 ms  ✓ SuiteSparse_jll
    423.9 ms  ✓ DocStringExtensions
    472.2 ms  ✓ Serialization
    506.8 ms  ✓ Requires
    302.8 ms  ✓ SimpleBufferStream
    321.2 ms  ✓ PtrArrays
    333.9 ms  ✓ PCRE2_jll
    325.4 ms  ✓ DataAPI
    320.1 ms  ✓ BitFlags
    412.5 ms  ✓ DelimitedFiles
    456.8 ms  ✓ URIs
    429.8 ms  ✓ TranscodingStreams
    386.6 ms  ✓ SortingAlgorithms
    894.8 ms  ✓ Format
    336.1 ms  ✓ PrecompileTools
    382.8 ms  ✓ LoggingExtras
   1078.9 ms  ✓ Unzip
    465.0 ms  ✓ Graphite2_jll
    400.9 ms  ✓ EpollShim_jll
    450.1 ms  ✓ Libmount_jll
    461.0 ms  ✓ LLVMOpenMP_jll
   1308.7 ms  ✓ Grisu
    467.6 ms  ✓ Bzip2_jll
    436.6 ms  ✓ Xorg_libpciaccess_jll
              ✗ MbedTLS
    456.3 ms  ✓ Xorg_libICE_jll
    457.7 ms  ✓ Xorg_libXau_jll
    442.8 ms  ✓ libpng_jll
    923.0 ms  ✓ UnicodeFun
    475.8 ms  ✓ libfdk_aac_jll
   1148.9 ms  ✓ StructUtils
    450.2 ms  ✓ LERC_jll
    428.8 ms  ✓ fzf_jll
    551.6 ms  ✓ LAME_jll
    452.8 ms  ✓ Ogg_jll
    454.5 ms  ✓ mtdev_jll
   1534.6 ms  ✓ IrrationalConstants
    525.4 ms  ✓ JpegTurbo_jll
    531.1 ms  ✓ XZ_jll
    456.0 ms  ✓ Xorg_libXdmcp_jll
    475.1 ms  ✓ x265_jll
   1965.4 ms  ✓ MacroTools
    459.8 ms  ✓ x264_jll
    499.2 ms  ✓ libaom_jll
    444.2 ms  ✓ Xorg_xtrans_jll
    452.4 ms  ✓ Opus_jll
    497.0 ms  ✓ Zstd_jll
    474.5 ms  ✓ Expat_jll
    455.3 ms  ✓ Libffi_jll
   1830.0 ms  ✓ FixedPointNumbers
    550.8 ms  ✓ libevdev_jll
    470.5 ms  ✓ FriBidi_jll
    465.3 ms  ✓ Libuuid_jll
    501.7 ms  ✓ eudev_jll
    530.8 ms  ✓ GettextRuntime_jll
    362.5 ms  ✓ AliasTables
    351.0 ms  ✓ CodecZlib
    531.2 ms  ✓ ConcurrentUtilities
    307.2 ms  ✓ Showoff
    453.5 ms  ✓ Pixman_jll
    469.5 ms  ✓ FreeType2_jll
    488.0 ms  ✓ libdrm_jll
    474.5 ms  ✓ Xorg_libSM_jll
    538.7 ms  ✓ JLFzf
    466.7 ms  ✓ LogExpFunctions
   1033.1 ms  ✓ NaNMath
    507.9 ms  ✓ libvorbis_jll
    981.7 ms  ✓ RecipesBase
   1032.5 ms  ✓ Missings
    480.5 ms  ✓ Libtiff_jll
    476.2 ms  ✓ Dbus_jll
    549.7 ms  ✓ Wayland_jll
   1250.3 ms  ✓ OpenSSL
    900.2 ms  ✓ Ghostscript_jll
    466.4 ms  ✓ libinput_jll
    608.1 ms  ✓ Glib_jll
    613.3 ms  ✓ Fontconfig_jll
   1047.4 ms  ✓ Xorg_libxcb_jll
    454.5 ms  ✓ Xorg_xcb_util_jll
    478.5 ms  ✓ Xorg_libX11_jll
   1555.2 ms  ✓ ColorTypes
    470.0 ms  ✓ Xorg_xcb_util_keysyms_jll
    472.9 ms  ✓ Xorg_xcb_util_image_jll
    480.4 ms  ✓ Xorg_xcb_util_renderutil_jll
    463.2 ms  ✓ Xorg_libXrender_jll
    491.2 ms  ✓ Xorg_xcb_util_wm_jll
   2586.9 ms  ✓ Test
    474.4 ms  ✓ Xorg_libxkbfile_jll
    477.3 ms  ✓ Xorg_libXext_jll
    485.7 ms  ✓ Xorg_libXfixes_jll
    408.2 ms  ✓ ColorTypes → StyledStringsExt
    417.9 ms  ✓ ExceptionUnwrapping
    477.2 ms  ✓ Xorg_xcb_util_cursor_jll
    464.0 ms  ✓ Xorg_xkbcomp_jll
    469.4 ms  ✓ Xorg_libXinerama_jll
    471.1 ms  ✓ Xorg_libXrandr_jll
    469.3 ms  ✓ Xorg_libXcursor_jll
    472.1 ms  ✓ Xorg_libXi_jll
   3113.6 ms  ✓ SparseArrays
    529.9 ms  ✓ libva_jll
    563.9 ms  ✓ Cairo_jll
    619.2 ms  ✓ Libglvnd_jll
    413.1 ms  ✓ Xorg_xkeyboard_config_jll
    405.5 ms  ✓ Statistics → SparseArraysExt
    549.6 ms  ✓ HarfBuzz_jll
   2505.9 ms  ✓ Latexify
    497.0 ms  ✓ xkbcommon_jll
   1788.1 ms  ✓ ColorVectorSpace
    486.9 ms  ✓ libass_jll
    546.0 ms  ✓ Pango_jll
    521.9 ms  ✓ Latexify → SparseArraysExt
    476.3 ms  ✓ Vulkan_Loader_jll
              ✗ HTTP
    503.5 ms  ✓ libdecor_jll
    729.3 ms  ✓ FFMPEG_jll
   2739.3 ms  ✓ Colors
    518.2 ms  ✓ GLFW_jll
    835.7 ms  ✓ Qt6Base_jll
    394.4 ms  ✓ FFMPEG
   2129.2 ms  ✓ StatsBase
    490.3 ms  ✓ Qt6ShaderTools_jll
    517.3 ms  ✓ Qt6Svg_jll
    655.4 ms  ✓ GR_jll
   5774.8 ms  ✓ Parsers
   1457.3 ms  ✓ Qt6Declarative_jll
   2463.1 ms  ✓ ColorSchemes
    516.1 ms  ✓ Qt6Wayland_jll
   2668.0 ms  ✓ JSON
   2752.8 ms  ✓ GR
   4930.4 ms  ✓ PlotUtils
   1838.2 ms  ✓ PlotThemes
   2294.6 ms  ✓ RecipesPipeline
  26532.6 ms  ✓ Plots
  139 dependencies successfully precompiled in 43 seconds. 38 already precompiled.
Precompiling packages...
    866.5 ms  ✓ QuartoNotebookWorkerJSONExt (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
   1706.0 ms  ✓ QuartoNotebookWorkerPlotsExt (serial)
  1 dependency successfully precompiled in 2 seconds

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

Complex Gaussian Saddle

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

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

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

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

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.

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)

(zc = 0.0 + 1.0im, H = 2.0 + 0.0im, takagi_direction = 1.0 - 0.0im, quadratic_coefficient = 2.0 + 0.0im)

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

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

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

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

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

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))

Precompiling packages...
    461.9 ms  ✓ OpenSpecFun_jll
   1685.1 ms  ✓ SpecialFunctions
  2 dependencies successfully precompiled in 2 seconds. 11 already precompiled.
Precompiling packages...
    391.1 ms  ✓ ColorVectorSpace → SpecialFunctionsExt
  1 dependency successfully precompiled in 0 seconds. 20 already precompiled.

(zplus = 0.0 + 1.0im, S_zplus = 0.6666666666666667 + 0.0im, zminus = 0.0 - 1.0im, S_zminus = -0.6666666666666667 - 0.0im, exact = 0.8500673223499207)

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

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_flow_rhs! (generic function with 1 method)

Airy: Runnable Finite Flow Cell

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)

Precompiling packages...
    302.5 ms  ✓ IteratorInterfaceExtensions
    301.9 ms  ✓ Future
    302.7 ms  ✓ SafeTestsets
    308.0 ms  ✓ CommonSolve
    319.0 ms  ✓ MuladdMacro
    325.6 ms  ✓ ConcreteStructs
    323.9 ms  ✓ CompositionsBase
    355.6 ms  ✓ EnumX
    360.4 ms  ✓ ExprTools
    371.8 ms  ✓ InverseFunctions
    376.5 ms  ✓ FastPower
    476.6 ms  ✓ EnzymeCore
    486.7 ms  ✓ ADTypes
    317.4 ms  ✓ FastClosures
    356.7 ms  ✓ Adapt
    356.6 ms  ✓ StaticArraysCore
    370.2 ms  ✓ SciMLPublic
    401.6 ms  ✓ CommonSubexpressions
    403.0 ms  ✓ TruncatedStacktraces
    412.5 ms  ✓ SciMLLogging
    536.4 ms  ✓ OrdinaryDiffEqRosenbrockTableaus
    483.6 ms  ✓ DiffRules
    328.3 ms  ✓ InverseFunctions → InverseFunctionsDatesExt
    307.2 ms  ✓ LogExpFunctions → LogExpFunctionsInverseFunctionsExt
    601.3 ms  ✓ LazyArtifacts
    313.7 ms  ✓ CompositionsBase → CompositionsBaseInverseFunctionsExt
    320.2 ms  ✓ ADTypes → ADTypesEnzymeCoreExt
   1116.9 ms  ✓ ConstructionBase
    411.6 ms  ✓ Adapt → AdaptSparseArraysExt
   1239.2 ms  ✓ FunctionWrappers
    308.1 ms  ✓ EnzymeCore → AdaptExt
    418.9 ms  ✓ GPUArraysCore
    329.4 ms  ✓ DiffResults
    310.8 ms  ✓ ConstructionBase → ConstructionBaseLinearAlgebraExt
    321.3 ms  ✓ ADTypes → ADTypesConstructionBaseExt
    760.6 ms  ✓ IntelOpenMP_jll
   1288.7 ms  ✓ RuntimeGeneratedFunctions
    793.1 ms  ✓ oneTBB_jll
   1824.5 ms  ✓ Distributed
   1097.2 ms  ✓ ArrayInterface
   1327.9 ms  ✓ DifferentiationInterface
    800.2 ms  ✓ Setfield
    304.7 ms  ✓ ArrayInterface → ArrayInterfaceStaticArraysCoreExt
    359.6 ms  ✓ ArrayInterface → ArrayInterfaceGPUArraysCoreExt
    304.3 ms  ✓ DifferentiationInterface → DifferentiationInterfaceGPUArraysCoreExt
    415.8 ms  ✓ FastBroadcast
    419.0 ms  ✓ ArrayInterface → ArrayInterfaceSparseArraysExt
    447.8 ms  ✓ MaybeInplace
   2044.5 ms  ✓ TimerOutputs
   1313.2 ms  ✓ FunctionWrappersWrappers
    843.0 ms  ✓ MKL_jll
    427.2 ms  ✓ DifferentiationInterface → DifferentiationInterfaceSparseArraysExt
    465.4 ms  ✓ FiniteDiff
    832.2 ms  ✓ PreallocationTools
    434.0 ms  ✓ MaybeInplace → MaybeInplaceSparseArraysExt
   1064.3 ms  ✓ SciMLStructures
    363.9 ms  ✓ DifferentiationInterface → DifferentiationInterfaceFiniteDiffExt
    416.5 ms  ✓ FiniteDiff → FiniteDiffSparseArraysExt
   1979.8 ms  ✓ Accessors
    381.2 ms  ✓ Accessors → LinearAlgebraExt
   2796.9 ms  ✓ ForwardDiff
   3397.4 ms  ✓ SparseMatrixColorings
    363.6 ms  ✓ FastPower → FastPowerForwardDiffExt
   4306.6 ms  ✓ Krylov
    557.6 ms  ✓ DifferentiationInterface → DifferentiationInterfaceSparseMatrixColoringsExt
    713.0 ms  ✓ PreallocationTools → PreallocationToolsForwardDiffExt
   1105.1 ms  ✓ DifferentiationInterface → DifferentiationInterfaceForwardDiffExt
   1463.0 ms  ✓ SymbolicIndexingInterface
   1885.5 ms  ✓ SciMLOperators
    359.7 ms  ✓ SciMLOperators → SciMLOperatorsStaticArraysCoreExt
    471.5 ms  ✓ SciMLOperators → SciMLOperatorsSparseArraysExt
   3029.1 ms  ✓ RecursiveArrayTools
    394.7 ms  ✓ RecursiveArrayTools → RecursiveArrayToolsStatisticsExt
    406.1 ms  ✓ RecursiveArrayTools → RecursiveArrayToolsFastBroadcastExt
    430.7 ms  ✓ RecursiveArrayTools → RecursiveArrayToolsForwardDiffExt
    520.1 ms  ✓ RecursiveArrayTools → RecursiveArrayToolsSparseArraysExt
   5986.5 ms  ✓ SciMLBase
    571.8 ms  ✓ SciMLBase → SciMLBaseDifferentiationInterfaceExt
   1399.8 ms  ✓ SciMLBase → SciMLBaseForwardDiffExt
   2498.1 ms  ✓ SciMLJacobianOperators
   4962.7 ms  ✓ LinearSolve
   1778.6 ms  ✓ LinearSolve → LinearSolveForwardDiffExt
   3831.6 ms  ✓ LineSearch
   2371.8 ms  ✓ LinearSolve → LinearSolveSparseArraysExt
   4952.3 ms  ✓ NonlinearSolveBase
    978.9 ms  ✓ LinearSolve → LinearSolveEnzymeExt
    713.9 ms  ✓ NonlinearSolveBase → NonlinearSolveBaseLineSearchExt
    835.6 ms  ✓ NonlinearSolveBase → NonlinearSolveBaseSparseMatrixColoringsExt
    880.2 ms  ✓ NonlinearSolveBase → NonlinearSolveBaseSparseArraysExt
   1442.2 ms  ✓ NonlinearSolveBase → NonlinearSolveBaseLinearSolveExt
   2582.1 ms  ✓ BracketingNonlinearSolve
   2952.9 ms  ✓ NonlinearSolveBase → NonlinearSolveBaseForwardDiffExt
    829.0 ms  ✓ BracketingNonlinearSolve → BracketingNonlinearSolveForwardDiffExt
   2174.6 ms  ✓ DiffEqBase
   4725.9 ms  ✓ NonlinearSolveSpectralMethods
   1030.5 ms  ✓ DiffEqBase → DiffEqBaseSparseArraysExt
    834.2 ms  ✓ NonlinearSolveSpectralMethods → NonlinearSolveSpectralMethodsForwardDiffExt
   2357.1 ms  ✓ DiffEqBase → DiffEqBaseForwardDiffExt
   3090.9 ms  ✓ OrdinaryDiffEqCore
   4918.1 ms  ✓ SimpleNonlinearSolve
    931.1 ms  ✓ OrdinaryDiffEqCore → OrdinaryDiffEqCoreSparseArraysExt
   7913.2 ms  ✓ NonlinearSolveQuasiNewton
   1046.7 ms  ✓ NonlinearSolveQuasiNewton → NonlinearSolveQuasiNewtonForwardDiffExt
   3801.6 ms  ✓ OrdinaryDiffEqTsit5
   3015.8 ms  ✓ OrdinaryDiffEqDifferentiation
   1010.6 ms  ✓ OrdinaryDiffEqDifferentiation → OrdinaryDiffEqDifferentiationSparseArraysExt
  15436.9 ms  ✓ NonlinearSolveFirstOrder
  12514.9 ms  ✓ OrdinaryDiffEqVerner
   8002.1 ms  ✓ OrdinaryDiffEqRosenbrock
  20220.4 ms  ✓ NonlinearSolve
   3374.8 ms  ✓ OrdinaryDiffEqNonlinearSolve
   7471.1 ms  ✓ OrdinaryDiffEqSDIRK
   5381.2 ms  ✓ OrdinaryDiffEqBDF
  11405.9 ms  ✓ OrdinaryDiffEqDefault
   1871.4 ms  ✓ OrdinaryDiffEq
   1909.8 ms  ✓ DifferentialEquations
  116 dependencies successfully precompiled in 93 seconds. 51 already precompiled.
Precompiling packages...
    324.7 ms  ✓ StructUtils → StructUtilsStaticArraysCoreExt
  1 dependency successfully precompiled in 0 seconds. 6 already precompiled.
Precompiling packages...
    628.0 ms  ✓ SparseMatrixColorings → SparseMatrixColoringsColorsExt
  1 dependency successfully precompiled in 1 seconds. 18 already precompiled.

(0.4055158488524138 + 1.0238351339876182im, 2.0553868474538546 + 0.23946842282017358im)

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

Airy: Plot the Flowed Contour

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

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]

4-element Vector{@NamedTuple{T::Float64, estimate::ComplexF64, error::Float64}}:
 (T = 0.0, estimate = 0.8186488540408526 - 6.938893903907228e-18im, error = 0.031418468309068026)
 (T = 0.1, estimate = 0.8417343457903285 + 0.0im, error = 0.008332976559592198)
 (T = 0.2, estimate = 0.8495218320187008 + 4.163336342344337e-17im, error = 0.0005454903312198667)
 (T = 0.35, estimate = 0.8500673221425528 + 6.938893903907228e-18im, error = 2.073679006286967e-10)

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)|. \]

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))

(acceptance = 0.6096666666666667, mean_tangent = -0.008843771899084503, std_tangent = 0.2986172803543729)

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))}. \]

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))

(average_phase = 0.9830826101470839 - 0.005260903257862508im, abs_average_phase = 0.9830966867387418)

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

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))

(mean = -0.003760121663699721, std = 0.9739459597389546, tau = 1.029268754983147, ess = 4857.817723303831)

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

MCMC Diagnostics: Plot the Chain

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))

(acceptance = 0.7085, mean = 1.3565169377133028, std = 1.8019059418886119)

MCMC Diagnostics: Trace, Histogram, ACF

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.

Geometry schematic redrawn from the construction in Fukuma-Matsumoto, arXiv: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

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)]

3-element Vector{@NamedTuple{label::String, Re_z::Float64, Im_z::Float64, Re_S::Float64, Im_S::Float64}}:
 (label = "left", Re_z = -1.1239, Im_z = -0.2961, Re_S = 0.1338, Im_S = -0.8407)
 (label = "middle", Re_z = 0.0, Im_z = 0.5923, Re_S = -0.0177, Im_S = 0.0)
 (label = "right", Re_z = 1.1239, Im_z = -0.2961, Re_S = 0.1338, Im_S = 0.8407)

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

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?

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

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

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

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

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.

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

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.

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.

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

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)]

5-element Vector{@NamedTuple{K::Int64, cond::Float64}}:
 (K = 1, cond = 8.886110520507984e6)
 (K = 2, cond = 2981.4582805786117)
 (K = 4, cond = 55.41958459646888)
 (K = 8, cond = 8.414658346931526)
 (K = 16, cond = 4.035065540210353)

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

Multiple Shooting: Conditioning Plot

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

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

(A_right_left = 0.0 + 1.6814478064077518im, stokes_angle_radians = 1.5707963267948966, stokes_angle_over_pi = 0.5)

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

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

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

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)

([(N = 0, estimate = 2.5066282746310002), (N = 1, estimate = 2.2810317299142104), (N = 2, estimate = 2.399469915890525), (N = 3, estimate = 2.282216111773973), (N = 4, estimate = 2.45369980029443), (N = 5, estimate = 2.121364411941785), (N = 6, estimate = 2.9239543748134227), (N = 7, estimate = 0.6021762679347563), (N = 8, estimate = 8.429470710749461)], 2.3467129309223336)

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). \]

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))

(exact = 2.3467129309223336, truncation_N8 = 8.429470710749461, borel_pade = 2.346714183716724)

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

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]))

(action = 1.9997142290558227, endpoints = (-0.999329299739067, 0.999329299739067))

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. \]

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))

(initial = 4.516666623999998, relaxed = 2.3328583894145565)

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

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)]

4-element Vector{@NamedTuple{N::Int64, phase::ComplexF64}}:
 (N = 4, phase = 0.7224016417296388 - 0.010651149814091493im)
 (N = 8, phase = 0.5257318763121992 - 0.013558221406814336im)
 (N = 16, phase = 0.280019015451821 + 0.00033492002894137946im)
 (N = 32, phase = 0.07119456932164227 - 0.00994114021796093im)

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

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))

(raw_error = 0.008062420900239342, shanks_error = 2.311610621807958e-6)

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

Spinfoam Data Bridge

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

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))

Precompiling packages...
   4097.1 ms  ✓ JSON3
  1 dependency successfully precompiled in 4 seconds. 10 already precompiled.
Precompiling packages...
    859.7 ms  ✓ QuartoNotebookWorkerJSON3Ext (serial)
  1 dependency successfully precompiled in 1 seconds

(tetrahedra = 5, angle_table_shape = (5, 5, 2), area_blocks = 2, first_angle_pair = [1.5707963267948966, 1.5707963267948966], kappa_blocks = 2)

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.