Asymptotic Bootstrap in Matrix Integrals

• Asymptotic Bootstrap Method is a non-perturbative technique that reconstructs observables in unitary matrix models using exact recursion relations and asymptotic tail conditions.
• It achieves high numerical precision by linearly propagating a finite boundary data set across infinitely many modes without relying on positivity constraints.
• The method efficiently solves banded linear systems to accurately compute key quantities, such as Wilson loops and instanton effects, even at complex couplings.

Searching arXiv for the cited work and closely related papers. The asymptotic bootstrap method is a non-perturbative computational procedure for reconstructing observables of matrix integrals from exact recursion relations supplemented by asymptotic information at large mode number. In the formulation developed for unitary matrix integrals at complex coupling, the method is a linear truncation-and-tail-control scheme that combines exact finite-difference recursions for Fourier moments with asymptotic control of their rapidly decaying tails, thereby producing very high numerical precision without relying on positivity or semidefinite programming (Berenstein et al., 20 Feb 2026). In this sense it is a “bootstrap” because it propagates a finite amount of boundary data across infinitely many modes by means of exact identities, but it differs fundamentally from positivity-based bootstrap programs: it remains applicable when the coupling is complex and positivity is absent (Berenstein et al., 20 Feb 2026). Closely related antecedents include the “shoestring bootstrap” for circle problems (Berenstein et al., 28 Feb 2025), while a parallel development for Hermitian matrix models shows analogous exponentially fast convergence in a cone of complex couplings (Berenstein et al., 29 Jun 2026).

1. Definition and conceptual role

In unitary matrix models, the asymptotic bootstrap estimate is designed to reconstruct Fourier moments of the one-eigenvalue weight and, from them, orthogonal polynomials on the unit circle (OPUC), partition functions, and Wilson loops (Berenstein et al., 20 Feb 2026). The method starts from exact Schwinger–Dyson or loop-equation recursions for the moments, truncates the infinite system at a large mode $M$, and closes the system by imposing an asymptotic tail condition informed by the exact large-mode decay rate (Berenstein et al., 20 Feb 2026). This yields a banded linear system whose solution determines the low modes with exponentially small truncation error.

Its distinctive feature is independence from positivity. Positivity or semidefinite programming (SDP) methods constrain moment matrices through positive semidefiniteness, but such constraints become unavailable or inapplicable at complex ’t Hooft coupling (Berenstein et al., 20 Feb 2026). The asymptotic bootstrap instead exploits the dichotomy between the physical decaying solution of the recursion and unphysical growing solutions. This suggests that the essential non-perturbative information is encoded in selecting the decaying branch by appropriate boundary conditions at infinity.

The method is therefore both numerical and structural. Numerically, it is a fast linear solver for moments and observables. Structurally, it identifies the physically admissible solution among all formal recursion solutions by enforcing the correct asymptotic sector (Berenstein et al., 20 Feb 2026). A plausible implication is that the method belongs to a broader class of bootstrap procedures in which asymptotic admissibility replaces positivity as the primary selection principle.

2. Matrix-integral setting and exact recursion structure

The unitary models considered are single-trace unitary matrix integrals with $U(N)$-Haar measure,

$$Z_N(\boldsymbol t) = \frac{1}{\mathrm{Vol}(U(N))}\int_{U(N)} dU\, \exp\!\Big(V(U)\Big),$$

with potential

$$V(U)= \sum_{k=1}^{K} \frac{t_k}{k}\,\Big(\mathrm{Tr}\, U^k + \mathrm{Tr}\, U^{-k}\Big),$$

where the couplings $t_k$ may be complex (Berenstein et al., 20 Feb 2026). Two emphasized examples are the Gross–Witten–Wadia model, $V(U)=t(\mathrm{Tr}\,U+\mathrm{Tr}\,U^{-1})$, and a quadratic single-trace model with first and second harmonics (Berenstein et al., 20 Feb 2026).

After diagonalization $U=\mathrm{diag}(e^{i\theta_1},\dots,e^{i\theta_N})$, the basic objects become Fourier moments of the one-eigenvalue weight

$$a_n=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\,e^{i n \theta}\, e^{V(e^{i\theta})},$$

which determine Toeplitz moment matrices and OPUC (Berenstein et al., 20 Feb 2026). With monic OPUC $p_n(z)$ defined by

$$\langle f|g\rangle=\int_{S^1}\frac{dz}{2\pi i z}\, e^{V(z)}\, g(z)\,f(z^{-1}), \qquad \langle p_m|p_n\rangle=h_n\,\delta_{mn},$$

