11institutetext: Potsdam Institute for Climate Impact Research (PIK), Member
of the Leibniz Association, 14473 Potsdam, Germany 22institutetext: University
of Potsdam, Institute of Geoscience, Karl-Liebknecht-Straße 32, 14476 Potsdam,
Germany

# Recurrence flow measure of nonlinear dependence

Tobias Braun 11<EMAIL_ADDRESS>K. Hauke Kraemer 11 Norbert Marwan 1122

###### Abstract

Couplings in complex real-world systems are often nonlinear and scale-
dependent. In many cases, it is crucial to consider a multitude of interlinked
variables and the strengths of their correlations to adequately fathom the
dynamics of a high-dimensional nonlinear system. We propose a recurrence based
dependence measure that quantifies the relationship between multiple time
series based on the predictability of their joint evolution. The statistical
analysis of recurrence plots (RPs) is a powerful framework in nonlinear time
series analysis that has proven to be effective in addressing many fundamental
problems, e.g., regime shift detection and identification of couplings. The
recurrence flow through an RP exploits artifacts in the formation of diagonal
lines, a structure in RPs that reflects periods of predictable dynamics. By
using time-delayed variables of a deterministic uni-/multivariate system,
lagged dependencies with potentially many time scales can be captured by the
recurrence flow measure. Given an RP, no parameters are required for its
computation. We showcase the scope of the method for quantifying lagged
nonlinear correlations and put a focus on the delay selection problem in time-
delay embedding which is often used for attractor reconstruction. The
recurrence flow measure of dependence helps to identify non-uniform delays and
appears as a promising foundation for a recurrence based state space
reconstruction algorithm.

##  1 Introduction

Measures of statistical dependence represent one of the cornerstones in the
analysis of empirical data. The study of time series measured from complex
real-world systems poses a broad variety of challenges in quantifying uni- and
multivariate data sets, e.g., lagged dependencies, non-stationarity, noise
contamination, uncertainties and the limited length of time series. The set of
tools to detect and quantify statistical dependencies ranges from standard
correlation analysis techniques Kwapień et al. (2015) over graph theoretical
approaches, such as complex networks Newman (2018); Feldhoff et al. (2012),
towards causal discovery algorithms Runge et al. (2015, 2019); Ramos et al.
(2017), and detecting critical transitions Dakos et al. (2008); Boettner et
al. (2021); Boers (2018). A notoriously challenging problem is to adequately
quantify non-monotonous and nonlinear relationships in stochastic and
deterministic systems. In this context, information theoretic measures have
established as an effective framework Balasis et al. (2013); Sun et al.
(2014); Porfiri et al. (2019). However, popular methods, such as the mutual
information (MI), are not designed to treat higher-dimensional data
appropriately (even though extensions have been suggested Pompe and Runge
(2011); Kraskov et al. (2004)). The nature and strength of links in complex
systems additionally often exhibits scale-dependence, i.e., multi-scale
behaviour Nawroth and Peinke (2006). This motivates the need of methods that
are capable of unravelling dependencies at a broad range of scales, e.g.,
wavelet based methods Maraun and Kurths (2004); Agarwal et al. (2018); Braun
et al. (2021). Few methods succeed to combine both capabilities of capturing
nonlinear dependencies at multiple time scales Vlachos and Kugiumtzis (2010);
Pecora et al. (2007).

A powerful method that captures both nonlinear and multi-scale properties of a
high-dimensional dynamical system is the recurrence plot (RP) Eckmann et al.
(1987). An RP is a mathematically simple yet effective tool that encodes the
tendency of a time series to recur to formerly visited states Marwan et al.
(2007). An RP is based on a binary recurrence matrix in which recurrences are
marked by value one, giving rise to intriguing and well-interpretable
structures in the RP. Various quantification measures can be applied to a
recurrence matrix and prove powerful in classifying differing systems Thiel et
al. (2004); Schinkel et al. (2008); Klimaszewska and Żebrowski (2009); Hirata
and Aihara (2010), identifying dynamical regime transitions Marwan et al.
(2002, 2021), and detecting non-linear correlations as well as synchronization
Ramos et al. (2017); Romano et al. (2007); Nkomidio et al. (2022). Recurrence
quantification analysis (RQA) based on diagonal lines in the RP not only
allows identification of periodic behaviour Zbilut and Marwan (2008); Braun et
al. (2021), but also helps to identify unstable periodic orbits in high-
dimensional chaotic systems Bradley and Mantilla (2002). The conceptual
simplicity of RPs allows for a broad range of real-world applications, also
for challenging data that is event-like or unevenly sampled in time Banerjee
et al. (2021); Ozdes and Eroglu (2022). Recurrence measures of dependence have
not only facilitated the study of synchronization in dynamical systems Romano
et al. (2004, 2005); Senthilkumar et al. (2008); Nkomidio et al. (2022) but
have also been extended to account for lagged and conditional dependencies
Romano et al. (2007); Zou et al. (2011); Goswami et al. (2013); Ramos et al.
(2017). Further concepts, including symbolic analyses of relationships, have
been conceived more recently Porfiri and Ruiz Marín (2019). Recurrence based
quantification of statistical dependencies, thus, bares high potential to meet
the combination of above mentioned challenges. Here, we propose a novel
recurrence based measure of dependence that uses delay coordinates from a
given observational time series. Since the measure is based on RPs, nonlinear
dependencies with multiple time lags can be quantified which makes the measure
applicable to the problem of non-uniform delay selection Vlachos and
Kugiumtzis (2010); Pecora et al. (2007). The proposed dependence measure,
thus, contributes to the challenge of characterizing complex real-world
interactions using RPs.

This work is structured as follows: in Sect. 2, we introduce the recurrence
flow as a measure of dependence along with a brief review of the RP method. We
showcase its scope in different numerical experiments in Sect. 3, covering the
characterization of lagged nonlinear dependence and delay selection for
uniform and non-uniform TDE. We conclude our findings in Sect. 4.

##  2 Recurrence Flow

We are interested in characterizing nonlinear dependencies in a deterministic,
high-dimensional system that is represented by $M$ observational time series
$\left\\{s_{n}(t)\,|\,n=1,\dots,M\right\\}$. In general, the relationships
between the different time series $s_{n}(t)$ and their coordinates do not need
to be instantaneous but are often associated with time delays
$\tau_{1},\tau_{2},\dots,\tau_{m}$. Consequently, we define the recurrence
flow measure of redundance to capture such lagged dependencies.

The key idea of the proposed measure is based on the existence of diagonal
lines in RPs. An RP is a two-dimensional matrix that encodes how a system
recurs to formerly visited states $\vec{v}_{i},\,i=1,\dots,N$. In general,
this representation can be computed for systems of any dimension $d$ and is
formally given by

