Analytical solution for time-integrals in diagrammatic expansions: application to real-frequency diagrammatic Monte Carlo

The past years have seen a revived interest in the diagrammatic Monte Carlo (DiagMC) methods for interacting fermions on a lattice. A promising recent development allows one to now circumvent the analytical continuation of dynamic observables in DiagMC calculations within the Matsubara formalism. This is made possible by symbolic algebra algorithms, which can be used to analytically solve the internal Matsubara frequency summations of Feynman diagrams. In this paper, we take a different approach and show that it yields improved results. We present a closed-form analytical solution of imaginary-time integrals that appear in the time-domain formulation of Feynman diagrams. We implement and test a DiagMC algorithm based on this analytical solution and show that it has numerous significant advantages. Most importantly, the algorithm is general enough for any kind of single-time correlation function series, involving any single-particle vertex insertions. Therefore, it readily allows for the use of action-shifted schemes, aimed at improving the convergence properties of the series. By performing a frequency-resolved action-shift tuning, we are able to further improve the method and converge the self-energy in a non-trivial regime, with only 3-4 perturbation orders. Finally, we identify time integrals of the same general form in many commonly used Monte Carlo algorithms and therefore expect a broader usage of our analytical solution.

The past years have seen a revived interest in the diagrammatic Monte Carlo (DiagMC) methods for interacting fermions on a lattice. A promising recent development allows one to now circumvent the analytical continuation of dynamic observables in DiagMC calculations within the Matsubara formalism. This is made possible by symbolic algebra algorithms, which can be used to analytically solve the internal Matsubara frequency summations of Feynman diagrams. In this paper, we take a different approach and show that it yields improved results. We present a closed-form analytical solution of imaginary-time integrals that appear in the time-domain formulation of Feynman diagrams. We implement and test a DiagMC algorithm based on this analytical solution and show that it has numerous significant advantages. Most importantly, the algorithm is general enough for any kind of single-time correlation function series, involving any single-particle vertex insertions. Therefore, it readily allows for the use of action-shifted schemes, aimed at improving the convergence properties of the series. By performing a frequency-resolved action-shift tuning, we are able to further improve the method and converge the self-energy in a non-trivial regime, with only 3-4 perturbation orders. Finally, we identify time integrals of the same general form in many commonly used Monte Carlo algorithms and therefore expect a broader usage of our analytical solution.
In Refs 42 and 51 it has been shown that a convenient transformation of the interaction-expansion series can be used to significantly improve its convergence and sometimes allows one to converge the electronic selfenergy with only a few perturbation orders where it would have otherwise been impossible. The method relies on a transformation of the action which affects the bare propagator at the cost of an additional expansion, i.e. more diagram topologies need to be taken into account. Alternatively, this transformation can be viewed as a Maclaurin expansion of the bare propagator with respect to a small chemical potential shift. The resulting convergence speed up comes from an increased convergence radius of the transformed series.
In a separate line of work, DiagMC methods have been proposed that are based on the Matsubara formalism that do not require an ill-defined analytical continuation 46 . Such methods have so far been implemented for the calculation of the self-energy 47,48 and the dynamical spin susceptibility 49 . The algorithms differ in some aspects but all rely on the symbolic algebra solution of the internal Matsubara frequency summations appearing in Feynman diagrams. However, this approach has some downsides. First, numerical regulators are needed to properly evaluate Bose-Einstein distribution functions and diverging ratios that appear in the analytical expressions, and also poles on the realaxis (effective broadening of the real-frequency results). In the case of finite cyclic lattice calculations, multiple precision algebra is needed in order to cancel divergences even with relatively large regulators. 47 Most importantly, in the Matsubara summation algorithm, applying the series transformation from Refs 42 and 51 would require a separate analytical solution for each of the additional diagram topologies, that are very numerous, and the calculation would become rather impractical. More generally, treating any distinct diagram requires that the Matsubara frequency summations be performed algorithmically beforehand. This makes it virtually impossible to devise MC sampling algorithms that go to indefinite perturbation orders because the arXiv:2011.08226v2 [cond-mat.str-el] 19 Nov 2020 upper limit of the perturbation orders must be set in advance.
In this paper, we show that it can be advantageous to start from the imaginary-time domain formulation of Feynman diagrams. A diagram contribution then features a multiple imaginary-time integral, rather than sums over Matsubara frequencies. The multiple integral can be solved analytically and we present a general solution. This analytical solution, although equivalent to the analytical Matsubara summation, has a simpler and a more convenient form that does not feature Bose-Einstein distribution functions or diverging ratios. As a result, numerical regulators are not needed and the need for multiple precision arithmetic may arise only at very high perturbation orders. The numerical evaluation yields a sum of poles of various orders on a uniform grid on the real axis. The ability to separate contributions of poles of different orders allows one to formally extract the real-frequency result without any numerical broadening. Finally, the analytical solution is general and applies to all diagram topologies that would appear in the transformed series proposed in Refs 42 and 51 or any other diagrammatic series for single-time correlation functions. This paves the way for real-frequency diagrammatic algorithms formulated in real-space that are not a priori limited to small perturbation orders (similarly to CTINT or CTAUX 41 ). In this work, we apply the analytical time-integral to the momentum-space DiagMC for the calculation of the selfenergy, and implement and thoroughly test the method. We reproduce the self-energy results from Ref. 51 and supplement them with real-axis results, free of the uncontrolled systematic error that would otherwise come from the analytical continuation. Furthermore, we show that even if a full convergence is not possible with a single choice of the action-tuning parameter, one can choose the optimal tuning parameter for each frequency independently 45 . Such a frequency-resolved resummation can be used to improve the solution and in some cases systematically eliminate the non-physical features that appear in the result due to the truncation of the series at a finite order. The paper is organized as follows. In Section I, we define the model and the basic assumptions of our calculations. In Section II, we introduce our method in details. First, in Section II A, we present the analytical solution of the general multiple-time integral that appears in time-domain formulation of Feynman diagrams and discuss the numerical evaluation of the final expression. Then in Section II B, we show the analytical solution for the Fourier transform of the Maclaurin expansion of the bare propagator, which is essential for our DiagMC algorithm. In Section II C, we discuss in detail how our analytical solutions can be applied in the context of Di-agMC for the self-energy. In Section III, we discuss our results and benchmarks and then give closing remarks in Section IV. Additional details of analytical derivations and further benchmarks and examples of calculations can be found in the appendices.

