The two-particle-irreducible vertex of the two-dimensional lattice $\phi^{4}$ model across the Ising transition

The single-particle propagator $G$ and self-energy $\Sigma$ are the standard objects of a correlated field theory [1, 16, 25], but they carry only part of the fluctuations: the response of the system to a change in $G$ is governed by the two-particle-irreducible (2PI) vertex $\Gamma=\delta\Sigma/\delta G$, the kernel of the Bethe-Salpeter equation [29]. It is the vertex, not the self-energy, that decides which collective channel becomes unstable and where the correlation length diverges. For the two-dimensional $\phi^{4}$ model, which is the paradigmatic realization of the Ising universality class [26], the transition is ferromagnetic, and one might expect the vertex to be dominated by the fully symmetric ($A_{1}$) channel at zero momentum transfer. Furthermore, the Green function and the self-energy depend only on a single momentum in a translation-invariant system and reside exclusively in $A_{1}$. However, as we will show, this expectation is incomplete: As fluctuations strengthen, the vertex develops structure in the other irreducible representations of the lattice point group, most notably a nematic $B_{1}$ component, followed by a diagonal nematic $B_{2}$ component.

In this work, we reconstruct the 2PI vertex $\Gamma(k,p;q)$, which follows by inversion from Monte Carlo measurements of the connected two-particle correlator. We resolve $\Gamma(k,p;q)$ in the irreducible representations of $C_{4v}$, track the leading Bethe-Salpeter eigenvalue as it approaches unity at the phase transition (also known as the Thouless point) and characterize the vertex both in momentum and in real space. To the best of our knowledge, such a channel- and transfer-resolved map of the 2PI vertex has not been reported in full for the two-dimensional Ising/$\phi^{4}$ problem. We will then extract the fully irreducible vertex of the parquet by stripping the reducible (crossed-channel ladder) contributions from the 2PI vertex, in the disordered phase up to the critical point. The fully irreducible vertex plays a central role in such theories as the parquet formalism [12, 11, 7], the $n-$PI theories [5, 9], the Schwinger-Dyson equations (SDE) [27] and the functional renormalization group (fRG) [24, 4]. This will allow us to establish D$\Gamma$A [28] as a very good approximation in the disordered phase, and to analyze the stability of the parquet equations in the critical region, where the physical solution turns into a repulsive fixed point.

We study the single-component $\phi^{4}$ theory on the square lattice with Euclidean action

$$S=-\beta\sum_{\langle ij\rangle}\phi_{i}\phi_{j}+m^{2}\sum_{i}\phi_{i}^{2}+\lambda\sum_{i}\phi_{i}^{4},$$

(1)

at $m^{2}=0$ and $\lambda=\tfrac{1}{2}$, sampled with the Brower–Tamayo cluster algorithm [8]. Monte Carlo does not give $\Gamma$ directly; what is accumulated is the connected two-particle correlator of the composite operator $\rho_{q}(k)=\phi_{k}^{*}\phi_{k+q}$,

$\displaystyle\chi(k,p;q)$

$\displaystyle=\langle\rho_{q}(k)\,\rho_{q}(p)^{*}\rangle-\langle\rho_{q}(k)\rangle\langle\rho_{q}(p)^{*}\rangle.$

(2)

The vertex is obtained by the Bethe–Salpeter equation (BSE) [29],

$$\Gamma=\chi_{0}^{-1}-\chi^{-1},$$

(3)

where $\chi$ is full susceptibility and the bare susceptibility $\chi_{0}$ is given by

$$\chi_{0}(k,p;q)=G(k)G(k+q)\,[\delta_{k,p}+\delta_{p,-k-q}],$$

(4)

