A stochastic PDE approach to large N problems in quantum field theory: A survey

In quantum field theory (QFT), there are important models where the “field” takes values in an N-dimensional space. A natural question is that what is the asymptotic behavior of such models as the dimensionality of the target space N → ∞. The perspective of investigating these questions in this Review is called stochastic quantization. Stochastic quantization is a procedure that formulates a quantum field theory in terms of a stochastic partial differential equations (SPDEs). More precisely, for a quantum field theory where the target space has dimension N, its stochastic quantization yields a system of N coupled SPDEs. This procedure, then, brings in many powerful analytical tools, including the recent developments of SPDE solution theories and their a priori estimates, and mean field theory techniques. The purpose of this Review is to survey the recent results along this direction.

For most of this Review, we will focus on the linear σ-model in two and three dimensions. However, we will first briefly review the origin of the large N problems in physics.

Large N methods (or “1/N expansions”) in theoretical physics are ubiquitous and are generally applied to models where the dimensionality of the target space is large. It was first noted by Stanley¹ for lattice spin models (rather than QFT) that they exhibit considerable simplification as N (the number of spin components) becomes large: at least at the level of free energy, as N → ∞, these models become the “spherical model,” which is a solvable model introduced by Berlin and Kac in 1952.² See Baxter’s book,³ Sec. 5, for an expository discussion on spherical model. Later, the large N convergence to the spherical model was established at the level of free energies⁴ and, then, of correlations.^5,6 Again, in the context of spin systems, the authors of Ref. 7 found a systematic way to expand in powers of 1/N and obtained 1/N order corrections to the spherical model.

In 1973, Wilson exploited this idea in quantum field theory.⁸ In particular, Wilson studied both bosonic (linear σ-model) and fermionic models. The linear σ-model is the N-component generalization of the $Φd4$ model, given by the (formal) measure

over $RN$ valued fields Φ = (Φ₁, Φ₂, …, Φ_N), where Z_N is a normalization constant and $DΦ$ is the formal “Lebesgue measure.” We will discuss below the physics predictions and how they are related with our rigorous results. Soon later, Coleman, Jackiw, and Politzer⁹ calculated the effective potential and studied spontaneous symmetry breaking of this model for d ⩽ 4; see also Ref. 10 for more extensive analysis in d = 4.¹¹ 

Around the same time, Gross and Neveu¹² studied the 1/N expansion for a two-dimensional fermionic model, where N Dirac fermions interact via a quartic term (“four-fermion interactions”): namely, the Lagrangian is given by $∑i=1Nψ̄i(i∂/)ψi+λN(∑i=1Nψ̄iψi)2$⁠. This model had also been touched in Wilson’s paper.⁸ It was observed that the model has asymptotic freedom and symmetry breaking (for a discrete chiral symmetry) and produces a fermion mass.

Another class of model considered from the 1970s is the nonlinear σ-models in two dimensions. For instance, see Ref. 13, which considered models where Φ takes values in a sphere in $RN$⁠. It can be viewed as the limit of $RN$ valued σ-model with a potential $K(Φ⋅Φ−a2)2$ in a formal limit K → ∞. The idea was also extended to CP^N−1-type nonlinear σ-models^14,15 (where more general homogeneous space G/H valued models were also formulated).

Matrix-valued field theories are considered to be more difficult. In 1974, t’Hooft¹⁶ gave the solution of the two-dimensional Yang–Mills model with the SU(N) gauge group as N → ∞. After this seminal work by t’Hooft, the large N behavior of gauge theories was extensively studied by physicists; see, e.g., Refs. 17–20. See also Ref. 21 (and the references therein) for a broad class of matrix models.

It is certainly impossible to list all the vast literature; instead, we refer to Ref. 22 for an edited comprehensive collection of articles on large N as applied to a wide spectrum of problems in quantum field theory and statistical mechanics; see also the review articles²³ and Ref. 24 or Coleman’s excellent book (Ref. 25, Chap. 8) for summaries of the progress.

A standard approach in quantum field theory is perturbation theory, namely, calculating physical quantities of a given model by perturbative expansion in a coupling constant of the model. In large N problems, terms in perturbation theory can have different orders in 1/N, and the “1/N expansion” is a re-organization of the series in the parameter 1/N, with each term typically being a (formal) sum of infinitely many orders of the ordinary perturbation theory.