$\displaystyle
R_{i,j}(\tau)\,=\,\Theta\left(\varepsilon-\|\vec{v}_{i}(\tau)-\vec{v}_{j}(\tau)\|\right)\
=\
\left\\{\begin{array}[]{lll}1&\text{if}&\|\vec{v}_{i}(\tau)-\vec{v}_{j}(\tau)\|\leq\varepsilon\\\
0&\text{if}&\|\vec{v}_{i}(\tau)-\vec{v}_{j}(\tau)\|>\varepsilon\\\
\end{array}\right..$ (3)

with two arbitrary times $i$ and $j$, the vicinity threshold $\varepsilon$ and
a suitable norm $\|\cdot\|$. The states denoted by $\vec{v}$ are either given
by the available components (state variables) of the system or, in case of
only limited access to the state variables, by delayed copies of the one (or
multiple) observational time series $s_{n}(t)$ of the studied system. In
particular, $\vec{v}$ is then obtained by stacking these copies on top of each
other as it is common practise in time delay embedding (TDE).

Diagonal lines of length $l_{d}$ in an RP resemble periods of enhanced
predictability as two trajectory segments at times $i$ and $j$ evolve in
parallel in an $\varepsilon\,$-tube for $l_{d}$ time instances. For a given
system, this may reveal time periods of continuously high determinism or
uncover abrupt regime shifts Westerhold et al. (2020). Properties of the
diagonal line length distribution of an RP are linked to dynamical invariants
of paradigmatic dynamical systems Thiel et al. (2003). However, spurious
artifacts are known to disrupt, lengthen or thicken diagonal lines due to
erroneous computation of an RP Marwan (2011). An inadequate choice of the
vicinity threshold $\varepsilon$ will disrupt diagonal lines, thus,
underestimating the system’s predictability. On the other hand, too high
values will artificially merge fundamentally distinct regions of phase space.
Sampling can also alter diagonal line structures: if the system is
undersampled, diagonal lines might not emerge continuously as deterministic
time intervals are not resolved sufficiently. On the other hand, oversampling
results in artificially thickened diagonal lines (tangential motion) Kraemer
and Marwan (2019). Erroneous time delay embedding of uni-/multivariate time
series can have several undesired effects on the formation of diagonal lines;
if the embedding dimension is chosen too high, diagonal lines are artificially
lengthened and can even emerge in absence of determinism for an uncorrelated
stochastic process due to correlations in the underlying distance matrix Thiel
et al. (2006). A non-optimal choice for the embedding delay will result in
diagonal lines that are perpendicular to the line of identity ($i=j$, LOI).
This indicates inclusion of erroneous time scales Marwan (2011). The formation
of perpendicular lines is caused by the ambiguity in the reconstruction that
is introduced by not eliminating the full serial dependence; this results in
close evolution of states both forward and backward in time:
$\|\vec{v}_{i}-\vec{v}_{j}\|<\varepsilon,\,\|\vec{v}_{i+1}-\vec{v}_{j-1}\|<\varepsilon$,
i.e., the trajectory segments closely evolve in parallel, but with opposite
time directions. On top of that, additional deformations to a diagonal line
can occur, e.g., in form of bowed diagonal lines indicating that the evolution
of states at different time intervals is similar but occurs with different
velocity or temporal resolution Marwan and Kurths (2005).

We utilize the formation of such diagonal line artifacts (DLA) to identify
time scales of the system that result in well-expressed diagonal lines. The
proposed method is, thus, based on the assumption that the studied system
exhibits (at least to some degree) deterministic dynamics which will result in
meaningful diagonal lines in an RP. An effective way of retrieving information
on the formation of DLA is given by scanning an RP diagonal-wise. We can
identify an index for each diagonal at which the first recurrence pixel is
located. It appears intuitive to regard these pixels as ‘obstacles’ to an
imaginary fluid that flows along each diagonal into the RP and is not allowed
to turn (Fig. 1C/E). The formation of DLA blocks the flow. As a basic example,
we consider a noisy sinusoidal time series (Fig. 1A) with $n=5,000$ and a
period of $T=100$. Formation of perpendicular diagonal lines for $\tau=nT$
(Fig. 1B) reduces the flow through the RP compared to
$\tau=\left.T\middle/4\right.$ (Fig. 1D). We use the symmetry of the RP by
only flooding the upper triangular matrix to save computation time. The
recurrence flow $\Phi(\tau)$ can be computed for varying delays $\tau$ and
encodes similar information as an inverse autocorrelation function, yielding a
continuous representation of the redundancy between the time series and its
delayed version (Fig. 1F). For continuous variations of $\tau$ from $\tau=0$
to $\tau=\left.T\middle/4\right.$, the perpendicular diagonal lines are
progressively eliminated. This reproduces the well-known result that a
sinusoidal signal needs to be shifted by (odd multiples of) a quarter of its
period against itself to minimize redundancy.

To quantify the flow through the RP, we define the recurrence flow matrix
$\mathbf{\phi}$ (Fig. 1C/E)

$\displaystyle\phi_{i,j}(\tau)\,=\,\Bigl{(}1-\Theta\bigl{(}\varepsilon-\|\vec{v}_{i}(\tau)-\vec{v}_{j}(\tau)\|\bigr{)}\Bigr{)}\,\Theta\left(\ell_{j}-i\right),\qquad
i,\,j=1,\,\dots,\,N$ (4)

with the length $\ell_{j}$ of the $j^{\text{th}}$ flooded diagonal, i.e., the
number of subsequent zeros up to the first one. $\mathbf{\phi}_{i,j}$ depends
on the vicinity threshold $\varepsilon$, i.e., the fraction of recurrences. We
fix $\varepsilon$ at some reasonable value that corresponds to a fixed
recurrence rate (RR).

We study the dependence on the time delay $\tau$ contained in the vector
$\vec{v}(\tau)$ similar as it is done in TDE where the delays
$\tau_{1},\tau_{2},\dots,\tau_{m}$ for the different coordinates are free
parameters and need to be chosen with respect to some notion of optimality
Kraemer et al. (2022). The recurrence flow $\Phi(\tau)$ is computed by summing
over the recurrence flow matrix $\mathbf{\phi}_{i,j}$ at given $\tau$ and
dividing by the number of non-recurrences (i.e., zeros in the RP):

$\displaystyle\Phi(\tau)\,=\,\frac{\sum_{i,j=1}^{N}\phi_{i,j}(\tau)}{\sum_{i,j=1}^{N}\bigl{(}1-R_{i,j}(\tau)\bigr{)}}.$
(5)

Figure 1: Recurrence flow for a noisy sinusoidal. (A) Time series of the noisy
sinusoidal and (B/C) RPs for two different embedding vectors $\vec{v}(t)$:
$\vec{v}(t)=[y(t),y(t+T)]$ and $\vec{v}(t)=[y(t),y(t+\frac{T}{4})]$ (left) and
flooded RPs/recurrence flow matrices (right). (D) Recurrence flow $\Phi(\tau)$
through the RP for varying delay $\tau$. The red dot marks the period $T$.
$\Phi(\tau)$ is maximized for delays that are odd multiples of $T/4$.

In a multivariate application, $\vec{v}$ can encompass time series from
different systems to study their cross-dependencies. In such a scenario, it is
more instructive to define the recurrence flow as a direct measure of
correlation/redundance. We, thus, define the recurrence flow measure of
redundance (RFMR) $\theta(\tau)$ as

$\displaystyle\theta(\tau)\,=\,1-\Phi(\tau).$ (6)

The significance of recurrence flow values can be tested against a random null
model based on uncorrelated white noise (App. B). Finally, it needs to be
noted that the idea of using RPs to identify optimal embedding parameters has
been considered before, but to our best knowledge has not been performed
systematically Zbilut and Webber Jr (1992); Atay and Altıntaş (1999).

##  3 Application to Model Examples

We now demonstrate the scope of the proposed method by highlighting two
different potential applications: the quantification of nonlinear correlations
(Sect. 3.1) and the identification of uniform embedding delays for TDE of
nonlinear signals (Sect. 3.2).

###  3.1 Nonlinear Dependence

We exemplify the efficacy of $\theta(\tau)$ as a nonlinear dependence measure
for deterministic systems with a simple bivariate system:

$\displaystyle\begin{split}x(t)\,&=\,\mathrm{sin}(2\pi
t/\omega)+\eta(t,\sigma_{1})\\\
y(t)\,&=\,ax(t-\tilde{\tau})^{2}+\eta(t,\sigma_{2})\end{split}$ (7)

with frequency $\omega=\left.2\pi\middle/T\right.$, time lag
$\tilde{\tau}=20$, and normal-distributed white noise processes
$\eta(t,\sigma)$ with standard deviations $\sigma_{1}$ and $\sigma_{2}$. This
system exhibits a sinusoidal cycle with frequency $\omega$ in its
$x$-component. The $y$-component is nonlinearly coupled to $x(t)$ and exhibits
a cycle with half of the period of $x(t)$. $y(t)$ follows $x(t)$ with a fixed
time lag $\tilde{\tau}$. Both components are superimposed by measurement noise
$\eta(t)$. We consider time series with $n=5,000$ samples (Fig. 2A). Due to
the specific coupling, the relationship between $x$ and $y$ is nonlinear (Fig.
2B).

We test whether we can detect the coupling and the corresponding time lag
$\tilde{\tau}$ by computing $\theta(\tau)$ for delays in the range
$\tau\in[-200,200]$ (Fig. 2C). In fact, we find that $\theta(\tau)$ reaches
local maxima at integer multiples of $\tilde{\tau}$, including
$\tau=\tilde{\tau}$ (red dashed line). Similar results are obtained if the
mutual information is used (Fig. 2C, black curve), confirming that nonlinear
relationships between deterministic time series can be captured by recurrence
flow in presence of measurement noise.

Figure 2: Detection of lagged, nonlinear dependence between nonlinearly
coupled sinusoidals. (A) The $y(t)$ time series (black) results from squaring
$x(t)$ (blue) and shifting it by a fixed time delay. (B) Nonlinear
relationship between $x(t)$ and $y(t)$. (C) Recurrence flow measure of
redundance $\theta(\tau)$ (blue) detecting the time delay between $x(t)$ and
$y(t)$ (red dashed line), confirmed by the cross-mutual information (black).

###  3.2 Uniform Time Delay Embedding

The ability of the recurrence flow to detect delayed dependencies between
multiple variables motivates its use in the delay selection problem faced in
TDE. Proficient delay selection must be based on a measure that captures the
redundancy in a (potentially large) set of correlated time series. A popular
solution is to use mutual information (MI). However, characterizing the
required joint probability density function $p(s_{1},s_{2},\dots,s_{m})$
becomes cumbersome for a large number of variables $m$ and is rendered
infeasible for many high-dimensional real-world systems. Compared to nonlinear
correlation measures like standard MI that are based on binning, recurrence
flow offers the advantage that computation times increase less rapidly: given
$k_{m}$ bins, an increase in dimensionality $m\rightarrow m+1$ results in
$(k_{m}-1)k_{m}$ additional bins while in the computation of $\Phi$, a
$k$-$d$-tree nearest neighbour search based RP computation increases only
linearly with the dimensionality of the system. It has yet to be noted that
more sophisticated nearest neighbour based approaches for MI computation do
not suffer from this drawback Kraskov et al. (2004).

Due to the popularity of this problem, other nonlinear correlation measures
that do not suffer from the curse of dimensionality have been conceived Ince
et al. (2017); Pecora et al. (2007). In order to validate the effectiveness of
the measure proposed here for selecting embedding delays, we do not only
compare it to the linear autocorrelation function (ACF) and the auto-mutual
information (MI), but also to the delay selection method proposed in Pecora et
al. (2007). The latter is based on a continuity statistic. In the following,
we estimate the optimal embedding dimension based on Cao’s method with a
threshold of $\Delta_{\mathrm{afn}}=0.2$ for the change of the number of
averaged false neighbors from $m\rightarrow m+1$ Cao (1997).

To illustrate the procedure, we consider a time series of the past 1 million
years of insolation on Earth at $52.39^{\circ}$ latitude (see App. A) Laskar
et al. (2004). The insolation depends on the Earth orbit and the Earth axis
tilt and precession, thus, varies in specific cycles (Milankovich cycles).
Cao’s method suggests an embedding dimension $m=4$. We present results for the
second and third component of the embedding vector. By using first-crossing of
ACF and first minimum of MI, both measures suggest a delay of
$\tau_{1}^{(\mathrm{ACF})}=\tau_{1}^{(\mathrm{MI})}=$6\text{\,}\mathrm{k}\mathrm{a}$$
(Fig. 3A). Since none of both measures is capable of selecting different
embedding delays for higher components of the embedding vector, this yields
the three-dimensional embedding vector
$\vec{v}(t)=\left[y(t),y(t-$6\text{\,}\mathrm{k}\mathrm{a}$),y(t-$12\text{\,}\mathrm{k}\mathrm{a}$)\right]$.
Next, we compute $\Phi(\tau)$ to check if the estimate of $\vec{v}(t)$ based
on the traditional TDE metrics is confirmed: we find
$\tau_{1}^{(\Phi)}=$5\text{\,}\mathrm{k}\mathrm{a}$$ in close agreement with
the ACF and MI-criteria (Fig. 3C). Moreover, the continuity statistic
$\langle\varepsilon^{*}\rangle(\tau)$ suggests the same embedding delay for
the first component as the auto-correlation and MI (first local maximum). It
does not indicate a global maximum at this delay, yielding a more ambiguous
choice of the optimal embedding delay than $\Phi(\tau)$. Finally, we examine
if both multi-dimensional measures suggest $\tau_{2}=2\tau_{1}$ for the second
embedding delay as expected for a traditional uniform time delay embedding
(UTDE). Interestingly, $\langle\varepsilon^{*}\rangle(\tau)$ offers only
limited information on an optimal embedding delay for the third component of
$\vec{v}(t)$ (Fig. 3D). Multiple local maxima offer a variety of choices with
no clear optimal value. We choose the global maximum (marked by star).
Conversely, $\Phi(\tau)$ once more provides a clear choice for the second
embedding delay with globally maximized flow for
$\tau_{2}^{(\Phi)}=2\tau_{1}^{(\Phi)}$.

We visually evaluate the quality of the resulting embeddings by comparing the
line structures in the corresponding RPs (Fig. 3E). The enlarged details of
the RPs illustrate how well deterministic intervals in the evolution of
insolation are resolved based on the phase space reconstructions yielded by
the ACF, continuity statistic, and recurrence flow. While the uniform
embedding vectors obtained from the ACF and $\Phi(\tau)$ result in well-
separated, undisturbed diagonal lines, multiple diagonal lines and the related
cycles are poorly expressed in the phase space suggested by
$\langle\varepsilon^{*}\rangle(\tau)$. The reconstructed phase space based on
the embedding vector obtained from the recurrence flow criterion reveals
several unstable periodic orbits (Fig. 3F), constituting a concentric spiral-
like phase space trajectory in three dimensions.

Many real-world systems allow taking only a relatively short series of
measurements for a single variable with high levels of superimposed
measurement noise. We study how well a known phase space of a paradigmatic
system can be reconstructed based on the four different measures considered
above with increasing noise strength. In particular, we generate $n=2,000$
samples of a Rössler system (see App. A) such that the resulting trajectory
only covers relatively few unstable periodic orbits. We reconstruct the known
three-dimensional phase space from the $y(t)$-component (Fig. 4A) with
superimposed uncorrelated white noise realizations (Fig. 4B).

Figure 3: Uniform delay selection for numerical insolation model. (A)
Insolation time series $y(t)$ in $\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$. (B)
Serial dependence of $y(t)$, measured in terms of ACF (red/blue) and auto-MI
(black) for the univariate time series, (C/D) as well as by the continuity
statistic $\langle\varepsilon^{*}\rangle(\tau)$ (gray) and the recurrence flow
$\Phi(\tau)$ (dark blue) for (C) two- and (D) three-dimensional embedding
vectors $\vec{v}(t)$. The optimal delay is marked by a
circle(/cross/star/vertical green line) for the ACF(/MI/continuity
statistic/recurrence flow), respectively. (E) Zoomed RPs for the ACF,
$\langle\varepsilon^{*}\rangle(\tau)$ and $\Phi(\tau)$ (from left to right).
(F) Three-dimensional phase space reconstruction based on recurrence flow
shows a spiral-like trajectory.

The noise strength (standard deviation of the noise) is varied in multiples of
the standard deviation $\sigma_{\mathrm{Roe}}$ of the undisturbed $y(t)$. Even
with only $10\%$ measurement noise, the original attractor is already
significantly less smooth (Fig. 4B). We compare the ACF, MI, continuity
statistic, and recurrence flow as delay selection measures while only uniform
time delays are considered, i.e., the optimal embedding delay is selected only
once for the step from a one- to a two-dimensional embedding. This ensures
that the two high-dimensional measures ($\Phi(\tau)$ and
$\langle\varepsilon^{*}\rangle(\tau)$) can be compared adequately to the
traditional measures.

To quantitatively evaluate the dependence of the reconstruction on noise
strength, we generate an RP for each reconstruction and for each of the four
delay selection methods. For each RP, we compute the joint recurrence rate
fraction (JRRF):

$\displaystyle\begin{split}JRRF\,&=\,\frac{\sum_{i,j}^{N}JR_{i,j}}{\sum_{i,j=1}^{N}R_{i,j}^{\mathrm{ref}}}\,,\qquad\mathrm{JRRF}\in[0,1]\\\
&\text{with}\\\
\mathbf{JR}\,&=\,\mathbf{R}^{\mathrm{ref}}\circ\mathbf{R}^{\mathrm{rec}}\end{split}$

from the RP of the (real) reference system $\mathbf{R}^{\mathrm{ref}}$ and the
RP of the respective reconstruction $\mathbf{R}^{\mathrm{rec}}$. We use it to
quantify the accordance of the real RP of the corresponding noisy Rössler
system to the reconstructions with respect to its recurrence structure (the
higher JRRF, the better the reconstruction).

We generate 50 individual noise realizations for each noise strength between
$0\%$ and $100\%$ of the original standard deviation $\sigma_{\mathrm{Roe}}$
and average the corresponding JRRF values (Fig. 4C). As expected, with
increasing noise level, the quality of the reconstruction decreases.

Figure 4: Sensitivity of recurrence flow against measurement noise, compared
to other dependence measures. (A) $y(t)$-component and reference phase space
of Rössler system without noise contamination; (B) the same, but with $10\%$
measurement noise. (C) Performance of ACF, auto-MI, continuity statistic
$\langle\varepsilon^{*}\rangle(\tau)$ and recurrence flow $\Phi(\tau)$ in
terms of JRRF for noise strengths from $0\%$ to $30\%$ (and $100\%$ inset).

For noise strengths between $10\%$ and $30\%$, the linear ACF performs worst
in terms of JRRF while the MI and continuity statistic perform equally well.
The almost perfect alignment of both might seem surprising but is due to the
discreteness immanent in the delay selection. While their agreement might be
interpreted in the way that this has to be the optimal delay (i.e., the
reconstructed system preserves most of the recurrence structure), the
recurrence flow $\Phi(\tau)$ shows superior performance with noise strengths
up to $30\%$, i.e., three times the noise level illustrated in Fig. 4B. Beyond
$30\%$, all four measures yield approximately the same performance (Fig. 4C
inset), as for JRRF$<0.5$, the alignment could be explained by random joint
recurrences.

###  3.3 Non-uniform Time Delay Embedding

For many real-world dynamical systems, it is not sufficient to consider only a
single characteristic time scale. Instead, multi-scale systems are governed by
a multitude of processes that imprint (quasi)-periodic cycles of various
lengths onto the measured time series. The selection of delays must account
for this complexity by considering non-uniform embedding delays (non-uniform
time delay embedding, NUTD). One of the most studied systems that exhibits
multi-scale dynamics is the El Niño-Southern Oscillation (ENSO). ENSO
represents a quasi-periodic climate pattern that is associated with spatio-
temporal variations of sea surface temperatures in the central and eastern
Pacific Ocean, oscillating between El Niño and La Niña events. We use the
delay differential ENSO model proposed in Ghil et al. (2008) to examine if the
recurrence flow can unveil distinct delays for a three-dimensional state space
reconstruction of model time series. The model is based on a nonlinear delay
differential equation and reproduces an abundance of key features of ENSO (see
App. A). We study two different solution types that are associated with
distinct dynamical regimes, i.e., a seasonal oscillation with superimposed
faster low-amplitude oscillations prior to a period-doubling and irregular
oscillations that are reminiscent of El Niño and La Niña events of random
magnitudes. A more detailed discussion of these solution types can be found in
Ghil et al. (2008).

The first solution type exhibits well-pronounced seasonal cycle and fast,
amplitude-modulated wiggles on top (Fig. 5A). Cao’s method suggests that this
solution type can be embedded in a $m=3$-dimensional embedding space. We find
that MI does not yield an unambiguous choice for an embedding delay while the
ACF suggests
$\tau_{1}^{(\mathrm{ACF})}=$0.25\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$$,
i.e., the expected value of a quarter of the seasonal cycle (Fig. 5B). The
same delay is identified with $\Phi(\tau)$ whereby
$\langle\varepsilon^{*}\rangle(\tau)$ yields a slightly higher optimal
embedding delay (Fig. 5C). While traditional UTDE now suggests
$\tau_{2}^{(\mathrm{ACF})}=2\tau_{1}^{(\mathrm{ACF})}$ for a three-dimensional
embedding, both $\Phi(\tau)$ and $\langle\varepsilon^{*}\rangle(\tau)$ instead
show that a different choice yields a superior phase space reconstruction in
terms of minimized redundancy (Fig. 5C/D). Both
$\langle\varepsilon^{*}\rangle(\tau)$ and $\Phi(\tau)$ effectively uncover the
faster cycle by means of local maxima. However, both also detect a delay that
is the sum of the seasonal and the fast cycle as a promising candidate. The
fact that the estimate of this conjoint cycle differs for both measures can be
explained by the different estimates on $\tau_{1}$. Since for $\Phi(\tau)$
both local maxima have the same height, we pick the first. The zoomed RPs
clearly express that both NUTD selection methods entail more coherent diagonal
lines with less perpendicular distortions (Fig. 5E). Despite the different
embedding vectors $\vec{v}^{(\Phi)}$ and $\vec{v}^{(\varepsilon^{*})}$, both
reconstructions give a convincing representation of the seasonal cycle in the
RP, respectively. The reconstructed phase space based on the delays selected
from the recurrence flow yields a clear visualization of the system’s periodic
oscillations (Fig. 5F).

Figure 5: Non-uniform delay selection for delay-differential ENSO model with
(A) periodic dynamics, analogously to Fig. 3. Serial dependence of $y(t)$,
measured in terms of (B) ACF (red/blue) and auto-MI (black) for the univariate
time series, as well as by (C/D) the continuity statistic
$\langle\varepsilon^{*}\rangle(\tau)$ (gray) and the recurrence flow
$\Phi(\tau)$ (dark blue) for (C) two- and (D) three-dimensional embedding
vectors $\vec{v}(t)$. The optimal delay is marked by a
circle(/cross/star/vertical green line) for the ACF(/MI/continuity
statistic/recurrence flow), respectively. (E) Zoomed RPs for the ACF,
$\langle\varepsilon^{*}\rangle(\tau)$ and $\Phi(\tau)$ (from left to right).
(F) Three-dimensional phase space reconstruction based on recurrence flow.

In the same manner, we study the second solution type of irregular ENSO-like
oscillations (Fig. 6). Cao’s method yields a 4-dimensional TDE. Both the ACF
and MI identify different but similar delays slightly larger than
$\tau_{1}=$0.35\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$$
(Fig. 6B). $\langle\varepsilon^{*}\rangle(\tau)$ and $\Phi(\tau)$ agree on an
embedding delay of
$\tau_{1}^{(\Phi)}=\tau_{1}^{(\varepsilon^{*})}=$0.33\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$$
which is slightly smaller than both $\tau_{1}^{(\mathrm{ACF})}$ and
$\tau_{1}^{(\mathrm{MI})}$ (Fig. 6C). Again, the estimates $\tau_{2}^{(\Phi)}$
and $\tau_{2}^{(\varepsilon^{*})}$ differ from the choice that would result
from UTDE and indicate that describing the dominant variability in the
observed irregular oscillations requires a multi-scale approach (Fig. 6D).
Slight deviations in both estimates give rise to minor discrepancies between
the emerging diagonal lines in the respective RPs (Fig. 6E). For none of the
three reconstructions, DLA are entirely removed which hints at an optimal
embedding dimension $m>3$ as identified by Cao’s method. This is supported by
the displayed attractor reconstruction that has a two-winged structure
reminiscent of the famous Lorenz attractor but potentially be unfolded further
(Fig. 6F). However, both NUTD selection measures once more provide a more
convincing result in terms of sparse, continuous diagonal lines than the ACF
and, thus, capture the system’s predictability more adequately.

Figure 6: Non-uniform delay selection for delay-differential ENSO model with
(A) irregular oscillations, analogously to Fig. 3. Serial dependence of
$y(t)$, measured in terms of (B) ACF (red/blue) and auto-MI (black) for the
univariate time series, as well as by (C/D) the continuity statistic
$\langle\varepsilon^{*}\rangle(\tau)$ (gray) and the recurrence flow
$\Phi(\tau)$ (dark blue) for (C) two- and (D) three-dimensional embedding
vectors $\vec{v}(t)$. The optimal delay is marked by a
circle(/cross/star/vertical green line) for the ACF(/MI/continuity
statistic/recurrence flow), respectively. (E) Zoomed RPs for the ACF,
$\langle\varepsilon^{*}\rangle(\tau)$ and $\Phi(\tau)$ (from left to right).
(F) Three-dimensional phase space reconstruction based on recurrence flow.

##  4 Conclusion

The nonlinearity and scale-dependence of relationships observed in high-
dimensional empirical data calls for appropriate and easily applicable
methods. Recurrence plots offer a mathematically simple yet effective
framework for the study of dependencies in high-dimensional dynamical systems
and are an established tool in applied nonlinear time series analysis. For
deterministic systems, diagonal lines in an RP yield valuable information on
the evolution of a system’s trajectory. We propose a novel recurrence based
dependence measure, called recurrence flow. It builds on the fact that an RP
can be computed from delayed copies of one (or multiple) time series that are
stacked on top of each other as it is common practise in time delay embedding.
The absence of spurious structures perpendicular to diagonal lines hints at a
correct choice of (embedding) delays between the embedding vector’s
coordinates. We exploit the emergence of these structures to characterize the
serial dependence in uni- and cross-dependence of multi-variate, high-
dimensional systems. We demonstrated that the recurrence flow $\Phi(\tau)$
captures nonlinear, lagged dependencies in presence of observational noise.
Due to its conceptual proximity to time delay embedding, we put a focus on the
delay selection problem that underlies attractor reconstruction. The
recurrence flow effectively pinpoints uniform and non-uniform embedding delays
in multi-dimensional nonlinear systems. It, thus, bares high potential to be
used in a recurrence based embedding algorithm which will be the focus of
future work. Compared to existing nonlinear dependence measures, it offers
several advantages: (i) it is not based on a binning procedure and does not
suffer from the curse of dimensionality (such as the mutual information), (ii)
it is able to detect non-uniform delays in high-dimensional multi-scale data,
(iii) it performs relatively well in the presence of observational noise, and,
(iv) for the considered examples, it yields well interpretable, unambiguous
maxima in the delay selection procedure (unlike, e.g., the continuity
statistic). The recurrence flow formalism appears promising for the
conceptualization of an automated recurrence based time delay embedding
algorithm.

## Code availability

Scripts and data used to perform this study are availabe via Zenodo.

## Acknowledgements

This research was supported by the Deutsche Forschungsgemeinschaft in the
context of the DFG project MA4759/11-1 “Nonlinear empirical mode analysis of
complex systems: Development of general approach and application in climate”.

## Conflict of interest

The authors declare that they have no conflict of interest.

##  Appendix A Appendix: Example Systems

To assess the performance of $\Phi(\tau)$ as a measure to select embedding
delays, the following systems are considered:

### Numerical Insolation Model

The complex superposition of gravitational forces between the earth and the
other planetary bodies in the solar system perturbs the earth’s rotation on an
elliptic orbit around the sun and its axial rotation. The cycles that manifest
due to these variations control the earth’s climate at time scales of millions
of years and are called Milankovich cycles. The nature of the underlying
perturbation renders the variations in insolation chaotic. The model proposed
in Laskar et al. (2004) considers all nine planets in the solar system and
describes the orbit of the moon separately. A Hamiltonian, consisting of an
integrable and perturbation component, is numerically integrated with a
symplectic integrator scheme (SABAC4). Several dissipative effects (tides,
core-mantle friction, climate friction) are included. The model returns time
series for the earth’s orbit’s eccentricity, climatic precession, obliquity,
and insolation, from which we study only the latter with $n=1,000$ samples.
For the computation of RPs, we fix the recurrence rate to $5\%$.

### Rössler System

The Rössler system is a three-dimensional, continuous dynamical system that
generates a strange attractor:

$\displaystyle\begin{split}\dot{x}\,&=\,-y-z,\\\ \dot{y}\,&=\,x+ay,\\\
\dot{z}\,&=\,b+z(x-c).\end{split}$

The time series only covers a few unstable periodic orbits with $3,000$
samples from which we discard $1,000$ as transients, resulting in $N=2,000$.
Uniform sampling intervals are fixed as $\Delta t=0.10$. We set
$a=0.02925,\,b=0.1$ and $c=8.5$, ensuring chaotic dynamics. Independent
realizations of uncorrelated white noise with different noise strengths
between $0\%$ and $100\%$ are superimposed on the $y(t)$-component to mimic
measurement noise. For the computation of RPs, we fix the recurrence rate to
$8\%$ regardless of the noise strength.

### Delay Differential ENSO Model

As an example for a system with multiple characteristic time scales, we
examine a delay differential model of ENSO. On top of the seasonal mode, ENSO
represents the predominant mode of sea surface temperature (SST) variability
in the tropical Pacific. The SSTs exhibit recurring variations with time
scales between 2 and 7 years between two regimes of well-distinguishable SST
anomalies: El Niño (warming phase) and La Niña (cooling phase). These
variations disturb large-scale air transport in the tropics and induce a
multitude of global climatic impacts, e.g., droughts and floods in Australia
or South America. Several conceptual models have reproduced key features of
this oscillation by including hypotheses on the mechanistic origins of ENSO,
including negative and positive feedbacks of temperature anomalies and
atmospheric circulation and potential resonance phenomena with the seasonal
forcing. The model studied here is taken from Ghil et al. (2008) and mimics
ENSO dynamics based on two key mechanisms, i.e., delayed negative feedback and
seasonal forcing

$\displaystyle\dfrac{\mathrm{d}y(t)}{\mathrm{d}t}\,=\,-\mathrm{tanh}\bigl{(}\kappa
y(t-\zeta)\bigr{)}+b\,\mathrm{cos}\left(2\pi\omega t\right).$

We set the frequency of the periodic forcing to seasonal forcing ($\omega=1$)
and fix $b=1$. Variations in the delay $\zeta$ and the parameter $\kappa$ give
rise to dynamically distinct time series. We study two solution types: a
regular solution with a seasonal cycle and fast, amplitude-modulated wiggles
$(\kappa=100,\,\zeta=0.025)$ and a solution of irregular ENSO-like
oscillations with stochastic amplitude variations $(\kappa=50,\,\zeta=0.42)$.
For both solution types, we generate $n=10,000$ values. For the computation of
RPs, we fix the recurrence rate to $5\%$.

##  Appendix B Appendix: Statistical Significance

For an uncorrelated white noise time series of infinite length, the presence
of recurrences along a diagonal can be described by a binomial distribution.
The flow along each diagonal can consequently be regarded as an idealized
sequence of Bernoulli trials for which a success is equivalent to a
recurrence, i.e., a black pixel. The probability of having $X$ unsuccessful
Bernoulli trials (no recurrences) until a trial succeeds is given by the
geometric distribution

$\displaystyle P(X=n)\,=\,(1-p)^{k-1}p$ (8)

with expectation value $\left.1\middle/p\right.$. The probability of success
$p$ for each Bernoulli trial is given by the recurrence rate and depends on
$\varepsilon$. It follows that we can derive the recurrence flow for an
idealized uncorrelated white noise time series of length $n$ as

$\displaystyle\theta(\varepsilon)\,=\,1-\frac{n}{p(\varepsilon)}.$ (9)

Given an observational time series of length $n$ and a suitable choice for the
vicinity threshold $\varepsilon$, $p(\varepsilon)$ can be identified with the
recurrence rate and eq. (7) can be used to test whether the recurrence flow of
the real signal can be distinguished from an uncorrelated random process.

## References

  * Kwapień et al. [2015] Jarosław Kwapień, Paweł Oświęcimka, and Stanisław Drożdż.  Detrended fluctuation analysis made flexible to detect range of cross-correlated fluctuations.  _Physical Review E_ , 92(5):052815, 2015.
  * Newman [2018] Mark Newman.  _Networks_.  Oxford university press, 2018.
  * Feldhoff et al. [2012] J. H. Feldhoff, R. V. Donner, J. F. Donges, N. Marwan, and J. Kurths.  Geometric detection of coupling directions by means of inter-system recurrence networks.  _Physics Letters A_ , 376(46):3504–3513, 2012\.  doi: 10.1016/j.physleta.2012.10.008.
  * Runge et al. [2015] Jakob Runge, Vladimir Petoukhov, Jonathan F Donges, Jaroslav Hlinka, Nikola Jajcay, Martin Vejmelka, David Hartman, Norbert Marwan, Milan Paluš, and Jürgen Kurths.  Identifying causal gateways and mediators in complex spatio-temporal systems.  _Nature communications_ , 6(1):1–10, 2015.
  * Runge et al. [2019] Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovic.  Detecting and quantifying causal associations in large nonlinear time series datasets.  _Science advances_ , 5(11):eaau4996, 2019.
  * Ramos et al. [2017] A. M. T. Ramos, A. Builes-Jaramillo, G. Poveda, B. Goswami, E. E. N. Macau, J. Kurths, and N. Marwan.  Recurrence measure of conditional dependence and applications.  _Physical Review E_ , 95:052206, 2017.  doi: 10.1103/PhysRevE.95.052206.
  * Dakos et al. [2008] Vasilis Dakos, Marten Scheffer, Egbert H van Nes, Victor Brovkin, Vladimir Petoukhov, and Hermann Held.  Slowing down as an early warning signal for abrupt climate change.  _Proceedings of the National Academy of Sciences_ , 105(38):14308–14312, 2008.
  * Boettner et al. [2021] Christopher Boettner, Georg Klinghammer, Niklas Boers, Thomas Westerhold, and Norbert Marwan.  Early-warning signals for cenozoic climate transitions.  _Quaternary Science Reviews_ , 270:107177, 2021.
  * Boers [2018] Niklas Boers.  Early-warning signals for dansgaard-oeschger events in a high-resolution ice core record.  _Nature communications_ , 9(1):1–8, 2018.
  * Balasis et al. [2013] Georgios Balasis, Reik V Donner, Stelios M Potirakis, Jakob Runge, Constantinos Papadimitriou, Ioannis A Daglis, Konstantinos Eftaxias, and Jürgen Kurths.  Statistical mechanics and information-theoretic perspectives on complexity in the earth system.  _Entropy_ , 15(11):4844–4888, 2013.
  * Sun et al. [2014] Jie Sun, Carlo Cafaro, and Erik M Bollt.  Identifying the coupling structure in complex systems through the optimal causation entropy principle.  _Entropy_ , 16(6):3416–3433, 2014.
  * Porfiri et al. [2019] Maurizio Porfiri, Raghu Ram Sattanapalle, Shinnosuke Nakayama, James Macinko, and Rifat Sipahi.  Media coverage and firearm acquisition in the aftermath of a mass shooting.  _Nature human behaviour_ , 3(9):913–921, 2019\.
  * Pompe and Runge [2011] Bernd Pompe and Jakob Runge.  Momentary information transfer as a coupling measure of time series.  _Physical Review E_ , 83(5):051122, 2011.
  * Kraskov et al. [2004] Alexander Kraskov, Harald Stögbauer, and Peter Grassberger.  Estimating mutual information.  _Physical review E_ , 69(6):066138, 2004.
  * Nawroth and Peinke [2006] Andreas P Nawroth and Joachim Peinke.  Multiscale reconstruction of time series.  _Physics Letters A_ , 360(2):234–237, 2006.
  * Maraun and Kurths [2004] Douglas Maraun and Jürgen Kurths.  Cross wavelet analysis: significance testing and pitfalls.  _Nonlinear Processes in Geophysics_ , 11(4):505–514, 2004.
  * Agarwal et al. [2018] Ankit Agarwal, Rathinasamy Maheswaran, Norbert Marwan, Levke Caesar, and Jürgen Kurths.  Wavelet-based multiscale similarity measure for complex networks.  _The European Physical Journal B_ , 91(11):1–12, 2018.
  * Braun et al. [2021] Tobias Braun, Vishnu R Unni, RI Sujith, Juergen Kurths, and Norbert Marwan.  Detection of dynamical regime transitions with lacunarity as a multiscale recurrence quantification measure.  _Nonlinear Dynamics_ , 104(4):3955–3973, 2021\.
  * Vlachos and Kugiumtzis [2010] Ioannis Vlachos and Dimitris Kugiumtzis.  Nonuniform state-space reconstruction and coupling detection.  _Physical Review E_ , 82(1):016207, 2010.
  * Pecora et al. [2007] Louis M Pecora, Linda Moniz, Jonathan Nichols, and Thomas L Carroll.  A unified approach to attractor reconstruction.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 17(1):013110, 2007.
  * Eckmann et al. [1987] J.-P. Eckmann, S. Oliffson Kamphorst, and D. Ruelle.  Recurrence Plots of Dynamical Systems.  _Europhysics Letters_ , 4(9):973–977, 1987.  doi: 10.1209/0295-5075/4/9/004.
  * Marwan et al. [2007] Norbert Marwan, M Carmen Romano, Marco Thiel, and Jürgen Kurths.  Recurrence plots for the analysis of complex systems.  _Physics reports_ , 438(5-6):237–329, 2007.
  * Thiel et al. [2004] M. Thiel, M. C. Romano, P. L. Read, and J. Kurths.  Estimation of dynamical invariants without embedding by recurrence plots.  _Chaos_ , 14(2):234–243, 2004.  doi: 10.1063/1.1667633.
  * Schinkel et al. [2008] S. Schinkel, O. Dimigen, and N. Marwan.  Selection of recurrence threshold for signal detection.  _European Physical Journal – Special Topics_ , 164(1):45–53, 2008.  doi: 10.1140/epjst/e2008-00833-5.
  * Klimaszewska and Żebrowski [2009] K. Klimaszewska and J. J. Żebrowski.  Detection of the type of intermittency using characteristic patterns in recurrence plots.  _Physical Review E_ , 80:026214, 2009.  doi: 10.1103/PhysRevE.80.026214.
  * Hirata and Aihara [2010] Y. Hirata and K. Aihara.  Devaney’s chaos on recurrence plots.  _Physical Review E_ , 82:036209, 2010.  doi: 10.1103/PhysRevE.82.036209.
  * Marwan et al. [2002] N. Marwan, N. Wessel, U. Meyerfeldt, A. Schirdewan, and J. Kurths.  Recurrence Plot Based Measures of Complexity and its Application to Heart Rate Variability Data.  _Physical Review E_ , 66(2):026702, 2002.  doi: 10.1103/PhysRevE.66.026702.
  * Marwan et al. [2021] Norbert Marwan, Jonathan F Donges, Reik V Donner, and Deniz Eroglu.  Nonlinear time series analysis of palaeoclimate proxy records.  _Quaternary Science Reviews_ , 274:107245, 2021.
  * Romano et al. [2007] M Carmen Romano, Marco Thiel, Jürgen Kurths, and Celso Grebogi.  Estimation of the direction of the coupling by conditional probabilities of recurrence.  _Physical Review E_ , 76(3):036211, 2007.
  * Nkomidio et al. [2022] A. M. Nkomidio, E. K. Ngamga, B. R. N. Nbendjo, J. Kurths, and N Marwan.  Recurrence-Based Synchronization Analysis of Weakly Coupled Bursting Neurons Under External ELF Fields.  _Entropy_ , 24(2):235, 2022.  doi: 10.3390/e24020235.
  * Zbilut and Marwan [2008] Joseph P Zbilut and Norbert Marwan.  The wiener–khinchin theorem and recurrence quantification.  _Physics Letters A_ , 372(44):6622–6626, 2008\.
  * Bradley and Mantilla [2002] Elizabeth Bradley and Ricardo Mantilla.  Recurrence plots and unstable periodic orbits.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 12(3):596–600, 2002.
  * Banerjee et al. [2021] Abhirup Banerjee, Bedartha Goswami, Yoshito Hirata, Deniz Eroglu, Bruno Merz, Jürgen Kurths, and Norbert Marwan.  Recurrence analysis of extreme event-like data.  _Nonlinear Processes in Geophysics_ , 28(2):213–229, 2021.
  * Ozdes and Eroglu [2022] Celik Ozdes and Deniz Eroglu.  Transformation cost spectrum for irregularly sampled time series.  _The European Physical Journal Special Topics_ , pages 1–12, 2022\.
  * Romano et al. [2004] M. C. Romano, M. Thiel, and J. Kurths.  Generalized Synchronization Indices based on Recurrence in Phase Space.  volume 742, pages 330–336, 2004.  doi: 10.1063/1.1846494.
  * Romano et al. [2005] M. C. Romano, M. Thiel, J. Kurths, I. Z. Kiss, and J. L. Hudson.  Detection of synchronization for non-phase-coherent and non-stationary data.  _Europhysics Letters_ , 71(3):466–472, 2005.  doi: 10.1209/epl/i2005-10095-1.
  * Senthilkumar et al. [2008] D. V. Senthilkumar, M. Lakshmanan, and J. Kurths.  Transition from phase to generalized synchronization in time-delay systems.  _Chaos_ , 18(2):023118, 2008.  doi: 10.1063/1.2911541.
  * Zou et al. [2011] Yong Zou, M Carmen Romano, Marco Thiel, Norbert Marwan, and Jürgen Kurths.  Inferring indirect coupling by means of recurrences.  _International Journal of Bifurcation and Chaos_ , 21(04):1099–1111, 2011.
  * Goswami et al. [2013] Bedartha Goswami, Norbert Marwan, Georg Feulner, and Jürgen Kurths.  How do global temperature drivers influence each other?  _The European Physical Journal Special Topics_ , 222(3):861–873, 2013.
  * Porfiri and Ruiz Marín [2019] Maurizio Porfiri and Manuel Ruiz Marín.  Transfer entropy on symbolic recurrences.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 29(6):063123, 2019.
  * Westerhold et al. [2020] Thomas Westerhold, Norbert Marwan, Anna Joy Drury, Diederik Liebrand, Claudia Agnini, Eleni Anagnostou, James SK Barnet, Steven M Bohaty, David De Vleeschouwer, Fabio Florindo, et al.  An astronomically dated record of earth’s climate and its predictability over the last 66 million years.  _Science_ , 369(6509):1383–1387, 2020.
  * Thiel et al. [2003] Marco Thiel, M Carmen Romano, and Jurgen Kurths.  Analytical description of recurrence plots of white noise and chaotic processes.  _arXiv preprint nlin/0301027_ , 2003.
  * Marwan [2011] Norbert Marwan.  How to avoid potential pitfalls in recurrence plot based data analysis.  _International Journal of Bifurcation and Chaos_ , 21(04):1003–1017, 2011.
  * Kraemer and Marwan [2019] K Hauke Kraemer and Norbert Marwan.  Border effect corrections for diagonal line based recurrence quantification analysis measures.  _Physics Letters A_ , 383(34):125977, 2019.
  * Thiel et al. [2006] Marco Thiel, M Carmen Romano, and Jürgen Kurths.  Spurious structures in recurrence plots induced by embedding.  _Nonlinear Dynamics_ , 44(1):299–305, 2006.
  * Marwan and Kurths [2005] N. Marwan and J. Kurths.  Line structures in recurrence plots.  _Physics Letters A_ , 336(4–5):349–357, 2005\.  doi: 10.1016/j.physleta.2004.12.056.
  * Kraemer et al. [2022] K Hauke Kraemer, Maximilian Gelbrecht, Induja Pavithran, RI Sujith, and Norbert Marwan.  Optimal state space reconstruction via monte carlo decision tree search.  _Nonlinear Dynamics_ , pages 1–21, 2022.
  * Zbilut and Webber Jr [1992] Joseph P Zbilut and Charles L Webber Jr.  Embeddings and delays as derived from quantification of recurrence plots.  _Physics letters A_ , 171(3-4):199–203, 1992\.
  * Atay and Altıntaş [1999] Fatihcan M Atay and Yiğit Altıntaş.  Recovering smooth dynamics from time series with the aid of recurrence plots.  _Physical Review E_ , 59(6):6593, 1999.
  * Ince et al. [2017] Robin AA Ince, Bruno L Giordano, Christoph Kayser, Guillaume A Rousselet, Joachim Gross, and Philippe G Schyns.  A statistical framework for neuroimaging data analysis based on mutual information estimated via a gaussian copula.  _Human brain mapping_ , 38(3):1541–1573, 2017\.
  * Cao [1997] Liangyue Cao.  Practical method for determining the minimum embedding dimension of a scalar time series.  _Physica D: Nonlinear Phenomena_ , 110(1-2):43–50, 1997.
  * Laskar et al. [2004] Jacques Laskar, Philippe Robutel, Frédéric Joutel, Mickael Gastineau, ACM Correia, and Benjamin Levrard.  A long-term numerical solution for the insolation quantities of the earth.  _Astronomy & Astrophysics_, 428(1):261–285, 2004.
  * Ghil et al. [2008] Michael Ghil, Ilya Zaliapin, and Sylvester Thompson.  A delay differential model of enso variability: parametric instability and the distribution of extremes.  _Nonlinear Processes in Geophysics_ , 15(3):417–433, 2008.


$\displaystyle\hat{T}_{00}^{\rm gf}={}$
$\displaystyle\Bigl{[}\hat{T}_{00}^{\rm gf}\Bigr{]}_{\scriptscriptstyle\rm
W}+a^{2-D}\,{\rm
Re}\Bigl{(}\bigl{[}\hat{K}_{2},\hat{A}_{0}\bigr{]}+\bigl{[}\partial_{k}\hat{K}_{1},\hat{A}_{k}\bigr{]}\Bigr{)}\,,$
(226) $\displaystyle\hat{T}_{0i}^{\rm gf}={}$
$\displaystyle\Bigl{[}\hat{T}_{0i}^{\rm gf}\Bigr{]}_{\scriptscriptstyle\rm
W}+a^{2-D}\,{\rm
Re}\Bigl{(}\bigl{[}\partial_{i}\hat{K}_{1},\hat{A}_{0}\bigr{]}+\bigl{[}\hat{K}_{2},\hat{A}_{i}\bigr{]}\Bigr{)}\,,$
(227) $\displaystyle\hat{T}_{ij}^{\rm gf}={}$
$\displaystyle\Bigl{[}\hat{T}_{ij}^{\rm gf}\Bigr{]}_{\scriptscriptstyle\rm
W}+a^{2-D}\,{\rm
Re}\Bigl{(}2\bigl{[}\partial_{(i}\hat{K}_{1},\hat{A}_{j)}\bigr{]}+\delta_{ij}\bigl{[}\hat{K}_{2},\hat{A}_{0}\bigr{]}-\delta_{ij}\bigl{[}\partial_{k}\hat{K}_{1},\hat{A}_{k}\bigr{]}\Bigr{)}\,.$
(228)

The Weyl-ordered parts above can easily be defined from the covariant
expression (221) directly, by simply symmetrizing the products,

$\displaystyle\Bigl{[}T_{\mu\nu}^{\rm gf}\Bigr{]}_{\scriptscriptstyle\rm
W}={}$
$\displaystyle-\frac{1}{\xi}\Bigl{\\{}\hat{A}_{(\mu},\Bigl{(}\nabla_{\nu)}+2\zeta
n_{\nu)}\Bigr{)}\Bigl{(}\nabla^{\rho}\hat{A}_{\rho}-2\zeta
n^{\rho}\hat{A}_{\rho}\Bigr{)}\Bigr{\\}}$ (229) $\displaystyle\hskip
5.69046pt+\frac{g_{\mu\nu}}{2\xi}\biggl{[}\Bigl{\\{}\hat{A}_{\rho},\Bigl{(}\nabla^{\rho}+2\zeta
n^{\rho}\Bigr{)}\Bigl{(}\nabla^{\sigma}\hat{A}_{\sigma}-2\zeta
n^{\sigma}\hat{A}_{\sigma}\Bigr{)}\Bigr{\\}}+\Bigl{(}\nabla^{\rho}\hat{A}_{\rho}-2\zeta
n^{\rho}\hat{A}_{\rho}\Bigr{)}^{\\!2}\,\biggr{]}\,.$

