****** Start of file apssamp.tex ******

# Learn one size to infer all: Exploiting translational symmetries in delay-
dynamical and spatio-temporal systems using scalable neural networks

Mirko Goldmann<EMAIL_ADDRESS>Claudio R. Mirasso  Ingo Fischer
<EMAIL_ADDRESS>Miguel C. Soriano<EMAIL_ADDRESS>Instituto
de Física Interdisciplinar y Sistemas Complejos (IFISC, UIB-CSIC), Campus
Universitat de les Illes Balears E-07122, Palma de Mallorca, Spain

###### Abstract

We design scalable neural networks adapted to translational symmetries in
dynamical systems, capable of inferring untrained high-dimensional dynamics
for different system sizes. We train these networks to predict the dynamics of
delay-dynamical and spatio-temporal systems for a single size. Then, we drive
the networks by their own predictions. We demonstrate that by scaling the size
of the trained network, we can predict the complex dynamics for larger or
smaller system sizes. Thus, the network learns from a single example and, by
exploiting symmetry properties, infers entire bifurcation diagrams.

##  I Introduction

Due to their many degrees of freedom and potentially multiple time scales, the
prediction and analysis of complex system dynamics represent challenging
tasks. Tackling unavailable analytical models, machine learning methods
emerged [1, 2, 3] that learn from data and forecast the dynamics of
biological[4, 5], climate[6], spatio-temporal[7, 8] and other complex
systems[9, 10, 11]. More recently, incorporating symmetries is considered to
guide such data-driven models towards preserving conservation laws [12],
improve prediction accuracy [13] and might yield efficient representations of
learned processes [14]. Nonetheless, machine learning demands large amounts of
data, and yet predictions are mainly restricted to the dynamical regime they
observed during training. In recent works[15, 16, 17, 18], parameterized
neural networks are being studied for the prediction of untrained dynamics.
Thereby, a neural network needs to be trained on several examples, often
covering different dynamical regimes. After training, changing the
parameterization renders the network the ability to generalize and predict
untrained dynamics, transformations of chaotic attractors, and close-by
bifurcations. However, in real-world applications, system parameters can often
not be easily changed, and consequently, certain parameter regimes are hardly
accessible. Here, the question arises - Is it possible to infer untrained
(size-dependent) dynamical regimes of a complex system while learning from one
example related to a certain system size only?

In this manuscript, we exploit symmetries in dynamical systems by designing
neural networks that exhibit i) excellent generalization properties allowing
us to infer untrained dynamical properties for a wide tuning of bifurcation
parameters and ii) high prediction capabilities. Here, we design scalable
neural networks that satisfy the same symmetries as the dynamical system under
study. We train these networks on a single time series of either delay-
dynamical or spatio-temporal systems with a fixed size. After training,
notably, the size of the trained network can be scaled up or down by adding or
removing neurons exploiting the translational symmetry in these systems.
Without further adaptations, the scaled network generalizes from learned to
untrained dynamics and enables far-reaching inferences revealing bifurcations
to various dynamical regimes. In the following, we first focus on delay-
dynamical systems and provide a single neural network that infers the entire
bifurcation diagram for changing the delay while learning from a single
example and exploiting temporal translational symmetry. Based on the analogies
between delay-dynamical and spatio-temporal systems, can interpret the
symmetry in delay-dynamical systems also as a quasi-spatial symmetry[19, 2,
1]. For comparison, in Section III, we apply our approach analogously to the
real spatial translational symmetry of a Kuramoto-Sivashinsky model with
periodic boundary conditions. Using a parallel network architecture we can
train the network on a single spatial size and then infer untrained spatial
extensions.

Figure 1: Scheme of the delayed echo state network. The reservoir contains
delayed internode weighted connections (green), where the delay $D$ of the
nodes connections can be adjusted.  Figure 2: Two-dimensional projection
(x-axis $y(t)$, y-axis $y(t-\tau)$) of attractors of the chaotic Mackey-Glass
system for different delay lengths in a), b), c), & d). e) Bifurcation diagram
generated using the Mackey-Glass delay system. Inferred attractors by the dESN
trained on a single example of the Mackey-Glass system with $\tau=100$ shown
in f), g), h) & i). j) Bifurcation diagram inferred by the dESN trained on
data of a Mackey-Glass system with $\tau=100$.

##  II Inference of Delay System Dynamics

Dynamical systems with delayed coupling of system variables play an important
role in many real-world contexts, such as climate systems[22], epidemiological
models[23], biological systems[24], control systems[25, 26, 27] and photonic
systems[28, 29, 30, 31]. Since the evolution of delay systems relies on a
continuous history function $h(t),\ t\in[0,-\tau]$, the phase space of these
systems is inﬁnite dimensional[32]. Depending on the length of the delay and
other system parameters, these systems either converge to fixed points, limit
cycles, or evolve on chaotic attractors[31, 33, 29, 30]. Furthermore, with
increasing delay length, these systems exhibit extensive chaotic behavior,
i.e., the maximal dimension of their chaotic attractors scales linearly with
the delay time [34].

In the following, we design echo state networks that incorporate a delay in
their topology, as illustrated in Fig. 1. Due to this delay, the network
exhibits the same temporal translational symmetry as the delay-dynamical
system it is built to predict[19]. During training, the network delays will be
tuned to match those of the target system, and the readout weights are adapted
to implement one-step-ahead prediction. After learning the time evolution
operator of the delay system, the network can be scaled up or down by
extending or shortening its inherent delay to infer the untrained dynamics
corresponding to shorter or longer delay, respectively. The dynamical
evolution of the delayed echo state network (dESN) is governed by the
following equation:

$\displaystyle\vec{x}(n{+}1)=$
$\displaystyle\alpha\vec{x}(n){+}\beta\tanh(\mathbf{W}\vec{x}(n{-}D)+\gamma\mathbf{W}_{in}s(n)+\mathbf{W}_{b})$
(1)

where $\vec{x}(n)\in\mathbb{R}^{K}$, $K$ is the network network size, $\alpha$
is the leak term, $\beta$ is a feedback gain, $\gamma$ is the input gain, and
$D$ is the delay length. The total size of the network scales linearly with
$D$, as all previous states $x(n)$ up to $x(n-D)$ must be stored. Thus, the
delay can be seen as a signal through a cue of $D$ hidden linear neurons along
each delayed connection. However, the read-out dimension of the network is
independent of the delay length since $\mathbf{W}_{out}\in\mathbb{R}^{K}$. The
randomly drawn matrix $\mathbf{W}_{in}\in\mathbb{R}^{K}$ gives the connection
between the input $s(n)\in\mathbb{R}$ and the network, while
$\mathbf{W}_{b}\in\mathbb{R}^{K}$ gives a random bias to each node. The
delayed connections are weighted via $\mathbf{W}\in\mathbb{R}^{K\times K}$,
the elements of which are randomly drawn from a uniform distribution
$\mathcal{U}[-1,1]$, with a sparsity of $1.5\%$.

The training data set used only contains a time series of a delay system with
a single fixed delay. Here, we consider inferring the dynamics of a Mackey-
Glass delay system[35], however, the same holds for other delay systems as we
demonstrate in Appendix D by performing a similar study for an Ikeda-type
delay system. The dynamics of the Mackey-Glass (MG) delay system is given by
the following equation:

$\dot{s}(t)=-\frac{s(t)}{T_{0}}+\dfrac{0.2s(t-\tau)}{1+s(t-\tau)^{10}}.$ (2)

We set the characteristic relaxation time $T_{0}=10$ and the system delay
$\tau=100$ and generate a time series $s(t)$ that is sampled with $\Delta t=1$
to obtain the training data set $s(n)$. For these parameters, the system
evolves along a chaotic attractor as depicted in Fig. 3.

Figure 3: a) Time series of the original Mackey Glass system (blue) and the
autonomously continued (green) chaotic attractor using a dESN. b) Divergence
rate of the chaotic system and the time series generated by the dESN (green)
as shown in a), the orange line indicates the average divergence rate of 20
different dESNs initialized at 20 different trajectories of the chaotic
Mackey-Glass system. The blue line indicates the divergence rate related to
the largest Lyapunov exponent $\lambda=0.009$ of the Mackey-Glass system with
$\tau=30$. In c) and d), two-dimensional projection of the chaotic attractor
of the Mackey-Glass system with a delay of $\tau=30$ and the dESN prediction
with $D=30$, respectively. Figure 4: The dynamics of the spatio-temporal
Kuramoto-Sivashinsky system $y(x,n)$ (top row), the predicted dynamics from a
parallel network architecture $\hat{y}(x,n)$ (mid row), and their difference
$y(x,n)-\hat{y}(x,n)$(bottom row). The parallel networks are trained with data
from the Kuramoto-Sivashinsky system of spatial extension $L=10\pi$ (red box),
by adapting the network architecture it generates the dynamics also for
smaller and larger spatial extensions (yellow box). For comparison both
systems were initialized with the same initial conditions.

Before training the reservoir on the MG time series, the delay $D$ of the dESN
must be adjusted to the delay of the target system. We find the optimal
performance when we match the delay of the reservoir connections directly with
the delay of the MG system, i.e. by setting $D=\tau/\Delta t=100$. In most
real-world systems, however, the delay is not known a priori and needs to be
extracted from the data. As we show in the Appendix A, we can optimize the
reservoir delay $D$ by determining the best performance in a one-step ahead
prediction that reveals the delay underlying the training data. All other
hyperparameters of the dESN are optimized using Bayesian optimization and are
given in table 1.

parameter | symbol | value
---|---|---
network size | $K$ | 1000
initial steps | $N_{init}$ | 5000
training steps | $N_{train}$ | 25000
feedback gain | $\beta$ | 0.176
input gain | $\gamma$ | 1.24
spectral radius | $\rho$ | 0.84
leak rate | $\alpha$ | 0.75
Table 1: List of parameters used for learning the Mackey-Glass system.

During training, the dESN is driven by the input $s(n)$ with added white
noise. The addition of noise during training was shown to improve the
stability of the reservoir in closed loop mode[36]. The network’s output layer
is optimized using linear regression to predict the next step in the training
time series. Once the output weights $W_{out}$ are computed, the reservoir is
decoupled from the input $s(n)$ and the prediction of the reservoir
$\hat{y}(n)=\mathbf{W}_{out}\vec{x}(n)$ is returned as the new input.
Consequently, the reservoir evolves autonomously in closed-loop mode[37, 38,
39]. Reservoir computing [2, 40, 41, 42] is applied here, as it offers fast
and efficient training; nevertheless, we expect that other training methods
yield comparable results.

In closed-loop mode, the reservoir evolves along the attractor on which it was
trained, as depicted in Figs. 2 a) and f). When resetting the delay $D$ of the
dESN connections, we do not need any additional data for initialization, no
further retraining, and the output weights can be kept fixed. As shown in
Figs. 2 g) to i), the dESN generates high-dimensional chaotic attractors at
$D=17,\ 200$ and a limit cycle at $D=39$. These inferred dynamical states are
similar to those of the original MG system with delay $\tau=17,\ 39$ and $200$
shown in Figs. 2 b) to d). Consequently, it is sufficient that the dESN
learned from a single time series of the MG system with delay $\tau=100$, to
infer the dynamics for other delay lengths. By scanning the delay of the dESN
from $D=1-120$ we generate a bifurcation diagram of the inferred dynamics as
shown in Fig. 2 j). Depending on $D$, the dESN exhibits chaotic behavior,
transitions to intermittent limit cycles, for example, at $D=34\ \&\ 39$,
period-doubling bifurcations at $D=14\ \&\ 16$, and transitions to a stable
fixed point at $D=5$. The observed transitions in the dESN coincide with those
found for the original MG system shown in Fig. 2 e). Furthermore, in Appendix
C we show that the dESN can infer multistabilities in the delay range
$\tau,D=70-80$. Röhm et al.[43] recently presented that reservoir computing is
capable of predicting multistabilities in a dynamical system even if the
reservoir is trained only around one of the attractors. Here, our dESN infers
infers coexisting attractors, while it was trained in a dynamical regime where
there was no multistability.

