# Heating and cooling in self-consistent many-body simulations

Yang Yu  Department of Physics, University of Michigan, Ann Arbor, MI 48109,
USA  Sergei Iskakov  Department of Physics, University of Michigan, Ann Arbor,
MI 48109, USA  Emanuel Gull  Department of Physics, University of Michigan,
Ann Arbor, MI 48109, USA

###### Abstract

We present a temperature extrapolation technique for self-consistent many-body
methods, which provides a causal starting point for converging to a solution
at a target temperature. The technique employs the Carathéodory formalism for
interpolating causal matrix-valued functions and is applicable to various
many-body methods, including dynamical mean field theory, its cluster
extensions, and self-consistent perturbative methods such as the self-
consistent GW approximation. We show results that demonstrate that this
technique can efficiently simulate heating and cooling hysteresis at a first-
order phase transition, as well as accelerate convergence.

##  I Introduction

When conducting experiments in condensed matter physics, it is common practice
to investigate the temperature dependence of observables while keeping other
parameters constant. A particular example are specific heat or transport
measurements, which are often used as a preliminary probe to identify
intriguing temperature-dependent behavior. For example, in a first-order
coexistence regime, heating and cooling curves may reveal history-dependent
hysteresis.

In theoretical calculations using self-consistent finite-temperature field
theories, changing the temperature of a system is generally not practical as
it causes a shift of Matsubara frequencies [1, 2]. Extrapolating the lowest
Matsubara frequency during cooling can be problematic, resulting in non-causal
solutions. Therefore, implementing heating and cooling protocols in self-
consistent finite-temperature field theories, in analogy to heating and
cooling measurements in experiment, is an unresolved issue.

Self-consistent finite-temperature simulation methods include non-perturbative
techniques, such as the dynamical mean field theory (DMFT) [3] and its cluster
variants [4, 5, 6, 7]; self-consistent perturbative methods, such as GW [8, 9,
10, 11, 12, 13], second-order perturbation theory [14, 15, 16, 17],
fluctuation-exchange or T-matrix approximations [18, 19], as well as bold-line
diagrammatic Monte Carlo methods [20, 21, 22, 23]; and combinations of
embedding theories with perturbation theory [24, 25, 26].

In these methods, solutions are obtained through an iterative process that
involves starting with an initial guess and continuing the process until self-
consistency is achieved. The number of iterations required to reach
convergence is directly related to how close the starting point is to the
iteration fixed point. A ‘good’ starting point can significantly reduce the
computational effort required for the simulation, whereas a ‘bad’ starting
point may result in iterations that diverge, iterate in limit cycles, or even
converge to unphysical fixed points [27]. In parameter regimes where first-
order coexistence occurs, multiple physical fixed points may exist [28, 29,
30, 31, 32].

This paper addresses the challenge of generating better starting points and
implementing heating and cooling protocols in self-consistent many-body
methods. It presents a solution that guarantees a causal starting point by
using a converged solution at a different temperature. The proposed method
relies on the Carathéodory formalism [33] for interpolating causal matrix-
valued functions. Originally developed for analytic continuation of matrix-
valued Matsubara functions to real frequencies, this formalism can be extended
to temperature extrapolation, which involves evaluating an interpolant at
different Matsubara frequencies. We demonstrate the effectiveness of this
approach in obtaining improved starting points in the context of DMFT and
real-materials perturbation theory. Additionally, we examine heating and
cooling hysteresis in the context of a first-order phase transition.

The paper proceeds as follows. In Sec. II we introduce the formalism. Sec. III
contains a brief description of a pedagogical implementation of the method
which we provide as a supplement. Sec. IV contains results for temperature
extrapolation, convergence acceleration, and first-order hysteresis. Finally,
Sec. V contains our conclusions.

Figure 1: Schematic representation of the interpolation procedure for a
matrix-valued finite-temperature Matsubara function $G(\mathrm{i}\omega_{n})$.
See text for detailed description.

##  II Method

The central object of this paper is a causal matrix-valued fermionic Matsubara
function which may represent a Green’s function, a self-energy, or a cumulant
[34]. The Matsubara function is expressed as a three-dimensional tensor,
$G_{ij}(\mathrm{i}\omega_{n})$, which associates a matrix (identified by the
indices $i$ and $j$) with every Matsubara frequency
$\omega_{n}=(2n+1)\pi/\beta$ ($n$ denotes an integer, $\beta$ the inverse
temperature). Its continuation to the complex plane [33] is denoted as
$G_{ij}(\tilde{z})$. $\mathrm{i}G(\tilde{z})$ is a Carathéodory function (up
to a convention-dependent minus sign), such that
$[\mathrm{i}G(\tilde{z})+(\mathrm{i}G(\tilde{z}))^{\dagger}]/2$ is a positive
semi-definite matrix for any $\tilde{z}$ in the upper half of the complex
plane.

The matrix structure of $G$ depends on the application. Matrices may be
scalar, diagonal, block-diagonal, or fully dense. The scalar case typically
appears in single-site single-orbital cases, such as single-site DMFT. The
diagonal case is often encountered in momentum-space simulations, such as
cluster DMFT [4, 5, 6, 7]. The general multi-orbital case with dense matrices
typically appears in real-materials ab-initio calculations [25, 26].

The frequency-dependence of $G$ is known at discrete Matsubara frequencies
$\omega_{n}$. In the simplest case, the frequencies are uniformly spaced
positive fermionic Matsubara frequencies $\omega_{n}=(2n+1)\pi/\beta$,
$n=0,1,\cdots N-1$. More generally, data is provided on a non-uniform
frequency grid using Chebyshev [35], intermediate representation [36],
Legendre [37], spline [38], or discrete Lehmann representation [39] schemes
with associated Fourier transforms [40, 41]. Bosonic Matsubara response
functions $\Pi(\tilde{z})$ are not directly related to Carathéodory functions.
However, $\Pi(\tilde{z})\tilde{z}$ and $\Pi(\tilde{z})/\tilde{z}^{*}$
correspond to Carathéodory functions [42], and Nogaki and Shinaoka [43]
recently investigated a related mapping of bosonic functions.

