Strong-coupling diagrammatic Monte Carlo technique for correlated fermions and frustrated spins

We describe a controllable and unbiased strong-coupling diagrammatic Monte Carlo technique that is applicable to a wide range of fermionic systems and spin models. Unlike previous strong coupling methods that generally rely on the Grassmannian Hubbard-Stratonovich transformation, our construction is based on Wick's theorem and a recursive procedure to group contractions into effective connected vertices that are non-perturbative in all local physics and can be calculated exactly. The resulting expansion is described by simple diagrammatic rules that make it suitable for systematic treatment via stochastic sampling. Benchmarks against numerical linked cluster expansion display excellent agreement.

Strongly correlated electrons and frustrated spin models are among the most challenging problems in condensed matter theory due to a combination of the sign problem, lack of a natural small parameter and the computational complexity of series expansions. Coincidentally, this topic is at the same time absolutely central to understanding the electronic structure of solids, and a wide range of numerical techniques have accordingly been devised to overcome these obstacles. A wellknown example of this is DMFT [1], with extensions based on diagram techniques, [2,3], cluster generalizations [4][5][6], and related methods [7,8]. Other examples include DMRG [9], wave function methods [10,11], and auxiliary-field quantum Monte Carlo [12][13][14][15].
Several of the aforementioned techniques are capable of producing states that seem highly relevant for cuprate superconductivity, including anti-ferromagnetism, stripes, pseudogap physics, and d-wave superconductivity. Nevertheless, it has been known for some time that notable discrepancies may appear both when comparing different techniques and when altering details of the implementation (such as discretization) of a given method [16]. This sensitivity that correlatedfermion models display to numerical protocol appears to have a physical origin, and be rooted in competition between different states situated very closely in terms of free energy [16,17].
More recently, systematic comparison of leading numerical protocols, applied to the Hubbard model, has demonstrated some of the very significant progress that has eventually been made in this field [18]. In a substantial part of the parameter space, key observables can now be obtained from completely different techniques with a high degree of consensus. The region that remains the most problematic corresponds to small but finite doping, and intermediate to large onsite repulsion. Incidentally, this is also a scenario that is highly relevant for cuprate superconductors.
The sensitivity that correlated systems display to perturbations suggests that a reliable solution to this problem requires numerical methods that can provide extremely accurate results in the strongly correlated regime, and this has proven to be a significant challenge. The unbiased methods-that are free of systematic errors beyond some form of truncation-are typically based on a series expansion of some form. However, for strongly correlated systems, conventional perturbation theory is not viable, as the interaction is far too strong. Instead, an alternative expansion parameter has to be found.
The perhaps most well-known formalism aimed at the correlated regime is strong-coupling expansion, where the nonlocal processes are treated as a perturbation while the unperturbed system corresponds to the atomic limit [19,20]. Thus far however, most of these works include only modest expansion orders, and they are primarily conducted at half-filling, or for actual spin systems [21], while the case of non-zero doping is technically far more challenging [22].
Numerical linked cluster expansion (NLCE) [23] has been applied successfully to both spin models [24] and itinerant fermionic theories like the t-J [25] and Hubbard models. For the latter, results exist at infinite onsite repulsion [26], and for finite interactions up to U/t = 16 [27], which is far into the strongly correlated regime. This method is based on exact diagonalization of small clusters, and allows convergence to macroscopic results to be observed with increasing cluster size.
More recently, the extremely correlated fermi liquid theory was developed specifically for Gutzwiller-projected models [28]. This framework allows a form of diagram technique to be employed on restricted Hilbert spaces [29], and currently published benchmarks, while limited to low expansion orders, appear encouraging [30].
Finally, a number of adaptions of diagrammatic Monte Carlo methods [31] have been made to address the strongly correlated regime. Universal fermionization [32] has opened a new analytical path where restrictions on the Hilbert space are encoded via non-Hermitian terms in the Hamiltonian, thus allowing Gutzwiller-projected systems to be treated within the framework of Wick's theorem. Via second fermionization, doubly occupied sites can then be reintroduced in the form of hardcore bosons, which are subsequently fermionized, thus allowing generic correlated systems to be addressed within this framework [32]. This technique suffered from poor convergence properties at its inception except at large doping, but this problem has since been overcome through spin-charge transformation, which essentially results in a representation involving fermionic carriers that propagate on a spin background [33]. Diagrammatically, these models can be treated as a spin system using Popov-Fedotov fermionization [34], where the spins are mapped onto fermions with an imaginary chemical potential. Results from spin-charge transformation arXiv:1909.00816v2 [cond-mat.str-el] 20 Apr 2021 diagrammatic Monte Carlo (SCT-DMC) indicate that the expansion converges quite rapidly, but that the complexity of the resulting theory limits the expansion order, making it hard to reach low temperatures [35].
Currently, several new analytical techniques are being evaluated to overcome the inherent problem of a large expansion parameter: By extracting the analytical structure of the selfenergy from weak-coupling data, it becomes possible to reconstruct it in the non-perturbative regime [36]. Homotopic action operates on the principle of altering the starting point of the expansion, as well as the expansion parameter, such that the point of interest falls within the convergence radius [37]. Taking advantage of determinant sampling techniques, [38,39], these methods give access to fairly strong interactions, up to U/t = 7 in the Hubbard model. While impressive, this is still short of U/t ∼ 12, which is relevant for the high-temperature superconductors.
Thus, despite the significant recent advances, much of the parameter space remains challenging to unbiased techniques, and it is the case of strong correlation that remains the most elusive. In this work, we will discuss how diagrammatic techniques can be adapted to the strongly correlated regime by an alternative series expansion that is based on a non-perturbative treatment of all the interactions, and expansion only in the nonlocal part of the Hamiltonian. This series expansion is computationally economical, possesses simple diagrammatic rules, and is free of any large expansion parameter, even for arbitrarily strong interactions.

