Polarons in spinless metals—a variational solution

We propose a simple variational solution for calculating one-particle spectral functions in lattice models of spinless metals with strong electron-phonon coupling. It is based on a generalization of the Momentum Average variational approximation for single polarons, combined with the assumption that the other fermions in the system are locked into an inert Fermi sea. We expect the method to be accurate for fermion addition spectral functions in metals with a small Fermi energy (nearly empty band), and for fermion removal spectral functions in metals with a large Fermi energy (nearly full band), provided that the characteristic phonon frequency is not too small. Both these regions are far from the region where the Migdal theorem holds, thus our results offer new insights into polaronic behavior in a largely unexplored part of the parameter space. Here, we show results for the Holstein coupling in one-dimension and present ways to gauge their accuracy, but ultimately this will need to be verified against numerical calculations. This variational method can be extended straightforwardly to higher dimensions and other forms of electron-phonon coupling.


Introduction
One of the important challenges in condensed matter physics is to fully understand the effects of electron-phonon coupling on the properties of quasiparticles in metals. To set the stage for the following discussion, we choose E F /Ω and λ to span the interesting parameter space, where E F is the Fermi energy of the bare carriers of concentration x, Ω is the characteristic energy of the phonon branch of interest, and λ is a dimensionless effective electron-phonon coupling defined so as to distinguish weak coupling λ ≪ 1 from strong coupling λ ⩾ 1. (Explicit expressions are provided below for our specific model). This choice means that other ingredients relevant in real materials-coupling to multiple phonons, presence of disorder, bare electron-electron interactions, etc-will be ignored hereafter, as a reasonable first step in attempting to solve this difficult problem.
Furthermore, we emphasize that we are only concerned here with metallic ground-states and their quasiparticles. We do not consider insulating ground-states due to a Peierls instability [1,2] and/or charge density waves [3], known to arise at special carrier concentrations x whose corresponding metallic Fermi surface would be nested.
Another complication we wish to ignore is an instability to a superconducting ground-state, whether of the well established Bardeen-Cooper-Schrieffer type at weak electron-phonon coupling [4], or a more exotic Bose-Einstein condensation of a bipolaron liquid that might occur at strong electron-phonon couplings [5,6]. To achieve this, we further restrict the following discussion to a liquid of spinless fermions (alternatively, one can think of fully spin-polarized charge carriers).
Even with all these multiple simplifications, the only regions of the parameter space where we have a good understanding of the properties of the resulting quasiparticles are (i) the single polaron limit corresponding to E F /Ω = 0, at any λ; and (ii) the Migdal region, where the Migdal theorem holds. The precise boundaries of the latter are not well established. Migdal argued that crossed diagrams can be ignored in the electronic self-energy for E F /Ω ≫ 1, because phase space considerations lead to a decrease of order Ω/E F of their contribution compared to that of non-crossed diagrams of the same order [7,8]. This argument becomes questionable when the number of relevant crossed diagrams vastly exceeds the number of uncrossed diagrams, as is the case at sufficiently strong couplings where one needs to sum up to high orders before convergence is achieved. This explains why E F /Ω ≫ 1 is supplemented by λ < 1, restricting the Migdal region to the top-left of the parameter space shown in figure 1. We 'tapered' this region all the way to E F /Ω → 0 for sufficiently small λ, because if the coupling is weak enough, then first order perturbation theory becomes sufficiently accurate and the self-energy reduces to the lowest order diagram. As a result, crossed diagrams can still be ignored (even though here they are not small compared to non-crossed diagrams of same order) because they only appear at higher orders. We note again that the precise shape of this region is not known. In fact, as we argue below, it likely depends on additional parameters.
Nevertheless, it is clear that so far, we do not have good knowledge regarding the quasiparticles of spinless metals with E F /Ω ⩽ 1 at medium and strong electron-phonon coupling λ ∼ 1. We begin to address this challenge here, by generalizing the momentum average (MA) approximation-a variational approximation known to work well at all λ in the single polaron limit [9][10][11]-to low concentrations x. The rough area where we believe its predictions to be accurate is colored cyan in figure 1, although the precise location of the top border depends on additional parameters.
The work presented here is limited to one-dimension (1D) and to the Holstein electron-phonon coupling, but can be straightforwardly extended to higher dimensions [9,10] and different electron-phonon couplings [12][13][14][15], as well as account for disorder [16,17], for multiple phonon modes [18] and for non-linear coupling [19], following the corresponding single polaron MA solutions.
The paper is organized as follows: section 2 briefly reviews the model, and its variational solution is discussed in section 3. Section 4 shows representative results, and section 5 discusses them and offers perspectives for further generalizations.