which also carries the exchange term $\delta_{p,-k-q}$ required by the fact that the field is real ($\phi_{-k}=\phi_{k}^{*}$). The BSE can be used, for instance, to describe collective excitations (such as plasmons and magnons), bound states (such as excitons and Cooper pairs), and thermodynamic response (such as viscosity and thermal conductivity). The 2PI version of the BSE can be shown to be derivable in second order from a Luttinger-Ward functional, which ensures that macroscopic conservation laws such as charge and energy are conserved and that non-perturbative physics can be studied [2, 3]. We note in passing that the results presented here stand on their own and do not rely on a Luttinger-Ward functional, although they do have implications for theories built on it.

We monitor the instability using the symmetrized kernel given by

$$K=\chi_{0}^{1/2}\Gamma\chi_{0}^{1/2}=1-\chi_{0}^{1/2}\chi^{-1}\chi_{0}^{1/2},$$

where the leading eigenvalue reaches unity at the Thouless point. The function $\rho_{q}(k)$ satisfies $\rho_{0}(k)=\rho_{0}(-k)$ at $q=0$ due to the fact that the field $\phi$ is real. Similar relations hold for other values of $q$. Consequently, the susceptibility $\chi$ is rank-deficient by construction. Throughout this work, we focus on its physical (non-null) subspace, where the condition number remains moderate (see Appendix  F). At $q=0$, the operator $\rho_{0}(k)=|\phi_{k}|^{2}$ is even under inversion, which means that only the four inversion-even irreducible representations $A_{1},A_{2},B_{1},B_{2}$ can appear.

Figure 1 presents the leading eigenvalues for each symmetry sector. In the disordered phase, all channels are weak and comparable. As we approach the transition, the fully symmetric $A_{1}$ channel begins to separate, with its leading eigenvalue increasing toward the Thouless value of unity. This indicates a ferromagnetic instability as observed through the 2PI vertex framework. Initially, the nematic $B_{1}$ channel grows alongside $A_{1}$, reaching just over half its value at the largest system size. The finite-size trends for the two leading channels are as follows: for $A_{1}$, the values are logarithmically diverging as $0.75,0.85,0.88,0.89$ at $L=8,16,32,64$; for $B_{1}$, the values are $0.34,0.44,0.51$ at $L=8,16,32$. Note that for $L=64$ we could only resolve the leading $A_{1}$ eigenvalue because of increasingly demanding statistical requirements. The $B_{2}$ channel behaves similarly to $B_{1}$ but exhibits a smaller amplitude, while $A_{2}$ remains negligible. Therefore, the soft sector is not characterized by a single mode; instead, it is a cooperating multiplet dominated by $A_{1}$. Although only $A_{1}$ ultimately diverges, $B_{1}$ and $B_{2}$ contribute through their off-diagonal coupling with $A_{1}$, and should therefore not be omitted. The leading eigenvalues peak slightly beyond the thermodynamic transition point before receding in the ordered phase: the maximal $A_{1}$ eigenvalue is observed near $\beta\simeq 0.70$, while the clean second-moment correlation-length crossing indicates $\beta_{c}\simeq 0.685$ (App. G). Therefore, the vertex signal slightly overestimates $\beta_{c}$ for finite system sizes; we adopt $\beta_{c}\simeq 0.685$ throughout our analysis. It is worth mentioning that stripe and even nematic phases are predicted for the Ising model, but they require additional terms, such as a sufficiently strong antiferromagnetic coupling (which leads to frustration) along the diagonal and the presence of an external magnetic field [18].

The $E$ channel is exactly zero at $q=0$ due to parity. However, it takes on a finite value at the corner transfer $M=(\pi,\pi)$, which describes a short-range antiferromagnetic fluctuation. This fluctuation increases gradually through the critical window without nearing instability.

Figure 2 illustrates the behavior of the vertex as a function of momentum transfer. In the disordered phase, the leading eigenvalue is small and relatively flat across different values of $|q|$, indicating that the vertex is weak and lacks structure. As we approach criticality, a prominent maximum at $q=0$ emerges, but the decline away from this point is gradual. This results in a band of small-momentum modes becoming softer, which serves as the momentum-space signature of the diverging correlation length.