In Fig. 2 we show that the dESN trained at the MG system with $\tau=100$ can
infer dynamics of the MG system with other delay lengths. In the following, we
further underline how precise the dESN infers the dynamics even far from the
training example. Therefore, in Fig. 3, we show the inference of the dESN with
$D=30$ compared to the original dynamics. We divide the prediction
capabilities into two regimes, which we term weather and climate[36]. Here,
the weather regime refers to the short-term behavior, when going from
externally driven operation to autonomous operation. As shown in Fig.3 a), the
dESN can precisely continue the trajectory of a MG system with $\tau=30$ for
around $1000$ steps, which corresponds to approximately nine Lyapunov times.
Due to the chaotic behavior of the MG system with $\tau=30$, two closely
initialized trajectories diverge in time, where the divergence rate is given
by the largest Lyapunov exponent. As depicted in Fig. 3 b), the divergence
rate of the predicted trajectory (green line) is similar to the largest
Lyapunov exponent of the MG system (blue line). Furthermore, the performance
that the dESN provides is robust against randomization. We find similar
results (orange line) for 20 different dESNs with varied connection matrices
and training data. The climate regime describes the long-term behavior of a
dynamic system, e.g. the evolution along the chaotic attractor. In Fig. 3 d),
we show the predicted trajectory after it diverges from the initialization
trajectory. The dESN reproduces the climate of the delay system by inferring
the chaotic attractor of the MG system with $\tau=30$ as shown in Fig. 3 c)
and is also the case for other delays $D$ shown in Fig. 2. Accordingly,
training the dESN only on a single example is sufficient to precisely infer
the dynamical properties of the MG system for various delay lengths.

As presented in Figs. 2 and 3, training the dESN with data from the MG system
with a delay of $\tau=100$ enables inferring untrained dynamics related to
shorter and even longer delays. A delay of $\tau=100$ places the system in the
long delay limit, where the delay is much longer than the characteristic
relaxation time $T_{0}$ in Eq. 2. Here, learning benefits from a clear time-
scale separation between response and delay time. We find that learning in the
short delay regime enables inferring towards relatively shorter delays,
whereas going towards larger delays, the inference quality deteriorates. For a
more detailed discussion, we refer to Appendix B.

##  III Inference of Spatio Temporal system dynamics

Taking into account the close relationship of delay and spatio-temporal
systems [2, 1, 19], we extend our approach to the spatial translational
symmetry of a homogeneous Kuramoto-Sivashinsky (KS) model with periodic
boundary conditions. By analogy, we design reservoirs that can infer
attractors for different system sizes when trained for a single size only. The
dynamics of the KS model are governed by:

$y_{t}=-yy_{x}-y_{xx}-y_{xxxx},$ (3)

where $y(x,t)$ is a scalar field. We further consider periodic boundary
conditions in the interval $[0,L)$. For the training, we generate a data set
using a spatial extension of $L=10\pi$. The generated data is sampled every
$\Delta t=0.25$ in the temporal domain and contains $Q=100$ equidistant
samples in the spatial domain leading to $y(n)\in\mathbb{R}^{100}$. To predict
the spatio-temporal evolution of the KS system, we construct a parallel
network architecture as described by Pathak et al.[8] and train this
architecture to perform a one-step-ahead prediction. The parallel architecture
contains $G=10$ subnetworks, each containing $K=1000$ neurons, respectively.
The evolution of the $g$-th subnetwork can be described by the following
equation:

$x^{g}(n+1)=\tanh(Wx^{g}(n)+\gamma W_{in}u^{g}(n)+0.2W_{b}),$ (4)

where $W\in\mathbb{R}^{N\times N}$ is a randomly drawn adjacency matrix with
spectral radius $\rho=1.3$, $W_{in}\in\mathbb{R}^{N\times M}$ are randomly
drawn input weights, the input gain $\gamma=0.001$, and
$W_{b}\in\mathbb{R}^{N}$ is a random bias. The input $u^{g}(n)$ of the
subnetworks is generated by dividing the data from the spatial domain into $G$
sections of the same size $Q/G=10$ (spatial size of $\pi$). Each of the $G$
subnetworks receives the 10 inputs of one particular section and additionally
the 3 closest inputs from both neighboring sections, leading to an input
dimension of $M=16$. Due to the spatial translational symmetry of the
homogeneous KS system along the spatial domain the subnetworks adjacency
matrix and hyperparameters can be chosen identical. Similar approaches are
mentioned in [8, 44]. During training, each subnetwork is trained to perform a
one-step-ahead prediction of the spatial domain section to which it
corresponds, resulting in 10 outputs per network computed using the output
weight matrix $W_{out}\in\mathbb{R}^{10\times 1000}$. Again, due to the
translational symmetry the output weights sets can be shared between the
subnetworks. After the training phase, we close the loop by feeding the
predicted state back to the reservoir. As shown in Figs. 4 a) and f), parallel
reservoirs autonomously predict the weather and climate of the chaotic KS
system with $L=10\pi$.

To infer spatial extensions $L$ that are different from the one used during
training, we take further advantage of the spatial translation symmetry and
similar to changing the delay in the dESN, we adapt the topology of the
reservoir. Therefore, either a subnetwork is removed from the parallel
architecture or a copy of a subnetwork is inserted into the architecture. By
varing the number of subnetworks $G$ in the architecture the output dimension
alters in steps of $10$ which effectively increases the predicted spatial
extensions. Due to the spatial length trained on and the used number of
subnetworks, we can vary the spatial extension of the predicted system in
units of $\pi$. In Fig. 4, we present the original KS dynamics (first row),
the predicted dynamics of the reservoir (middle row), and the difference
between both (bottom row) for different spatial extensions, respectively.
Thereof, we initialized both using the same initial conditions. In analogy to
the results for the MG delay-dynamical system, we observe that the once-
trained reservoir is able to infer the untrained dynamics at significantly
shorter and larger spatial extensions. Therefore, the parallel network
architecture takes advantage of the spatial translational invariance of the
homogeneous KS model with periodic boundary conditions.

##  IV Conclusion

We demonstrated that exploiting symmetries of dynamical systems by designing
neural networks obeying these symmetries improves their prediction ability and
further enable far-reaching inference. We particularly showed this for
translational symmetry in delay-dynamical and spatio-temporal systems.
Transformation of the trained networks along this symmetry enabled learning
from a single example to infer an entire bifurcation diagram. Thus, we obtain
minimal requirements for training data and gain a single model that can infer
a wide variety of dynamical behaviors. This represents a very efficient use of
resources, such as the required network size, training data, and energy.
Therefore, it is also a step in the direction of more sustainable machine
learning. In addition, the provided method might be used to analyze real-world
systems for which certain parameter settings might not be accessible.
Recently, Tegmark et al.[45] presented a machine learning method to manifest
hidden symmetries in physical systems. The discovery of hidden symmetries and
the use of our approach to exploit them for far-reaching inferences while
learning from single examples could give rise to powerful predictive models
for a variety of high-dimensional systems, including complex networks[9].

##  V appendix

##  Appendix A Delay estimation from time series using delayed echo state
networks

Figure 5: Estimating the delay a from time series of the chaotic Mackey-Glass
system for three different delay lengths $\tau$. The red line marks the delay
$\tau$ that underlies the data set a) $\tau=30$, b) $\tau=60$ and c)
$\tau=100$. The blue dots mark the one-step-ahead prediction accuracy of dESNs
with different delays $D$ (x-axis) and different sets of hyperparameters. The
accuracy is determined by the NRMSE (normalized root mean square error) of the
one step ahead prediction.

As mentioned in section II, during training, we match the delay of the delayed
echo state network (dESN) with the delay of the Mackey-Glass system to enhance
the prediction abilities of the dESN. As the Mackey-Glass time series is
sampled with $\Delta t=1$ this leads to $D=\tau$. In many real-world systems,
however, the underlying delay is not known a priori and needs to be determined
from a sampled time series of the dynamical system under study. There exist
several methods to estimate delays from time series based on, e.g.,
autocorrelation function, delayed mutual information, local linear fitting in
a low-dimensional subspace[46], and even deep learning-based methods [47]. In
the following, we show how to estimate the delay that underlies a time series
using the delayed echo state network (dESN) given in Eq. 1. Thus, we seek for
high one-step-ahead prediction accuracy of dESNs by scanning the delay $D$ of
the dESN and optimizing the reservoir hyperparameters using Bayesian
optimization respectively. In Fig. 5, we show the results of this scan for
three different delays $\tau=30,60,100$ of the Mackey-Glass system that
generated the training data. The delay scan of the dESN reveals an increased
performance (reduced NRMSE) when its delay is in resonance or equal to the
delay $D=\tau$ of the Mackey-Glass system. As it can be seen in Fig. 5,
identifying these resonances in turn allows to extract the delay that
underlies the data. There can be relative offsets of $\pm 1$ in the optimal
estimated delay, which here is the smallest possible offset related to the
sampling of the Mackey-Glass time series used here ($\Delta t=1$). In this
manner, the optimal setting of the dESN delay might depend on the frequency
used for sampling the real-world system.

##  Appendix B Dependence of inference capabilities on the trained chaotic
attractor

In section II, we show that the dESN with delay $D=100$ trained only on the
Mackey-Glass time series with $\tau=100$ can predict the entire bifurcation
diagram of the Mackey-Glass system by scaling its size (parameter $D$) after
training. This means that the scalable dESN is able to infer dynamics of much
shorter but also longer delays while being trained on data of a single delay
length only. In the following, we evaluate how these prediction capabilities
depend on the properties of the chaotic attractor used for the training of the
dESN. To do this, we use the absolute difference between the autocorrelation,
denoted $\Delta_{ACF}$, of the Mackey-Glass attractor and the attractor
predicted by the respective dESN to quantify the quality of the prediction.
Here, we computed the autocorrelation (ACF) of the original and predicted time
series for 1100 steps into the past, as shown in Fig. 6. Subsequently, the
absolute difference $\Delta_{ACF}$ between the original and predicted ACF is
calculated and summed in the range $n\in[0,1100]$.

Figure 6: Comparison of the autocorrelation function ACF of the Mackey-Glass
system (blue) with a) $\tau=30$ and b) $\tau=100$ and the predicted time
series using the dESN (green). The dESN used in this plot was trained on a
chaotic Mackey-Glass system with $\tau=100$.

For comparison, we train three networks where $D_{T}=17,\ 30\ \&\ 100$ defines
the delay during training, respectively, and the corresponding training data
set is given by the Mackey-Glass time series with similar corresponding delay
$\tau=17,\ 30\ \&\ 100$. In Fig. 7 we show the autocorrelation difference
$\Delta_{ACF}$ for two ranges: a) $\tau,D=17-100$ and b) $\tau,D=100-900$.
While all trained dESN can predict towards smaller delays, the predicted
attractors from the dESNs trained at $D_{T}=17$ (blue) and $D_{T}=30$ (green)
start to deviate from the original Mackey-Glass attractor at a delay $D>30$
and $D>50$, respectively. In contrast, the one trained at $D=100$ (red) does
not show an increase in the difference between the predicted and original
autocorrelation function $\Delta_{ACF}$ even if the delay $D$ of the dESN
becomes much larger than the one used during training (deviations in the range
$\tau,D=75-95$ are caused by multistability, see appendix C). These results
illustrate that the dESN $D_{T}=100$ trained at $\tau=100$ has the highest
capabilities and can predict even attractors with much longer delays.

We relate this ability of the dESN with $D_{T}=100$ to the properties of the
dynamics observed during training. As mentioned in the section II, a delay of
$\tau=100$ sets the Mackey-Glass system in the long delay limit, where a clear
separation of the local and the delayed dynamics appear. This can be indicated
by using the autocorrelation function of the time series, as shown in Fig. 6,
wherein panel a) a short delay $\tau=30$ does not show a clear separation
between delay and local dynamics. In contrast, in panel b) of Fig. 6, one can
observe this separation showing a decay of local dynamics until $n=50$ and the
delayed dynamics indicated by the peak around the delay time $n=100$.

In conclusion, exploiting translational symmetries enables one to infer
untrained dynamics while learning from a chaotic time series related to a
single system size only. The prediction ability of the trained scalable neural
network further depends on the (size-dependent) dynamical regimes on which it
was trained. Thereby, the prediction ability of the dESN is the strongest if
the network is trained in the long delay limit of delay-dynamical systems.