To explain this formalism a bit more, we first recall the usual perturbation theory for the Φ⁴ model with N = 1 and m = 0 (for simplicity) in (1.1). This could be found in a standard QFT textbook, e.g., Peskin and Schroeder.²⁶ Loosely speaking, the perturbative calculation of, for instance, a two-point correlation of Φ is given by a sum of Feynman graphs with two external legs and degree-4 internal vertices that represents the quartic interaction Φ⁴,

For the vector model (1.1) with N > 1, which was studied by Wilson,⁸ the perturbative calculation of the two-point correlation of Φ_i (with i fixed) is also given by the sum of Feynman graphs with two external legs and degree-4 internal vertices, but now, each vertex carrying two distinct summation variables and a factor 1/N that represents the interaction $1N∑j,kΦj2Φk2$ term in (1.1),

Heuristically, graph (a) is of order $1N∑j≈O(1)$ and graph (b) is of order $1N2∑j≈O(1N)$⁠. The philosophy of Ref. 8 is that graphs with “self-loops” such as (a) get canceled by Wick renormalization,²⁷ and all other graphs with internal vertices such as (b) are at least of order O(1/N) and, thus, vanish, so the theory would be asymptotically Gaussian free field. (We rigorously prove this in Secs. VI and VIII C in 2D and 3D.)

Since model (1.1) is invariant under any O(N) rotation of Φ, it is more natural to study the O(N)-invariant observables than the field itself, such as $∑j=1NΦj2$⁠. One needs to find a suitable scaling for such observables to see interesting fluctuations. In the “law of large numbers behavior” scaling $1N∑j=1NΦj2$⁠, one would expect a deterministic limit. Interesting fluctuation, then, arises in the next scale.

Following the above heuristic, we assume that the model is Wick renormalized, and so the perturbation series does not contain self-loops, and we consider $1N∑j=1N:Φj2:$⁠. Note that here $:Φj2:$ is the Wick square. The bubbled diagrams of the following form are all the O(1) contributions to its two-point correlation:

In fact, it is known in physics literature that these “bubble chains” are dominant (after Wick renormalization) contributions for this model: for instance, if we consider another observable $1N:∑i=1NΦi22:$⁠, the leading contribution to its one-point correlation would, then, be the “bubble loops,” i.e., the above type of diagrams with x, y identified.

Besides directly examining the perturbation theory, alternative (and more systematic) methodologies of analyzing such expansion were discovered in physics. We briefly recall some of them.

In Ref. 28, Schwinger–Dyson equations were used to study models of the type (1.1). These are essentially a type of integration by parts formulas. They yield a hierarchy of relations for correlation functions of different orders. This hierarchy does not close, so it is somewhat similar to perturbation theory, but one advantage of using Schwinger–Dyson equations is that one could again write a truncation remainder in terms of correlations of the field. This will be useful for rigorous studies (see Ref. 29), and we will come back to this point in Sec. VII.

Schwinger–Dyson equations are particularly powerful when we move from vector valued models of the type (1.1) to matrix valued models, such as Yang–Mills model. In the context of Yang–Mills model, the authors of Ref. 19 discovered relations between Wilson loop observables among different loops, now called Makeenko–Migdal equations or master loop equations. These are also broadly used in the study of random matrices, for instance, see the book³⁰ by Guionnet.

The problem of the perturbative argument discussed in Sec. I B (besides being non-rigorous) is that keeping track of all the diagrams and proving that there is no other diagrams appearing at the leading order are a combinatoric challenge; this challenge also exists when iteratively applying Schwinger–Dyson equations. A better way is to introduce a “dual” field or “auxiliary” field, which we briefly review now. A very nice exposition of this idea can be found in Coleman’s book (Ref. 25, Chap. 8) [see also the original paper⁹ or the review (Ref. 24, Sec. 2)]. The exact ways of implementing this idea vary among the literature, but roughly speaking, one will shift the action in (1.1) by a quadratic term involving a new field σ when calculating an expectation with respect to the measure dν^N,

This is because the integral over σ is simply a Gaussian integral, which yields a constant that can be absorbed into the partition function. We, then, obtain $exp(−S)$⁠, where

In particular, the quartic term is canceled. Still considering two-point correlation of Φ_i, for instance, one has the following type of diagrams:

Here, the dashed line represents the covariance of the “Gaussian measure” $exp(−∫N2λσ(x)2dx)$⁠, which is simply $λNδ$⁠. Each internal vertex is degree 3 because of the $σΦj2$ interaction term.

