# On the Positive Definiteness of Response Functions in the Time Domain

Alexander F. Kemper<EMAIL_ADDRESS>Department of Physics, North Carolina
State University, Raleigh, North Carolina 27695, USA  Chao Yang<EMAIL_ADDRESS>Computational Research Division, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA  Emanuel Gull<EMAIL_ADDRESS>Department of Physics,
University of Michigan, Ann Arbor, Michigan 48109, USA

###### Abstract

Response functions of quantum systems, such as electron Green’s functions,
magnetic, or charge susceptibilities, describe the response of a system to an
external perturbation. They are the central objects of interest in field
theories and quantum computing and measured directly in experiment. Response
functions are intrinsically causal. In equilibrium and steady-state systems,
they correspond to a positive spectral function in the frequency domain. This
article shows that response functions define an inner product on a Hilbert
space and thereby induce a positive definite function. The properties of this
function can be used to reduce noise in measured data and, in equilibrium and
steady state, to construct positive definite extensions for data known on
finite time intervals, which are then guaranteed to correspond to positive
spectra.

_Introduction._ Response functions are critical for the understanding of
physics in a wide variety of contexts. Given a state, they are the natural
framework for considering the properties of said state, and in particular the
_dynamical_ properties that involve excitation out of the state[1, 2, 3, 4].
Response functions are also what is measured experimentally, ranging from low-
frequency THz conductivity to magnetic susceptibilities to photoemission
spectroscopy. A large body of literature has been devoted to the study of
response functions, and a number of “field” theoretical approaches avoid the
calculation of eigenstates entirely, and instead cast the formalism entirely
in terms of response functions — called correlation functions in this context
— typically the electronic Green’s function. These include embedding
techniques such as the dynamical mean field theory (DMFT) [5, 6, 7] and its
cluster extensions [8, 9, 10, 11], Monte Carlo approaches for lattice [12] and
impurity [13, 14, 15, 16, 17] models, self-consistent partial summation
methods for real materials and model systems [18, 19, 20], and non-equilibrium
Green’s function methods [4, 21, 22, 23, 24, 25, 26, 27, 28].

An important characteristic of correlation functions lies in their analytical
properties. The retarded correlation functions have no content at negative
times due to causality[1, 2, 3]. In the complex frequency domain, they
correspond to so-called Nevanlinna functions [29] whose poles are restricted
to the lower half of the complex plane. This analytical framework has long
been utilized to evaluate integrals that emerge in many-body theory [1, 2, 3].
More recently, it was used in order to perform analytic continuation from a
Wick-rotated frame to a standard frame. Typical approaches that rely e.g. on
the maximum entropy method [30] have significant uncertainty in the final
result, resulting in washed out spectra. However, explicitly enforcing the
highly constraining analytical properties of Nevanlinna functions results in a
drastic reduction of the uncertainty, leading to sharp spectral functions [29,
31]. It is clear from these examples that encoding the mathematical structure
of correlation functions into numerical algorithms can be used to great
benefit.

In this Letter, we derive and analyze a fundamental property of correlation
functions in the time domain: several correlation functions of interest are
_positive definite functions_ of their time arguments. This property arises
from viewing correlation functions as an inner product in the vector space of
operators, combined with the fact that the time translation operator is a
unitary representation of the time translation group. This positive
definiteness sets a strong constraint on the correlation function, similar to
its other analytic properties — in fact, some of these directly follow from
the positive definiteness[32]. Moreover, as we will show, insisting that a
correlation function is positive definite enables both the extension of
numerical correlation functions to later times, and the extraction of clean
spectra from noisy data such as that obtained from Monte Carlo and Quantum
Computing approaches.

We consider a quantum system in a quantum state $\rho$. In a second quantized
formalism, quantum states can be described as states in Fock space [1, 2, 3].
The energetics (and dynamics) of the system is described by a Hamiltonian
$\mathcal{H}$; in the Heisenberg picture, the time-evolution of an operator
$A$ is given by $A(t)=e^{i\mathcal{H}t}Ae^{-i\mathcal{H}t}$. Expectation
values of operators for the state $\rho$ are computed as $\langle
A(t)\rangle=\text{Tr}[\rho A(t)]$. In a system with time-translational
invariance, $\mathcal{H}$ commutes with $\rho$, and in thermal equilibrium
$\rho=e^{-\beta\mathcal{H}}/\mathcal{Z}$, with $\beta$ denoting the inverse
temperature and $\mathcal{Z}\equiv\mathrm{Tr}e^{-\beta\mathcal{H}}$.

We are primarily interested in time-dependent correlation functions of the
type

$\displaystyle G_{AB}(t,t^{\prime})=\langle
A^{\dagger}(t)B(t^{\prime})\rangle.$ (1)

In a system with creation operators $c_{i}^{\dagger}$ creating particles in
state $i$, the so-called “lesser” Green’s function [4]
$G_{ij}^{<}(t,t^{\prime})=i\text{Tr}\left[\rho
c_{j}^{\dagger}(t^{\prime})c_{i}(t)\right]$ and the “greater” Green’s function
$G_{ij}^{>}(t,t^{\prime})=-i\text{Tr}\left[\rho
c_{i}(t)c_{j}^{\dagger}(t^{\prime})\right]$, as well as the corresponding
charge and spin correlation functions are related to correlation functions of
this type. Below, we will consider these functions without their usual $\pm i$
prefactors.