Figure 7: Absolute difference of the original and the predicted
autocorrelation $\Delta_{ACF}$ of the chaotic attractors a) in the range
$\tau,D\in[17,100]$ and b) $\tau,D\in[100,900]$. The colors refer to different
delays of the Mackey-Glass system used in the training of the dESN, $\tau=17$
(blue), $\tau=30$ (green) and $\tau=100$ (red). Solid lines indicate the
average over 100 random seeds used to generate the dESN and initialize the
Mackey-Glass system that generated the training data. The shaded areas
indicate the standard deviation, respectively.

##  Appendix C Inferring coexisting limit cycle and chaotic attractor

As mentioned in the manuscript, the Mackey-Glass system exhibits a
multistability leading to coexistence of a limit cycle and a chaotic attractor
in the region $\tau\in[79,95]$. In Fig. 8 a), we show an overlay of the limit
cycle and the chaotic attractor for a delay of $\tau=85$. During the
computation of the bifurcation diagram, depending on their initial conditions,
the trajectory of the Mackey-Glass system and the dESN converge either to one
or the other of the attractors. This effect causes deviations in the
computation of the autocorrelation comparison shown in Fig. 7. In Fig. 8 b),
we quantify the frequency with which the Mackey-Glass system and dESN end up
in the limit cycle, respectively. Therefore, we generate time series of both
systems starting from 100 different initial conditions. Using the
autocorrelation function of the time series generated, we discriminate if the
observed attractor is either chaotic or a limit cycle. As shown in Fig. 8 b),
both the Mackey-Glass system and the dESN exhibit comparable probabilities of
ending in the limit cycle. It is worth mentioning that the dESN was trained
with data of a Mackey-Glass system with $\tau=100$ where there exists no
multistability. The prediction of attractors in a multistable system that was
not part of the training using RC was shown recently by Röhm et al.[43]. The
fact that the dESN can predict these multistabilities even if it is trained on
a Mackey-Glass system with a different delay illustrates how strong the
generalization ability of the dESN is.

Figure 8: a) 2D representation of the coexisting limit cycle and chaotic
attractor generated by the Mackey-Glass system at a delay of $\tau=85$. b)
Number of occurrences of a limit cycle in the delay range $\tau,D\in[70,100]$
for 100 different initial conditions. The blue colored line indicates the
occurrences in differently initialized Mackey-Glass systems, while green
indicates occurrences in the differently initialized dESN systems. The dESN is
trained at $D=100$ on a Mackey-Glass system with $\tau=100$ and afterwards the
delay D was scaled down.

##  Appendix D Inferring Ikeda delay system dynamics

Figure 9:  Two-dimensional projection (x-axis $y(t)$, y-axis $y(t-\tau)$) of
attractors of the chaotic Ikeda system for different delay lengths in a), b),
c), & d). e) Bifurcation diagram generated using the Ikeda delay system.
Inferred attractors by the dESN trained on a single example of the Ikeda
system with $\tau=100$ shown in f), g), h) & i). j) Bifurcation diagram
inferred by the dESN trained on data of a Ikeda system with $\tau=100$.

In the main text, we show the autonomous continuation of a Mackey-Glass delay
system for short, medium, and long delays and the prediction of the
bifurcation diagram for arbitrarily long delays. However, the Mackey-Glass
system is only a single example of a delay dynamical system, and there exists
a variety of such delay systems featuring other nonlinearities. In the
following, we show that the proposed method of using a dESN to predict the
bifurcation diagram is not only restricted to the Mackey-Glass system, but can
be applied to other delay systems. Another well-known delay system, called the
Ikeda delay system, was investigated in the field of optics because it can be
generated using a Mach-Zehnder interferometer and a cavity[29]. Here, we
attempt to learn and infer the dynamics of an Ikeda-type delay system with a
sine-square nonlinearity:

$\dot{s}(t)=-1/T_{0}s(t)+\beta sin^{2}(s(t-\tau)).$ (5)

The state of the delay system is given by $s(t)$, the delay is indicated by
$\tau$, the characteristics relaxation time is fixed at $T_{0}=10$ and the
feedback gain is $\beta=0.4$ in the following. The data set is generated by
integrating the delay differential equation in Eq. (5) and the time series is
sampled with $\Delta t=1$. Due to the delay, this Ikeda system can evolve in
different dynamical regimes including chaos [29, 30]. In the following, we
will train a dESN to continue the time series of the Ikeda delay system with a
delay of $\tau=100$. Afterward, the trained dESN will be used to infer the
entire bifurcation of the Ikeda system.

As in the case of the Mackey-Glass system, the delay $D$ of the dESN and the
delay $\tau$ of the Ikeda system to be learned are matched. Afterwards, the
other hyperparameters of the dESN are optimized using a Bayesian optimization
approach. The found optimal parameters are given in Tab. 2. The optimal leak
rate to predict the Ikeda delay system is found at $\alpha=0$, which means
that the leak term in the equation of the dESN can be neglected.

parameter | symbol | value
---|---|---
network size | $K$ | 1000
initial steps | $N_{init}$ | 1000
training steps | $N_{train}$ | 20000
feedback gain | $\beta$ | 0.1
input gain | $\gamma$ | 1.71429
spectral radius | $\rho$ | 0.6975
leak rate | $\alpha$ | 0
Table 2: Hyperparameters used to continue the time series of the Ikeda delay
system.

In Fig. 9 f), we show the results of the dESN trained on an Ikeda system with
$\tau=100$ and compare it to the original attractor in Fig. 9 a). We find that
the dESN properly reproduces the chaotic attractor. Similarly as for the
Mackey-Glass system, we now use the dESN trained on the Ikeda system with
delay $\tau=100$ to infer the entire bifurcation diagram. In Fig. 9, we show a
comparison of the original bifurcation diagram in e) and the one predicted by
the dESN in j). The dESN is able to predict stable fixed points and limit
cycles in the short delay range as well as the route to chaos and the chaotic
attractors for longer delays. Furthermore, it precisely infers intermittent
limit cycles in the chaotic regime for delays $\tau=23-25$. The inferred
bifurcation diagram is again very similar; we even reproduce the fine details
in the bifurcation diagram. Nevertheless, we observe small deviations between
the inferred and the original attractor. Similarly, as for the Mackey-Glass
system, these deviations are caused by the appearance of multistabilities in
the range of $\tau,D=50-55$. Depending on the initial conditions of the dESN
and Ikeda system, respectively, they either evolve on a limit cycle or on a
chaotic attractor.

## ACKNOWLEDGEMENTS

We thank André Röhm for helpful discussions. Authors acknowledge the support
of the Spanish State Research Agency, through the Severo Ochoa and María de
Maeztu Program for Centers and Units of Excellence in R&D (MDM-2017-0711
funded by MCIN/AEI/10.13039/501100011033) and through the QUARESC project
(PID2019-109094GB-C21 and -C22/ AEI / 10.13039/501100011033) and the DECAPH
project (PID2019-111537GB-C21 and -C22/ AEI / 10.13039/501100011033). MG
acknowledges financial support by the European Union’s Horizon 2020 research
and innovation program under the Marie Skłodowska-Curie grant agreement No.
860360 (POST DIGITAL). This work has been partially funded by the European
Union’s Horizon 2020 research and innovation programme under grant agreement
No. 899265 (ADOPD).

## References

  * LeCun  _et al._ [2015] Y. LeCun, Y. Bengio, and G. Hinton, Deep learning, Nature 521, 436 (2015).
  * Jaeger [2001] H. Jaeger, The “echo state” approach to analysing and training recurrent neural networks-with an erratum note, Bonn, Germany: German National Research Center for Information Technology GMD Technical Report 148, 13 (2001).
  * Gauthier  _et al._ [2021] D. J. Gauthier, E. Bollt, A. Griffith, and W. A. Barbosa, Next generation reservoir computing, Nat. Commun. 12, 1 (2021).
  * Elmarakeby  _et al._ [2021] H. A. Elmarakeby, J. Hwang, R. Arafeh, J. Crowdis, S. Gang, D. Liu, S. H. AlDubayan, K. Salari, S. Kregel, C. Richter, _et al._ , Biologically informed deep neural network for prostate cancer discovery, Nature 598, 348–352 (2021).
  * Alber  _et al._ [2019] M. Alber, A. Buganza Tepole, W. R. Cannon, S. De, S. Dura-Bernal, K. Garikipati, G. Karniadakis, W. W. Lytton, P. Perdikaris, L. Petzold, _et al._ , Integrating machine learning and multiscale modeling—perspectives, challenges, and opportunities in the biological, biomedical, and behavioral sciences, NPJ Digital Medicine 2, 1 (2019).
  * Dijkstra  _et al._ [2019] H. A. Dijkstra, P. Petersik, E. Hernández-García, and C. López, The application of machine learning techniques to improve el niño prediction skill, Frontiers in Physics , 153 (2019).
  * Pathak  _et al._ [2017] J. Pathak, Z. Lu, B. R. Hunt, M. Girvan, and E. Ott, Using machine learning to replicate chaotic attractors and calculate lyapunov exponents from data, Chaos 27, 121102 (2017).
  * Pathak  _et al._ [2018] J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach, Phys. Rev. Lett. 120, 024102 (2018).
  * Srinivasan  _et al._ [2022] K. Srinivasan, N. Coble, J. Hamlin, T. Antonsen, E. Ott, and M. Girvan, Parallel machine learning for forecasting the dynamics of complex networks, Phys. Rev. Lett. 128, 164101 (2022).
  * Tang  _et al._ [2020] Y. Tang, J. Kurths, W. Lin, E. Ott, and L. Kocarev, Introduction to focus issue: When machine learning meets complex systems: Networks, chaos, and nonlinear dynamics, Chaos 30, 063151 (2020).
  * Lu  _et al._ [2018] Z. Lu, B. R. Hunt, and E. Ott, Attractor reconstruction by machine learning, Chaos 28, 061104 (2018).
  * Karniadakis  _et al._ [2021] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Physics-informed machine learning, Nat. Rev. Phys. 3, 422 (2021).
  * Barbosa  _et al._ [2021] W. A. Barbosa, A. Griffith, G. E. Rowlands, L. C. Govia, G. J. Ribeill, M.-H. Nguyen, T. A. Ohki, and D. J. Gauthier, Symmetry-aware reservoir computing, Phys. Rev. E 104, 045307 (2021).
  * Higgins  _et al._ [2022] I. Higgins, S. Racanière, and D. Rezende, Symmetry-based representations for artificial and biological general intelligence, Frontiers in Computational Neuroscience , 28 (2022).
  * Kim  _et al._ [2021] J. Z. Kim, Z. Lu, E. Nozari, G. J. Pappas, and D. S. Bassett, Teaching recurrent neural networks to infer global temporal structure from local examples, Nat. Mach. Intell. 3, 316 (2021).
  * Kong  _et al._ [2021a] L.-W. Kong, H.-W. Fan, C. Grebogi, and Y.-C. Lai, Machine learning prediction of critical transition and system collapse, Phys. Rev. Res. 3, 13090 (2021a), 2012.01545 .
  * Kong  _et al._ [2021b] L.-W. Kong, H. Fan, C. Grebogi, and Y.-C. Lai, Emergence of transient chaos and intermittency in machine learning, JPhys Complexity 2, 035014 (2021b).
  * Klos  _et al._ [2020] C. Klos, Y. F. K. Kossio, S. Goedeke, A. Gilra, and R.-M. Memmesheimer, Dynamical learning of dynamics, Phys. Rev. Lett. 125, 088103 (2020).
  * [19] See Supplemental Material at [URL will be inserted by publisher] for a discussion on the symmetries of delay-dynamical and spatio-temporal systems and the analogies between both including two illustrative video sequences.
  * Arecchi  _et al._ [1992] F. Arecchi, G. Giacomelli, A. Lapucci, and R. Meucci, Two-dimensional representation of a delayed dynamical system, Phys. Rev. A 45, R4225 (1992).
  * Giacomelli and Politi [1996] G. Giacomelli and A. Politi, Relationship between delayed and spatially extended dynamical systems, Phys. Rev. Lett. 76, 2686 (1996).
  * Keane  _et al._ [2017] A. Keane, B. Krauskopf, and C. M. Postlethwaite, Climate models with delay differential equations, Chaos 27, 114309 (2017).
  * Bauer  _et al._ [2017] L. Bauer, J. Bassett, P. Hövel, Y. N. Kyrychko, and K. B. Blyuss, Chimera states in multi-strain epidemic models with temporary immunity, Chaos 27, 114317 (2017).
  * Shi  _et al._ [2017] X. Shi, Y. Kuang, A. Makroglou, S. Mokshagundam, and J. Li, Oscillatory dynamics of an intravenous glucose tolerance test model with delay interval, Chaos 27, 114324 (2017).
  * Orosz  _et al._ [2010] G. Orosz, R. E. Wilson, and G. Stépán, Traffic jams: dynamics and control, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 368, 4455 (2010).
  * Stepan [2009] G. Stepan, Delay effects in the human sensory system during balancing, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367, 1195 (2009).
  * Zhu  _et al._ [2019] Q. Zhu, H. Ma, and W. Lin, Detecting unstable periodic orbits based only on time series: When adaptive delayed feedback control meets reservoir computing, Chaos 29, 093125 (2019).
  * Appeltant  _et al._ [2011] L. Appeltant, M. C. Soriano, G. Van Der Sande, J. Danckaert, S. Massar, J. Dambre, B. Schrauwen, C. R. Mirasso, and I. Fischer, Information processing using a single dynamical node as complex system, Nat. Commun. 2, 466 (2011).
  * Larger  _et al._ [2013] L. Larger, B. Penkovsky, and Y. Maistrenko, Virtual chimera states for delayed-feedback systems, Phys. Rev. Lett. 111, 054103 (2013).
  * Erneux  _et al._ [2004] T. Erneux, L. Larger, M. W. Lee, and J.-P. Goedgebuer, Ikeda hopf bifurcation revisited, Physica D 194, 49 (2004).
  * Fischer  _et al._ [1994] I. Fischer, O. Hess, W. Elsä$\beta$er, and E. Göbel, High-dimensional chaotic dynamics of an external cavity semiconductor laser, Phys. Rev. Lett. 73, 2188 (1994).
  * Erneux [2009] T. Erneux,  _Applied delay differential equations_ , Vol. 3 (Springer Science & Business Media, 2009).
  * Heiligenthal  _et al._ [2011] S. Heiligenthal, T. Dahms, S. Yanchuk, T. Jüngling, V. Flunkert, I. Kanter, E. Schöll, and W. Kinzel, Strong and weak chaos in nonlinear networks with time-delayed couplings, Phys. Rev. Lett. 107, 234102 (2011).
  * Le Berre  _et al._ [1986] M. Le Berre, E. Ressayre, A. Tallet, and H. M. Gibbs, High-dimension chaotic attractors of a nonlinear ring cavity, Phys. Rev. Lett. 56, 274 (1986).
  * Mackey and Glass [1977] M. C. Mackey and L. Glass, Oscillation and chaos in physiological control systems, Science 197, 287 (1977).
  * Estébanez  _et al._ [2019] I. Estébanez, I. Fischer, and M. C. Soriano, Constructive role of noise for high-quality replication of chaotic attractor dynamics using a hardware-based reservoir computer, Phys. Rev. App. 12, 034058 (2019).
  * Jaeger and Haas [2004] H. Jaeger and H. Haas, Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication, Science 304, 78 (2004).
  * Antonik  _et al._ [2018] P. Antonik, M. Gulina, J. Pauwels, and S. Massar, Using a reservoir computer to learn chaotic attractors, with applications to chaos synchronization and cryptography, Phys. Rev. E 98, 012215 (2018).
  * Haluszczynski and Räth [2019] A. Haluszczynski and C. Räth, Good and bad predictions: Assessing and improving the replication of chaotic attractors by means of reservoir computing, Chaos 29, 103143 (2019).
  * Maass  _et al._ [2002] W. Maass, T. Natschläger, and H. Markram, Real-time computing without stable states: A new framework for neural computation based on perturbations, Neural Comput. 14, 2531 (2002).
  * Nakajima and Fischer [2021] K. Nakajima and I. Fischer, eds.,  _Reservoir Computing: Theory, Physical Implementations, and Applications_, Natural Computing Series (Springer, Singapore, 2021).
  * Lukoševičius  _et al._ [2012] M. Lukoševičius, H. Jaeger, and B. Schrauwen, Reservoir computing trends, KI-Künstliche Intelligenz 26, 365 (2012).
  * Röhm  _et al._ [2021] A. Röhm, D. J. Gauthier, and I. Fischer, Model-free inference of unseen attractors: Reconstructing phase space features from a single noisy trajectory using reservoir computing, Chaos 31, 103127 (2021).
  * Barbosa and Gauthier [2022] W. A. Barbosa and D. J. Gauthier, Learning spatiotemporal chaos using next-generation reservoir computing, arXiv preprint arXiv:2203.13294 (2022).
  * Liu and Tegmark [2022] Z. Liu and M. Tegmark, Machine learning hidden symmetries, Phys. Rev. Lett. 128, 180201 (2022).
  * Hegger  _et al._ [1998] R. Hegger, M. J. Bünner, H. Kantz, and A. Giaquinta, Identifying and modeling delay feedback systems, Phys. Rev. Lett. 81, 558 (1998).
  * Gao  _et al._ [2021] X. Gao, W. Zhu, Q. Yang, D. Zeng, L. Deng, Q. Chen, and M. Cheng, Time delay estimation from the time series for optical chaos systems using deep learning, Opt. Express 29, 7904 (2021).