One could study this theory of (σ, Φ), but as explained in Ref. 25, since the above action is quadratic in Φ, one could “integrate out” (Φ₁, …, Φ_N) and obtain an “effective action” for σ. This Gaussian integral would, however, yield a non-local interaction because of the term |∇Φ_j|², namely, a determinant $det(−Δ+m−σ2)N/2$ or $det(1−12(−Δ+m)−1σ)N/2$ if we factor out the determinant of the massive Laplacian. Many physical properties can be deduced from this dual field theory for σ, for instance, one could exponentiate this determinant to obtain a measure of the form e^NF(σ), where F has the form tr log(⋯) so that large N limit becomes a question of critical point analysis for F. The 1/N expansion for correlations will also turn out to be a simple perturbative expansion in the dual picture.

A rigorous study of large N in mathematical physics was initiated by Kupiainen^29,31 around 1980. See also his review.³² The literature mostly related to the present article is Ref. 29, which studied the linear sigma model in continuum in d = 2 (with Wick renormalization) and proved that the 1/N expansion of the pressure (i.e., vacuum energy or log of partition per area) is asymptotic, and each order in this expansion can be described by sums of infinitely many Feynman diagrams of certain types. This paper uses both Dyson–Schwinger equations and the dual field method, together with constructive field theory methods, such as chessboard estimates³³ based on reflection positivity. Borel summability of the 1/N expansion of Schwinger functions for this model was discussed in Ref. 34.

Most of the other rigorous results are proved on lattices and, thus, ignore the small scale singularity (which is one of the key difficulties in our work). We also review some of these results here even if they are not necessarily directly related to our results. In Ref. 31, Kupiainen considered the N-component nonlinear sigma model (fields take values in a sphere in $RN$⁠) on the d-dimensional lattice with fixed lattice spacing. As mentioned above, this simplifies to “spherical model” as N → ∞, and for the spherical model, it is known that the critical temperature is 0 when d ⩽ 2 and is strictly positive when d > 2. The author of Ref. 31 proved that the large N expansion of correlation functions and free energy is asymptotic above the spherical model critical temperature (meaning all temperature if d ⩽ 2), and a mass gap was established for these temperatures when N is sufficiently large. The main idea therein was the dual field representation (for which the heuristic is discussed in Sec. I C), which yields nonlocal expectations, and this technical difficulty was solved using random walk expansion ideas from Brydges and Federbush.³⁵ Asymptoticity was later extended to Borel summability by Ref. 36.

In the 1990s, Kopper, Magnen, and Rivasseau³⁷ proved the mass gap for fermions in the Gross–Neveu model when N is sufficiently large, and the model has at least two pure phases. They again work with the dual field as illustrated above (which ends up with a determinant of Dirac operator) with suitable ultraviolet regularization and apply cluster and Mayer expansion techniques. Note that the Gross–Neveu model action has a (discrete) chiral symmetry, which prohibits a mass, but the mass is acquired in a symmetry breaking. The authors of Ref. 38 (on lattice Z²)³⁹ [on $R2$ but with a certain ultraviolet cutoff and a “soft constraint” K(Φ · Φ − 1)² in the potential], then, revisited the nonlinear σ model and showed a mass gap for N sufficiently large, and they both work with dual fields and determinants using cluster expansions.

A large N limit and expansion for the Yang–Mills model have drawn much attention, following t’Hooft’s 1974 work.¹⁶ In two space dimensions (continuum plane), Lévy⁴⁰ proved the convergence of Wilson loop observables to the master field in the large N limit, with orthogonal, unitary, and symplectic structure groups. The master field is a deterministic object, which was described conjecturally by Singer in 1995⁴¹ who pointed out its connections with free probability. Lévy’s work⁴⁰ also rigorously establishes integration by parts, i.e., Schwinger–Dyson equations both for finite N and in the large N limit for the expectations of Wilson loops, which in this context are called the Makeenko–Migdal equations. Note that the two-dimensional Yang–Mills model has special integrability structures, which allows one to make sense of the random Wilson loop observables in continuum—this was established by Lévy in Refs. 42 and 43, which Ref. 40 relies on. Around the same time, Anshelevich and Sengupta⁴⁴ also gave a construction of the master field in the large N limit via a different approach, which is based on the use of free white noise and of free stochastic calculus; these methods were developed earlier by Sengupta in Refs. 45 and 46. On the lattice of any dimension, we refer to Ref. 47 (respectively, Ref. 48) for computation of correlations of the Wilson loops in the large N limit (respectively, 1/N expansion), which relates to string theory.