_Mathematical exposition._ In this section, we will show that
$G_{AB}(t,t^{\prime})$ are _positive definite_ functions when $A=B$. First, we
observe that the correlation function can be viewed as an inner product of the
linear operators that act on Fock space $V$ (we denote this vector space as
$\mathcal{L}(V)$). That is, we define an inner product on $\mathcal{L}(V)$ as

$\displaystyle\langle A,B\rangle:=\mathrm{Tr}\left[\rho A^{\dagger}B\right],$
(2)

which is conjugate symmetric, linear in the second component, and positive for
$\langle A,A\rangle$ when $A\neq 0$. The properties of conjugate symmetry and
linearity are straightforward to verify. The positivity property $\langle
A,A\rangle\geq 0$ can be established by noting that (i) if a square matrix $M$
equals the multiplication of a matrix with its Hermitian transpose, i.e.,
$M=AA^{\dagger}$, then $M$ is a Hermitian positive semi-definite matrix; and
(ii) the trace of the product of two positive semi-definite matrices $M$ and
$N$ is always greater than or equal to zero, i.e., $\text{Tr}(MN)\geq 0$. We
further restrict our consideration to the cases where $A,B$ are not orthogonal
to $\rho$, which would yield a zero expectation value and is thus not
physically relevant for correlation functions. $\mathcal{L}(V)$ together with
the inner product of Eq. 2 forms a Hilbert space.

The correlation functions are defined as complex-valued two-time functions
($\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{C}$):

$\displaystyle G_{AA}(t,t^{\prime})=\mathrm{Tr}\left[\rho
A(t)^{\dagger}A(t^{\prime})\right].$ (3)

Given the inner product, we can show that the Green’s functions arising in the
study of equilibrium and non-equilibrium dynamics of many-body physics [1, 2,
3, 4] are positive definite functions [33] of two variables. That is, given
any set of arbitrarily spaced time points $t$, the eigenvalues of the matrix
that result from evaluating the Green’s function at those points,
$\left[G_{AA}(t_{i},t_{j})\right]$, are all positive or zero 111Note that in
the standard terminology that positive definite functions include zero
eigenvalues [33].

We prove this by using an equivalent definition of positive definiteness,
which is that, given the aforementioned set of time points $t$,

$\displaystyle\sum_{ij}G_{AA}(t_{i},t_{j})$
$\displaystyle\lambda_{i}^{*}\lambda_{j}\geq 0$ (4)

for any set of $\lambda_{i}$. This follows from

$\displaystyle\sum_{ij}$ $\displaystyle
G_{AA}(t_{i},t_{j})\lambda_{i}^{*}\lambda_{j}=\sum_{i,j}\mathrm{Tr}\left[\rho
A(t_{i})^{\dagger}A(t_{j})\right]\lambda_{i}^{*}\lambda_{j}$ (5)
$\displaystyle=\mathrm{Tr}\left[\rho\left(\sum_{i}\lambda_{i}A(t_{i})\right)^{\dagger}\left(\sum_{j}\lambda_{j}A(t_{j})\right)\right]$
(6)
$\displaystyle=\bigg{\langle}\sum_{i}\lambda_{i}A(t_{i}),\sum_{j}\lambda_{j}A(t_{j})\bigg{\rangle}\geq
0.$ (7)

The correlation functions (and Green’s function as well as self-energies [35,
36]) are therefore positive definite functions. This holds for any operator
$A\in\mathcal{L}(V)$, and in particular it holds for fermionic
creation/annihilation operators $c^{\dagger}/c$, for densities $n$, for
magnetization operators $n_{\uparrow}-n_{\downarrow}$, and trivially for the
identity. In multi-orbital systems, the diagonal components of electronic
Green’s functions $G^{\lessgtr}_{ii}$ are positive semi-definite. Off-diagonal
components can be constructed from linear combinations of $\langle
c_{i}+c_{j},c_{i}+c_{j}\rangle$, $\langle c_{i}+ic_{j},c_{i}+ic_{j}\rangle$,
and the diagonal components, all of which are positive definite.

In the presence of time-translation invariance, (i.e. for steady-state and
equilibrium systems), the Green’s function becomes a function of a single time
argument corresponding to the time difference: $G_{AA}(t-t^{\prime})$, i.e. a
positive definite function of a single variable where

$\displaystyle\sum_{ij}G_{AA}(t_{i}-t_{j})$
$\displaystyle\lambda_{i}^{*}\lambda_{j}\geq 0.$ (8)

In order for time translation to hold, the set of time points $t$ must form a
group under addition such as $\mathbb{R}$ or $\mathbb{Z}$. With this
additional structure, and as long as $\rho$ commutes with the time evolution
operator,

$\displaystyle\langle A(t),B(t^{\prime})\rangle=\langle
A(t-t^{\prime}),B\rangle=\langle A,B(t^{\prime}-t)\rangle$ (9)

follows from the cyclicity of the trace. For a more mathematical discussion of
the above points, see the supplementary information).

The positive definiteness of the Green’s function sets a strong constraint on
the function and has important practical consequences, which we will explore
and exploit in the remainder of the paper. In short, we will show that (i) we
can improve the signal-to-noise ratio in noisy data from experiment and
theory, and (ii) we can construct causal extensions (with a positive spectrum)
from short time data.

