A Clique-Based Method for Improving Motif Scanning Accuracy, v3.01

Braslav Rabar, Keti Nižetić and Pavle Goldstein

University of Zagreb, Faculty of Science, Mathematics Department

## Abstract

### Background

Motif scanning is a very common method in bioinformatics. Its objective is to
detect motifs of sufficient similarity to the query, which is then used to
determine familiy membership, or structural or functional features or
assignments. Considering a variety of uses, accuracy of motif scanning
procedures is of great importance.

### Results

We present a new approach for improving motif scanning accuracy, based on
analysis of in-between similarity. Given a set of motifs obtained from a
scanning process, we construct an associated weighted graph. We also compute
the expected weight of an edge in such a graph. It turns out that restricting
results to the maximal clique in the graph, computed with respect to the
expected weight, greatly increases precision, hence improves accuracy of the
scan. We tested the method on an ungapped motif-characterized protein family
from five plant proteomes. The method was applied to three iterative motif
scanners - PSI-BLAST, JackHMMer and IGLOSS - with very good results.

### Conclusions

We presented a method for improving protein motif scanning accuracy, and have
successfully applied it in several situations. The method has wider
implications, for general pattern recognition and feature extraction
strategies, as long as one can determine the expected similarity between
objects under consideration.

##  1 Background

Motif scanning - or local similarity search - is a very important part of
sequence analysis. It can be used for various purposes - protein family
assignment ([2]), secondary structure prediction ([8]) and similar. Motif
scanning methods are typically based on a local alignment algorithm - Smith-
Waterman or Viterbi algorithm - with various modifications added, such as
approximations and variations in the scoring function or model building.

Motif scanning procedures normally take an instance of a motif - or even a
profile - as an input - that is a query \- and search for similar patterns in
a set of sequences. The output consists of a set of sufficiently similar
matches, and these form a set of positives, or - what we call it here - the
response.

In this paper, we are concerned with accuracy of motif scanning procedures.
Namely, the aim of motif scanning applications is to detect as many
significant motifs as possible - these are the true positives, while keeping
the number of wrong assignments - or false positives \- to the minimum. In
other words, accuracy is measured by how closely the response matches the set
of biologically relevant sequences in the sample.

Typically, motif scanning methods are based on a scoring function - usually a
log-likelihood ratio - and the response is generated in two steps: first, by
ranking all elements of the sample by their similarity to the query, and,
second, by considering all candidates with a similarity score above a certain
threshold. Both steps can be a source of errors - an inaccurate ranking
scheme, with a low threshold, will generate a huge response, causing a large
type I-like error, whereas a high threshold, while testing positive on
strongest examples, might miss the candidates with a slightly weaker signal.
Various applications deal with these problems in various ways - for example,
JackHMMer ([6]) uses - effectively - several scoring functions, each with its
own ranking and threshold, while IGLOSS ([9]) uses detailed parameter
estimation to minimize both issues.

In this paper, we explore an alternative approach for improving accuracy,
based on pairwise similarity. Namely, given a large response, containing -
presumably - a large percentage of false positives, we search for subsets in
which each pair of elements is sufficiently similar. We then consider the
largest of these subsets as the new, modified response. Since two true
positives are more likely to be similar than a true and a false positive or
two false positives, this is a sensible strategy, and it also turns out to be
very robust. Furthermore, we determine what “sufficiently similar” means -
that is, we compute the appropriate similarity threshold, in terms of the
expected conservation and length of the motif. The algorithm is presented in
the framework of graph theory, where the new response is obtained as the
maximal clique in a certain graph, that is, in turn, derived from the response
graph, modified by applying the similarity threshold. As already mentioned,
using the maximal clique as the new response greatly increases accuracy of the
search.

##  2 Methods

###  2.1 Response Graph, Derived Graph and Maximal Clique Algorithm

The starting point for our analysis is a set of positives from a motif-
scanning process. We assume that our motifs are ungapped, so this set, called
response, is just a collection of k-mers. We form an undirected, weighted
graph $\Gamma=(V,E)$, with vertices given by the elements of the response, and
weights $w(e,f)$, $e,f\in V$, given by

$w(e,f)=\sum_{i=1}^{k}B(e_{i},f_{i}),$

where $B(\cdot,\cdot)$ is the BLOSUM50-score for corresponding amino acids.

We now form an undirected $\\{0,1\\}$-graph $\bar{\Gamma}=(V,\bar{E})$ using
the expected weight $st$ from Section 2.3: for $e,f\in V$ and the unordered
pair $\\{e,f\\}$ we have

$\\{e,f\\}\in\bar{E}\Leftrightarrow w(e,f)-st>0.$

Put simply, two vertices in $\bar{\Gamma}$ are connected if and only if their
similarity score is above the threshold.

It is now straightforward to apply a maximal-clique algorithm to
$\bar{\Gamma}$. We used the standard Bron-Kerbosh algorithm ([3]), implemented
in Python.

###  2.2 Accuracy Measures

In this subsection, we establish notation and define relevant accuracy
measures, and we closely follow ([9, Supplementary material, Section 5]) in
the presentation. Sequences in the sample that have been annotated as
belonging to test families are marked as condition positive and their number
is denoted as $|CP|$, while the rest are marked as condition negative (CN).
Now, each application or combination of applications under consideration
produces - for a specified similarity level - a list of hits, and their
respective sequences are denoted as positive (P) - with the rest of the sample
being negative (N) - while $|P|$ and $|N|$ denote the corresponding set sizes.
We then have true positives (TP) and false positives (FP) as

$TP=P\cap CP,\ FP=P\cap CN,$ (1)

and likewise for true negatives (TN) and false negatives (FN)

$TN=N\cap CN,\ FN=N\cap CP.$ (2)

The usual way of assessing diagnostic ability of an application would be to
compare sensitivity or true positive rate

$TPR=|TP|/|CP|$ (3)

and false positive rate

$FPR=|FP|/|CN|.$ (4)

However, in the present context, there is a serious imbalance between the
sizes of the condition positive and condition negative sets: CN is several
orders of magnitude larger than CP, so, for any reasonable test outcome,
$|FPR|$ will be close to $0$. Consequently, we consider precision or positive
predictive value

$PPV=|TP|/|P|,$ (5)

and use $PPV$ and $TPR$ as accuracy measures. Finally, we combine these two by
considering their harmonic mean, called the F1-score, hence

$F1=2\cdot\frac{PPV\cdot TPR}{PPV+TPR},$ (6)

and plot threshold-F1 diagrams to assess accuracy.

###  2.3 Similarity Threshold

We compute the similarity threshold as the expected BLOSUM ([7]) score of two
$k$-mers, sampled from a certain distribution. Namely, given
$x=(x_{1},\ldots,x_{k})$ and $y=(y_{1},\ldots,y_{k})$, let

$s(x,y)=\sum_{i}B(x_{i},y_{i})$

be their BLOSUM score and $\mathbb{E}[s(x,y)]$ the expected value. Then,
somewhat informally,

$\mathbb{E}[s(x,y)]=\sum\mathbb{E}[B(x_{i},y_{i})],$

and, averaging over the whole length of the motif,

$\mathbb{E}[s(x,y)]=k\mathbb{E}[B(x_{0},y_{0})],$

for some “average” amino acids $x_{0}$ and $y_{0}$. Now, let

$e_{i}=(0,\ldots,1,0,\ldots,0),$

with $1$ on $i$-th position, let $\alpha\in(0,1)$, and let $bg=bg(i)$ be the
average distribution of amino acids. Let

$f_{i}=\alpha\cdot e_{i}+(1-\alpha)bg$ (7)

be an $e_{i}-bg$-mixture of distributions. Here, $\alpha$ is a “conservation
parameter”, representing the percentage of the dominant amino acid $i$ in an
alignment column sampled from $f_{i}$. Then

$\sum_{j,k=1}^{20}B(j,k)f_{i}(j)f_{i}(k)$

is the expected BLOSUM score for two amino acids sampled from $f_{i}$. Setting
$\alpha=0.68$ and averaging over the distribution $bg$, we get

$\sum_{i=1}^{20}bg(i)\sum_{j,k=1}^{20}B(j,k)f_{i}(j)f_{i}(k)=2.522,$ (8)

and, for a motif length $k$, we take the similarity threshold $st$ to be

$st=k\cdot 2.5$ (9)

We further comment on this in Section 3.2.3.

##  3 Results and Discussion

###  3.1 Tests and Results

In order to test the method, we applied it to responses from three iterative
motif scanners - PSI-BLAST (PB) ([1]), JackHMMER (JH) ([6]) and IGLOSS (IG)
([9]) - and compared the maximal clique with the original response. As in
([9]), scanners were applied to five plant proteomes - Arabidopsis thaliana
(AT, v. TAIR9), Oryza sativa (OS, v. MSU v7), Solanum tuberosum (ST, v.
ITAG1), Solanum lycopersicum (SL, v. ITAG2.3) and Beta vulgaris (BV, v.
KWS2320) - where we searched for members of an extensively studied, motif
characterized protein family - GDSL lipases.

GDSL lipases belong to lipid hydrolyzing enzymes that exhibit a GDSL motif.
Proteins in this family display fairly low overall sequence similarity, but
are reasonably well described by the presence of conserved residues in four
conserved blocks (I, II, III, and V) ([11]). Block I contains the main
characteristic motif (PROSITE:PS01098) ([10]) from which the main search query
of 10 amino acids was constructed. As in [9], the condition positive set was
determined by processing the information from GoMapMan resource [12].

Altogether, we performed approximately $900$ tests. We used three search
queries:

  * •

FVFGDSLSDA - consensus query, defined above

  * •

FVFnDSLSDA - a single mutation, at a highly conserved site

  * •

vfFGDSLSDn - three substitutions

This was done with all three scanners, for approximately $20$ threshold levels
each, and all five proteomes.

The average gain in the F1-score was around $0.20$ in the most interesting
threshold region (around 1/3 of the x-axis). However, gains can sometimes be
spectacular - Table 1 shows results of scans where we more than doubled
F1-score, albeit at a fairly low threshold level. Another feature of our tests
that we should comment on are peaks-and-troughs present in the Figure 3.
Namely, the corresponding responses are diffuse and inhomogeneous, as a result
of a “wrong” query, so the clique algorithm hardly improves PPV at certain
threshold levels - hence a “jagged” PPV-curve. We further analyzed this
situation, and detected features that indicate that the clique recognizes this
“wrongness”. We discuss this in Section 3.2.4.

Figures in the last section were obtained by merging the results from all five
proteomes. Here, responses are matched by their size (i.e. the number of
positives), and the average threshold - scale (for IG) or negative logarithm
of the e-value (for PB and JH) - was assigned to the x-axis, with averaged
PPV, TPR or F1-score on the y-axis. For the tables, responses were again
matched by their sizes, with the last column reporting cumulative results for
positives and true positives, and the average threshold.

###  3.2 Discussion

####  3.2.1 Overview

Put somewhat abstractly, the aim of our method is to detect - using pairwise
comparison - an optimal subset of true positives in a (fairly) large set of
hits (positives). This is achieved with the help of a pairwise similarity
threshold and a maximal clique algorithm. We have developed this concept in
the context of ungapped motif scanning - hence the present method - and tested
it in conjunction with three different iterative motif scanners. As stated
before, these principles can be applied more generally, provided suitable
assumptions are fulfilled.