Extrapolating in temperature requires the transformation of one Matsubara grid
to another one corresponding to a different temperature. We propose to
construct a causal matrix-valued interpolant through $G(\mathrm{i}\omega_{n})$
and evaluating it at the Matsubara points corresponding to the changed
temperature grid to realize the temperature extrapolation. The method’s
algorithmic steps follow those presented for analytic continuation in Ref.
[33], with the only distinction that the interpolant is evaluated on the
imaginary axis, rather than just above the real axis. The major steps of the
algorithm are described below.

For a Carathéodory function $\tilde{F}(\tilde{z})=\mathrm{i}G(\tilde{z})$,
$\tilde{z}$ in the upper half of the complex plane, we assume the values of
$\tilde{F}$ are known at a set of points
$\\{\tilde{z}_{n}|n=0,1,\cdots,N-1\\}$. A Möbius transform

$z=\frac{\tilde{z}-\mathrm{i}}{\tilde{z}+\mathrm{i}}$ (1)

maps the upper half of the complex plane to the unit disk, the real axis to
the unit circle, and the Matsubara points to points onto the real axis. Since
a Carathéodory function $F$ is a function defined on an open subset of the
complex plane and exhibits the property $(F+F^{\dagger})/2$ being positive
semi-definite, $C(z)=\tilde{F}(\tilde{z})$ then defines a Carathéodory
function on the unit disk $|z|<1$. The function $C(z)$ is known at a set of
points $\\{z_{n}|n=0,1,\cdots,N-1\\}$ with the values of

$\displaystyle C_{n}\equiv C(z_{n})=\tilde{F}(\tilde{z}_{n}).$ (2)

The Cayley transform [44] and its inverse, defined as

$\displaystyle S(z)$
$\displaystyle=\left[I-C\left(z\right)\right]\left[I+C\left(z\right)\right]^{-1},$
(3) $\displaystyle C(z)$
$\displaystyle=\left[I+S\left(z\right)\right]^{-1}\left[I-S\left(z\right)\right],$
(4)

map the Carathéodory function $C(z)$ to a matrix-valued Schur function $S(z)$
[45] and vice versa. A Schur function is defined as a function $S(z)$ on the
unit disk with the property $||S(z)||\leq 1$, where the matrix norm
$||\cdot||$ is defined as the largest singular value of $S$, or equivalently
the largest eigenvalue of $(SS^{\dagger})^{1/2}$.

The problem of interpolating the Carathéodory function $C(z)$ is thus
converted into the problem of determining an interpolant for the Schur
function $S(z)$, where $S(z)$ is known at $N$ discrete points, given as
$S_{n}\equiv S(z_{n})$, with $n=0,1,\cdots,N-1$.

To identify the interpolant of the Schur function $S(z)$, we proceed to find a
set of Schur functions $S^{0}(z),S^{1}(z),\cdots,S^{m}(z),\cdots,S^{N}(z)$,
where each Schur function $S^{m}(z)$ is known at $N-m$ points
$z=z_{m},z_{m+1},\cdots,z_{N-1}$ with the values of $S^{m}_{n}\equiv
S^{m}(z_{n})$. $S^{N}(z)$ is an arbitrary Schur function without any
constraints. We require $S^{0}_{n}=S_{n}$ such that $S^{0}(z)$ serves as the
desired interpolant for the Schur function $S(z)$.

To establish a connection between the two Schur functions $S^{m}(z)$ and
$S^{m+1}(z)$, a function $B^{m}(z)$ is introduced, enabling a progressive
determination of $S^{0}(z)$ starting from $S^{N}(z)$ [33]. The relationship
between $B^{m}(z)$ and $S^{m}(z)$ is given by:

$\displaystyle\frac{|z_{m}|(z_{m}-z)}{z_{m}(1-z_{m}^{*}z)}B^{m}(z)$
$\displaystyle=L^{m}[S^{m}(z)-S^{m}_{m}]$ (5)
$\displaystyle\times[I-{S_{m}^{m}}^{\dagger}S^{m}(z)]^{-1}R^{m}$

with

$\displaystyle L^{m}$
$\displaystyle=[I-S^{m}_{m}{S^{m}_{m}}^{\dagger}]^{-\frac{1}{2}},$ (6)
$\displaystyle R^{m}$
$\displaystyle=[I-{S^{m}_{m}}^{\dagger}S_{m}^{m}]^{\frac{1}{2}}.$ (7)

An equivalent form of this relation is

$S^{m}(z)=[I+V^{m}(z){S^{m}_{m}}^{\dagger}]^{-1}[V^{m}(z)+S^{m}_{m}]$ (8)

with

$V^{m}(z)=\frac{|z_{m}|(z_{m}-z)}{z_{m}(1-z_{m}^{*}z)}(L^{m})^{-1}B^{m}(z)(R^{m})^{-1}.$
(9)

The relation between $B^{m}(z)$ and $S^{m+1}(z)$ reads

$\displaystyle S^{m+1}(z)$
$\displaystyle=[I-K^{m}{K^{m}}^{\dagger}]^{-\frac{1}{2}}[B^{m}(z)-K^{m}]$ (10)
$\displaystyle\cdot[I-{K^{m}}^{\dagger}B^{m}(z)]^{-1}[I-{K^{m}}^{\dagger}K^{m}]^{\frac{1}{2}},$

where $K^{m}$ is an arbitrary matrix with $||K^{m}||<1$.