one has

$U(N)$0

so the moment sequence is the computational backbone of the full finite-$U(N)$1 problem (Berenstein et al., 20 Feb 2026).

The exact recursions arise from integration by parts on the circle. Writing

$U(N)$2

the moments satisfy, for all integers $U(N)$3,

$U(N)$4

(Berenstein et al., 20 Feb 2026). Thus $U(N)$5 and the $U(N)$6 seeds $U(N)$7 determine the full sequence. For $U(N)$8,

$U(N)$9

and for $$Z_N(\boldsymbol t) = \frac{1}{\mathrm{Vol}(U(N))}\int_{U(N)} dU\, \exp\!\Big(V(U)\Big),$$0,

$$Z_N(\boldsymbol t) = \frac{1}{\mathrm{Vol}(U(N))}\int_{U(N)} dU\, \exp\!\Big(V(U)\Big),$$1

(Berenstein et al., 20 Feb 2026). These identities are exact and do not invoke positivity.

3. Large-mode asymptotics and tail control

The decisive input is the asymptotic structure of the recursion at large mode number. For nondegenerate highest coupling $$Z_N(\boldsymbol t) = \frac{1}{\mathrm{Vol}(U(N))}\int_{U(N)} dU\, \exp\!\Big(V(U)\Big),$$2, the physical decaying solution satisfies

$$Z_N(\boldsymbol t) = \frac{1}{\mathrm{Vol}(U(N))}\int_{U(N)} dU\, \exp\!\Big(V(U)\Big),$$3

where the coefficients $$Z_N(\boldsymbol t) = \frac{1}{\mathrm{Vol}(U(N))}\int_{U(N)} dU\, \exp\!\Big(V(U)\Big),$$4 do not grow with $$Z_N(\boldsymbol t) = \frac{1}{\mathrm{Vol}(U(N))}\int_{U(N)} dU\, \exp\!\Big(V(U)\Big),$$5 (Berenstein et al., 20 Feb 2026). The decay is therefore faster than exponential. By contrast, a basis of growing solutions behaves as

$$Z_N(\boldsymbol t) = \frac{1}{\mathrm{Vol}(U(N))}\int_{U(N)} dU\, \exp\!\Big(V(U)\Big),$$6

(Berenstein et al., 20 Feb 2026). This separation between rapidly decaying and rapidly growing sectors is the central asymptotic mechanism.

Because a small seed error $$Z_N(\boldsymbol t) = \frac{1}{\mathrm{Vol}(U(N))}\int_{U(N)} dU\, \exp\!\Big(V(U)\Big),$$7 excites the growing solution with amplification roughly proportional to

$$Z_N(\boldsymbol t) = \frac{1}{\mathrm{Vol}(U(N))}\int_{U(N)} dU\, \exp\!\Big(V(U)\Big),$$8

the boundary condition imposed at large $$Z_N(\boldsymbol t) = \frac{1}{\mathrm{Vol}(U(N))}\int_{U(N)} dU\, \exp\!\Big(V(U)\Big),$$9 determines low-mode accuracy (Berenstein et al., 20 Feb 2026). Two truncation strategies are contrasted. A “cheap truncation” merely enforces boundedness of the tail, for example $$V(U)= \sum_{k=1}^{K} \frac{t_k}{k}\,\Big(\mathrm{Tr}\, U^k + \mathrm{Tr}\, U^{-k}\Big),$$0, which leads to low-mode error scaling

$$V(U)= \sum_{k=1}^{K} \frac{t_k}{k}\,\Big(\mathrm{Tr}\, U^k + \mathrm{Tr}\, U^{-k}\Big),$$1

The asymptotic bootstrap, or “shoestring,” instead sets

$$V(U)= \sum_{k=1}^{K} \frac{t_k}{k}\,\Big(\mathrm{Tr}\, U^k + \mathrm{Tr}\, U^{-k}\Big),$$2

matching the exponentially small physical tail, and produces a low-mode error parametrically equal to the square of the cheap one,

$$V(U)= \sum_{k=1}^{K} \frac{t_k}{k}\,\Big(\mathrm{Tr}\, U^k + \mathrm{Tr}\, U^{-k}\Big),$$3

(Berenstein et al., 20 Feb 2026).

This is the formal basis for the observed gain in precision. The low-mode estimates typically acquire about twice as many correct digits as under the cheap truncation (Berenstein et al., 20 Feb 2026). Pre-asymptotic drift can occur because subleading recursion terms are only suppressed by fractional inverse powers of $$V(U)= \sum_{k=1}^{K} \frac{t_k}{k}\,\Big(\mathrm{Tr}\, U^k + \mathrm{Tr}\, U^{-k}\Big),$$4, but the exponential rate predicted by the large-mode analysis remains correct (Berenstein et al., 20 Feb 2026).

