# Hypothesis testing of scientific Monte Carlo calculations

Markus Wallerberger, Emanuel Gull  Department of Physics, University of
Michigan, Ann Arbor, MI 48109

###### Abstract

The steadily increasing size of scientific Monte Carlo simulations and the
desire for robust, correct, and reproducible results necessitates rigorous
testing procedures for scientific simulations in order to detect numerical
problems and programming bugs. However, the testing paradigms developed for
deterministic algorithms have proven to be ill suited for stochastic
algorithms. In this paper we demonstrate explicitly how the technique of
statistical hypothesis testing, which is in wide use in other fields of
science, can be used to devise automatic and reliable tests for Monte Carlo
methods, and we show that these tests are able to detect some of the common
problems encountered in stochastic scientific simulations. We argue that
hypothesis testing should become part of the standard testing toolkit for
scientific simulations.

##  I Introduction

Scientific computing, i.e., the process of obtaining numerical results from
scientific theories using algorithms, relies on correct and reproducible
implementations of computer programs. In condensed matter and statistical
physics, these computer programs were traditionally small, often implemented
by a single researcher, and tested and debugged by hand until no more problems
could be found.

Over time, the size and complexity of programs in this field has grown
rapidly. For example, computer programs for complex many-body problems, such
as finding the ground state energy of an interacting solid Kim  _et al._
(2012) or evaluating response functions of correlated quantum impurity models
Gaenko  _et al._ (2017); Parcollet  _et al._ (2015), now span hundreds of
thousands of lines that are developed and maintained by large and constantly
changing teams. For such programs, manual testing becomes inefficient and
expensive.

This challenge is not unique to scientific computing, and software engineering
has responded by establishing automated testing practices. The corresponding
arsenal of methods includes, in order of increasing granularity: contract
programming, where invariants in the program state are verified continuously
during execution Meyer (1992); unit tests, which ensure the correctness of
small sections of the code Binder (2000); as well as integration and system
tests, which check that implementations yield correct non-trivial results for
predefined benchmark problems British Computer Society Working Group on
Testing (1986).

These techniques have permeated scientific software engineering Dubois (2005),
and they are by now standard in many computational science packages. Combined
with continuous testing, i.e., the automatic execution of tests after a change
to the code base, they have led to a massive improvement of the quality and
resilience of scientific software Sanders and Kelly (2008).

Nevertheless, there is a large part of computational and statistical physics
where such tests were so far not practical, namely the field of stochastic
Monte Carlo simulations. In this domain, results make use of random or pseudo-
random number generators, and are therefore intrinsically stochastic in
nature. Agreement with a reference result has to be “within error bars” only.

As far as we are aware, most practitioners of these techniques therefore
either enforce a deterministic procedure (e.g., a simulation with a fixed seed
of the pseudo-random number generator or an otherwise fixed sequence of
updates on a given configuration) or resort to “visual inspection” of the
results to determine agreement between simulation and reference, neither of
which is optimal. The former breaks whenever the sequence or ratio of updates
are changed, and therefore is prone to false negatives, i.e., failed tests
even though the results are correct. The latter relies on human intervention
and is therefore neither reliable nor automatable.

In this paper, we show how tools of statistics Stuart  _et al._ (2010), known
for more than a century and in wide use in many fields, should be used to
construct automated tests for physics simulations. Our formulations are
general and applicable to any stochastic simulation. While we are not aware of
applications to physics so far, we emphasize that similar applications have
been pioneered both in the field of image synthesis Ševčíková  _et al._ (2006)
and urban simulations Subr and Arvo (2007).

In the remainder of this paper we will introduce the concept of statistical
testing or “hypothesis testing” in Sec. II and III, with applications to the
two-dimensional Ising model. Sec. IV shows an application to the Anderson
impurity model, and Sec. V summarizes our conclusions.

##  II Scalar Tests

###  II.1 One-sample test for the mean

The basic idea of statistical hypothesis testing in the context of Monte Carlo
is straight-forward: one first chooses a model for which an exact benchmark
result $y$ exists. The null hypothesis, $H_{0}$, is that there is no
significant difference between this reference result and the expectation value
$\mathrm{E}[\hat{X}]$ of a simulation with the estimator $\hat{X}$ Subr and
Arvo (2007). The alternate hypothesis, $H_{1}$, is that this is not the case:

$\displaystyle H_{0}$ $\displaystyle:\mathrm{E}[\hat{X}]=y$ (1a)
$\displaystyle H_{1}$ $\displaystyle:\mathrm{E}[\hat{X}]\neq y.$ (1b)

We first discuss the scalar case. Let $\hat{X}$ be a “simple” Monte Carlo
estimator, i.e., an average $\langle X\rangle$ over $N$ independent random
variables identically distributed according to $X$. (In the case of sampling
on a Markov chain, one has to correct the number $N^{\prime}$ of Monte Carlo
samples by the integrated autocorrelation time:
$N=N^{\prime}/\tau_{\mathrm{int},X}$.) We then find:

$\frac{\langle X\rangle-y}{\sigma_{X}/\sqrt{N}}\sim t_{N-1},$ (2)

where $\sim$ is shorthand for “is distributed according to”, $t_{\nu}$ is
Student’s $t$ distribution for $\nu$ degrees of freedom and $\sigma_{X}^{2}$
is the variance of $X$.

Following standard practice Neyman and Pearson (1933), we turn Eq. (2) into a
likelihood estimate for $H_{0}$, known as Student’s $t$ test. We compute the
two-sided $p$-value as $p=2P^{-1}(-|z|)$, where $P^{-1}$ is the inverse of the
cumulative distribution function of the right-hand side and $z$ is the
observed left-hand side in Eq. (2). Finally, we compare $p$ with a
significance level $\alpha\in(0,1)$ and reject the null hypothesis (1a) if
$p<\alpha$. In other words, the $p$ value is the probability of observing $z$
or a “more unlikely” event given $H_{0}$, and we reject the $H_{0}$ if that
probability becomes smaller than $\alpha$.

Figure 1: Scalar one-sample test against $\langle m\rangle=0$ for single spin-
flip updates (red curves) and Wolff cluster updates (blue curves) in a
classical two-dimensional Ising model with length $L=16$: (a) result for
$\langle m\rangle$ from $N^{\prime}=10^{6}$ Monte Carlo sweeps (a sweep is
either a set of $L^{2}$ single spin flips or a single cluster update); (b)
$|t|$ score as the left-hand side of Eq. 2; and (c) $p$ values from a two-
tailed test of the $t$ score against the Student $t_{N-1}$ distribution (the
shaded area indicates $p<\alpha=0.01$ and thus a failed test).

Let us illustrate the procedure with a simple example, the ferromagnetic Ising
model Landau and Binder (2005)

$\mathcal{H}=-\sum_{\langle ij\rangle}\sigma_{i}\sigma_{j},$ (3)

where $\langle ij\rangle$ runs over all pairs of directly neighboring Ising
spins $\sigma_{i}\in\\{1,-1\\}$ on a $L\times L$ square lattice with periodic
boundary conditions and $L=16$. Since the system is finite and there is no
external magnetic field, $\langle m\rangle=0$. We perform a Markov chain Monte
Carlo simulation Metropolis  _et al._ (1953) for Eq. (3) for two different
types of updates: (a) a set of single spin flips $\sigma_{i}\to-\sigma_{i}$,
and (b) Wolff cluster updates Wolff (1989). In both cases, the magnetization
estimator is constructed as $\hat{m}=\langle\sum_{i}\sigma_{i}\rangle/L^{2}$.

Fig. 1 shows the temperature-dependent magnetization curves obtained by the
simulation. From Fig. 1(a), we immediately see that the single spin flip
updates (red curve) produce a spurious spin polarization at low temperature
for the parameters chosen. This is to be expected, since in order to restore
$\langle m\rangle=0$, all spins must be flipped, which due to the
exponentially divergent autocorrelation time $\tau\propto\exp(L)$ requires far
more updates than performed in our test. Figure 1(b) shows the $|t|$ score or
deviation in units of the standard error computed from Eq. 2. Figure 1(c)
shows the $p$ value as result of a two-tailed test with the Student
distribution, which amounts to $p=2P^{-1}(-|t|)$. If we choose a significance
level of $\alpha=0.01$, we see that the test fails for all temperatures below
the critical temperature, $T<2.2$. In contrast, the Wolff updates, which
circumvent the problem of divergent autocorrelation times by updating clusters
of spins, pass the test for all temperatures.

