Marjan Maček

Philipp T Dumitrescu
email: pdumitrescu@flatironinstitute.org

Corentin Bertrand

Bill Triggs

Olivier Parcollet

Xavier Waintal
email: xavier.waintal@cea.fr

Quantum Quasi-Monte Carlo Technique for Many-Body Perturbative Expansions

High order perturbation theory has seen an unexpected recent revival for controlled calculations of quantum many-body systems, even at strong coupling. We adapt integration methods using low-discrepancy sequences to this problem. They greatly outperform state-of-the-art diagrammatic Monte Carlo. In practical applications, we show speed-ups of several orders of magnitude with scaling as fast as 1/N in sample number N ; parametrically faster than 1/ √ N in Monte Carlo. We illustrate our technique with a solution of the Kondo ridge in quantum dots, where it allows large parameter sweeps.

The exponential complexity of quantum many-body systems is at the heart of many remarkable phenomena. Advances in correlated materials and recently developed synthetic quantum systems -e.g. atomic gases [START_REF] Gross | Quantum simulations with ultracold atoms in optical lattices[END_REF], trapped ions [START_REF] Blatt | Quantum simulations with trapped ions[END_REF], and nanoelectronic devices [START_REF] Goldhaber-Gordon | Kondo effect in a single-electron transistor[END_REF][START_REF] Goldhaber-Gordon | From the Kondo regime to the mixed-valence regime in a single-electron transistor[END_REF][START_REF] Cronenwett | A tunable Kondo effect in quantum dots[END_REF][START_REF] Iftikhar | Tunable quantum criticality and superballistic transport in a "charge" Kondo circuit[END_REF] -have allowed many-body states to be characterized and controlled with unprecedented precision. The latest of these systems, quantum computing chips, are highly engineered out-ofequilibrium many-body systems, where the interacting dynamics performs computational tasks [START_REF] Bernien | Probing many-body dynamics on a 51-atom quantum simulator[END_REF]. However, our understanding of these many-body systems is limited by their intrinsic complexity. While uncontrolled approximations can give insight into possible behaviors, there is a growing effort to develop controlled, high-precision methods [START_REF] Leblanc | Solutions of the two-dimensional Hubbard model: Benchmarks and results from a wide range of numerical algorithms[END_REF], especially ones that apply far from equilibrium [START_REF] Cohen | Taming the dynamical sign problem in real-time evolution of quantum many-body problems[END_REF][START_REF] Rosario | Quantum Monte Carlo for correlated out-of-equilibrium nanoelectronic devices[END_REF][START_REF] Bertrand | Reconstructing nonequilibrium regimes of quantum many-body systems from the analytical structure of perturbative expansions[END_REF]. These allow us to make quantitative predictions about the physics of many-body systems and to uncover qualitatively new effects at strong coupling.

Among theoretical approaches, perturbative expansions in the interaction strength have seen an unexpected recent revival, in particular using a family of "diagrammatic" Quantum Monte Carlo (DiagQMC) methods [START_REF] Rosario | Quantum Monte Carlo for correlated out-of-equilibrium nanoelectronic devices[END_REF][START_REF] Bertrand | Reconstructing nonequilibrium regimes of quantum many-body systems from the analytical structure of perturbative expansions[END_REF][START_REF] Prokof | Polaron problem by diagrammatic quantum Monte Carlo[END_REF][START_REF] Prokof | Bold diagrammatic Monte Carlo: A generic sign-problem tolerant technique for polaron models and possibly interacting many-body problems[END_REF][START_REF] Mishchenko | Comprehensive study of Fröhlich polaron[END_REF][START_REF] Van Houcke | Feynman diagrams versus Fermi-gas Feynman emulator[END_REF][START_REF] Wu | Controlling Feynman diagrammatic expansions: Physical nature of the pseudogap in the two-dimensional Hubbard model[END_REF][START_REF] Rossi | Determinant diagrammatic Monte Carlo algorithm in the thermodynamic limit[END_REF][START_REF] Chen | A combined variational and diagrammatic quantum Monte Carlo approach to the many-electron problem[END_REF][START_REF] Bertrand | Quantum Monte Carlo algorithm for out-of-equilibrium Green's functions at long times[END_REF][START_REF] Moutenet | Cancellation of vacuum diagrams and the long-time limit in out-of-equilibrium diagrammatic quantum Monte Carlo[END_REF][START_REF] Rossi | Renormalized perturbation theory at large expansion orders[END_REF]. Using various techniques [START_REF] Bertrand | Reconstructing nonequilibrium regimes of quantum many-body systems from the analytical structure of perturbative expansions[END_REF][START_REF] Prokof | Polaron problem by diagrammatic quantum Monte Carlo[END_REF][START_REF] Chen | A combined variational and diagrammatic quantum Monte Carlo approach to the many-electron problem[END_REF][START_REF] Rossi | Renormalized perturbation theory at large expansion orders[END_REF], it is now possible to sum perturbative series beyond their radius of convergence and thus access strongly correlated regimes. The effects of strong interactions have been studied in diverse systems, including unitary quantum gases [START_REF] Van Houcke | Feynman diagrams versus Fermi-gas Feynman emulator[END_REF], polarons [START_REF] Prokof | Polaron problem by diagrammatic quantum Monte Carlo[END_REF], quantum dots [START_REF] Rosario | Quantum Monte Carlo for correlated out-of-equilibrium nanoelectronic devices[END_REF][START_REF] Bertrand | Reconstructing nonequilibrium regimes of quantum many-body systems from the analytical structure of perturbative expansions[END_REF][START_REF] Bertrand | Quantum Monte Carlo algorithm for out-of-equilibrium Green's functions at long times[END_REF], and pseudo-gap metals [START_REF] Wu | Controlling Feynman diagrammatic expansions: Physical nature of the pseudogap in the two-dimensional Hubbard model[END_REF].

DiagQMC is currently the preferred strategy for computing series coefficients at large perturbation order n, as this involves integrals of dimension proportional to n (practically around 5 -30). High dimensional integration is notoriously difficult, and Monte Carlo provides a robust and flexible solution with errors that scale as 1/ √ N independently of the dimension; here N is the number of sample points.

Nonetheless, there has been tremendous progress in integration methods for problems that lie in-between traditional quadrature (very low dimensions) and Monte Carlo (high dimensions). In intermediate dimensions (typically 5-200), 'Quasi-Monte Carlo' methods have become well established. These sample the integrand in a deterministic and structured way that ensures improved uniformity and better convergence rates. In favorable cases they can achieve error scalings of 1/N or even 1/N 2 , far outperforming traditional Monte Carlo [START_REF] Dick | Highdimensional integration: The quasi-Monte Carlo way[END_REF][START_REF] Nuyens | The construction of good lattice rules and polynomial lattice rules[END_REF][START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF][START_REF] Pierre | Randomized quasi-Monte Carlo: An introduction for practitioners[END_REF].

In this paper we show how to apply these integration techniques to perturbative expansions for quantum many-body systems. Our "Quantum Quasi-Monte Carlo" (QQMC) approach is broadly applicable. It can be formulated for both equilibrium and non-equilibrium cases and extended to various lattices and dimensions. Here we demonstrate it on a quantum dot model and show computational accelerations of several orders of magnitude compared to state-of-the-art DiagQMC [START_REF] Rosario | Quantum Monte Carlo for correlated out-of-equilibrium nanoelectronic devices[END_REF][START_REF] Bertrand | Reconstructing nonequilibrium regimes of quantum many-body systems from the analytical structure of perturbative expansions[END_REF] (Fig. 1). A crucial ingredient of QQMC is the warping of the in- tegral. This is a multi-dimensional change of variables constructed from a model function that approximates the integrand. We show that a simple model already yields remarkable results and propose directions for future optimizations. We demonstrate convergence as fast as 1/N in a high-precision benchmark against an exact Bethe Ansatz solution, to order n = 12. To illustrate the power of QQMC, we calculate the finite-bias current through a quantum dot in the Kondo regime, sweeping electrostatic gating and interaction strength as parameters. This experimentally relevant calculation was computationally unfeasible for previous techniques.

Formalism. In perturbative calculations, an observable F (U ) such as a current or susceptibility is expressed as a power series in the interaction U :

F (U ) = ∞ n=0 F n U n , (1)
where the coefficients F n are n-dimensional integrals