Model and diagrammatic expansion
As a starting point for the derivation of the new diagrammatic description, let us assume a Hamiltonian of a form that encapsulates the processes generally found in models of twocomponent lattice fermions and fermionized spin systems: Here,μ is assumed to be local and bilinear, i.e., a chemical potential. The termÛ is a contact interaction that is local and non-bilinear. The operatorĴ describes a nonlocal interaction that is mediated by a Boson, like super-exchange or the nonlocal part of a Coulomb interaction term. Finally,t is nonlocal and fermionic, generally corresponding to hopping. In principle, we can treat the model (1) through expansion in H 1 , and due to bi-linearity of H 0 , the contractions can be evaluated using standard Matsubara formalism based on Wick's theorem [40]. Furthermore, we note that the unperturbed theory is also local, so that all contractions of operators that are separated in space vanish a priori, i.e.
At this stage, we aim to exploit the combination of bi-linearity and locality of the unperturbed theory. Thus, we first recall that Wick's theorem allows us to obtain an answer from the series of connected diagrams by cancellation of disconnected contributions and the partition function [40]. Secondly, we notice, in accordance with (3), that all calculations are carried out in the atomic limit, where the full expectation value of an operator is generally trivial to obtain, and does not even require the evaluation of diagrams. In particular, this allows the evaluation of certain classes of terms up to infinite order, thus paving the way for non-perturbative treatment of contact interactions, for example. When using these properties in conjunction, we do however face a fundamental problem in that the full contraction of a set of operators contains a mix of connected and disconnected topologies, which runs very much contrary to the concept of diagrammatic expansions. This complication is further bolstered because connectivity of a set of contractions is a nonlocal property. The principal solution to this problem is to divide the set of contractions on a lattice site i into groups according to their connectivity, which effectively gives rise to a set of connected vertices that form the basis for an alternative diagrammatic technique.