The non-vanishing commutators, that are just $c$-numbers, can be evaluated as
coincidence limits of position space commutators from (190)–(198). This leads
to the following form for the gauge-fixing part of the energy-momentum tensor
operator,

$\displaystyle\hat{T}_{\mu\nu}^{\rm gf}={}$
$\displaystyle\Bigl{[}\hat{T}_{\mu\nu}^{\rm gf}\Bigr{]}_{\scriptscriptstyle\rm
W}-{\rm
Re}\Biggl{\\{}\biggl{[}2\Bigl{(}\nabla_{(\mu}\nabla^{\prime}_{\nu)}+2\zeta
n_{(\mu}\nabla^{\prime}_{\nu)}\Bigr{)}$ $\displaystyle\hskip 71.13188pt-
g_{\mu\nu}\Bigl{(}\nabla^{\rho}\nabla^{\prime}_{\rho}+2\zeta
n^{\rho}\nabla^{\prime}_{\rho}\Bigr{)}\biggr{]}\biggl{[}\Bigl{(}\frac{a^{\prime}}{a}\Bigr{)}^{\\!\zeta}i\bigl{[}\tensor*[^{\scriptscriptstyle\\!-\\!}]{\Delta}{{}^{\scriptscriptstyle\\!+\\!}}\bigr{]}(x;x^{\prime})\biggr{]}_{x^{\prime}\to
x}\Biggr{\\}}\,.$ (230)