The freedom in choosing arbitrary $K^{0}$, $K^{1}$, $\cdots$, $K^{N-1}$, along
with $S^{N}(z)$, allows for a full coverage of all possible interpolants for
$S(z)$. The seemingly undetermined quantities in the aforementioned equations
are the values of $\set{S^{m}_{m}}{m=0,1,\cdots,N-1}$. Starting from the given
$\set{S^{0}_{n}}{n=0,1,\cdots,N-1}$, as $S^{m+1}_{n}$ are entirely determined
by $S^{m}_{n}$ and $S^{n}_{n}$ by proceeding from Eq. (5) to Eq. (10), all
$S^{m}_{m}$ can be acquired in the process of calculating
$\set{S^{m}_{n}}{m=0,1,\cdots N-1;n=m,m+1,\cdots,N}$ from $m=0$ to $m=N-1$.
After obtaining all $S^{m}_{m}$ and the corresponding $L^{m}$ and $R^{m}$, one
can reverse the procedure, transitioning from Eq. (10) to Eq. (8), to compute
$S^{m}(z)$ from $S^{m+1}(z)$ with a chosen $S^{N}(z)$ as the starting point,
until the desired interpolant $S^{0}(z)$ is attained. Finally, an interpolant
for $C(z)$ can be derived using Eq. (4), and by multiplying $-\mathrm{i}$ and
substituting $z$ with $\tilde{z}$ through the Möbius transform, the
interpolant for $G(\tilde{z})$ is obtained.

A sketch of the algorithm is shown in Fig. 1 (starting top left, following the
arrows, terminating bottom left). In practice, we take
$\\{\tilde{z}_{n}\\}=\\{\mathrm{i}\omega_{n}|n=0,1,\cdots,N-1\\}$ along with
the corresponding $G(\mathrm{i}\omega_{n})$ as the input where
$\omega_{n}=(2n+1)\pi/\beta$. Setting $\tilde{F}=\mathrm{i}G$, and with Eqs.
(1), (2), and (3), we obtain the constraints $\\{(z_{n},S_{n})\\}$ for the
interpolation problem of a Schur function $S(z)$. By taking $K^{m}=0$ for all
$m$s, from Eq.(10), we obtain $S^{m+1}(z)=B^{m}(z)$ which leads to

$S^{m+1}_{n}=\frac{z_{m}(1-z_{m}^{*}z_{n})}{|z_{m}|(z_{m}-z_{n})}L^{m}[S^{m}_{n}-S^{m}_{m}][I-{S_{m}^{m}}^{\dagger}S^{m}_{n}]^{-1}R^{m}$
(11)

according to Eq. (5). This equation offers a concrete way to calculate all
$S^{m}_{n}$, starting from $S^{0}_{n}=S_{n}$. Subsequently, we can obtain
$S^{m}_{m}$ and the corresponding $L^{m}$ and $R^{m}$. Assume the desired new
Matsubara points for a different temperature $\beta^{\prime}$ are
${\mathrm{i}\omega^{\prime}_{n}=\mathrm{i}(2n+1)\pi/\beta^{\prime}}$. Using
the Möbius transform, the corresponding points in the domain of $S(z)$ are
obtained as ${z^{\prime}_{n}}$. By setting $S^{N}(z)=I$, and using the
obtained $\\{S^{m}_{m},L^{m},R^{m}\\}$, the values of ${S_{n^{\prime}}\equiv
S^{0}(z^{\prime}_{n})}$ can be determined by proceeding from Eq. (10) to Eq.
(8). With the inverse Cayley transform in Eq. (4) and
$G=-\mathrm{i}\tilde{F}$, the desired values of the causal interpolant
${G(\mathrm{i}\omega^{\prime}_{n})}$ at ${\mathrm{i}\omega^{\prime}_{n}}$ are
obtained as output.

##  III Implementation

We provide a straightforward pedagogical implementation of the algorithm as
supplement to this paper 111See Supplemental Material at [URL will be inserted
by publisher] for Python implementation of the algorithm, documentation, and
examples.. The code is written in the python programming language, using no
dependencies other than numpy. Unlike related real-frequency analytic
continuation codes [33, 47], we found that an implementation in standard
double precision was sufficient to perform all calculations, and that the
temperature extrapolation of noisy Monte Carlo data presented no difficulties.

The implementation requires a Green’s function (or, equivalently, self-energy
or cumulant) $G[n,i,j]=G_{ij}(\mathrm{i}\omega_{n})$ as a three-dimensional
tensor and $w[n]=\mathrm{i}\omega_{n}$ as a one-dimensional vector. After
reading the data, all preprocessing steps indicated by the red arrows in Fig.
1 are performed, and $\\{S^{m}_{m},L^{m},R^{m}\\}$ are obtained and stored.

The code also requires a set of frequencies
$wp[n]=\mathrm{i}\omega^{\prime}_{n}$, provided as a one-dimensional vector.
These frequencies should correspond to the Matsubara points at the new
temperature $\beta^{\prime}$. Using this data, the code evaluates the
interpolated Green’s function at the new Matsubara points, according to the
blue arrows in Fig. 1 and returns
$G[n^{\prime},i,j]=G_{ij}(\mathrm{i}\omega^{\prime}_{n})$ as a three-
dimensional tensor.

The supplemental materials contain detailed instructions and usage examples
for the code.

##  IV Results

We showcase results for the temperature extrapolation technique applied to a
range of typical self-consistent finite-temperature methods. Sec. IV.1
investigates the accuracy and convergence properties of the extrapolation
technique. Sec. IV.2 examines the convergence acceleration resulting from an
enhanced starting point provided by the extrapolation technique. This study is
conducted specifically within the context of single-site DMFT and cluster DMFT
calculations for the Hubbard model, as well as the self-consistent GW
approximation for nickel oxide (NiO) real material calculations. Sec. IV.3
demonstrates the occurrence of hysteresis during heating and cooling processes
at a first-order phase transition, utilizing the extrapolation technique.