Strong-coupling vertices
Let us start by dividing the second part of the Hamiltonian (1) into local and inter-site terms according to where i refers to lattice sites. Then, we proceed to introduce the following shorthand notation for the normalized timeordered integration which appears in expansions of the form (2): with the generalization We can now write the expansion in H 1 as where the subscript of U refers to lattice site and the string of operators H n 1 are assumed to depend on τ 1 ...τ n . Next, we introduceŌ i to denote the set of operators on the site i that are attributable to nonlocal terms in the Hamiltonian (i.e., H I ) or the measured operatorô. Expressing the expansion (2) in this new language, we find Here, n is the expansion order of the inter-site terms H I , whilē x denotes their spatial degrees of freedom, which are summed over accordingly. The subscript c implies connected topologies.
Since the bare Greens function is local (3), it follows that all contractions in (9) may be carried out locally also. However, the problem that remains is that we are interested in topologies that are connected, and this very property is nonlocal. To overcome this difficulty, we cannot simply compute the local trace; instead, we have to sort all local contractions according to their connectivity. To do so, we begin by breaking out the local terms on the site i that are not connected to any of the operatorsŌ i : Here, the subscript µ,e denotes the subset of contractions such that all diagrammatic elements are connected to at least one external line (i.e., an operator that is attributable to a nonlocal process). Discarding the disconnected topologies of (10) we may write (9) as where subscript c once again implies the subset of these topologies that are connected. However, this is a nonlocal property that depends both on the contractions on a site and the inter-site processes that connect different sites. To make further progress, we have to sort the local contractions according to the connectivity of the set of operatorsŌ i . Thus, we define the connected set of contractions denoted by ... µ,c as the subset of ... µ,e for which all elements ofŌ i are connected by local contractions. As an example, we may consider the case whenŌ has only two elements: Thus, we have sorted the contractions into two parts: Those where O 1 , O 2 are connected, which defines the Connected vertex, and those where they are disconnected. In principle, we can generalize this procedure to the case of an arbitrary number of elements ofŌ by constructing a recursion that is reminiscent of Determinant diagrammatic techniques [38,39]: For a given set of operatorsŌ, we take as our starting point a set of local contractions of the form The set of connected topologies may be obtained by subtracting those that are disconnected. To list these, we begin by sorting them according to which of the operatorsŌ are connected to O 1 , and denote this set by A. The set of contractions for which the operators in A are connected to each other, but not to the remaining operators (Ô \ A) is by definition given by where ξŌ ,A is a fermionic sign given by where c is the number of fermionic commutations associated with the reordering In the next stage, we have to sum over all possible choices of A. Here, we have two restrictions: Firstly, since A is the set of operators connected to O 1 it follows that A must contain O 1 . Secondly, since we are only interested in the set of disconnected topologies, it follows thatŌ \ A must be a nonempty set-otherwise, all operators are connected. Thus, A is a proper subset ofŌ. Summing over all choices of A, we obtain Using (17) we can then construct a recursive relation for the connected vertex on the form such that it can be directly computed from terms of the form (13). We then recognise that (10) can be written Here it should be pointed out that for the Hubbard model, the summation (20) has a finite convergence radius with respect to U [41]. It is also known that this model possesses singularities in the irreducible vertex functions at strong interactions [42]. However, the divergence of (20) can be removed by applying second fermionization [32] to the Hubbard model. A proof of this is given in appendix A.
Noting that the summation (20) may be rendered convergent and equated to an expectation value taken in the atomic limit, it can be computed exactly, and so the construction of the connected vertex (18) is an exactly solvable problem. This allows the contact interactions to be treated completely nonperturbatively.
We also note that any set of the form (13) may be decomposed into a sum of sets of connected vertices: Inserting this into (11) we find that all topologies in the expansion may be expressed in terms of connected vertices. Likewise, any set of connected vertices may be expressed in terms of contractions of the form (11). Therefore, we can expand directly in connected topologies of connected vertices: where η is the normalisation. In the next section, we will establish rules for the expansion (22).