Figure 1: Real part of the on-site Hubbard dimer Green’s function
$G^{>}_{11\uparrow}(t)$. Red triangles: Data polluted by Gaussian noise with
$\sigma=0.1$ (left panel), $\sigma=0.05$ (middle), and $\sigma=0.01$ (right).
Green filled circles: Projection to nearest causal function. Black circles:
Exact solution. Other parameters are: temperature $T=1/10$, hopping $t=1$,
interaction $U=5$, level energy $\epsilon=2.3$.

_Denoising correlation functions._ A first, and natural application of the
mathematics presented above, is to take a correlation function from a source
that has some inherent noise, and to use the positive definite property to
project its values to the nearest positive definite correlation function; in
effect, denoising the data. Noisy correlation functions can arise from Monte
Carlo evaluations, from experimental measurements, or from simulations on
quantum computers.

We consider a discretized, time-translation invariant correlation function on
a regularly spaced time axis $t$, $G_{AA}(t_{i}-t_{j})$. If we label the
elements of the correlation function as $G_{AA}(t_{i}-t_{j})\rightarrow
f_{i-j}$, then the resulting matrix $\underline{G}$ with entries
$\underline{G}_{ij}=f_{i-j}$ is a positive semi-definite Hermitian Toeplitz
matrix [33],

$\displaystyle\underline{G}$
$\displaystyle=\begin{pmatrix}f_{0}&f_{1}&f_{2}&\cdots&f_{n}\\\
f_{1}^{*}&f_{0}&f_{1}&\cdots&f_{n-1}\\\
f_{2}^{*}&f_{1}^{*}&f_{0}&\cdots&f_{n-2}\\\ \vdots&&&\ddots&\vdots\\\
f_{n}^{*}&f_{n-1}^{*}&f_{n-2}^{*}&\cdots&f_{0}\end{pmatrix}.$ (10)

$\underline{G}$ is commonly known as the Gramian or “Gram” matrix.

In the presence of noise in the correlation function data ($f_{j}$), this
matrix is not positive semi-definite (PSD). However, alternating projections
to the nearest PSD matrix [37] followed by projections to the nearest Toeplitz
matrix [38] and an enforcement of the value at time zero (which is typically
known precisely, e.g. from equal-time measurements or sum rules) results in
quick convergence to a positive definite function. (Note that more elaborate
schemes, see e.g. [39, 40] and [41], may be substantially more efficient; we
have not tried them here). We illustrate this in Fig. 1 for a positive
definite Green’s function of the Hubbard dimer which has been polluted with
Gaussian noise (for details see supplement). Enforcing positive definiteness
by the procedure above shows a dramatic improvement of the quality of the
data. A second application to state-of-the-art quantum Monte Carlo data with a
continuous spectrum [28] is presented in the supplement.

_Extending correlation functions to longer times._ A second consequence of
positive definiteness is that positive definite extensions of response
functions to longer times exist, from which spectral measures can be
constructed. We make use of two well-known mathematical facts for positive
definite functions. First, functions that are only known on a subset of their
domain (e.g. $\mathbb{Z_{+}}$ or $\mathbb{R}$) have at least one extension to
the full domain [42, 43, 44]. This is a consequence of the extension theorems
of Kreĭn (in the case of $\mathbb{R}$) and Carathéodory (in the case of
$\mathbb{Z}$). Second, according to Bochner’s theorem [32], the Fourier
transform of a positive definite function (over $\mathbb{Z}$ or $\mathbb{R}$)
is guaranteed to have positive real part; moreover, the inverse Fourier
transform of a positive spectral function is guaranteed to be positive
definite. In the discrete case, these extensions may not be unique; in the
continuum case, uniqueness is guaranteed by the analyticity of the Green’s
function.

Conceptually, this realization offers a straightforward methodology to obtain
causal spectral functions from finite-time data: in a first step, the function
is extended from short to long times by taking advantage of an extension
theorem. In a second step, the uniquely defined and positive Fourier transform
of the extended data is computed.

Numerically, we can use extension theorems to predict Green’s function values
at later times from known values at short times [33]. Supplementing a known
time series $f_{0},\cdots,f_{n}$ by a single unknown element $f_{n+1}$, the
Gramian of Eq. 10 is a Toeplitz matrix with a single unknown complex number
$f_{n+1}=\underline{G}_{0,n+1}=(\underline{G}_{0,n+1})^{*}$. The real and
imaginary parts of $f_{n+1}$ can then be obtained in a two-dimensional search
in the complex plane for the region where the lowest eigenvalue of $G$ is zero
or larger, i.e. $G$ is positive definite; such a value must exist [42] and
$|f_{n+1}|\leq f_{0}$ due to the positive definiteness 222We thank Z. Sasvári
for suggesting this approach. In practice, we find that in systems with only
few well-defined excitations, and consequentially a low-rank Gram matrix, the
search converges quickly to a unique solution, whereas systems with continuous
spectra may present ambiguities in the extensions. Fig. 2 demonstrates the
approach for the Green’s function of the Hubbard dimer (for details see
supplement). Once the data is extended to all times, it yields a unique
spectral function [32].