We remark that while the Schwinger–Dyson equations for vector models (as we will discuss more in Sec. VII) are easy to derive, Makeenko–Migdal equations for Yang–Mills models are more difficult to derive (rigorously). The first rigorous version was established for the two-dimensional Yang–Mills model in continuum in Ref. 40, and alternative proofs and extensions were given in Refs. 49–52 on plane or surfaces. On lattice, for finite N and in any dimension, these were derived by Chatterjee (Refs. 47 and 53) and, then, in Ref. 54 using stochastic ordinary differential equations (ODEs) on Lie groups.

Stochastic quantization refers a procedure to turn Euclidean quantum field theories to singular SPDEs and, thus, bring stochastic analysis and PDE methods into the study. They were introduced by Parisi and Wu in Ref. 55. Given an action $S(Φ)$⁠, which is a functional of Φ, one considers a gradient flow of $S(Φ)$ perturbed by space-time white noise ξ,

Here, $δS(Φ)δΦ$ is the variational derivative of the functional $S(Φ)$⁠; for instance, when $S(Φ)=12∫(∇Φ)2dx$ is the Dirichlet form, $δS(Φ)δΦ=−ΔΦ$ and (2.1) boils down to the stochastic heat equation,

Note that Φ can also be multi-component fields, with ξ being the likewise multi-component. Another example is the “dynamical Φ⁴ equation,”

and it arises from the gradient of $∫12|∇Φ|2+λ4Φ4dx$⁠, which is model (1.1) with m = 0, N = 1.

The significance of these “stochastic quantization equations” (2.1) is that given an action $S(Φ)$⁠, the formal measure that defines an Euclidean quantum field theory

is formally an invariant measure⁵⁶ for Eq. (2.1). Here, DΦ is the formal Lebesgue measure and Z is a “normalization constant.” We emphasize that (2.3) are only formal measures because, among several other reasons, there is no “Lebesgue measure” DΦ on an infinite-dimensional space, and it is a priori not clear at all if the measure can be normalized. The task of constructive quantum field theory is to give precise meaning or constructions to these formal measures. This area of study has yielded numerous deep results, which we do not mean to give a full list of literature; we only refer to the book.^57,58

Regarding the well-posedness of (2.2), in two dimensions, two classical works are Ref. 59 where martingale solutions are constructed and Ref. 60 where local strong solutions are obtained, as well as the more recent approach to global well-posedness in Ref. 61. In three dimension, the local strong solution is obtained as an application of regularity structures⁶² and, then,⁶³ as an application of paracontrolled analysis.

Given the very recent progress of SPDE, a new approach to construct the measure of the form (2.3) is to construct the long-time solution to the stochastic PDE (2.1) and average the distribution of the solution over time. This approach has shown to be successful for the Φ⁴ model in d ⩽ 3 in a series of very recent studies, which starts with Ref. 64 on the 3D torus where a priori estimates were obtained to rule out the possibility of finite time blow-up. Then, in Refs. 65 and 66, the authors established a priori estimates for solutions on the full space $R3$⁠, yielding the construction of Φ⁴ quantum field theory on the whole $R3$ and verification of some key properties that this invariant measure must satisfy as desired by physicists, such as reflection positivity. See also Ref. 67. Similar uniform a priori estimates are obtained by Ref. 68 using the maximum principle.

Large N problems in the stochastic quantization formalism have also been discussed (on a heuristic level) in the physics literature, for instance, Refs. 69 and 70 and Ref. 71, Sec. 8. The review on large N by Moshe and Zinn-Justin (Ref. 24, Sec. 5.1) is close to our setting; it makes an “ansatz” that $1N∑j=1NΦj2$ in the equation would self-average in the large N limit to a constant; our results justified this ansatz and in the non-equilibrium setting generalize it.