The spurious spin polarization is already obvious from a fleeting inspection
of Fig. 1(a), and a formal verification of Eq. (1b) may seem superfluous.
However, we emphasize that the formal procedure can easily be turned into an
automated test and run as part of an automated test suite. This extends the
test coverage from the deterministic parts of the algorithm to the stochastic
updates and the magnetization estimator and its autocorrelation effects.

The choice of significance level $\alpha$ is a trade-off between the
probability of two kinds of errors:

$\displaystyle\alpha$
$\displaystyle=P(H_{0}\mathrm{\,rejected}\,|\,H_{0}\mathrm{\,is\,true})$ (4a)
$\displaystyle\beta$
$\displaystyle=P(H_{0}\mathrm{\,accepted}\,|\,H_{0}\mathrm{\,is\,false}),$
(4b)

known as type-I and type-II errors, or false positives and false negatives,
respectively. We empirically find that the rather conservative $\alpha\approx
0.01$ provides such a good trade-off for a single test. In the case of a test
suite of $K$ tests, one can either substitute
$\alpha\to\alpha^{\prime}\approx\alpha/K$ to keep the probability of a type-I
error constant or keep the threshold as-is to keep the probability of a type-
II error constant. The former scheme is suited for automatized stochastic unit
tests, the latter strategy is advantageous when combined with test refinement.
In such a scheme, we choose a window $p\in[\alpha,\beta)$ corresponding to
ambiguous test results and re-run these cases with double the number of
samples until they are either accepted or rejected.

###  II.2 Two-sample test; biased estimator

In many cases, exact benchmark results may not be available or cumbersome to
obtain. In these cases, we can also compare two stochastic results: the
estimator $\hat{X}$ to be tested and a trusted estimator $\hat{Y}$. This
corresponds to replacing $y$ with $\mathrm{E}[\hat{Y}]$ in Eqs. 1:

$\displaystyle H_{0}$ $\displaystyle:\mathrm{E}[\hat{X}]=\mathrm{E}[\hat{Y}]$
(5a) $\displaystyle H_{1}$
$\displaystyle:\mathrm{E}[\hat{X}]\neq\mathrm{E}[\hat{Y}]$ (5b)

In the scalar case, with $\hat{X}$ and $\hat{Y}$ averages over $N_{X}$ and
$N_{Y}$ independent random variables distributed according to $X$ and $Y$, we
find the analogue of Eq. (2),

$\frac{\langle X\rangle-\langle Y\rangle}{\sigma/N_{\mu}}\sim
t_{N_{X}+N_{Y}-2},$ (6)

with $N_{\mu}^{-1}=N_{X}^{-1}+N_{Y}^{-1}$ and the pooled variance

$\sigma^{2}=\frac{(N_{X}-1)\sigma_{X}^{2}+(N_{Y}-1)\sigma_{Y}^{2}}{N_{X}+N_{Y}-2}.$
(7)

The rest of the test proceeds exactly the same as for the one-sample test
(Sec. II.1).

Figure 2: Scalar two-sample test for the Binder cumulant $U_{4}$ in the Ising
model obtained from Markov chain Monte Carlo with Wolff cluster updates at
temperature $T=2.3$ for system length $L=32$. (a) Sample mean of $U_{4}$ of an
erroneous implementation for different simulation times $N$, and the result of
a correct reference simulation; (b) corresponding $p$ values computed using
Eq. (6); (c) illustration of the $p$ value for $N^{\prime}=10^{3}$ ($N=123$
after the removal of autocorrelation and variance pooling) as the shaded area
under the probability density function (p.d.f.) of the corresponding $t$
distribution.

As an example, we reexamine data from the Ising model, this time on a
$32\times 32$ square lattice. We verify the estimator for the Binder cumulant
Binder (1981)

$\hat{U}_{4}=\frac{\langle m^{4}\rangle}{1-3\langle m^{2}\rangle^{2}}.$ (8)

The Student $t^{2}$ test is sensitive to non-Gaussian distributed errors,
which occur in the computation of Eq. (8) due to non-linear error propagation.
To remedy this, we use the jackknife resampling procedure, which replaces Eq.
(8) with a simple average $\langle U_{4}^{\prime}\rangle$ over pseudovalues
$U_{4}^{\prime}$, removing the linear order of the bias and restoring the
validity of Eq. (6) Efron (1982). Alternatively, one could abandon the the
Student test altogether in favor of the parametric bootstrap method Efron
(1982). However, we will see that the jackknife method suffices in our case.

To simulate a common programming error, we have artificially broken periodic
boundary conditions on the corners of the lattice (they are reduced to having
two neighbors each). Fig. 2(a) compares this erroneous implementation (green
curve) with a simulation result where the error is not present. As evident
from Fig. 2(b), as we increase the number of Monte Carlo sweeps $N^{\prime}$,
the error bars shrink and the null hypothesis (5a) is rejected more and more
strongly.

##  III Data series tests

###  III.1 Tests for the mean

While tests for scalar quantities (Sec. II) are useful, we empirically find
that it is often easier to identify problems when comparing functions and data
series. In the case of a one-sample test, this corresponds to the benchmark
result $y$ being a vector of $n$ elements rather than a scalar. Consequently,
the Monte Carlo estimator $\hat{X}$ is vector-valued. Again assuming
independent and identically distributed results, Eq. (2) is replaced by
Hotelling (1931)

$N(\langle X\rangle-x_{0})^{\mathrm{T}}\Sigma_{X}^{-1}(\langle X\rangle-
x_{0})\sim\frac{n(N-1)}{N-n}F_{n,N-n}\;,$ (9)

where $\langle X\rangle$ is the sample mean, $\Sigma_{X}$ is the sample
covariance matrix, and $F_{a,b}$ is the Fisher–Snedecor distribution with
parameters $a$, $b$. One proceeds in a similar way to the Student’s $t$-test.
The observed left-hand side of Eq. (9) is again used as the test statistic and
checked against the right-hand side distribution. However, since the $F$
distribution is not symmetric for low $n$ (cf. Fig. 3(c)), one uses two one-
sided tests instead of a two-sided test and subsequently obtains two
$p$-values, which we will call $p_{<}$ and $p_{>}$. This is known as
Hotelling’s $T^{2}$ test.

In the case where we compare the estimator to a trusted result $\langle
Y\rangle$, we proceed similar as in Sec. II.2 and replace Eq. (9) with:

$\begin{split}&N_{\mu}(\langle X\rangle-\langle
Y\rangle)^{\mathrm{T}}\Sigma^{-1}(\langle X\rangle-\langle Y\rangle)\\\
&\qquad\sim\frac{n(N_{X}+N_{Y}-2)}{N_{X}+N_{Y}-n-1}F_{n,N_{X}+N_{Y}-n-1}\;,\end{split}$
(10)

where $\Sigma$ is the pooled covariance obtained by replacing all variances
$\sigma_{a}^{2}$ with covariance matrices $\Sigma_{a}$ in Eq. (7).

Figure 3: Vector-valued two-sample test on the spin correlation function
$\chi_{x,0}$ (cf. Eq. (11)). (a) Simulation result for Wolff updates (black
curve) and single spin-flip updates (red curve); (b) deviation of spin-flip
from Wolff update (the shaded region are the Wolff result error bars); (c)
Projected deviations and errors (numerator and denominator in the l.h.s. of
Eq. (14)); (d) correlation matrix (12); (e) p.d.f. of the corresponding $F$
distribution in Eq. (14) with the mean as vertical blue line and $p_{<}$-value
as blue shaded area to the left of the observed $T^{2}$ score (blue dot) as
well as test failure threshold as red-shaded area.

In order to illustrate the procedure, we revisit our Ising model example for
$L=32$ and $T=2.3$ (close to the critical temperature) and examine the spin
correlation function

$\begin{split}\chi_{x,y}&=\langle\sigma_{0,0}\sigma_{x,y}\rangle\\\
&=\frac{1}{L^{2}}\langle\sum_{k,q}\sum_{x^{\prime},y^{\prime}}\mathcal{F}_{x,y;k,q}^{-1}|\mathcal{F}_{k,q;x^{\prime},y^{\prime}}\sigma_{x^{\prime},y^{\prime}}|^{2}\rangle,\end{split}$
(11)