Figure 2: Extension of the on-site Hubbard dimer Green’s function
$G^{>}_{11\uparrow}(t)$. Shown are real (black) and imaginary (red) parts of
input data up to $t=2$ (circles), the numerically computed extension from
$t=2$ to $t=10$ (squares), and the analytically known Green’s function up to
$t=10$ (lines).

More generally, given that the Gram matrix Eq. 10 is a Toeplitz matrix of size
$n$, and assuming a rank of size $r$, a classical result by Carathéodory and
Fejér [46] guarantees the existence of the decomposition $T=APA^{\dagger}$,
where $A$ is a $n\times r$ Vandermonde matrix and $P$ is a $r\times r$
positive diagonal matrix. The columns of $A$ can be interpreted as uniformly
sampled oscillation frequencies, and the entries of $P$ as positive ‘pole
strengths’. This approach therefore gives both access to the spectral function
in terms of a series of $r$ discrete poles at given frequencies, and a method
to construct extensions for all times by enlarging the matrix $A$ with
additional frequencies. The decomposition also establishes a connection to
signal processing and control theory, where PSD covariance matrices of time-
invariant processes are analyzed using this decomposition in super-resolution
algorithms such as the multiple signal classification (MUSIC) [47], and where
thereby approximate extensions from noisy data can be constructed. We will
explore this connection, together with connections to reproducing kernel
Hilbert spaces [48], in a future paper.

_Application to Quantum Computing._ As a final illustration, we apply the
denoising and extension implications of Eq. 8 to a noisy correlation function
measured on IBM’s quantum computer ibm_auckland. We have measured the
momentum-space greater Green’s function for the empty state $G_{k}(t):=\langle
0|c^{\dagger}_{k}(t)c_{k}|0\rangle$ for an 8-site Su-Schrieffer-Heeger model,
which is a model for free electrons with a hopping parameter that alternates
with an amplitude $\delta$

$\mathcal{H}=-V_{nn}\sum_{i}\left[1+(-\delta)^{i}\right]c^{\dagger}_{i}c_{i+1}+\mathrm{H.c}-\mu\sum_{i}c^{\dagger}_{i}c_{i},$
(11)

where we set $V_{nn}=1$. The raw data (originally published in Ref. 49), is
shown in Fig. 3, in the gapped phase with $\delta=0.4$ at $k=\pi/2$. We report
details of the calculation in the Supplementary information. Because the model
is translationally invariant and has 2 bands, the momentum basis Green’s
functions should exhibit at most 2 frequencies. However, due to the hardware
noise from the quantum computer, there is significant noise in the time domain
signal. This is similarly reflected in the Fourier transform, where peaks can
be identified, but only in the power spectrum — the spectral function (shown
inFig. 3) does not show the requisite analytic structure (a consistent sign
across all frequencies). We de-noise the data by asserting the positive
definiteness and performing a point-by-point optimization on the Green’s
function, where the cost function is the square of the negative eigenvalues,
and obtain an improved spectrum. As can be seen in the Figure, much of the
broad spectrum noise has disappeared, the peaks are clearly visible in the
spectral function, and it is starting to approach analyticity. We detail the
procedure and show the result on additional data from different $k$ points in
the supplement; we find that this iterative procedure is more suitable in the
presence of extreme noise than the alternating projection scheme used for
Monte Carlo data. The final step is to extend the data as discussed above,
which results in a long positive definite signal, with sharp peaks in the
spectrum and almost no non-analytic part.

Note that applying a PSD projection as done here is conceptually distinct from
recent efforts using advanced signals processing approaches to extract a
cleaner signal from noisy quantum simulation [50, 51, 52], and in fact can be
used as a complementary tool.

Figure 3: Main panel: greater Green’s function $G_{k}^{>}(t)$ at $k=\pi/2$ of
an 8-site Su-Schieffer-Heeger with $\delta=0.4$ and $\mu=-3$. The figure shows
the raw data obtained from the ibm_auckland quantum computer (blue), as well
as the PSD denoised data (orange), and PSD denoised and extended data (green).
Inset: Corresponding spectral functions. The Fourier transform used a damping
factor $\tau=100$. The vertical dashed lines indicate the expected analytic
frequencies.

_Discussion._ Positive semidefinite response functions are ubiquitous in many-
body physics[1, 4, 53], as is the desire to reduce systematic or stochastic
noise and to obtain corresponding spectral functions. The theory and
algorithms presented here are broadly applicable to problems ranging from the
analysis of experimental measurements to simulations of quantum systems on
classical and quantum hardware.

In this Letter, we have demonstrated the efficiency of removing noise with a
PSD projection in the examples of synthetic data and real-world quantum
computing data. We believe that a PSD noise filter should be applied much more
broadly to any positive definite response function, simulated or measured,
before data is analyzed and/or spectra are computed, including data from
quantum Monte Carlo, tensor networks, and time-resolved experiment.

While the examples in this work focused on Green’s functions, a PSD projection
is not limited to this particular use case since the operators in Eq. 4 are
rather general. For example, the identity operator is included, and thus
quantities of the form $\langle\psi_{0}|e^{-i\mathcal{H}t}|\psi_{0}\rangle$
are also amenable to denoising in the same way. This covers a range of novel
quantum algorithms that make use of such quantities, including subspace
expansions [54, 55] and energy minimization based on quantum phase estimation
[56, *ding2023simultaneous], as well as Loschmidt echo measurements [58].