In the ordered phase, the scenario at the origin reverses: at $\beta=1.2$, the eigenvalue at $q=0$ drops to zero. The ferromagnetic weight that facilitated the transition has condensed into the uniform order parameter (or a collective coordinate for a finite system), which is one-particle-reducible and consequently absent from the 2PI vertex. However, the loss of weight is not limited to the single point $q=0$; rather, the entire small-$q$ region becomes depleted. At $\beta=1.2$, the eigenvalue in the ordered phase is reduced to approximately two-thirds of its critical value at the smallest momentum resolved on the $L=8$ lattice ($|q|=\pi/4$), while the curves converge near the zone boundary.

For the larger $L=16$ lattice, we are able to access smaller values of $|q|=\pi/8$, where the suppression intensifies, approaching one-half (data not shown). Thus, the short-wavelength structure of the vertex remains largely unchanged as we enter the ordered phase, whereas the long-wavelength component is suppressed. The gradual depletion of weight as $|q|\to 0$ is physically expected, as ordering suppresses long-wavelength fluctuations while preserving short-wavelength fluctuations.

The comparison with the disordered side is quite striking. For $\beta=1.2$, the absolute value of the relative detuning from $\beta_{c}$, calculated as $\left|\beta-\beta_{c}\right|/\beta_{c}\approx 0.71$, is significantly greater than that for $\beta=0.6$, where $\left|\beta-\beta_{c}\right|/\beta_{c}\approx 0.14$. However, the finite-$q$ eigenvalues for $\beta=1.2$ are consistently larger by about a factor of two at the smallest momenta, measuring $0.37$ compared to $0.17$ for $L=8$.

This indicates that the finite-$q$ fluctuation structure gaps out remarkably slowly on the ordered side. While the order parameter condenses sharply in the $q=0$ channel, the short-wavelength fluctuations are only slightly reduced and retain most of their critical strength well into the ordered phase. In contrast, these fluctuations have not yet developed on the disordered side.

This observation does not contradict the fact that the connected correlation length $\xi_{c}$ on the ordered side, as shown in Figure 7 and Appendix G, decreases much more rapidly as a function of detuning compared to the disordered side. For instance, $\xi_{c}\approx 2.3$ at the shoulder ($\beta=0.6$) drops to only $\xi_{c}\approx 0.4$ at $\beta=1.2$.

In the decomposition $K=\chi_{0}^{1/2}\Gamma\chi_{0}^{1/2}$, the vertex $\Gamma(q)$ remains essentially flat across all temperatures, while its overall magnitude increases significantly—by roughly a factor of fifty—from the shoulder into the ordered phase. Therefore, all the transfer dependence of $K$ is contained within the bubble $\chi_{0}(k;q)=G(k)G(k+q)$, whose profile is determined by the connected correlation length $\xi_{c}$: for long $\xi_{c}$, it is strongly $q$-dependent (enhanced at small transfer), while for short $\xi_{c}$, it is weakly $q$-dependent.

In Figure 3, we analyze the locality of the 2PI vertex in real space (refer to App. B for a discussion on the center-of-mass coordinate). The $x$-axis represents the relative coordinate conjugate to the momentum difference $k-p$, evaluated at zero transfer ($q=0$) and summed over the center-of-mass momentum $k+p$.

When away from criticality, the kernel shows a sharp on-site spike, which has effectively decayed by one lattice spacing. This indicates that the vertex is ultra-local, primarily influenced by its $r=0$ component, consistent with the nearly structureless momentum dependence depicted in Fig. 2. As we approach the transition, the weight begins to redistribute: the on-site component decreases, and the maximum shifts outward to $|r|=1$, with a faint tail extending to larger distances.