I. MODEL
We solve the Hubbard model given by the Hamiltonian (1) where σ ∈ {↑, ↓}, i, j enumerate lattice sites, t ij is the hopping amplitude between the sites i and j, U is the onsite coupling constant, and µ is the chemical potential. We only consider the Hubbard model on the square lattice with the nearest-neighbor hopping t and next-nearest-neighbor hopping t . The bare dispersion is given by We define D = 4t, which will be used as the unit of energy unless stated otherwise. We restrict to thermal equilibrium and paramagnetic phases with full lattice symmetry.

II. METHODS
The idea of DiagMC algorithms is to stochastically compute the coefficients of a perturbation series describing some physical quantity. We will focus on expansions in the coupling constant U and a shift in the chemical potential δµ. The calculation of each coefficient involves the evaluation of many Feynman diagrams expressed in terms of the bare propagator, in our case taken as a function of momentum and two imaginary times. The evaluation of a diagram then boils down to a sum over multiple momentum variables and a multiple imaginary time integral that is always of the same generic form. The goal of this Section is to find a general analytical solution for these time integrals and reformulate the perturbation series as a function of a complex frequency z.

A. Analytical solution of time-integrals
We are interested in solving analytically integrals of the form with τ K+1 = β, and iΩ η any Matsubara frequency, either fermionic or bosonic, depending on η (iΩ η=−1 ≡ iω ≡ i(2m + 1)πT and iΩ η=1 ≡ iν ≡ 2imπT , with m ∈ Z). The parameters ω i may be complex, l i and K are integers, l i ≥ 0, K ≥ 2; r can be any integer. δ x,y is the Kronecker delta (it will be used throughout this paper, also in the shortened version δ x ≡ δ x,0 ). The reason for our choice to label times starting from 2 will become clear later.
The main insight is that upon applying the innermost integral, one gets a number of terms, but each new integrand has the same general form ∼ τ n e τ z . The solution therefore boils down to a recursive application of The number of terms obtained after each integration is apparently 1 + (1 − δ z )(n + 1), and we can enumerate all terms obtained after the full integration by a set of integers {k i }, where k i ≥ 0 denotes the choice of the term of the integral i (over dτ i ). For a given choice of k i , the propagation of exponents (n and z in Eqs 4 and 5) across successive integrals can be fully described by a simple set of auxiliary quantities. The exponent of e in the integration i we denote asz i and it is given bỹ where we introduced b i ≡ b ni,ki , and n i is specified below. Alternatively one can obtain the exponentz i as and The exponent of τ that will be carried over from integration i to integration i + 1 depends on the choice of the term from the integral i, and is given by Pos(n i −k i ) where Pos denotes the positive part of the number (Pos(x) = (x + |x|)/2). n i denotes the maximum exponent that can be carried over from integration i, and is obtained as: and In the case of Eq.4, the maximal exponent that can be carried over to the next integration coincides with the exponent that entered the integral (the integral Eq.4 does not raise the power of τ ), so the definition of n i coincides with the meaning of n in Eq.4. In the case of integral Eq.5, n i rather denotes the exponent after the integration, i.e. n + 1.
The solution for the integral can be in the end continued to the whole of the complex plane iΩ η → z, and can be written down as (introducing shorthand notation X = (r, {l 2 ...l K }, {ω 2 ...ω K })): Note that we have rewritten the sum over {k i } as a sum over {b i } and a partial (inner) sum over {k i }. This is not necessary, being that b i is a function of k i . Each b i is fully determined by k i , but not the other way around, so the inner sum over k i in Eq.13 goes over values that are allowed by the corresponding b i . We present this form of Eq.13 to emphasize that the factor e b K βω K depends only on {b i }, and can thus be pulled out of the inner {k i } sum. The notation "k i : b i = 1" means that we sum only over k i with i such that b i = 1. The remaining k i are fixed to n i + 1, which is the only possibility if b i = 0. The notation is applied analogously in other products over i.
The only remaining step is to expand the product of poles in Eq.13 into a sum of poles Note also that if r / ∈ [2, K] (no Matsubara frequency appearing in any exponent), the result of the integral is a number, rather than a frequency dependent quantity. In that case, the integral can be straightforwardly generalized to the case of real time, where integrations go to some externally given time t (instead of β), and the resulting expression is a function of that time. The step Eq.14 is then not needed. See Appendix A for details.