where $(x,y)$ denote row and column of the lattice site, and $\mathcal{F}$
denotes the discrete Fourier transform used in the actual estimator. Fig. 3(a)
shows $\chi_{x0}$ for the Wolff cluster update (black curve), which we take as
the trusted result, and for a set of spin-flip updates (red curve). The inset
Fig. 3(b) shows the deviation of the red curve from the black one, where the
shaded region marks the error bars of the Wolff update result. We see
significant correlation of the error bars, which underscores the importance of
a proper treatment of the covariance matrix (cf. Fig. 3(d)).

###  III.2 Cross-correlated data

A common complication with the $T^{2}$ test are perfect correlation or anti-
correlation within the dataset (duplicates), which implies a singular
covariance matrix in Eq. (9). In our example, the symmetry of the system
implies $\chi_{x,y}=x_{L-x,y}$, thus half of the points yielded by the
estimator (11) are just copies of the other half. We can confirm this by
examining the correlation matrix:

$\mathrm{Corr}[\chi_{x,0},\chi_{x^{\prime},0}]=\frac{\mathrm{Cov}[\chi_{x,0},\chi_{x^{\prime},0}]}{\sqrt{\mathrm{Var}[\chi_{x,0}]\mathrm{Var}[\chi_{x^{\prime},0}]}},$
(12)

plotted in Fig. 3(d), which is one on the anti-diagonal.

This can be solved by first diagonalizing $\Sigma$ and retaining only the non-
zero eigenvalues (a relative threshold of $10^{-14}$ seems to be practical for
most cases we studied):

$\Sigma=\mathcal{P}\;\mathrm{diag}(s_{1}^{2},\ldots,s_{m}^{2})\;\mathcal{P}^{\mathrm{T}},$
(13)

where $\mathcal{P}$ is the $n\times m$ projection to the non-zero eigenvalues
$s_{1}^{2},\ldots,s_{m}^{2}$. This is shown in Fig. 3(c), where there is sharp
drop of $s_{n}$ (red curve) in magnitude after $m=15$. Eq. (9) is then amended
to:

$\sum_{i=1}^{m}\frac{|\sum_{k=1}^{n}\mathcal{P}_{ki}(\langle X_{k}\rangle-
y_{k})|^{2}}{s_{i}^{2}/N}\sim\frac{m(N-1)}{N-m}F_{m,N-m}.$ (14)

Note the reduction in the degrees of freedom from $n$ to $m$, which
corresponds to discarding the $n-m$ correlated data points. Note also that for
$n=m$ Eq. (9) and Eq. (14) are equivalent, such that in practical
calculations, one can always use Eq. (14). Finally, we perform a $T^{2}$ test
against the appropriate $F$ distribution and find that the null hypothesis is
accepted with $p\approx 0.2$ (Fig. 3(e)).

###  III.3 Error bars

By using the sample mean and covariance as input rather than the individual
samples, one can interpret Hotelling’s $t^{2}$ test as statistical test on the
error bars $\sigma_{0}$ :

$\displaystyle H_{0}$ $\displaystyle:\sigma=\sigma_{0}$ (15a) $\displaystyle
H_{1}^{-}$ $\displaystyle:\sigma<\sigma_{0}$ (15b) $\displaystyle H_{1}^{+}$
$\displaystyle:\sigma>\sigma_{0}$ (15c)

For a scalar estimator (Sec. II), we can distinguish $H_{0}$ from $H_{1}^{-}$:
error bars being “too small” (15b) is equivalent to the result being
inconsistent with the benchmark (Eq. (1b)). However, we cannot test against
$H_{1}^{+}$, since we may have accidentally hit the benchmark accurately.
Using a data series, we can also distinguish it from $H_{1}^{+}$, formalizing
the rule that “roughly two-thirds of the data should fall within one-sigma
error-bars”. This is reflected in the fact that for $n>1$, the $F$
distribution turns from a one-tailed to a two-tailed distribution, and becomes
more symmetric around $1$ as $n$ gets larger. We can make use of this by
testing the lower tail as score for $H_{1}^{+}$ and the upper tail as score
for $H_{1}^{-}$.

This procedure is illustrated at the example of the estimator for $\chi_{i,0}$
(Sec. III). If we ignore the cross-correlation (Fig. 3d) and interpret the
error bars in Fig. 3b as uncorrelated errors, it is evident from visual
inspection that they are too large. We can confirm this numerically by
(erroneously) plugging the diagonal elements $\Sigma_{ii}$ of the covariance
matrix instead of its eigenvalues $s_{i}$ into Eq. (14). We then find a
$T^{2}$ score of $0.03$ and an acceptance of the lower alternate hypothesis
$H_{1}^{+}$ (15b) with $p=1-10^{-16}$.

##  IV Example: Anderson impurity model

Figure 4: (a) Green’s function $G(\tau)$ for AIM parameters as in main text:
Monte Carlo result (red) and exact result with an artificially introduced
shift of the half bin size, i.e., $G(\tau-0.02)$, (black) modeling a binning
error; (b) deviation from the exact result in multiples of the standard error
as well as a linear regression (black dashed line) $0.539-0.098\tau$ with a
goodness of fit of $R^{2}\approx 0.065$.

To illustrate our method on a research example, we examine the single-orbital
Anderson impurity model Anderson (1961) (AIM) which characterizes a few
discrete and potentially correlated impurity states coupled to a non-
interacting bath. The model is in wide use in nano- and transport science
Hanson  _et al._ (2007); Brako and Newns (1981) and as an auxiliary model in
the dynamical mean field theory Georges  _et al._ (1996), and in many
parameter regimes quantum Monte Carlo methods are the standard tools for
obtaining its properties Gull  _et al._ (2011). Its Hamiltonian is

$\begin{split}\mathcal{H}=&\;Uc_{\uparrow}^{\dagger}c_{\downarrow}^{\dagger}c_{\downarrow}c_{\uparrow}-\mu\sum_{\sigma}c_{\sigma}^{\dagger}c_{\sigma}\\\
&+\sum_{p\sigma}(V_{p\sigma}f_{p\sigma}^{\dagger}c_{\sigma}+\mathrm{h.c.})+\sum_{p\sigma}\epsilon_{p}f_{p\sigma}^{\dagger}f_{p\sigma}.\end{split}$
(16)

Here, $c_{\sigma}$ annihilates a fermion of spin-$\sigma$ on the impurity and
$f_{p\sigma}$ annihilates a bath fermion of momentum $p$ and spin $\sigma$.
Impurity interactions are characterized by $U,$ $\mu$ denotes a chemical
potential, $V$ a spin- and momentum dependent hybridization term, and
$\epsilon_{p}$ a momentum-dependent bath dispersion. In the context of the
AIM, a truncation of the bath to a few states and subsequent exact
diagonalization of the finite system is particularly suitable for testing.
While the complexity of solving the model with Monte Carlo methods is the same
as for a model without bath truncation, one empirically finds that the
truncated model shares much of the physics of the AIM __ and can thus be used
to generate non-trivial, analytically accessible test cases for Monte Carlo
simulations.

Our example consists of two momenta and correspondingly two bath sites with
energies of $\epsilon_{p}=\pm 0.5$ and a hybridization strength $V=1$, as well
as $U=5$, $\mu=U/2$, and temperature $T=1/10$. Stochastic results were
obtained using continuous-time quantum Monte Carlo in the hybridization
expansion Werner and Millis (2006); Gull  _et al._ (2011).

The imaginary time Green’s function $G(\tau)=-\langle
Tc(\tau)c^{\dagger}(0)\rangle,$ which is the fundamental quantity of interest
in this model and which is directly related to the interacting spectral
function, is shown in Fig. 4. In order to mimic the effect of a typical
binning programming error, we have shifted the exact result (black) by half a
bin $G(\tau)\to G(\tau-0.02)$. The top panel shows that the Monte Carlo result
(red) is still consistent with the exact result in this case when gauged by
visual inspection. This is reinforced by the bottom panel, where the deviation
from the exact result in multiples of the standard error is plotted (red).
Overall we find the expected result, even though a linear fit of the data
(shown as black dashed line) shows a slight downward slope indicative of a
problem.

However, Hotelling’s $T^{2}$ test finds a test statistic of $T^{2}\approx
1.28$ and therefore a rejection of the null hypothesis in favor of $H_{1}^{+}$
with $p\approx 0.0026$ (about three sigma). This is because by using all
$n=250$ data points, the test becomes sensitive to a small increase of the
values outside of error bars. Systematically increasing the statistics would
eventually expose the deviation to visual inspection.