###  IV.1 Temperature extrapolation

Figure 2: Imaginary part of a cluster DMFT [7] self-energy,
$\operatorname{Im}\Sigma(\mathrm{i}\omega_{n})$, as a function of $\omega_{n}$
for the 2D Hubbard model at the antinodal point $(\pi,0)$ with $U=6t$ (see
text for calculation setup). Converged data at $\beta t=18$ (red crosses) and
$\beta t=20$ (black dots) along with temperature extrapolation (red circle)
from $\beta t=18$ to $\beta t=20$. The labels of the extrapolation results
indicate the initial inverse temperature used for extrapolation.  Figure 3:
Difference
$|\Delta\text{Im}\Sigma(\mathrm{i}\omega_{n})/\text{Im}\Sigma(\mathrm{i}\omega_{n})|$
between the converged self-energy at $\beta t=20$ and the extrapolation from
the converged self-energy at higher temperature, shown as a function of the
extrapolation inverse temperature for the three lowest Matsubara frequencies
$\omega_{n}$, with $n=0,1,2$ (see text for calculation setup).

In this section, we provide an illustration of temperature extrapolation in a
cluster DMFT [7] calculation, as depicted in Figure 2. Specifically, we
present the self-energy at the antinodal point, $(\pi,0)$, for a simulation of
the two-dimensional Hubbard model [48]. The simulation was performed on a
16-site ($4\times 4$) cluster, and the continuous-time auxiliary field
impurity solver [49, 50] was employed. The system was solved at half-filling,
in the paramagnetic state, with no next-nearest neighbor hopping, and at an
interaction strength of $U=6t$, which corresponds to a Mott insulating regime.
The self-energy in this regime exhibits a strong temperature dependence, and
various phases such as metallic, insulating, superconducting, and pseudo-
gapped phases are in close proximity [51, 31, 52]. For a comprehensive
discussion of the physics of this self-energy, we refer the reader to the
extensive literature on this system, including the reviews [7, 48] and
references therein.

Fig. 2 presents the fully converged, self-consistent imaginary part of the
self-energy, $\Sigma(\mathrm{i}\omega_{n})$, at $\beta t=18$ (red crosses) and
at $\beta t=20$ (black dots). Red circles illustrate the extrapolation of the
$\beta=18t$ self-energy to $\beta=20t$ using the Carathéodory formalism. The
extrapolated data aligns with the converged points for all frequency points,
except for the lowest one. The discrepancy between converged data and
extrapolation data at $\beta=20t$ serves as an indicator of additional
correlations emerging in the system as it cools down.

Fig. 3 demonstrates the deviation of the extrapolated self-energy at the three
lowest Matsubara frequencies at $\beta t=20$. The extrapolation is conducted
using converged data at $\beta t=10,11,...,19$. Notably, while the lowest
frequency of the self-energy exhibits significant deviations when extrapolated
by a factor of two from $\beta t=10$ to $\beta t=20$, these deviations rapidly
diminish if the extrapolation interval (the difference between the initial
temperature and targeted temperature) is reduced.

Figs. 2 and 3 demonstrate that the Carathéodory structure of many-body
functions can be effectively utilized to extrapolate data in temperature.

Figure 4: Cluster DMFT convergence of the imaginary part of the local Green’s
function $G_{\text{loc}}$ as a function of iteration for the three lowest
Matsubara frequencies at $\beta t=20$ (see text for calculation setup).
Convergence to the fixed point is shown for the non-interacting starting point
(filled symbols) and the starting point extrapolated from the converged
solution at $\beta t=18$ (open symbols). The labels of the extrapolation
results indicate the initial inverse temperature used for extrapolation.

###  IV.2 Starting point and Convergence

####  IV.2.1 Convergence of the extrapolated starting point away from phase
transitions

A temperature extrapolation of this type can be used to improve the starting
point for the convergence of self-consistent many-body methods, such as DMFT
[53, 54, 28] and its cluster extensions [4, 5, 6, 7]. The DMFT equations are
solved using a fixed-point iteration scheme, in which an initial estimate of
the self-energy (commonly chosen as $\Sigma(i\omega_{n})\equiv 0$) serves as
an initial guess for the iteration. The equations typically converge rapidly
if the initial starting point of the iteration is near the stationary point.
Since the computational effort is proportional to the number of iterations
required, a starting point closer to the stationary point reduces the
calculation time.

Fig. 4 displays the imaginary part of the lowest three frequencies of the
local Green’s function, $G_{\text{loc}}(\mathrm{i}\omega_{n})$, for a typical
cluster DMFT calculation away from phase transitions, using non-interacting
starting points (filled symbols) and starting points derived from temperature
extrapolation (open symbols), as a function of iteration. The calculation
setup is identical to the one in Sec. IV.1 (a 16-site cluster in the
paramagnetic state with $\beta t=20$, $U=6t$, and at half filling). A naive
non-interacting starting point converges in approximately seven iterations.
Moreover, a few additional iterations are necessary to confirm that
convergence has been achieved. In contrast, a starting point generated by
extrapolating from higher to lower temperatures almost immediately yields a
converged result.

####  IV.2.2 DMFT convergence near a phase transition

Figure 5: DMFT convergence of $-G(\beta/2)$ on a half-filled inifinite-
coordination Bethe lattice with $U=4.7t$ and bandwidth $W=4t$ as a function of
iteration. The top and lower panels display the insulating ($\beta t=20.31$)
and metallic ($\beta t=21.56$) phases, respectively. Atomic starting point:
gray dots. Non-interacting starting point: black dots. Extrapolated starting
points: red circles and blue circles. Inset: phase diagram adapted from Ref.
[55]. The labels of the extrapolation results indicate the initial inverse
temperature used for extrapolation.