Numerical evaluation and relation to other algorithms
The implementation of Eq.13 is rather straightforward and much simpler than the algorithmic Matsubara summations in our previous work Ref. 47. Indeed, most of the calculations just require the numerical evaluation of an analytical expression and it is not necessary to implement a dedicated symbolic algebra to manipulate the expressions. The only exception is the last step Eq.14. This transformation was the centerpiece of the algorithm in Ref. 47 and was applied recursively many times, leading to complex book-keeping and data structures. Ultimately, the result was a symbolic expression that was stored, and a separate implementation was needed for the comprehension and numerical evaluation of such a general symbolic expression. In the present context, however, Eq.14 is applied only once to produce numbers, and is simple to implement. The other important point is that we treat analytically cases with δz i = 1 by employing Eq.5. With the frequency-summation algorithms 47,48 , one cannot take into account possible cancellations of ω i terms in Eq.9 without computing a large number of separate analytical solutions. When untreated, these cancellations yield diverging ratios in the final expressions, which need to be regularized. On the contrary, in Eq.13 the ratio 1/ω ki+bi i cannot have a vanishing denominator and its size will in practice be limited by the energy resolution. This will also allow us to have the final result in the form of a sum of poles on an equidistant grid on the real-axis, and extract the real-axis results without any numerical pole-broadening (see Section II C 2 and Appendix B).

B. Expansion of the bare propagator
The central quantity is the Green's function defined in Matsubara formalism as: The non-interacting Green's function (or the bare propagator) in the eigenbasis of the non-interacting Hamiltonian has a very simple general form and for the plane-wave k, the propagator is As we will discuss below, the diagrammatic series for the self-energy will in general be constructed from different powers of the bare propagator: Indeed, these powers naturally arise after expanding the bare propagator in a Maclaurin series 1 z+x = ∞ n=0 (−x) n z n+1 around a small chemical potential shift This series converges (for all iω) if δµ is smaller in amplitude than the first Matsubara frequency: |δµ| < πT . Nevertheless, this expression will become a part of a larger series with additional expansion parameters, which may result in a modified convergence radius of the overall series with respect to δµ. We anticipate that the Feynman diagrams will be formulated in the imaginary-time domain, so it is essential to work out the Fourier transform of G l 0 (ε, iω). We present the full derivation in Appendix E and here only write the final solution with s τ,τ = sgn(τ − τ ). In our notation, l in G l 0 is a superscript index, rather than the power of G 0 (although these meanings coincide in the case of G l 0 (ε, iω)). The Fermi function is defined as n F (ε) = 1/(e βε + 1) and the coefficients that go with τ ζ τ ς terms are and Here we make use of binomial coefficients n k = n! k!(n−k)! and the Stirling number of the second kind n k = k i=0 C. Application to DiagMC