As we will see, the methods in the proofs of our main theorems borrow some ingredients from mean field limit theory (MFT). To the best of our knowledge, the study of mean field problems originated from the work of McKean.⁷² Typically, a mean field problem is concerned with a system of N particles interacting with each other, which is often modeled by a system of stochastic ordinary differential equations, for instance, driven by independent Brownian motions. A prototype of such systems has the form $dXi=1N∑jf(Xi,Xj)dt+dBi$⁠; see, for instance, the classical reference by Sznitman [Ref. 73, Sec. I(1)], and in the N → ∞ limit, one could obtain decoupled SDEs each interacting with the law of itself: we will briefly review this calculation below.

In simple situations, the interaction f is assumed to be “nice,” for instance, globally Lipschitz;⁷² much of the literature aims to prove such limits under more general assumptions on the interaction; see Ref. 73 for a survey.⁷⁴ 

We note that mean field limits are studied under much broader frameworks or scopes of applications, such as the mean field limit in the context of rough paths (e.g., Refs. 75–77), mean field games (e.g., survey⁷⁸), and quantum dynamics (e.g., Ref. 79 and the references therein). We do not intend to have a comprehensive list, but rather refer to survey articles^80,81 and the book (Ref. 82, Chap. 8) in addition to Ref. 73.

The study of the mean field limit for SPDE systems also has precursors; see, for instance, the book (Ref. 83, Chap. 9) or Ref. 84. However, these results make strong assumptions on the interactions of the SPDE systems, such as linear growth and globally Lipschitz drift, and certainly do not cover the singular regime where renormalization is required as in our case.

Before we delve into the mean field limit for singular SPDE systems, we briefly review a very simple example regarding SDEs, which is taken from the work of Sznitman [Ref. 73, Sec. I(1)]. We will see that some of the simple ideas will be reflected when we move on to the much more technical SPDE setting. Consider

where B_i (i = 1, …, N) are independent Brownian motions. One can show that as N → ∞, the solution converges to that of (“McKean–Vlasov limit”)

where $ut=Law(X̄i(t))$⁠. Note that this is a decoupled system now. To prove the above result, we subtract the two systems and get

Here, the terms in the second line are put “by hand” and they obviously sum to zero. We, then, compare the first and third terms on the RHS and the fourth and fifth terms and the second and the last term; we show that the differences all vanish as N → ∞. The most important case is as follows:

Note that the nontriviality here is that the left-hand side “appears to be” 1/N² times two summations, which would be O(1) as N → ∞, but actually, we have a decay of 1/N. To show this, note that $f(X̄i,X̄j)−∫Rf(X̄i,y)ut(y)dy$ is mean-zero, so for j₁ ≠ j₂,

since $X̄j1$ and $X̄j2$ are independent.

To bound the other terms on the RHS of (3.2), we need some assumptions on f, for instance, we assume that f is globally Lipschitz. Then, we can easily bound the other terms on the RHS of (3.2) by $|Xj−X̄j|$⁠. Upon taking expectation and by symmetry of law, these are just $E|Xi−X̄i|$⁠. We could, then, apply Gronwall’s lemma to show that $Esupt∈[0,T]|Xi(t)−X̄i(t)|∼1/N$⁠.

The key to the above argument is the decay in (3.3), which relies on the mechanism demonstrated in (3.4) for L² products of independent and mean-zero random variables. One also needs some way to “move” other terms from the RHS to the LHS using Gronwall’s lemma or some other methods. We will see that the estimates in singular SPDEs will be much more complicated, but these basic mechanisms will still be present. For instance, our Theorem 5.1 can be viewed as a result of this flavor, in an SPDE setting, and in fact, the starting point of our proof is, indeed, close in spirit to Ref. 73, Sec. I(1) where one subtracts X_i from Y_i to cancel the noise and, then, bound a suitable norm of the difference.

Before discussing large N problems, we first briefly review the solution theory for the stochastic quantization of Φ⁴ on T², i.e., the case N = 1. For simplicity, here and in Secs.. V–VII, we set λ = 1. Equation (2.2) is only formal in d = 2 since Φ is expected to be distributional. To give a meaning to it, one should regularize and insert renormalization terms (also often called “counter-terms” in the context of quantum field theory⁸⁵). Namely, we take a sequence of mollified noises ξ_ɛ → ξ (for instance, convolving ξ with smooth functions χ_ɛ, which approximate Dirac distribution as ɛ → 0) and consider the mollified equation,

where C_ɛdiverges as ɛ → 0 at a suitable rate, with some initial condition Φ_ɛ(0) = ϕ.⁸⁶ If the sequence of constants C_ɛ is suitably chosen, the sequence of smooth solutions Φ_ɛ of (4.1) will converge to a nontrivial limit Φ = lim_ɛ→0Φ_ɛ.