Diagrammatic rules
To obtain connected topologies, the connected vertices defined by the recursion (18) have to be connected via external lines that originate in the nonlocal operators, i.e.t,Ĵ. However, a complication that remains is determining the sign of a contribution, which generally requires the construction of diagrammatic rules that govern the expansion. Since this derivation is essentially based on Wick's theorem, we first have to make contact with Feynman type diagrammatics in order to derive the corresponding principles for the strong-coupling expansion.
In standard literature, the expansion is typically conducted in conventional two-body interactions, and the overall sign of a diagram is generally expressed in terms of the number of fermionic loops [40]. Proceeding to more general models that for example include projected hopping, the resulting Feynman rules must generally be obtained from Wick's theorem. A convenient way of doing this is to introduce a reference contraction: Specifically, we understand that we can write a fermionic theory on a form where creation and annihilation operators form natural pairs, whose contraction corresponds to an infinitesimally backward propagating fermion. Thus, for an expression of the form U n i↑ n i↓ Jn iσ n iσ tn j↑ c † j↓ c k↓ n k↑ ... n iσ = c † iσ c iσ , (23) every operator is contracted with its natural partner to form a diagrammatic element as shown in Fig. (1, a), for which the fermionic sign is positive. Swapping the operators being contracted (Fig 1, b) gives rise to a fermionic sign, and so all diagram topologies can be characterized by whether they are related to the reference by an even or an odd number of such swaps.
Adapting this idea to the strong-coupling expansion, the first natural stage is to define a reference contraction for the To confirm that this diagram indeed carries a positive fermionic sign, we simply note that from the underlying operators, we can form the Feynman reference contractions of the type (1, a) without commuting any of them. Thus, if we, for example, assume that the external lines in Fig. (2, a) are fermionic, then we obtain an operator product of the form Summing up all contractions of (24) in accordance with the underlying Feynman series (including disconnected topologies), this is equivalent to the expectation value of the operators, corresponding to a positive fermionic sign.
As illustrated in Fig. 2, we require two basic updates to generate arbitrary diagrams from a set of reference contractions, namely swapping the connectivity of two external lines, and commuting operators within a vertex. For fermionic lines or operators, particle statistics suggest that these operations are odd, and this is indeed what transpires from the underlying Feynman type diagrammatics: Firstly, swapping the connectivity of two fermionic lines is equivalent to changing the connectivity of an odd number of fermionic propagators, which according to Wick's theorem, is an odd operation. As an example, we may consider the operation (a → b) in Fig. 1, which could alternatively be realized by a swap of the fermionic operators or the t−lines. Secondly, the process (c→ d) in Fig.  2 can be achieved with either a swap or a commute, implying that these operations have the same parity. By the same logic, operations on bosonic external lines do not give rise to a sign. Figure 2. Diagrammatic rules. The reference contraction of a connected vertex (a) is obtained by taking the external lines to be non-intersecting, have no time-overlap, and be forward-propagating in the case of fermions. In our notation, the horizontal line corresponds to imaginary time, and so forward propagation implies that the external line is traveling from left to right. To generate arbitrary diagram topologies from reference vertices, we require two basic updates: Swapping the connectivity of two external lines, we can go from (a) to (b), whilst commuting the operator order takes us from (a) to (c). The diagrams (c) and (d) are related by alternatively a swap or a commute operation. Swapping external lines also allows us to connect different vertices, such as going from (e) to (f).