Hartree-shifted series
In this section we discuss the construction of the selfenergy series, where all tadpole-like insertions are omitted in the topologies of the diagrams. Rather, the full Hartree shift is absorbed in the bare propagator. The diagrams are therefore expressed in terms of the Hartreeshifted bare propagator with the Hartree-shifted dispersion defined as where nσ is the average site-occupation per spin. After constructing the tadpole-less topologies, we are free to expand all propagators that appear in the diagrams according to Eq.19: In frequency domain, this step can be viewed as introducing new topologies: we now have diagrams with any number of single-particle-vertex (δµ) insertions on any of the propagator lines. Each arrangement of these additional single-particle vertices on the diagram does require a separate solution by the symbolic algebra algorithm as presented in Refs. 47 and 48. Nevertheless, as a δµ-vertex cannot carry any momentum or energy, the formal effect of it is that it just raises the power l of the propagator that passes through it. In the imaginarytime domain, it turns out that the contribution of δµdressed diagrams is readily treatable by the analytical expression Eq.13, and we no longer have to view the δµ-insertions as changes to topology, but rather as additional internal degrees of freedom to be summed over. This is illustrated in Fig.1.
Up to the Hartree shift, the self-energy expansion can now be made in powers of the interaction U and the small chemical-potential shift δµ where j enumerates the propagators, of which there are N prop = 2N − 1, N is perturbation order in U , each l j goes from 1 to ∞, Υ N enumerates distinct topologies of the diagram at order N (without any δµ or Hartree insertions), and D is the contribution of the diagram. The general form of the diagram contribution is We denote N bub the number of closed fermion loops in the diagram; τ 1 ...τ N are internal times, and we fix τ i=1 = 0; τ is the external time; k is the external momentum, k 1 ..k N are the independent internal momenta; j indexes the propagator lines, andk are corresponding linear combinations of the momentã k j ≡ N λ=0s jλ k λ , wheres jλ ∈ {−1, 0, 1} and we index with 0 the external momentum k 0 ≡ k.τ j andτ j are outgoing and incoming times for the propagator j, and take values in {τ 1 ...τ N }, where we denote with index N the external time τ N ≡ τ . The coefficientss jλ , times τ j ,τ j and the number N bub are implicit functions of the topology Υ N . Throughout the paper, we assume normalized k-sums, k ≡ 1 N k k , where N k is the number of lattice sites. We can perform the Fourier transform of the external time, to obtain the contribution of the diagram in the Matsubara-frequency domain: The Green's function G l 0 (ε, τ − τ ) is discontinuous at τ = τ , so to be able to perform the τ -integrations analytically, we first need to split the integrals into ordered parts where P denotes all (N − 1)! permutations of time indices. p labels the permutation and p i is the permuted index of vertex i. Let us rewrite the contribution of the diagram, with propagators written explicitly using the expression Eq.20: where J i/o (i) is the set of incoming/outgoing propagators j of the vertex i, which depends on the topology Υ N . We also introduced shorthand notation s j = sτ j ,τ j . Practically, s j depends on whether p(i(j)) > p(i (j))) or the other way around, where i(j)/i (j) is the outgoing/incoming vertex of propagator j in the given permutation p. The total number of forward-facing propagators is N fwd (p) = j δ −1,sj , which depends on the permutation and the topology. The products of δ ζj and δ ςj are there to ensure that the time τ 1 = 0 is not raised to any power other than 0, as such terms do not contribute.
Now we can apply the analytic solution for the time integrals (Eq.13) to arrive at the final expression: where i(p i ) is the vertex index i of the permuted index p i and we have introduced a new expansion variable L = j (l j − 1) and a convenient variablel j = l j − 1, so that which is the series we implement and use in practice.