We note that, similar to the extension techniques presented in this letter,
memory kernel approaches [59, 60, 61, 62, 63, 64] pioneered in transport
physics aim to extend data measured for short times to longer times. The
approach presented is complementary to those techniques and may be combined
with them.

The connection highlighted here between positive definite response functions
and the Toeplitz Vandermonde decomposition of the Gram matrix in addition
allows to make a connection to algorithms developed for signal processing and
control theory, such as MUSIC [47], which are used to analyze positive
semidefinite covariance matrices. This analysis will be the topic of a future
paper.

###### Acknowledgements.

We thank Vadym Adamyan for making us aware of extension theorems, Zoltan
Sasvári for pointing out the possibility of considering a Gram matrix extended
by a single point, Efekan Kökcu for suggestions regarding non-equilibrium
Green’s functions, and Michael Ragone for feedback on the mathematical
exposition. AFK was supported by the Department of Energy, Office of Basic
Energy Sciences, Division of Materials Sciences and Engineering under grant
no. DE-SC0023231. EG and CY were supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research and Office
of Basic Energy Science, Scientific Discovery through Advanced Computing
(SciDAC) program under Award Number(s) DE-SC0022088. The data for the figures
are submitted as a supplement to this paper.

## References

  * Alekseevich  _et al._ [1975] A. Alekseevich, L. P. Gorkov, and I. E. Dzyaloshinski,  _Methods of quantum field theory in statistical physics_ , Dover Books on Physics Series (Dover Publications, 1975).
  * Mahan [2010] G. D. Mahan,  _Many Particle Physics_ (Springer, New York, NY10013, USA, 2010).
  * Bruus and Flensberg [2004] H. Bruus and K. Flensberg,  _Many-body quantum theory in condensed matter physics: an introduction_ (OUP Oxford, 2004).
  * Stefanucci and van Leeuwen [2013] G. Stefanucci and R. van Leeuwen,  _Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction_ , 1st ed. (Cambridge University Press, 2013).
  * Metzner and Vollhardt [1989] W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
  * Georges and Kotliar [1992] A. Georges and G. Kotliar, Phys. Rev. B 45, 6479 (1992).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Reviews of Modern Physics 68, 13 (1996).
  * Hettler  _et al._ [1998] M. H. Hettler, A. N. Tahvildar-Zadeh, M. Jarrell, T. Pruschke, and H. R. Krishnamurthy, Phys. Rev. B 58, R7475 (1998).
  * Lichtenstein and Katsnelson [2000] A. I. Lichtenstein and M. I. Katsnelson, Phys. Rev. B 62, R9283 (2000).
  * Kotliar  _et al._ [2001] G. Kotliar, S. Y. Savrasov, G. Pálsson, and G. Biroli, Phys. Rev. Lett. 87, 186401 (2001).
  * Maier  _et al._ [2005] T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Reviews of Modern Physics 77, 1027 (2005).
  * Blankenbecler  _et al._ [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
  * Hirsch and Fye [1986] J. E. Hirsch and R. M. Fye, Phys. Rev. Lett. 56, 2521 (1986).
  * Rubtsov  _et al._ [2005] A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein, Phys. Rev. B 72, 035122 (2005).
  * Werner  _et al._ [2006] P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006).
  * Gull  _et al._ [2008] E. Gull, P. Werner, O. Parcollet, and M. Troyer, Europhysics Letters 82, 57003 (2008).
  * Gull  _et al._ [2011a] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011a).
  * Hedin [1965] L. Hedin, Phys. Rev. 139, A796 (1965).
  * Bickers  _et al._ [1989] N. E. Bickers, D. J. Scalapino, and S. R. White, Phys. Rev. Lett. 62, 961 (1989).
  * Yeh  _et al._ [2022] C.-N. Yeh, S. Iskakov, D. Zgid, and E. Gull, Phys. Rev. B 106, 235104 (2022).
  * Mühlbacher and Rabani [2008] L. Mühlbacher and E. Rabani, Phys. Rev. Lett. 100, 176403 (2008).
  * Werner  _et al._ [2009] P. Werner, T. Oka, and A. J. Millis, Phys. Rev. B 79, 035320 (2009).
  * Gull  _et al._ [2011b] E. Gull, D. R. Reichman, and A. J. Millis, Phys. Rev. B 84, 085134 (2011b).
  * Aoki  _et al._ [2014] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Rev. Mod. Phys. 86, 779 (2014).
  * Schüler  _et al._ [2020] M. Schüler, D. Golež, Y. Murakami, N. Bittner, A. Herrmann, H. U. Strand, P. Werner, and M. Eckstein, Computer Physics Communications 257, 107484 (2020).
  * Dong  _et al._ [2022] X. Dong, E. Gull, and H. U. R. Strand, Phys. Rev. B 106, 125153 (2022).
  * Núñez Fernández  _et al._ [2022] Y. Núñez Fernández, M. Jeannin, P. T. Dumitrescu, T. Kloss, J. Kaye, O. Parcollet, and X. Waintal, Phys. Rev. X 12, 041018 (2022).
  * Erpenbeck  _et al._ [2023] A. Erpenbeck, E. Gull, and G. Cohen, Phys. Rev. Lett. 130, 186301 (2023).
  * Fei  _et al._ [2021a] J. Fei, C.-N. Yeh, and E. Gull, Physical Review Letters 126, 056402 (2021a).
  * Jarrell and Gubernatis [1996] M. Jarrell and J. Gubernatis, Physics Reports 269, 133 (1996).
  * Fei  _et al._ [2021b] J. Fei, C.-N. Yeh, D. Zgid, and E. Gull, Phys. Rev. B 104, 165111 (2021b).
  * Bochner [1932] S. Bochner,  _Vorlesungen über Fouriersche Integrale_ (Leipzig : Akademische Verlagsgesellschaft, 1932).
  * Sasvári [1994] Z. Sasvári,  _Positive Definite and Definitizable Functions_ (Berlin: Akademie Verlag, 1994.).
  * Note [1] Note that in the standard terminology that positive definite functions include zero eigenvalues [33].
  * Balzer and Eckstein [2014] K. Balzer and M. Eckstein, Phys. Rev. B 89, 035148 (2014).
  * Gramsch and Potthoff [2015] C. Gramsch and M. Potthoff, Phys. Rev. B 92, 235135 (2015).
  * Higham [1988] N. J. Higham, Linear Algebra and its Applications 103, 103 (1988).
  * Eberle and Maciel [2003] M. G. Eberle and M. C. Maciel, Computational & Applied Mathematics 22, 1–18 (2003).
  * Bauschke and Borwein [1996] H. H. Bauschke and J. M. Borwein, SIAM Review 38, 367 (1996), https://doi.org/10.1137/S0036144593251710 .
  * Bauschke  _et al._ [2002] H. H. Bauschke, P. L. Combettes, and D. R. Luke, Journal of the Optical Society of America. A, Optics, image science, and vision 19 7, 1334 (2002).
  * Lewis  _et al._ [2009] A. S. Lewis, D. R. Luke, and J. Malick, Foundations of Computational Mathematics 9, 485 (2009).
  * Carathéodory [1907] C. Carathéodory, Mathematische Annalen 64, 95 (1907).
  * Krein [1940] M. Krein, in  _Comptes Rendus de l’Académie des Sciences de l’URSS_ , Vol. 46 (1940) pp. 17–22.
  * Sasvári [2006] Z. Sasvári, The extension problem for positive definite functions. a short historical survey, in  _Operator Theory and Indefinite Inner Product Spaces: Presented on the occasion of the retirement of Heinz Langer in the Colloquium on Operator Theory, Vienna, March 2004_, edited by M. Langer, A. Luger, and H. Woracek (Birkhäuser Basel, Basel, 2006) pp. 365–379.
  * Note [2] We thank Z. Sasvári for suggesting this approach.
  * Caratheodory and Fejér [1911] C. Caratheodory and L. Fejér, Rendiconti del Circolo Matematico di Palermo (1884-1940) 32, 218 (1911).
  * Schmidt [1986] R. Schmidt, IEEE Transactions on Antennas and Propagation 34, 276 (1986).
  * Paulsen and Raghupathi [2016] V. I. Paulsen and M. Raghupathi,  _An Introduction to the Theory of Reproducing Kernel Hilbert Spaces_, Cambridge Studies in Advanced Mathematics (Cambridge University Press, 2016).
  * Kökcü  _et al._ [2023] E. Kökcü, H. A. Labib, J. Freericks, and A. F. Kemper, arXiv preprint arXiv:2302.10219 (2023).
  * Cruz and Magano [2022] D. Cruz and D. Magano, arXiv preprint arXiv:2210.04919 (2022).
  * Steckmann  _et al._ [2023] T. Steckmann, T. Keen, E. Kökcü, A. F. Kemper, E. F. Dumitrescu, and Y. Wang, Physical Review Research 5, 023198 (2023).
  * Shen  _et al._ [2023] Y. Shen, D. Camps, S. Darbha, A. Szasz, K. Klymko, N. M. Tubman, R. Van Beeumen, _et al._ , arXiv preprint arXiv:2306.01858 (2023).
  * Uhrig  _et al._ [2019] G. S. Uhrig, M. H. Kalthoff, and J. K. Freericks, Physical Review Letters 122, 130604 (2019).
  * Cortes  _et al._ [2022] C. L. Cortes, A. E. DePrince III, and S. K. Gray, Physical Review A 106, 042409 (2022).
  * Klymko  _et al._ [2022] K. Klymko, C. Mejuto-Zaera, S. J. Cotton, F. Wudarski, M. Urbanek, D. Hait, M. Head-Gordon, K. B. Whaley, J. Moussa, N. Wiebe, _et al._ , PRX Quantum 3, 020323 (2022).
  * Ding and Lin [2023a] Z. Ding and L. Lin, PRX Quantum 4, 020331 (2023a).
  * Ding and Lin [2023b] Z. Ding and L. Lin, arXiv preprint arXiv:2303.05714 (2023b).
  * Schuckert  _et al._ [2022] A. Schuckert, A. Bohrdt, E. Crane, and M. Knap, arXiv preprint arXiv:2206.01756 (2022).
  * Shi and Geva [2003] Q. Shi and E. Geva, The Journal of Chemical Physics 119, 12063 (2003), https://pubs.aip.org/aip/jcp/article-pdf/119/23/12063/10853520/12063_1_online.pdf .
  * Cohen and Rabani [2011] G. Cohen and E. Rabani, Phys. Rev. B 84, 075150 (2011).
  * Cohen  _et al._ [2013] G. Cohen, E. Gull, D. R. Reichman, A. J. Millis, and E. Rabani, Phys. Rev. B 87, 195108 (2013).
  * Cerrillo and Cao [2014] J. Cerrillo and J. Cao, Phys. Rev. Lett. 112, 110401 (2014).
  * Pollock and Modi [2018] F. A. Pollock and K. Modi, Quantum 2, 76 (2018).
  * Pollock  _et al._ [2022] F. A. Pollock, E. Gull, K. Modi, and G. Cohen, SciPost Phys. 13, 027 (2022).
  * Cohen  _et al._ [2015] G. Cohen, E. Gull, D. R. Reichman, and A. J. Millis, Phys. Rev. Lett. 115, 266802 (2015).