F n = d n u f n (u 1 , u 2 , . . . , u n ). (2)
The integrands f n (u) are time-ordered correlators expressed in terms of 2 n determinants (Wick's theorem), in both Schwinger-Keldysh [START_REF] Rosario | Quantum Monte Carlo for correlated out-of-equilibrium nanoelectronic devices[END_REF] and Matsubara formalisms [START_REF] Rossi | Determinant diagrammatic Monte Carlo algorithm in the thermodynamic limit[END_REF]. The exponential complexity of evaluating f n (u) leads us to seek fast integration methods. Here the u i specify the locations of interaction vertices in space and time. We present the formalism generally and will specialize to a concrete application later. We will perform the integral Eq. ( 2) by direct sampling using quasi-Monte Carlo. The crucial step is to warp the integral, i.e. to make a change of variables u(x) that maps the hypercube x ∈ [0, 1] n onto the u domain. The integral Eq. ( 2) becomes

F n = [0,1] n d n x f n [u(x)] ∂u ∂x , (3)
where |∂u/∂x| is the associated Jacobian.

The most important property of the warping is to make the function fn (x) = f n [u(x)] |∂u/∂x| as smooth as possible in the new variables x. If f n were positive, the perfect change of variables would make fn constant and thus trivial to integrate with a single sample. That would be tantamount to ideal sampling from the distribution f n (u) and it is as challenging as the original integration. Instead, a judicious warping must provide sufficient smoothing while remaining efficiently computable.

Mathematically, convergence theorems can only be established for fn (x) that belong to specific smooth function spaces, or whose Fourier coefficients have rapid asymptotic decay properties [START_REF] Dick | Highdimensional integration: The quasi-Monte Carlo way[END_REF][START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF]. Although we cannot prove that our warped integrands satisfy assumptions of this kind, in practice we find that the change of variables are good enough to provide excellent error scaling.

To warp the integral, we consider a positive model function p n (u), which should be viewed as an approximation of |f n |. The inverse change of variables x(u) is then defined by (for 1

≤ m ≤ n) x m (u m , u m+1 , . . . , u n ) = u m 0 du m ∞ 0 m-1 i=1 du i p n (u) ∞ 0 du m ∞ 0 m-1 i=1 du i p n (u)
(4) Here we adopt a case where u i is defined on the interval [0, ∞). Since x m (u) only depends on u m , . . . , u n , the Jacobian is |∂u/∂x| = [ du p n (u)]/p n (u) (see Appendix D). In quasi-Monte Carlo, the integral Eq. ( 3) is approximated by a sum over the first N points of a lowdiscrepancy sequence xi . This is a deterministic sequence of points with specific properties that uniformly samples the hypercube [START_REF] Dick | Highdimensional integration: The quasi-Monte Carlo way[END_REF][START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF]. We have

F n ≈ F n (N ) = C N N i=0 f n [u(x i )] p n [u(x i )] (5)
where C = du p n (u) is a constant. Here we use a Sobol' sequence [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF][START_REF] Kuo | Application of quasi-Monte Carlo methods to elliptic pdes with random diffusion coefficients: A survey of analysis and implementation[END_REF] to obtain xi . The model function p n (u) should have two key properties. First, it should approximate |f n (u)| well. Second, its form should be simple enough for the partial integrals Eq. ( 4) to be evaluated exactly and quickly. This allows the reciprocal function u(x) to be computed by first inverting the one-dimensional function x n (u n ), then inverting x n-1 (u n , u n-1 ) for fixed u n , and so on (see Appendix D).

Many classes of model functions are possible, as discussed later. This paper applies the method to impurity models, using a real-time Schwinger-Keldysh formalism, in which the u i are the times of the interaction vertices. We consider the simple form

p n (u) = n i=1 h (i) u i-1 -u i . (6)
with 0 < u n < u n-1 < . . . < u 1 < u 0 . Here u 0 = t is defined to be the measurement time and the h (i) are positive scalar functions. (They may depend on n, but we omit this index). The factored structure allows Eq. ( 4) to be inverted rapidly (see Appendix D).

Anderson Impurity. We illustrate our method on the Anderson impurity model coupled to two leads. This is the canonical model for a quantum dot with Coulomb repulsion and the associated Kondo effect. It has been realized in many nanoelectronic experiments [START_REF] Goldhaber-Gordon | Kondo effect in a single-electron transistor[END_REF][START_REF] Goldhaber-Gordon | From the Kondo regime to the mixed-valence regime in a single-electron transistor[END_REF][START_REF] Cronenwett | A tunable Kondo effect in quantum dots[END_REF][START_REF] Iftikhar | Tunable quantum criticality and superballistic transport in a "charge" Kondo circuit[END_REF]. Importantly, some quantities including the electron occupation on the dot Q can be computed analytically in the universal limit with the Bethe ansatz [START_REF] Tsvelick | Exact results in the theory of magnetic alloys[END_REF][START_REF] Okiji | Thermodynamic properties of the Anderson model (invited)[END_REF]. This provides us with a high-precision benchmark for QQMC at any perturbation order n.

We consider an infinite one-dimensional chain with the impurity at site i = 0. The non-interacting Hamiltonian is

H 0 = i,σ (γ i c † i,σ c i+1,σ + H.c.) + ε d σ c † 0σ c 0σ , where
Q n (N )/Q Bethe n -1 n = 4 n = 8 n = 12 n = 4 n = 8 n = 12
H int = U c † 0↑ c 0↑ c † 0↓ c 0↓ .
The electron tunneling between the leads and dot is γ 0 = γ -1 = γ. All other γ i = D/2, corresponding to hopping within the leads; the lead halfbandwidth D is a constant. We perform the perturbative expansion in powers of U (see Appendix E).

Benchmark. To validate the QQMC method, we consider the special case solved by the Bethe Ansatz. For this, we set temperature T = 0, capacitive gate ε d = 0, and half-bandwidth D → +∞ such that Γ = 4γ 2 /D = 1 is the unit of energy. The measurement time t = 30/Γ is sufficiently long that the system reaches steady-state. We compute the expansion of the occupation number

Q(U ) = c † 0↑ c 0↑ + c † 0↓ c 0↓ .
The system is particle-hole symmetric for ε d = -U/2 so the non-interacting case is Q 0 = 1. For higher-order Q n , particle-hole symmetry is broken, but the expansion stays in the symmetric regime (U + 2ε d ) √ U Γ [START_REF] Tsvelick | Exact results in the theory of magnetic alloys[END_REF]. Figure 2 shows the relative error between Q n (N ) using QQMC and the exact result Q Bethe n (see Appendix F), as a function of the number of integrand evaluations N . Following an initial transient, we enter an asymptotic regime in which there is rapid convergence: for n = 4 this is consistent with pure 1/N while for n = 8, 12 it is 1/N δ with δ 0.9, 0.8. These calculations used the product model function Eq. ( 6) with a single exponential h (i) (v i ) = exp(-v i /τ ), where τ = 0.95. The same set-up was used in Fig. 1. The level of precision that we obtained revealed limitations in the conventional evaluation of the non-interacting Green functions, which warranted special consideration (see Appendix E).

It is expected that the convergence rate gradually slows as n increases. First, the quality of the warping decreases as the disparity between the increasingly-severe requirements of convergence theory and the behavior of our integrands grows. This can be mitigated by constructing more expressive model functions, which we discuss below. Second, for larger n the integrands generally become more oscillatory. The model functions Eq. ( 4) were not designed to handle cases with massive cancellation, and this may become a limiting factor. We will see this effect below for calculations with ε d /U > 0.5, although in practice enough orders can be computed accurately to obtain the desired physical results (see Appendix E).

In Quasi-Monte Carlo methods, a standard technique to estimate errors is to perform computations using Eq. ( 5) with several 'randomized' low-discrepancy sequences [START_REF] Dick | Highdimensional integration: The quasi-Monte Carlo way[END_REF][START_REF] Nuyens | The construction of good lattice rules and polynomial lattice rules[END_REF][START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF][START_REF] Pierre | Randomized quasi-Monte Carlo: An introduction for practitioners[END_REF] and we use this method below (see Appendix G).

Having made these technical points, let us reiterate the lessons of Fig. 1 and Fig. 2: (i) QQMC provides a dramatic speed-up with better asymptotic error scaling than DiagQMC; (ii) the speed-up persists up to at least order n = 12, which is what is needed for practical applications.

Coulomb Diamond. We now apply QQMC to solve a topical physics problem. We explore the current-voltage characteristic I(V ) across the quantum dot for finite bias and varying U . Since quantum dots are considered promising platforms for building qubit systems, it is of primary importance to understand how many-body effects influence their properties, especially the phase coherence.

Quantum dots can be in three different experimentally accessible regimes [START_REF] Sander | Individual single-wall carbon nanotubes as quantum wires[END_REF][START_REF] Nygård | Kondo physics in carbon nanotubes[END_REF][START_REF] Liang | Fabry -Perot interference in a nanotube electron waveguide[END_REF][START_REF] Roch | Quantum phase transition in a single-molecule quantum dot[END_REF]: Fabry-Pérot (small U ), Kondo (intermediate U ) and Coulomb blockade (large U ). The Fabry-Pérot and Coulomb blockade limits are well described by, respectively, non-interacting and semi-classical theories; the out-of-equilibrium Kondo regime is more challenging. Two controlled approaches have recently appeared, but both are too slow for some applications: the Schwinger-Keldysh DiagQMC used in Figs. 1 and 4 [10, 11] and the real-time inchworm algorithm [START_REF] Cohen | Taming the dynamical sign problem in real-time evolution of quantum many-body problems[END_REF][START_REF] Ridley | Numerically exact full counting statistics of the energy current in the Kondo regime[END_REF][START_REF] Krivenko | Dynamics of Kondo voltage splitting after a quantum quench[END_REF]. QQMC provides the speed and precision to allow large parameter sweeps, which is mandatory to make good contact with experiments. In [START_REF] Bertrand | Reconstructing nonequilibrium regimes of quantum many-body systems from the analytical structure of perturbative expansions[END_REF], some of us studied the Kondo ridge close to ε d = -U/2. QQMC allows us to present results scanning the entire (U, ε d ) phase diagram, including slowly converging regions with even numbers of electrons or near the degeneracy points.

Figure 3 (inset) shows a cartoon of the differential conductance for varying (ε d , V ) as predicted by Coulomb blockade theory [START_REF] Beenakker | Theory of Coulomb-blockade oscillations in the conductance of a quantum dot[END_REF] and seen experimentally at low temperatures and large U [START_REF] Hofheinz | Capacitance enhancement in Coulomb blockade tunnel barriers[END_REF]. At small bias, the Coulomb blockade forbids current flows except at two special points: ε d = 0, where the dot energies for Q = 0 and Q = 1 electrons are degenerate, and ε d = -U (likewise for Q = 1, 2). At intermediate U , the Kondo effect changes this picture drastically: the zero-bias Kondo resonance forms in the 'forbidden' region of odd Q and enables current flow.

Figure 3 shows the current I versus gate voltage ε d for V = U/7 and temperature T = 0. We choose Sweeping the interaction U/Γ shows several regimes. For U/Γ = 1.75, 3.5 a current plateau emerges in the local moment regime (Q = 1) due to Kondo resonance formation. The current develops new local maxima seen for U/Γ = 5.25, 6.00. These grow toward the Coulomb blockade limit at larger U (black lines); at the same time, current around ε d /U = -0.5 reaches a maximum and decreases. This is a competition between resonance formation and narrowing. At small U , the Kondo temperature T K is much larger than the bias V and we are in the linear response regime. In this regime near ε d = -U/2 we approach perfect transmission I = V /π; see Ref. [START_REF] Bertrand | Reconstructing nonequilibrium regimes of quantum many-body systems from the analytical structure of perturbative expansions[END_REF]. At larger U 4Γ, T K decreases exponentially with U and become smaller than V , leaving the linear response regime. Throughout, as U increases, the already-small current in the side regions (Q = 0, 2) is increasingly suppressed.

Q = 2 Q = 1 Q = 0 V = U/7 I/Γ U/Γ = +∞ U/Γ = 1.75 U/Γ = 3.5 U/Γ = 5.25 U/Γ = 6.0 -U 0 ε d -U 0 V =U
Model Function. Let us reexamine the importance of integral warping and model functions. Figure 4 shows the convergence of Q 8 (N ) using different integration methods; parameters are identical to Fig. 2. When the integral is evaluated using Sobol' points without warping ('Sobol' only') the convergence is poor, showing that naively applying low-discrepancy sequences provides little benefit for these integrands. Next, contrast regular DiagQMC with the warped integrand using pseudo-random numbers.