Numerical evaluation
The expression Eq.30 is very convenient for numerical evaluation.
First, we restrict the values ofε k to a uniform grid on the real-axis with the step ∆ω (ε k = j∆ω). These appear in ω 2 , ..., ω K as terms with integer coefficients, which means that {ω i } entering I X will also be restricted to the same uniform grid. The final result there-fore has the form: This form allows us to reinterpret the finite-lattice results as that of the thermodynamic limit and extract D Υ N ,k,L,δµ (ω + i0 + ) without any numerical broadening (see Appendix B for details).
In our present implementation we perform a flat-weight (uniform) MC sampling over internal momenta {k i }, and do a full summation of all the other sums, and accumulate the amplitudes A j,p . There are, however, other options. For example, one may sample {k i }, {p i }, {b i } and use P ≡ j n F (s jεk j )e b N βω N as the weighting function. We have checked thoroughly that the factor P correlates closely with the contribution to A j,p coming from a given choice of {k i }, {p i }, {b i } variables (with other variables summed over), and thus P could be a good choice for a weighting function. However, this requires additional operations related to move proposals and trials, and we have not yet been able to make such an algorithm more efficient than the flat-weight MC. Nevertheless, it is apparent that our approach offers more flexibility than the algorithmic Matsubara summations (AMS). In AMS no convenient weighting function can be defined for the Monte Carlo, so one either does the flat-weight summation 47 or uses the whole contribution to the result as the weight, which comes at the price of having to repeat the calculation for each frequency of interest 46 (on the contrary, in Ref. 47, as well as in this paper, the entire frequency dependence of self-energy is obtained in a single MC run). Concerning floating point arithmetic, it is important that the factor e b N βω N stemming from I X can always be absorbed into the product of n F functions in the second row of Eq.30. This can be understood as follows. A givenεk j can at most appear twice as a term inω N , once with sign +1 and once with sign −1, corresponding to the incomingτ j and outgoingτ j ends of the propagator j. In that case, the exponent cancels. The other possibility is that it appears only once, in which case it must correspond to the later time in the given permutation. If the later time is the outgoing end of the propagator, then the propagator is forward facing, the sign in front is s = −1; if it's the incoming end, then the propagator is backward facing, and the sign in front is s = 1. In both cases we can make use of e sβε n F (sε) = n F (−sε) Therefore, no exponentials will appear in the final expression. A product of n F functions is at most 1, and the coefficients c are not particularly big. Then, the size of the pole amplitudes that come out of Eq.13 is determined by the energy resolution (1/∆ω) and temperature (β n K +1−b K −k K ). In our calculations so far, the amplitudes remain relatively small. Our approach ensures that we do not have very large canceling terms, like we had in Ref. 47. Indeed, we have successfully implemented Eq.30 without the need for multiple-precision floating-point types.
It is interesting to compare the computational effort for the numerical evaluation of our analytical solution to the straightforward numerical integration. In the most straightforward integration algorithm, one would discretize the imaginary-time interval [0, β] with N τ times, and then perform the summation which has the complexity O(N N −1 τ ) for each external τ , so overall O(N N τ ). With our algorithm we do not have to go through all configurations of internal times, but we do need to go through all possible permutations of the internal times, and for each permutation there is at least 2 N −1 terms to be summed over. So the number of terms one has to sum grows at least as O((N − 1)!2 N −1 ). At sufficiently high N , this number is bound to outgrow the exponential N N τ , whatever the N τ . This will happen, however, only at very large N . For example, if N τ = 30, the analytical solution becomes slower at around N = 40. Moreover, one actually needs a much larger N τ , especially at low temperature. In any case, the additional computational effort can be understood as coming from the difference in the information content of the result, which is a lot more substantial in the case of the analytical solution. At orders N < 6, we find that the implementation of our new algorithm is significantly more efficient than our current implementation of the Matsubara summations from Ref. 47, and at N = 6 they are about equally efficient. However, we anticipate that further optimizations will be possible at the level of Eq.13 and Eq.30.