High-statistics data at $L=16$ indicate that the tail is best described by a power law of the form $\sim r^{-a}$. The effective exponent softens as it approaches the transition, with values of $a\simeq 2.4,2.0,1.8$ at $\beta=0.64,0.66,0.68$, respectively. This decay is significantly faster than the power-law decay of the propagator, which follows $r^{-\eta}$ with $\eta=1/4$ at the critical point.

The three cooperating channels illustrated in Fig. 1 are not an arbitrary selection. For those familiar with conformal field theory, these channels are the lattice analogs of the low-lying even (spin-neutral) operators of the Ising CFT: the energy density ($A_{1}$) and the two components of the stress-energy tensor ($B_{1}$, $B_{2}$). For readers acquainted with diagrammatic expansions, the channel $B_{1}$ is initially produced in the 2PI vertex by amputating the sunset diagram, while $B_{2}$ is first generated in higher order.

After determining the exact vertex through Monte Carlo simulations, we now turn our attention to its implications for diagrammatic theories. The vertex we have measured is irreducible in the particle-hole ($ph$) channel, but it does incorporate the resummation of the ladder diagrams in the two crossed channels: the particle-particle ($pp$) and the transverse particle-hole ($\overline{\text{ph}}$) channels. The fully irreducible vertex $\Lambda$ is the key object in the parquet formalism, which is irreducible in all channels [12, 11, 7]. Therefore, it is natural to consider what remains once these ladders are removed.

Since the field is real and the interaction involves a single on-site $\phi^{4}$ term, the connected four-point function maintains full crossing symmetry on the lattice (see App. C). The vertices from the three channels reduce to just one function, which is evaluated at three different momentum transfers: $q$ (ph), $k+p$ (pp), and $k-p$ ($\overline{\text{ph}}$). We have directly verified these identities on the connected correlator $\chi_{c}=\chi-\chi_{0}$ and found that they hold within the Monte Carlo data’s error bars.

We construct the vertex function $\Lambda$ from the same data using the parquet equation (our conventions for the parquet formalism are explained in App. D):

$$\Lambda=\Gamma_{\text{ph}}+\Gamma_{\text{pp}}+\Gamma_{\overline{\text{ph}}}-2F,$$

(5)

where each vertex is obtained through a Bethe-Salpeter inversion, given by $\Gamma_{r}=\chi_{0,r}^{-1}-\chi_{r}^{-1}$ for $r\in(\text{pp},\text{ph},\overline{\text{ph}})$. Here, $F=\chi_{0}^{-1}(\chi-\chi_{0})\chi_{0}^{-1}$ represents the fully amputated vertex. In the case of a finite Euclidean lattice system, the Luttinger-Ward functional is proven to be unique [23]. Consequently, the susceptibility $\chi$ is positive definite and possesses the properties of a covariance. This implies that both $\chi$ and each $\Gamma_{r}$ are invertible, ensuring that $\Lambda$ is finite by construction. It is worth noting that for fermionic systems, the Luttinger-Ward functional can become multivalued [22, 30, 20, 19], which adds complexity to a similar analysis.

The locality of $\Lambda$ is illustrated in the inset of Fig. 3. The relative-coordinate weight at one lattice spacing $(r_{1})$ is only 0.2% to 3% of its on-site value $(r_{0})$ across the entire window. We cannot discern the weight at larger distances due to noise. By stripping the Bethe-Salpeter ladders, we remove this tail as expected: $\Lambda$ is 20 to 110 times more local than the particle-hole vertex at $L=8$ and, within the resolution of our construction, behaves like a pure contact term. This observation is not a result of finite-size artifacts. For the $L=16$ data at $\beta=0.6$, we find that $\Lambda(r_{1})/\Lambda(r_{0})=0.0029\pm 0.0030$, indicating a pure contact term within $1\sigma$. At $\beta=0.64$, this ratio is $\Lambda(r_{1})/\Lambda(r_{0})=0.0150\pm 0.0056$, which remains a pure contact term within $3\sigma$. At $\beta=0.68$, we observe $\Lambda(r_{1})/\Lambda(r_{0})=0.0250\pm 0.0064$, again a pure contact term within $4\sigma$. While it is likely that the fully irreducible vertex develops a short-ranged but low-amplitude contribution for $r\neq 0$ at very large system sizes, establishing this with certainty is challenging, even with high-precision Monte Carlo data. All local mean values remain consistent with the error bars observed for $L=8$. Additionally, we note that removing the ladders also resolves the center-of-mass issue discussed in App. B. The center of mass (COM) structure of the 2PI vertex resides in the ladders, while $\Lambda$ carries less than 1% of it, increasing only slightly as we approach $\beta_{c}$ in the $A_{1}$ channel. Our findings place the fRG/parquet expectations and methods, such as D$\Gamma$A (see App E) [31, 28], on solid ground for $\phi^{4}$ models.