Supplemental Material:
Learn one size to infer all: Exploiting translational symmetries in delay-
dynamical and spatio-temporal systems using scalable neural networks

Mirko Goldmann, Claudio R. Mirasso, Ingo Fischer, Miguel C. Soriano
Instituto de Física Interdisciplinar y Sistemas Complejos (IFISC, UIB-CSIC),
Campus Universitat de les Illes Balears E-07122, Palma de Mallorca, Spain

I. Symmetries in delay-dynamical and spatio-temporal systems

In the manuscript, we use the knowledge of symmetries in dynamical systems to
design neural networks, train them for a single size, and subsequently infer
the untrained dynamics of the system for different sizes. In particular, we
infer various untrained dynamical regimes and even bifurcations by exploiting
the temporal translational symmetry in delay-dynamical systems and the spatial
and temporal translational symmetries in spatio-temporal systems. The symmetry
in the delay-dynamical systems can also be interpreted as a quasi-space
translational symmetry and a discrete temporal translational symmetry applying
the spatio-temporal interpretation of delay systems.

In the following, we describe the analogy between delay and spatio-temporal
dynamical systems and comment on the symmetries they incorporate. To highlight
the analogy, we visualize in Fig. 10 a) a time series of a Mackey-Glass delay
system given by $\dot{y}(t)=-y(t)+3y(t-100)/(1+y(t-100)^{10})$[1] in the
spatial representation introduced by Arecchi et al.[2]. Hereby, the y-axis
represents the length of the history interval $[0,\tau)$ (quasi-space) and the
x-axis is the time evolution $k$ in steps of the delay length $\tau$. Using
this representation, we can obtain the discrete time evolution operator of the
delay dynamical system along the x-axis mapping from the interval
$[t_{0},t_{0}+\tau)$ to the interval $[t_{0}+\tau,t_{0}+2\tau)$. The delay-
dynamical system thereby reveals quasi-space translational symmetry and a
discrete temporal translational symmetry, meaning that the time evolution
operator of the system is invariant under shift in quasi-space, temporal
shifts in the interval $[0,\tau)$ and for discrete temporal shifts in
multiples of $\tau$. Both symmetries are due to the fact that $t_{0}$ can be
set arbitrarily.

Figure 10: a) Spatial representation of a chaotic time series generated by a
Mackey Glass delay system as given in [1]. In this representation delay
intervals $[k\tau,(k+1)\tau)_{k\in\mathbb{N}}$ are plotted along the y-axis
whereas the x-axis gives the time evolution in $k$-multiples of the delay
$\tau$. As suggested in [1], in this panel every second delay interval is
plotted. b) Chaotic time series of a Kuramoto-Sivashinsky model with a spatial
size of $L=20\pi$.

In Fig. 10 b), we present a time series of the homogeneous Kuramoto-
Sivashinsky model with a size of $L=20\pi$. In analogy to the delay dynamical
system, on the y-axis, we present the spatial domain of the spatio-temporal
system. Due to its periodic boundary condition, the spatio-temporal Kuramoto-
Sivashinsky model reveals, besides spatial and temporal translational symmetry
also a spatial reﬂectional symmetry. Spatial translational symmetry is
satisfied since the time evolution operator of the KS system is invariant
under spatial shifts in the domain $x\in[x_{o},x_{o}+L)$ and accordingly,
$x_{0}$ can be set arbitrarily. For an illustration of translational
symmetries in the spatio-temporal and delay-dynamical systems, we refer to the
video files in the supplemental material. In the videos, we show the evolution
of both systems along a cylinder surface in order to emphasize commonalities
and visualize the corresponding translational symmetries. Both, the dynamics
of a constant delay system (CDS) and a homogenous spatio-temporal system
(HSTS) can be projected onto a cylinder with the spatial (HSTS) or quasi-
spatial (CDS) variable shown in the angular position and the temporal
evolution in the vertical direction. For the HSTS the cylinder is filled
during temporal evolution by horizontal rings of spatial profiles from bottom
to top, while for the CDS it is spiral, filling the cylinder from bottom to
top. The rotation of the cylinder illustrates the equivalence of all spatial
and quasi-spatial positions, respectively. For full chaotic dynamics, the
statistics of all spatial (quasi-spatial) positions are the same. And for
sufficiently large radii (system sizes), the statistics even does not change
with changing the radii (system sizes).

In conclusion, both delay-dynamical and homogeneous spatio-temporal systems
contain translational symmetries. In the manuscript, we exploit these
symmetries using scalable neural networks capable of far-reaching inference
beyond the dynamical regime they were trained on. As discussed in the
manuscript, while the symmetries allow for the scaling of the neural networks,
their prediction ability will further depend on the dynamical regime (system
size) used during training as we further discuss in detail in Appendix B.

## References

  * Giacomelli and Politi [1996] G. Giacomelli and A. Politi, Relationship between delayed and spatially extended dynamical systems, Phys. Rev. Lett. 76, 2686 (1996).
  * Arecchi  _et al._ [1992] F. Arecchi, G. Giacomelli, A. Lapucci, and R. Meucci, Two-dimensional representation of a delayed dynamical system, Phys. Rev. A 45, R4225 (1992).


∎

11institutetext: Arzu Kurt
11email<EMAIL_ADDRESS>

1 Department of Physics, Bolu Abant İzzet Baysal University, Bolu 14030,
Turkey.

# Two-time correlation functions beyond quantum regression theorem: Effect of
external noise

Arzu Kurt1

(Received: date / Accepted: date)

###### Abstract

We present the results of a study of the two-time correlation functions of a
dichotomously driven two-level system in contact with a thermal bath by using
corrections beyond quantum regression theorem. In the strong system-
environment coupling regime, it is found that the noise parameters could be
tuned to control the magnitude of corrections at low environmental
temperatures. Motional averaging and narrowing effect of the external noise
were observed on the absorption and emission spectra of the two-state system.
Furthermore, effects similar to the destruction of tunneling and noise
enhancement of transport are observed in the dynamics of the two-time
correlations.

###### Keywords:

Quantum regression theorem telegraph noise non-Markovian dynamics

##  1 Introduction

Developments in spectroscopy techniques have increased the need to develop
effective methods to compute the multi-time correlation functions which encode
large amounts of crucial information on the measured spectra Scully97 (1, 2,
3, 4, 5, 6), that might not be obtained from the expectation values of the
system operators. Use of two-time correlation functions to model measured
spectra range from the optical spectrum of the semiconductor quantum dot
Bertelot2006 (7) and optical emission in quantum-dot-cavity systems Winger2009
(8) to Mollow triplet spectra of a resonantly driven quantum dot in a
microcavity Ulrich2011 (9). Jorgensen et. al.’s work on the phonon emission
spectrum of the spin-boson modelJorgensen2020 (10) is a recent example of
modeling the spectral response of open quantum systems by using two-time
correlation functions (TTCFs) which has a long history. TTCFs have, also, been
used in inferring direct quantum dynamics in experimental investigation of the
long-lived electronic quantum coherence for a photosynthetic pigment-protein
structure known as Fenna-Matthews-Olsen (FMO) complex Brixner2005 (11).
Besides, the authors Li2013 (4) reported a study of motional averaging and
narrowing in the situation where the energy levels of qubit system are
modulated by a random telegraph noise in experimental set-up.