Model
As already mentioned, we study a 1D metal of spinless fermions interacting with an Einstein branch of phonons through a Holstein coupling: Here is the hopping of the spinless fermions, c i being the annihilation operator for a fermion at site i of the chain with N → ∞ sites. We set the lattice constant a = 1 and also take ℏ = 1. For a concentration x = N f /N, the Fermi wavevector is k F = πx and the Fermi energy is Ignoring the zero-point energy, the phonon energy is H ph = Ω i b † i b i , where b i is an annihilation operator for a phonon at site i. The Holstein coupling [20,21] has the well-known form ). In the single particle limit, it is customary to define the dimensionless effective coupling as: We will use this definition, and we will analyze its suitability to characterize the finite-x region below. At finite concentrations, it is convenient to rewrite the Hamiltonian in terms of dressed boson operators B i = b i + gx/Ω, so as to account for the average distortion due to the uniform fermion density. We note that this is an exact reformulation: The mean-field solution corresponds to settingH e-ph → 0, i.e. to ignoring the small fluctuations of both the fermion concentration, and of the dressed phonon displacements. The mean-field ground-state is then is the fermion Fermi sea, and| describes the uniform coherent distortion that is the vacuum for the dressed phonon operators, B i| 0⟩ = 0. The corresponding mean-field energy is: Thus, at the mean-field level, the fermions' dispersion is not affected by the electron-phonon coupling. The only effect of the latter is a lowering of the total energy due to a uniform lattice distortion driven by the uniform fermion density.

Variational solution for one-particle propagators at finite carrier concentrations
To move beyond the underwhelming mean-field solution, we use a variational method to calculate fermion addition and removal propagators, to try to better capture the effects of the electron-phonon coupling on the formation of quasiparticles and their dispersion. Because we work in a canonical formulation, we treat the two propagators separately as they characterize Hilbert spaces with different numbers of fermions N f ± 1.

Fermion addition
The fermion addition propagator is defined as: where |GS⟩ is the ground-state of the system with N f fermions, E GS is its energy,Ĝ A (z) = [z − H] −1 is the resolvent and z = ω + iη is the energy plus an artificial broadening η → 0. The first obvious obstacle in trying to calculate this propagator is that we do not know |GS⟩ and E GS . To make progress, we limit our calculation to small concentrations x = N f /N, where the probability of two fermions to meet is already low. Additionally, we note that for spinless fermions we expect the phonon-mediated effective interactions to be repulsive because Pauli's principle prevents two fermions from efficiently sharing the same distortion cloud. As a result, it is a reasonable starting point to assume that the lattice distortion dressing the additional fermion is entirely due to its coupling to phonons, while the other fermions act simply to block states in the Fermi sea and thus limit its scattering phase-space.
Thus, we impose the variational constraint that N f of the fermions are always in the Fermi sea state |FS⟩ and only allow the additional fermion to turn into a quasiparticle: This propagator can be calculated similarly to the MA solution for the single carrier limit, as we show next. The MA variational space is characterized by: (a) the size of the phonon cloud, i.e. the number of consecutive sites over which it is allowed to spread. There are no constraints on how many phonons are at any of these cloud sites, nor on the distance between the cloud and the fermion that is creating it. In the single carrier limit we have implemented calculations with up to a three-site cloud [13]. For the Holstein model, a one-site cloud is very accurate for Ω/t ⩾ 0.5 or so, and a three-site cloud further improves accuracy down to Ω/t ⩾ 0.1 or so [10,11,15]. Here we will implement both one-site and three-site cloud solutions, and only show results for Ω/t ⩾ 0.5. (b) the number n of free phonons which are not part of the polaron cloud defines the version MA (n) of the approximation [11]. Here we implement both n = 0 and n = 1 solutions. The former describes well the bottom of the polaron band but misplaces the location of the polaron+one-phonon continuum defining the top of the band, as those states imply the existence of at least one free phonon.
We now briefly review the main steps in the derivation of the one-site MA (0) solution, to show that it mirrors the single polaron result, except that the momentum-averaged free propagator is replaced with the one for the given fermion concentration x. The same holds true for all other MA versions, hence their corresponding expressions can be taken directly from previous work.
To generate an equation of motion (EOM) we use Dyson's identityĜ Next we use: where the Heaviside function Θ(x) insures that the additional fermion's momentum is above the Fermi sea, iH e-ph |mf ⟩ ≡ 0 if we do not allowH e-ph to excite particle-hole pairs from the Fermi sea. This leads to the first EOM: where we defined the generalized propagators: Their EOM are generated similarly. For a one-site MA (0) variational approximation, once the first phonon is created at a site i we do not allow phonons to be created elsewhere in the system. Combining this with always freezing N f of the fermions in the Fermi sea, leads to the recurrence equation: where we define the real-space free propagator: These 1D free propagators have analytic expressions. For completeness, these are presented in appendix A, where we also present additional arguments for the expected accuracy of the MA solution. These EOM are thus formally identical to those obtained for a single polaron, and reduce to them as . The solution is then formally similar: where A n (z) are continuous fractions defined by and A M (z) = 0 for a sufficiently large M ≫ 1. More details on the meaning of this variational approximation, as well as arguments for its expected quantitative accuracy at small x, are provided in appendix B.
The one-site MA (1) solution and the three-site MA (0) solutions proceed similarly but within the corresponding enlarged variational spaces. The corresponding expressions for the self-energy are like those for the single polaron counterparts, upon replacing g 0 (z) → g A 0 (z) [11,15].