As the tail of $\Lambda$ is obscured by noise, we must explore alternative methods to determine its relevance. To achieve this, we replace $\Lambda$ with its effective contact value, defined as $\lambda_{\text{eff}}=\frac{1}{N^{2}}\sum_{Q,k,p}\Lambda(k,p;Q)$, and then we close the parquet loop. It is worth noting that the term $Q=0$ predominantly contributes to this summation; however, it exceeds $\lambda_{\text{eff}}$ by an unacceptable 2-3% when $\beta\geq 0.6$, resulting in significantly poorer outcomes that are not presented here. We will proceed with our analysis in several steps.

Firstly, as a preliminary test, we close the loop at the measured propagator. By inserting the single value $\lambda_{\rm eff}$ into the parquet approach at the Monte Carlo $G$ and propagating through Eq. (15) and the Dyson equation, we can reproduce the measured self-energy to better than a tenth of a percent throughout the convergent window. This confirms that the physical solution obtained from Monte Carlo methods is a fixed point of the parquet formalism, even in the critical region (though excluding the thermodynamic phase transition). From a physical perspective, the fully irreducible vertex functions as a screened local contact, which is only mildly renormalized from the bare interaction by a factor of approximately 0.85 (see Fig. 6) in the disordered phase, and this changes to an enhancement of approximately 1.15 at criticality. This indicates that the strong correlations contributing to the self-energy are primarily contained in the ladder diagrams that the parquet formalism resums, rather than in $\Lambda$ itself.

Secondly, our observation raises the natural question within the D$\Gamma$A framework: how does the measured contact compare with the local vertex that a practitioner would actually use, specifically the fully irreducible vertex of a dynamical mean-field impurity? In classical $\phi^{4}$ field theory, the impurity corresponds to a single site, which is solved self-consistently within a Gaussian bath through quadrature (this approach is equivalent to the DMFT approximation). The fully irreducible vertex, denoted as $\Lambda_{\rm DMFT}$, is derived from Eq. (5) specialized for one mode (refer to App. E for a detailed derivation). This leads to the D$\Gamma$A approximation, and its vertex is compared to $\lambda_{\rm eff}$ in Fig. 6. It is important to note that due to the simplicity of the theory, the fully irreducible vertex is frequency-independent, which is why we denote it as $\Lambda_{\rm DMFT}$.

Figure 4 compares the local self-energy to the one obtained through $\lambda_{\rm eff}$. This provides a fair test because local theories, such as DMFT and D$\Gamma$A, are expected to effectively capture local quantities. However, there is a notable difference compared to the previous paragraph, which analyzed the parquet equation for the measured propagator: here, we need to solve the full parquet equations (known as closure) without prior knowledge of the exact $G$. This can lead to stability issues (see below); in fact, the physical solution transitions from an attractive fixed point to a repelling one, notably throughout the entire critical region. For future reference, we refer to this procedure as a parquet–Schwinger-Dyson–Dyson sweep, or PSD; the earlier approach of parquet–Schwinger-Dyson sweeps at a fixed, exact $G$ will be referred to as PS.