Many different approaches to calculate the TTCFs have been developed over the
yearsOnsager31a (12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23). Quantum
regression theorem (QRT) provides a direct relation between the dynamics of
single-time expectation values of system operators and their two-time
correlation functions in Markovian systems and is used extensively in the
earlier modelling attempts Onsager31a (12, 13). Validity of QRT in formulation
of TTCFs for non-Markovian systems was questioned by various groups Ford1996
(24, 25, 26, 27) which lead to various proposals to amend it to account for
the memory effects Alonso2005 (14, 15, 16, 17, 18, 19, 20, 21, 22, 23). de
Vega and Alonso deVega2008 (22) used the stochastic Schrödinger equation
approach to derive time evolution equations for single and two-time
correlation functions of arbitrary system operators and used it to calculate
the spectral properties of a two level atom weakly interacting with a highly
non-Markovian electromagnetic environment. This model was criticized as being
valid only for a limited class of system-environment coupling operators and
the zero temperature environment by Goan et. al. goan11 (18) who used a second
order master equation approach to systematically improve the QRT evolution
equations by including terms that account for the memory effects. They have
found considerable differences in time-evolution of the TTCFs obtained in QRT
and the improved QRT models even for the weak system-bath coupling regime. Jin
et. al. derived time-local equations of motion for the two-time non-Markovian
correlation function of arbitrary system operators in open quantum systems and
used them to demonstrate that they could account for the typical non-Markovian
effects in current fluctuation spectra of single and double quantum dots.
Jin2016 (23). Cosacchi et. al. Cosacchi2018 (28) have improved the iterative
path integral algorithm for the numerically exact computation of multi-time
correlation functions introduced by Shao and Makri Shao2002 (29) and used it
to study spectra of a two-level system simultaneously in contact with
Markovian and non-Markovian environments. Their method was found to reproduce
the non-Markovian features in the emission of a quantum dot in contact with
acoustic phonons. Recent attempts to include multi-time measurement scenario
without applying QRT, also exist Pollock2018 (30, 31, 10). Jorgensen and
Pollock suggest, for example, an alternative formulation based on time-
evolving matrix product operator(TEMPO) algorithm to both effectively
calculate multi-time correlation functions for more general systems and
indicate to enhance the efficiency of the method by directly computing the the
emission spectra in non-Markovian case for the spin-boson model beyond QRT
Jorgensen2019 (31). In order to incorporate multi-time correlation functions
Jorgensen and Pollock developed a powerful protocol based on the transfer-
tensor framework Pollock2018 (30) by examining full of steady-state
correlations for the spin-boson model without computing inhomogeneous terms.
Also, it is studied the dynamics of system-bath correlations in the spirit of
the non-Markovianity relative entropy measurement Jorgensen2020 (10).

Moreover, various authors suggested the possibility of using the violation of
QRT as an indication of non-Markovianity of the system dynamics Guarnieri2014
(32, 33, 34, 28). Guarnieri et. al. Guarnieri2014 (32) compared the validity
of QRT with CP-divisibility and distinguishability non-Markovianity measures
for the exactly solvable pure dephasing spin boson model and showed that QRT
might be violated even when the measures indicate Markovian dynamics. Chen and
Manirul Ali have used the non-Markovian TTCFs developed by Goan et. al. goan11
(18) to study the violation of Leggett-Garg inequality for a two-level system
in a non-Markovian dephasing environment and reported that the violation
decreases with increasing environmental coupling Chen2014 (33). Also, Kurt
Kurt2020 (35) has considered the violation of QRT for the spin-boson model in
the strong coupling regime by using the method developed in Ref. goan11 (18)
and found that beyond QRT corrections were small for the range of studied
parameters.

Most of the studies on the non-Markovian two-time correlation functions deals
with systems weakly interacting with a non-Markovian thermal environment or
Markovian plus another non-Markovian environment. However, some realistic
systems, such as FMO light harvesting complexes, are in contact with a thermal
environment which might be Markovian or non-Markovian as well as driven by
external disturbances which can not be described as a collection of harmonic
oscillators Jang2018 (36). Aim of the present work is to investigate the
behaviour of two-time correlation functions for a two-state system (TSS) in a
thermal bath when the TSS energy gap is modulated by an external random
telegraph signal in the strong coupling regime. Noise averaged time-
convolutionless master equation with QRT plus corrections introduced by Goan
et. al. goan11 (18) is used to examine the effect of the external noise on the
QRT violations as a function of system, environment and noise parameters.

The paper is organized as the following. We present the model in the polaron
frame and the time evolution of a single time as well as two-time correlations
functions in non-Markovian case in Sec. 2. In Sec. 3, we display the result
including the averaged over the noise realization for the TTCFs for the non-
Markovian system quantities. We briefly summarize the paper in Sec. 4.

##  2 Model

We consider a two-state system (TSS) that is driven by a random telegraph
noise (RTN), which in turn is in contact with a structured bath. The
Hamiltonian describing the closed system (TSS + bath) is:

$H=H_{S}(t)+H_{B}+H_{I},$ (1)

where $H_{S}(t)$ is time dependent system Hamiltonian, $H_{B}$ is the
Hamiltonian of the thermal bath which consists of a set of independent
harmonic oscillators, and interaction Hamiltonian between these two is
described by $H_{I}$. Using the units with $\hbar=1$, the total system
Hamiltonian is:

$H=\frac{\epsilon(t)}{2}\sigma_{z}+\frac{V}{2}\sigma_{x}+\sum_{k}\omega_{k}\,b^{\dagger}_{k}\,b_{k}+\sigma_{z}\sum_{k}g_{k}\left(b^{\dagger}_{k}+b_{k}\right),$
(2)

where $V$ is the tunneling matrix element between the two states, $\sigma_{i}$
are the Pauli spin matrices, $b^{\dagger}_{k}(b_{k})$ are creation
(annihilation) operators of the $k^{th}$ bath mode, $g_{k}$ is the interaction
strength between the system and the $k$th mode of the environmental
oscillators. Here, $\epsilon(t)=\epsilon_{0}+\Omega\,\alpha(t)$ where
$\epsilon_{0}$ is a bare energy difference between the states of the TSS and
the time-dependent term accounts for the random telegraph noise modulation of
the state energies. RTN has two states with noise amplitude $\pm\Omega$. Its
character is described by zero average ($\langle\alpha(t)\rangle=0$) and
exponentially decaying auto-correlation function
($\langle\alpha(t)\alpha(t^{\prime})\rangle=e^{-\nu|t-t^{\prime}|}$) with
noise frequency $\nu$ that is a measure of number of flippings of the noise
signal in unit time interval.

In the present study, we consider strong interaction between the TSS and the
bath and use polaron transformation to express the interaction Hamiltonian in
a form that lets one use it in perturbative master equation. In the polaron
frame, the total time-dependent Hamiltonian of Eq. (1) can be expressed as:

$\displaystyle H^{{}^{\prime}}$ $\displaystyle=$ $\displaystyle
H_{S}^{{}^{\prime}}(t)+H_{B}^{{}^{\prime}}+H_{I}^{{}^{\prime}},$ (3)
$\displaystyle=$
$\displaystyle\frac{\epsilon(t)}{2}\sigma_{z}+V_{r}\,\sigma_{x}+\sum_{k}\omega_{k}b^{\dagger}_{k}b_{k}+\sigma_{+}\,B_{-}+\sigma_{-}\,B_{+},$

where the superscript $"\,^{\prime}\,"$ means that the operator
$O^{{}^{\prime}}$ is in the polaron frame. After this point, we will drop the
superscript notation, for simplicity. $\sigma_{\pm}$ are the flip operators of
the TSS. $B_{\pm}$ are the bath correlation operators and $V_{r}$ is the
reduced tunneling matrix element, and given as a function of system-bath
coupling strength $g_{k}$ with the $k^{th}$ oscillator mode respectively:

$\displaystyle B_{\pm}$ $\displaystyle=$ $\displaystyle\langle
e^{\mp\sum_{k}\frac{2\,g_{k}}{\omega_{k}}\left(b_{k}^{\dagger}-b_{k}\right)}\rangle_{R},$
(4) $\displaystyle V_{r}$ $\displaystyle=$ $\displaystyle
V\exp\left(-\frac{1}{4\pi}\int_{0}^{\infty}\frac{J(\omega)}{\omega^{2}}\coth\left(\beta\omega\right)\right)$
(5)

where $\langle\dots\rangle_{R}$ implies averaging over the environmental
degrees of freedom. It is significant to note that the reduced tunneling rate
$V_{r}$ is equal to zero for the bath spectral densities of which frequency
exponent is less than two which is the case in the current study.

The evolution equation of the reduced density matrix of the TSS in the
interaction picture can be written as

$\frac{d}{d\,t}\rho_{S}(t)=-\int_{0}^{t}dt_{1}\,\mathrm{Tr}_{R}{\left[H_{I}(t),\left[H_{I}(t_{1}),\rho_{S}\otimes\rho_{R}\right]\right]}.$
(6)

where $\mathrm{Tr}_{R}$ indicates partial trace over the environmental modes.
In the Schrödinger picture, the evolution equations of any system operator $A$
in the second order can be given by

$\displaystyle\frac{d}{d\,t_{1}}\langle A(t_{1})\rangle$ $\displaystyle=$
$\displaystyle i\,\mathrm{Tr}_{S\otimes
R}\left(\\{\left[H_{S},A\right]\\}(t_{1})\,\rho_{T}(0)\right)$ (7)
$\displaystyle+$ $\displaystyle\int_{0}^{t_{1}}\,d\tau\,\mathrm{Tr}_{S\otimes
R}\left(\\{\tilde{H}_{I}(\tau-
t_{1})\left[A,H_{I}\right]\\}(t_{1})\rho_{T}(0)\right.$
$\displaystyle\left.+\\{\left[H_{I},A\right]\tilde{H}_{I}(\tau-
t_{1})\\}(t_{1})\rho_{T}(0)\right).$

Here,
$\tilde{H}_{I}(t)=V\left(\sigma_{-}(t)\,B_{+}(t)+\sigma_{+}(t)\,B_{-}(t)\right)$
which describes the time evolution in the interaction picture in the polaron
frame. Tr with the symbol $S\otimes R$ refers to a trace over the Hilbert
space of the total system. Curly brackets in Eqs. (7) and (8) indicate that
the expression should be evaluated in the Heisenberg picture and its time
should be taken as given in the post bracket. It is straightforward to express
the TTCFs in non-Markovian regime by using the result of Refs.goan11 (18, 21):

$\displaystyle\frac{d}{d\,t_{1}}\langle A(t_{1})B(t_{2})\rangle$
$\displaystyle=$ $\displaystyle i\,\mathrm{Tr}_{S\otimes
R}\left(\\{\left[H_{S},A\right]\\}(t_{1})B(t_{2})\,\rho_{T}(0)\right)$ (8)
$\displaystyle+$ $\displaystyle\int_{0}^{t_{1}}\,d\tau\,\mathrm{Tr}_{S\otimes
R}\left(\\{\tilde{H}_{I}(\tau-
t_{1})\left[A,H_{I}\right]\\}(t_{1})B(t_{2})\rho_{T}(0)\right.$
$\displaystyle\left.+\\{\left[H_{I},A\right]\tilde{H}_{I}(\tau-
t_{1})\\}(t_{1})B(t_{2})\rho_{T}(0)\right)$ $\displaystyle+$
$\displaystyle\int_{0}^{t_{2}}d\tau\,\mathrm{Tr}_{S\otimes
R}\left(\\{\left[H_{I},A\right]\\}(t_{1})\\{\left[B,\tilde{H}_{I}(\tau-
t_{2})\right]\\}(t_{2})\,\rho_{T}(0)\right).$

Here, the first two terms on the right hand side of Eq. (8) belong to the QRT
terms while the last term accounts for the corrections for the non-Markovian
effects. The last integral term in Eq. (8) is the source of violation of the
quantum regression theorem. By using the set of TTCFs in non-Markovian
approximation in Eq.(8) we obtain the time evolution of TTCFs of system
operators of dichotomously driven spin-boson model in non-Markovian case and
average those equations over the large number of noise trajectories in the
following section.

##  3 Results