The following argument is originally given by Da Prato and Debussche.⁶⁰ Write $L=∂t−Δ$⁠, and let Z_ɛ be the stationary solution to the mollified linear stochastic heat equation $LZε=ξε$⁠. The key observation is that the most singular part of Φ is Z, so if we decompose

we can expect the remainder Y_ɛ to converge in a space of better regularity. Substituting (4.2) into (4.1) yields

Here, we have defined

We will choose $Cε=E[Zε2(0,0)]$⁠, and the reason will be clear below. Equation (4.3) now looks more promising since noise ξ_ɛ (whose limit is rough) has dropped out. This manipulation has not solved the problem of multiplying distributions since the limit of Z_ɛ is still a distribution, which is why we need renormalization when define its powers. However, Z_ɛ is a rather concrete object since it is Gaussian distributed. This makes it possible to study the renormalization and convergence of $Zε2$ and $Zε3$ via probabilistic methods. For example, consider the expectation

where $χε(2)$ is the convolution of the mollifier χ_ɛ with itself and P is the heat kernel. Due to the singularity of the heat kernel P at the origin, this integral diverges like O(log ɛ) as ɛ → 0 in two spatial dimensions. The constant C_ɛ defined above simply subtracts this divergent mean from $Zε2$ and similarly for the other term $Zε3$⁠. This amounts to considering the renormalized Φ⁴ equation [Eq. (4.1)].

It can be shown that Z_ɛ, , and do converge to nontrivial limits in probability as space-time distributions of regularity C^−κ for any κ > 0.⁸⁷ In fact, thanks to Gaussianity of Z_ɛ, given a smooth test function f, one can explicitly compute any probabilistic moment of (f) and prove its convergence. By choosing f from a suitable set of wavelets or Fourier basis, one can apply a version of Kolmogorov’s theorem to prove that Z_ɛ, , and converge in C^−κ. We denote these limits by Z, , and .

Passing (4.3) to the limit, we get

We can prove the local well-posedness of this equation as a classical PDE by a standard fixed point argument. For this, we use a classical result in harmonic analysis under the name “Young’s theorem” (e.g., Ref. 62, Proposition 4.14 and Ref. 88, Theorems 2.47 and 2.52), which states that if $f∈Cα$⁠, $g∈Cβ$⁠, and α + β > 0, then $f⋅g∈Cmin(α,β)$⁠. Thus, if we assume that $Y∈Cβ$ for β ∈ (0, 2), then the worst term in the parentheses in (4.5) has regularity −κ. By the classical Schauder estimate, which states that the heat kernel improves regularity by 2, the fixed point map

is well defined on C^β since C^−κ+2 ⊂ C^β. With a bit of extra effort, one can show that over short time interval, the fixed point map is contractive and, thus, has a fixed point in C^β, which is the solution. To conclude, one has Φ = Z + Y, which is the local solution to the renormalized Φ⁴ equation on T²; the above renormalization is well known as “Wick renormalization,” and a commonly used notation is to write the equation as

The global (long-time) bound is first established by Mourrat and Weber.⁶¹ The start point is a PDE energy estimate—here, it will be an L^p estimate for Y. Namely, we multiply (4.5) by Y^p−1 and integrate over space,

where ⟨, ⟩ denotes the L² product. Note that the second term on the LHS arises from ΔY and the third term on the LHS arises from −Y³ in (4.5) (this negative sign is crucial since it yields a positive quantity on the LHS).

One, then, needs to control the terms on the RHS. We demonstrate the basic idea with the first term, namely, −⟨Y^p−1, 3ZY²⟩. Standard Besov space inequalities allow us to show that

Recall that Z ∈ C^−κ (which is essentially the same as the Besov space $B∞,∞−κ$⁠), so $‖Z‖B∞,∞−α$ is well bounded.

One, then, has the following interpolation inequality:

We would like to use Young’s inequality to write the product on the RHS into a sum and, then, “absorb” them by the second and third terms on the LHS of (4.8). To this end, some manipulation is needed,