####  3.2.2 Robustness

Before we discuss the way we computed the similarity threshold - in the next
section - let us first note that the threshold is rather robust. Namely, as
seen from the Figure 10(a), $10\%$ change in the level of the threshold
produces only minor changes in terms of accuracy (measured by the F1-score).
Hence, gradual changes in the threshold level will affect accuracy gradually.
On the other hand, this is in contrast to the iterative approach - as
mentioned in the Introduction - where small changes of the threshold (e-value
or scale) might produce very different results. This is due to the different
nature of two approaches - there are many more comparisons carried out in the
pairwise case, so the effect of the change is smoothed out. Furthermore, to
phrase this in the graph theory framework: the pairwise threshold acts on the
edges of the response graph, while we are looking for a suitable subset of the
vertices; so, while a small change in the threshold level might affect many
edges, that will - eventually - add or remove only a couple of vertices from
the maximal clique.

####  3.2.3 Similarity Threshold Comments

The aim of the similarity threshold is to distinguish between true and false
positives by means of the pairwise score. Namely, protein motifs are
chracterized by a specific substitution pattern, so it is to be expected that
the pairwise score between two true positives will be higher than the one
between true and false positives, let alone two false positives (note that
this naturally leads to the maximal clique approach). However, knowledge of
the substitution pattern amounts to a detailed description of the motif under
consideration, in which case the type of analysis that we study here - single-
query, iterative approximation - becomes almost redundant.

Hence, we derived the threshold a priori, as an abstract, average score, and
dependent on a single parameter $\alpha$, so-called conservation coefficient.
More precisely, the threshold was obtained as the expected BLOSUM score of two
$k$-mers, sampled from an $\alpha$-convex combination of distributions, as
well as assuming average (i.e. background) distribution of amino acids across
the length $k$. The parameter $\alpha$ should be thought of as the relative
frequency of the dominant amino acid in the hypothetical motif profile,
averaged across the length $k$. Incidently, in five proteomes that we
considered, profile conservation varied from $68\%$ to $72\%$.

There are other strategies to obtain the threshold, as long as appropriate
conservation is maintained. For example, one could use a suitable power of the
PAM matrix ([4]) instead of distributions $f_{i}$ in the Equation 7. So,
taking PAM120 - which yields the average diagonal value around $2/3$ \- and
repeating the procedure from Section 2.3, gives the value $st=2.58$ \-
approximately the same. Likewise, using uniform distribution - instead of
background - in the Equation 8, one obtains $st=2.4$. Furthermore, other
similarity measures, rather than the BLOSUM50 matrix, could be used - other
BLOSUM matrices, or even other, non-BLOSUM substitution matrices. However,
that would involve setting up a new response graph and a new threshold, in
parallel. Considering the comments above, one should expect marginal changes,
or no changes at all.

Finally, we should mention that, although the threshold was obtained as an
expected score, it is used as the “minimal allowed” score. This is because
assumptions for the conservation and, especially, amino acid composition are
rather weak. It is possible to tighten these assumptions, and then set the
“minimal score” to, say, $-2\sigma$ (i.e. two standard deviations) from the
mean. However, tightening the assumptions would again amount to a description
of the motif under consideration - hence, not an a priori approach - so we
decided not to explore this further.

####  3.2.4 Peaks-and-Troughs

As mentioned in Section 3.2.4, peaks-and-troughs in the Figure 3 were caused
by the incorrect query, resulting in an inhomogenous, diffuse response. This,
in turn, produced the maximal clique containing a large number of false
positives, yielding a fairly low PPV for some threshold levels. Obviously,
this is in contrast with Figure 1, where the “right” - i.e. consensus - query
produces more homogenous response, and the clique yields consistent
improvement. We analyzed this a bit futher, in order to detect features that
might differentiate between these two situations.

First of all, note that most of the scan descriptions used above - “low PPV,
varying F1-score”, and figures such as Figure 3 and Figure 1 \- are available
only a posteriori. In other words, this information is not available in an
exploratory setting, where one wants to asses the validity of the response
without knowing the desired outcome. So, we should be looking for some other -
“a priori” \- features.

Homogeneity of the response is a possible candidate, and it will distinguish
between these two sets of scans, but, again, homogeneity is a feature that is
best measured a posteriori, when some information regarding the variability of
the condition positive set is available. So, we looked at stability, that is,
relationship between responses and their respective cliques, for neighbouring
threshold levels.

A sequence of threshold levels, from low to high, should - in principle -
produce a descending family of responses. This is, clearly, the case in a non-
iterative setting, where a simple scan results in a fixed ranking scheme - an
ordering of the sample with respect to the similarity to the query. In an
iterative situation, the similarity function is being optimized, which might
produce a different ranking scheme from one iteration to another, and result
in a different ranking scheme from one threshold to another. Consequently, a
smaller response - produced with a higher threshold - might not be the subset
of a larger one, obtained with a lower threshold. Furthermore, a significant
deviation from this stability principle indicates, in general, problems with
either the sample or the query.

Let us analyze from this point of view the series of scans from Figure 3 (we
will focus on Arabidopsis thaliana here; some of the scan results - for
neighbouring e-values - are presented in Tables 4 and 3; for a complete set of
tables, consult the server web-page): as already mentioned, we have scanned
with a non-consensus query; quite surprisingly, responses have shown to be
rather stable, with sizes $98$ and $81$, and the size of the intersection
$80$; however, the corresponding cliques - with sizes $26$ and $24$ \- have a
single element in the intersection, which should be considered as a
significant deviation from the stability principle. Hence, our approach
provides another method to assess validity of a scan, with the clique as a new
stability criterion.

####  3.2.5 Pairwise Similarity vs EM-algorithm

All the iterative motif scanners that we combined our clique-method with use
some form of the expectation maximization algorithm (EM-algorthm) to find the
“optimum” - the optimal set of positives, for a given significance level. On
the other hand, the maximal clique algorithm also provides the optimal
solution for the set of positives - the maximal clique itself. A natural
question arises: are these two optima the same? More precisely, given the
right parameters - similarity and significance threshold - will these two
approaches return the same, or very similar, response?

The answer appears to be yes - optimal solutions will be more-or-less the same
for the right choice of parameters, with a couple of outliers added or
subtracted. This can already be inferred from figures in the next Section,
where we see stability in TPR across all threshold ranges, and F1-convergence,
as the threshold becomes higher. In the opposite direction, iterative scanners
will, invariably, accept the maximal clique as the optimal solution, again
with minor changes.

All this means that these two approaches are complementary and
interchangeable, at least in the present context. How general can such
framework be is, at the moment, unclear. The underlying algorithms - the EM-
algorithm and maximal-clique approach - are very different, so this agreement
should not be considered a rule. However, in a very tractable situation - an
$n$-dimensional Euclidean space $\mathbb{E}^{n}$\- this can be made more
precise, as follows: we have to show that sets of positives from the two
approaches - the EM-algorithm and the maximal clique - are identical; it is
well known that in $\mathbb{E}^{n}$ the k-means clustering algorithm is a form
of EM-algorithm (see [5]), with clusters (i.e. sets of positives) given by
$n$-balls, where

$B(x_{0};r)=\\{x\in\mathbb{E}^{n};|x-x_{0}|\leq r\\}$ (10)

is an $n$-ball around the point $x_{0}$, with the radius $r$; on the other
hand, set $d=2r$, and note that $B=B(x_{0};r)$ can be realized as the largest
subset of $\mathbb{E}^{n}$ such that

$B=\\{x\in\mathbb{E}^{n};|x-y|\leq d,\forall y\in B\\},$ (11)

which yields a clique-like object.

Consequently, we see that the considerable improvement in accuracy that we
have recorded is not a question of a superior, but complementary approach.
Namely, our queries - deliberately - consist of a single string, sometimes
even not a consensus query, and iterative process reaches a local optimum - a
set of positives with a fairly low F1-score and a blurred signal. And this is
a situation where maximal clique approach offers greatest gains, by filtering
the response and providing a better foundation for the next step in analysis.

####  3.2.6 Availability

http://compbioserv.math.hr/igloss/index.html?clique

##  4 Tables and Figures

###  4.1 BLAST

Figure 1: BLAST FVFGDSLSDA

Figure 2: BLAST FVFNDSLSDA

Figure 3: BLAST VFFGDSLSDN

###  4.2 JackHMMER

Figure 4: HMMER FVFGDSLSDA

Figure 5: HMMER FVFNDSLSDA

Figure 6: HMMER VFFGDSLSDN

###  4.3 IGLOSS

Figure 7: IGLOSS FVFGDSLSDA

Figure 8: IGLOSS FVFNDSLSDA

Figure 9: IGLOSS VFFGDSLSDN

###  4.4 Examples

(a) robustness BLAST

(b) robustness IGLOSS

Figure 10: F1-threshold level

###  4.5 Tables