In the present study, we will analyze the dynamics of
$\langle\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle$ for a TSS which is in
contact with a thermal bath at temperature $T$. We will consider the strong
system-bath coupling regime and use a structured spectral density for
describing the interaction. The energy levels of the TSS is assumed to be
driven stochastically with an RTN signal. By using the Hamiltonian Eq. (3) in
the master equation Eq. (8) goan11 (18) with $A=\sigma_{z}$ and
$B=\sigma_{z}$, one can obtain following set of four coupled differential
equations for various TTCFs and $\langle\sigma_{z}(t)\rangle$ which include
beyond the quantum regression theorem contributions:

$\displaystyle\frac{d}{d\,t_{1}}\langle\sigma_{z}(t_{1})\,\sigma_{z}(t_{2})\rangle$
$\displaystyle=$
$\displaystyle-\Gamma_{1}(t_{1})\langle\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle-\Gamma_{2}(t_{1})\langle\sigma_{z}(t_{2})\rangle$
(9)
$\displaystyle+4\,\Gamma_{3}(t_{1},t_{2})\langle\sigma_{-}(t_{1})\sigma_{+}(t_{2})\rangle+4\,\Gamma_{4}(t_{1},t_{2})\langle\sigma_{+}(t_{1})\sigma_{-}(t_{2})\rangle$
$\displaystyle\frac{d}{d\,t_{1}}\langle\sigma_{z}(t_{1})\rangle$
$\displaystyle=$
$\displaystyle-\Gamma_{1}(t_{1})\langle\sigma_{z}(t_{1})\rangle-\Gamma_{2}(t_{1}),$
(10)
$\displaystyle\frac{d}{d\,t_{1}}\langle\sigma_{+}(t_{1})\sigma_{-}(t_{2})\rangle$
$\displaystyle=$
$\displaystyle\left[i\,\left(\epsilon_{0}+\Omega\,\alpha(t_{1})\right)-\Gamma_{5}(t_{1})\right]\langle\sigma_{+}(t_{1})\sigma_{-}(t_{2})\rangle$
(11)
$\displaystyle+\Gamma_{3}(t_{1},t_{2})\langle\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle,$
$\displaystyle\frac{d}{d\,t_{1}}\langle\sigma_{-}(t_{1})\sigma_{+}(t_{2})\rangle$
$\displaystyle=$
$\displaystyle-\left[i\,\left(\epsilon_{0}+\Omega\,\alpha(t_{1})\right)+\Gamma_{5}(t_{1})^{*}\right]\langle\sigma_{-}(t_{1})\sigma_{+}(t_{2})\rangle$
(12)
$\displaystyle+\Gamma_{4}(t_{1},t_{2})\langle\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle,$

We should note that derivation of Eqs. (9)-(12) is simplified by the fact that
the reduced tunneling rate $V_{r}$ in Eq. (3) equals to zero for the chosen
spectral density. All the terms in Eqs. (9)-(12) that contain $\Gamma_{i}$
coefficients with $i=3$ and 4 are due to the corrections to the QRT (the last
line of Eq. (8)) while those with $i=1,2,5$ and $6$ account for the quantum
regression theorem. The environmental correlation functions $\Gamma_{i}(t)$
are derived from thermal averages of various system and bath operators as:

$\displaystyle\Gamma_{1}(t_{1})$ $\displaystyle=$ $\displaystyle
4\,V^{2}\int_{0}^{t_{1}}d\tau\,e^{-Q_{2}\left(t_{1}-\tau\right)}\cos{\left[Q_{1}(t_{1}-\tau)\right]}\,\cos{\left[f(t_{1},\tau)\right]},$
$\displaystyle\Gamma_{2}(t_{1})$ $\displaystyle=$ $\displaystyle
4\,V^{2}\int_{0}^{t_{1}}d\tau\,e^{-Q_{2}\left(t_{1}-\tau\right)}\,\sin{\left[Q_{1}(t_{1}-\tau)\right]}\,\sin{\left[f(t_{1},\tau)\right]},$
$\displaystyle\Gamma_{3}(t_{1},t_{2})$ $\displaystyle=$ $\displaystyle
V^{2}\int_{0}^{t_{2}}d\tau\,e^{-Q_{2}\left(t_{1}-\tau\right)+i\,Q_{1}\left(t_{1}-\tau\right)}\,e^{i\,f(t_{2},\tau)},$
$\displaystyle\Gamma_{4}(t_{1},t_{2})$ $\displaystyle=$ $\displaystyle
V^{2}\int_{0}^{t_{2}}d\tau\,e^{-Q_{2}\left(t_{1}-\tau\right)+i\,Q_{1}\left(t_{1}-\tau\right)}\,e^{-i\,f(t_{2},\tau)},$
$\displaystyle\Gamma_{5}(t_{1})$ $\displaystyle=$ $\displaystyle
2\,V^{2}\int_{0}^{t_{1}}d\tau\,e^{-Q_{2}\left(t_{1}-\tau\right)}\,\cos{\left[Q_{1}(t_{1}-\tau)\right]}\,e^{i\,f(t_{1},\tau)},$

where

$\displaystyle Q_{1}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2\,\pi}\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega^{2}}\sin(\omega
t),$ (13) $\displaystyle Q_{2}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2\,\pi}\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega^{2}}\coth\left(\frac{\beta\,\omega}{2}\right)(1-\cos(\omega
t)),$ (14) $\displaystyle f(t,\tau)$ $\displaystyle=$
$\displaystyle\epsilon_{0}(t-\tau)+\Omega\,\int_{t}^{\tau}d\zeta\,\alpha(\zeta),$
(15)

where $Q_{1,2}(t)$ are the imaginary and the real parts of the bath
correlation function. $f(t,\tau)$ accounts for the static as well as
stochastically fluctuating system energy gap. Spectral density $J(\omega)$
characterizes the interaction between the TSS and the environment as a
function of environmental modes. We choose $J(\omega)$ as the structured
spectral density. In this model, the environment is described by single
harmonic oscillator (with frequency $\omega_{0}$) which is in contact with an
Ohmic thermal bath that broadens its energy levels Garg1985 (37):

$J(\omega)=8\,\kappa^{2}\frac{\gamma\,\omega_{0}\,\omega}{(\omega^{2}-\omega_{0}^{2})^{2}+4\gamma^{2}\,\omega^{2}},$
(16)

where $\kappa$ is the magnitude of system-bath coupling, $\omega_{0}$ is the
frequency of the central harmonic oscillator, $\gamma$ is the broadening of
the oscillator levels due to interaction with the environment. A rough measure
of the strength of the interaction between the TSS and its environment is the
reorganization energy of the bath
($E_{r}=\int_{0}^{\infty}d\omega\,J(\omega)/\omega$). The reorganization
energy for the chosen $J(\omega)$ is equal to $\kappa^{2}/\omega_{0}$.

###  3.1 Stochastic Averaging

The RTN signal $\alpha(t)$ appears in the coupled differential equations that
describe the dynamics of the two-time correlation functions in Eqs.(9)-(12) as
bare coefficients as well as its integral in definition of $\Gamma_{i}(t)$
coefficients. One can either use an ensemble averaging approach to find the
noise averaged TTCFs by solving Eqs. (9)-(12) for an appropriately controlled
number of different realizations of RTN noise and average the solutions or use
probability density of RTN signal to average the coupled differential
equations and solve those averaged equations. Thanks to theorems due to
Bourret-Frisch Bourret73 (38) which provide a factorization for the multi-time
correlation functions of dichotomous process and Loginov-Shapiro Shapiro78
(39) which relates the average of time evolution operator of the noise with
that of the system operators, Eqs. (9)-(12) can be averaged exactly over the
RTN Magazzu2017 (40). The averaging process effectively doubles the number of
correlation functions by coupling the correlations of the noise signal
$\alpha(t)$ with the system operator correlations. We will adopt the single
average notation to describe both the correlation and noise averaging, i. e.,
$\langle\sigma_{z}(t)\rangle$ denotes noise averaged correlation of
$\sigma_{z}$ at times $t_{1}$ and $t_{2}$ in the rest of the paper. The set of
RTN averaged coupled differential equations for
$\langle\sigma_{z}(t_{1})\rangle$ and
$\langle\alpha(t_{1})\sigma_{z}(t_{2})\rangle$ are:

$\displaystyle\frac{d}{d\,t_{1}}\langle\sigma_{z}(t_{1})\rangle$
$\displaystyle=$
$\displaystyle-\Gamma_{1}(t_{1})\langle\sigma_{z}(t_{1})\rangle+\Gamma_{2}(t_{1})\langle\alpha(t_{1})\,\sigma_{z}(t_{1})\rangle-\Gamma_{3}(t_{1}),$
$\displaystyle\frac{d}{d\,t_{1}}\langle\alpha(t_{1})\sigma_{z}(t_{1})\rangle$
$\displaystyle=$
$\displaystyle-(\nu+\Gamma_{1}(t_{1}))\,\langle\alpha(t_{1})\sigma_{z}(t_{1})\rangle+\Gamma_{2}(t_{1})\langle\sigma_{z}(t_{1})\rangle-\Gamma_{4}(t_{1}).$
(17)

The noise averaged coupled differential equations for the six two-time
correlation functions can be expressed as:

$\frac{d}{dt_{1}}Y=A(t_{1},t_{2})\cdot Y+b(t_{1},t_{2})$ (18)

where $Y$ is a column vector whose elements are the different correlations:

$\displaystyle
Y=\left(\langle\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle,\langle\alpha(t_{1})\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle,\langle\sigma_{+}(t_{1})\sigma_{-}(t_{2})\rangle,\langle\alpha(t_{1})\sigma_{+}(t_{1})\sigma_{-}(t_{2})\rangle,\right.$
$\displaystyle\left.\langle\sigma_{-}(t_{1})\sigma_{+}(t_{2})\rangle,\langle\alpha(t_{1})\sigma_{-}(t_{1})\sigma_{+}(t_{2})\rangle\right)^{T}$

and $A$ is the time-dependent coefficient matrix of the noise-averaged coupled
differential equations:

$\displaystyle
A=\left(\begin{array}[]{cccccc}-\Gamma_{1,1}&\Gamma_{1,2}&-4\Gamma_{3,1}&-4\Gamma_{3,2}&4\Gamma_{4,1}&4\Gamma_{4,2}\\\
\Gamma_{1,2}&-\left(\nu+\Gamma_{1,1}\right)&4\Gamma_{3,2}&4\Gamma_{3,1}&-4\Gamma_{4,2}&4\Gamma_{4,1}\\\
\Gamma_{4,1}&-\Gamma_{4,2}&(i\epsilon_{0}-\Gamma_{5,1})&(i\Omega+\Gamma_{5,2})&0&0\\\
-\Gamma_{4,2}&\Gamma_{4,1}&(i\Omega+\Gamma_{5,2})&-(\nu-i\epsilon_{0}+\Gamma_{5,1})&0&0\\\
\Gamma_{3,1}&\Gamma_{3,2}&0&0&-(i\epsilon_{0}+\Gamma_{6,1})&-(i\Omega+\Gamma_{6,2})\\\
\Gamma_{3,2}&\Gamma_{3,1}&0&0&-(i\Omega+\Gamma_{6,2})&-(\nu+i\epsilon_{0}+\Gamma_{6,1})\end{array}\right)$

where all $\Gamma$s in $A$ depend on time $t_{1}$ except $\Gamma_{3,}$ and
$\Gamma_{4,}$ which have double time arguments $(t_{1},t_{2})$. The
$b(t_{1},t_{2})$ term in Eq. (18) contributes to the change in
$\langle\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle$ and
$\langle\alpha(t_{1})\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle$ terms and is
found to be:

$\displaystyle
b=\left(-\Gamma_{2,1}(t)g_{1}(t_{2})-e^{-\nu|t-t_{2}|}\Gamma_{2,2}(t)g_{2}(t_{2}),-\Gamma_{2,2}(t)g_{1}(t_{2})-e^{-\nu|t-t_{2}|}\Gamma_{2,1}(t)g_{2}(t_{2}),0,0,0,0\right)^{T}$

where $g_{1}(t)=\langle\sigma_{z}(t)\rangle$ and
$g_{2}(t)=\langle\alpha(t)\sigma_{z}(t)\rangle$. $\Gamma_{i,j}$s of $A$
involve the time evolution operator
$S(t,t^{\prime})=\exp\left[-i\,\Omega\int_{t^{\prime}}^{t}d\tau\alpha(\tau)\right]$
of the Kubo oscillator and are given as:

$\displaystyle\Gamma_{1,1}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}d\tau\,\mathcal{E}_{cc}(t-\tau)S_{0}(t-\tau)$
$\displaystyle\Gamma_{1,2}(t)$ $\displaystyle=$ $\displaystyle
i\int_{0}^{t}d\tau\,\mathcal{E}_{cs}(t-\tau)S_{1}(t-\tau)$
$\displaystyle\Gamma_{2,1}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}d\tau\,\mathcal{E}_{ss}(t-\tau)S_{0}(t-\tau)$
$\displaystyle\Gamma_{2,2}(t)$ $\displaystyle=$ $\displaystyle
i\int_{0}^{t}d\tau\,\mathcal{E}_{sc}(t-\tau)S_{1}(t-\tau)$
$\displaystyle\Gamma_{3,1}(t,t_{2})$ $\displaystyle=$
$\displaystyle\int_{0}^{t_{2}}d\tau\,\mathcal{E}_{f+}(t-\tau)e^{i\epsilon(t_{2}-\tau)}S_{0}(t_{2}-\tau)$
$\displaystyle\Gamma_{3,2}(t,t_{2})$ $\displaystyle=$
$\displaystyle\int_{0}^{t_{2}}d\tau\,\mathcal{E}_{f+}(t-\tau)e^{i\epsilon(t_{2}-\tau)}S_{1}(t_{2}-\tau)$
(20) $\displaystyle\Gamma_{4,1}(t,t_{2})$ $\displaystyle=$
$\displaystyle\int_{0}^{t_{2}}d\tau\,\mathcal{E}_{f-}(t-\tau)e^{-i\epsilon(t_{2}-\tau)}S_{0}(t_{2}-\tau)$
$\displaystyle\Gamma_{4,2}(t,t_{2})$ $\displaystyle=$
$\displaystyle\int_{0}^{t_{2}}d\tau\,\mathcal{E}_{f-}(t-\tau)e^{-i\epsilon(t_{2}-\tau)}S_{1}(t_{2}-\tau)$
$\displaystyle\Gamma_{5,1}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}d\tau\,\mathcal{E}_{c-}(t-\tau)S_{0}(t-\tau)$
$\displaystyle\Gamma_{5,2}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}d\tau\,\mathcal{E}_{c-}(t-\tau)S_{1}(t-\tau)$
$\displaystyle\Gamma_{6,1}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}d\tau\,\mathcal{E}_{c+}(t-\tau)S_{0}(t-\tau)$
$\displaystyle\Gamma_{6,2}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}d\tau\,\mathcal{E}_{c+}(t-\tau)S_{1}(t-\tau)$

where
$\mathcal{E}_{cs}(t)=e^{-Q_{2}(t)}\mathrm{(c)os}\left(Q_{1}(t)\right)\,\mathrm{(s)in}\left(\epsilon_{0}\,t\right)$,
$\mathcal{E}_{c\pm}(t)=e^{-Q_{2}(t)}\mathrm{cos}\left(Q_{1}(t)\right)\,e^{\pm\,i\,\epsilon_{0}\,t}$,
and $\mathcal{E}_{f\pm}(t)=e^{-Q_{2}(t)}e^{i(Q_{1}(t)\pm\epsilon_{0}\,t)}$.
$S_{0}(t)$ and $S_{1}(t)$ in kernel definitions are, also, dichotomous noise
propagators Goychuk95 (41) which can be expressed as:

$\displaystyle S_{0}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2\eta}\left(\nu_{+}e^{-t\,\nu_{-}/2}-\nu_{-}e^{-t\,\nu_{+}/2}\right),$
$\displaystyle S_{1}(t)$ $\displaystyle=$ $\displaystyle
i\frac{\Omega}{\nu}\left(e^{-t\,\nu_{+}/2}-e^{-t\,\nu_{-}/2}\right)$ (21)

where $\eta=\sqrt{\nu^{2}-4\Omega^{2}}$, $\nu_{+}=\nu+\eta$, and
$\nu_{-}=\nu-\eta$.

One should note that the $\Gamma_{i,j}$ that depends on a single time argument
stem from the non-Markovian quantum regression theorem, while the ones with
double time argument arise as corrections to the QRT. They will be referred to
as QRT$+$ in the rest of the current paper. Also, when there is no external
noise, $S_{0}(t)=1$ and $S_{1}(t)=0$, so all the $\Gamma_{i,2}$ are zero.

---
(a) $\epsilon_{0}=1,\,\beta=0.02,\,\nu=0.01$

---
(b) $\epsilon_{0}=0,\,\beta=50,\,\nu=0.01$

---
(c) $\epsilon_{0}=1,\,\beta=0.02,\,\nu=1$

---
(d) $\epsilon_{0}=0,\,\beta=50,\,\nu=1$

Figure 1: Dynamics of the real part of two-time correlations
$\langle\sigma_{z}(t)\sigma_{z}(t_{2})\rangle$,
$\langle\sigma_{+}(t)\sigma_{-}(t_{2})\rangle$ and
$\langle\sigma_{-}(t)\sigma_{+}(t_{2})\rangle$ for the QRT (dashed line) and
QRT+ (solid line) approaches for various combinations of the TSS bias, the
noise frequency and temperature of the environment. Inset shows the power
spectra of $\langle\sigma_{+}(t)\sigma_{-}(t_{2})\rangle$ (solid line) and
$\langle\sigma_{+}(-)\sigma_{+}(t_{2})\rangle$ (dashed line) correlations
calculated with the QRT+ corrections. The noise amplitude $\Omega=0.75$ for
all the figures.

For the present study, we consider the strong system-environment coupling
regime with $\omega_{0}=1$, $\kappa=2$ along with a relatively large
dispersion $\gamma=0.5$ and investigate the low ($\beta=50$) and the high
temperature ($\beta=0.02$) as well as the biased ($\epsilon_{0}=1$) and the
non-biased ($\epsilon_{0}=1$) cases. With these coupling parameters, the short
time approximationGarg1985 (37) can be used to write the bath correlation
functions in a simpler form as $Q_{1}(t)=E_{r}t$ and $Q_{2}(t)=-\xi t^{2}$
where

$\xi=E_{r}/\beta+\kappa^{2}/\left(\pi\sqrt{\omega_{0}^{2}-\gamma^{2}}\right)\omega_{0}\,\mathrm{Im}\left[\Psi\left(0,1+\beta\left(\gamma+i\sqrt{\omega_{0}^{2}-\gamma^{2}}\right)/(2\pi)\right)\right].$

For this choice of the bath correlation functions, the $\Gamma_{i,j}$
coefficients in Eq. (3.1) can be expressed as sum of error functions of
various sum and difference combinations of the TSS bias, noise frequency and
the bath reorganization energy. But, the expressions are too cumbersome and
will not be displayed here.

We start the discussion of the results with presentation of the various
qualitatively different forms of time-dependence of the two-time correlation
functions calculated in non-Markovian QRT and QRT+ approximations in Fig.
1a-d. While $\langle\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle$ is found to
decay exponentially with a rate that depends on the temperature, the bias and
the properties of the noise for the strong system-bath coupling,
$\langle\sigma_{+}(t_{1})\sigma_{-}(t_{2})\rangle$ and
$\langle\sigma_{-}(t_{1})\sigma_{+}(t_{2})\rangle$ display richer dynamics. It
can be seen from the figures that the difference between the QRT and QRT+
approximations for all three correlation functions at low temperature is
higher than the difference at high temperature. The difference between the two
will be discussed in more detail below.

---
(a) $\beta=50$

---
(b) $\beta=0.02$

Figure 2: Noise amplitude and color dependence of absorption spectra at low
and high temperature for the non-biased TSS.

The dynamics of a two state system driven by RTN can be deduced intuitively at
the asymptotic limits with the help of the concept of the noise color or the
Kubo number $K=\Omega/\nu$. As the noise jumps of magnitude $\Delta
E=\hbar\Omega$ occur with average time spacing $\Delta t=1/\nu$, the energy-
time uncertainty relation ($\Delta E\Delta t=\hbar\Omega/\nu\geq\hbar/2$)
would make resolving those fluctuations impossible at the fast jumping limit.
Resulting motional averaging leads to emission or absorption at dynamically
averaged $\epsilon_{0}$ frequency although the system is at one of the states
$\epsilon_{0}\pm\Omega$ at all times. At the slow noise limit, the mean time
between the jump events is long and the dynamics of the system can be
considered as the static average at frequencies $\epsilon_{0}\pm\Omega$.
Fourier transform of $\langle\sigma_{-}(t+\tau)\sigma_{+}(t)\rangle$ and
$\langle\sigma_{+}(t+\tau)\sigma_{-}(t)\rangle$ are measures of absorption and
emission spectrum, respectively. The insets in Fig. 1a-d display the power
spectrum of the correlations
$\langle\sigma_{+}(t_{1})\sigma_{-}(t_{2})\rangle$ and
$\langle\sigma_{-}(t_{1})\sigma_{+}(t_{2})\rangle$. The first row where $K=75$
is the slow noise (strongly colored noise) limit and one can observe emission
and absorption phenomena at frequencies $\epsilon_{0}\pm\Omega$ clearly. The
plots in the second row display the dynamics of the correlation at
intermediate noise color with $K=0.75$ for which the motional averaging starts
to dominate the power spectrum as the displayed spectra tend to broaden and
become a single peak at frequency $\epsilon_{0}$. We display the $\Omega$, $K$
and frequency dependence of the absorption spectrum at high and low
temperatures in Fig.2. As can be seen from the plots, the single peak to
double peak behaviour start emerging around $K\approx 0.7$ independent of the
bath temperature.

---
(a) $\langle\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle$

---
(b) $\langle\sigma_{-}(t_{1})\sigma_{+}(t_{2})\rangle$

Figure 3: Noise frequency $\nu$ and amplitude $\Omega$ dependence of the
relative difference between the QRT and QRT+ ($\epsilon_{0}=0,\,\beta=50$).
Thick line indicates $\nu=2\Omega$.

There has been some discussion on the relevance of QRT$+$ terms on the non-
Markovianity of the dynamics of the system, which suggests that the difference
can be used as a measure of non-Markovianity Ali2015 (34, 28, 32, 33). As a
rough measure of the difference between the two dynamics, we adopt a heuristic
approach and define

$\Delta_{c}=100\left|\frac{I_{\textrm{QRT}}(c)-I_{\textrm{QRT+}}(c)}{I_{\textrm{QRT}}(c)}\right|,$
(22)

where $I_{QRT}(c)$ is the integral of the absolute value of the correlation
function $c$ obtained as solution of the QRT formulation. $\Delta_{c}$ is
found to be highly sensitive to the environmental temperature for the system
and the other environment parameters considered in the present study. It is
found to be close to zero ($<0.25$) at high temperature ($\beta=0.02$) for all
three correlation functions. Figure 3 displays $\Delta_{zz}$ and $\Delta_{-+}$
as function of the noise frequency and amplitude at low temperature
($\beta=50$) for the non-biased system. $\Delta_{zz}$ and $\Delta_{-+}$ are
found to show similar qualitative and quantitative dependence on noise
parameters. Interestingly, for both correlations, the difference display a
death and revival cycle as function of the noise amplitude. The cycle has a
weak dependence on the noise frequency which is peculiar because one would
expect effect of noise to be dependent on its color or Kubo number.

---
(a) $\epsilon_{0}=0,\,\beta=50$

---
(b) $\epsilon_{0}=1,\,\beta=50$

---
(c) $\epsilon_{0}=0,\,\beta=0.02$

---
(d) $\epsilon_{0}=1,\,\beta=0.02$

Figure 4: Noise frequency $\nu$ and amplitude $\Omega$ dependence of the
inverse life time of the $\langle\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle$.