This tail-control viewpoint was foreshadowed by the “shoestring bootstrap” for distributions on the circle, which similarly used the asymptotic exponential smallness of high Fourier modes to bypass a full SDP search (Berenstein et al., 28 Feb 2025). The unitary-matrix asymptotic bootstrap turns that observation into a systematic finite-$$V(U)= \sum_{k=1}^{K} \frac{t_k}{k}\,\Big(\mathrm{Tr}\, U^k + \mathrm{Tr}\, U^{-k}\Big),$$5 algorithm (Berenstein et al., 20 Feb 2026).

4. Algorithmic realization

The implementation is entirely linear. One fixes the number $$V(U)= \sum_{k=1}^{K} \frac{t_k}{k}\,\Big(\mathrm{Tr}\, U^k + \mathrm{Tr}\, U^{-k}\Big),$$6 of nonzero couplings and a large cutoff $$V(U)= \sum_{k=1}^{K} \frac{t_k}{k}\,\Big(\mathrm{Tr}\, U^k + \mathrm{Tr}\, U^{-k}\Big),$$7. The procedure begins with $$V(U)= \sum_{k=1}^{K} \frac{t_k}{k}\,\Big(\mathrm{Tr}\, U^k + \mathrm{Tr}\, U^{-k}\Big),$$8 and the symmetry $$V(U)= \sum_{k=1}^{K} \frac{t_k}{k}\,\Big(\mathrm{Tr}\, U^k + \mathrm{Tr}\, U^{-k}\Big),$$9, imposes asymptotic tail boundary conditions $t_k$0, and solves the exact recursion over modes $t_k$1 as a banded linear system (Berenstein et al., 20 Feb 2026). Because the recursion bandwidth is $t_k$2, the solve scales as $t_k$3 or $t_k$4, with small constants (Berenstein et al., 20 Feb 2026).

Once the moments are obtained, one forms the Toeplitz matrix $t_k$5, with $t_k$6, and computes the monic OPUC by Gram–Schmidt (Berenstein et al., 20 Feb 2026). This determines $t_k$7, the three-term recurrences if desired, and Wilson loops through the overlap formula. The naive complexity of OPUC construction up to degree $t_k$8 is $t_k$9, though Toeplitz or OPUC structure can improve this (Berenstein et al., 20 Feb 2026).

Convergence is monitored by increasing $V(U)=t(\mathrm{Tr}\,U+\mathrm{Tr}\,U^{-1})$0 until changes in $V(U)=t(\mathrm{Tr}\,U+\mathrm{Tr}\,U^{-1})$1 or in $V(U)=t(\mathrm{Tr}\,U+\mathrm{Tr}\,U^{-1})$2 fall below a target tolerance (Berenstein et al., 20 Feb 2026). The asymptotic error model predicts that the required precision is typically reached once $V(U)=t(\mathrm{Tr}\,U+\mathrm{Tr}\,U^{-1})$3 (Berenstein et al., 20 Feb 2026). Stability diagnostics include consistency under $V(U)=t(\mathrm{Tr}\,U+\mathrm{Tr}\,U^{-1})$4, comparison to the asymptotic error envelope

$V(U)=t(\mathrm{Tr}\,U+\mathrm{Tr}\,U^{-1})$5

and conditioning estimates from eigenvalue ratios of cumulative transfer matrices $V(U)=t(\mathrm{Tr}\,U+\mathrm{Tr}\,U^{-1})$6 (Berenstein et al., 20 Feb 2026).

The method is particularly amenable to high-precision arithmetic. In the worked example $V(U)=t(\mathrm{Tr}\,U+\mathrm{Tr}\,U^{-1})$7, the computation achieved more than $V(U)=t(\mathrm{Tr}\,U+\mathrm{Tr}\,U^{-1})$8 correct digits at $V(U)=t(\mathrm{Tr}\,U+\mathrm{Tr}\,U^{-1})$9 (Berenstein et al., 20 Feb 2026). The paper also notes the practical use of rational couplings $U=\mathrm{diag}(e^{i\theta_1},\dots,e^{i\theta_N})$0 to avoid roundoff (Berenstein et al., 20 Feb 2026).

5. Wilson loops, instantons, and non-perturbative resolution