Self-consistent many-body methods are notoriously slow to converge in the
vicinity of phase transitions. This can be illustrated at the example of a
single-site DMFT calculation on an infinite coordination number Bethe lattice
with bandwidth $W=4t$. The single-site DMFT calculation exhibits a first-order
phase transition between a paramagnetic metal at weak interaction and low
temperature and a paramagnetic Mott insulator at large interaction and higher
temperature [54, 56, 57, 58, 59, 60, 55], as depicted in the inset of Fig. 5.
At an interaction strength $U=4.7t$, the system is insulating at high
temperature, metallic at low temperature, and both metallic and insulating
solutions coexist in an intermediate temperature regime. We show the
convergence of the imaginary-time Green’s function in the middle of the
imaginary time interval, $-G(\beta/2)$, in Fig. 5, which is related to the
spectral density at the Fermi surface in the low-temperature limit as
$A(\omega=0)\sim-\beta G(\beta/2)$, at temperatures close to the two
boundaries of the coexistence region.

At a temperature of $\beta t=20.31$, the system exhibits Mott insulating
behavior. The upper panel of Fig. 5 shows that hundreds of iterations are
required to achieve convergence when initiating from the non-interacting limit
(black dots). Similarly, starting with an extrapolation from a lower
temperature metallic solution in the coexistence region (blue circles) results
in slow convergence toward the insulating fixed point. Conversely, the
extrapolation from a higher temperature insulating phase converges rapidly
(red circles). A starting point derived from the atomic limit also leads to
relatively fast convergence (gray dots).

Conversely, at a temperature of $\beta t=21.56$ (Fig. 5 lower panel), the
system is in a metallic state. Convergence from an ‘insulating’ starting point
requires hundreds of iterations to reach convergence (red circles for the
starting point extrapolated from a higher temperature insulating solution, and
gray dots for the atomic limit starting point). Convergence from a metallic
solution, like the non-interacting limit solution (black dots), is
substantially faster, and a starting point extrapolated from a converged
metallic solution at lower temperature leads to even faster convergence (blue
circles).

Fig. 5 therefore shows that even though all starting points, including atomic,
non-interacting, and extrapolated, converge to the identical solution,
extrapolation efficiently accelerates convergence, provided that the
extrapolated solution is in the same phase as the solution at the target
temperature. For extrapolation originating from a different phase (at the
boundary of phase transitions, solutions with distinct properties become
possible for two neighboring temperatures), the extrapolated starting points
still perform better than the simplistic starting points corresponding to the
‘wrong’ phase, specifically, the non-interacting solution or the atomic
solution.

Figure 6: Cluster DMFT convergence of the squared magnetic moment on the half-
filled 3D Hubbard model (using a 34-site cluster with a doubled unit cell) at
$U=6t$ and $\beta t=3.3$ (upper panel) and $3.5$ (lower panel) as a function
of iteration. Dots: non-interacting starting point. Circles: temperature
extrapolation. Inset: phase diagram adapted from Ref. [61]. $M^{2}$ at the
first and last iteration indicated by numbers. $\overline{M^{2}}$ corresponds
to the values at the last iteration of red or blue circles, which serves as a
reference for converged values. The labels of the extrapolation results
indicate the initial inverse temperature used for extrapolation.

####  IV.2.3 Cluster DMFT convergence near a phase transition

At second order phase transitions, critical slowing down, rather than
coexistence and hysteresis, is expected. This phenomenon leads to a slow
convergence of the fixed point iteration. We illustrate this at the example of
a half-filled three-dimensional Hubbard model treated within cluster DMFT [62,
63, 64], calculated on a 34-site cluster with a doubled unit cell. This model
exhibits a phase transition between a paramagnetic (PM) state at high
temperature and an antiferromagnetic (AFM) insulator at low temperature, as
depicted in the inset of Fig. 6. The model has been studied extensively within
the context of ultracold atomic gases. In Fig. 6, we examine two points,
$\beta t=3.3$ and $\beta t=3.5$, in the ordered phase at $U=6t$, near the AFM
phase transition (the critical point is at about $\beta t=2.7$). A starting
point from temperature extrapolation cannot completely overcome the slowing
down (since the self-energy itself exhibits strong temperature dependence) but
does lead to an acceleration of the convergence by at least a factor of two,
as illustrated by the convergence of the squared magnetic moment shown in Fig.
6.

Figure 7: GW convergence of the squared magnetic moment (left) and energy
(right) near the AFM phase transition of NiO at $\beta=30\;\mathrm{Ha}^{-1}$
(upper panel) and $\beta=31\;\mathrm{Ha}^{-1}$ (lower panel) as a function of
iteration. Calculations with unrestricted Hartree-Fock self-energy starting
point are represented by black dots, while calculations with temperature-
extrapolated self-energy are denoted by red circles (indicating extrapolation
from higher temperature solutions) or blue circles (indicating extrapolation
from lower temperature solutions). The values of $M^{2}$ and $E$ at the first
and last iterations are indicated by the numbers. The values of
$\overline{M^{2}}$ and $\overline{E}$ are determined as the values at the last
iteration of red or blue circles. The Hartree atomic units are used. The
labels of the extrapolation results indicate the initial inverse temperature
used for extrapolation. Figure 8: GW convergence of the squared magnetic
moment (left) and energy (right) deep within the AFM phase of NiO at
$\beta=100\;\mathrm{Ha}^{-1}$ (upper panel) and $\beta=200\;\mathrm{Ha}^{-1}$
(lower panel) as a function of iteration. Calculations with unrestricted
Hartree-Fock self-energy starting point are represented by black dots, while
calculations with temperature-extrapolated self-energy are denoted by red
circles (indicating extrapolation from higher temperature solutions) or blue
circles (indicating extrapolation from lower temperature solutions). The
values of $M^{2}$ and $E$ at the first and last iterations are indicated by
the numbers. The values of $\overline{M^{2}}$ and $\overline{E}$ are
determined as the values at the last iteration of red or blue circles. The
Hartree atomic units are used. The labels of the extrapolation results
indicate the initial inverse temperature used for extrapolation.