Fermion removal
The fermion removal propagator is defined as: Its calculation proceeds similarly to that of the fermion addition propagator, with some care needed with details such as the changed Dyson identityĜ R (z) =Ĝ 0 R (z) −Ĝ 0 R (z)H e-phĜR (z) and the relevant free fermion propagator: whose expression is provided in appendix A. The one-site MA (0) solution is found to be: where these continued fractions are defined asÃ This approximation and its more accurate variants are expected to hold for x → 1, i.e. for a small concentration of holes away from a full band. This is because one can use particle-hole symmetry and map the removal propagator at a concentration x onto the addition propagator at the concentration 1 − x, upon also changing the energy ω → −ω. This is why the accuracy considerations mirror those for the fermion addition propagator (see appendix B), and why in the following we only show results for electron addition spectra at a small x.

Fermion addition spectral weights
We start by considering the spectral weights for fermion addition A A (k, ω) = − 1 π ImG A (k, ω + iη) for low fermion concentrations and at weak effective couplings. Representative results obtained with MA In both cases, we see a polaron band at low energies-this is a discrete eigenstate (here broadened by η) located below a much broader polaron + phonon continuum that starts at Ω above the bottom of the polaron band. At first sight, the only difference appears to be that at finite x, the spectral weight for fermion addition vanishes for k < k F = πx.
However, a more careful look shows a qualitative difference for the finite x results as compared to the single polaron results, related to the end point of the polaron band. It is well-known that in 1D, the single polaron band extends across the entire Brillouin zone for any λ > 0, i.e. there is a discrete polaron eigenstates below the continuum at any k. Its quasiparticle weight decreases very fast with increasing k, which is why this state is difficult to see on the scale of figure 2. To make it more apparent, in figure 3 we plot instead tanh[10A A (k, ω)] for λ = 0.1. The left panel shows the x = 0 result, and indeed reveals the entire polaron band, flattened just below the continuum but extending across the entire Brillouin zone. In contrast, the right panel shows that at finite x, the polaron state instead 'joins' the continuum, i.e. there is a k max so that for k > k max there is no discrete polaron state, the lowest feature being the continuum.
In other words, for low λ at any finite x we observe the famous 'kink' that is textbook knowledge in the Migdal limit. As already discussed, for these parameters the Migdal theorem does not hold, however at such small λ the lowest self-energy diagram is dominant so that Σ A (k, ω) ≈ g 2 g A 0 (z − Ω). This agrees with the Migdal result for g → 0. For x = 0, there is a divergence at ω = −2t in the density of states (DOS) of the free fermion band. The DOS is proportional to the imaginary part of g 0 (z), so the DOS divergence is responsible for a divergence in the real part of g 0 (z) at ω = −2t, through the Kramers-Krönig relations; the latter divergence, in turn, explains why there is a discrete low-energy pole at the polaron energy E P (k) satisfying E P (k) = ϵ k + Σ A (k, E P (k)) for any k. For any finite x, on the other hand, the bottom of the free fermion band (energies below ϵ(k F )) is blocked out and there is no singularity in the DOS at ϵ(k F ). As a result, the real part  of g A 0 (z) has only a finite peak at ω = ϵ(k F ), therefore the equation ω = ϵ k + Σ A (k, ω) only has a discrete solution at ω = E P (k) below the continuum for a restricted range k F ⩽ k ⩽ k max .
Spectral weights for larger values λ = 0.4, 0.7 and 1 are shown in figure 4 for the single polaron (top row) and at x = 0.1 (bottom row). Here the polaron band becomes fully visible in the single polaron case, but we can also observe it extending over the full Brillouin zone for the finite x. This is a consequence of the fact that k max increases with λ, as can be discerned from figure 2, and eventually reaches the Brillouin zone boundary.
For Ω/t = 1, x = 0.1, this happens somewhere between λ = 0.15 and λ = 0.4, i.e. still at weak couplings. Insofar as we use the existence of a full polaron band as a sign of 'polaronic physics' , it is clear that in 1D and for E F /Ω < 1, this sets in fast. This is to be contrasted with the Migdal limit E F /Ω ≫ 1 where 'polaronic physics' is expected to become relevant only for λ ∼ 1 [8], although that analysis is in 3D. Furthermore, it is worth noting that in 3D there is a critical value of λ even for the single Holstein polaron band to extend over the whole Brillouin zone (there is no singularity at the bottom of the 3D DOS, unlike in 1D) [22]. Thus, it is not clear how much of this difference is due to dimensionality, and how much to having very different E F /Ω ratios.
The other observations apparent from figure 4 are that at finite x, the polaron bandwidth is broader (even in the absence of the k < k F region) and the quasiparticle weight at the bottom of the band is larger than for the corresponding single polaron band (see different scales for λ = 1 contours). We can understand qualitatively the former observation as follows. As mentioned, the continuum starts at min k⩾kF E P (k) + Ω. For the single polaron, the minimum is at k = 0 whereas at finite x it is as k = k F , leaving more 'space' for a broader dispersion in the latter case. The latter observation is encouraging as it bodes well for the internal self-consistency of this variational scheme even at such strong couplings. Apart from the variational assumptions about the structure of the cloud dressing the additional fermion (whose validity improves with increasing λ), we were also forced to assume that the other fermions are undressed and locked in an inert Fermi sea. Of course, in reality each fermion will create its own cloud. As a result, even if the resulting object has a well defined momentum k < k F , with some probability the fermion could have any momentum k ′ with the cloud carrying k − k ′ . In other words, in the presence of these clouds, our assumption that the original fermions only block the states in the Fermi sea becomes suspect. The quasiparticle weight equals the probability that there is no cloud so the higher the former, the lower the latter, and the more valid is the inert Fermi sea approximation. We further discuss these issues in the last section.