##  Appendix A Application to quantum computer data

Fig. S1 shows the greater Green’s function obtained from IBM quantum computers
for the Su-Schrieffer-Heeger (SSH) model

$\displaystyle\mathcal{H}=-\sum_{\langle
i,j\rangle}\left[V_{nn}+(-1)^{i}\delta/2\right]c^{\dagger}_{i}c_{j}-\mu\sum_{i}c^{\dagger}_{i}c_{i}.$
(12)

This data was originally presented in Kökcü et al. [49]. The figure shows the
data for $\delta=0.4$, $\mu=-3$, and momenta
$k\in\\{0,1,2,3,4\\}\frac{2\pi}{8}$. Note that while the Fourier transform
power spectrum shows peaks near the expected values (indicated as vertical
dashed lines), there is a substantial amount of noise present.

Figure S1: Top panel: Raw data for an 8-site Green’s function
$\mathrm{Re}G_{k}(t)$ and corresponding power spectrum for the SSH model with
$\delta=0.4$ and $\mu=-3$ obtained from IBM’s quantum computer ibm_auckland
(originally from Ref. [49]). Bottom panel: Data after PSD denoising as
discussed below. Individual curves are $k$ points. The dashed lines indicate
the expected frequencies obtained by analytically solving Eq. 12.