As expected for pure Monte Carlo approaches, both show 1/ √ N convergence. Nonetheless, sampling the warped integrand still converges faster than DiagQMC, despite Comparison of convergence of Q8 for different methods of integration: evaluating unwarped integrand with a Sobol' sequence (cyan), DiagQMC (red), warped integral sampled with Mersenne Twister pseudo-random numbers (green) or Sobol' sequence (blue). For the warped cases, we used Eq. ( 6) with h (i) (vi) = exp(-vi/τ ), τ = 0.95. After an initial warping with exponential functions τ = 1.1, we can apply an additional warping obtained by projection (orange); see Appendix I. For visibility, data (except Sobol' and DiagQMC) has been smoothed in the same way as in Fig. 2.

the fact that the latter uses importance sampling via the Metropolis algorithm. As anticipated, QQMC using Sobol' points and the model function Eq. ( 6) based on exponential h (i) converges even more rapidly.

How can the model function Eq. ( 6) with simple h (i) provide such dramatic convergence improvements? Our integrands describe physical correlators that are highly structured and have decaying exponential or power-law tails; see Appendix I and Refs. [START_REF] Rosario | Quantum Monte Carlo for correlated out-of-equilibrium nanoelectronic devices[END_REF][START_REF] Bertrand | Quantum Monte Carlo algorithm for out-of-equilibrium Green's functions at long times[END_REF]. The tail contributions become ever more important as the dimension increases. The model function properly describes the long-time asymptotics (see Appendix I). We also emphasize the importance of a well-chosen coordinate system in the model function: the differences of closest times

v i = u i-1 -u i used to parametrize the h (i) .
Optimization of the model function should allow further performance gains, particularly at higher orders n. One possibility is to better adapt the functions h (i) to f n . To illustrate this, we apply a second warping constructed by sampling points from the first warping. These samples are projected along the dimensions of v space and smoothed; see Appendix I for details. As shown in Fig. 4, this optimization reduces the error by a factor of 2. More importantly, it automatically gives robust convergence without the need to manually optimize the τ parameter.

Finally, other families of model functions exist beyond Eq. ( 6), that provide versatile and expressive approximations while still allowing for fast inversion of Eq. ( 4). One such family is Matrix Product States (MPS) or functional tensor-trains [START_REF] Schollwöck | The density-matrix renormalization group in the age of matrix product states[END_REF][START_REF] Glasser | Expressive power of tensor-network factorizations for probabilistic modeling, with applications from hidden Markov models to quantum machine learning[END_REF], of which Eq. ( 6) is just the simplest case:

p n (u) = h (1) a (v 1 )h (2) ab (v 2 ) • • • h (n-1) cd (v n-1 )h (n) d (v n ). (7)
Here

h (i)
ab are matrices and repeated indices are summed. Another promising family is

p n (u) = n-1 i=1 h(i) (v i+1 , v i ).

Conclusion.

We have shown how to use sampling techniques based on low-discrepancy sequences to compute high orders of many-body perturbation theory. Although we cannot show that the integrands obey the assumptions of formal Quasi-Monte Carlo convergence theory, practical scaling as fast as 1/N is still achievable. This success was possible due to the warping of the integral based on a model function. Using benchmarks on exactly solvable quantities in the Anderson impurity model, we unambiguously validated the convergence of this 'Quantum Quasi-Monte Carlo' (QQMC) method at high-precision. This calculation was about ∼ 10 4 times faster than the DiagQMC equivalent.

We can apply the techniques established here to models with interesting strongly correlated physics in all dimensions, for equilibrium and especially non-equilibrium situations. For continuum models, the integrands are smooth and warping should be particularly simple. For lattice models, the discrete summation may degrade convergence, although this may be addressed with sufficiently good model functions. QQMC can also be applied to other diagrammatic expansions, e.g. in hybridization [START_REF] Cohen | Taming the dynamical sign problem in real-time evolution of quantum many-body problems[END_REF]. Constructing more expressive model functions should further increase speed and accuracy and is an ideal application for recent machine learning techniques in quantum systems.

Finally, we have shown in our calculations that the simple model function Eq. ( 6) captures the behavior of perturbation theory integrands in asymptotic large-coordinate regions. This is not accidental, but reveals a simplifying structure of the correlation functions arising from Wick's theorem that was not previously appreciated in diagrammatic numerical simulations. For the real-time Schwinger-Keldysh calculations, the contour index means that the MPS structure Eq. ( 7) is the natural approximation for generic many-body systems. It can be used as a starting point to efficiently compute and integrate these functions, even beyond the Monte Carlo or QQMC sampling discussed here.

The integration technique for perturbation theory integrals we presented in the main manuscript is based on so-called 'quasi-Monte Carlo' methods, which use sampling with low-discrepancy sequences. Despite the naming, such quasi-random sequences are highly structured and deterministic, unlike the pseudo-random number sequences used in conventional Monte Carlo sampling (see Fig. 5). Here we briefly summarize the history of this field of mathematics as well as key features of these sequences.

The roots of low-discrepancy and quasi-random sequences lie in the ancient field of Diophantine Approximation -the theory of approximating sets of real numbers by rational ones, especially in its modern form pioneered by Roth [START_REF] Roth | On irregularities of distribution[END_REF] -and in Weyl's early efforts to characterize uniformity of distributions [START_REF] Weyl | Über die Gleichverteilung von Zahlen mod. Eins[END_REF]. These methods were applied to practical multidimensional integration in the late 1950's and 1960's, with contributions from, among many others, Korobov [START_REF] Korobov | Approximate calculation of repeated integrals by number-theoretical methods[END_REF] and Sobol' [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF] in the USSR and Hammersley [START_REF] Hammersley | Monte Carlo methods for solving multivariable problems[END_REF] (building on work by van der Corput [START_REF] Van Der Corput | Verteilungsfunktionen[END_REF]), Halton [START_REF] Halton | On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals[END_REF] and Haselgrove [START_REF] Haselgrove | A method for numerical integration[END_REF] in the UK. There was even an unsuccessful attempt to use them as a component of Markov Chain Monte Carlo simulation for multi-point physical integrals in 1951 [START_REF] Richtmyer | The evaluation of definite integrals, and a quasi-Monte-Carlo method based on the properties of alge-braic numbers[END_REF]. Niederreiter [START_REF] Niederreiter | Quasi-Monte Carlo methods and pseudo-random numbers[END_REF] surveys much of the early work in this area. Since that time, quasi-Monte Carlo methods have become standard in many fields including, for example, computational chemistry [START_REF] Conroy | Molecular Schrödinger equation. VIII. A new method for the evaluation of multidimensional integrals[END_REF]. Indeed, convergence rate improvements analogous to, but smaller than, ours are seen in molecular excitation computations [START_REF] Berblinger | Monte Carlo integration with quasi-random numbers: some experience[END_REF]. Most of the early applications had relatively low dimensions 10. However, the field was rejuvenated in the late 1990's and early 2000's by the discovery that quasi-Monte Carlo works unexpectedly well for certain high-dimensional financial integrals stemming from discretized stochastic partial differential equations [START_REF] Spassimir | Faster Valuation of Financial Derivatives[END_REF]. This was followed by the application of Reproducing Kernel Hilbert Space theory to explain this success (see e.g. Ref. [START_REF] Wang | Why are highdimensional finance problems often of low effective dimension?[END_REF]), and the advent of 'fast component-by-component' construction techniques Quasi-Random (Sobol') Pseudo-Random (MT) to create optimized quasi-random generators for these problems [START_REF] Nuyens | Fast algorithms for component-by-component construction of rank-1 lattice rules in shift-invariant reproducing kernel Hilbert spaces[END_REF]. A survey of this more recent work appears in Ref. [START_REF] Dick | Highdimensional integration: The quasi-Monte Carlo way[END_REF].

The main strength of quasi-Monte Carlo relative to Markov Chain Monte Carlo is its greatly accelerated convergence rate in applications with well-behaved integrands. Its main disadvantage is its reduced flexibility. Quasi-Monte Carlo methods require highly symmetric integration domains -usually hypercubes -and careful preparation of the integrand. In our approach we implement this through a change of coordinates and warping. However, perhaps most importantly, they are intrinsically non-adaptive during a calculation: changing the distribution during the main sampling run destroys its carefully constructed uniformity properties, which typically reduces the convergence rate to Monte Carlo scaling O(1/ √ N ). To handle generic integrands whose properties are unknown in advance, separate runs and analysis are needed before the main sampling run.

Finally, we note that the theory of quasi-Monte Carlo typically gives worst-case error bounds, but not averagecase ones. This is because no stochastic averaging is done. An exception is when we perform averaging of random shifts. These worst-case bound are theoretically tight but often rather pessimistic compared to observed behavior. For pure Monte Carlo methods, the opposite situation is true and we typically have average-case but not worst-case bounds.

Appendix B: Summary of the QQMC algorithm Below, we summarize the QQMC algorithm. More details on the building blocks are given in other sections of this supplementary material. Note that unlike usual diagrammatic Monte Carlo techniques, the coefficients F n for different orders n are calculated separately.

Preprocessing steps.

• Calculate non-interacting Green functions in real time. These functions form the basic elements from which the integrands are calculated. They can be obtained through Fourier transform of the real frequency Green functions which can be obtained analytically or numerically.

• Calculate Warping. Once a model function p n (u) has been selected and possibly adapted to the integrand (see projection method below), construct the mapping x(u) by computing the partial cumulative functions of the model function p n (u).

• For each component of the mapping x(u), invert it to obtain the inverse mapping u(x). This can be done, for example, be interpolating the original mapping x(u) on a sufficiently fine linear mesh and inverting.

• Initialize r = 0.

Main computing loop. For each sample point i ∈ 1 . . . N ,

• Generate the low-discrepancy quasi-random vector x i . These vectors span uniformly the hypercube [0, 1] n .

• Calculate the corresponding point u i = u(x i ) in the original integration space.

• Calculate the corresponding value of the integrand f (i) = f n (u i ) and the value of the model function

p (i) = p n (u i ).
• Update the result r → r + f (i) /p (i)

Returns. The final estimate of the observable is F n ≈ r/N .

Appendix C: Benchmarking at higher orders

The main manuscript focuses on data for perturbation theory orders n 12, where QQMC enables calculations We used Eq. ( 6) with h (i) (vi) = exp(-vi/τ ) with τ = 0.95 for n ≤ 15 and τ = 0.90 for n ≥ 16. The parameters are the same as in Fig. 1.

with an unprecedented precision of up to six or seven digits. In this section, we show additional data for a smaller number of samples N . This enables us to calculate much larger orders up to n = 22, well beyond previous diagrammatic quantum Monte Carlo techniques. We recall that the a single evaluation of the integrand f n (u) has an exponential complexity ∼ 2 n , which has limited previous calculations to order n ≤ 15 with less than two digits precisions. While N = 10 5 is not yet deep into the asymptotic regime, we find that we could reach two to three digits of accuracy.

The data are presented in Fig. 6 respectively for the relative error (top panel) and absolute error (lower panel). Although the error deteriorates with the order n, the speed up provided by QQMC more than compensates for the imperfection of the model function. We find that the absolute error actually decreases with n as the Q n gets smaller at large n.

Appendix D: Model Function Properties

Here we expand on the properties of the change of variables x(u) arising from the model p n (u). These were defined in the main text as:

x m (u m , . . . , u n ) = u m 0 du m ∞ 0 m-1 i=1 du i p n (u) ∞ 0 du m ∞ 0 m-1 i=1 du i p n (u) (4 )
To understand the structure of Eq. ( 4) and how it can be useful when sampling from p n (u), let us consider the explicit transformation for small orders n.

For n = 1, Eq. ( 4) is simply the normalized cumulative function x 1 (u 1 ) = u1 0 dū 1 p 1 (ū 1 )/ ∞ 0 dū 1 p 1 (ū 1 ) which upon differentiation gives,

dx 1 = p 1 (u 1 )du 1 ∞ 0 dū 1 p 1 (ū 1 )
.

(D1)

This means that uniformly sampling x 1 leads to the sampling of p(u 1 ). In practice, one need to invert the cumulative distribution x 1 (u 1 ) which can be done through interpolation techniques.

Next, let us consider the procedure for n = 3 using a model function p 3 (u 1 , u 2 , u 3 ). The reverse coordinate transform x(u) is

x 1 (u 1 , u 2 , u 3 ) = u1 0 dū 1 p 3 (ū 1 , u 2 , u 3 ) ∞ 0 dū 1 p 3 (ū 1 , u 2 , u 3 ) x 2 (u 2 , u 3 ) = u2 0 dū 2 ∞ 0 dū 1 p 3 (ū 1 , ū2 , u 3 ) ∞ 0 dū 2 ∞ 0 dū 1 p 3 (ū 1 , ū2 , u 3 ) (D2) x 3 (u 3 ) = u3 0 dū 3 ∞ 0 dū 2 ∞ 0 dū 1 p 3 (ū 1 , ū2 , ū3 ) ∞ 0 dū 3 ∞ 0 dū 2 ∞ 0 dū 1 p 3 (ū 1 , ū2 , ū3 )
The consecutive coordinate integration, gives the coordinate transformation a special structure: x 3 (u 3 ) does not depend on u 1 , u 2 and x 2 (u 3 , u 2 ) does not depend on u 1 . This means that the Jacobian matrix for the reverse coordinate transformation x(u) has an upper-triangular form

∂x ∂u =   (∂x 1 /∂u 1 ) (∂x 1 /∂u 2 ) (∂x 1 /∂u 3 ) 0 (∂x 2 /∂u 2 ) (∂x 2 /∂u 3 ) 0 0 (∂x 3 /∂u 3 )   (D3)
so that the Jacobian determinant is simply

∂x ∂u = ∂x 1 ∂u 1 • ∂x 2 ∂u 2 • ∂x 3 ∂u 3 . ( D4
)
Differentiating Eq. (D2), cancelling common factors and using |∂u/∂x| = 1/ |∂x/∂u| gives

∂u ∂x = ∞ 0 du 3 ∞ 0 du 2 ∞ 0 du 1 p 3 (u 1 , u 2 , u 3 ) p 3 (u 1 , u 2 , u 3 ) . (D5)
The same procedure straightforwardly generalizes to an arbitrary number of dimensions n. This reproduces the result quoted in main text:

∂u ∂x = ∞ 0 n i=1 du i p n (u) p n (u) . ( D6
)
In practice, we uniformly sample the hypercube [0, 1] 3 using Sobol' sequence to obtain (x 1 , x 2 , x 3 ). From x 3 one obtains u 3 by inverting the one dimensional equation x 3 (u 3 ). With the obtained value of u 3 , the equation x 2 (u 2 , u 3 ) becomes a one dimensional function of u 2 which can be inverted. Last with the obtained (u 2 , u 3 ), one can invert x 1 (u 1 , u 2 , u 3 ) to obtain u 1 .

Product Model. While the coordinate transform described above is very general, it is only useful if the multiple integrals in Eq. ( 4) can be performed efficiently. Otherwise, it is as or more costly than the actual integral of the perturbation series coefficient that we wish to compute.

The product model function

p n (u) = n i=1 h (i) u i-1 -u i , (6 )
is particularly efficient. It is simpler to view this as a composition of two transforms. First, we change variables v i = u i-1u i , which has Jacobian |∂u/∂v| = 1. Second, in the v i variables, the coordinate transform Eq. ( 4) separates entirely, so that each x m only depends on a single variable v m :

x m (v m ) = vm 0 dv m h (m) (v m ) ∞ 0 dv m h (m) (v m ) (D7)
These one-dimensional integrals can be integrated quickly and precisely analytically or using quadrature algorithms.

We then invert x m (v m ) numerically to find the coordinate transform v m (x m ). In practice, it is possible to completely pre-compute these integrals on a fine mesh for fast evaluation during calculation.

We note that the product form Eq. ( 6) is known for importance sampling in Monte Carlo applications, e.g. as part of the VEGAS algorithm [START_REF] Lepage | A new algorithm for adaptive multidimensional integration[END_REF][START_REF] Lepage | VEGAS -an adaptive multi-dimensional integration program[END_REF]. The choice of coordinate system v i as compared to u i affects the quality of model function and is an important physical consideration.

MPS Model. The model function Eq. ( 7) can also be efficiently computed using the above algorithm and standard MPS techniques [START_REF] Schollwöck | The density-matrix renormalization group in the age of matrix product states[END_REF]. Unlike the product model, the integrations for different v i have to be performed in sequence, with a matrix-vector multiplication at each step.

Appendix E: Explicit expressions of the integrands

In this appendix, we will describe the Anderson impurity model as well as the perturbation expansion formalism of our calculations in more detail.

Model

The Hamiltonian of the Anderson impurity model is H = H 0 + H int θ(t), with the non-interacting term:

H 0 = i,σ γ i c † i,σ c i+1,σ + H.c. + E d σ c † 0σ c 0σ . (E1)
The hopping parameters are all γ i = D/2, except at the impurity (i = 0) where γ 0 = γ -1 = γ. The interaction term is

H int = U (c † 0↑ c 0↑ -α)(c † 0↓ c 0↓ -α), ( E2
)
where α is a quadratic shift to the perturbation. This shift means that the U expansion is performed about a different starting point, and is commonly used to improve perturbation series convergence [START_REF] Rosario | Quantum Monte Carlo for correlated out-of-equilibrium nanoelectronic devices[END_REF][START_REF] Wu | Controlling Feynman diagrammatic expansions: Physical nature of the pseudogap in the two-dimensional Hubbard model[END_REF][START_REF] Rubtsov | Continuous-time quantum Monte Carlo method for fermions: Beyond auxiliary field framework[END_REF]. Note that the energy of a single electron localized on the impurity is ε d = E d -αU . A symmetric voltage bias V is applied between the two leads. As is standard in the Keldysh formalism, H int is turned on at time t = 0 and observables are computed after a large time t when the stationary regime has been reached. By integrating out the leads, their effect on the dot is represented by a retarded hybridization function ∆(ω). The non-interacting dot retarded Green function is g R (ω) ≡ 1/(ω -E d -∆(ω)). By symmetry, all Green functions are the same for spin up and down.

The density of states of the leads is semi-circular with half-bandwidth D. An important parameter of the noninteracting model is the tunneling rate from the impurity to the leads at the (equilibrium) Fermi level Γ = 4γ 2 /D. In terms of D and Γ, the hybridization function is

∆(ω) = Γ D ×      ω + √ ω 2 -D 2 for ω < -D ω -i √ D 2 -ω 2 for -D < ω < D ω - √ ω 2 -D 2 for ω > D (E3)
The Bethe ansatz provides results only in the universal regime, where D → +∞ with Γ fixed. In this limit, the density of states becomes independent of energy (flat band). The hybridization function is simply ∆(ω) = -iΓ. When comparing results to the Bethe ansatz, we we will always work in this regime.

Expressions for the integrands of series expansions

To obtain the number of electrons Q or the current I, we compute a perturbation series in U for the equaltime lesser Green function

G < 0i (t, t) = -i c † 0↑ (t)c i↑ (t)
, where the creation and annihilation operators are in the Heisenberg picture. We take i = 0 (on-site) to obtain Q, and i = 1 (dot-lead) to obtain I [START_REF] Meir | Landauer formula for the current through an interacting electron region[END_REF].

For α = 0, this series can be written compactly using the "Wick determinant" notation introduced in Ref. [START_REF] Bertrand | Quantum Monte Carlo algorithm for out-of-equilibrium Green's functions at long times[END_REF]:

G < 0i (t, t) = n≥0 i n U n n! t 0 n k=1 du k × a1,...,an (-1) a k (0, t, 0), U 1 , . . . , U n (i, t, 1), U 1 , . . . , U n U 1 , . . . , U n U 1 , . . . , U n , (E4)
where a k ∈ {0, 1} are Keldysh indices, and U k = (0, u k , a k ) represents a point on the Keldysh contour composed of a site index (here 0 for the impurity), a time u k and a Keldysh index a k . The Wick determinant . . . is defined, for A 1 , . . . , A m and B 1 , . . . , B m any set of points on the Keldysh contour, in the case α = 0, the Wick determinant

A 1 , . . . , A m B 1 , . . . , B m = g(A 1 , B 1 ) . . . g(A 1 , B m ) . . . . . . . . . g(A m , B 1 ) . . . g(A m , B m ) , (E5)
where

g[(x, u, a), (x , u , a )] = g T xx (t, t ) g < xx (t, t ) g > xx (t, t ) g T xx (t, t ) aa , ( E6
) is the non-interacting one-particle Keldysh Green function. Here g T , g T , g < and g > are respectively the time-ordered, anti-time-ordered, lesser and greater Green functions. In Eq. (E4), the determinant on the left is from spin up operators, while the one on the right is from spin down operators. Nevertheless, by spin symmetry their elements share the same Green functions g.

The case α = 0 is similar, but the diagonal terms of the Wick determinants in Eq. (E4) must be shifted by -iα, except the one connecting to the measurement point (involving the Green function at time t) [START_REF] Rosario | Quantum Monte Carlo for correlated out-of-equilibrium nanoelectronic devices[END_REF].

Provided that the non-interacting Green function g is known as a function of time, Eq. (E4) explicitly defines the integrand that we refer to in the main text of this article. We compute the time domain Green functions by Fourier transform of the Green functions g < (ω) and g > (ω). These can be derived in the Schwinger-Keldysh formalism [START_REF] Meir | Landauer formula for the current through an interacting electron region[END_REF]:

g < 00 (ω) =2i|g R (ω)| 2 Im[∆(ω)] (E7) × n F ω - V 2 + n F ω + V 2 g > 00 (ω) =g < 00 (ω) -4i|g R (ω)| 2 Im[∆(ω)] (E8) γg < 01 (ω) = ∆(ω) 2 g < 00 (ω) -2i n F ω + V 2 Im[∆(ω)]g R (ω) * (E9)
Here n F (ω) is the Fermi function. Note that the function g > 0i is not used in Eq. (E4).

3. Precision calculation of g < (t) and g > (t)

Care has to be taken when performing the Fourier transform to obtain g < and g > in the time domain. As integration methods get increasingly precise, the accuracy of the integrand becomes more important. In order to provide benchmarks with relative error of ∼ 10 -6 at order n ∼ 10, and to rule out any bias due to inexact integrands, we need to refine the calculation of g(t). In particular, using a Fast Fourier Transform (FFT) algorithm produces an error which decreases too slowly with the number of samples for functions with sharp features or power law tails, such as the ones we encounter here.

One approach to high precision is to compute the Fourier transform using adaptive quadrature methods. This is precise enough if a system has a finite bandwidth and the integrand is proportional to the density of states of the leads and therefore has bounded support. In general, however, the integrand decays slowly and oscillates at high frequencies, which renders direct integration methods inaccurate. Alternatively, when the tails are dominated by simple poles, it is possible to separate them out analytically and perform the finite remainder using a FFT. In our case, however, the tails are dominated by the Fermi functions, which have an essential singularities at |ω| = ∞, and we must resort to other methods.

To circumvent this problem, we deform the integration path in the ω complex plane to find a more favorable integrand, and apply an adaptive quadrature method. We show that a path can be found for a generic class of problems which improves the decay rate and eliminates oscillations near infinity. We consider the general case of finite temperature, and denote the inverse temperature by β. We will work out the case for g < , but g > can be treated equally by first applying a change of variable ω → -ω. Specifically, at time t, the Fourier Transforms we are interested in can always be decomposed in a sum of integrals of the form:

dω ζ(ω )n F (ω )e -iω t (E10)
with ω = ω ± V /2. The function ζ depends on g R and ∆. Its exact form does not affect the choice of a new path, as long as it has no singularity at |ω | = ∞. We will further assume that one can bound the complex singularities of ζ and n F inside a vertical band ω -< Re[ω ] < ω + .

For t = 0, we compute the integral Eq. (E10) along a new path parametrized by x and defined as:

ω (x) =      ω --(x -ω -)it, for x < ω -, x, for ω -< x < ω + , ω + + (x -ω + )(β -it), for x > ω + .
(E11) The new path is made of three pieces, joined together at ω = ω ± . The central one is simply a segment of the real axis, left unchanged to prevent crossing singularities. The other two are straight lines at an angle with the real axis, which have been chosen so that the integrand becomes asymptotically proportional to a decaying, oscillation-free, exponential. The points ω ± can be moved away from one another, in particular to avoid the x < ω -piece being too close to singularities. For β = +∞, the integrand is zero on the half-plane Re[ω ] > 0, hence the x > ω + piece of the path can be ignored. In the case t = 0, the integrand is simply g < 0i (ω), which is free of oscillations, and deforming the integration path would make some appear. Hence the integration path is left untouched in this case. This technique is easily generalized to more complex impurity systems. However, it relies on an analytical continuation and knowledge of singularities. This may not be easily accessible for numerically computed Green functions.

Appendix F: Bethe Ansatz Comparison

Here we briefly discuss the Bethe ansatz solution for the Anderson impurity model and how specifically we extract the coefficients Q Bethe n from the general solution (see [START_REF] Tsvelick | Exact results in the theory of magnetic alloys[END_REF][START_REF] Okiji | Thermodynamic properties of the Anderson model (invited)[END_REF] and references therein). We are interested in the case were ε d = 0 and we then perform a perturbative expansion in U . This always corresponds to the so-called

symmetric limit U/2 + ε d √ U Γ.
We use the results of [START_REF] P B Wiegmann | Exact solution of the Anderson model: I[END_REF]. In the symmetric limit, the occupation number Q on the quantum dot is given by a series Here:

Q = 1 - ∞ n=0 √ 2 π (-1) n (2n + 1) G (+) [iπ(2n + 1)] • ∞ -∞ dk ∆(k)e -π(2n+1)[g(k)-Λ] . (F1) n Occupation Number Q Bethe n 0 1 1.000000000000000000000000 1 -1 π -0.3183098861837906715377675 2 1 π 2 0.1013211836423377714438795 3 π 2 -9 3π 3 0.009348692094998422959323018 4 90-11π
∆(k) = Γ π • 1 (k -ε d ) 2 + Γ 2 , (F2) g(k) = (k -ε d -U/2) 2 2U Γ , (F3)
G (+) (ω) = √ 2π Γ 1 2 -1 2π iω -iω + 0 2πe -iω/2π (F4)
and Λ is the energy cutoff. All other quantities are the same as in Appendix E. The cutoff Λ is given implicitly by the series: ) [iπ(2n + 1)]e πΛ(2n+1) (2n + 1) 3/2 (F5)

U/2 + ε d √ U Γ = 1 2 • 2 π ∞ n=0 (-1) n G (+
We emphasize that Eq. (F5) differs from the expression in [START_REF] P B Wiegmann | Exact solution of the Anderson model: I[END_REF] by a factor of 1/2. The constraint on Λ arises from imposing number conservation on the distribution function, which we take to be of the form used by [START_REF] Okiji | Thermodynamic properties of the Anderson model (invited)[END_REF].

We now set ε d = 0, so that our expansion parameter is U/Γ 1. To extract the coefficients Q n we use the following procedure. First, write Eq. (F5) as a polynomial in x = e Λπ up to order N cutoff = 30. Second, perform polynomial inversion to find x as a function of U/Γ, repeatedly using the smallness of U/Γ. Third, evaluate the k integral in Eq. (F1) analytically in an asymptotic expansion in U/Γ up to order N cutoff , (see e.g. [START_REF] Horvatić | Equivalence of the perturbative and Bethe-Ansatz solution of the symmetric Anderson Hamiltonian[END_REF]). Fourth, substitute the expansion of x, collecting terms of the same order in U/Γ. The final result is a power series expansion in U/Γ, with analytic coefficients. These are evaluated numerically with high precision arithmetic and shown in Table I. For orders n = 0 -5, the analytic expressions match the results calculated explicitly in perturbation theory in [START_REF] Horvatić | Equivalence of the perturbative and Bethe-Ansatz solution of the symmetric Anderson Hamiltonian[END_REF]. Absolute error of occupation number Q6(N ), with curves matching three cases of Fig. 4. Here we perform calculations with K = 25 shifted sequences. The pale curves show the convergence of the average over these runs to the Bethe ansatz result. The solid curves show the corresponding uncertainty estimates computed from the standard error of the mean over the shifted sequence values. We see that this procedure gives good estimates of the error. a system away from particle-hole symmetry (E d ≈ 0.9Γ, α ≈ 0.27, which in Fig. 3 corresponds to the point ε d /U ≈ -0.12 and U = 6Γ). We observe a scaling slightly below 1/N , which deteriorates with increasing order. Also, at large orders, the final scaling regime is reached at a later N .

Note that the data of Figs. 2, 4 and 8 have been smoothed for visibility: for N > 100, we show the maximum of the error in a moving window around N of fixed size in log-space (5% of N ). This smoothing generates an upper bound of the error.

For each integration, digitally shifted [START_REF] Kuo | Application of quasi-Monte Carlo methods to elliptic pdes with random diffusion coefficients: A survey of analysis and implementation[END_REF] generators are used to produce 10 different Sobol' sequences from which we take the average to obtain the final result I n . The error is estimated by taking the standard deviation of the 10 results and dividing by √ 10. An example of a series computed with its estimated error is shown in Fig. 9. The error (black dots) on the coefficients (colored dots) is low enough so that only the truncation of the series limits the resummation accuracy.

The convergence radii of the series are about 3-4 Γ, so the series are resummed to obtain answers at larger interaction U . We use the conformal transform technique of Ref. [START_REF] Bertrand | Reconstructing nonequilibrium regimes of quantum many-body systems from the analytical structure of perturbative expansions[END_REF] with the so-called parabola transform W =tan 2 ( U/p). Here p is a real negative parameter optimized for each series. Depending on their analytical structure, the series I(U ) or its inverse 1/I(U ) is resummed, whichever gives smaller error. The integration error is propagated through the resummation process and added to the truncation error, the latter being estimated from the convergence radius of the resummed series. Unlike in Ref. [START_REF] Bertrand | Reconstructing nonequilibrium regimes of quantum many-body systems from the analytical structure of perturbative expansions[END_REF], no Bayesian inference is used.

At large ε d , the integrands become more difficult to This data corresponds to a typical point of Fig. 3, away from particle-hole symmetry (ε d /U ≈ -0.12, α ≈ 0.27, U = 6Γ). At low order (blue and orange lines) the relative error scales as 1/N 0.9 (black plain line). At larger orders (n = 7 and 10, green and red lines), the convergence slows down and scales only as 1/N 0.7 (dashed black line). Finally at order n = 10, the final scaling starts at a larger number of function evaluations N than for lower orders. For visibility, the data has been smoothed as described in Appendix H (see main text). For all orders calculated, up to n = 10, we obtained at least two digits of accuracy.

integrate, but at the same time the summation of the series requires less orders for the same precision. When ε d /U > 0.5 (or ε d /U < -1.5 by symmetry) we only computed and summed the series up to order n = 5. Finally, it is worth noting that small bias V setups reach a stationary regime after a longer time. Hence performing the integration away from the linear response regime, at large V , is actually less expensive and less prone to sign problem. Physically, the voltage reduces the coherence of the system, which explains the reduction of almost exact cancellations in the numerical integration.

Exponential form for h (i)

Our second change of variable v → w will be based on a model function with a simple analytic form for h (i) , designed to correctly describe the asymptotics of the integrand. This asymptotic region becomes increasingly important at large order n, especially in the long time limit t → ∞. We choose a simple form, independent of i h

(i) exp (v) = e -v/τ , ∀i. (I2)
This choice is motivated by a direct study of the integrand. In Fig. 10, we plot the absolute value of the integrand

|f 5 (v 1 , v 2 , v 3 , v 4 , v 5 )| of Q 5
with the parameters of our benchmark along various directions of the 5 dimensional v-space. The different colors correspond to different directions: red corresponds to (v, δ, δ, δ, δ) where δ is fixed to δ = 0.5 and v is varied. The 5 different curves correspond to different permutations of the v with respect to the δ. Green curves correspond to the 10 different permutations of (v, v, δ, δ, δ) and so on. The dashed line corresponds to Eq. (I2) with τ = 0.95. Remarkably, such a simple ansatz with a single parameter, already captures the integrand asymptotics well, in various directions in 5 dimensions. For a fixed number of integrand evaluations N , the error made in the calculation of Q n is very sensitive to the choice of the parameter τ as shown in Fig. 11. For N = 10 6 , Fig. 11 shows the relative error as a function (v, δ, δ, δ, δ) of τ for various orders n = 5, 6, 7 and 10. The error possesses a sharp minimum around τ = 0.85 (note the log scale). If τ is too small, we under-sample the tails of the function leading to potentially incorrect results. If τ is too large, the calculation is correct, but less efficient since the sampling puts a lot of points in regions which contribute little to the result. We will now see how to make the computation more robust regarding the choice of τ .

(v, v, δ, δ, δ) (v, v, v, δ, δ) (v, v, v, v, δ) (v, v, v, v, v) |f 5 |

Learning h (i) from the integrand

In this section, we make an additional change of variable w → x. To approximate the new integrand fn (w) in the w variable, we search for a new model function pn (w), again having the form:

pn (w) = n i=1 h(i) w i . (I3)
If we assume that the integrand fn (w) is well approximated by such a simple form, we estimate the functions h(i) by projecting in each dimension

h(i) proj (y) ≡ d n w fn (w)δ(w i -y). (I4)
In practice, we calculate M values of the integrand f (w α ), α ∈ 1...M generated by sampling the w space using the Sobol' sequence, and bin them in each dimension with N b bins. The function that we obtain is rather noisy due to The red solid curve is the new error after optimization of h (i) by projection (see section I 2), at order n = 7. The parameters are the same as in Fig. 1.

the binning, so in a second step, we use a Gaussian kernel smoothing

h(i) proj (w) → N b l=1 K λ w, l N b h(i) proj l N b N b l=1 K λ w, l N b (I5)
where

K λ (w, w ) = e -(w-w ) 2 /λ 2 . (I6)
The upper right plot of Fig. 12 shows an example of the smoothing procedure. The change of variable based on h(i) proj was used to compute the continuous line in Fig. 11 for order n = 7 (N b = 100, M = 10 5 , λ = 0.05). The error is improved at the optimum point τ = 0.95 only by approximately a factor 2. However, it remains largely independent of the value τ selected for the v → w change of variable, showing that the method has become much more robust. The projection automatically fits the exponential tails, without manual adjustment (as long as τ is not too small, to avoid under-sampling as explained above). Note that in Fig. 11, only M = 10 5 values of the integrand were used for the learning step. This is 1% of the total number of function calls, hence negligible. If τ is increased, the sampling of the integrand decreases in quality, as the tails are over-sampled. Correcting this in the learning step becomes increasingly more demanding. Therefore the change of variable v → w with a good initial guess for τ is essential to the success of the projection.

Because we only used model functions of the form Eq. ( 6), according to Eq. (D7), x m only depends on w m , which itself only depends on v m . In this special situation, the change of variable v → x can be represented by a model function of the form Eq. ( 6), and compared to the simpler v → w. As a result, Fig. 12 shows h (i) proj (v i ) along with the initial exponential guess, for all i and n ≤ 5, and for two values of τ = 0.95 and τ = 1.5. h (i) proj computed from two different initial guess for τ are indistinguishable (plain lines, on top of each other), showing that the result is independent of the initial choice of τ . At small v, they are quite different from a pure exponential. Finally, we note that the projected model functions vary only slightly with n for most values of i. This could be turned into an advantage by reusing the h functions from a lower dimension to a higher one in future developments.

Comparison to VEGAS algorithm

The form of model function Eq. ( 6) and the projection procedure described above are related to the VEGAS algorithm [START_REF] Lepage | A new algorithm for adaptive multidimensional integration[END_REF][START_REF] Lepage | VEGAS -an adaptive multi-dimensional integration program[END_REF]. However, simply applying this algorithm to perform the perturbation theory integrals produces poor results. Here we describe the essential differences between VEGAS and our approach.

The standard VEGAS algorithm uses random numbers to sample the integrand. It generates a weight function by projecting samples onto coordinate axes, which is analogous to the procedure we described. During the sampling steps, VEGAS uses the information it obtained to continuously improve the weight function.

To have a meaningful comparison, we will add physics information in the form of the u → v coordinate mapping (see above). The form Eq. ( 6) is only a good approximation to our integrand in the v variables. Although in our case VEGAS converges poorly in both u and v variables, it is far worse in u.

Our approach still differs from VEGAS in two essential ways. First, a key feature of our approach is to correctly capture the decaying asymptotic structure in the model function analytically. As we have shown in Fig. 11, a pure projection method is prohibitively inefficient in sampling the long tails and leads to poor convergence. It is only useful as a small correction to a good starting point. Second, the continuous change of model function of VEGAS cannot be used together with low-discrepancy sequences, since it breaks their special properties. For QQMC to work and achieve 1/N convergence, it is essential to construct the model function entirely before the full calculation is performed. Samples cannot be "reused" as in Monte Carlo sampling. This again requires an numerically efficient way to construct the model function which is not a priori provided by VEGAS.

Finally, we note that the idea of the model function is more general than the product form of Eq. ( 6) and VEGAS. Although they were sufficient for the Anderson impurity model calculations presented here, more general model function such as an MPS Eq. ( 7) are important for more complex many-body systems.

8 QQQMC, n = 6 QQMC, n = 8 DiagQMC, n = 6 DiagQMC, n = 8 Figure 1 .

 868681 Figure 1. Comparison of the convergence rates for QQMC and DiagQMC. Here Qn(N ) is the expansion coefficient of the occupation number of the Anderson impurity model at order n as a function of the number of integrand evaluations N . Each result is normalized to the exact analytic result Q Bethe n .

Figure 3 .

 3 Figure 3. Current at finite bias voltage through the Anderson impurity at T = 0, sweeping through several interaction regimes. Each point is a different QQMC calculation up to order n = 10, including series resummation [11]. The error bars are a combination of integration error and truncation error of the resummation; the latter dominates. By construction, the data is symmetric with respect to the particle-hole symmetric point ε d = -U/2. Inset: Coulomb diamond in the Coulomb blockade picture (large U ). Regions where current can flow are shaded grey. The dashed line indicates the scan shown in the main plot (varying ε d for fixed V /U = 1/7).

Figure 4 .

 4 Figure 4.Comparison of convergence of Q8 for different methods of integration: evaluating unwarped integrand with a Sobol' sequence (cyan), DiagQMC (red), warped integral sampled with Mersenne Twister pseudo-random numbers (green) or Sobol' sequence (blue). For the warped cases, we used Eq. (6) with h (i) (vi) = exp(-vi/τ ), τ = 0.95. After an initial warping with exponential functions τ = 1.1, we can apply an additional warping obtained by projection (orange); see Appendix I. For visibility, data (except Sobol' and DiagQMC) has been smoothed in the same way as in Fig.2.

Figure 5 .

 5 Figure 5. Sampling the two-dimensional square [0, 1] 2 with quasi-random numbers from a Sobol' sequence (left) and pseudo-random numbers from a Mersenne-Twister sequence (right). The Sobol' sequence gives a far more uniform sampling.

Figure 6 .

 6 Figure 6. Relative error (upper) an absolute error (lower) of Qn versus perturbation theory order n for our benchmark calculation against Bethe ansatz results, up to n = 22. Different curves correspond to different numbers of calculated points N . We used Eq. (6) with h (i) (vi) = exp(-vi/τ ) with τ = 0.95 for n ≤ 15 and τ = 0.90 for n ≥ 16. The parameters are the same as in Fig.1.

  Figure 7.Absolute error of occupation number Q6(N ), with curves matching three cases of Fig.4. Here we perform calculations with K = 25 shifted sequences. The pale curves show the convergence of the average over these runs to the Bethe ansatz result. The solid curves show the corresponding uncertainty estimates computed from the standard error of the mean over the shifted sequence values. We see that this procedure gives good estimates of the error.

7 Figure 8 .

 78 Figure 8. Convergence of the current In with number of samples N at different orders n. This data corresponds to a typical point of Fig.3, away from particle-hole symmetry (ε d /U ≈ -0.12, α ≈ 0.27, U = 6Γ). At low order (blue and orange lines) the relative error scales as 1/N 0.9 (black plain line). At larger orders (n = 7 and 10, green and red lines), the convergence slows down and scales only as 1/N 0.7 (dashed black line). Finally at order n = 10, the final scaling starts at a larger number of function evaluations N than for lower orders. For visibility, the data has been smoothed as described in Appendix H (see main text).

Figure 9 .

 9 Figure 9. Perturbation series of the current for the same parameters as in Fig. 8. The upper line (colored symbols)gives the absolute values of the series coefficients |In|: positive coefficients are depicted in blue, and negative ones in red. The lower line (black symbols) is the estimated absolute error ∆In. For all orders calculated, up to n = 10, we obtained at least two digits of accuracy.

Figure 10 .

 10 Figure 10. Comparison between the absolute value of the order 5 integrand |f5(v)| in Eq. (2) (colored lines) and the model function p5(v) of Eq. (6 ) (black dashed lines, τ = 0.95) along various directions in v-space. Each color corresponds to v = (v, . . . , δ) as indicated in the figure. Lines of the same color correspond to different permutations within a given direction (see main text). The parameters are the same as in Fig. 1.

7 Figure 11 .

 711 Figure 11. Exact relative error of Qn computed with N = 10 6 points as a function of the model function parameter τ . The dashed lines are for a model defined by a single exponential function h (i) (vi) = exp(-vi/τ ), at different orders.The red solid curve is the new error after optimization of h (i) by projection (see section I 2), at order n = 7. The parameters are the same as in Fig.1.

Figure 12 .

 12 Figure 12. Examples of functions h (i) (vi) for various n and 1 ≤ i ≤ n: exponential form (dashed lines) and after projection (plain lines, on top of each other). We used N b = 250 bins, M = 20000 samples and λ = 0.05; see text. An example of binned values and smoothing for h

  Figure 2. Expansion coefficients Qn for the Anderson impurity occupation number relative to the analytic result Q Bethe =↑, ↓ is the electronic spin and ε d represents a capacitive gate coupled to the dot. The local Coulomb repulsion is

	n	.
	QQMC converges at rates close to 1/N with the number of
	integrand evaluations N . For visibility, the data has been
	smoothed (see Appendix H). The black lines indicate exact 1/N (dotted) and 1/ √ N (dashed) convergence. Each run
	was performed with one Sobol' sequence. Inset: Cartoon of
	quantum dot set-up.

σ

Table I .

 I Perturbation coefficients for the occupation number from the Bethe Ansatz.

	2 5 -9450-1175π 2 +18π 4 6π 4 90π 5 6 945 27π 4 + 137 60π 2 π 6 -3185 7 -10395 π 7 + 11690 9π 5 -2653 100π 3 + 1 7π 8 135135 π 8 -912065 54π 6 + 524881 140π 2 1500π 4 -363 9 -2027025 π 9 + 41046005 162π 7 -118576073 22500π 5 + 623873 14700π 3 -1 9π 10 34459425 π 10 -348898550 81π 8 + 5045327287 56250π 6 -761337511 1029000π 4 + 7129 2520π 2 11 -654729075 π 11 + 19887342475 243π 9 -287681226833 168750π 7 + 306048256943 21609000π 5 -4785253 79380π 3 + 1 11π 12 13749310575 π 12 -2505809831525 1458π 10 + 90628717412233 2531250π 8 -84543422632097 283618125π 6 + 131145705977 100018800π 4 -83711 27720π 2	-0.03176576952402088647648727 0.01428453978718527764087931 0.003296943760155087659102276 -0.007417285598070720865765088 0.003074980697930490893903131 0.001257173321207511996877863 -0.002096420663210080104411625 0.0007382700406215484131765644 0.0004844375173694184755334739

ACKNOWLEDGMENTS

We thank N. Andrei, L. Greengard, E.M. Stoudenmire, N. Wentzell and especially A.H. Barnett for helpful discussions. The algorithms in this paper were implemented using code based on the TRIQS library [40] and the QMCgenerators library [27]. The Flatiron Institute is a division of the Simons Foundation. XW and MM acknowledge funding from the French-Japanese ANR QCONTROL, E.U. FET UltraFastNano and FLAG-ERA Gransport. (2) proj for order n = 3, 6 is shown in the upper right inset. The parameters are the same as in Fig. 1.

Appendix G: Error Calculation in Quasi-Monte Carlo

If we have reached the asymptotic regime for a single sequence Q(N ), we can estimate an error by fitting the approach to convergence. This is similar to error estimation for traditional quadrature.

A more robust estimate of the error can be achieved by reintroducing a random component to the methodso-called Randomized Quasi-Monte Carlo [START_REF] Dick | Highdimensional integration: The quasi-Monte Carlo way[END_REF][START_REF] Nuyens | The construction of good lattice rules and polynomial lattice rules[END_REF][START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF][START_REF] Pierre | Randomized quasi-Monte Carlo: An introduction for practitioners[END_REF]. Here we repeat the calculations with K 'randomly shifted' sequences giving a distribution of values Q k , from which we obtain an estimate of the mean and error as for conventional Monte Carlo. Typically one chooses only a moderate number of K ∼ 10 -100 as it is advantageous, for fixed computational time, to maximize N . In Fig. 7, we show this approach in practice for different methods of Fig. 4.

Appendix H: Details on the current I calculation

This appendix gives details of the calculations made to obtain Fig. 3.

In this application, each lead has centered, semi-circular density of states with half-bandwidth D = 2. The coupling to the leads is chosen so that Γ = 0.1. A symmetric bias voltage V is applied between the two leads. For each parameter set (ε d , V, U ), E d and α are chosen so that the first order (Hartree term) of the perturbation series vanishes at U = 7Γ, to improve its convergence radius. The series, defined in Appendix E, is computed at a time t = 10/Γ after switching the interaction on.

Due to the finite bandwidth D, the integrand decays polynomially at large times, but exponentially at intermediate times. We used an exponential warping with τ = 1.5/Γ, which was enough to capture the general shape of the integrand up to the observation time t.

Figure 8 shows the convergence of the calculation for In this appendix, we explain how the h (i) functions are constructed in the model function defined in the main text,

for 0 < u n < u n-1 < . . . < u 2 < u 1 and u 0 = t is defined to be the measurement time.

We are going to use successive changes of variable in this section:

First, we change to the v i variables defined as v i ≡ u i-1u i > 0 which are natural since the integrand only depends on time differences. The corresponding Jacobian is one. Note that the integration on v is performed on [0, ∞) n . When going back to u space, this generates extra points that are not in the original u domain. The value of the integrand for these points is simply zero so that they induce no extra computational cost.


Marion Le Guellec

Florent Lebreton

Laure Ramond

Abibatou Ndiaye

Thierry Gervais

Guillaume Bernard-Granger
email: guillaume.bernard-granger@cea.fr

Sintering investigations of a UO2-PuO2 powder synthesized using the freeze-granulation route

Keywords: Sintering, ceramics, oxides, microstructure, MOX

2

MOX fuels (UO2-PuO2) are used in light water nuclear reactors and are potential candidates for fast neutron reactors. Industrial MOX pellets are manufactured by a dry-route process (grinding, pressing and sintering). Sintering at high temperature of MOX fuels was investigated in the past in function of the plutonium oxide content, in different kinds of atmospheres and sometimes using hot-pressing [1][2][3][4][5][6]. However the obtained results were not very conclusive in identifying the mechanisms controlling densification and microstructure development. Recently, the freeze-granulation route was shown to be a promising approach to manufacture MOX fuel pellets [7]. Highly flowable, dustless and easy-to-press MOX granules were elaborated [7]. Sintering of green compacts gives highly dense and defect-free pellets exhibiting a very homogeneous U-Pu spatial distribution [7]. The purpose of this paper is to rigorously investigate the sintering behavior of powder compacts made of such granules.

A homogeneous UO2-PuO2 powder, containing 11 wt% PuO2 was prepared using the freeze-granulation route optimized by La Lumia [7]. The UO2 powder used was synthesized through a liquid route from uranyl nitrate solution. The PuO2 powder was obtained by oxalic precipitation of plutonium (IV) nitrate and calcination. Supplementary material shows the morphology of some individual granules observed using SEM (scanning electron microscopy, JCM-6000Plus, Jeol). They have a spherical shape and are devoid of any central cavity. Laser particle size distribution measured in dry mode (Mastersizer 3000, Malvern Panalytical, Mie configuration, air pressure fixed to 0.1 bar) is also shown in supplementary material. The volume distribution appears monomodal. The volume mean (D [4;3]) and median (Dv(50)) diameters are 915 and 775 µm, respectively. The Pu/(U+Pu) content in the freeze-granulated powder was determined to be 10.5 at% using thermal ionization mass spectrometry (TIMS, VG-54 magnetic sector mass spectrometer, Isotopx).

12/7/2020 3 Cylindrical samples (diameter of 5.4 mm and a height of 7 mm) were obtained by uniaxial pressing the freeze-granulated powder with a compaction pressure set to 450 MPa (manual hydraulic press, Atlas 25, Specac). By taking a theoretical volume mass of 11.034 g/cm 3 at room temperature, the relative green density was around 54-55 %.

The green pellets were sintered from room temperature to 1700 °C in a high temperature dilatometer (DIL402C, Setaram) under Ar/4 vol% H2. Dioxygen partial pressure of the atmosphere was monitored and controlled with an oxygen pump (Gen'air, Setnag) in order to have a water concentration of 1200 vpm during all the sintering runs (the oxygen partial pressure, calculated using the Wheeler approach [8], is fixed at 5.0x10 -26 bar at 650 °C, which gives 5.8x10 -11 bar at 1700 °C). Different heating rates: 2, 5 and 10 °C/min were used. All the runs have been completed using a heating rate fixed to 2 °C/min, from room temperature to 600 °C, in order to eliminate the organic compounds (dispersant and binder) used to elaborate the freeze-granulated powder. Because the mass of the samples does not change significantly between the green and fired states (small weight loss accounting for the organic compounds and adsorbed humidity present in the granules), the following relation is obtained [9][10]:

𝐷 (𝑇) (%) = 𝐷 0 1 𝜉 2 [1 + ∆𝐿 (𝑇) 𝐿 0 ] 3 𝑒 3𝛼(𝑇-𝑇 0 ) (1)
where:𝐷 (𝑇) is defined as the sample instantaneous relative density, 𝐷 0 is the relative green density, 𝜉is representative of the shrinkage anisotropy of the sample and is given by 𝜉 = Fig. 1a shows how the relative density is progressing in function of temperature.

Whatever the heating rate, samples start to densify around 650 °C. The lower the heating rate, the higher the relative density for a given temperature. Such a behavior is well known and was previously observed for numerous nonradioactive oxides as zirconia, alumina and spinel [9][10][11][12]. As proposed by Sato, the lower the heating rate, the longer the exposure time above the sintering onset temperature and the higher the shrinkage and relative density reached at a given temperature [11]. Fig. 1b shows the variation of the instantaneous densification rate as a function of temperature. The higher the heating rate, the higher the instantaneous densification rate whatever the temperature. This trend is also well known and was previously reported for non-radioactive oxides [9][10]13]. It has also to be outlined that the higher the heating rate, the higher the temperature of 1050°C, the samples are more homogeneous and can be described as a polycrystalline body constituted by pores and grains of U1-xPuxO2±x solid solutions, rather than UO2 and PuO2 grains. We are aware that the previous hypothesis needs to be corroborated in the future using, for example, X-ray diffraction (XRD) experiments at high temperature and interrupted sintering tests around 830 °C followed by careful microstructure characterizations.

Indeed, other reasons may lead to a maximum in the densification rate for low temperatures, as for example stepwise particle rearrangement at early sintering offset. Interesting is also the tendency given by Fig. 1d showing the variation of the absolute instantaneous densification rate in function of 𝐷(𝑇). All curves exhibit a maximum value around 1.2x10 -5 /K when the instantaneous relative density is in the range 67-70%, which corresponds more or less to the end of the step related to the necks formation between the elemental crystallites and/or the individual freeze-dried granules. In addition, the three curves overlap. This means that the densification mechanism is the same for the three heating rates used.

Figs exhibit an intergranular fracture surface (Figs. 2a and2b). Conversely, samples with a relative density above 90% give an intragranular fracture surface (Fig. 2c). For samples giving an intergranular fracture, as is also the case for green compacts, it is easy to determine the average grain size by the intercept method (at least 300 grains are considered, correction factor set to 1.56 even if such a value is usually used for planar surfaces which is not strictly the case of a 12/7/2020 6 fracture surface) [16]. For samples with an intragranular fracture, the grain size is only approximated by finding zones where the elementary grains are observed unequivocally (around 40 grains only, same correction factor). Knowing the relative density and the grain size of each sintered sample, it is possible to construct the "grain size versus relative density" trajectory named the "sintering path", as shown on Fig. 2d. The experimental points belong to a single trajectory (dashed curve on Fig. 2d), independent of sintering temperature, soak time and heating rate. Accordingly, the grain size can be described as a monotonous function of the relative density and a multitude of sintering experimental conditions led to the same microstructure, as was previously reported for other non-radioactive oxides [9][10]13]. From Fig. 2d, it also appears that grain growth is moderate when the relative density is lower than 90% and is on contrary pronounced above. has been determined to be 10.5 at% using TIMS (same value as in granules).

Philipponneau proposed a relationship linking the lattice parameter, 𝑎, to the Pu/(U+Pu) content, 𝑦, and the deviation from stoichiometry, 𝑥, for U1-yPuyO2-x solid solutions [17]:

𝑎(𝑝𝑚) = 547 -7.4𝑦 + 32𝑥 (2)
According to relation (2) we obtain a 𝑥 value close to 1x10 -3 . Thus the average composition of the sintered polycrystal is U0.895Pu0.105O1.999. The O/M value is almost 2, which is compatible to the one required for MOX fuels dedicated to nuclear light water reactors.

All the experimental points constituting the sintering path shown on Fig. 2d were fitted with theoretical expressions linking the grain size to the relative density [18]. The best determination coefficient value amongst these fits is obtained when grain growth is controlled by the grain boundaries and densification is controlled by volume diffusion (see fit on Fig. 3a).

Let us now determine the apparent activation energy for the mechanism controlling densification supposed to be volume diffusion. The equation for the densification rate can be separated into temperature-dependent, grain sizedependent, density-dependent quantities and rearranged as follows [10,12,19]:

𝐿𝑛 [𝑇 𝑑𝐷 𝑑𝑇 𝑑𝑇 𝑑𝑡 ] = - 𝑄 𝑑 𝑅𝑇 + 𝐿𝑛{𝐹[𝐷]} + 𝐿𝑛(𝐴) -𝑛𝐿𝑛(𝐺) (3)
where 𝐴 is a constant, 𝐹[𝐷] is a function only of density, 𝑄 𝑑 is the apparent activation energy for the mechanism controlling densification, 𝑅 is the universal gas constant, 𝑇 is the absolute temperature, 𝐺 is the grain size, 𝑛 is the grain size exponent whose value depends on whether the densification rate is controlled by lattice diffusion (n = 3) or grain-boundary diffusion (n = 4) and 𝜁 = 𝑑𝑇 𝑑𝑡 ⁄ is the heating rate that is held constant during anisothermal sintering experiments.

A plot of the left-hand side of relation ( 3) versus 1/𝑇 would give a value for 𝑄 𝑑 provided that the data points are taken at a constant value of 𝐷 and 𝐺. A second method to determine the apparent activation energy for the mechanism controlling densification is to construct the Master Sintering Curve (MSC). By extending the analysis of sintering beyond the confined segments described by the individual stage models [20][21], it was shown, for isotropic shrinkage (true for the freeze-granulated powder investigated) and if there exists only one dominant diffusion mechanism (either volume or grain boundary diffusion), that [20]:

𝑅 3𝛾 𝑠𝑣 ΩΦ 0 ∫ (𝐺[𝐷]) 𝑛 𝐷Γ[𝐷] 𝐷 𝐷 0 𝑑𝐷 = ∫ 𝑒 -𝑄 𝑑 𝑅𝑇 𝑇 𝑑𝑡 𝑡 0 = Θ 𝑑 (ℎ/𝐾) = 1 𝜁 ∫ 𝑒 -𝑄 𝑑 𝑅𝑇 𝑇 𝑑𝑇 𝑇 𝑇 0 ( 4
)
where 𝛾 𝑠𝑣 is the solid/vapor surface tension, Φ0 = Φv0 and n = 3 for volume diffusion, Φ0 = bΦb0with b the grain boundary thickness and n = 4 for grain boundary diffusion (Φ0 is the pre-exponential factor of the diffusion coefficient),

Γ[𝐷]and 𝐺[𝐷]

(the grain size) are functions only of the relative density, 𝑇 0 is the temperature at which densification is effectively starting (assumed to be 1000 °C at the end of solid solution formation) and the other parameters have been described previously. Because the left-hand site of Equation ( 4) is a function only of relative density, the MSC is simply defined as the relationship between the relative density 𝐷 and Θ 𝑑 [20]. To construct the MSC curve under anisothermal sintering conditions, we used the results from the dilatometer runs, the trapezoidal approximation and a curve-fitting procedure to determine the value of 𝑄 𝑑 by minimizing the residual sum of squares (RSS) between the experimental and calculated values, assuming that all points must line up to form a single curve. Such an analysis leads to the results shown in Figs. 3c and3d.

The apparent activation energy for the mechanism controlling densification is then calculated to be around 600 kJ/mol.

12/7/2020 9 Thereby, an average apparent activation energy equal to 630 ± 25 kJ/mol, associated with the densification mechanism, is therefore determined from the CRH and MSC methods (average value and standard deviation calculated by taking into account the five values obtained from the CRH method and the single value calculated from the MSC analysis). Tab. 1 summarizes the volume diffusion coefficients and corresponding activation energies for uranium/plutonium cations and oxygen anions for U1-yPuyO2±x [22][23][24][25].

Plutonium cations are the slowest species and the activation energies measured by Noyau [24] and Lamb [25] are very close to the one measured during our sintering experiments. Accordingly, we postulate that volume diffusion of plutonium cations is controlling densification of the freeze-granulated powder investigated.

Finally, a polished cross section from the sample sintered for 4 h at 1700 °C (heating rated fixed to 2 °C/min) was characterized using electron probe microanalysis (SX100, Cameca). Measurements performed at 20 kV at the U-Mα and Pu-Mβ lines allowed the U-Pu distribution to be assessed and the Pu content was quantified by applying the ZAF method and using UO2 and PuO2 samples as standards for U, Pu and O. The resulting chemical mapping in Fig. 4 shows that the Pu repartition inside the pellet is very homogeneous, notably compared to pellets produced by conventional powder metallurgy processes [START_REF] Gervais | Proceedings of the International Conference on Fast Reactors and Related Fuel Cycles: Next Generation Nuclear Systems for Sustainable Developments[END_REF][START_REF] Oudinet | [END_REF]. MOX fuels with a very homogeneous U/Pu distribution are thought to have a better behavior under irradiation in reactors by preventing the formation of the high burnup structure [28][29]. 0.0E+00 Tab.1

Figure captions

  the maximum in densification rate: 1421 °C for 2 °C/min, 1461 °C for 5 °C/min and 1547 °C for 10 °C/min. Fig.1cis a zoom of Fig.1bin the 600-1200 °C temperature range. Whatever the heating rate (and even if it is more delicate to observe for the lowest heating rate), the instantaneous densification rate exhibits a first maximum value around 830 °C. It means that densification starts around 650 °C and slows-down from 830 °C. Similar results were reported by several authors investigating the sintering behavior in different kinds of atmospheres of UO2-20 wt% PuO2[1], UO2-30 wt% PuO2[14] fuels and UO2-30 wt% CeO2 surrogates[15]. For them, the densification slow-down correlates with the onset of a UO2-PuO2 or UO2-CeO2 solid solution formation. Even if additional experiments are required, we postulate that a solid solution forms in the densifying samples around 830 °C, whatever the heating rate used. Because the solid solution is formed by the interdiffusion of Pu 4+ cations into UO2 lattice and U 4+ cations into PuO2 lattice it decreases the densification rate temporarily until 12/7/2020 5 completion and then densification restarts with a polycrystalline sample made of U1-xPuxO2±x grains. Accordingly, when densification resumes around 950-

  Fig. 2d displays that each value of relative density corresponds to a single value of grain size, regardless of the sintering conditions. Then, points for a constant value of 𝐷 are generated by changing the heating rate. The measurements led to values of 𝑄 𝑑 at different values of 𝐷. This formalism is known as the CRH method (Constant Rates of Heating). Using such an approach leads to Fig. 3b and an 12/7/2020 8 almost constant apparent activation energy of 635 ± 25 kJ/mol is calculated for the mechanism controlling densification.

Fig. 1 :Fig. 2 :Fig. 3 :

 123 Fig. 1: a) Relative density in function of temperature; b) Densification rate in function of temperature; c) zoom of b) in the 600-1200 °C temperature range; d) Absolute densification rate in function of temperature.

  Fig.1

Acknowledgements:

The authors warmly thank Gauthier Jouan, Romain Lauwerier, Romain Vauchy and Patrice Signoret for their precious help during the experiments carried out within the framework of this study.

12/7/2020