| AT | OS | ST | SL | BV | ALL
---|---|---|---|---|---|---
CP | 118 | 155 | 123 | 108 | 86 | 590
| IGLOSS |
scale | 5.0 | 6.0 | 5.0 | 4.9 | 4.6 | 5.1
TP/P | 105/421 | 111/321 | 95/389 | 94/408 | 64/413 | 469/1952
PPV | 0.2494 | 0.3458 | 0.2442 | 0.2304 | 0.1550 | 0.2403
TPR | 0.8898 | 0.7161 | 0.7724 | 0.8704 | 0.7442 | 0.7949
F1 | 0.3896 | 0.4664 | 0.3711 | 0.3643 | 0.2565 | 0.3690
| IGLOSS + CLIQUE |
TP/P | 101/115 | 106/120 | 87/97 | 87/92 | 60/67 | 441/491
PPV | 0.8783 | 0.8833 | 0.8969 | 0.9457 | 0.8955 | 0.8982
TPR | 0.8559 | 0.6839 | 0.7073 | 0.8056 | 0.6977 | 0.7475
F1 | 0.8670 | 0.7709 | 0.7909 | 0.8700 | 0.7843 | 0.8159
Table 1: GDSL lipases IGLOSS vs IGLOSS+CLIQUE (FVFGDSLSDA) | AT | OS | ST | SL | BV | ALL
---|---|---|---|---|---|---
CP | 118 | 155 | 123 | 108 | 86 | 590
| IGLOSS |
scale | 6.6 | 7.6 | 6.6 | 6.5 | 6.2 | 6.7
TP/P | 80/180 | 118/176 | 92/174 | 93/174 | 63/156 | 446/860
PPV | 0.4444 | 0.6705 | 0.5287 | 0.5345 | 0.4038 | 0.5186
TPR | 0.6780 | 0.7613 | 0.7480 | 0.8611 | 0.7326 | 0.7559
F1 | 0.5369 | 0.7130 | 0.6195 | 0.6596 | 0.5207 | 0.6152
| IGLOSS + CLIQUE |
TP/P | 75/92 | 116/124 | 87/95 | 86/90 | 59/66 | 423/467
PPV | 0.8152 | 0.9355 | 0.9158 | 0.9556 | 0.8939 | 0.9058
TPR | 0.6356 | 0.7484 | 0.7073 | 0.7963 | 0.6860 | 0.7169
F1 | 0.7143 | 0.8315 | 0.7982 | 0.8687 | 0.7763 | 0.8004
Table 2: GDSL lipases IGLOSS vs IGLOSS+CLIQUE (FVFGDSLSDA) | AT | OS | ST | SL | BV | ALL
---|---|---|---|---|---|---
CP | 118 | 155 | 123 | 108 | 86 | 590
| BLAST |
ev | 226 | 226 | 226 | 226 | 226 | 226
TP/P | 26/98 | 32/51 | 26/ 109 | 32/108 | 35/87 | 151/453
PPV | 0.2653 | 0.6275 | 0.2385 | 0.2963 | 0.4023 | 0.3333
TPR | 0.2203 | 0.2065 | 0.2114 | 0.2963 | 0.4070 | 0.2559
F1 | 0.2407 | 0.3107 | 0.2241 | 0.2963 | 0.4046 | 0.2895
| BLAST + CLIQUE |
TP/P | 24/35 | 31/38 | 10/40 | 12/41 | 34/37 | 111/191
PPV | 0.6857 | 0.8158 | 0.2500 | 0.2927 | 0.9189 | 0.5812
TPR | 0.2034 | 0.2000 | 0.0813 | 0.1111 | 0.3953 | 0.1881
F1 | 0.3137 | 0.3212 | 0.1227 | 0.1611 | 0.5528 | 0.2843
Table 3: BLAST vs BLAST+CLIQUE (VFFGDSLSDN) | AT | OS | ST | SL | BV | ALL
---|---|---|---|---|---|---
CP | 118 | 155 | 123 | 108 | 86 | 590
| BLAST |
ev | 203 | 203 | 203 | 203 | 203 | 203
TP/P | 24/81 | 29/44 | 24/98 | 29/93 | 32/80 | 138/396
PPV | 0.2963 | 0.6591 | 0.2449 | 0.3118 | 0.4000 | 0.3485
TPR | 0.2034 | 0.1871 | 0.1951 | 0.2685 | 0.3721 | 0.2339
F1 | 0.2412 | 0.2915 | 0.2172 | 0.2886 | 0.3855 | 0.2799
| BLAST + CLIQUE |
TP/P | 1/30 | 29/35 | 10/40 | 12/39 | 30/33 | 82/177
PPV | 0.03333 | 0.8286 | 0.2500 | 0.3077 | 0.9091 | 0.4633
TPR | 0.00847 | 0.1871 | 0.0813 | 0.1111 | 0.3488 | 0.1390
F1 | 0.01351 | 0.3053 | 0.1227 | 0.1633 | 0.5042 | 0.2138
Table 4: BLAST vs BLAST+CLIQUE (VFFGDSLSDN)

## References

  * [1] Altschul, Stephen F and Madden, Thomas L and Schäffer, Alejandto A and Zhang, Jinghui and Zhang, Zheng and Miller, Webb and Lipman, David J  Gapped BLAST and PSI-BLAST:a new generation of protein database search programs.  Nucleic Acids Research. 1997;25(17):3389-402.
  * [2] Bateman A, Coin L, Durbin R, Finn RD, Hollich V, Griffiths-Jones S, et al.  The Pfam protein families database.  Nucleic Acids Research. 2004;32(DI):D138-41.
  * [3] Bron C, Kerbosch J.  Algorithm 457: finding all cliques of an undirected graph;  Communications of the ACM. 1973;16(9):575-77.
  * [4] Dayhoff M, Schwartz R.  A Model of Evolutionary Change in Proteins.  In Atlas of protein sequence and structure. 1978:345-52.
  * [5] Hastie, T et al,  Elements of Statistical Learning,  Springer, 2009
  * [6] Finn RD, Clements J, Arndt W, Miller BL, Wheeler TJ, Schreiber F, et al.  HMMER: web server 2015 update.  Nucleic Acids Research. 2015;43(W1):W30-8.
  * [7] Henikoff S, Henikoff JG.  Amino acid substitution matrices from protein blocks.  Proc Natl Acad Sci U S A. 1992;89(22):10915-9.
  * [8] Pirovano W. and Heringa J.  Protein Secondary Structure Prediction.  Methods in Molecular Biology. 2010;609.
  * [9] Rabar B, Zagorščak M, Ristov S, Rosenzweig M, Goldstein P.  IGLOSS: iterative gapless local similarity search.  Bioinformatics. 2019;35(18):3491-2.
  * [10] Sigrist CJA, de Castro E, Cerutti L, Cuche BA, Hulo N, Bridge A, et al.  New and continuing developments at PROSITE.  Nucleic Acids Research. 2013;41(D1):D344-7.
  * [11] Vujaklija I, Bielen A, Paradžik T, Bidin S, Goldstein P, Vujaklija D.  An effective approach for annotation of protein families with low sequence similarity and conserved motifs: Identifying GDSL hydrolases across the plant kingdom.  BMC Bioinformatics. 2016;17(1):1-17.
  * [12] Živa Ramšak, Špela Baebler, Rotter A, Korbar M, Mozetič I, Usadel B, et al.  GoMapMan: integration, consolidation and visualization of plant gene annotations within the MapMan ontology.  Nucleic Acids Research. 2014;42(D1):D1167-75.


Notes on a non-thermal fluctuation-dissipation relation in quantum Brownian
motion

Xinyi<EMAIL_ADDRESS>
NORDITA
Hannes Alfvéns väg 12, SE-106 91 Stockholm, Sweden

Abstract

We review how unitarity and stationarity in the Schwinger-Keldysh formalism
naturally lead to a (quantum) generalized fluctuation-dissipation relation
(gFDR) that works beyond thermal equilibrium. Non-Gaussian loop corrections
are also presented. Additionally, we illustrate the application of this gFDR
in various scenarios related to quantum Brownian motion and the generalized
Langevin equation.

###### Contents

  1. 1 Introduction
  2. 2 Effective SK theories
    1. 2.1 Open quantum systems
    2. 2.2 Stationary regime
    3. 2.3 Generating functional for stationary states
    4. 2.4 Quadratic effective action
    5. 2.5 Green’s function
    6. 2.6 Loop corrections
    7. 2.7 Thermalization condition
  3. 3 Quantum Brownian motion
    1. 3.1 Bath in thermal equilibrium
    2. 3.2 Multiple baths
    3. 3.3 Out-of-equilibrium bath
    4. 3.4 Anomalous diffusion
    5. 3.5 Anharmonic oscillator
  4. 4 Conclusion
  5. A Schwinger-Keldysh
    1. A.1 Density operator
    2. A.2 Time-evolution kernel
    3. A.3 Generating functional
      1. A.3.1 Properties
    4. A.4 Perturbation theory
    5. A.5 Keldysh basis
    6. A.6 Green’s functions
  6. B Fluctuation-Dissipation Theorem
  7. C Generalized Langevin equation
    1. C.1 Solution
    2. C.2 Stationary limit
  8. D Fourier transform

##  1 Introduction

The Schwinger-Keldysh (SK) formalism [1, 2] is a well-known path integral
framework to study quantum non-equilibrium many-body systems [3]. In the
classical limit, it reduces to the Martin-Siggia-Rose (MSR) path integral,
which is also equivalent to the stochastic Langevin description.

One important aspect of the SK approach is that symmetries become a powerful
tool to derive important relations in both quantum and classical regimes. For
example, the thermodynamic equilibrium can be formulated as a symmetry of the
SK action [4], and leads to the celebrated fluctuation-dissipation theorem
(FDT). Another example is the breaking of the time-reversal symmetry giving
rise to fluctuation theorems in non-equilibrium physics [5, 6]. Besides, the
path integral approach allows a systematic and perturbative way to deal with
non-linear theories beyond the scope of the standard Langevin approach.

These notes focus on the unitarity constraint on the SK two-point correlation
functions in the stationary limit, which can be viewed as a generalized
fluctuation-dissipation relation (gFDR) that works beyond thermal stationary
states and Markov approximation. While the classical limit of this gFDR is
certainly known in the study of aging glassy systems, where it is used to
define an effective temperature [7], its use in the broader classical
stochastic community does not seem widespread, in contrast to the use of the
non-Markovian generalization of the Langevin equation (GLE).

The objective of these notes is not to provide an extensive review of existing
literature222The present work was conducted mostly independently and with
limited awareness of the literature on glassy systems. It is part of my
journey to learn about classical and quantum non-equilibrium frameworks..
Instead, the aim is to present a self-contained derivation of the gFDR for
non-driven non-equilibrium systems and illustrate its application through
various examples.

The paper is structured as follows. In section 2, we review the SK framework
for open quantum systems with one effective degree of freedom. Then, we impose
the stationary condition at the generating functional level. Next we
perturbatively derive the gFDR in Fourier space. In section 3, the gFDR is
applied to GLE and quantum Brownian motion (QBM). Finally, we conclude in
section 4. In the appendix, we provide additional review material on the SK
formalism in section A, FDT in section B, and GLE in section C.

##  2 Effective SK theories

###  2.1 Open quantum systems

Consider a quantum system whose action is separable:

$S=S_{\text{subsystem}}+S_{\text{environment}}+S_{\text{interaction}}.$ (1)

By integrating out the environment degrees of freedom, we obtain an SK
effective theory for the subsystem, which is equivalent to the Feynman-Vernon
theory [8]. The SK generating functional, see A.3, can be rewritten in terms
of the reduced density matrix $\rho_{r}$ and the reduced time evolution kernel
$J_{r}$:

$\displaystyle Z(t,t_{0};\phi,\phi^{\prime}]$
$\displaystyle=\int_{-\infty}^{\infty}dx_{0}\int_{-\infty}^{\infty}dx_{0}^{\prime}\rho_{r}(x_{0},x_{0}^{\prime},t_{0})\int_{-\infty}^{\infty}dx\,J_{r}(x,x,t,x_{0},x_{0}^{\prime},t_{0};\phi,\phi^{\prime}]$
(2)

with

$\displaystyle
J_{r}(x,x^{\prime},t,x_{0},x_{0}^{\prime},t_{0};\phi,\phi^{\prime}]$
$\displaystyle=\int_{q(t_{0})=x_{0}}^{q(t)=x}\mathcal{D}q\int_{q^{\prime}(t_{0})=x_{0}^{\prime}}^{q^{\prime}(t)=x^{\prime}}\mathcal{D}q^{\prime}\,\exp{\left(\frac{i}{\hbar}A_{\text{CG}}[q,q^{\prime};\phi,\phi^{\prime}]\right)},$
(3)

where the subindex CG emphasizes a coarse-grained effective action, and
$\phi,\phi^{\prime}$ are the external sources. The curly/square parenthesis
notation refers to function/functional of the left/right parameters inside,
separated by a semi-colon. This action includes the action that affects only
the subsystem, and the influence phase $\Phi[q,q^{\prime}]$, which summarizes
all the contributions from the environment, its initial state and its
interaction with the subsystem:

$\displaystyle A_{\text{CG}}[q,q^{\prime};\phi,\phi^{\prime}]$
$\displaystyle=S_{\text{subsystem}}[q]-S_{\text{subsystem}}[q^{\prime}]+\Phi[q,q^{\prime}]$
$\displaystyle+\int_{t_{0}}^{t}ds\,(q(s)\phi(s)-q^{\prime}(s)\phi^{\prime}(s)).$
(4)

Here, it is also assumed that the initial density matrix factorizes as a
product of the reduced density matrix and that of the environment. It does not
matter, though, as we are interested only in the stationary regime.

###  2.2 Stationary regime

Let us define the _stationary regime_ as the one with the correlators being

  1. 1.

independent on the initial conditions and distribution
$\rho(x_{0},x_{0}^{\prime},t_{0})$,

  2. 2.

time-translation invariant.

###  2.3 Generating functional for stationary states

Consider a quadratic effective action
$A_{0}[q,q^{\prime};\phi,\phi^{\prime}]$. Its corresponding time evolution
kernel is a double Gaussian path integral. Therefore, it is solvable and
splits into the classical path contribution and fluctuations around the
classical path:

$\displaystyle J_{r}(x,x,t,x_{0},x_{0}^{\prime},t_{0};\phi,\phi^{\prime}]$
$\displaystyle=\int_{q(t_{0})=x_{0}}^{q(t)=x}\mathcal{D}q\int_{q^{\prime}(t_{0})=x_{0}^{\prime}}^{q^{\prime}(t)=x}\mathcal{D}q^{\prime}\,\exp{\left(\frac{i}{\hbar}A_{0}[q,q^{\prime};\phi,\phi^{\prime}]\right)}$
(5)
$\displaystyle=C(t,t_{0})\exp{\left(\frac{i}{\hbar}S_{\mathrm{cl}}(x,x_{0},x_{0}^{\prime};\phi,\phi^{\prime}]\right)},$
(6)

where the integration over quantum fluctuations hidden in $C(t,t_{0})$, does
not depend on the external sources.

The classical action is separable into a term independent of the initial
conditions and a term that is not:

$\displaystyle S_{\mathrm{cl}}(x,x_{0},x_{0};\phi,\phi^{\prime}]$
$\displaystyle=S_{\mathrm{cl}}(x,0,0;\phi,\phi^{\prime}]+\delta
S(x,x_{0},x_{0}^{\prime},\phi,\phi^{\prime}]$ (7) $\displaystyle\approx
S_{\mathrm{cl}}^{\infty}(x,0,0;\phi,\phi^{\prime}],\quad(t\longrightarrow\infty)$
(8)

where we imposed the condition of independency on the initial distribution, in
the stationary limit. The time-translational invariance condition allows us to
take the initial and final times to be minus and plus infinities. Finally, the
generating functional for stationary states is simplified to:

$\displaystyle Z_{0}^{\infty}[\phi,\phi^{\prime}]$ $\displaystyle\equiv
Z_{0}(\infty,-\infty;\phi,\phi^{\prime}]$ (9)
$\displaystyle=C(\infty,-\infty)\int_{-\infty}^{\infty}dx\,\exp{\left(\frac{i}{\hbar}S_{\mathrm{cl}}^{\infty}(x,0,0;\phi,\phi^{\prime}]\right)},$
(10)

since the integrations over the initial distribution is just one.

For quadratic theories perturbed by an external potential $V(q,q^{\prime})$,
see appendix A.4, we can do perturbation theory as usual:

$\displaystyle Z_{V}^{\infty}[\phi,\phi^{\prime}]$
$\displaystyle=\exp\left(-\frac{i}{\hbar}\int_{-\infty}^{\infty}dt\,V\left(\frac{\hbar}{i}\frac{\delta}{\delta\phi},\frac{\hbar}{i}\frac{\delta}{\delta\phi^{\prime}}\right)\right)Z_{0}^{\infty}[\phi,\phi^{\prime}].$
(11)

###  2.4 Quadratic effective action

In the stationary limit, the most general SK quadratic effective action (the
external currents are set to zero here for simplicity) in Keldysh basis is as
follows:

$\displaystyle A_{0}[q_{r},q_{a}]=$
$\displaystyle\int_{-\infty}^{\infty}dt^{\prime}\,q_{a}(t^{\prime})\left\\{(D_{ar}*q_{r})(t^{\prime})+i(D_{aa}*q_{a})(t^{\prime})\right\\},$
(12)

where the star operation is the Fourier convolution. Note that $q_{r}-q_{r}$
term is prohibited by unitarity, see appendix A. The kernel $D_{ar}$ contains
time derivative operators, and $D_{ar}$ is non-zero and satisfies

$D_{aa}(-t)=D_{aa}(t).$ (13)

The action above can be identified with the stochastic Martin-Siggia-Rose
action, where $q_{a}(t)$ plays the role of the auxiliary field (see also
appendix of A of [9]). Therefore, it is also equivalent to a generalized
stationary Langevin equation with a Gaussian noise $F(t)$:

$\displaystyle-(D_{ar}*q_{r})(t)$ $\displaystyle=F(t)$ (14)
$\displaystyle\frac{1}{2}\braket{\\{F(t),F(t^{\prime})\\}}$
$\displaystyle=D_{aa}(t-t^{\prime}).$ (15)

where the curly parenthesis represents the anticommutator.

###  2.5 Green’s function

Because the action (12) is time-translation invariant, we can Fourier
transform it:

$\displaystyle\tilde{A}_{0}[\tilde{q}_{r},\tilde{q}_{a}]=$
$\displaystyle\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega\,\tilde{q}_{a}(-\omega)\left\\{\tilde{D}_{ar}(\omega)\tilde{q}_{r}(\omega)+i\tilde{D}_{aa}(\omega)\tilde{q}_{a}(\omega)\right\\},$
(16)

from which two Feynman diagrams can be read off, see Fig. 1. The $q_{a}-q_{a}$
diagram is the interaction vertex, and the $q_{r}-q_{a}$ diagram is the
retarded Green’s function, which is:

$\tilde{G}^{R}(\omega)=-\frac{1}{\tilde{D}_{ar}(\omega)},$ (17)

since it is the solution of the homogeneous equation:

$-(D_{ar}*G^{R})(t)=\delta(t).$ (18)

Of course, the retarded Green’s function must satisfy the causality condition

$\displaystyle G^{R}(t<0)=0,$ (19)

which implies that the lower half plane in the frequency $\omega$-space must
be analytic.

The autocorrelator, or the symmetric Green’s function (100), is a simple
composite diagram (see Fig. 2) in this limit:

$\tilde{G}^{S}(\omega)=\left|\tilde{G}^{R}(\omega)\right|^{2}\tilde{D}_{aa}(\omega).$
(20)

For complex retarded Green functions, we use the following identity

$\displaystyle\mathrm{Im}\tilde{G}^{R}(\omega)=\left|\tilde{G}^{R}(\omega)\right|^{2}\mathrm{Im}\tilde{D}_{ar}(\omega),$
(21)

to obtain:

$\boxed{\frac{\tilde{G}^{S}(\omega)}{\mathrm{Im}\tilde{G}^{R}(\omega)}=\frac{\tilde{D}_{aa}(\omega)}{\mathrm{Im}\tilde{D}_{ar}(\omega)}}.$
(22)

Given that $\tilde{G}^{S}(\omega)$ is even in $\omega$ space, the inverse
Fourier transform can be written as:

$\displaystyle G^{S}(t)$
$\displaystyle=\frac{1}{\pi}\int_{0}^{\infty}d\omega\,\tilde{G}^{S}(\omega)\cos(\omega
t)$ (23)
$\displaystyle=\frac{1}{\pi}\int_{0}^{\infty}d\omega\frac{\tilde{D}_{aa}(\omega)}{\mathrm{Im}\tilde{D}_{ar}(\omega)}\mathrm{Im}\tilde{G}^{R}(\omega-i0^{+})\cos(\omega
t),\quad(t\longrightarrow\infty)$ (24)

where the latter expression (we added the prescription for the poles of the
retarded Green’s function) is valid only if the stationarity conditions in
section 2.2 are fulfilled.

We will see some examples, but before that, let us generalize the above
relation for non-quadratic theories.

Figure 1: Feynman diagram for the retarded Green’s function $G_{ra}$ on the
left and the noise kernel $D_{aa}$ on the right. The solid line represents
$q_{r}$ and the dashed line is for $q_{a}$.  Figure 2: Feynman diagram for the
symmetric Green’s function $G_{rr}$.

###  2.6 Loop corrections

Small non-quadratic terms added to the action can be treated perturbatively
according to (11). For the retarded Green’s function, also known as the
propagator, the loop corrections (aka self-energy) which we call
$\delta\tilde{D}_{ar}(\omega)$, form a geometric sum as shown in Fig. 3:

$\displaystyle\tilde{G}_{r}^{R}(\omega)$
$\displaystyle=\tilde{G}^{R}(\omega)(1+\delta\tilde{D}_{ar}(\omega)\,\tilde{G}^{R}(\omega)+\delta\tilde{D}_{ar}(\omega)^{2}\,\tilde{G}^{R}(\omega)^{2}+\ldots)$
(25)
$\displaystyle=\frac{\tilde{G}^{R}(\omega)}{1-\delta\tilde{D}_{ar}(\omega)\,\tilde{G}^{R}(\omega)}.$
(26)

The equation above is known as the Dyson equation. In terms of the kernel
$\tilde{D}_{ar}(\omega)$, which is the negative inverse of the propagator, the
above expression simplifies to:

$\displaystyle\tilde{D}_{ar,r}(\omega)$
$\displaystyle=\tilde{D}_{ar}(\omega)+\delta\tilde{D}_{ar}(\omega).$ (27)

Note that we use an additional subindex $r$ to stand for
renormalized/corrected values, and it should not be confused with the $r$ from
the Keldysh variables.

The correction to the other kernel, shown as the black circle in Fig. 4, is
simply additive:

$\tilde{D}_{aa,r}(\omega)=\tilde{D}_{aa}(\omega)+\delta\tilde{D}_{aa}(\omega).$
(28)

Now, new interaction vertices connecting $q_{r}$ with $q_{r}$ can emerge, and
we call it $\delta Q(\omega)$ and represent it as the diamond in Fig. 4. This
type of interaction breaks the tree level composition rule (22). Nevertheless,
we can still write down a generic loop-corrected expression for the symmetric
Green’s function in a closed form, which corresponds to Fig. 5:

$\displaystyle\tilde{G}_{r}^{S}(\omega)=\tilde{D}_{aa,r}(\omega)\left|\tilde{G}_{r}^{R}(\omega)\right|^{2}(1+\delta
Q(\omega)\tilde{D}_{aa,r}(\omega)\left|\tilde{G}_{r}^{R}(\omega)\right|^{2}).$
(29)

Replacing the modulus squared by the imaginary part:

$\displaystyle\mathrm{Im}\tilde{G}^{R}_{r}(\omega)=\left|\tilde{G}_{r}^{R}(\omega)\right|^{2}\mathrm{Im}\tilde{D}_{ar,r}(\omega),$
(30)

we can rewrite the above expression as:

$\displaystyle\frac{\tilde{G}_{r}^{S}(\omega)}{\mathrm{Im}\tilde{G}^{R}_{r}(\omega)}$
$\displaystyle=\frac{\tilde{D}_{aa,r}(\omega)}{\mathrm{Im}\tilde{D}_{ar,r}(\omega)}\left(1+\delta
Q(\omega)\frac{\tilde{D}_{aa,r}(\omega)}{\mathrm{Im}\tilde{D}_{ar,r}(\omega)}\mathrm{Im}\tilde{G}^{R}_{r}(\omega)\right).$
(31)

This expression is useful when we Taylor-expand the right-hand-side (RHS) in
terms of the unperturbed quantities and the loop contributions of different
orders. For the leading order correction, we get

$\boxed{\frac{\tilde{G}_{r}^{S}(\omega)}{\mathrm{Im}\tilde{G}^{R}_{r}(\omega)}\approx\frac{\tilde{G}^{S}(\omega)}{\mathrm{Im}\tilde{G}^{R}(\omega)}\left(1+\frac{\delta\tilde{D}_{aa}(\omega)}{\tilde{D}_{aa}(\omega)}-\frac{\mathrm{Im}\delta\tilde{D}_{ar}(\omega)}{\mathrm{Im}\tilde{D}_{ar}(\omega)}+\delta
Q(\omega)\,\tilde{G}^{S}(\omega)\right)},$ (32)

where the three corrections are shown in Fig. 6.

Figure 3: Renormalized retarded Green’s function. The self-energy
$\delta\tilde{D}_{ar}(\omega)$ is the white square.  Figure 4: Corrections to
the interaction vertices. The solid line represents $q_{r}$ and the dashed
line $q_{a}$. The black circle is $\delta\tilde{D}_{aa}(\omega)$ and the
diamond is $\delta Q(\omega)$. Figure 5: The renormalized symmetric Green’s
function is composed of the renormalized retarded Green’s function, see Fig.
3, and the vertex corrections in Fig. 4. Figure 6: The 1-loop contributions
(the white geometric shapes are symbolic only) to the renormalized symmetric
Green’s function.

###  2.7 Thermalization condition

If thermalization happens, both (22) and (31) must satisfy FDT (116):

$\displaystyle\frac{\tilde{G}^{S}(\omega)}{\mathrm{Im}\tilde{G}^{R}(\omega)}$
$\displaystyle=\hbar\coth\frac{\hbar\beta\omega}{2},$ (33)
$\displaystyle\frac{\tilde{G}_{r}^{S}(\omega)}{\mathrm{Im}\tilde{G}^{R}_{r}(\omega)}$
$\displaystyle=\hbar\coth\frac{\hbar\beta_{r}\omega}{2}.$ (34)

Therefore

$\displaystyle\coth\frac{\hbar\beta_{r}\omega}{2}$
$\displaystyle\approx\coth\frac{\hbar\beta\omega}{2}\left(1+\Delta(\omega)\right)$
(35)

where we use $\Delta(\omega)$ as a short-hand notation for the loop
contributions to the ratio of the Green’s functions, for example, at the
leading order, see (32), it is:

$\Delta(\omega)=\frac{\delta\tilde{D}_{aa}(\omega)}{\tilde{D}_{aa}(\omega)}-\frac{\mathrm{Im}\delta\tilde{D}_{ar}(\omega)}{\mathrm{Im}\tilde{D}_{ar}(\omega)}+\delta
Q(\omega)\,\tilde{G}^{S}(\omega).$ (36)

Note that $\Delta(\omega)$ can depend on the temperature, but its dependence
is not made explicit here.

On the other hand, FDT is a non-perturbative result. Therefore, we must
systematically expand the hyperbolic cotangent in terms of the corrections to
the inverse temperature, namely $\beta_{r}=\beta+\delta\beta$, as
$\delta\beta$ shall absorb the loop contributions in the thermal case. Then,
the expansion is to be compared order by order with the one from (35). Let us
show how it works for the leading order. First:

$\displaystyle\coth\frac{\hbar\beta_{r}\omega}{2}$
$\displaystyle=\coth\frac{\hbar(\beta+\delta\beta)\omega}{2}$ (37)
$\displaystyle\approx\coth\frac{\hbar\beta\omega}{2}\left(1-\hbar\delta\beta\omega\operatorname{csch}\hbar\beta\omega\right).$
(38)

Then, comparing the RHS of (35) with (38), we see that thermalization
condition at this order implies:

$\delta\beta\overset{!}{=}-\frac{\Delta(\omega)}{\hbar\omega\operatorname{csch}\hbar\beta\omega}.$
(39)

Since $\delta\beta$ must be a constant, there are three possible scenarios for
thermalization:

$\displaystyle\Delta(\omega)\begin{cases}\propto\omega\operatorname{csch}\hbar\beta\omega&\text{Quantum
thermalization}\\\ =0&\text{Quantum thermalization: }\beta_{r}=\beta\\\
=\text{constant}&\text{Classical/high-T thermalization}\\\ \end{cases}$ (40)

Again, this a perturbative result, hence the conclusion for (quantum or
classical) thermalization is only valid up to the level studied. This means
higher order terms could potentially drive the system away from
thermalization. On the other hand, the negative statement about thermalization
is sufficient with only perturbative information.

##  3 Quantum Brownian motion

In this section, we focus on one of the simplest integrable SK effective
theories, that is, a quantum Brownian motion. This is a quantum harmonic
oscillator coupled linearly to a harmonic bath. The classical action is the
sum of the following actions:

$\displaystyle S[x]=$
$\displaystyle\int^{t}_{0}ds\,\frac{1}{2}M\left(\dot{x}(s)^{2}-\Omega_{R}^{2}x(s)^{2}\right)$
(41) $\displaystyle S_{\mathrm{bath}}[\\{x_{n}\\}]=$
$\displaystyle\int^{t}_{0}ds\sum_{n}\frac{1}{2}m_{n}\left(\dot{x}_{n}(s)^{2}-\omega_{n}^{2}x_{n}(s)^{2}\right)$
(42) $\displaystyle S_{\mathrm{int}}[x,\\{x_{n}\\}]=$
$\displaystyle-\int^{t}_{0}ds\sum_{n}c_{n}x_{n}(s)x(s)$ (43)

where

$\Omega_{R}^{2}=\Omega^{2}+\delta\Omega^{2},\quad\delta\Omega^{2}=\sum_{n}\frac{c_{n}^{2}}{2m_{n}\omega_{n}^{2}}.$
(44)

In the continuum limit333The finite version is also integrable, but
thermalization happens only in the continuum limit., the bath is fully
characterized by a spectral density function $I(\omega)$, that is related to
the discrete frequencies as follow:

$I(\omega)=\sum_{n}\frac{c_{n}^{2}}{2m_{n}\omega_{n}}\delta\left(\omega-\omega_{n}\right).$
(45)

Traditionally, power law distributions $I(\omega)\propto\omega^{\alpha}$ were
studied, where $\alpha=1$ is called the Ohmic bath, and when $\alpha$ is
smaller or bigger than 1, sub-Ohmic or supra-Ohmic baths.

The Feynman-Vernon/SK approach gives rise to an effective quadratic action
(12) with the kernels [10, 11]:

$\displaystyle D_{ar}(t)$
$\displaystyle=-M\delta(t)(\partial^{2}_{t}+\Omega_{R}^{2})-\mu(t)$ (46)
$\displaystyle D_{aa}(t)$ $\displaystyle=\nu(t)$ (47)

where $\mu$ and $\nu$ are known as the dissipation and the noise kernel,
respectively. In this case of a linearly coupled bath, the dissipation kernel
in Fourier space is fully determined by the spectral density of the bath, [12,
13]444Note that [13] absorbs the $\pi$ in the definition of the spectral
density.:

$\tilde{\mu}(\omega)=i\pi I(\omega).$ (48)

This problem is solvable for all time [14], and it is equivalent to the
stationary generalized Langevin equation, see appendix C. The Fourier
transform of the retarded Green’s function is:

$\tilde{G}^{R}(\omega)=\frac{1/M}{-\omega^{2}+\Omega_{R}^{2}+\tilde{\mu}(\omega)/M},$
(49)

and the symmetric Green’s function in the stationary limit is obtained from
(22):

$\frac{\tilde{G}^{S}(\omega)}{\mathrm{Im}\tilde{G}^{R}(\omega)}=\frac{\tilde{\nu}(\omega)}{\mathrm{Im}\tilde{\mu}(\omega)}.$
(50)

###  3.1 Bath in thermal equilibrium

If the bath is in initial thermal equilibrium, the noise kernel is:

$\displaystyle\nu(t)$
$\displaystyle=\int_{0}^{\infty}d\omega\,\hbar\coth\left(\frac{\hbar\beta\omega}{2}\right)I(\omega)\cos\omega
t,$ (51)

which in Fourier space is just the FDT for the bath (sometimes referred as
2FDT in the classical literature):

$\frac{\tilde{\nu}(\omega)}{\mathrm{Im}\tilde{\mu}(\omega)}=\hbar\coth\left(\frac{\hbar\beta\omega}{2}\right)\xrightarrow{\hbar\rightarrow
0}\frac{2}{\beta\omega}.$ (52)

If the bath were to be modeled by a scalar field, $\nu$ and $\mu$ would
correspond precisely to the symmetric and retarded Green’s functions of the
field, see [9].

Now, the gFDR (50) implies:

$\frac{\tilde{G}^{S}(\omega)}{\mathrm{Im}\tilde{G}^{R}(\omega)}=\hbar\coth\left(\frac{\hbar\beta\omega}{2}\right)\xrightarrow{\hbar\rightarrow
0}\frac{2}{\beta\omega}.$ (53)

In other words, when the system eventually thermalizes, the final temperature
is the same as the one of the initial bath. This case is also studied and
emphasized in [9]. The FDT above is sometimes called 1FDT in the classical
literature, and it will differ from 2FDT in non-equilibrium scenarios.

###  3.2 Multiple baths

The sum of Gaussian random variables are Gaussian, hence, when there are many
harmonic baths, the gFDR (50)is easily generalized to:

$\frac{\tilde{G}^{S}(\omega)}{\mathrm{Im}\tilde{G}^{R}(\omega)}=\frac{\sum_{i}\tilde{\nu}_{i}(\omega)}{\sum_{i}\mathrm{Im}\tilde{\mu}_{i}(\omega)}.$
(54)

Even when the baths are in thermal equilibrium (i.e. satisfying (52)), it is
clear from the above expression that, in general, the system does not
thermalize. It was also shown in [5] that the multi-bath setting breaks the
thermal symmetry. Only in the classical limit, when all the baths have the
same spectral density, then, the particle thermalizes at an average
temperature of the baths:

$\frac{\tilde{G}^{S}(\omega)}{\mathrm{Im}\tilde{G}^{R}(\omega)}=\frac{2k_{B}\bar{T}}{\omega}.$
(55)

This case has been numerically studied for realistic thermal baths in [15].

The above expression can be generalized to effective temperatures for
classical out-of-equilibrium baths, see [16].

###  3.3 Out-of-equilibrium bath

The generalized Langevin equation is used to model classical glassy systems
with a slow relaxation dynamics [17]. The gFDR (50) naturally appears in this
context, where an effective temperature [7] defined for the bath encodes the
non-equilibrium properties.