Here we can recognize the FP ghost two-point function (202) in the $c$-number
contribution. This contribution matches exactly the contribution that would
descend from the FP ghost action (205). Thus we see that even in the linear
Abelian gauge theory, and even without being explicitly introduced, the FP
ghosts naturally arise as commutators accounting for the difference between
proper operator ordering of operators and Weyl ordering of operators. Either
form of the definition we adopt, (224) or (230), they both lead to the
vanishing expectation value (225).

Had we ignored the question of operator ordering and defined (229) as the
observable, its expectation value would no longer vanish identically. This is
what is behind the conclusion in BeltranJimenez:2008enx ;
BeltranJimenez:2009oao ; BeltranJimenez:2011nm  that the gauge-fixing
contribution to the energy-momentum tensor of the photon engenders a
nonvanishing cosmological constant contribution. This conclusion, however,
does not hold up. Ref. Zhang:2022csl  attempted to address the problem in de
Sitter space by considering Weyl ordering of operators, but employing
adiabatic subtraction to obtain (225). That approach is not appropriate
either, as it leaves the option for the gauge-fixing part of the energy-
momentum tensor to produce a physical contribution in some spacetime other
than de Sitter. Moreover, their conclusion that FP ghosts are not necessary
when operators are Weyl ordered contradicts the results of this section, as
well as contradicting consistent analogous for Stueckelberg vector fields
Belokogne:2015etf ; Belokogne:2016dvd .