##  V Conclusions

In this paper, we have shown how hypothesis testing can be used to develop
tests for code correctness of Monte Carlo codes in statistical and condensed
matter physics. We also have shown how these tests are sensitive to different
types of simulation problems, and how they can therefore be used as diagnostic
tools to ensure the correctness of simulations.

The mathematical framework for hypothesis testing has been known for over 100
years and statistical tests are in wide use across many scientific fields.
Despite this, the technique is not used on a routine basis for testing
scientific simulation results. With the advent of automatic testing and unit
test frameworks, which have permeated most of software engineering and to some
extent also scientific computing, our techniques add to the testing toolkits
that can be used to systematically ensure correctness and reproducibility of
stochastic physics simulations. These tests integrate well into existing
testing frameworks and can validate parts of the programs that are otherwise
difficult to test.

Hypothesis testing allows to gain and keep trust in complex codes as they
undergo modifications, and to uncover problems that are difficult to uncover
by other means, e.g. manual visual inspection. This both increases the speed
of scientific software development and the trust in results produced by
complex computer programs.

In our opinion hypothesis testing should be widely adopted in statistical
simulation codes and should become a standard tool in scientific software
development. While implementing these tests carries a small overhead, we argue
that rigorous, frequent, and automatic testing is necessary for today’s codes,
especially in light of the replication crisis Baker (2016) observed in other
fields of science.

The code for the stochastic solvers and the hypothesis testing post-processing
scripts are available from the authors upon request. An open-source software
implementation of hypothesis testing is scheduled for inclusion in the
upcoming version of the ALPS core libraries.Gaenko  _et al._ (2017)

###### Acknowledgements.

The authors would like to thank Alexander Gaenko for fruitful discussions. MW
was supported by the Simons Foundation via the Simons Collaboration on the
Many-Electron Problem. EG was supported by DOE grant no. ER46932. This
research 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.