Now, we see that the terms on the LHS of (4.8) show up here. By Young’s inequality $ab⩽aq1+bq2$ with a proper choice of q₁, q₂ such that $q1−1+q2−1=1$⁠, we can, then, bound $⟨Yp−1,ZY2⟩$ by (say) $110$ times the second and third terms on the LHS of (4.8) plus a finite constant. The other terms on the RHS of (4.8) are bounded in the same way. Putting these bounds together, one obtains global control on the L^p norm.

We emphasize that an energy identity of the form (4.8) followed by various bounds will be the start-point of our analysis below in the large N problems.

In this section, we consider the N-vector generalization of (4.7), which is the stochastic quantization of the linear σ-model in 2D,

Note that the solution Φ_i depends on N, but we omit this dependence in our notation. The solution theory in Refs. 60 and 61 explained in Sec. IV corresponds to the case N = 1. For a fixed N, these well-posedness results are easy to generalize to the present setting. Our primary goal in this section is to explain a priori bounds, which are stable with respect to N, and we will, then, send N to infinity to obtain the following mean field limit:

The precise meaning of this equation will be given in (5.3) and (5.10).

This equation is reminiscent of the mean field limit SDE (3.1) discussed in Sec. III in the sense that the equation is self-consistently defined, which involves the law of the solution (here, it depends on the second moment of Ψ_i). On the formal level, this equation arises naturally: assuming that the initial conditions ${ϕi}i=1N$ are exchangeable,⁸⁹ the components ${Φi}i=1N$ will have identical laws so that replacing the empirical average $1N∑j=1NΦj2$ by its mean and re-labeling Φ as Ψ lead us to (5.2).

In two space dimensions, (5.2) is a singular SPDE where the ill-defined non-linearity also requires a renormalization. As far as we know, this is the first example of singular “McKean–Vlasov SPDE,” which requires renormalization.

The above result is precisely given in the next theorem.

Note that (5.3) is still “formal,” i.e., its meaning needs interpretation; see (5.10). We could actually prove this convergence result under more general conditions for initial data (see Ref. 90, Sec. 4).

Following the trick of Da-Prato and Debussche discussed in Sec. IV, we consider the decomposition Φ_i = Z_i + Y_i, where Z_i is a solution to the linear SPDE,

and Y_i solves the remainder equation,

with some initial condition Y_i(0) = y_i. Here, we assume that the initial datum satisfies $E‖zi‖C−κp≲1$ for κ > 0 small enough and every p > 1, and $E‖yi‖L22≲1$⁠, uniformly in i, N. (5.8) has similar structure as (4.3): it has on the RHS terms cubic in Y and in Z and cross-terms between Y and Z.

The notations : Z_iZ_j:, $:Zj2:$⁠, and $:ZiZj2:$ denote the Wick products. It is easier to explain the Wick products for $Z̃i$⁠, which is the stationary solution to $LZ̃i=ξi$⁠. As in Sec. IV A and using the independence of $(Z̃i)i=1N$⁠, we should define

where $aε=E[Z̃i,ε2(0,0)]$ is a divergent constant independent of i (in fact, it is equal to C_ɛ in Sec. IV A). As in Sec. V, the above limits exist in C^−κ for κ > 0. The Wick products : Z_iZ_j:, $:Zj2:$⁠, and $:ZiZj2:$ are just suitable modifications by additional terms involving the initial condition z of Z (also subtracting the same a_ɛ; see Ref. 90, Sec. 2.1 for details).

A key step in Proof of Theorem 5.1 is a new uniformity in N bounds through suitable energy estimates on the remainder equation [Eq. (5.8)]. This is inspired by the approach of a priori bounds in Ref. 61 discussed in Sec. IV C, but subtleties arise as we track carefully the dependence of the bounds on N. Indeed, the approach to obtain global in time bounds for N = 1 is to exploit the damping effect from the term −Y³ in (4.5), which corresponds to the term $−Yj2Yi$⁠. However, the extra factor 1/N makes this effect weaker as N becomes large. In fact, the moral is that we cannot exploit the strong damping effect at the level of a fixed component Y_i; rather, we are forced to consider aggregate quantities, and ultimately, we focus on the empirical average of the L²-norm (squared) instead of the L^p-norm for arbitrary p, as in Sec. IV. We skip the details in this expository note and jump into the discussion on the mean-field equation and the proof of convergence.

We now discuss a bit more the solution theory for the mean-field SPDE (5.3). While the notion of solution we use is again via the Da-Prato/Debussche trick, the well-posedness theory requires more care than that for $Φ24$ discussed in Sec. IV since here the equation depends on its own law and we cannot proceed by pathwise arguments alone. Again, we understand (5.3) via the decomposition Ψ_i = Z_i + X_i with X_i satisfying