##  9 Discussion

In this work we considered the canonical quantization of the photon (massless
vector field) in spatially flat FLRW spacetimes in the two-parameter family of
linear gauges (1). We used this framework to demonstrate that observables with
appropriate operator ordering respect expectations set by the correspondence
principle. In particular, we had demonstrated how the gauge-fixing term (1)
does not contribute to the energy-momentum tensor expectation value for any
physical state.

Along the way we had elucidated how the canonical formulation is the
appropriate framework for the Gupta-Bleuler quantization, and had derived the
subsidiary condition on the physical state from the first-class constraint
structure of the classical theory. This subsidiary condition translates into
subsidiary conditions on two-point functions derived in Sec. 7.1. The form of
these subsidiary conditions depends on the two gauge-fixing parameters $\xi$
and $\zeta$ from (1), as in fact does the two-point function itself. Two-point
functions with free gauge-fixing parameters are particularly useful for
computations, as they allow for explicit checks of gauge independence of
observables, that should not depend on the gauge-fixing parameters. The
construction of inflationary gauge-independent quantum observables in still in
its early stages Miao:2017feh .

Different two-point functions considered in Sec. 7 are basic building blocks
of non-equilibrium loop computations in the Schwinger-Keldysh formalism
appropriate for early universe cosmology. We believe that the framework set up
in this paper will facilitate the construction of photon propagators in
realistic inflationary spacetimes, that will in turn allow the investigation
of slow-roll corrections to large infrared effects found when photons interact
with spectator scalars and gravitons in de Sitter Prokopec:2002jn ;
Prokopec:2002uw ; Prokopec:2003bx ; Prokopec:2003iu ; Kahya:2005kj ;
Kahya:2006ui ; Prokopec:2006ue ; Prokopec:2007ak ; Prokopec:2008gw ;
Leonard:2012si ; Leonard:2012ex ; Chen:2016nrs ; Chen:2016uwp ; Chen:2016hrz ;
Kaya:2018qbj ; Popov:2017xut ; Glavan:2019uni ; Leonard:2012fs ;
Leonard:2013xsa ; Glavan:2013jca ; Wang:2014tza ; Glavan:2015ura ;
Wang:2015eaa ; Glavan:2016bvp ; Miao:2018bol . We had expressed the photon
two-point function in Sec. 7.2 in terms of several scalar mode functions. The
equations of motion that these scalar mode functions satisfy are collected in
Sec. 5.3. In that form it should be considerably easier to identify convenient
choices for gauge-fixing parameters that lead to simpler propagators for given
inflationary backgrounds. Constructing such two-point functions will be the
subject of future work Domazet:2023 . Furthermore, collecting subsidiary
conditions for two-point functions in Sec. 7.1, both in the double derivative
and the single derivative form, will facilitate consistency checks of photon
two-point functions constructed in future studies.