####  IV.2.4 GW convergence of realistic many-body simulations

Finally, we turn to realistic simulations within a weak-coupling framework. We
show examples for periodic solid NiO treated within the so-called GW
approximation [8, 9, 10, 11, 12, 13]. The GW approximation takes into account
screening processes via a renormalized, frequency-dependent interaction, but
neglects the second-order exchange diagram. It is therefore mostly used for
weakly correlated systems such as semiconductors. In the calculation for NiO,
a 4x4x4 cluster with a doubled unit cell along the $[1,1,1]$ direction and a
fixed lattice constant $a=4.1705$ Å [65] is utilized. The gth-dzvp-molopt-sr
basis [66] along with the gth-pbe pseudopotential [3] is employed. For the
density fitting of Coulomb integrals, the def2-svp-ri basis set is chosen as
the auxiliary basis [67]. Finite-size errors in the GW exchange diagram are
corrected using the Ewald probe-charge approach [68, 69]. The Coulomb
integrals and non-interacting matrix elements are obtained from PYSCF [70]. In
order to decrease the number of frequencies utilized in the computation, the
intermediate representation grid [36] is employed. A comprehensive description
of the methods and implementation, in conjunction with the computational setup
for evaluating NiO, is extensively detailed in Refs. [71, 12].

When applied to the antiferromagnetic materials NiO, the GW method shows a
continuous transition to an ordered state with a non-zero staggered
magnetization at low temperature. Within the GW framework, the transition
temperature is situated near $\beta=25\;\mathrm{Ha}^{-1}$. As illustrated in
Fig. 7, the convergence to this ordered state, initiated from an unrestricted
Hartree-Fock solution (black dots), is relatively slow. This is evident in the
convergence of both the squared magnetic moment for Ni and the total energy.
However, when starting from extrapolated starting points (blue or red
circles), the convergence occurs at a significantly faster pace, with an
initial starting point akin to iteration 30 of the Hartree-Fock-initialized
convergence.

At considerably lower temperatures, deep within the AFM phase, the
extrapolation of the self-energy results in a starting point that is
essentially exact because the self-energy is only weakly temperature-
dependent, similar to the 2D Hubbard model case shown at the beginning of Sec.
IV.2. Fig. 8 displays the temperature extrapolation from $\beta t=100$ to
$\beta t=200$, and vice versa, in comparison with the Hartree-Fock starting
points.

###  IV.3 Heating and Cooling

The possibility of obtaining a starting point for fixed-point iteration by
extrapolating from a nearby temperature point allows to smoothly change
temperature in subsequent simulation, in analogy to ‘heating’ and ‘cooling’
measurements in experiments. This capability is especially important at first-
order phase transitions, where multiple stable fixed points, corresponding to
the different phases, exist.

Such a heating and cooling process is shown in Fig. 9. Shown is a coexistence
region and hysteresis curve at the example of the single-site DMFT on a Bethe
lattice (the same calculation setup as the one mentioned in Sec. IV.2). We
plot the double occupancy, which is directly related to the potential energy,
as a function of inverse temperature $\beta$. With starting points
extrapolated from higher temperature converged results (red open circles), we
find a transition from an insulating phase with small double-occupancy to a
metallic phase with large double-occupancy around $\beta t=21.56$. Conversely,
with starting points extrapolated from lower temperature converged results
(blue open circles), we find a transition from a metallic phase to an
insulating phase around $\beta t=20.31$. Between those temperatures, both
metal and insulator are stable solutions of the fixed-point equations,
indicating phase coexistence [54, 56, 57, 58, 59, 60, 55, 72].

Figure 9: DMFT analysis of a half-filled Bethe lattice with $U=4.7t$ and a
bandwidth of $4t$, displaying the hysteresis near $\beta t=21$. Points on the
red curve utilize extrapolation from the nearest higher temperature converged
results as the starting points for the calculation, while points on the blue
curve employ extrapolation from the nearest lower temperature converged
results as the starting points for the calculation. Inset: phase diagram
adapted from Ref. [55].

##  V Conclusions

The Carathéodory temperature extrapolation technique, introduced in this
paper, provides an effective way to accelerate the convergence of self-
consistent many-body calculations for both model Hamiltonians and realistic
systems. It allows for smooth temperature variation in simulations, making it
suitable for studying heating and cooling processes in many-body systems. The
starting points provided by this method are generally superior, especially in
systems with phase transitions and convergence issues. The Carathéodory
technique therefore offers a versatile and efficient approach for studying
temperature-dependent properties in self-consistent many-body calculations
that should be adapted in any self-consistent finite-temperature many-body
simulation.

We note that ‘convergence acceleration’ techniques such as DIIS [73, 74, 75]
and Anderson acceleration [76, 77, 78, 79] are complementary to the
Carathéodory temperature extrapolation method employed in this study.
Convergence acceleration techniques work by extrapolating from a set of
initial iterations. The combination of such techniques with the better
starting points provided by Carathéodory temperature extrapolation is
therefore straightforward.

###### Acknowledgements.

We thank André Erpenbeck for helpful discussions. This material is based upon
work supported by the National Science Foundation under Grant No. NSF DMR
2001465. Real material calculations used resources of the National Energy
Research Scientific Computing Center, a DOE Office of Science User Facility
supported by the Office of Science of the U.S. Department of Energy under
Contract No. DE-AC02-05CH11231 using NERSC award BES-ERCAP0023196.