A central application is the computation of Wilson loop expectation values with sensitivity to exponentially small instanton corrections. In the strong-coupling limit with fixed $U=\mathrm{diag}(e^{i\theta_1},\dots,e^{i\theta_N})$1 as $U=\mathrm{diag}(e^{i\theta_1},\dots,e^{i\theta_N})$2, perturbation theory predicts

$U=\mathrm{diag}(e^{i\theta_1},\dots,e^{i\theta_N})$3

equivalently $U=\mathrm{diag}(e^{i\theta_1},\dots,e^{i\theta_N})$4 in the ’t Hooft scaling $U=\mathrm{diag}(e^{i\theta_1},\dots,e^{i\theta_N})$5 (Berenstein et al., 20 Feb 2026). The moments reconstructed by asymptotic bootstrap and then fed into OPUC reproduce these large-$U=\mathrm{diag}(e^{i\theta_1},\dots,e^{i\theta_N})$6 limits with residuals decaying as $U=\mathrm{diag}(e^{i\theta_1},\dots,e^{i\theta_N})$7, enabling direct numerical access to instanton sectors (Berenstein et al., 20 Feb 2026).

In the ungapped phase, the partition function has perturbative part

$U=\mathrm{diag}(e^{i\theta_1},\dots,e^{i\theta_N})$8

and non-perturbative completion

$U=\mathrm{diag}(e^{i\theta_1},\dots,e^{i\theta_N})$9

with $$a_n=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\,e^{i n \theta}\, e^{V(e^{i\theta})},$$0 labeling tunneling eigenvalue or anti-eigenvalue pairs (Berenstein et al., 20 Feb 2026). For the Gross–Witten–Wadia model in the ungapped regime $$a_n=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\,e^{i n \theta}\, e^{V(e^{i\theta})},$$1, the two saddle points are

$$a_n=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\,e^{i n \theta}\, e^{V(e^{i\theta})},$$2

and the instanton action is

$$a_n=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\,e^{i n \theta}\, e^{V(e^{i\theta})},$$3

(Berenstein et al., 20 Feb 2026). The leading one-instanton contribution obeys

$$a_n=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\,e^{i n \theta}\, e^{V(e^{i\theta})},$$4

with explicit coefficients $$a_n=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\,e^{i n \theta}\, e^{V(e^{i\theta})},$$5 given in the paper (Berenstein et al., 20 Feb 2026).

Wilson loops inherit a trans-series expansion,

$$a_n=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\,e^{i n \theta}\, e^{V(e^{i\theta})},$$6

and the numerical bootstrap confirms this structure by the ratio

$$a_n=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\,e^{i n \theta}\, e^{V(e^{i\theta})},$$7

for which $$a_n=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\,e^{i n \theta}\, e^{V(e^{i\theta})},$$8 as $$a_n=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\,e^{i n \theta}\, e^{V(e^{i\theta})},$$9 increases, with systematic improvement upon inclusion of successive $p_n(z)$0 corrections (Berenstein et al., 20 Feb 2026). The same pattern persists in the quadratic model, even when multiple instanton channels interfere (Berenstein et al., 20 Feb 2026).

This instanton sensitivity is one of the defining achievements of the method. Because the truncation error is itself exponentially small, the numerics can resolve exponentially tiny non-perturbative contributions rather than merely reproduce perturbative data (Berenstein et al., 20 Feb 2026).

6. Complex coupling, phase structure, and OPUC zeros

The absence of positivity becomes especially consequential in the complex ’t Hooft plane. The asymptotic bootstrap remains applicable there because it never required positivity in the first place (Berenstein et al., 20 Feb 2026). In this setting the ungapped-phase instanton action $p_n(z)$1 defines two distinguished loci. Anti-Stokes lines satisfy $p_n(z)$2, and Stokes lines satisfy $p_n(z)$3 (Berenstein et al., 20 Feb 2026). Anti-Stokes lines track where eigenvalue transport becomes energetically neutral and accurately follow phase boundaries, while Stokes lines signal changes in Lefschetz-thimble decomposition. The unitary-matrix analysis finds that some Stokes lines also serve as useful proxies for additional phase boundaries, especially near critical points (Berenstein et al., 20 Feb 2026).

Phase diagrams are mapped numerically by the zeros of $p_n(z)$4 at large $p_n(z)$5. These zeros accumulate on the planar support of the eigenvalue spectral density in the ’t Hooft limit, and plotting them reveals the topology of the spectral support, such as the number of cuts (Berenstein et al., 20 Feb 2026). The OPUC–density correspondence is rigorous in Hermitian settings and is reported to work empirically and robustly for the unitary models studied (Berenstein et al., 20 Feb 2026).