## References

  * Kim  _et al._ (2012) J. Kim, K. P. Esler, J. McMinis, M. A. Morales, B. K. Clark, L. Shulenburger, and D. M. Ceperley, Journal of Physics: Conference Series 402, 012008 (2012).
  * Gaenko  _et al._ (2017) A. Gaenko, A. Antipov, G. Carcassi, T. Chen, X. Chen, Q. Dong, L. Gamper, J. Gukelberger, R. Igarashi, S. Iskakov, M. Konz, J. LeBlanc, R. Levy, P. Ma, J. Paki, H. Shinaoka, S. Todo, M. Troyer, and E. Gull, Computer Physics Communications 213, 235 (2017).
  * Parcollet  _et al._ (2015) O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. Krivenko, L. Messio, and P. Seth, Computer Physics Communications 196, 398 (2015).
  * Meyer (1992) B. Meyer, Computer 25, 40 (1992).
  * Binder (2000) R. V. Binder,  _Testing Object-oriented Systems: Models, Patterns, and Tools_ (Addison-Wesley, 2000).
  * British Computer Society Working Group on Testing (1986) British Computer Society Working Group on Testing,  _Testing in Software Development_ , edited by M. Ould and C. Unwin (Cambridge University Press, 1986).
  * Dubois (2005) P. F. Dubois, Comput. Sci. Eng. 7, 80 (2005).
  * Sanders and Kelly (2008) R. Sanders and D. Kelly, IEEE Softw. 25, 21 (2008).
  * Stuart  _et al._ (2010) A. Stuart, K. Ord, and S. Arnold,  _Kendall’s Advanced Theory of Statistics – Classical Inference and the Linear Model_ , Kendall’s Advanced Theory of Statistics (Wiley, 2010).
  * Ševčíková  _et al._ (2006) H. Ševčíková, A. Borning, D. Socha, and W.-G. Bleek, in  _Proceedings of the 2006 International Symposium on Software Testing and Analysis_ , edited by L. Pollock and M. Pezzè (Association for Computing Machinery, New York, NY, 2006) p. 106.
  * Subr and Arvo (2007) K. Subr and J. Arvo, in  _Proceedings of the 15 th Pacific Conference on Computer Graphics and Applications_, edited by M. Alexa, S. Gortler, and Tao Ju (IEEE Computer Society, Los Alamitos, CA, 2007) p. 215.
  * Neyman and Pearson (1933) J. Neyman and E. S. Pearson, Phil. Trans. R. Soc. A 231, 289 (1933), http://rsta.royalsocietypublishing.org/content/231/694-706/289.full.pdf .
  * Landau and Binder (2005) D. Landau and K. Binder,  _A Guide to Monte Carlo Simulations in Statistical Physics_ (Cambridge University Press, New York, NY, USA, 2005).
  * Metropolis  _et al._ (1953) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, The Journal of Chemical Physics 21, 1087 (1953), http://dx.doi.org/10.1063/1.1699114 .
  * Wolff (1989) U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
  * Binder (1981) K. Binder, Zeitschrift für Physik B Condensed Matter 43, 119 (1981).
  * Efron (1982) B. Efron,  _The Jackknife, the Bootstrap and Other Resampling Plans_ (Society for Industrial and Applied Mathematics, 1982).
  * Hotelling (1931) H. Hotelling, Ann. Math. Statist. 2, 360 (1931).
  * Anderson (1961) P. W. Anderson, Phys. Rev. 124, 41 (1961).
  * Hanson  _et al._ (2007) R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, and L. M. K. Vandersypen, Rev. Mod. Phys. 79, 1217 (2007).
  * Brako and Newns (1981) R. Brako and D. M. Newns, Journal of Physics C: Solid State Physics 14, 3065 (1981).
  * Georges  _et al._ (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  * Gull  _et al._ (2011) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
  * Werner and Millis (2006) P. Werner and A. J. Millis, Phys. Rev. B 74, 155107 (2006).
  * Baker (2016) M. Baker, Nature 533, 452 (2016).


# Non-thermal superconductivity in photo-doped multi-orbital Hubbard systems

Sujay Ray  Department of Physics, University of Fribourg, Fribourg-1700,
Switzerland  Yuta Murakami  Center for Emergent Matter Science, RIKEN, Wako,
Saitama 351-0198, Japan  Philipp Werner  Department of Physics, University of
Fribourg, Fribourg-1700, Switzerland

###### Abstract

Superconductivity in laser-excited correlated electron systems has attracted
considerable interest due to reports of light-induced superconducting-like
states. Here we explore the possibility of non-thermal superconducting order
in strongly interacting multi-orbital Hubbard systems, using non-equilibrium
dynamical mean field theory. We find that a staggered $\eta$-type
superconducting phase can be realized on a bipartite lattice in the high
photo-doping regime, if the effective temperature of the photo-carriers is
sufficiently low. The $\eta$ superconducting state is stabilized by Hund
coupling – a positive Hund coupling favors orbital-singlet spin-triplet $\eta$
pairing, whereas a negative Hund coupling stabilizes spin-singlet orbital-
triplet $\eta$ pairing.

Introduction: The study of non-equilibrium quantum systems has become a
frontier topic of research in condensed matter physics and quantum field
theory int_qm1 ; int_qm2 ; int_qm3 ; int_qm4 . A particularly interesting
topic is the exploration of hidden phases of matter lat_struc ; mang_hid ,
that is, long-lived metastable states which can only be accessed via a
nonthermal pathway. In strongly correlated electron systems, which are often
characterized by competing low-energy states Dagotto2005 , non-equilibrium
protocols offer the intriguing opportunity to manipulate material properties
lat_struc ; mang_hid ; exciton_insulator ; mag_switch ; opt_exint ; opt_magord
; ind_sc_cuprate ; phind_sc  and to create novel states with non-thermal
electronic orders nat_rev ; AFM_mott ; Ono2017 ; KK_jiajun ; yuta_exciton ;
Kaneko2019 ; eta_philipp ; Tindall2019 ; oneband_jiajun ; Werner2021 ;
Murakami2023 ; Takashi2002 ; Rosch2008 ; Murakami2022 . A promising route for
realizing such non-equilibrium states is laser driving. This allows to create
Floquet states Oka2009 ; Wang2013  or photo-doped states Iwai2003 , which
under suitable circumstances can favor electron pairing and superconductivity.
In photo-doped large-gap Mott insulators, the charge excitations survive for a
long time Sensarma2010 ; Eckstein2011 ; Strohmaier2010 ; Lenarcic2014 . Rather
than relaxing back to a (heated) equilibrium state, the charge excitations
thermalize within the Hubbard bands, resulting in the formation of a
metastable non-thermal state which can be characterized by effective
temperatures of the charge carriers.

A broad range of theoretical tools has been employed to simulate photo-induced
non-equilibrium states Maeshima2005 ; Golez2012 ; Shinjo2017 ; Sandri2013 ;
Eckstein2013 ; Tacogne2018 ; Rincon2018 ; Zawadzki2019 ; Schueler2020 ;
Bittner2021 ; Perfetto2022 . One of the challenges in laser-driven systems is
the heating effect, which tends to suppress the emergence of electronic
orders. Sophisticated protocols have been developed to realize effectively
cold metastable systems. Entropy cooling techniques entropy_cooling ;
eta_philipp  have proven useful in overcoming relaxation bottlenecks and
enabled the preparation of effectively cold non-equilibrium states of photo-
doped Mott systems. Alternatively, such states can be realized by the coupling
to suitable fermionic NESS_martin  or bosonic baths Eckstein2013 . Both
techniques have been successfully used to simulate nonthermal superconducting
(SC) states in the single band repulsive Hubbard model eta_philipp ;
oneband_jiajun ; tri_chiral . SC states with staggered order parameter can be
stabilized over a wide photo-doping region on bipartite lattices
oneband_jiajun , while chiral SC order with a $120^{\circ}$ phase twist has
been induced on a triangular lattice tri_chiral . Here, we extend these
studies to multi-orbital Hubbard systems, more specifically the two-orbital
model with strong repulsive interactions.

Nonthermal orders in the two-orbital repulsive Hubbard model have been
previously studied. These investigations include nonthermal spin and orbital
orders in the three-quarter filled (Kugel-Khomskii) model KK_jiajun ;
Kugel_Khomskii ; Kugel_Khomskii_1 , and excitonic orders in the model with
crystal field splitting yuta_exciton . A photo-induced enhancement of a spin-
triplet pairing susceptibility has also been reported enhanced_pairing , but
due to computational limitations, this study could not reach a cold enough
state for symmetry breaking, and it did not consider the staggered order that
we will identify below as the dominant nonthermal phase.

In this letter, we employ the recently developed non-equilibrium steady state
(NESS) technique NESS_martin  to explore the photo-doped two-orbital repulsive
Hubbard model. With this method we are able to achieve low enough effective
temperatures to demonstrate the emergence of nonthermal SC states which are
specific to multi-orbital systems.

Model and Method: We consider the two-orbital Hubbard model with Hamiltonian

$\displaystyle
H=-t_{\text{hop}}\sum_{\left<ij\right>,\sigma}\sum_{\alpha=1,2}c^{\dagger}_{i,\alpha\sigma}c_{j,\alpha\sigma}$
$\displaystyle+U\sum_{i}\sum_{\alpha=1,2}n_{i,\alpha\uparrow}n_{i,\alpha\downarrow}-\mu\sum_{i}\sum_{\alpha=1,2}\left(n_{i,\alpha\uparrow}+n_{i,\alpha\downarrow}\right)$
$\displaystyle+(U-2J)\sum_{i,\sigma}n_{i,1\sigma}n_{i,2\bar{\sigma}}+(U-3J)\sum_{i,\sigma}n_{i,1\sigma}n_{i,2\sigma}$
(1)

on a bipartite lattice. Here, $\alpha$ and $\sigma$ denote the orbital and
spin indices respectively, $t_{\text{hop}}$ is the nearest neighbor hopping
amplitude between sites $i$ and $j$, $U$ is the intra-orbital Hubbard
repulsion, $J$ the Hund coupling, and $\mu$ the chemical potential. The
equilibrium ground state of the half-filled model at large $U$ is a Mott
insulator with antiferro spin ($J>0$) or orbital ($J<0$) order. Away from
half-filling, for $J>0$, orbital-singlet spin-triplet superconductivity
appears over a wide range of dopings between quarter filling ($n=1/4$) and
three-quarter filling ($n=3/4$) Hoshino ; Hoshino_2 ; enhanced_pairing , and
it is interesting to ask if a nonthermal manifestation of this order exists.

To create an effectively cold photo-doped state, we employ the NESS formalism
NESS_martin ; oneband_jiajun , where the system is weakly coupled to cold
fermion baths at each site. By choosing the densities of states (DOS)
$\rho_{b}(\omega)$ and chemical potentials $\mu_{b}$ of these baths
appropriately, we can inject (remove) electrons into (from) the upper (lower)
Hubbard band, and stabilize a photo-doped state. Here, we use fermion baths
with a semi-elliptic DOS with a coupling strength $\Gamma$. The hybridization
function $D_{b}(\omega)$ of these baths is given by $D_{b}(\omega)=$
$\Gamma\rho_{b}(\omega)=\Gamma\sqrt{W_{b}^{2}-(\omega-\omega_{b})^{2}}$, where
$W_{b}$ denotes the half-bandwidth and $\omega_{b}$ indicates the center of
the energy spectrum of the bath. We use four fermion baths with
$\omega_{b}=\pm U/2,\pm 3U/2$. The two baths with $\omega_{b}=\pm U/2$ are
chosen to coincide with the lower and upper Hubbard bands, to create a large
population of triply occupied sites (triplons) and singly occupied sites
(singlons). The other two baths, a completely empty one with $\omega_{b}=3U/2$
and a completely filled one at $\omega_{b}=-3U/2$, are added to prevent the
production of high-energy fully occupied (quadruplon) and empty (holon)
states, see Fig. 1(a).

The NESS solution of the two-orbital Hubbard model is obtained on an
infinitely connected Bethe lattice using non-equilibrium Dynamical Mean-Field
theory (DMFT) Aoki2014  and a non-crossing approximation (NCA) impurity solver
Eckstein2010 . To study uniform or staggered SC orders, we use the Nambu
Keldysh formalism to write the DMFT self-consistency condition as
$\Delta(t,t^{\prime})=t_{0}^{2}\gamma G(t,t^{\prime})\gamma+D(t,t^{\prime})$.
Here, $t_{0}=W/4$ is the quarter-bandwidth of the Bethe-lattice DOS, which is
used as the unit of energy. The hybridization function $\Delta$, impurity
Green’s function $G$ and hybridization function of the fermion baths
$D=\sum_{b}D_{b}$ can be represented by matrices in the Nambu basis, see
Supplemental Material (SM), while the choice of the $\gamma$ matrix determines
whether the solution corresponds to staggered or uniform pairing eta_philipp .

Figure 1: Spectra of the photo-doped two-orbital model with (a) $U=20$, $J=0$,
(b) $U=20$, $J=-1$, and (c) $U=20$, $J=1$. The occupied DOS is indicated by
the filled area. The four fermion baths and their occupation are indicated
below the spectra in (a) by the gray dashed lines and filled blue areas,
respectively. The insets in panel (a) show a zoom of the region indicated by
the gray box and a fit to the Fermi distribution function with $\beta=67.0$
(black dashed line).

In the calculations, we use $U=20$, $J=\pm 1$, $\Gamma=0.016$-$0.032$ and
$W_{b}=5$. A typical spectral function $A(\omega)$ for a photo-doped system
with triplon density $=0.474$ is shown in Fig. 1(a). The occupied density of
states $A^{<}(\omega)$ is indicated by the filled area. An effective
temperature $T_{\text{eff}}$ of the system can be defined by comparing the
ratio $A^{<}(\omega)/A(\omega)$ to a (shifted) Fermi distribution function
$f_{\beta}(\omega)$ for inverse temperature $\beta$. The inset of Fig. 1(a)
shows that $A^{<}(\omega)/A(\omega)$ can be very well fitted with
$f_{\beta=67.0}(\omega)$, which indicates that the triplons (and singlons) in
the non-equilibrium state have a well-defined and low effective temperature.
We also see sharp quasi-particle peaks at the edges of the Hubbard bands,
which are associated with the remaining doublon states. In equilibrium, at low
temperature, a single doublon state dominates and the energy separation
between the upper and lower Hubbard bands is roughly $U+J$ ($J>0$) or $U-5J$
($J<0$). In a photo-doped system, all doublon states contribute to the
spectrum, which exhibits peaks with a separation between $U-5J$ and $U+J$.
Hence for $J>0$ ($J<0$), the separation between the Hubbard bands decreases
(increases), see Figs. 1(b,c) and SM. A large enough gap ensures a slow
relaxation across the gap Sensarma2010 ; Eckstein2011  and is necessary for
the stability of the photo-doped state, which sets limits to values of $J$.

Results: In the two-orbital model we can define six types of SC orders,
depending on the spin and orbital symmetry. The three orbital-singlet spin-
triplet orders are measured by Hoshino

$\displaystyle p_{i,s\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\sum_{\alpha\sigma\alpha^{\prime}\sigma^{\prime}}c^{\dagger}_{i,\alpha\sigma}\epsilon_{\alpha\alpha^{\prime}}(\sigma^{\nu}\epsilon)_{\sigma\sigma^{\prime}}c^{\dagger}_{i,\alpha^{\prime}\sigma^{\prime}},$
(2)

and the three spin-singlet orbital-triplet orders by

$\displaystyle p_{i,o\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\sum_{\alpha\sigma\alpha^{\prime}\sigma^{\prime}}c^{\dagger}_{i,\alpha\sigma}(\sigma^{\nu}\epsilon)_{\alpha\alpha^{\prime}}\epsilon_{\sigma\sigma^{\prime}}c^{\dagger}_{i,\alpha^{\prime}\sigma^{\prime}},$
(3)

where $\sigma^{\nu}$ denotes the three Pauli matrices with $\nu=x,y,z$ and the
antisymmetric tensor is defined as $\epsilon=i\sigma^{y}$. In analogy to
antiferromagnetic order in magnetically ordered systems, a staggered SC order
($\eta$ order) can be defined, where $p_{i,\alpha}(\alpha=sx,sy,sz,ox,oy,oz)$
has opposite signs on the two sublattices $A$ and $B$, such that

$\displaystyle\eta^{+}_{i,\alpha}=\delta_{i}^{A/B}p_{i,\alpha},$ (4)

with $\delta_{i}^{A(B)}=1(-1)$ for $i$ on the $A(B)$ sublattice.

Without Hund coupling, the local Hamiltonian is fourfold degenerate (spin
degeneracy and orbital degeneracy) and at half-filling, for sufficiently large
$U$, each site is doubly occupied in the ground state. The photo-excitation
creates singlon-triplon pairs. In the extreme photo-doping limit, all doubly
occupied sites are converted to singlons and triplons, and if the fluctuations
between these states become large, we may expect the emergence of
superconductivity with order parameters given by Eqs. (2) or (3), similar to
the $\eta$ pairing states in the photo-doped single-band Hubbard model
eta_philipp ; oneband_jiajun ; Murakami2023 ; Rosch2008 . An important
distinction however is that the triplons and singlons in the two-orbital model
have spin and orbital degrees of freedom. The Hund coupling $J$ lifts the
degeneracy between the SC orders in the half-filled system and favors a spin
($J>0$) or orbital ($J<0$) moment at each site. It also stabilizes the SC
pairing, so that an $\eta$ SC order can be realized at sufficiently large
photo-doping.

Figure 2: SC order parameter $\eta^{x}_{sx}$ and inverse susceptibility
$1/\chi_{SC}$ as a function of the triplon density for $U=20$, $J=1$,
$\Gamma=0.02$ and $P_{\text{seed}}=10^{-4}$. While $\chi_{SC}$ for $\eta$
pairing is as high as $10^{3}$ for large triplon densities, for uniform
pairing it remains of order $1$.

Effective model: Physical properties of photo-doped states can be well
captured by an effective model derived from the Schrieffer-Wolff
transformation (SWT) or 2nd order perturbation theory oneband_jiajun . For
$J=0$ and half filling, the effective model for complete doping, where only
tripons and singlons exist, can be derived by projecting the transformed
Hamiltonian onto the singlon and triplon subspace (see SM) SWT ,

$\displaystyle H_{\text{eff}}$ $\displaystyle=$ $\displaystyle
H_{t_{\text{hop}}^{2}/U}^{st}+H_{t_{\text{hop}}^{2}/U}^{tt}+H_{t_{\text{hop}}^{2}/U}^{ss},$
(5)

where the low-energy effective interaction between triplons and singlons is
divided into three terms representing the singlon-triplon, triplon-triplon and
singlon-singlon interaction. Each term can be expressed in terms of spin,
orbital and $\eta$-spin operators as

$\displaystyle H_{t_{\text{hop}}^{2}/U}^{st}$ $\displaystyle=$
$\displaystyle-\frac{4t_{\text{hop}}^{2}}{U}\sum_{\left<ij\right>,\alpha}\left(\eta_{i\alpha}^{+}\eta_{j\alpha}^{-}+\eta_{i\alpha}^{-}\eta_{j\alpha}^{+}\right)$
(6)
$\displaystyle-\frac{4t_{\text{hop}}^{2}}{3U}\sum_{\left<ij\right>}\left(\frac{1}{2}+2\bm{s_{i}}.\bm{s_{j}}\right)\left(\frac{1}{2}+2\bm{\tau_{i}}.\bm{\tau_{j}}\right),$
$\displaystyle H_{t_{\text{hop}}^{2}/U}^{ss(tt)}$ $\displaystyle=$
$\displaystyle\frac{2t_{\text{hop}}^{2}}{U}\sum_{\left<ij\right>}\left(\frac{1}{2}+2\bm{s_{i}}.\bm{s_{j}}\right)\left(\frac{1}{2}+2\bm{\tau_{i}}.\bm{\tau_{j}}\right),$
(7)

where $\eta^{+}_{i\alpha}$ is defined in Eq. (4), $\eta^{-}_{i\alpha}$ is the
Hermitian conjugate $(\eta^{+}_{i\alpha})^{\dagger}$,
$\bm{s}_{i}=\sum_{\alpha}c^{\dagger}_{i,\alpha\sigma}\bm{\sigma}_{\sigma\sigma^{\prime}}c_{i,\alpha\sigma^{\prime}}$
the spin moment, and
$\bm{\tau}_{i}=\sum_{\sigma}c^{\dagger}_{i,\alpha\sigma}\bm{\sigma}_{\alpha\alpha^{\prime}}c_{i,\alpha^{\prime}\sigma}$
the orbital moment. As can be seen from Eqs. (6) and (7), in addition to
$\eta$-spin interaction terms, there are also spin interactions, orbital
interactions and spin-orbital composite terms. The temperature scales for all
types of $\eta$ orders are set by the coefficients which are $\sim
t_{\text{hop}}^{2}/U$. Nonzero Hund coupling lifts the degeneracy between
these $\eta$ orders, because it affects the energies of the doublon states
involved in the second order hopping processes. The exchange coupling for the
$\eta$ pairing channels $\alpha$ (the first term in Eq. (6)) is modified as
$t_{\text{hop}}^{2}/U\rightarrow t_{\text{hop}}^{2}/\tilde{U}_{\alpha}$, where
$\tilde{U}_{\alpha}$ is a function of $U$ and $J$. For $J>0$, $\tilde{U}_{sx}$
takes the smallest value and favors spin-triplet $\eta_{sx}$ superconducting
order (see SM).

A mean-field (MF) decoupling may be applied to the effective model to obtain a
MF Hamiltonian (see SM). This MF calculation suggests that an $\eta$ SC state
can be achieved at low $T_{\text{eff}}$, and that the Hund coupling enhances
spin-triplet or orbital-triplet superconductivity, depending on $J$. However,
compared to the non-equilibrium DMFT results, MF overestimates the SC
transition temperature (red triangle in Fig. 3). This is because MF neglects
temporal fluctuations. Moreover, the effective model is valid only in the
large-$U$ limit and the DMFT solution is affected by the bath couplings.

Figure 3: Phase diagram of the two-orbital model with $U=20$, $J=1$ and
$\Gamma=0.02$ in the space of $\beta_{\text{eff}}$ and triplon density. The
AFM phase appears in the undoped and weakly photo-doped region, while the
$\eta$ SC state is stabilized in the high photo-doping region. The red
triangle marks the MF transition temperature at triplon density $0.5$. The
right panel shows the color map of the order parameter Re[$\eta^{+}_{sx}$] in
the high photo-doping region.

DMFT results: To detect a SC state, we apply a small seed field
$P_{\text{seed}}$ and couple it to the SC order we are interested in. For
example, if we look for the $p_{sx}$ SC order, we add the term
$H_{\text{pair}}=P_{\text{seed}}\sum_{i}(p_{i,sx}+h.c.)$ to the Hamiltonian
(1) and measure the order parameter $\left<(p_{sx}+h.c.)\right>$. For a $\eta$
SC state, we use a staggered local seed field with opposite signs on the two
sublattices and measure the SC order parameter
$\left<\eta^{x}_{\alpha}\right>$ and SC susceptibility $\chi_{SC}$ as

$\displaystyle\left<\eta^{x}_{\alpha}\right>=\text{Re}[\left<\eta^{+}_{\alpha}\right>]=\frac{\left<\eta^{+}_{\alpha}+\eta^{-}_{\alpha}\right>}{2},\quad\chi_{SC}=\frac{\left<\eta^{x}_{\alpha}\right>}{P_{\text{seed}}}.$
(8)

The effective temperature and band filling can be controlled by tuning the
temperature and chemical potentials of the fermion baths, which allows us to
map out the phase diagram of the non-equilibrium photo-doped state in the
space of triplon density and $T_{\text{eff}}$. For $J=0$, all $\eta$ pairing
phases are degenerate and the highest transition temperature is
$T_{\text{eff}}=1/62.6$. A nonzero positive (negative) Hund coupling favors
the spin triplet (orbital triplet) $\eta$ phases, and the corresponding
transition temperatures increase.

Figure 2 plots the SC susceptibility $\chi_{SC}$ and order parameter
$\left<\eta^{x}_{sx}\right>$ as a function of the triplon density for $U=20$
and $J=1$. The $\eta$ SC state appears at high triplon density ($\gtrsim
0.46$), as indicated by the large increase in $\chi_{SC}$. We observe an order
parameter of the order of $10^{-1}$, which for $P_{\text{seed}}=10^{-4}$
corresponds to $\chi_{SC}\sim 10^{3}$. Importantly, this order parameter is
almost unchanged if the seed is reduced to $P_{\text{seed}}=10^{-5}$, which
demonstrates that the order is due to a spontaneous symmetry breaking (see
SM). For comparison, we also plot in Fig. 2 the susceptibility for the uniform
order $p_{sx}$, which remains small.

The $\beta_{\text{eff}}$ vs triplon density phase diagram for the same
interaction parameters is shown in Fig. 3. At half filling, without any photo-
doping, we have an antiferromagnetic (AFM) state, which is quickly suppressed
with increasing triplon density. Over a wide region of intermediate photo-
dopings, the system is in a normal paramagnetic state. As we approach the
extreme doping limit, an $\eta$ spin-triplet SC state is realized below a
sufficiently low effective temperature ($\beta_{\text{eff}}\gtrsim 52.0$).

As mentioned earlier, in the two-orbital system there can be various SC orders
with different spin and orbital symmetries. Thus, in Fig. 4, we also show the
phase diagram for fixed $U=20$ and $\beta_{\text{eff}}=60$ in the space
spanned by the triplon density and Hund coupling. The data for the orbital-
triplet order $\eta^{x}_{ox}$ and spin-triplet order $\eta^{x}_{sx}$, on which
this phase diagram is based, are shown in the SM. Increasing $|J|$ favors the
$\eta$ SC order in both cases with a dominant $\eta^{x}_{sx}$ pairing for
positive $J$. Insights into the $J<0$ case can be obtained by noting that the
$J>0$ model can be mapped (qualitatively) to the $J<0$ model by the
transformation $c_{i,1\downarrow}\rightarrow c_{i,2\uparrow}$ and
$c_{i,2\uparrow}\rightarrow c_{i,1\downarrow}$ Hoshino_2 . This exchanges the
roles of spin and orbital and it is thus natural that $\eta^{x}_{ox}$ pairing
becomes dominant for negative $J$. The smaller region of $\eta^{x}_{ox}$
pairing compared to $\eta^{x}_{sx}$ is due to the fact that the relevant
interaction $\tilde{U}_{\alpha}$ is larger for $J<0$ than for $J>0$, which
reduces the critical temperature ($\propto
t_{\text{hop}}^{2}/\tilde{U}_{\alpha}$) for $J<0$.

Figure 4: Nonquilibrium phase diagram in the space of triplon density and Hund
coupling $J$, for $U=20$ and $\beta_{\text{eff}}=60$.

Optical conductivity: We finally calculate the optical conductivity and study
how the optical response changes as we approach the SC phase from the normal
phase. The steady state optical conductivity $\sigma(\omega)$ is calculated
for the Bethe lattice from the local Green’s function, following Ref.
oneband_jiajun .

Figure 5: Optical conductivity results for $U=20$, $J=1$ and for triplon
density 0.473. (a) Order parameter $\eta^{x}_{sx}$ as a function of effective
temperature $T_{\text{eff}}$, (b) imaginary part and (c) real part of the
optical conductivity as a function of $\omega$. The different colors in (b)
and (c) correspond to the corresponding data points in (a). The inset of (b)
shows the log-log plot of (shifted) Im$[\sigma(\omega)]$. (d) Square root of
the SC phase stiffness $D$ as a function of the order parameter
$\eta^{x}_{sx}$.

Figure 5(a) shows the order parameter $\eta^{x}_{sx}$ as a fuction of
$T_{\text{eff}}$ for the triplon density $0.473$, $U=20$ and $J=1$. The SC
transition temperature is $T_{c}\sim 0.019$. In both the normal and SC state,
in addition to a narrow Drude peak, which we associate with the doublon
population, the real part of $\sigma(\omega)$ shows a range of negative
spectral weight around $\omega\sim 15$ (panel (c)). The latter comes from
triplon-singlon recombination processes in the population inverted non-
equilibrium system, which are associated with local energy changes in the
range from $U-5J$ to $U+J$. The imaginary part, however, shows a distinct low-
frequency behavior in the SC and non-SC states (panel (b)). In the non-SC
state (blue), Im$[\sigma(\omega)]$ decreases to $0$ as $\omega\rightarrow 0$.
The onset of superconductivity (green dot in panel (a)) is marked by a
diverging Im[$\sigma(\omega)$] as $\omega\rightarrow 0$, as shown by the green
curve in panel (b). All the states which have a lower $T_{\text{eff}}$ than
this state, and hence are superconducting, exhibit a $1/\omega$ divergence at
low $\omega$. The inset of Fig. 5(b) shows the log-log plot of appropriately
shifted Im[$\sigma(\omega)$] and we can clearly observe the $1/\omega$
behavior below $\omega\approx 10^{-1}$ for the SC states (red and yellow). As
$T_{\text{eff}}$ increases, the energy region over which the $1/\omega$
scaling can be observed decreases and eventually disappears in the non-SC
state (blue). We note that all the curves in Fig. 5(b) show an approximate
power-law decay for $\omega>10^{-1}$. However, the normal state behavior is
different from a $1/\omega$ decay, as can be seen from the inset. The
$1/\omega$ scaling at low frequencies thus provides evidence for the
realization of a nonthermal SC state.

The $1/\omega$ divergence of Im[$\sigma(\omega)$] implies that we can define
the SC phase stiffness $D$ from the relation

$\lim\limits_{\omega\to 0}\text{Im}[\sigma(\omega)]=\frac{D}{\omega}.$

In Fig. 5(d) we plot $\sqrt{D}$ as a function of the order parameter
$\eta^{x}_{sx}$ for fixed $U=20$, $J=1$ and observe a linear scaling. This
$D\propto|\eta|^{2}/U$ scaling is reminiscent of the two-fluid model and holds
in the entire SC region oneband_jiajun .

Conclusions: We demonstrated that upon photo-doping, staggered $\eta$ SC
states can be realized in a small region near triplon density $0.5$ in the
half-filled two-orbital repulsive Hubbard model. Unlike in the single band
model, the relevant local states have spin and orbital degrees of freedom,
which potentially allows to realize different SC states. We showed that a
positive Hund coupling enhances orbital-singlet spin-triplet $\eta$ SC pairing
whereas a negative Hund coupling favors orbital-triplet spin-singlet $\eta$
pairing. Due to the large $U$ needed to stabilize the non-equilibrium states,
the highest $T_{\text{eff}}$ for realizing multi-orbital $\eta$ pairing states
are about an order of magnitue lower than in the single band case
oneband_jiajun . However, in the realistic case where the energy unit is
$\sim$0.5 eV, the SC state can still be realized up to high effective
temperatures of about 100 K. It will be interesting to explore similar $\eta$
pairing states in other photo-doped multi-orbital systems and to check if SC
states coexisting with spin or orbital orders can be realized.

Acknowledgements We thank J. Li and M. Eckstein for providing their NESS DMFT
solver and J. Li and S. Hoshino for helpful discussions. The calculations have
been run on the beo05 and beo06 clusters at the University of Fribourg. We
acknowledge support from ERC Consolidator Grant No. 724103 and from SNSF Grant
No. 200021-196966.

## References

  * (1) D. N. Basov, R. D. Averitt, and D. Hsieh , Nat. Mat. 16, 1077 (2017).
  * (2) C. Giannetti, M. Capone, D. Fausti, M. Fabrizio, F. Parmigiani, and D. Mihailovic, Advances in Physics 65, 58-238 (2016).
  * (3) T. Oka and S. Kitamura, Annual Review of Condensed Matter Physics 10, 387 (2019).
  * (4) A. dela Torre, D. M. Kennes, M. Claassen, S. Gerber, J. W. McIver, and M. A. Sentef, Rev. Mod. Phys. 93, 041002 (2021).
  * (5) H. Ichikawa, S. Nozawa, T. Sato, A. Tomita, K. Ichiyanagi, M. Chollet, L. Guerin, N. Dean, A. Cavalleri, S.-i. Adachi, et al., Nat. Mater. 10, 101 (2011).
  * (6) L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D. Svetin, S. Brazovskii, and D. Mihailovic, Science 344, 177 (2014).
  * (7) S. Mor, M. Herzog, D. Golez, P. Werner, M. Eckstein, N. Katayama, M. Nohara, H. Takagi, T. Mizokawa, C. Monney, and J. Stahler, Phys. Rev. Lett. 119, 086401 (2017).
  * (8) E. Dagotto, Science 309, 257 (2005).
  * (9) T. Li, A. Patz, L. Mouchliadis, J. Yan, T. A. Lograsso, I. E. Perakis, and J. Wang, Nature 496, 69 (2013).
  * (10) R. Mikhaylovskiy, E. Hendry, A. Secchi, J. H. Mentink, M. Eckstein, A. Wu, R. Pisarev, V. Kruglyak, M. Katsnel-son, T. Rasing, et al., Nat. Commun. 6, 8190 (2015).
  * (11) A. Kirilyuk, A. V. Kimel, and T. Rasing, Rev. Mod. Phys. 82, 2731 (2010).
  * (12) D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, Science 331, 189 (2011).
  * (13) M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S. Lupi, P. Di Pietro, D. Pontiroli, M. Ricco, S. R. Clark, et al., Nature 530, 461 (2016).
  * (14) D. Basov, R. Averitt, and D. Hsieh, Nat. Mater. 16, 1077 (2017).
  * (15) P. Werner, N. Tsuji, and M. Eckstein, Phys. Rev. B 86, 205101 (2012).
  * (16) A. Ono and S. Ishihara, Phys. Rev. Lett. 119, 207202 (2017).
  * (17) J. Li, H. U. R. Strand, P. Werner and M. Eckstein , Nat. Com. 9, 4581 (2018).
  * (18) P. Werner, and Y. Murakami, Phys. Rev. B 102, 241103(R) (2020).
  * (19) T. Kaneko, T. Shirakawa, S. Sorella, and S. Yunoki. Phys. Rev. Lett. 122, 077002 (2019).
  * (20) P. Werner, J. Li, D. Golez, and M. Eckstein, Phys. Rev. B 100, 155130 (2019).
  * (21) J. Tindall, B. Buca, J. R. Coulthard, and D. Jaksch, Phys. Rev. Lett. 123, 030603 (2019).
  * (22) J. Li, D. Golez, P. Werner, and M. Eckstein, Phys. Rev. B 102, 165136 (2020).
  * (23) P. Werner and Y. Murakami, Phys. Rev. B 104, L201101 (2021).
  * (24) Y. Murakami, S. Takayoshi, T. Kaneko, A. M. Läuchli, and P. Werner, Phys. Rev. Lett. 130, 106501 (2023).
  * (25) A. Takahashi, H. Gomi, and M. Aihara, Phys. Rev. B 66, 115103 (2002).
  * (26) A. Rosch, D. Rasch, B. Binz, and M. Vojta, Phys. Rev. Lett. 101, 265301 (2008).
  * (27) Y. Murakami, S. Takayoshi, T. Kaneko, Z. Sun, D. Golež, A. J. Millis, and P. Werner, Communications Physics 5, 23 (2022).
  * (28) T. Oka and H. Aoki. Phys. Rev. B 79, 081406(R) (2009).
  * (29) Y. H. Wang, H. Steinberg, P. Jarillo-Herrero, and N. Gedik, Science 342, 453 (2013).
  * (30) S. Iwai, M. Ono, A. Maeda, H. Matsuzaki, H. Kishida, H. Okamoto, and Y. Tokura, Phys. Rev. Lett. 91, 057401 (2003).
  * (31) R. Sensarma, D. Pekker, E. Altman, E. Demler, N. Strohmaier, D. Greif, R. Jördens, L. Tarruell, H. Moritz, and T. Esslinger, Phys. Rev. B 82, 224302 (2010).
  * (32) M. Eckstein and P. Werner, Phys. Rev. B 84, 035122 (2011).
  * (33) N. Strohmaier, D. Greif, R. Jordens, L. Tarruell, H. Moritz, T. Esslinger, R. Sensarma, D. Pekker, E. Alt- man, and E. Demler, Phys. Rev. Lett. 104, 080401 (2010).
  * (34) Z. Lenarcic and P. Prelovsek, Phys. Rev. B 90, 235136 (2014).
  * (35) N. Maeshima, and K. Yonemitsu, J. Phys. Soc. Jpn. 74, 2671 (2005).
  * (36) D. Golez, J. Bonca, L. Vidmar, and S. A. Trugman, Phys. Rev. Lett. 109, 236402 (2012).
  * (37) K. Shinjo and T. Tohyama, Phys. Rev. B 96, 195141 (2017).
  * (38) M. Sandri and M. Fabrizio, Phys. Rev. B 88, 165113 (2013).
  * (39) M. Eckstein and P. Werner, Phys. Rev. Lett. 110, 126401 (2013).
  * (40) N. Tancogne-Dejean, M. A. Sentef, and A. Rubio, Phys. Rev. Lett. 121, 097402 (2018).
  * (41) J. Rincon, E. Dagotto, A. E. Feiguin, Phys. Rev. B 97, 235104 (2018).
  * (42) K. Zawadzki and A. E. Feiguin, Phys. Rev. B 100, 195124 (2019).
  * (43) M. Schüler, U. De Giovannini, H. Hübener, A. Rubio, M. A. Sentef, T. P. Devereaux, and P. Werner, Phys. Rev. X 10, 041013 (2020).
  * (44) N. Bittner, D. Golez, M. Casula, and P. Werner, Phys. Rev. B 104, 115138 (2021).
  * (45) E. Perfetto, Y. Pavlyukh, and G. Stefanucci, Phys. Rev. Lett. 128, 016801 (2022).
  * (46) P. Werner, M. Eckstein, M. Müller, and G. Refael, Nat. Comm. 10, 5556 (2019).
  * (47) J. Li and M. Eckstein, Phys. Rev. B 103, 045133 (2021).
  * (48) J. Li, M. Müller, A. J. Kim, A. Läeuchli, and P. Werner, Phys. Rev. B 107, 205115 (2023).
  * (49) K. I. Kugel and D. I. Khomskii, Sov. Phys. Usp. 25, 231 (1992).
  * (50) G. Khaliullin and V. Oudovenko, Phys. Rev. B 56, R14243(R) (1997).
  * (51) P. Werner, H. U. R. Strand, S. Hoshino, Y. Murakami, and M. Eckstein, Phys. Rev. B 97, 165119 (2018).
  * (52) S. Hoshino and P. Werner, Phys. Rev. B 93, 155161 (2016).
  * (53) K. Steiner, S. Hoshino, Y. Nomura, and P. Werner, Phys. Rev. B 94, 075107 (2016).
  * (54) H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Rev. Mod. Phys. 86, 779 (2014).
  * (55) M. Eckstein and P. Werner, Phys. Rev. B 82, 115115 (2010).
  * (56) S.-S. B. Lee, J. v. Delft, and A. Weichselbaum, Phys. Rev. B 96, 245106 (2017).