The two-time correlation function
$\langle\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle$ decays exponentially for
the system and environment parameters chosen for the present study as can be
seen from Figs. 1a-d. We fit the time dependence of
$\langle\sigma_{z}(t_{1})\sigma_{z}(t_{2})\rangle$ to a function of the form
$p+(1-p)e^{-kt}$ where $p$ is the limit of
$\langle\sigma_{z}(t_{1})\sigma_{z}(1)\rangle$ as $t_{1}\to\infty$ and $k$ is
the decay parameter of the correlation. We display the noise frequency and
amplitude dependence of this exponential decay parameter $k$ at different
temperature and TSS bias parameters in Fig. 4 a-d. It can be seen from the
sub-figures that $k$ decreases towards zero at high noise amplitude
irrespective of the bath temperature, TSS bias and the noise frequency which
is reminiscent of the noise induced destruction of tunneling Berry95 (42).
Effect of $\nu$ and $\Omega$ on $k$ is found to be observably different at low
and high temperatures. At high temperature, $k$ is almost independent of the
bias and the noise frequency and it decreases with increasing noise amplitude
monotonously while one can observe an enhancement of the transport especially
for the low-frequency noise for both the biased and non-biased TSS at low
temperature. This finding is another example of so called ENAQT Rebentrost2009
(43) which describes noise assisted enhancement of quantum transport.

---
(a) $\epsilon_{0}=0,\,\beta=50$

---
(b) $\epsilon_{0}=1,\,\beta=50$

---
(c) $\epsilon_{0}=0,\,\beta=0.02$

---
(d) $\epsilon_{0}=1,\,\beta=0.02$

Figure 5: Noise frequency $\nu$ and amplitude $\Omega$ dependence of the
decoherence rate for $\langle\sigma_{-}(t_{1})\sigma_{+}(t_{2})\rangle$.

We found that the real part of the correlation
$\langle\sigma_{-}(t_{1})\sigma_{+}(t_{2})\rangle$ could be fitted to a
function of the form $e^{-\Lambda
t}\left(a_{1}\cos{\left(\omega_{1}t\right)}+a_{2}\cos{\left(\omega_{2}t\right)}\right)$
which could in return be used to extract the exponential decay rate $\Lambda$
which is similar to the decoherence rate of the $\langle\sigma_{x}\rangle$.
Figure 5 indicates the noise amplitude and frequency dependence of $\Lambda$
at high and low temperatures for the biased and non-biased TSS. The $\Lambda$
value at low $\Omega$ might be considered as the $1/T_{2}$ Wang2005 (44) of
the system. As expected, its value at high temperature is higher than at low
temperature. Decoherence rate $\Lambda$ is found to be high at large $\nu$ and
$\Omega$ values irrespective of the environmental temperature and the bias of
the TSS which is expected because dynamics is dominated by the noise. An
interesting decrease of decoherence rate can be observed at low temperatures
for the slow noise for both the biased and non-biased TSS (Fig. 5a and b).
Also, for the very weak noise strength, the decoherence rate displays a
resonance structure as function of the noise frequency as can be seen in all
four sub-plots in Fig. 5.

##  4 Conclusions

We have investigated the effect of random telegraph noise on the two-time
correlation functions of a strongly coupled spin-boson model with structured
spectral density. Both the quantum regression theorem and corrections to the
QRT terms were included in the formulation of the master equations for the
correlations. One of the important findings of the study is that the QRT+
corrections are found to be qualitatively insignificant and quantitatively
small. As the temperature of the bath or the strength of the noise is
increased the difference tends to zero with some structure in the noise
dependence at low temperature which is an indication of noise induced non-
Markovianity. The motional averaging and narrowing caused by the random
telegraph noise driving of the TSS was observed in the emission and absorption
spectra of the system as function of the noise color with a transition region
around $K=0.7$. The decoherence rate of the spectra were found to strongly
depend on the noise parameters as they are expected to be directly
proportional to the noise frequency $\nu$ and the noise strength
$\Omega^{2}/\nu$ in the slow and fast jumping limits. We have, also, studied
the decay rate of the two-time correlation of $\sigma_{z}$ which can be an
indicator of the transport and found that slow noise can enhance transport at
low temperature. The reported findings of this study might be beneficial in
the open quantum systems when the physical properties corresponding to the
two-time correlation functions are of interests and the system-environment
interaction might not be described by using a single thermal bath.

##  5 Acknowledgments

The author would like to acknowledge many useful comments and discussions with
Prof. Dr. Resul Eryiğit.

## References

  * (1) M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * (2) H. J. Carmichael,  _Statistical Methods in Quantum Optics_ (Springer Press, 1999).
  * (3) C. W. Gardiner and P. Zoller,  _Quantum Noise_ (Springer-Verlag, 2000).
  * (4) J. Li, M.P. Silveri, K.S. Kumar, J.-M Pirkkalainen, A. Vepsäläinen, W.C. Chien, J. Tuorila, M.A. Sillanpää, P.J. Hakonen, E.V. Thuneberg, and G.S. Paraoanu, “Motional averaging in a superconducting qubit,” Nat. commun. 4, 1420 (2013).
  * (5) L. Li, M. J.W. Hall, and H. M. Wisemanv, “Concepts of quantum non-markovianity: A hierarchy,” Phy. Rep. 759, 1–51 (2018).
  * (6) I. de Vega and D. Alonso, “Dynamics of non-markovian open quantum systems,” Rev. Mod. Phys. 89, 015001 (2017).
  * (7) A. Berthelot, I. Favero, G. Cassabois, C. Voisin, C. Delalande, Ph. Roussignol, R. Ferreira, and J. M. Gérard, “Unconventional motional narrowing in the optical spectrum of a semiconductor quantum dot,” Nat. Phys. 2, 759–764 (2006).
  * (8) M. Winger, T. Volz, G. Tarel, S. Portolan, A. Badolato, K. J. Hennessy, E. L. Hu, A. Beveratos, J. Finley, V. Savona, and A. Imamoğlu, “Explanation of photon correlations in the far-off-resonance optical emission from a quantum-dot–cavity system,” Phys. Rev. Lett. 103, 207403 (2009).
  * (9) S. M. Ulrich, S. Ates, S. Reitzenstein, A. Löffler, A. Forchel, and P. Michler, “Dephasing of triplet-sideband optical emission of a resonantly driven $\mathrm{InAs}/\mathrm{GaAs}$ quantum dot inside a microcavity,” Phys. Rev. Lett. 106, 247402 (2011).
  * (10) Mathias R. Jørgensen and Felix A. Pollock, “A discrete memory-kernel for multi-time correlations in non-markovian quantum processes,” Phys. Rev. A 102, 052206 (2020).
  * (11) T. Brixner, J. Stenger, H. M. Vaswani, M. Cho, R. E. Blankenship, and G. R. Fleming, “Two-dimensional spectroscopy of electronic couplings in photosynthesis,” Nature 434, 625–628 (2005).
  * (12) L. Onsager, “Reciprocal relations in irreversible processes. i.” Phys. Rev. 37, 405 (1931).
  * (13) M. Lax, “Formal theory of quantum fluctuations from a driven state,” Phys. Rev. 129, 2342–2348 (1963).
  * (14) D. Alonso and I. de Vega, “Multiple-time correlation functions for non-markovian interaction: Beyond the quantum regression theorem,” Phys. Rev. Lett. 94, 200403 (2005).
  * (15) I. de Vega and D. Alonso, “Non-markovian reduced propagator, multiple-time correlation functions, and master equations with general initial conditions in the weak-coupling limit,” Phys. Rev. A 73, 022102 (2006).
  * (16) D. Alonso and I. de Vega, “Hierarchy of equations of multiple-time correlation functions,” Phys. Rev. A 75, 052108 (2007).
  * (17) A. A. Budini, “Operator correlations and quantum regression theorem in non-markovian lindblad rate equations,” J. Stat. Phys. 131, 51–78 (2008).
  * (18) H.-S. Goan, P.-W. Chen, and C.-C. Jian, “Non-markovian finite-temperature two-time correlation functions of system operators: Beyond the quantum regression theorem,” J. Chem. Phys. 134, 124112 (2011).
  * (19) P.-W. Chen, C.-C. Jian, and H.-S. Goan, “Non-markovian dynamics of a nanomechanical resonator measured by a quantum point contact,” Phys. Rev. B 83, 115439 (2011).
  * (20) H.-S. Goan, C.-C. Jian, and P.-W. Chen, “Non-markovian finite-temperature two-time correlation functions of system operators of a pure-dephasing model,” Phys. Rev. A 82, 012111 (2010).
  * (21) D. P. S. McCutcheon, “Optical signatures of non-markovian behavior in open quantum systems,” Phys. Rev. A 93, 022119 (2016).
  * (22) I. de Vega and D. Alonso, “Emission spectra of atoms with non-markovian interaction: Fluorescence in a photonic crystal,” Phys. Rev. A 77, 043836 (2008).
  * (23) J. Jin, C. Karlewski, and M. Marthaler, “Non-markovian correlation functions for open quantum systems,” New J. Phys. 18, 083038 (2016).
  * (24) G. W. Ford and R. F. O’Connell, “There is no quantum regression theorem,” Phys. Rev. Lett. 77, 798–801 (1996).
  * (25) G. W. Ford and R. F. O’Connell, “Calculation of correlation functions in the weak coupling approximation,” Ann. Phys. 276, 144–151 (1999).
  * (26) G. W. Ford and R. F. O’Connell, “Driven systems and the lax formula,” Opt. Commun. 179, 451–461 (2000).
  * (27) Ya.M. Blanter and M. Büttiker, “Shot noise in mesoscopic conductors,” Phys. Rep. 336, 1–166 (2000).
  * (28) M. Cosacchi, M. Cygorek, F. Ungar, A. M. Barth, A. Vagov, and V. M. Axt, “Path-integral approach for nonequilibrium multitime correlation functions of open quantum systems coupled to markovian and non-markovian environments,” Phys. Rev. B 98, 125302 (2018).
  * (29) J. Shao and N. Makri, “Iterative path integral formulation of equilibrium correlation functions for quantum dissipative systems,” J. Chem. Phys. 116, 507 (2002).
  * (30) F. Pollock, C. A. Rodriguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, “Non-markovian quantum processes: Complete framework and efficient characterization,” Phys. Rev. A 97, 012127 (2018).
  * (31) M. R. Jørgensen and F. A. Pollock, “Exploiting the causal tensor network structure of quantum processes to efficiently simulate non-markovian path integrals,” Phys. Rev. Lett. 123, 240602 (2019).
  * (32) G. Guarnieri, A. Smirne, and B. Vacchini, “Quantum regression theorem and non-markovianity of quantum dynamics,” Phys. Rev. A 90, 022110 (2014).
  * (33) P.-W. Chen and Md. M. Ali, “Investigating leggett-garg inequality for a two level system under decoherence in a non-markovian dephasing environment,” Sci. Rep. 4, 6165 (2014).
  * (34) Md. M. Ali, P.-Y. L., M. W.-Y. Tu, and W.-M. Zhang, “Non-markovianity measure using two-time correlation functions,” Phys. Rev. A 92, 062306 (2015).
  * (35) A. Kurt, “Non-Markovian Corrections to Quantum Regression Theorem for the Strong Coupling Spin-Boson Model,” SAUJS 24, 596–604 (2020).
  * (36) S. J. Jang, M. Benedetta, “Delocalized excitons in natural light-harvesting complexes,” Rev. Mod. Phys. 90, 035003 (2018).
  * (37) A. Garg, J. N. Onuchic, and V. J. Ambegaokar, “Effect of friction on electron transfer in biomolecules,” J. Chem. Phys. 83, 4491–4503 (1985).
  * (38) R. C. Bourret, U. Frisch, and A. Pouquet, “Browian motion of harmonic oscillator with stochastic frequency,” Physica 65, 303–320 (1973).
  * (39) V. E. Shapiro and V. M. Loginov, “”formulae of differentiation” and their use for solving stochastic equations,” Physica A 91, 563 (1978).
  * (40) L. Magazzù, P. Hänggi, B. Spagnolo, and D. Valenti, “Quantum resonant activation,” Phys. Rev. E 95, 042104 (2017).
  * (41) I. A. Goychuk, E. G. Petrov, and V. May, “Dynamics of the dissipative two-level system driven by external telegraph noise,” Phys. Rev. E 52, 2392–2400 (1995).
  * (42) M. Berry, “Two-state quantum asymptotics,” Ann. N. Y. Acad. Sci. 755, 303 (1995).
  * (43) P. Rebentrost, M. Mohseni, I. Kassal, S. Lloyd, and A. Aspuru-Guzik, “Environment-assisted quantum transport,” New J. Phys. 11, 033003 (2009).
  * (44) X. R. Wang, Y. S. Zheng, and Sun Yin, “Spin relaxation and decoherence of two-level systems,” Phys. Rev. B 72, 121303 (2005).