Bare series
We are also interested in constructing a bare series where tadpole insertions are present in diagram topologies. Tadpole (or Hartree) insertions are static and an evaluation of their amplitudes can be done relatively simply by various means. Since Hartree insertions have the same effect on the underlying diagram topology as the δµ-insertions, our expression Eq.30 is already general enough: we do not actually have to consider Hartree dressed topologies. Inclusion of the Hartree insertions can be entirely accounted for in the resummation of the D Υ N ,k,L,δµ (z) contributions from the previous section, with the replacementε (i.e. full Hartree shift excluded). We will redefine the meaning of the expansion parameter K to be the total number of Hartree and δµ in- sertions, and consider the series only up to K = 5. In that case, we can have at most three interaction vertices in a Hartree insertion, which leaves only 5 such diagrams (Fig.2). We can evaluate these 5 amplitudes with very little effort, by making use of spatial and temporal Fourier transforms. Note that the expansion of the propagators in δµ is performed in Hartree insertions as well, so we need to account for possible additional δµ insertions inside the Hartree diagrams. We first define the bare density and the real-space propagator: We will also need the bare polarization χl 1 ,l2 0,q=0 (iν = 0) = r dτ χl 1 ,l2 0,r (τ ) (38) and the second-order self-energy which can be Fourier transformed to yield Σl 1,l2 ,l3 We can now calculate the amplitudes of the possible Hartree insertions with a number L of δµ insertions on them, in any arrangement The series can be now resummed as: where Ω(L, {M L i }) is the combinatorial prefactor which counts all possible ways the selected single-particle vertices δµ,{D i } can be arranged. This corresponds to the number of permutations of multisets We emphasize that Eq. 45 is fully general, but at orders K ≥ 5, additional Hartree insertions D (compared to  need to be considered.

A. Convergence speed-up with δµ expansion in the bare series
Here we focus on supplementing the results from Ref.51 with real-frequency self-energies calculated without any numerically ill-defined analytical continuation. The model parameters are t = −0.3t, µ = 0, U = 1.0D, T = 0.125D, and n σ = 0.3625. In Ref. 51, the calculation was performed with the Hartree shifted series with δµ = 0, as well as with the bare series, with two values of δµ, namely 0.15D and 0.3825D. We repeat these calculations with our method. We use lattice size 32 × 32, and project the dispersion on a uniform energy grid as described in Ref. 47, and discussed in Section II C 2. In Fig.3 we show our results, and compare them with the results of Ref. 51.
In the upper row of Fig.3 are the real-frequency selfenergies calculated up to order K ≤ 5. We are keeping a finite broadening η = 0.3D, to smoothen the curves. As discussed in Appendix B, in our method, numerical pole-broadening is not a formal necessity. However, there is still a significant amount of statistical noise in our real-frequency result (although the imaginaryfrequency result is already very well converged). It is important to note that some of the noisy features in our real-frequency result may be artifacts of the finite lattice size that would not vanish with increasing number of MC steps. However, by comparing the result with a 256 × 256 lattice calculation (Appendix C) we check that at already at η = 0.2D, no such artifact should be visible. Our method appears rather insensitive to the lattice size, i.e. the lattice size can be increased arbitrarily at no apparent cost: at η = 0.2D, the 256 × 256 lattice calculation appears equally well converged as the 32 × 32 lattice calculation with the equal number of MC steps, and yields a result that is on top of the 32 × 32 calculation.
In the bottom row of Fig.3  The excellent agreement with the results from Ref. 51 serves as a stringent test of our implementation. In the δµ = 0.3825D calculation, even on the real-axis, the self-energy does appear well converged by order K = 5, although there is some discrepancy between K = 4 and K = 5 at around ω = 1.5D.

B. ω-resolved resummation
We can now go one step further by resumming the series presented in Fig.3a and Fig.3c for each ω individually, using an ω-dependent optimal shift δµ * (ω).
The results are shown in Figs.4 and 5.
We determine the optimal δµ * (ω) by minimizing the spread of the ImΣ(ω + iη) results between orders K = 3 and K = 5. This spread as a function of ω and δµ is color-plotted in Figs.4 and 5. We have results for a discrete set of δµ ∈ {δµ i }, so the optimal δµ * (ω) is a priori a discontinuous curve. As this is clearly nonsatisfactory, we smoothen the curve (shown with the blue line on the top panels in Figs.4 and 5). Clearly, we do not have results for each precise value of this optimal δµ * (ω). One could take for each ω the available δµ i that is closest to δµ * (ω), but this would, again, result in a discontinuous curve. To avoid this, we average the available results as The results of the averaging around the optimal δµ * (ω) are shown in the middle and bottom panels of Figs.4 and 5. In both cases, the ω-resolved resummation helps to converge the result. In the case of the bare series, the convergence is now almost perfect, and already order K = 3 is on top of the exact result. In the case of the Hartree shifted series, the results are not perfectly converged at ω < 0, yet the K = 5 calculation is practically on top of the exact result on the imaginary axis, and presents an improvement to the δµ = 0 series in Fig.3a. Note that the improvement in convergence is seen on the imaginary axis, as well.

C. Removing non-physical features
In this section we focus on the parameters case discussed in Ref. 47. We calculate the Hartree-shifted series with parameters of the model t = 0, µ − U n σ = −0.1D, T = 0.1D, and employ various δµ shifts. The lattice size is again 32 × 32 and we focus on the self-energy at k = (0, π).
The results are presented in Fig.6 for three values of U .
At low U , the series is well converged by K = 5, and the result is entirely insensitive to the choice of δµ, as expected.
At intermediate and high U , the result can be strongly δµ sensitive. The δµ dependence of the result, however, strongly varies with ω. It appears that for a given ω, there are ranges of δµ value where the result (at fixed order K) is insensitive to the precise choice of δµ. This presents an alternative way of choosing an optimal δµ (similar idea was employed in a different context in Ref. 57).
The striking feature at large U is the causality violations at |ω| ≈ 2 that were previously discussed in Ref. 47 (note that the broadening somewhat masks the extent of the problem). The dips in the self-energy spectrum appear to happen only at certain values of δµ: at ω = −2, the problem is present at δµ large and negative, and at ω = 2 at δµ large and positive. In particular at ω = 2, the result appears to vary uniformly with δµ, and one cannot select an optimal δµ based on sensitivity of the result to the δµ value. We therefore repeat the procedure from the previous section and select the optimal δµ * (ω) based on level of convergence between orders K = 4 and K = 5. The spread of results and a smooth choice of δµ * (ω) are presented in Fig.7. In Fig.8, the results of the averaging are shown and compared to δµ = 0 results at the highest available orders K = 4 and K = 5, at three values of U . The conver- gence is visibly better around our δµ * than with δµ = 0 at problematic frequencies |ω| ≈ 2. More importantly, the non-physical features are clearly absent. At U = 1, in the δµ = 0 calculation, the causality is not yet violated, but the dip at ω = 2 is already starting to appear, which is clearly an artifact of the series truncation which should be removed systematically. It is important that the intermediate frequency behavior we obtained by averaging results around the optimal δµ is indeed the cor- rect one, and it will not change much further with increasing orders. We show in the top panels the K = 6, the δµ = 0 result which has been benchmarked against a fully converged imaginary-axis result in Fig.9 (the converged result was obtained with the ΣDet method 58,59 at order 8). Clearly, the improved convergence between orders 4 and 5 that we have achieved by choosing δµ appropriately, does indeed mean an improved final result. However, our procedure does not improve the result at around ω = 0 where the optimal δµ does appear to be close to 0. The K = 6, δµ = 0 result shown in the upper panels of Fig.8 is still a bit different from K = 5,   Fig.4a, for the parameters of the model corresponding to Fig.6. The blue line is the optimal δµ * , to be used in Fig.8.
δµ ≈ δµ * (ω) results around ω = 0. In the case of U = 1D, it is interesting that a large negative δµ does bring the ω ≈ 0 result at order K = 5 much closer to the the exact value. This can be anticipated from Fig.6 where we show corresponding results for U = 0.8D and U = 1.2D. Also, by looking at the color plot in Fig.7, we see that at ω = 0, there is indeed a local minimum in the spread at around δµ = −0.2 which could be used as the optimal δµ * . This minimum, however, cannot be continuously connected with the other minima that we observe at ω < 0, so we chose a different trajectory in the (ω, δµ)-space. It would be interesting for future work to inspect the behavior at even more negative δµ, where another continuous trajectory δµ * (ω) might be found.

IV. DISCUSSION, CONCLUSIONS AND PROSPECTS
In this paper we have derived an analytical solution for the multiple-time integral that appears in the imaginary-time Feynman diagrams of an interaction series expansion. The solution is general for any diagram with a single external time, or no external times. We find this generality to be a great advantage compared to the recently proposed algorithmic solutions of the corresponding Matsubara-frequency summations. Our analytical solution allowed us to develop a very flexible DiagMC algorithm that can make use of the possibility to optimize the series with shifted actions. As a result we were able to almost perfectly converge a realfrequency self-energy in just 3-4 orders of perturbation, in a non-trivial regime and practically in the thermody- namic limit.
More importantly, the fact that one does not have to prepare a solution for each diagram topology individually opens the possibility to develop algorithms more akin to CTINT and allow the MC sampling to go to indefinite perturbation orders. In fact, upon a simple inspection of CTINT and segment-CTHYB equations 41 , it becomes clear that our solution can in principle be applied there, so as to reformulate these methods in real-frequency. This would, however, come at the price of having to break into individual terms the determinant that captures all the contributions to the partition function at a given perturbation order. In turn, this may lead to a more significant sign problem, and an effective cap on the perturbation orders that can be handled in practice. On the other hand, it is not entirely clear how much of the sign problem comes from summing the individual terms, and how much from the integration of the internal times, and we leave such considerations for future work. In any case, DiagMC algorithms based on hybridization expansion have been proposed before (see Ref. 23, 28, and 60), where our analytical solution may be applied.
Our solution also trivially generalizes to real-time integrals and may have use in Keldysh and Kadanoff-Baym 9 calculations, where the infamous dynamical sign problem arises precisely due to oscillating time integrands. There have been recent works 61,62 with imaginary-time propagation of randomized walkers where our solution may also find application.

Appendix A: Real-time integration
Let's consider the following special case of the integral Eq.3, which is relevant for real-time integrations featuring integrands of the form e itE : with t K+1 ≡ t. This corresponds to the case r / ∈ [2, K] in Eq.3, and ω i = iE i , and we will defineẼ i analogously toω i . The result is then obtained straightforwardly from Eq.13 which has the following general form Appendix B: Extracting real-axis results without pole-broadening In this section we show how the results on the realaxis can be extracted without any numerical broadening of the poles. Rather, we make use of the pole amplitudes, by interpreting the result as being representative of the thermodynamic limit, where poles on the realaxis merge into a branch cut, and we consider that the pole amplitude is then a continuous function of the realfrequency. We extract the imaginary part of the contribution, and then the Hilbert transform can be used to reconstruct the real part. In case of simple poles, the procedure is trivial, as the pole amplitude at ω is simply proportional to its contribution to the imaginary part at ω + i0 + . When we also have higher order poles, we must use the general expression: ImD(j∆ω + i0 + ) = − π ∆ω p (−1) p−1 p!∂ p−1 A j,p (B1) where∂ p−1 is the numerical (p − 1)-th derivative with respect to j∆ω. The energy resolution is therefore a Diagram used is the second order diagram, with L = 2 (illustrated at the top). Top three panels are contributions from 1st, 2nd and 3rd order poles, respectively. Bottom panel is the total result. Lines with η > 0 are obtained with numerical broadening. The crosses on the η = 0 result denote the available frequencies (in between we assume linear interpolation).
measure of the systematic error made in this procedure. Furthermore, it is necessary that the pole amplitudes form a reasonably smooth function of real-frequency. To avoid statistical noise and noisy features coming from the finite size of the lattice (see next section), we test our method on the example of a N = 2, L = 2 diagram, which we can solve with full summation of Eq.30, on a lattice of the size 96 × 96. This diagram produces poles up to order 3. The result is shown in Fig.10. In the first three panels we show contribution from poles of each order, and in the bottom panel we show the total result.

Appendix C: Convergence with lattice size
In this section we discuss the convergence of the result with respect to the lattice size.
In Fig.11 we compare the results for a single N = 5, L = 0 diagram on the lattices of size 32 × 32 and 256 × 256. We observe that the result is almost exactly the same at broadening level η = 0.2, which brings further confidence in the results in the main part of the paper. In Fig.12 we illustrate how the size of the lattice determines the highest energy-resolution one can have, under requirement that the results form a continuous curve on the real axis and are, therefore, representative of the thermodynamic limit. We perform the full summation for the second order diagram with L = 0, with various sizes of the lattice and various resolutions. Clearly, the bigger the lattice, the higher the energy-resolution one can set without affecting the smoothness of the results.
The numerical parameters of the calculation are therefore the size of the lattice, the energy resolution and the broadening (resolution and the broadening can be tuned a posteriori), and one can tune them to get the optimal ratio between performance and the error bar. If the pole amplitudes A jp are a relatively smooth function of j, no broadening is then needed at all.