Nevertheless, the figure clearly shows that away from criticality, $\Lambda_{\rm DMFT}$ and $\lambda_{\rm eff}$ agree within a few percent. Thus, we validate the assumption that the fully irreducible vertex is local and can be well approximated by a local impurity against the exact vertex. However, as we approach $\beta_{c}$, the measured contact is enhanced while the D$\Gamma$A vertex remains nearly constant (and incorrectly always screens). Next, we examine the non-local properties, which are not directly controlled by DMFT/D$\Gamma$A. As long as the (connected) correlation length $\xi_{c}$ remains less than one lattice spacing, DMFT/D$\Gamma$A remains virtually exact. However, as $\xi_{c}$ increases, the limitations of the local vertex become evident, as illustrated in Fig. 5. While the correlation function $G(r)$ obtained with the exact vertex aligns with Monte Carlo results at all distances, the D$\Gamma$A local vertex overshoots the long-range propagator by up to 10–20%, and this deviation increases with the distance $|r|$. This discrepancy highlights the non-local corrections that a momentum-independent $\Lambda$ cannot accommodate.

Thirdly, we return to the issue of convergence. The loss of convergence observed in the parquet self-consistency equations is not attributed to a physical singularity; rather, it is a result of damped fixed-point iterations acquiring a Jacobian eigenvalue of unit modulus or greater. This causes the physical fixed point to become repulsive [14].

This picture is supported by the Monte Carlo data: when inserting the exact solution at a fixed propagator, the spectral radius of the map increases smoothly from $0.90$ at $\beta=0.60$ to $0.94$ at $\beta=0.64$, crossing unity between $\beta=0.64$ and $0.66$. This places the linear stability threshold at approximately $\beta\simeq 0.65$. The self-consistent iteration, which couples this map to the Dyson update, loses convergence significantly earlier, at $\beta\simeq 0.62$, because the Dyson feedback is itself destabilizing (see below and App. H). Anderson-accelerated iterations [32, 15] and a $\beta$-annealed continuation do not help; they seek an attractor and stop at the same point where stability is lost. However, a Jacobian-free Newton-Krylov root-finder [21], which can converge to repulsive fixed points, does manage to extend the solution up to $\beta=0.65$—the onset of the single-mode repeller—but stalls just above this point at the mode proliferation phase ($\beta\geq 0.66$; see App. H). Therefore, the extent to which such methods can be applied to the two-dimensional lattice up to $\beta_{c}$ remains an open question. The scenario presented in Ref. [14] suggests a procedure to stabilize the equations in the non-convergent regime, particularly when dealing with a large eigenvalue. We reference App. H, where we analyze this in detail, but summarize the two key results here: First, at its onset, there is a single mode (specifically, the $A_{1}$ zero-transfer (order-parameter) direction), so the fixed-$G$ (PS) map remains an attractor up to its wall ($\beta\approx 0.64$–$0.65$). However, closing the Dyson loop (PSD) feeds that same soft mode back through the propagator, which makes it repulsive earlier (Eq. (20)): self-consistency is destabilizing, and the PSD map fails first. Second, this wall is not a fixed barrier below $\beta_{c}$; rather, its onset is a finite-size feature that recedes toward $\beta_{c}$ with increasing system size. Near criticality, the unstable set proliferates into many complex modes, and single-mode stabilization is no longer sufficient (App. H). In summary, solving the parquet equations remains a notoriously challenging and open problem, even with knowledge of the physical irreducible vertex. We emphasize again that this is not a limitation of the parquet equations themselves; they remain valid in the critical region.

We have reconstructed the 2PI vertex for lattice $\phi^{4}$ theory using Monte Carlo measurements of the connected two-particle correlator, which can serve as a benchmark object against which various theories such as fRG, parquet, D$\Gamma$A, and Schwinger–Dyson can be tested. We then analyzed the 2D $\phi^{4}$ Ising transition from the vertex perspective. Our main results can be summarized as follows:

First (see Figs. 1 and 2), the soft sector at the transition is multidimensional. The $A_{1}$ irreducible representation drives the instability, but the $B_{1}$ and $B_{2}$ representations cooperate and strengthen as the system size increases. Consequently, the minimal reduced description of the near-critical vertex spans these channels rather than relying solely on the $A_{1}$ mode. The contribution from $B_{2}$ is tied to $B_{1}$ because of the emergent rotational invariance at the critical point.

Second, the two vertices differ in range (see Figs. 3 and 6). The fully irreducible vertex $\Lambda$ is spatially compact—a near-contact object with approximately 40 times more amplitude at $r=0$ than at $r=1$ even at criticality—while the particle-hole irreducible (2PI) vertex develops a power-law tail at criticality and shows non-monotonic behavior at small $r$ in the relative coordinate, along with a complex structure in the center-of-mass coordinate. In the disordered phase, both vertices are localized within $r\approx 1$; here, $\Lambda$ screens the bare contact by approximately 10%, whereas toward criticality, the effective contact is instead enhanced (as shown in Fig. 6).

Third, the D$\Gamma$A approximation for $\phi^{4}$ theory accurately reproduces the fully irreducible vertex and the local self-energy (see Fig. 4), but begins to fail when the correlation length exceeds one lattice spacing. This failure is first evident in the closure of the correlation function (see Fig. 5) and can also be observed, albeit to a lesser extent, in local quantities (see Fig. 4).

Fourth, in the critical region, the physical solution of the parquet-SDE-Dyson (PSD) equations turns into a repulsive fixed point, and we clarify the structure of this instability (see App H). This instability begins as a single order-parameter ($A_{1}$) mode, which contributes an eigenvalue to the Jacobian with a magnitude greater than one. The fixed-propagator map (PS) remains an attractor up to a certain boundary that self-consistency (PSD) pushes to a lower wall. This boundary is a finite-size feature that approaches $\beta_{c}$ as the system size increases. Following Ref. [14], we stabilize that single mode, enhancing convergence; however, nearer to the transition, the unstable set proliferates into many largely complex modes, and a rank-one flip is no longer sufficient. This is a stability, not an existence, issue—it is not a property of the PSD equations themselves, nor is it related to the uniqueness of the Luttinger-Ward functional (which is guaranteed to be unique in our case); only the solver’s Jacobian becomes unstable.

Collectively, these results establish the directly measured vertex as a controlled, statistically valid benchmark for two-particle-based theories—such as parquet, D$\Gamma$A, fRG, and Schwinger–Dyson—on the disordered side of lattice field theory, and they precisely locate where these closures begin to fail.

In principle, this approach could be extended to any $N$-component $\phi^{4}$ theory on any $D$-dimensional finite lattice in the disordered phase, provided an efficient cluster algorithm exists for the two- and four-point correlators. However, we caution that this is not a scalable, general-purpose method: measuring the fully connected four-point function accurately enough to extract the Bethe–Salpeter ladders is costly in both statistics and memory, and these costs grow rapidly with the number of components, dimensionality, and correlation length. Therefore, we expect this method to remain practical only for a few comparatively simple models. The most natural candidates include the 3D Ising model, the $O(N)$ vector models (with the $N=2$ quasi-long-range-ordered critical phase being a particularly challenging test), and Heisenberg $O(3)$ spin models, whose vertex already couples nearest neighbors. In each case, the disordered side is accessible, but solving the closure relation (PSD) near criticality remains a formidable and likely model-dependent problem; we see no reason to expect it to be easier in that region than it was in our case. Future work will address the ordered phase and extensions to bosonic and fermionic many-body systems, where the vertex is significantly more complex and further conceptual complications arise, including the multivalued nature of the Luttinger-Ward functional in the fermionic case.