Analytic structure of the connected vertices
Whilst the recursion (18) provides a principal definition of the diagrammatic elements of the expansion, computing, and also storing these objects in memory is only possible given an efficient representation. In particular, for a vertex with N external lines, the naive description yields N −1 imaginary-time differences and equally many dimensions of the mathematical object to be constructed and stored, so that the task quickly becomes intractable.
The solution to this problem can be found by noting that in the recursion (18), the vertex is expressed in expectation values of the form (20), which are taken with respect to the entire local part of the Hamiltonian, i.e. H L =μ +Û . If we express the nonlocal part H I in an operator basis which possesses a trivial time-evolution with respect to H L , then the time-dependence of the entire vertex becomes equally simple. The specific representation which allows this to be achieved, and which notably also forms the starting point for derivation of the t-J model [43], takes the form whereσ = −σ, while the corresponding time-dependence with respect to H L is given by Given a set of operators of the form (25) that are evaluated with respect to H L , we can use (26) to divide it into a scalar part, which consists of an analytic function and an operator part which only depends on the order of the terms, according to Since the recursion (18) implies that the connected vertex can be expressed in terms of expectation values of the form (20), it follows that we can break out the scalar part from this expression, and thus obtain an object of the form where f is an analytic function, while is a constant which only depends on the order of the operators, and correspondingly may be stored as a single floating point. Furthermore, we may note that using the basis (25) and exploiting the property (28), (20) essentially corresponds to the expectation value of a projection operator, which can be calculated exactly, and so the connected vertex is naturally obtained to machine precision. Finally, let us comment on the prelusive question about the feasibility of storing the vertices in lookup tables: For the Hubbard model, the basis (25) gives a total of 8 operators, implying that the number of vertices scales as 8 N where N is the number of external lines or legs. At N = 10, this gives ∼ 10 9 vertices, which translates to approximately 8 GB at double precision. For the Heisenberg model, which can be described by only 4 operators, we can afford to store all vertices up to N = 15 with the same resources. Generally, Gutzwillerprojected systems will perform better than the Hubbard model in this respect. Exploiting symmetries and the fact that most vertices actually vanish due to particle and spin conservation, it might be possible to store somewhat larger objects.

Observables
In diagrammatic Monte Carlo, the extraction of observables is typically achieved using a measuring line as illustrated in Fig. 4, see also [31]. One of the lines is then tagged and treated as an entrance and exit of a particle from the system, while the remains of the diagram are interpreted as a contribution to the self energy or the polarization, depending on the line type being considered. In the strong-coupling expansion, the particle propagators are hidden inside the connected vertices, and we only have access to the external lines that originate in the nonlocal processes. Therefore, the Greens function must be obtained from the polarization of the t-line, as At an expansion order Nt, the largest vertex that can be constructed has 2Nt external lines or legs. Using fermionization techniques and conventional diagrammatics, the number of topologies at the same expansion order can be estimated to ∼ 10 6 . opposed to via Dysons equation: where Π is the polarization operator of the t−line. In spin models, two-point correlations can be computed from the polarization of the J-line, while access to further observables that do not correspond to any specific external line can in principle, be obtained by constructing appropriate operators solely for the purpose of measuring. Using multiple measuring lines, many-point correlators can be obtained.