Before embarking to compute the photon propagators in power-law and slow-roll
inflation it is instructive to check the existing two-point functions in de
Sitter space ($\epsilon\\!=\\!0$) from the literature versus subsidiary
conditions for the Wightman function (186) and for the Feynman propagator
(187). The checks are summarized in Table 1.

de Sitter propagator | $\ \ \bigl{\langle}\hat{\Pi}_{0}\hat{\Pi}_{0}\bigr{\rangle}\ \ $ | $\ \bigl{\langle}\hat{\Pi}_{0}\partial_{i}\hat{\Pi}_{i}\bigr{\rangle}\ $ | $\bigl{\langle}\partial_{i}\hat{\Pi}_{i}\partial_{j}\hat{\Pi}_{j}\bigr{\rangle}$
---|---|---|---
Allen-Jacobson Allen:1985wd  ($\xi\\!=\\!1,\zeta\\!=\\!0,D$) | $\boldsymbol{\times}$ | $\checkmark$ | $\checkmark$
Tsamis-Woodard Tsamis:2006gj  ($\xi\\!=\\!0,\zeta\\!=\\!0,D$) | $\checkmark$ | $\checkmark$ | $\checkmark$
Youssef Youssef:2010dw  ($\xi,\zeta\\!=\\!0,D\\!=\\!4$) | $\boldsymbol{\times}$ | $\checkmark$ | $\checkmark$
Fröb-Higuchi Frob:2013qsa  ($\xi,\zeta\\!=\\!0,D$) | $\boldsymbol{\times}$ | $\checkmark$ | $\checkmark$
Woodard Woodard:2004ut  ($\xi\\!=\\!1,\zeta\\!=\\!1,D$) | $\checkmark$ | $\checkmark$ | $\checkmark$
Table 1:  Checks of subsidiary conditions for massless photon two-point
functions in de Sitter reported in the literature. The checks refer to both
the Wightman functions and Feynman propagators. The first four entries with
$\zeta\\!=\\!0$ correspond to covariant gauges, while the fifth entry with
$\zeta\\!=\\!1$ corresponds to a conformal gauge in $D\\!=\\!4$.