For the GWW model, the method reproduces the known complex-coupling phase diagram: an ungapped phase, two one-cut phases centered at $p_n(z)$6, and a two-cut phase (Berenstein et al., 20 Feb 2026). Anti-Stokes lines arise from explicit actions

$p_n(z)$7

and

$p_n(z)$8

(Berenstein et al., 20 Feb 2026). For the quadratic model, the observed phases include ungapped, one-, two-, three-, and four-cut regions, inferred from OPUC zero distributions guided by anti-Stokes and complementary Stokes lines (Berenstein et al., 20 Feb 2026).

A plausible implication is that the asymptotic bootstrap furnishes a non-perturbative mapping from exact recursion data to spectral geometry in coupling space, with OPUC zeros serving as the practical interface between finite-$p_n(z)$9 computation and planar phase structure.

7. Relation to other bootstrap traditions, advantages, and limitations

The term “asymptotic bootstrap method” has multiple meanings across disciplines, but the matrix-integral version is technically distinct from statistical resampling bootstraps in econometrics or semiparametric inference. Here the “bootstrap” refers to recursive propagation of exact identities closed by asymptotic data, not to sampling with replacement. Its nearest relatives are circle and Hermitian moment bootstrap methods that likewise use asymptotic tail behavior to replace or supplement positivity-based feasibility conditions (Berenstein et al., 28 Feb 2025, Berenstein et al., 29 Jun 2026).

Relative to positivity or SDP bootstrap approaches, the asymptotic bootstrap for unitary matrix integrals is linear, extremely fast, and applicable at complex couplings (Berenstein et al., 20 Feb 2026). It produces high-precision estimates rather than certified bounds, although positivity checks on Toeplitz matrices can be performed a posteriori when the couplings are real (Berenstein et al., 20 Feb 2026). Relative to exact integrability or resurgence methods, it is numerical rather than symbolic: saddle-point and instanton data are used for validation and interpretation, not for the primary computation (Berenstein et al., 20 Feb 2026).

The method also has explicit scope conditions. The analysis exploits the symmetry $$\langle f|g\rangle=\int_{S^1}\frac{dz}{2\pi i z}\, e^{V(z)}\, g(z)\,f(z^{-1}), \qquad \langle p_m|p_n\rangle=h_n\,\delta_{mn},$$0, which implies $$\langle f|g\rangle=\int_{S^1}\frac{dz}{2\pi i z}\, e^{V(z)}\, g(z)\,f(z^{-1}), \qquad \langle p_m|p_n\rangle=h_n\,\delta_{mn},$$1 and permits a one-sided tail condition (Berenstein et al., 20 Feb 2026). Fully generic complex potentials with $$\langle f|g\rangle=\int_{S^1}\frac{dz}{2\pi i z}\, e^{V(z)}\, g(z)\,f(z^{-1}), \qquad \langle p_m|p_n\rangle=h_n\,\delta_{mn},$$2 would require simultaneous control of positive and negative tails, a case not implemented in the paper (Berenstein et al., 20 Feb 2026). Precision also depends on entering the asymptotic regime $$\langle f|g\rangle=\int_{S^1}\frac{dz}{2\pi i z}\, e^{V(z)}\, g(z)\,f(z^{-1}), \qquad \langle p_m|p_n\rangle=h_n\,\delta_{mn},$$3 and on avoiding extreme hierarchies among the couplings $$\langle f|g\rangle=\int_{S^1}\frac{dz}{2\pi i z}\, e^{V(z)}\, g(z)\,f(z^{-1}), \qquad \langle p_m|p_n\rangle=h_n\,\delta_{mn},$$4 (Berenstein et al., 20 Feb 2026).

Finally, the method should not be conflated with positivity-free asymptotic bootstrap schemes in unrelated contexts unless the underlying mechanism is specified. In the matrix-model literature, the defining mechanism is the imposition of asymptotically correct decaying tails on exact moment recursions. That mechanism first appeared in simplified circle bootstrap problems (Berenstein et al., 28 Feb 2025), was developed here for unitary matrix integrals (Berenstein et al., 20 Feb 2026), and has since been mirrored for Hermitian matrix models, where convergence is asymptotically exponential in a coupling cone determined by the first subleading coupling (Berenstein et al., 29 Jun 2026). This suggests a broader methodological principle: when exact recursions admit both physical decaying and unphysical growing branches, asymptotic selection can replace positivity as the organizing principle of non-perturbative bootstrap computation.