The Gramian of the Green’s function, i.e. the Toeplitz matrix constructed from
the Green’s function, must follow the structure of Eq. 10 and be positive
semidefinite. In the presence of noise, the Gramian is a Toeplitz matrix but
is not positive definite. Defining a cost function

$\displaystyle\mathcal{C}=\sum_{i}\left[\lambda_{i}\left(\mathrm{sgn}(\lambda_{i})-1\right)\right]^{2}$
(13)

where $\lambda_{i}$ are the eigenvalues of $G$, we cycle through the
individual entries in $G$ and numerically optimize them with respect to this
cost function. In cases with very large noise, such as with data from quantum
computers, we find this procedure to be more reliable than the alternating
projection described in the main text.

We apply this approach to the data set from Ref. [49] with the worst quality,
$\delta=0.4$; the results are shown in Fig. S1. For most of the momenta, the
assertion that the Green’s function should have the structure of Eq. 10 yield
quite good results — the algorithm correctly picks out the correct peak(s) in
the Fourier transform. The notable exception is the $k=\pi$, which is because
the original data is particularly poor.

##  Appendix B Application to steady-state Monte Carlo data

Figure S2: De-noised Quantum Monte Carlo data. Main panel: inchworm [65]
quantum Monte Carlo of a steady-state transport problem (for data and
parameters see Fig. 1 of [28]), real part (orange) and imaginary part (blue).
Denoised data, real part (black) and imaginary part(red). Inset: Zoom into the
region near t=2 highlighting the effect of denoising.

In this section, we present de-noised data from a real-time simulation with
current state-of-the-art quantum Monte Carlo technology. We use the steady-
state formulation [28] of the Inchworm [65] Monte Carlo real-time propagation
scheme. The method overcomes the dynamical sign problem inherent to many real-
time Monte Carlo methods [17, 21] and is therefore able to simulate transport
data for long times.

The data presented has Monte Carlo noise whose size differs as a function of
time but is largest in the region near $t=2$. Visible are also data ‘glitches’
near $t=9.1$ which are likely due to ergodicity issues in the simulation and
which are clearly unphysical.

We employ the alternate projection described in the main text, while keeping
the norm fixed to the known analytical value of 1. It is evident that the
causal projection onto a positive definite function eliminates much of the
Monte Carlo noise and results in a real-time function that is much smoother.
Unphysical outliers in the data are eliminated. Obtaining data of the same
quality as the de-noised data would require a substantially larger investment
of computer time.

##  Appendix C The Hubbard Dimer – Details

For the data in Fig. 1, we consider a Hubbard dimer (a two-site Hubbard model)
with the following Hamiltonian:

$\displaystyle H$ $\displaystyle=H_{0}+H_{I},$ (14) $\displaystyle H_{0}$
$\displaystyle=-\epsilon\sum_{i=1,2}\left(n_{i\uparrow}+n_{i\downarrow}\right)-v\sum_{\sigma=\uparrow,\downarrow}c_{0\sigma}^{\dagger}c_{1\sigma}+c_{1\sigma}^{\dagger}c_{0\sigma},$
(15) $\displaystyle H_{I}$
$\displaystyle=U\sum_{i=1,2}\left(n_{i\uparrow}n_{i\downarrow}-\frac{1}{2}(n_{i\uparrow}+n_{i\downarrow})\right)$
(16)