Fermi speed and effective mass renormalization
A more quantitative way to characterize the quasiparticle (polaron) dispersion and its renormalization from the free fermion band is by considering the values of the polaron speed at the Fermi momentum v * F = dE P (k) dk k=kF (18) and corresponding effective mass: The ratio defining the effective mass is well defined for k F > 0. In the limit k F = 0 (for a single polaron), applying the L'Hôpital rule gives the usual link between the single polaron effective mass and the inverse of the 2nd derivative of the energy. Below, we calculate the ratios between these quantities and their bare fermion counterparts v F = 2t sin k F and m(k , it suffices to show only one of them. We choose to present the effective mass renormalization as it is a quantity more studied in the single polaron limit. The derivative dE P (k)/dk can be calculated directly, by extracting the energy E P (k) of the lowest pole over a small momentum range above k F , and fitting it to extract the slope at k F . Because the MA self-energies are momentum-independent at the levels implemented here (for the Holstein coupling, weak momentum dependence only appears in the self-energy for MA (2) and higher levels), a simpler approach is to use the fact that for a local self-energy, v * is the quasiparticle weight. It can be extracted by fitting A A (k F , ω) near ω ≈ E P (k F ) with a Lorentzian of width η, to find the corresponding weight Z(k F ). We implemented both approaches and found agreement to better that 0.1% accuracy. In the following we show results obtained with the second method, i.e. we use:  [24], also see [25]. In all cases t = Ω = 1.  and Ω/t = 1, as predicted by the three versions of MA implemented here. The largest variation between the three MA predictions is observed at λ = 1 in panel (a). This is the single polaron limit where numerical results are available at these parameters, for example in [25]. The red circle shows the bold diagrammatic Monte Carlo (BDMC) result [24] (error bars are smaller than the symbol size). Both the 1-site MA (1) and the 3-site MA (0) predictions are within less than 2% of this BDMC value. The agreement between the three MA results improves with increasing x, showing that these variational results are basically converged for these fermion concentrations.
These results allow us to revisit the issue of how to define a 'proper' effective coupling at finite x. In metals with a large E F , it is customary to define it as λ eff (k F ) = m * (k F )/m(k F ) − 1. In the low-λ limit we obtain a linear scaling with λ (as expected from perturbation theory), so if we followed that approach we would find λ eff (k F ) ∝ λ. In 3D and the Migdal limit, is the free electron DOS at energy ϵ [8]. This result is obtained after assuming that which is a good approximation if E F ≫ Ω but is clearly inappropriate in the limit E F < Ω analyzed here. The analytic relationship valid in our limit can likely be found with some good will, but its usefulness is doubtful considering that nonlinear corrections set in at rather low λ values, as shown in figure 5. This is why defining an effective coupling from the mass renormalization is less useful when stronger couplings are considered, and why instead we continue to use the single polaron effective coupling λ.
These results confirm that indeed the mass renormalization is reduced (i.e. the polaron bandwidth and its quasiparticle weight are increased) as x increases. This is clearly demonstrated in figure 6, where we plot the mass renormalization ratio versus x. The MA prediction is that especially at stronger couplings, the mass renormalization reduces significantly and quickly as the concentration of carriers increases. This is in qualitative agreement with 2D results reported recently in figure 1 of [26], which shows an increasing Z(k F ) with increasing E F /Ω at small values of the latter parameter (due to a small concentration x). Those results also show an eventual upturn as the Migdal limit is being reached, but we do not expect that regime to be well described by our variational approximation.
As a final comment, we note that we could also plot the results of figure 6 against E F /Ω = 2t [1 − cos(πx)]/Ω. Clearly, however, this 'composite' parameter is not sufficient to fully characterize the behavior in the small E F /Ω limit, as evident from the fact that in the limit x → 0 the results depend strongly on both λ and the adiabaticity ratio Ω/t. The situation might improve if a better definition for λ eff was identified, nevertheless it is likely that one needs three dimensionless parameters, e.g. λ, x, Ω/t to fully characterize the parameter space at small x, far from the Migdal limit.