Benchmarks for the Hubbard model
To evaluate the strong-coupling method outlined above, we will now compare it to results obtained with two other state of the art numerical protocols. For the Hubbard model, NLCE can produce unbiased results in the strong-coupling regime, including U = ∞ [26]. This technique is exact in the limit of infinite cluster size, and correspondingly it is also controllable, as convergence of the result with respect to cluster size can be readily checked. A second method that is also applicable in this case is SCT-DMC [35], which is based on spin-charge transformation and a skeleton expansion in the hopping integral t. This also has an additional benefit: Since strongcoupling theory and SCT-DMC share the same expansion parameter and rely on identical skeleton schemes, they are comparable on an order by order basis.
When calculating the strong-coupling expansion, there are two principal computational efforts: Firstly, the connected vertices have to be obtained from the recursion (18). In practice, this set has to be truncated at some given vertex size. We were able to obtain all vertices possessing up to 16 operators attributable to nonlocal processes, which translates to as many external lines, or legs. Since the nonlocal terms ∼ t ij c † i c j each posses two operators, the largest vertex that can be constructed from N t nonlocal terms has 2N t legs, see also Fig.  (3). Secondly, the observables have to be extracted from an expansion in t using the connected vertices. Here, truncation of the total number of nonlocal terms N t is also necessary. We used a standard worm protocol [31] that is very similar to that of [35], and were able to reach an order of up to N t = 9. At this order it is in principle possible to create a connected vertex with as many as 18 legs, and the fact that we had to truncate the vertex size thus affects the last term. Summation was conducted for a chemical potential of µ/t = 2 and an infinite onsite repulsion U = ∞ in a temperature range 1/4 ≤ T /t ≤ 1.
A summary of the result is given in Fig. (5). Comparing these to SCT-DMC data, we obtain a first confirmation that strong-coupling treatment indeed reproduces results from Feynman type diagrammatics. The two series are in good agreement in all cases where results can be obtained to sufficiently high accuracy. What also transpires from this comparison is the disparity in efficiency of the two methods. We are now able to reach N t = 9, while the SCT-DMC results are limited N t = 4 or less. This improvement is very significant considering that the computational complexity scales factorially with expansion order.
To examine the impact of truncating the vertex size, we also compute the equation of state at an order of N t = 7, and vary the maximal vertex size in the span 6 ≤ N legs ≤ 14. At T /t = 1, we estimate that truncation at 8 legs gives an error of less than 3×10 −5 , though no lower bound could be set, as this is within error bars. At T /t = 1/4, we estimate that truncation at 10 legs gives an error of less than 10 −3 , with no lower bar. This implies that the truncation error present at N t = 9 in the results presented here is small compared to the statistical noise. It does also indicate that the current results are limited by the expansion order rather than achievable vertex size.
Extracting estimates for the equation of state from strongcoupling theory, we can compare these to NLCE, see Fig. (6). For temperatures T /t ≥ 1/2, we observe excellent agreement that is within moderate error bars. At the lowest temperature, the uncertainty increases and the NLCE data also begins to diverge, but the results still agree within error bars.
Further improvement to this method can be made by tailoring new sampling protocols specifically to this expansion. Very significant gains can also be made by altering the struc- In (a-d) the filling factor is given as a function of expansion order, in the temperature range 1/4 ≤ T /t ≤ 1 for µ/t = 2 and U/t = ∞. The blue bars correspond to strong-coupling theory (this work), while the red bars were obtained from conventional diagrammatic treatment of the spin-charge transformed Hubbard model [35]. Since both these methods rely on a bold expansion in t we expect them to provide identical results, and this holds true within or almost within error bars. The shaded region gives an estimate of the density at N = ∞, though a precise extrapolation to infinity is not possible at this stage. The vertex size is truncated at 16 legs, which affects the last term when N = 9. In (e-f) we examine the effect of truncating the vertex size at T /t = 1 and T /t = 1/4 respectively. We thus set the expansion order in t to N = 7, and observe how the predicted carrier density varies with the maximal number of vertex legs. At T /t = 1, the corrections beyond 8 legs vanish within the error bars, which are of the order ∼ 3 × 10 −5 . At T /t = 1/4, the corrections beyond 10 legs falls within the error bar of ∼ 10 −3 . Hence, at N = 9, the error due to truncation of the vertex size can be expected to be vanishingly small compared to statistical noise. ture of the series itself. Shifted action [44] or generalizations thereof [37] can be employed to improve the rate of convergence of the series, both with respect to overall expansion order, as well as vertex size.