The operators $c_{i\sigma}$ describe four fermion annihilation operators that
annihilate particles on site $i=1,2$ with spin $\sigma=\uparrow,\downarrow$.
Their adjoints $c_{i\sigma}^{d}agger$ create the corresponding particles on
site $i$ with spin $\sigma$. These operators correspond to matrices of size
$16\times 16$ in Fock space and satisfy the fermion commutation rules.

The parameter $U$ describes the on-site repulsion and is set to $U=5$.
$\epsilon$ is an on-site level energy term ($\epsilon=0$ corresponds to half
filling), and $v$ describes the hopping between sites $1$ and $2$. At
temperature $T=1/\beta=0.1,$ the choice of $\epsilon=2.3$ corresponds to a
density of $0.710556$ on all spin-orbitals.

The Green’s function discussed in Fig. 1 is defined as

$\displaystyle
G_{11\uparrow}(t)=\frac{1}{Z}\text{Tr}\left[e^{(-\beta+it)H}c_{1\uparrow}e^{-itH}c_{1\uparrow}^{\dagger}\right],$
(17)

and is a positive definite function that is symmetric in the real part and
anti-symmetric in the imaginary part, with $G_{11\uparrow}(t=0)=0.289444.$
(note that the system is degenerate in site and spin indices). The positive
definite function corresponds, up to a factor of $-i$, to the spin-up on-site
greater Green’s function of site $1$ of the system.

##  Appendix D Additional Mathematical notes

In an attempt to make the presentation in the main text accessible to most
physicists, we have eliminated some of the more formal mathematical aspects.
We would like to highlight the following additional connections:

###  D.1 Inner products induce positive definite functions

We reframe the time dependence of the operators as a map on $\mathcal{L}(V)$.
Note that the time argument is an element of the group of real numbers
$\mathbb{R}$ in the continuous case. When time is discretized (e.g. due to a
numerical sampling on an equidistant time grid), it is indexed by an element
the group of integers $\mathbb{Z}$; we will generically use $\mathcal{G}$ to
refer to the group. The time evolution of the operators in the Heisenberg
picture is a map $\phi_{t}$ that acts on elements of $\mathcal{L}(V)$. That
is, for an operator $A(t^{\prime})\in\mathcal{L}(V)$ and a time argument
$t\in\mathbb{R}$,

$\displaystyle\phi_{t}:\ $
$\displaystyle\mathcal{L}(V)\rightarrow\mathcal{L}(V)$ (18)
$\displaystyle\phi_{t}A(t^{\prime})=A(t+t^{\prime})$ (19)

It is straightforward to show that $\phi_{t}$ preserves that group structure
of the time argument: $\phi_{t}\phi_{t^{\prime}}=\phi_{t+t^{\prime}}$. Thus,
it is a representation of $\mathcal{G}$. Moreover, for the inner product in
Eq. 2, this map is unitary. This can be seen from the relations

$\displaystyle G_{AB}(t-t^{\prime})=\langle
A(t),B(t^{\prime})\rangle=\langle\phi_{t}A,\phi_{t^{\prime}}B\rangle$ (20)
$\displaystyle=\langle\phi_{t-t^{\prime}}A,B\rangle=\langle
A,\phi_{t^{\prime}-t}B\rangle$ (21)

which follow from the cyclicity of the trace. Combining these two points, and
because $\phi$ is a unitary representation of $\mathcal{G}$, we can show the
following theorem.

######  Theorem 1.

Consider a correlation function of two operators on Fock space $\langle
A(t_{i}),A(t_{j})\rangle$. This correlation function is a positive definite
function.

###### Proof.

For a positive definite function evaluated on an arbitrary set of time points
$t$, $f(t_{i},t_{j})$, the relation

$\displaystyle\sum_{i,j}f(t_{i},t_{j})\lambda_{i}^{*}\lambda_{j}\geq 0$ (22)

holds for any complex vector $\lambda$. Expanding out the definition of
$\langle A(t_{i}),A(t_{j})\rangle$, it is straightforward to show that

$\displaystyle\langle
A(t_{i}),A(t_{j})\rangle=\mathrm{Tr}\left[\rho\left(\sum_{i}\lambda_{i}\phi_{i}(A)\right)^{\dagger}\left(\sum_{j}\lambda_{j}\phi_{j}(A)\right)\right],$
(23)

which is the inner product of an element of $\mathcal{L}(V)$ with itself and
thus must be $\geq 0$. ∎

Thus, $G_{AA}(t-t^{\prime})$ is a positive definite function.

###  D.2 Connection to Reproducing Kernel Hilbert Spaces

Inner products induce reproducing kernel Hilbert spaces. Given a Hilbert space
$\mathcal{L}$ with inner product $\langle\cdot,\cdot\rangle$, the function
$K(x,y)=\langle x,y\rangle$ is a kernel function and determines a
corresponding reproducing kernel Hilbert space [48] on the space of functions
$\mathbb{R}\rightarrow\mathbb{C}$. Reproducing kernel Hilbert spaces form an
important tool in statistics and machine learning, with connections to
interpolation and approximation theory. This connection will be further
examined in a follow-up paper.