Discussions and conclusion
We proposed a variational method for studying the effects of strong electron-phonon coupling on polaron formation in spinless metals. One key assumption is to lock the original N f fermions in an inert Fermi sea, and only allow the additional fermion/hole to create a polaron cloud. This assumption is only likely to be reasonable at small x (nearly empty band) for fermion addition, and large x (nearly full band) for fermion removal. Our specific model has particle-hole symmetry so we only showed results for the fermion addition problem at small x. We continue to focus our discussion on this case but note that it carries over to fermion removal spectral weights if x → 1 − x.
The second key assumption is that the dressing cloud accompanying the quasiparticle at small x is not all that different from that in the x = 0 limit, i.e. from that of a single polaron. For single polarons, MA has been shown to be quite accurate at all coupling strengths λ, if Ω/t is not too small. The particular MA versions implemented here allow for a rather compact phonon cloud, of up to three sites; for the 1D single Holstein polaron the results are accurate if Ω/t > 0.5 or so (in the strongly adiabatic limit, the extent of the phonon cloud increases and the variational space must be modified accordingly) [10,15].
The validity of the second assumption can be gauged by considering the variation of the results upon increasing the cloud size. For the parameters shown here, this leads to changes of few percent, giving us some confidence that the variational space for cloud configurations is quite sufficient (appendix B provides additional insight into this).
The validity of the first approximation is harder to gauge, especially at larger λ where the quasiparticle weights become quite small, i.e. where the fermions are dressed by quite large clouds. This may limit this variational approximation to only being quantitatively accurate at small λ. However, our results predict an increase in the spectral weight with increasing x at all λ, so perhaps the method is applicable up to strong(er) couplings. It will be necessary to compare these MA predictions against numerical results in order to fully understand its range of validity.
In this context, it is important to mention that BDMC results for the 2D version of this problem have been recently published [26], so there is hope that 1D results will become available soon either from BDMC and/or from 1D specific numerical methods such as density matrix renormalization group. Taken together, we hope that they will provide a good starting point for understanding the properties of these quasiparticles in this part of the parameter space that lies far from the Migdal limit.
If the comparison with numerical results finds that this finite-x MA extension performs poorly, then it may be possible to use an iterative approach to improve it, for example by replacing the bare fermions' dispersion with the calculated polaron dispersion E P (k).
If the comparison with numerical results will validate this finite-x MA extension in some region of the parameter space, it is then straightforward to extend MA to higher dimensions and to other couplings. In particular, MA is known to become more accurate in higher dimensions for single polarons [10]. The main challenge in going to higher dimensions is to calculate the real-space propagators of the bare fermions. Here we presented 1D results because we obtained analytical expressions for the 1D bare propagators, see appendix A. In higher D, numerical calculations will be necessary to get them; we plan to investigate this issue next. Studying other couplings in 1D is straightforward and can be done immediately, but is of dubious value before the method is validated.
Another obvious future direction is to consider spinful fermions. Here, there are two qualitatively different regimes to consider. The first is when the bare electron-electron repulsion is rather weak. The main new physics expected here is the formation of bipolarons, so it would make better sense to consider propagators for the addition of a singlet or triplet of two particles (again, in the presence of an inert Fermi sea). This is a lattice generalization of the famous Cooper pair calculation [27], except it would again be in the less studied limit of small E F /Ω and moderate-to-large λ. We note that MA has already been used successfully to study a single bipolaron [28] so this extension can be implemented straightforwardly as well. The other limit is when the bare electron-electron repulsion is large enough to prevent the binding of bipolarons in real-space. In this case, the presence of the opposite-spin carriers is irrelevant because they would neither block states for the additional carrier, nor have a tendency to share a cloud with it. As a result, in this case the approximation presented here can be used upon rescaling x → x/2.
Finally, it is worth commenting on the perspectives of using this approximation to study actual materials. An interesting possibility is related to the recent proposal of spin polaron formation in hole-doped spin-polarized moiré superlattices [29]. The method proposed here can be straigthforwardly extended to study finite hole concentrations, given that carrier-magnon coupling can be treated similarly to carrier-phonon coupling, as demonstrated in the context of hole-doped cuprate layers [30,31]. For spinful systems, the most intriguing option are STiO 3 and other quantum paraelectrics, which are known to have substantial electron-phonon coupling and to superconduct at surprisingly low carrier concentrations x [32][33][34][35].
To summarize, we believe that this work opens up a range of interesting possibilities for making progress in the study of this old problem.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors. larger n, but their amplitude is still small. The decay is even faster as the energy argument becomes z − mΩ, m > 1 in higher order diagrams.
The one-site MA (0) is equivalent to ignoring all contributions with n ̸ = 0 from all diagrams (not just 2nd order ones). As discussed in the main text, for Ω = t = 1 this is an excellent approximation at any λ, and the results in figure B1 suggest that the accuracy should remain similarly good for x < 0.15 or so.
It is important to note that the n = 0 contribution of the crossed and non-crossed diagrams are precisely equal. This is a consequence of the fact that the phonons are emitted at the same site and thus indistinguishable, so the order of re-absorbing them cannot possibly matter. Thus, the Migdal argument is clearly invalid in this small x limit, where we can ignore the n ̸ = 0 contributions.
For completeness, we mention that the increase of the variational space selectively adds more contributions to the self-energy diagrams. For example, the one-site MA (1) sets g A n (z − mΩ) → δ n,0 g A 0 (z − mΩ) only for m ⩾ 2, so Σ 2,nc (k, z) is calculated without further approximations. The three-site MA (1) keeps all contributions with n ⩽ 2 exactly and only ignores free propagators that move the free fermion by n ⩾ 3 sites, etc. Comparing their results thus gives us an estimate for how important are these ignored propagators, and therefore of the accuracy of the MA prediction.