Unexpectedly, they reveal that the general covariant gauge two-point functions
satisfy all the required subsidiary conditions only in the exact transverse
gauge limit $\xi\\!\to\\!0$. For non-vanishing $\xi$ the double divergence of
the covariant gauge two-point functions fails to vanish off-coincidence,

$\nabla^{\mu}\nabla^{\prime\nu}i\bigl{[}\tensor*[_{\mu}^{\scriptscriptstyle-\\!}]{\Delta}{{}_{\nu}^{\scriptscriptstyle\\!+}}\bigr{]}(x;x^{\prime})=-\frac{\xi
H^{D}\,\Gamma(D)}{(4\pi)^{\frac{D}{2}}\Gamma\bigl{(}\frac{D}{2}\bigr{)}}=\nabla^{\mu}\nabla^{\prime\nu}i\bigl{[}\tensor*[_{\mu}^{\scriptscriptstyle+\\!}]{\Delta}{{}_{\nu}^{\scriptscriptstyle\\!+}}\bigr{]}(x;x^{\prime})\,,\qquad
x\neq x^{\prime}\,,$ (231)

violating subsidiary conditions (186a) and (187a). This points to
inconsistencies of known results even for the relatively simple case of the
maximally symmetric de Sitter space. These inconsistencies will be examined in
more detail and addressed elsewhere Glavan:2022dwb ; Glavan:2022nrd .

## Acknowledgements

I am grateful to Chun-Yen Lin, Danijel Jurman, and Jose Beltrán Jiménez for
discussions on the topic, and especially to Tomislav Prokopec for the critical
reading of the manuscript. This work was supported by the Czech Science
Foundation (GAČR) grant 20-28525S.