This follows by a bit algebra, noting that $Zi2$ and $Z̃i2$ in (5.3) “basically” cancel, except that the latter one is stationary, while the former one is not—this is the term (⋯) that we omit in the following discussion. What is more important is that here we actually introduced an independent copy (X_j, Z_j) of (X_i, Z_i) inside E, which turns out to be useful for both the local and global well-posedness of (5.3). Indeed, one point is that the term E[X_jZ_j]Z_i in (5.10) cannot be understood in a classical sense; however, we can view it as a conditional expectation E[X_jZ_jZ_i|Z_i] and use properties of the Wick product Z_iZ_j [such as in (5.9)] to give a meaning to this. In fact, after taking expectation, $E[Xj2]Xi$ in (5.10) also plays the role of the damping mechanism, which helps us to obtain uniform bounds on the mean-squared L²-norm of X_i.

Assume the corresponding decomposition for initial datum ψ_i = z_i + η_i. We denote $LT2L2=L2([0,T],L2)$⁠. We have the following.

Here, R contains terms that only depend on renormalized powers of Z and are well bounded. It has the following explicit form where i ≠ j and s > 0 is a sufficiently small parameter:

Recall in Sec. IV that the start point of the global estimate is an energy identity (4.8). The proof of Lemma 5.3 is again the same as in Sec. IV based on an energy identity,

where

We, then, proceed by bounding I^1,2,3 by $110E‖∇Xi‖L22+‖EXi2‖L22$ plus terms that only depend on Z (which eventually constitute R). We refer to Ref. 90, Sec. 3.2 for these details.

We now explain the basic ideas in the Proof of Theorem 5.1. Recall that Φ_i = Y_i + Z_i and Ψ_i = X_i + Z_i, and we define (in the same spirit as we did in Sec. III)

Now, we have the following energy identity:

where

In the definition of $I3N$⁠, to have a compact formula, we slightly abuse the notation for the contribution of the diagonal part i = j, where we understand Z_iZ_j to be $:Zi2:$ and $:Zj2:Zi$ to be $:Zi3:$⁠. The above could be viewed as a much more complicated version of (4.8).

Next, we bound $I1,2,3N$ one by one. This requires substantial calculations, but we illustrate a key point with $I3N$ in the following. Define

It is not very hard [see Ref. 90, (4.18)] to prove the following bound for s > 0 small:

with $Λ=(1−Δ)12$ and

A key point is that $Gj(k)$ are centered. In addition, components of X are independent, and so are components of Z. Therefore, although the terms in $R̄N$⁠, $R̄N′$ appear to be “order O(N),” we can actually show that upon taking expectation that they are bounded uniformly in N. Essentially, this boils down to the following fact: for mean-zero independent random variables U₁, …, U_N taking values in a Hilbert space H, we have (“gaining a factor 1/N”)

This is ultimately similar to the mechanism (3.4) in Sec. III for the ODE case. After considerably technical estimates, we can show that $supt∈[0,T]‖vi(t)‖L22→0$ in probability, as N → ∞, and Theorem 5.1 holds (see Ref. 90, Sec. 4 for details).

In this section, we are only interested in the stationary setting, and we simply write Z for the stationary solution to the linear equation $LZ=ξ$⁠.

We now study the invariant measure for the mean field Eq. (5.3), namely,

with E[Ψ² − Z²] = E[X²] + 2E[XZ] for X = Ψ − Z. The remainder X satisfies

An immediate observation is that Z is a stationary solution to (6.1). This follows since the unique solution to (6.2) starting from zero is identically zero [or one could simply set Ψ = Z in (6.1)]. Since the (massive) Gaussian free field $N(0,12(−Δ+m)−1)$ is invariant measure for the Z equation, we see that Gaussian free field is an invariant measure for (6.1).

Below, we use a semigroup $Pt*ν$ to denote the law of Ψ(t) with the initial condition distributed according to a measure ν. A natural question is whether (6.1) has a unique invariant measure. Since the equation involves its own law, it is unclear if the general ergodic theory can be applied directly in this setting. Fortunately, it has a strong damping property in the mean-square sense, which comes to our rescue and allows us to proceed directly by a priori estimates.