In this section, we apply the gFDR to a relatively recent example of a
semiclassical time glass model in [18], that was solved using fractional
calculus. The microscopic model is characterized by the following noise kernel
and the spectral density (hence the dissipation kernel through (48))555The
paper [18] defines their spectral density $J(\omega)$ that differs a factor
$\pi$ from our definition, i.e. $J(\omega)=\pi I(\omega)$.:

$\displaystyle\nu(t)$
$\displaystyle=2\int_{0}^{\infty}d\omega\frac{(t_{s}\omega)^{1-s}}{\beta\omega}I(\omega)\cos(\omega
t)$ (56) $\displaystyle I(\omega)$
$\displaystyle=\frac{\eta}{\pi}\sin\left(\frac{\pi}{2}s\right)\omega^{s},$
(57)

where $0<s<1$. The integration gives:

$\displaystyle\nu(t)$
$\displaystyle=\frac{2}{\pi}\frac{t_{s}^{1-s}}{\beta}\eta\sin\left(\frac{\pi}{2}s\right)\int_{0}^{\infty}d\omega\cos(\omega
t)$ (58)
$\displaystyle=\frac{2t_{s}^{1-s}}{\beta}\eta\sin\left(\frac{\pi}{2}s\right)\delta(t).$
(59)

The prefactor before the Dirac delta function is the Fourier transform of the
noise kernel, therefore, together with (48), the relation (50) reduces to:

$\displaystyle\frac{\tilde{G}^{S}(\omega)}{\mathrm{Im}\tilde{G}^{R}(\omega)}=\frac{2t_{s}^{1-s}}{\beta\omega^{s}}.$
(60)

The effective frequency-dependent inverse temperature is then:

$\beta_{\text{eff}}(\omega)=(\omega t_{s})^{s-1}\beta.$ (61)

When $s=1$, which corresponds to an Ohmic bath, we recover the classical FDT
(118). According to [18], only when $s<1$, i.e. in the sub-Ohmic regime, the
non-equilibrium time glassy behavior appears.

###  3.4 Anomalous diffusion

The diffusion of a classical Brownian particle is characterized by the late-
time behavior of the mean square displacement (MSD):

$MSD(t)=\left\langle[q(t)-q(0)]^{2}\right\rangle.$ (62)

Ford and O’Connell in [19] related MSD to FDT and studied anomalous diffusion
for non-Ohmic baths. Their relation for the time derivative of MSD can be
readily generalized to:

$MSD^{\prime}(t)=\int_{0}^{\infty}d\omega\,\tilde{G}^{S}(\omega)\omega\sin(\omega
t)$ (63)

where the symmetric Green’s function can be determined by the gFDR (50).

Let us apply this method to the non-equilibrium time glass model in the
subsection above, and we will reproduce the asymptotic results in the
supplemental material of [18]. Recall (60) and the retarded Green’s function
(49) (with free particle as done in [18]):

$\displaystyle\tilde{G}_{S}(\omega)$
$\displaystyle=\frac{2t_{s}^{1-s}}{\beta\omega^{s}}\frac{\pi
I(\omega)M^{-2}}{\omega^{4}+\pi^{2}I(\omega)^{2}M^{-2}}$ (64)
$\displaystyle\propto\frac{1}{\omega^{4}+A^{2}\omega^{2s}}$ (65)
$\displaystyle\propto\frac{1}{\omega^{2s}(\omega^{4-2s}+A)}$ (66)
$\displaystyle\propto\frac{1}{\omega^{2s}},\quad(\omega\approx 0)$ (67)

where $A=M^{-1}\eta\sin\left(\frac{\pi}{2}s\right)$, and we took the small
frequency limit in order to obtain the late-time behavior for the time
derivative of the MSD (63), which is

$\displaystyle MSD^{\prime}(t)$
$\displaystyle\propto\int_{0}^{\infty}d\omega\,\frac{\sin(\omega
t)}{\omega^{2s-1}}\propto t^{2s-2}$ (68)

hence, for MSD, we recover the following asymptotic results:

$\displaystyle MSD(t)$ $\displaystyle\propto\log(t),\quad s=1/2,$ (69)
$\displaystyle MSD(t)$ $\displaystyle\propto t^{2s-1},\quad s\neq 1/2.$ (70)

###  3.5 Anharmonic oscillator

Now, let us show an example for the loop-corrected gFDR (32).

Instead of the harmonic oscillator, let us consider a quartic anharmonic
oscillator subject to the following SK potential:

$\displaystyle V(q,q^{\prime})=\frac{\lambda}{4!}(q^{4}-q^{\prime 4})$ (71)

which in Keldysh basis becomes

$\displaystyle
V(q_{r},q_{a})=\frac{\lambda}{4!}(q_{r}q_{a}^{3}-3q_{r}^{3}q_{a}),$ (72)

where the two contributions are diagrammatically represented in Fig. 7. With
these diagrams, we can build the leading order loop corrections shown in Fig.
8, which correspond to the $\delta\tilde{D}_{ar}$, $\delta\tilde{D}_{aa}$ and
$\delta Q$ term in (31), respectively. However, the loop integral with the
retarded Green’s function vanishes because of causality. Furthermore, the loop
integral with the symmetric Green’s function is real. That means, at leading
order in $\lambda$, the qFDT is _not_ corrected. Hence, at the leading order
perturbation, the anharmonic oscillator thermalizes with the same temperature
as the harmonic oscillator, as discussed in subsection 2.7. The loop
computation was done explicitly by Hsiang et al [20], and that is precisely
their result. Furthermore, the same authors have a non-perturbative argument
for thermalization discussed in [21]. Then, it would be interesting to compute
the next-to-leading order corrections, where the Feynman diagrams are shown in
Fig. 9, because these do not seem to vanish and could potentially violate the
thermalization condition.

Figure 7: The quartic interaction vertices Figure 8: Leading order
corrections. However, the diagrams with retarded Green’s function loop (the
two right ones) are vanishing because of causality. Their corresponding
complex conjugate ones are not shown. Figure 9: Next-to-leading order
corrections. Their corresponding complex conjugate ones are not shown.

##  4 Conclusion

In summary, we presented the generalized fluctuation-dissipation relation
(gFDR) (22) for the quadratic Schwinger-Keldysh (SK) path integral and its
loop generalization (31). The gFDR generalizes the fluctuation-dissipation
theorem (FDT) to non-thermal baths, providing valuable insight into FDT itself
and its breakdown in various scenarios. We were also able to reproduce the
asymptotic behavior of the correlation function and the anomalous diffusion
for the time glass model in [18]. Besides, the loop gFDR (31) can be used to
determine thermalization if the quadratic system is perturbed.

Beyond the presented cases, equations (22) and (31) apply to other effective
SK theories in the stationary limit, such as the SK effective field theory for
diffusion [22, 23]. This theory, however, is thermal by construction, leading
to the reduction of the gFDR to FDT for all loop orders, as one can check
[24]. Another example of an SK effective theory where these equations apply is
a non-linear SK diffusion model without an external driven force, which was
proposed to model a quantum time crystal state in [25]. Unitarity and
stationarity can also be used to generalize the non-linear quantum FDTs
presented in [26] to non-linear gFDRs. Finally, for driven systems, there is
some work done towards a FDR for non-equilibrium steady states in [27].

### Acknowledgements

We would like to thank K. Zarembo, S. Krishnamurthy and E. Aurell for
discussion and reviewing the manuscript. We would also like to thank L.
Cugliandolo for pointing out some relevant existing work. This work was
supported by the Knut and Alice Wallenberg Foundation.

##  Appendix A Schwinger-Keldysh

In this section, we briefly review the Schwinger-Keldysh formalism in the
coordinate respresentation. The reader can learn more from [3, 28, 29].

###  A.1 Density operator

Consider a time-evolving quantum statistical system characterized by the
density operator:

$\displaystyle\rho(t)$
$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}p_{n}\ket{\psi_{n}(t)}\bra{\psi_{n}(t)},$
(73)

such that the statistical weights $p_{n}$ sum to one, and the evolution of the
quantum states $\ket{\psi_{n}(t)}$ is governed by the Schroedinger equation:

$i\hbar\frac{d}{dt}\ket{\psi_{n}(t)}=H\ket{\psi_{n}(t)},$ (74)

whose solution is

$\ket{\psi_{n}(t)}=U(t-t_{0})\ket{\psi_{n}(t_{0})}$ (75)

with the evolution operator being

$U(t)=e^{-itH/\hbar}.$ (76)

Therefore the evolution of the density operator is:

$\rho(t)=U\left(t-t_{0}\right)\rho(t_{0})U^{\dagger}\left(t-t_{0}\right).$
(77)

###  A.2 Time-evolution kernel

The density operator in the coordinate basis is:

$\displaystyle\rho(x,x^{\prime},t)$
$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}p_{n}\,\psi_{n}(x,t)\psi_{n}(x^{\prime},t)^{*}$
(78)
$\displaystyle=\int_{-\infty}^{\infty}dx_{0}\,\int_{-\infty}^{\infty}dx_{0}^{\prime}\,J(x,x^{\prime},t,x_{0},x_{0}^{\prime},t_{0})\rho(x_{0},x_{0}^{\prime},t_{0}),$
(79)

where the time-evolution kernel is a double copy of the quantum mechanical
time-evolution kernel:

$\displaystyle J(x,x^{\prime},t,x_{0},x_{0}^{\prime},t_{0})$
$\displaystyle\equiv\langle x|U(t-t_{0})|x_{0}\rangle\langle
x^{\prime}|U(t-t_{0})|x^{\prime}_{0}\rangle^{*}$ (80)
$\displaystyle=\int_{q(t_{0})=x_{0}}^{q(t)=x}\mathcal{D}q\int_{q^{\prime}(t_{0})=x_{0}^{\prime}}^{q^{\prime}(t)=x^{\prime}}\mathcal{D}q^{\prime}\,\exp{\left(\frac{i}{\hbar}A[q,q^{\prime}]\right)},$
(81)

where the action is:

$A[q,q^{\prime}]=S[q]-S[q^{\prime}].$ (82)

Notice that the time-evolution kernel is time-translation invariance, because
of the unitary evolution.

###  A.3 Generating functional

We can probe the system by minimally coupling it to an external source
$\phi(t)$, that means:

$\displaystyle J(x,x,t,x_{0},x_{0}^{\prime},t_{0};\phi,\phi^{\prime}]$
$\displaystyle=\int_{q(t_{0})=x_{0}}^{q(t)=x}\mathcal{D}q\int_{q^{\prime}(t_{0})=x_{0}^{\prime}}^{q^{\prime}(t)=x^{\prime}}\mathcal{D}q^{\prime}\,\exp{\left(\frac{i}{\hbar}A[q,q^{\prime};\phi,\phi^{\prime}]\right)}$
(83) $\displaystyle A[q,q^{\prime};\phi,\phi^{\prime}]$
$\displaystyle=S[q]-S[q^{\prime}]+\int_{t_{0}}^{t}ds\,(q(s)\phi(s)-q^{\prime}(s)\phi^{\prime}(s)).$
(84)

Then the partition function of this system is the generating functional for
the correlation functions of the original system:

$\displaystyle Z[\phi,\phi^{\prime}]$
$\displaystyle=\,\int_{-\infty}^{\infty}dx_{0}\,\int_{-\infty}^{\infty}dx_{0}^{\prime}\,\rho(x_{0},x_{0}^{\prime},t_{0})\int_{-\infty}^{\infty}dx\,J(x,x,t,x_{0},x_{0}^{\prime},t_{0};\phi,\phi^{\prime}].$
(85)

In particular, its logarithm

$W[\phi,\phi^{\prime}]\equiv\ln Z[\phi,\phi^{\prime}]$ (86)

generates connected correlators.

####  A.3.1 Properties

The unitary evolution of the density matrix implies the following properties
for the SK generating functional:

  1. 1.

Normalization condition:

$Z[\phi,\phi]=1\longrightarrow W[\phi,\phi^{\prime}]=0$ (87)

  2. 2.

Reflection symmetry:

$Z[\phi^{\prime},\phi]=Z[\phi,\phi^{\prime}]^{*}\longrightarrow
W[\phi^{\prime},\phi]=W[\phi,\phi^{\prime}]^{*}$ (88)

###  A.4 Perturbation theory

For a quadratic action $A_{0}[q,q^{\prime};\phi,\phi^{\prime}]$ perturbed with
a non-quadratic potential $V(q,q^{\prime})$:

$\displaystyle
A_{V}[q,q^{\prime};\phi,\phi^{\prime}]=A_{0}[q,q^{\prime};\phi,\phi^{\prime}]-\int_{t_{0}}^{t}ds\,V(q(s),q^{\prime}(s)),$
(89)

the corresponding generating functional of the resulting action
$A_{V}[q,q^{\prime};\phi,\phi^{\prime}]$ can be expressed in terms of the
unperturbed one:

$\displaystyle Z_{V}[\phi,\phi^{\prime}]$
$\displaystyle=\exp\left(-\frac{i}{\hbar}\int_{t_{0}}^{t}ds\,V\left(\frac{\hbar}{i}\frac{\delta}{\delta\phi},\frac{\hbar}{i}\frac{\delta}{\delta\phi^{\prime}}\right)\right)Z_{0}[\phi,\phi^{\prime}]$
(90)
$\displaystyle\approx\left(1-\frac{i}{\hbar}\int_{t_{0}}^{t}ds\,V\left(\frac{\hbar}{i}\frac{\delta}{\delta\phi},\frac{\hbar}{i}\frac{\delta}{\delta\phi^{\prime}}\right)+\ldots\right)Z_{0}[\phi,\phi^{\prime}],$
(91)

where in the last step, we expanded the exponential.

###  A.5 Keldysh basis

In practice, it is more convenient to work in the so-called Keldysh basis,
i.e.:

$\displaystyle q_{r}$ $\displaystyle=\frac{1}{2}(q+q^{\prime}),\quad
q_{a}=q-q^{\prime},$ (92)

that also applies to the external currents. Then, the action (84) becomes:

$\displaystyle
A[q_{r},q_{a};\phi_{r},\phi_{a}]=S[q_{r},q_{a}]+\int_{t_{0}}^{t}ds\,(q_{r}(s)\phi_{a}(s)+q_{a}(s)\phi_{r}(s)).$
(93)

###  A.6 Green’s functions

The two-point (path-ordered) correlation functions in Keldysh basis are
defined below:

$\displaystyle G_{rr}(t_{1},t_{2})$
$\displaystyle\equiv\braket{\mathcal{P}q_{r}(t_{1})q_{r}(t_{2})}=\left(\frac{\hbar}{i}\right)^{2}\left.\frac{\delta^{2}W[\phi_{r},\phi_{a}]}{\delta\phi_{a}(t_{1})\phi_{a}(t_{2})}\right|_{\phi_{r,a}=0}$
(94) $\displaystyle\frac{\hbar}{i}G_{ra}(t_{1},t_{2})$
$\displaystyle\equiv\braket{\mathcal{P}q_{r}(t_{1})q_{a}(t_{2})}=\left(\frac{\hbar}{i}\right)^{2}\left.\frac{\delta^{2}W[\phi_{r},\phi_{a}]}{\delta\phi_{a}(t_{1})\phi_{r}(t_{2})}\right|_{\phi_{r,a}=0}$
(95) $\displaystyle\frac{\hbar}{i}G_{ar}(t_{1},t_{2})$
$\displaystyle\equiv\braket{\mathcal{P}q_{a}(t_{1})q_{r}(t_{2})}=\left(\frac{\hbar}{i}\right)^{2}\left.\frac{\delta^{2}W[\phi_{r},\phi_{a}]}{\delta\phi_{r}(t_{1})\phi_{a}(t_{2})}\right|_{\phi_{r,a}=0}$
(96) $\displaystyle\left(\frac{\hbar}{i}\right)^{2}G_{aa}(t_{1},t_{2})$
$\displaystyle\equiv\braket{\mathcal{P}q_{a}(t_{1})q_{a}(t_{2})}=\left(\frac{\hbar}{i}\right)^{2}\left.\frac{\delta^{2}W[\phi_{r},\phi_{a}]}{\delta\phi_{r}(t_{1})\phi_{r}(t_{2})}\right|_{\phi_{r,a}=0}=0.$
(97)

The latter one vanishes due to the normalization condition (87). In fact, (87)
implies that all

$G_{a\ldots a}=0.$ (98)

This identity is sometimes called _Schwinger-Keldysh collapse rule_ in the
literature.

The reflection symmetry (88) implies:

$G_{ra}(t_{1},t_{2})=G_{ar}(t_{2},t_{1}),$ (99)

i.e. it relates the retarded and the advanced Green’s functions, since the
non-vanishing two-point correlators are identified with the symmetric,
retarded and advanced Green’s functions:

$\displaystyle G_{rr}(t_{1},t_{2})$
$\displaystyle=G^{S}(t_{1},t_{2})\equiv\frac{1}{2}\langle\\{q(t_{1}),q(t_{2})\\}\rangle$
(100) $\displaystyle G_{ra}(t_{1},t_{2})$
$\displaystyle=G^{R}(t_{1},t_{2})\equiv\frac{i}{\hbar}\,\theta(t_{1}-t_{2})\Delta(t_{1},t_{2})$
(101) $\displaystyle G_{ar}(t_{1},t_{2})$
$\displaystyle=G^{A}(t_{1},t_{2})\equiv-\frac{i}{\hbar}\,\theta(t_{2}-t_{1})\Delta(t_{1},t_{2}),$
(102)

where

$\Delta(t_{1},t_{2})=\langle[q(t_{1}),q(t_{2})]\rangle.$ (103)

Hence, for a quadratic theory, the influence phase is simply:

$\displaystyle W_{0}[\phi_{r},\phi_{a}]$
$\displaystyle=\frac{i}{\hbar}\frac{1}{2}\int_{t_{0}}^{t}dt_{1}\int_{t_{0}}^{t}dt_{2}\,\Phi(t_{1})^{T}\mathbf{G}(t_{1},t_{2})\Phi(t_{2})$
(104)

where

$\displaystyle\Phi$ $\displaystyle=\left(\begin{array}[]{l}\phi_{r}\\\
\phi_{a}\end{array}\right),\quad\mathbf{G}=\left(\begin{array}[]{cc}0&G^{A}\\\
G^{R}&\frac{i}{\hbar}G^{S}\end{array}\right).$ (109)

##  Appendix B Fluctuation-Dissipation Theorem

When the density matrix is thermal, i.e. a Gibbs state:

$\displaystyle\rho_{\beta}=\frac{1}{Z_{\beta}}e^{-\beta H},\quad
Z_{\beta}=\mathrm{Tr}\,\left(e^{-\beta H}\right),$ (110)

then the generating functional satisfies the so-called Kubo-Martin-Schwinger
(KMS) condition [30, 31]. In the case of thermal two-point correlators, the
KMS condition translates into the well-known quantum fluctuation-dissipation
theorem (FDT). Let us show how to derive it. First,

$\displaystyle\langle O_{i}(t_{1})O_{j}(t_{2})\rangle$
$\displaystyle=\frac{1}{Z_{\beta}}\mathrm{Tr}\,[e^{-\beta
H}U^{\dagger}(t_{1}-t_{0})O_{i}U(t_{1}-t_{2})O_{j}U(t_{2}-t_{0})]$
$\displaystyle=\frac{1}{Z_{\beta}}\mathrm{Tr}\,[e^{-\beta
H}U^{\dagger}(t_{1}-t_{0})O_{i}$ $\displaystyle\qquad\qquad
U(t_{1}-t_{0})U(t_{0}-t_{2})O_{j}U(t_{2}-t_{0})]$
$\displaystyle=\frac{1}{Z_{\beta}}\mathrm{Tr}\,[e^{-\beta H}e^{\beta
H}U(t_{0}-t_{2})O_{j}U(t_{2}-t_{0})e^{-\beta H}U^{\dagger}(t_{1}-t_{0})$
$\displaystyle\qquad\qquad O_{i}U(t_{1}-t_{0})]$
$\displaystyle=\frac{1}{Z_{\beta}}\mathrm{Tr}\,[e^{-\beta
H}U(t_{0}-t_{2}+i\hbar\beta)O_{j}U(t_{2}-i\hbar\beta-
t_{1})O_{i}U(t_{1}-t_{0})]$ $\displaystyle=\langle
O_{j}(t_{2}-i\hbar\beta)O_{i}(t_{1})\rangle.$ (111)

Now, applying the following operator identities:

$\displaystyle AB$ $\displaystyle=\frac{1}{2}\\{A,B\\}+\frac{1}{2}[A,B]$ (112)
$\displaystyle BA$ $\displaystyle=\frac{1}{2}\\{A,B\\}-\frac{1}{2}[A,B]$ (113)

and the definitions of the Green’s functions in (100), the identity (111)
becomes

$\displaystyle G^{S}_{ij}(t_{1}-t_{2})+\frac{1}{2}\Delta_{ij}(t_{1}-t_{2})$
$\displaystyle=G^{S}_{ij}(t_{1}-t_{2}+i\hbar\beta)-\frac{1}{2}\Delta_{ij}(t_{1}-t_{2}+i\hbar\beta).$
(114)

Its Fourier transform is simply:

$\tilde{G}^{S}_{ij}(\omega)=\frac{\hbar}{2}\coth{\frac{\hbar\beta\omega}{2}}\tilde{\Delta}_{ij}(\omega),$
(115)

which, in terms of the retarded Green’s function, reduces to the familiar
quantum FDT:

$\displaystyle\tilde{G}^{S}_{ij}(\omega)$
$\displaystyle=\hbar\coth{\frac{\hbar\beta\omega}{2}}\mathrm{Im}\tilde{G}^{R}_{ij}(\omega)$
(116)

since

$\mathrm{Im}\tilde{G}^{R}(\omega)=\frac{\tilde{\Delta}(\omega)}{2\hbar}.$
(117)

Note that the factor 2 dividing is due to the retarded Green’s function
defined only for positive time, see D.

Finally, in the classical or high temperature limit, we obtain the classical
FDT:

$\displaystyle\tilde{G}^{S}_{ij}(\omega)=\frac{2}{\beta\omega}\mathrm{Im}\tilde{G}^{R}_{ij}(\omega),\quad(\hbar\rightarrow
0\;\text{and/or}\;\beta\rightarrow 0).$ (118)

##  Appendix C Generalized Langevin equation

A quantum or classical harmonic oscillator coupled linearly to a harmonic bath
is an exactly solvable model [14, 12], and, as discussed in section 3, can be
described in terms the generalized Langevin equation:

$\displaystyle
M\ddot{q}(t)+M(\Omega^{2}+\delta\Omega^{2})q(t)+\int_{0}^{t}d\tau\mu(t-\tau)q(\tau)=F(t)$
(119)

with the initial conditions:

$q(0)=q_{0},\quad\dot{q}(0)=\dot{q}_{0}.$ (120)

and a Gaussian noise:

$\displaystyle\langle F(t)\rangle$ $\displaystyle=0$ (121)
$\displaystyle\frac{1}{2}\braket{\\{F(t),F(t^{\prime})\\}}$
$\displaystyle=\nu(t-t^{\prime}).$ (122)

###  C.1 Solution

It is common in the literature to find the generalized Langevin equation in
terms of the damping kernel $\gamma$, aka the memory function, which is
defined as:

$\mu(t-\tau)=M\frac{\partial}{\partial t}\gamma(t-\tau)$ (123)

then:

$\ddot{q}(t)+\Omega^{2}q(t)+\int_{0}^{t}d\tau\gamma(t-\tau)\dot{q}(\tau)+\gamma(t)q(0)=F(t)/M$
(124)

The generalized Langevin equation (124) can be solved using Laplace transform,
which gives:

$(s^{2}+\Omega^{2}+s\hat{\gamma}(s))\hat{q}(s)=sq_{0}+\dot{q}_{0}+\hat{F}(s)/M$
(125)

from which the solution is:

$\hat{q}(s)=M(s\,q_{0}+\dot{q}_{0})\hat{G}(s)+\hat{F}(s)\hat{G}(s)$ (126)

where the Laplace transform of the (retarded) Green’s function is

$\hat{G}(s)=\dfrac{1/M}{s^{2}+\Omega^{2}+s\hat{\gamma}(s)},$ (127)

which, in time space, must satisfy the initial conditions

$G(0)=0,\quad\dot{G}(0)=\frac{1}{M}$ (128)

that fully determine it.

The solution back to the time space is:

$q(t)=M\left(q_{0}\dot{G}(t)+\dot{q}_{0}G(t)\right)+\int_{0}^{t}d\tau
G(t-\tau)F(\tau).$ (129)

###  C.2 Stationary limit

The general solution (129) simplifies in the stationary limit, reducing to:

$\displaystyle q_{s}(t)=\int_{0}^{\infty}d\tau
G(t-\tau)F(\tau),\quad(t\rightarrow\infty),$ (130)

that means for non-vanishing initial conditions, we require
$G(t\longrightarrow\infty)=0$ and $\dot{G}(t\longrightarrow\infty)=0$.
However, if $q_{0}=0$ then $G(t\longrightarrow\infty)=0$ is sufficient to
reach stationarity.

The causality condition of the Green’s function (19) allows the integration to
be extended to the full real axis. Hence, upon Fourier transform, the
stationary solution (130) is just a product:

$\tilde{q}_{s}(\omega)=\tilde{G}(\omega)\tilde{F}(\omega).$ (131)

The gFDR (50) can be readily obtained using the definition of the correlation
function and the stationary solution (131).

##  Appendix D Fourier transform

Consider the following definitions of Fourier transform:

$\displaystyle\tilde{f}(\omega)$
$\displaystyle\equiv\int_{-\infty}^{\infty}dt\,e^{-i\omega t}f(t)$ (132)
$\displaystyle\tilde{f}_{+}(\omega)$
$\displaystyle\equiv\int_{0}^{\infty}dt\,e^{-i\omega t}f(t)$ (133)
$\displaystyle\tilde{f}_{-}(\omega)$
$\displaystyle\equiv\int^{0}_{-\infty}dt\,e^{-i\omega t}f(t)$ (134)

such that

$\displaystyle\tilde{f}(\omega)$
$\displaystyle=\tilde{f}_{-}(\omega)+\tilde{f}_{+}(\omega).$ (135)

It follows that:

  1. 1.

For functions even in time $f(t)=f(-t)$:

$\tilde{f}_{-}(\omega)=\tilde{f}_{+}(-\omega)\Longrightarrow\tilde{f}(\omega)=\tilde{f}(-\omega)$
(136)

  2. 2.

For functions odd in time $f(t)=-f(-t)$:

$\tilde{f}_{-}(\omega)=-\tilde{f}_{+}(-\omega)\Longrightarrow\tilde{f}(\omega)=-\tilde{f}(-\omega)$
(137)

  3. 3.

If $\tilde{f}_{-}(\omega)=\tilde{f}_{+}(\omega)$, then
$\tilde{f}(\omega)=2\tilde{f}_{+}(\omega)$.

  4. 4.

If $g(t)=f(t)\theta(t)$, then $\tilde{g}(\omega)=\tilde{f}_{+}(\omega)$.

## References

  * [1] J. Schwinger, Brownian motion of a quantum oscillator, Journal of Mathematical Physics 2, 407–432, 1961.
  * [2] L. V. Keldysh, Diagram technique for nonequilibrium processes, Sov. Phys. JETP 20, 1018–1026, 1965.
  * [3] A. Kamenev, _Field theory of non-equilibrium systems_.  Cambridge University Press, 2023.
  * [4] L. M. Sieberer, A. Chiocchetta, A. Gambassi, U. C. Täuber and S. Diehl, Thermodynamic equilibrium as a symmetry of the schwinger-keldysh action, 92, Physical Review B, 2015.
  * [5] C. Aron, G. Biroli and L. F. Cugliandolo, Symmetries of generating functionals of langevin processes with colored multiplicative noise, Journal of Statistical Mechanics: Theory and Experiment 2010, P11018, 2010.
  * [6] C. Aron, G. Biroli and L. F. Cugliandolo, (non) equilibrium dynamics: a (broken) symmetry of the keldysh generating functional, SciPost Physics 4, 008, 2018.
  * [7] L. F. Cugliandolo, The effective temperature, Journal of Physics A: Mathematical and Theoretical 44, 483001, 2011.
  * [8] R. Feynman and F. Vernon, The theory of a general quantum system interacting with a linear dissipative system, Annals of Physics 24, 173, 1963.
  * [9] J. T. Hsiang, C. H. Chou, Y. Subaşı and B. L. Hu, Quantum thermodynamics from the nonequilibrium dynamics of open systems: Energy, heat capacity, and the third law, Phys. Rev. E 97, 012135, 2018, [arXiv:1703.04970 [quant-ph]].
  * [10] A. O. Caldeira and A. J. Leggett, Path integral approach to quantum Brownian motion, Physica A 121, 587–616, 1983.
  * [11] B. L. Hu, J. P. Paz and Y.-h. Zhang, Quantum Brownian motion in a general environment: I. Exact master equation with nonlocal dissipation and colored noise, Phys. Rev. D 45, 2843–2861, 1992.
  * [12] G. W. Ford, J. T. Lewis and O’Connell, Quantum langevin equation, Phys. Rev. A 37, 4419, 1988.
  * [13] A. O. Caldeira, _An introduction to macroscopic quantum phenomena and quantum dissipation_.  Cambridge University Press, 2014.
  * [14] C. Fleming, A. Roura and B. Hu, Exact analytical solutions to the master equation of quantum brownian motion for a general environment, Annals of Physics 326, 1207–1258, 2011, [arXiv:1004.1603 [quant-ph]].
  * [15] H. Ness, A. Genina, L. Stella, C. Lorenz and L. Kantorovich, Nonequilibrium processes from generalized langevin equations: Realistic nanoscale systems connected to two thermal baths, Physical Review B 93, 174303, 2016.
  * [16] F. Zamponi, F. Bonetto, L. F. Cugliandolo and J. Kurchan, A fluctuation theorem for non-equilibrium relaxational systems driven by external forces, Journal of Statistical Mechanics: Theory and Experiment 2005, P09013, 2005.
  * [17] L. F. Cugliandolo, J. Kurchan and G. Parisi, Off equilibrium dynamics and aging in unfrustrated systems, Journal de Physique I 4, 1641–1656, 1994.
  * [18] R. C. Verstraten, R. F. Ozela and C. M. Smith, Time glass: A fractional calculus approach, Phys. Rev. B 103, L180301, 2021, [arXiv:2006.08786 [cond-mat]].
  * [19] G. Ford and R. O’Connell, Anomalous diffusion in quantum brownian motion with colored noise, Phys. Rev. A 73, 032103, 2006.
  * [20] J.-T. Hsiang and B.-L. Hu, Nonequilibrium nonlinear open quantum systems: Functional perturbative analysis of a weakly anharmonic oscillator, Phys. Rev. D 101, 125002, 2020, [arXiv:1912.12803 [hep-th]].
  * [21] J.-T. Hsiang and B.-L. Hu, Fluctuation-Dissipation Relation from the Nonequilibrium Dynamics of a Nonlinear Open Quantum System, Phys. Rev. D 101, 125003, 2020, [arXiv:2002.07694 [hep-th]].
  * [22] M. Crossley, P. Glorioso and H. Liu, Effective field theory of dissipative fluids, JHEP 09, 095, 2017, [arXiv:1511.03646 [hep-th]].
  * [23] X. Chen-Lin, L. V. Delacrétaz and S. A. Hartnoll, Theory of diffusive fluctuations, Phys. Rev. Lett. 122, 091602, 2019, [arXiv:1811.12540 [hep-th]].
  * [24] A. Jain, P. Kovtun, A. Ritz and A. Shukla, Hydrodynamic effective field theory and the analyticity of hydrostatic correlators, 2021, Journal of High Energy Physics, 2021.
  * [25] T. Hayata and Y. Hidaka, Diffusive Nambu-Goldstone modes in quantum time-crystals, 2018, [arXiv:1808.07636 [hep-th]].
  * [26] E. Wang and U. W. Heinz, A Generalized fluctuation dissipation theorem for nonlinear response functions, Phys. Rev. D 66, 025008, 2002, [arXiv:hep-th/9809016].
  * [27] J.-T. Hsiang and B.-L. Hu, Fluctuation-dissipation relation for open quantum systems in a nonequilibrium steady state, Phys. Rev. D 102, 105006, 2020, [arXiv:2007.00906 [quant-ph]].
  * [28] Y. BenTov, Schwinger-Keldysh path integral for the quantum harmonic oscillator, 2021, [arXiv:2102.05029 [hep-th]].
  * [29] H. Liu and P. Glorioso, Lectures on non-equilibrium effective field theories and fluctuating hydrodynamics, PoS TASI2017, 008, 2018, [arXiv:1805.09331 [hep-th]].
  * [30] R. Kubo, Statistical mechanical theory of irreversible processes. I. General theory and simple applications in magnetic and conduction problems, J. Phys. Soc. Jap. 12, 570–586, 1957.
  * [31] P. C. Martin and J. Schwinger, Theory of many-particle systems. i, Phys. Rev. 115, 1342–1373, 1959.