Summary and outlook
In conclusion, we have derived a diagrammatic technique based on connected strong-coupling vertices, that can be applied to lattice fermions and quantum spin models. This method allows large, or even infinite, contact interactions to be treated non-perturbatively, thus overcoming a longstanding obstacle for diagrammatic methods in the strongly correlated regime. For the Hubbard model, we are able to obtain self consistent solutions up to an expansion order of N t = 9, that display good agreement with results from NLCE.
Experimental progress with strongly correlated systemsusing ultracold atomic gases-is now rapid. With quantum gas microscopy [45,46], key features of the doped Mott insulator can be observed at the single particle level, and at temperatures where spin-correlations become significant [47][48][49][50][51]. Strong coupling Diagrammatic Monte Carlo can be used to extract virtually exact correlators at low temperatures, that can be directly compared to these experiments. frastructure for Computing (SNIC) at the National Supercomputer Centre in Linköping, Sweden. The author would like to thank Marcos Rigol for providing NLCE results for benchmarking as well as Kun Chen, Boris Svistunov and Nikolay Prokof'ev for important input and discussions. nal line, is given by Eq. (20), i.e. n Γ n U n iŌ μ,e = Ō μ+Û .
We begin by noting that the set of operatorsŌ may be expressed in the operator basis (25) as follows: whereŌ α is a set of operators of the form (25). For a finite setŌ, we furthermore note that the range of α is also finite.
Using (26) we obtain where f ({τ i }) gives the time dependence in accordance with (28), and H is expressed in units of temperature. For the Hubbard model, (32) is analytic on the real axis, but not in the entire complex plane due to zeros of the partition function that occur for complex values of U , and so the convergence radius is finite when expanding in contact interactions.
To solve this problem, we use second fermionization to construct a dual representation which is free of a large expansion parameter, and thus possesses a convergent series regardless of model parameters. A detailed discussion of fermionization techniques is given in [32], but we will here recapitulate the central ideas of this approach: Essentially, the goal is to remove the doublons from the trace, and then reintroduce them as hard core bosons that are subsequently fermionized. The end results of this procedure is that the contact interaction becomes a bilinear term in the Hamiltonian.
First, we thus remove the doublons entirely by introducing a projection operator p G and an auxiliary fermionic field with the number operator n A : H = −µ(n e ↑ + n e ↓ ) + p G , p G = n e ↑ n e ↓ iπn A , (33) where n e σ are electron number operators. When we trace over n A = 0, 1, the configurations for which n e ↑ n e ↓ = 1 obtain an imaginary energy shift of 0 or iπ respectively which in turn give them opposite sign in the trace, such that the contribution vanishes.
We then proceed to reintroduce the doublon in the form of a hard core boson, with an energy U − 2µ. The boson can in turn be fermionized, and thus gives rise to two fermionic components with number operators given by n d 0 , n d 1 . The state space correspondence is given by The remaining states in the construction (34) which correspond to n d ↑ + n d ↓ = 1 has no physical counterpart, and are thus removed from the trace by the introduction of a Popov-Fedotov projection term [33] of the form such that the contribution from n d = 0, 2, obtain a complex phase in the in the trace and cancel. Finally, we are required to project out configurations where n e ↑ + n e ↓ = 1, n d 1 = 1, as this has no correspondence in the original state space. This can be achieved by Including also the energy scale of the doublon, we thus arrive at a dual description of the local Hamiltonian according to where E D = U − 2µ is the doublon energy. The partition function of (37) is given by which is indeed the partition function of the Hubbard model in the atomic limit, except for a trivial factor 2 which we obtain when tracing over the auxiliary field. In (37), the contract interaction is described by a bilinear term, and expansion is instead conducted in the projection operators p G , p H . To examine the analyticity of the density matrix as a function of the expansion parameter, we parameterize the expansion terms p G , p H → ξp G , ξp H such that ξ = 1 corresponds to the fully projected system. For convergence of the series, we then require analyticity within the unit circle |ξ| ≤ 1, regardless of model parameters. Next, we recall that the density matrix takes the form For finite model parameters, W i and Z are analytic, implying that the density matrix is also analytic for non-vanishing Z. Correspondingly, demonstrating convergence of the series translates to ruling out zeros of Z(ξ) within the unit circle |ξ| ≤ 1, which we will now do: We begin by expressing the partition function in terms of ξ Z(ξ) = a + be −iπξ + ce iπξ (40) with a = 2e µ−E D +e 2µ−E D +2e −E D +4e µ +e 2µ +2, (41) b = e µ−E D + e 2µ−E D + e 2µ , c = e µ−E D , (42) where in particular we note that Then we observe that since the exponents in (40) are real on the imaginary axis. Furthermore we note that on the real axis, the exponentials in (40) only provide a phase, which together with (43) implies and so there are no poles on the real axis either.
Away from the axes, the partition function is generally complex. For the imaginary part to vanish, we require be πξ I sin(πξ R ) = ce −πξ I sin(πξ R ), ξ = ξ R + iξ I . (46) This equation has two types of solutions: Firstly, we have ξ R = 0, ξ R = ±1, but these lie on the axes since |ξ| ≤ 1.