## References

  * Abrikosov  _et al._ [1975] A. A. Abrikosov, L. P. Gor’kov, and I. E. Dzi͡aloshinskiĭ,  _Methods of quantum field theory in statistical physics_ , rev. english ed ed. (Dover Publications, New York, 1975).
  * Mahan [2000] G. D. Mahan,  _Many-Particle Physics_ , 3rd ed., Physics of Solids and Liquids (Kluwer Academic/Plenum Publishers, New York, 2000).
  * Goedecker  _et al._ [1996] S. Goedecker, M. Teter, and J. Hutter, Physical Review B 54, 1703 (1996).
  * Hettler  _et al._ [1998] M. H. Hettler, A. N. Tahvildar-Zadeh, M. Jarrell, T. Pruschke, and H. R. Krishnamurthy, Physical Review B 58, R7475 (1998).
  * Lichtenstein and Katsnelson [2000] A. I. Lichtenstein and M. I. Katsnelson, Physical Review B 62, R9283 (2000).
  * Kotliar  _et al._ [2001] G. Kotliar, S. Y. Savrasov, G. Pálsson, and G. Biroli, Physical Review Letters 87, 186401 (2001).
  * Maier  _et al._ [2005] T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Reviews of Modern Physics 77, 1027 (2005).
  * Hedin [1965] L. Hedin, Physical Review 139, A796 (1965).
  * Aryasetiawan and Gunnarsson [1998] F. Aryasetiawan and O. Gunnarsson, Reports on Progress in Physics 61, 237 (1998).
  * Kutepov  _et al._ [2009] A. Kutepov, S. Y. Savrasov, and G. Kotliar, Physical Review B 80, 041103 (2009).
  * Kutepov [2020] A. L. Kutepov, Computer Physics Communications 257, 107502 (2020).
  * Yeh  _et al._ [2022a] C.-N. Yeh, S. Iskakov, D. Zgid, and E. Gull, Physical Review B 106, 235104 (2022a).
  * Yeh  _et al._ [2022b] C.-N. Yeh, A. Shee, Q. Sun, E. Gull, and D. Zgid, Physical Review B 106, 085121 (2022b).
  * Dahlen and van Leeuwen [2005] N. E. Dahlen and R. van Leeuwen, The Journal of Chemical Physics 122, 164102 (2005).
  * Phillips and Zgid [2014] J. J. Phillips and D. Zgid, The Journal of Chemical Physics 140, 241101 (2014).
  * Rusakov and Zgid [2016] A. A. Rusakov and D. Zgid, The Journal of Chemical Physics 144, 054106 (2016).
  * Iskakov  _et al._ [2019] S. Iskakov, A. A. Rusakov, D. Zgid, and E. Gull, Physical Review B 100, 085112 (2019).
  * Bickers  _et al._ [1989] N. E. Bickers, D. J. Scalapino, and S. R. White, Physical Review Letters 62, 961 (1989).
  * Iskakov and Gull [2022] S. Iskakov and E. Gull, Physical Review B 105, 045109 (2022).
  * Prokof’ev and Svistunov [2007] N. Prokof’ev and B. Svistunov, Physical Review Letters 99, 250201 (2007).
  * Gull  _et al._ [2010a] E. Gull, D. R. Reichman, and A. J. Millis, Physical Review B 82, 075109 (2010a).
  * Gull  _et al._ [2011a] E. Gull, D. R. Reichman, and A. J. Millis, Physical Review B 84, 085134 (2011a).
  * Cohen  _et al._ [2014] G. Cohen, D. R. Reichman, A. J. Millis, and E. Gull, Physical Review B 89, 115139 (2014).
  * Biermann  _et al._ [2003] S. Biermann, F. Aryasetiawan, and A. Georges, Physical Review Letters 90, 086402 (2003).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Reviews of Modern Physics 78, 865 (2006).
  * Zgid and Gull [2017] D. Zgid and E. Gull, New Journal of Physics 19, 023047 (2017).
  * Kozik  _et al._ [2015] E. Kozik, M. Ferrero, and A. Georges, Physical Review Letters 114, 156402 (2015).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Reviews of Modern Physics 68, 13 (1996).
  * Kotliar and Vollhardt [2004] G. Kotliar and D. Vollhardt, Physics Today 57, 53 (2004).
  * Macridin  _et al._ [2006] A. Macridin, M. Jarrell, and Th. Maier, Physical Review B 74, 085104 (2006).
  * Gull  _et al._ [2009] E. Gull, O. Parcollet, P. Werner, and A. J. Millis, Physical Review B 80, 245102 (2009).
  * Walsh  _et al._ [2019] C. Walsh, P. Sémon, D. Poulin, G. Sordi, and A.-M. S. Tremblay, Physical Review Letters 122, 067203 (2019).
  * Fei  _et al._ [2021a] J. Fei, C.-N. Yeh, D. Zgid, and E. Gull, Physical Review B 104, 165111 (2021a).
  * Stanescu and Kotliar [2006] T. D. Stanescu and G. Kotliar, Physical Review B 74, 125110 (2006).
  * Gull  _et al._ [2018] E. Gull, S. Iskakov, I. Krivenko, A. A. Rusakov, and D. Zgid, Physical Review B 98, 075127 (2018).
  * Shinaoka  _et al._ [2017] H. Shinaoka, J. Otsuki, M. Ohzeki, and K. Yoshimi, Physical Review B 96, 035147 (2017).
  * Boehnke  _et al._ [2011] L. Boehnke, H. Hafermann, M. Ferrero, F. Lechermann, and O. Parcollet, Physical Review B 84, 075145 (2011).
  * Kananenka  _et al._ [2016] A. A. Kananenka, A. R. Welden, T. N. Lan, E. Gull, and D. Zgid, Journal of Chemical Theory and Computation 12, 2250 (2016).
  * Kaye  _et al._ [2022] J. Kaye, K. Chen, and O. Parcollet, Physical Review B 105, 235115 (2022).
  * Li  _et al._ [2020] J. Li, M. Wallerberger, N. Chikano, C.-N. Yeh, E. Gull, and H. Shinaoka, Physical Review B 101, 035144 (2020).
  * Kaltak and Kresse [2020] M. Kaltak and G. Kresse, Physical Review B 101, 205145 (2020).
  * Fei [2021] J. Fei,  _A Probe into Propagators_ , Thesis (2021).
  * Nogaki and Shinaoka [2023] K. Nogaki and H. Shinaoka, Journal of the Physical Society of Japan 92, 035001 (2023).
  * Cayley [1846] A. Cayley,  1846, 119 (1846).
  * Schur [1918] J. Schur, Journal für die reine und angewandte Mathematik 148, 122 (1918).
  * Note [1] See Supplemental Material at [URL will be inserted by publisher] for Python implementation of the algorithm, documentation, and examples.
  * Fei  _et al._ [2021b] J. Fei, C.-N. Yeh, and E. Gull, Physical Review Letters 126, 056402 (2021b).
  * Qin  _et al._ [2022] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, Annual Review of Condensed Matter Physics 13, 275 (2022).
  * Gull  _et al._ [2008] E. Gull, P. Werner, O. Parcollet, and M. Troyer, Europhysics Letters 82, 57003 (2008).
  * Gull  _et al._ [2011b] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Reviews of Modern Physics 83, 349 (2011b).
  * Werner  _et al._ [2009] P. Werner, E. Gull, O. Parcollet, and A. J. Millis, Physical Review B 80, 045120 (2009).
  * Gull  _et al._ [2010b] E. Gull, M. Ferrero, O. Parcollet, A. Georges, and A. J. Millis, Physical Review B 82, 155101 (2010b).
  * Metzner and Vollhardt [1989] W. Metzner and D. Vollhardt, Physical Review Letters 62, 324 (1989).
  * Georges and Kotliar [1992] A. Georges and G. Kotliar, Physical Review B 45, 6479 (1992).
  * Blümer [2003] N. Blümer,  _Mott-Hubbard metal insulator transition and optical conductivity in high dimensions_ , Ph.D. thesis, Universität Augsburg (2003).
  * Jarrell [1992] M. Jarrell, Physical Review Letters 69, 168 (1992).
  * Rozenberg  _et al._ [1992] M. J. Rozenberg, X. Y. Zhang, and G. Kotliar, Physical Review Letters 69, 1236 (1992).
  * Zhang  _et al._ [1993] X. Y. Zhang, M. J. Rozenberg, and G. Kotliar, Physical Review Letters 70, 1666 (1993).
  * Rozenberg  _et al._ [1999] M. J. Rozenberg, R. Chitra, and G. Kotliar, Physical Review Letters 83, 3498 (1999).
  * Bulla [1999] R. Bulla, Physical Review Letters 83, 136 (1999).
  * Staudt  _et al._ [2000] R. Staudt, M. Dzierzawa, and A. Muramatsu, The European Physical Journal B - Condensed Matter and Complex Systems 17, 411 (2000).
  * Kent  _et al._ [2005] P. R. C. Kent, M. Jarrell, T. A. Maier, and Th. Pruschke, Physical Review B 72, 060411 (2005).
  * Fuchs  _et al._ [2011a] S. Fuchs, E. Gull, L. Pollet, E. Burovski, E. Kozik, T. Pruschke, and M. Troyer, Physical Review Letters 106, 030401 (2011a).
  * Fuchs  _et al._ [2011b] S. Fuchs, E. Gull, M. Troyer, M. Jarrell, and T. Pruschke, Physical Review B 83, 235113 (2011b).
  * Bartel and Morosin [1971] L. C. Bartel and B. Morosin, Physical Review B 3, 1039 (1971).
  * VandeVondele and Hutter [2007] J. VandeVondele and J. Hutter, The Journal of Chemical Physics 127, 114105 (2007).
  * Hättig [2005] C. Hättig, Physical Chemistry Chemical Physics 7, 59 (2005).
  * Paier  _et al._ [2005] J. Paier, R. Hirschl, M. Marsman, and G. Kresse, The Journal of Chemical Physics 122, 234102 (2005).
  * Sundararaman and Arias [2013] R. Sundararaman and T. A. Arias, Physical Review B 87, 165122 (2013).
  * Sun  _et al._ [2018] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, WIREs Computational Molecular Science 8, e1340 (2018).
  * Iskakov  _et al._ [2020] S. Iskakov, C.-N. Yeh, E. Gull, and D. Zgid, Physical Review B 102, 085105 (2020).
  * Erpenbeck  _et al._ [2023] A. Erpenbeck, W.-T. Lin, T. Blommel, L. Zhang, S. Iskakov, L. Bernheimer, Y. Núñez-Fernández, G. Cohen, O. Parcollet, X. Waintal, and E. Gull, A Tensor Train Continuous Time Solver for Quantum Impurity Models (2023), arxiv:2303.11199 [cond-mat] .
  * Pulay [1980] P. Pulay, Chemical Physics Letters 73, 393 (1980).
  * Pulay [1982] P. Pulay, Journal of Computational Chemistry 3, 556 (1982).
  * Pokhilko  _et al._ [2022] P. Pokhilko, C.-N. Yeh, and D. Zgid, The Journal of Chemical Physics 156, 094101 (2022).
  * Bacskay [1981] G. B. Bacskay, Chemical Physics 61, 385 (1981).
  * Walker and Ni [2011] H. F. Walker and P. Ni, SIAM Journal on Numerical Analysis 49, 1715 (2011).
  * Evans  _et al._ [2020] C. Evans, S. Pollock, L. G. Rebholz, and M. Xiao, SIAM Journal on Numerical Analysis 58, 788 (2020).
  * Pollock and Rebholz [2021] S. Pollock and L. G. Rebholz, IMA Journal of Numerical Analysis 41, 2841 (2021).