## References

  * (1) T. Prokopec, O. Tornkvist and R. P. Woodard, “Photon mass from inflation,” Phys. Rev. Lett. 89 (2002), 101301 [arXiv:astro-ph/0205331 [astro-ph]].
  * (2) T. Prokopec, O. Tornkvist and R. P. Woodard, “One loop vacuum polarization in a locally de Sitter background,” Annals Phys. 303 (2003), 251-274 [arXiv:gr-qc/0205130 [gr-qc]].
  * (3) T. Prokopec and R. P. Woodard, “Vacuum polarization and photon mass in inflation,” Am. J. Phys. 72 (2004), 60-72 [arXiv:astro-ph/0303358 [astro-ph]].
  * (4) T. Prokopec and R. P. Woodard, “Dynamics of superhorizon photons during inflation with vacuum polarization,” Annals Phys. 312 (2004), 1-16 [arXiv:gr-qc/0310056 [gr-qc]].
  * (5) E. O. Kahya and R. P. Woodard, “Charged scalar self-mass during inflation,” Phys. Rev. D 72 (2005), 104001 [arXiv:gr-qc/0508015 [gr-qc]].
  * (6) E. O. Kahya and R. P. Woodard, “One Loop Corrected Mode Functions for SQED during Inflation,” Phys. Rev. D 74 (2006), 084012 [arXiv:gr-qc/0608049 [gr-qc]].
  * (7) T. Prokopec, N. C. Tsamis and R. P. Woodard, “Two Loop Scalar Bilinears for Inflationary SQED,” Class. Quant. Grav. 24 (2007), 201-230 [arXiv:gr-qc/0607094 [gr-qc]].
  * (8) T. Prokopec, N. C. Tsamis and R. P. Woodard, “Stochastic Inflationary Scalar Electrodynamics,” Annals Phys. 323 (2008), 1324-1360 [arXiv:0707.0847 [gr-qc]].
  * (9) T. Prokopec, N. C. Tsamis and R. P. Woodard, “Two loop stress-energy tensor for inflationary scalar electrodynamics,” Phys. Rev. D 78 (2008), 043523 [arXiv:0802.3673 [gr-qc]].
  * (10) K. E. Leonard, T. Prokopec and R. P. Woodard, “Covariant Vacuum Polarizations on de Sitter Background,” Phys. Rev. D 87 (2013) no.4, 044030 [arXiv:1210.6968 [gr-qc]].
  * (11) K. E. Leonard, T. Prokopec and R. P. Woodard, “Representing the Vacuum Polarization on de Sitter,” J. Math. Phys. 54 (2013), 032301 [arXiv:1211.1342 [gr-qc]].
  * (12) X. Chen, Y. Wang and Z. Z. Xianyu, “Loop Corrections to Standard Model Fields in Inflation,” JHEP 08 (2016), 051 [arXiv:1604.07841 [hep-th]].
  * (13) X. Chen, Y. Wang and Z. Z. Xianyu, “Standard Model Background of the Cosmological Collider,” Phys. Rev. Lett. 118 (2017) no.26, 261302 [arXiv:1610.06597 [hep-th]].
  * (14) X. Chen, Y. Wang and Z. Z. Xianyu, “Standard Model Mass Spectrum in Inflationary Universe,” JHEP 04 (2017), 058 [arXiv:1612.08122 [hep-th]].
  * (15) A. Kaya, ‘Superhorizon Electromagnetic Field Background from Higgs Loops in Inflation,” JCAP 03 (2018), 046 [arXiv:1801.02032 [gr-qc]].
  * (16) F. K. Popov, “Debye mass in de Sitter space,” JHEP 06 (2018), 033 [arXiv:1711.11010 [hep-th]].
  * (17) D. Glavan and G. Rigopoulos, “One-loop electromagnetic correlators of SQED in power-law inflation,” JCAP 02 (2021), 021 [arXiv:1909.11741 [gr-qc]].
  * (18) K. E. Leonard and R. P. Woodard, “Graviton Corrections to Maxwell’s Equations,” Phys. Rev. D 85 (2012), 104048 [arXiv:1202.5800 [gr-qc]].
  * (19) K. E. Leonard and R. P. Woodard, “Graviton Corrections to Vacuum Polarization during Inflation,” Class. Quant. Grav. 31 (2014), 015010 [arXiv:1304.7265 [gr-qc]].
  * (20) D. Glavan, S. P. Miao, T. Prokopec and R. P. Woodard, “Electrodynamic Effects of Inflationary Gravitons,” Class. Quant. Grav. 31 (2014), 175002 [arXiv:1308.3453 [gr-qc]].
  * (21) C. L. Wang and R. P. Woodard, “Excitation of Photons by Inflationary Gravitons,” Phys. Rev. D 91 (2015) no.12, 124054 [arXiv:1408.1448 [gr-qc]].
  * (22) D. Glavan, S. P. Miao, T. Prokopec and R. P. Woodard, “Graviton Loop Corrections to Vacuum Polarization in de Sitter in a General Covariant Gauge,” Class. Quant. Grav. 32 (2015) no.19, 195014 [arXiv:1504.00894 [gr-qc]].
  * (23) C. L. Wang and R. P. Woodard, “One-loop quantum electrodynamic correction to the gravitational potentials on de Sitter spacetime,” Phys. Rev. D 92 (2015), 084008 [arXiv:1508.01564 [gr-qc]].
  * (24) D. Glavan, S. P. Miao, T. Prokopec and R. P. Woodard, “One loop graviton corrections to dynamical photons in de Sitter,” Class. Quant. Grav. 34 (2017) no.8, 085002 [arXiv:1609.00386 [gr-qc]].
  * (25) S. P. Miao, T. Prokopec and R. P. Woodard, “Scalar enhancement of the photon electric field by the tail of the graviton propagator,” Phys. Rev. D 98 (2018) no.2, 025022 [arXiv:1806.00742 [gr-qc]].
  * (26) J. Beltran Jimenez and A. L. Maroto, “Cosmological electromagnetic fields and dark energy,” JCAP 03 (2009), 016 [arXiv:0811.0566 [astro-ph]].
  * (27) R. P. Woodard, “de Sitter breaking in field theory,” published in “Deserfest: A Celebration of the Life and Works of Stanley Deser” (Eds. J. T. Liu, M. J. Duff, K. S. Stelle, R. P. Woodard), World Scientific (2006), [arXiv:gr-qc/0408002 [gr-qc]].
  * (28) B. Allen and T. Jacobson, “Vector Two Point Functions in Maximally Symmetric Spaces,” Commun. Math. Phys. 103 (1986), 669
  * (29) N. C. Tsamis and R. P. Woodard, “A Maximally symmetric vector propagator,” J. Math. Phys. 48 (2007), 052306 [arXiv:gr-qc/0608069 [gr-qc]].
  * (30) A. Youssef, “Infrared behavior and gauge artifacts in de Sitter spacetime: The photon field,” Phys. Rev. Lett. 107 (2011), 021101 [arXiv:1011.3755 [gr-qc]].
  * (31) M. B. Fröb and A. Higuchi, “Mode-sum construction of the two-point functions for the Stueckelberg vector fields in the Poincaré patch of de Sitter space,” J. Math. Phys. 55 (2014), 062301 [arXiv:1305.3421 [gr-qc]].
  * (32) D. Glavan, A. Marunović, T. Prokopec and Z. Zahraee, “Abelian Higgs model in power-law inflation: the propagators in the unitary gauge,” JHEP 09 (2020), 165 [arXiv:2005.05435 [gr-qc]].
  * (33) E. Huguet and J. Renaud, “Two-point function for the Maxwell field in flat Robertson-Walker spacetimes,” Phys. Rev. D 88 (2013) no.12, 124018 [arXiv:1310.7333 [hep-th]].
  * (34) I. I. Cotaescu and C. Crucean, “The Quantum theory of the free Maxwell field on the de Sitter expanding universe,” Prog. Theor. Phys. 124 (2010), 1051-1066 [arXiv:0806.2515 [gr-qc]].
  * (35) S. N. Gupta, “Theory of longitudinal photons in quantum electrodynamics,” Proc. Phys. Soc. A 63 (1950) 681.
  * (36) K. Bleuler, “A New method of treatment of the longitudinal and scalar photons,” Helv. Phys. Acta 23 (1950) 567.
  * (37) F. Finster and A. Strohmaier, “Gupta–Bleuler Quantization of the Maxwell Field in Globally Hyperbolic Space-Times,” Annales Henri Poincare 16 (2015) no.8, 1837-1868 [arXiv:1307.1632 [math-ph]].
  * (38) Y. Zhang and X. Ye, “Maxwell field with gauge fixing term in de Sitter space: Exact solution and stress tensor,” Phys. Rev. D 106 (2022) no.6, 065004 [arXiv:2208.08057 [gr-qc]].
  * (39) S. P. Miao, N. C. Tsamis and R. P. Woodard, “Transforming to Lorentz Gauge on de Sitter,” J. Math. Phys. 50 (2009) 122502 [arXiv:0907.4930 [gr-qc]].
  * (40) D. Glavan and T. Prokopec, “Photon propagator in de Sitter space in the general covariant gauge,” JHEP 05 (2023), 126 [arXiv:2212.13982 [gr-qc]].
  * (41) D. Glavan and T. Prokopec, “Even the photon propagator must break de Sitter symmetry,” Phys. Lett. B 841 (2023), 137928 [arXiv:2212.13997 [hep-th]].
  * (42) N. A. Chernikov and E. A. Tagirov, “Quantum theory of scalar fields in de Sitter space-time,” Ann. Inst. H. Poincare Phys. Theor. A 9 (1968), 109.
  * (43) T. S. Bunch and P. C. W. Davies, “Quantum Field Theory in de Sitter Space: Renormalization by Point Splitting,” Proc. Roy. Soc. Lond. A 360 (1978), 117-134.
  * (44) J. S. Schwinger, “Brownian motion of a quantum oscillator,” J. Math. Phys. 2 (1961), 407-432
  * (45) K. T. Mahanthappa, “Multiple production of photons in quantum electrodynamics,” Phys. Rev. 126 (1962), 329-340
  * (46) P. M. Bakshi and K. T. Mahanthappa, “Expectation value formalism in quantum field theory. 1.,” J. Math. Phys. 4 (1963), 1-11
  * (47) P. M. Bakshi and K. T. Mahanthappa, “Expectation value formalism in quantum field theory. 2.,” J. Math. Phys. 4 (1963), 12-16
  * (48) L. V. Keldysh, “Diagram technique for nonequilibrium processes,” Zh. Eksp. Teor. Fiz. 47 (1964), 1515-1527
  * (49) K. c. Chou, Z. b. Su, B. l. Hao and L. Yu, “Equilibrium and Nonequilibrium Formalisms Made Unified,” Phys. Rept. 118 (1985), 1-131
  * (50) E. Calzetta and B. L. Hu, “Closed Time Path Functional Formalism in Curved Space-Time: Application to Cosmological Back Reaction Problems,” Phys. Rev. D 35 (1987), 495
  * (51) R. D. Jordan, “Effective Field Equations for Expectation Values,” Phys. Rev. D 33 (1986), 444-454
  * (52) M. Henneaux and C. Teitelboim, “Quantization of Gauge Systems,” (Princeton University Press, Princeton, New Jersey, USA, 1994)
  * (53) P. A. M. Dirac, “Lectures on Quantum Mechanics,” (1964. Reprint, Dover Publications, Mineola, New York, 2011).
  * (54) D. N. Vollick, “Lorentz gauge quantization in a cosmological space-time,” Phys. Rev. D 86 (2012), 084057 [arXiv:1301.0032 [gr-qc]].
  * (55) J. C. Cresswell and D. N. Vollick, “Lorenz gauge quantization in conformally flat spacetimes,” Phys. Rev. D 91 (2015) no.8, 084008 [arXiv:1504.05914 [gr-qc]].
  * (56) J. Beltran Jimenez and A. L. Maroto, “The electromagnetic dark sector,” Phys. Lett. B 686 (2010), 175-180 [arXiv:0903.4672 [astro-ph.CO]].
  * (57) J. Beltran Jimenez and A. L. Maroto, “The Dark Magnetism of the Universe,” Mod. Phys. Lett. A 26 (2011), 3025-3039 [arXiv:1112.1106 [astro-ph.CO]].
  * (58) A. Belokogne and A. Folacci, “Stueckelberg massive electromagnetism in curved spacetime: Hadamard renormalization of the stress-energy tensor and the Casimir effect,” Phys. Rev. D 93 (2016) no.4, 044063 [arXiv:1512.06326 [gr-qc]].
  * (59) A. Belokogne, A. Folacci and J. Queva, “Stueckelberg massive electromagnetism in de Sitter and anti–de Sitter spacetimes: Two-point functions and renormalized stress-energy tensors,” Phys. Rev. D 94 (2016) no.10, 105028 [arXiv:1610.00244 [gr-qc]].
  * (60) S. P. Miao, T. Prokopec and R. P. Woodard, “Deducing Cosmological Observables from the S-matrix,” Phys. Rev. D 96 (2017) no.10, 104029 [arXiv:1708.06239 [gr-qc]].
  * (61) S. Domazet, D. Glavan and T. Prokopec, in progress (2023).