# Vacuum birefringence and dichroism in a strong plane-wave background

I. A. Aleksandrov  Department of Physics, Saint Petersburg State University,
Universitetskaya Naberezhnaya 7/9, Saint Petersburg 199034, Russia  Ioffe
Institute, Politekhnicheskaya street 26, Saint Petersburg 194021, Russia  V.
M. Shabaev  Department of Physics, Saint Petersburg State University,
Universitetskaya Naberezhnaya 7/9, Saint Petersburg 199034, Russia  National
Research Centre “Kurchatov Institute” B.P. Konstantinov Petersburg Nuclear
Physics Institute, Gatchina, Leningrad district 188300, Russia

###### Abstract

In the present study, we consider the effects of vacuum birefringence and
dichroism in strong electromagnetic fields. According to quantum
electrodynamics, the vacuum state exhibits different refractive properties
depending on the probe photon polarization and one also obtains different
probabilities of the photon decay via production of electron-positron pairs.
Here we investigate these two phenomena by means of several different
approaches to computing the polarization operator. The external field is
assumed to be a linearly polarized plane electromagnetic wave of arbitrary
amplitude and frequency. Varying the probe-photon energy and the field
parameters, we thoroughly examine the validity of the locally-constant field
approximation (LCFA) and techniques involving perturbative expansions in terms
of the external-field amplitude. Within the latter approach, we develop a
numerical method based on a direct evaluation of the weak-field Feynman
diagrams, which can be employed for investigating more complex external
backgrounds. It is demonstrated that the polarization operator depends on two
parameters: classical nonlinearity parameter $\xi$ and the product
$\eta=\omega q_{0}/m^{2}$ of the laser field frequency $\omega$ and the photon
energy $q_{0}$ ($m$ is the electron mass). The domains of validity of the
approximate techniques in the $\xi\eta$ plane are explicitly identified.

##  I Introduction

According to quantum electrodynamics (QED), the physical vacuum state contains
quantum fluctuations of the electromagnetic and electron-positron fields,
which can be viewed as spontaneous creation and annihilation of electron-
positron pairs interacting with each other via virtual photons. Although these
virtual particles are not observable themselves, their existence can manifest
itself while interacting with external fields and real particles giving rise
to a number of remarkable nonlinear phenomena such as light-by-light
scattering euler_kockel ; heisenberg_euler ; weisskopf ; karplus_pr_1950_1951
, Sauter-Schwinger effect heisenberg_euler ; sauter_1931 ; schwinger_1951 ,
and so on (for review, see, e.g., Refs. dipiazza_rmp_2012 ; xie_review_2017 ;
fedotov_review ). In this investigation, we consider propagation of a probe
photon in vacuum in the presence of a strong external background. The latter
polarizes the physical vacuum, so the probe photon effectively interacts with
a nonlinear medium, which leads to the phenomena of vacuum birefringence and
dichroism Toll:1952 ; Baier ; Baier2 ; baier_1976 ; becker_1975  which are in
the focus of the present study (we note that the nontrivial properties of the
vacuum state in the presence of real photons give also rise to recently
discussed stimulated photon emission aleksandrov_antonino_2022 ).

Observing these processes in the laboratory represents currently an intriguing
and challenging task. There are mainly two different approaches to probing
vacuum birefringence. First, one can rely on unprecedented accuracy of
experimental measurements in the optical domain, i.e., in the regime of
relatively low probe-photon energies (see, e.g., Refs. di_piazza_prl_2006 ;
heinzl_2006 ; dinu_prd_2014 ; karbstein_zepf_2015 ; felix_mosman_prd_2020 ;
felix_prl_2022 ). From the theoretical viewpoint, this domain allows one to
employ local approximations, i.e., to treat the external (laser) field as a
locally constant background. The corresponding locally-constant field
approximation (LCFA) has basically two different implementations based either
on employing the exact expressions for the Heisenberg-Euler effective
Lagrangian felix_rashid_2015  or on using the local values of the polarization
operator derived in constant crossed fields meuren_2013 ; bragin_2017 . The
second approach to vacuum birefringence involves high-energy probe photons
king_elkina_2016 ; nakamiya_2017 ; bragin_2017 . The advantage of this
technique appears due to large probabilities of the corresponding quantum
processes resulting in large values of the experimental signal. On the other
hand, it is significantly more difficult to perform measurements in the high-
energy domain as, e.g., the Heisenberg-Euler approximation is only valid in
the low-energy domain. To properly assess the feasibility of the corresponding
scenarios, one has to obtain accurate and reliable theoretical predictions.

In order to avoid approximate local treatment of the external electromagnetic
field, one can model it with a plane-wave background allowing one to deduce
explicit analytical expressions for the polarization tensor baier_1976 ;
becker_1975 ; meuren_2013 . On the other hand, this simplified setup may not
properly reflect the properties of real experimental conditions.

In the present study, we have two primary aims. First, we will thoroughly
examine the plane-wave scenario by means of analytical nonperturbative
expressions derived in Refs. baier_1976 ; becker_1975 ; meuren_2013 . We will
compute the polarization tensor in a wide range of physical parameters
governing the process under consideration: laser-field amplitude, laser
frequency, and probe-photon energy. Expanding the nonperturbative result in
powers of the external-field amplitude, we will assess the accuracy of the
calculations based on perturbation theory (PT). Besides, we will
quantitatively analyze the validity of the LCFA in the two forms described
above. Second, the polarization tensor will be directly evaluated via the
corresponding Feynman diagrams. This approach is very important since it can
allow one to consider other field configurations, which differ from a simple
plane-wave scenario. In what follows, we will benchmark our direct
computational procedures and also provide an additional insight into the
analytical properties of the integrands involved in the Feynman diagrams. For
instance, it will be demonstrated that the overlap between the branch cuts
that appears for sufficiently high photon energies is closely related to the
decay of the probe photon via production of electron-positron pairs. We also
mention that $e^{+}e^{-}$ pairs can be produced directly by a classical strong
field, i.e., via the Sauter-Schwinger mechanism. The validity of the LCFA in
this context was recently examined in Refs. aleksandrov_prd_2019_1 ;
sevostyanov_prd_2021 ; aleksandrov_symmetry ; aleksandrov_sevostyanov_2022 .

The paper has the following structure. In Sec. II we describe the setup under
consideration involving a probe photon and external plane-wave background. In
Sec. III we present nonperturbative expressions which we employ in our
numerical computations. In Sec. IV we calculate the leading-order contribution
with respect to the external-field amplitude. Section V is devoted to the
description of the two possible implementations of the LCFA. In Sec. VI we
discuss how one can directly evaluate the leading-order Feynman diagrams.
Section VII contains our numerical results obtained by means of the various
techniques. Finally, we conclude in Sec. VIII.

Throughout the text, we employ the units $\hbar=c=1$,
$\alpha=e^{2}/(4\pi)\approx 1/137$.

##  II Setup and notation

We assume that the external plane wave is polarized along the $x$ axis and
propagates in the $z$ direction, i.e., it depends on $\varphi=\omega
n^{\mu}x_{\mu}=\omega(t-z)$, where $\omega$ is the laser frequency. The null
vector $n$ obeys $n_{0}=1$, $n^{2}=0$. The corresponding vector potential has
the following form:

$\displaystyle\boldsymbol{A}(x)$ $\displaystyle=$
$\displaystyle\mathcal{A}(\omega(t-z))\boldsymbol{e}_{x},$ (1)
$\displaystyle\mathcal{A}(\varphi)$ $\displaystyle=$
$\displaystyle\frac{E_{0}}{\omega}\sin\varphi,$ (2)

where $E_{0}$ is the field strength amplitude. We also introduce a
dimensionless parameter $\xi=|eE_{0}|/(m\omega)$. The initial photon momentum
$\boldsymbol{q}$ points in the opposite direction to
$\boldsymbol{n}=\boldsymbol{e}_{z}$,
$\boldsymbol{q}=-q_{0}\boldsymbol{e}_{z}$. Accordingly, the initial 4-momentum
of the photon is $q^{\mu}=q_{0}(1,0,0,-1)^{\text{t}}$. The final momentum will
be denoted by $k^{\mu}$. In what follows, we will also employ the light-cone
components which for arbitrary 4-vector $v^{\mu}$ read

$\displaystyle v_{+}$ $\displaystyle=$
$\displaystyle\frac{v_{0}+\boldsymbol{n}\boldsymbol{v}}{2},$ (3)
$\displaystyle v_{-}$ $\displaystyle=$ $\displaystyle
v_{0}-\boldsymbol{n}\boldsymbol{v},$ (4) $\displaystyle\boldsymbol{v}_{\perp}$
$\displaystyle=$
$\displaystyle\boldsymbol{v}-(\boldsymbol{n}\boldsymbol{v})\boldsymbol{n}.$
(5)

The scalar product of two vectors can be evaluated via

$vw\equiv
v^{\mu}w_{\mu}=v_{+}w_{-}+v_{-}w_{+}-\boldsymbol{v}_{\perp}\boldsymbol{w}_{\perp}.$
(6)

For instance, $n_{+}=1$, $n_{-}=0$, $\boldsymbol{n}_{\perp}=0$, and
$\varphi=\omega x_{-}$.

The amplitude $\mathcal{S}(q,k)$ of the process described by the diagram in
Fig. 1 involves two photon wavefunctions defined as

$f^{\mu}_{q}(x)=\frac{1}{\sqrt{2q_{0}}}\mathrm{e}^{-iqx}\varepsilon^{\mu}(q),$
(7)

where $\varepsilon^{\mu}(q)$ is the polarization 4-vector. The amplitude can
be represented in the form

$\mathcal{S}(q,k)=\frac{1}{\sqrt{4q_{0}k_{0}}}\varepsilon_{\mu}(q)i\big{[}\Pi_{0}^{\mu\nu}(q,k)+\Pi^{\mu\nu}(q,k)\big{]}\varepsilon^{*}_{\nu}(k).$
(8)

Here $\Pi_{0}^{\mu\nu}(q,k)$ denotes the zero-field contribution to the
polarization operator, which corresponds to the diagram with the free-electron
Green’s functions describing vacuum polarization in the absence of external
fields. This contribution diverges and requires a usual renormalization
procedure. Since this term does not affect the processes of vacuum
birefringence and dichroism, our task is to compute the field-dependent part
$\Pi^{\mu\nu}(q,k)$, which is finite.

Figure 1: Feynman diagram describing the leading-order contribution to the
photon polarization operator. The amplitude of the process is proportional to
the fine-structure constant $\alpha$ and exactly takes into account the
interaction with the classical external background (double lines represent the
dressed Green’s functions).

In what follows, we will evaluate $\Pi^{\mu\nu}(q,k)$ by means of several
different techniques mentioned above. As will be seen below, the polarization
operator involving $\xi$, $\omega$, and $q_{0}$ depends, in fact, only on
$\xi$ and the product $\omega q_{0}$. We will consider $\xi$ and
$\eta\equiv\omega q_{0}/m^{2}$ as two independent dimensionless parameters
governing the processes of vacuum birefringence and dichroism. We will also
introduce the so-called quantum nonlinearity parameter $\chi=2\xi\eta$ which
will be considered as a derived quantity $\chi(\xi,\eta)$.

##  III Nonperturbative analytical formulas

In the case of a plane-wave external background, it is possible to compute the
polarization tensor analytically. In Ref. baier_1976  it was done by means of
the operator approach. In Ref. becker_1975  the calculations were performed in
the case of a monochromatic plane wave. Recently, in Ref. meuren_2013  the
results of Refs. baier_1976 ; becker_1975  were confirmed by direct
computations of the Feynman diagram in Fig. 1 with the aid of the exact
Green’s functions, which can be constructed from the Volkov solutions.

Here we will first employ the general expressions presented in Refs.
baier_1976 ; becker_1975 ; meuren_2013 . Due to the symmetry of the external
plane-wave field, it can only change the $q_{+}$ component of the photon
momentum, so the amplitude corresponding to the Feynman diagram in Fig. 1
contains
$\delta(k_{-}-q_{-})\delta(\boldsymbol{k}_{\perp}-\boldsymbol{q}_{\perp})$. It
turns out that the cumbersome expressions for the amplitude derived in Refs.
baier_1976 ; becker_1975 ; meuren_2013  become relatively simple in the
particular case of a circularly polarized plane-wave background. Due to the
helicity conservation, the momentum component $q_{+}$ can change only by $\pm
2\omega$ or remain the same. It is not the case if the external field has a
linear polarization since such a plane wave does not possess a well-defined
helicity quantum number. Accordingly, the $q_{+}$ momentum component of the
photon may change by an arbitrary integer number of $\omega$. The general
expression for the setup described above has the following form:

$\Pi^{\mu\nu}(q,k)=-\frac{4\pi^{2}\alpha}{\omega}\delta(k_{-}-q_{-})\delta(\boldsymbol{k}_{\perp}-\boldsymbol{q}_{\perp})\int\limits_{-1}^{1}dv\int\limits_{0}^{\infty}\frac{d\tau}{\tau}\int\limits_{-\infty}^{\infty}d\varphi\,\mathrm{e}^{i\Phi}\begin{pmatrix}c&0&0&0\\\
0&b+\Delta b&0&0\\\ 0&0&b&0\\\ 0&0&0&c\end{pmatrix},$ (9)

where

$\displaystyle b$ $\displaystyle=$
$\displaystyle\Big{(}\frac{i}{\tau}+\frac{1}{2}kq\Big{)}(1-\mathrm{e}^{i\tau\beta})+\frac{2m^{2}\tau\xi^{2}}{\mu}\,\mathrm{e}^{i\tau\beta}\sin^{2}(\mu\omega
q_{0})\cos^{2}\varphi,$ (10) $\displaystyle\Delta b$ $\displaystyle=$
$\displaystyle 2m^{2}\xi^{2}\Big{[}\mathrm{sinc}^{2}(\mu\omega
q_{0})\sin^{2}\varphi-2\,\mathrm{sinc}(2\mu\omega
q_{0})\sin^{2}\varphi-\sin^{2}(\mu\omega
q_{0})+\sin^{2}\varphi\Big{]}\mathrm{e}^{i\tau\beta},$ (11) $\displaystyle c$
$\displaystyle=$
$\displaystyle\frac{k_{0}q_{0}\mu}{\tau}\,(1-\mathrm{e}^{i\tau\beta}),$ (12)
$\displaystyle\mu$ $\displaystyle=$ $\displaystyle\frac{1}{2}\tau(1-v^{2}),$
(13) $\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\frac{k_{+}-q_{+}}{\omega}\varphi+\frac{1}{2}\mu kq-m^{2}\tau,$
(14) $\displaystyle\beta$ $\displaystyle=$ $\displaystyle
m^{2}\xi^{2}\bigg{[}\mathrm{sinc}^{2}(\mu\omega
q_{0})\sin^{2}\varphi-\frac{1}{2}+\frac{1}{2}\,\mathrm{sinc}(2\mu\omega
q_{0})\cos 2\varphi\bigg{]}.$ (15)

In what follows, we will be interested only in the elastic process, where
$k_{+}=q_{+}$ as the other channels are significantly suppressed (actually,
they rather represent reactions involving photon merging or splitting than the
phenomenon of birefringence). To extract the particular process of elastic
scattering, one has to isolate the zeroth-order Fourier harmonics with respect
to $\varphi$ dependence in the functions $b$, $\Delta b$, and $c$, so the
integration of $\mathrm{exp}(i\Phi)$ yields the necessary delta-function. This
can be straightforwardly attained with the aid of the Jacobi-Anger identity.
The result reads

$\Pi_{\text{elastic}}^{\mu\nu}(q,k)=-(2\pi)^{3}\alpha\delta(k-q)\int\limits_{-1}^{1}dv\int\limits_{0}^{\infty}\frac{d\tau}{\tau}\,\mathrm{e}^{-im^{2}\tau}\begin{pmatrix}\tilde{c}&0&0&0\\\
0&\tilde{b}+\Delta\tilde{b}&0&0\\\ 0&0&\tilde{b}&0\\\
0&0&0&\tilde{c}\end{pmatrix},$ (17)

where

$\displaystyle\tilde{b}$ $\displaystyle=$ $\displaystyle\frac{i}{\tau}[1-\Xi
J_{0}(A)]+\frac{m^{2}\tau\xi^{2}}{\mu}\sin^{2}(\mu\omega
q_{0})\Xi[J_{0}(A)+iJ_{1}(A)],$ (18) $\displaystyle\Delta\tilde{b}$
$\displaystyle=$ $\displaystyle m^{2}\xi^{2}\Xi\big{\\{}-2\sin^{2}(\mu\omega
q_{0})J_{0}(A)+[\mathrm{sinc}^{2}(\mu\omega q_{0})-2\,\mathrm{sinc}(2\mu\omega
q_{0})+1][J_{0}(A)-iJ_{1}(A)]\big{\\}},$ (19) $\displaystyle\tilde{c}$
$\displaystyle=$ $\displaystyle\frac{q^{2}_{0}\mu}{\tau}\,[1-\Xi J_{0}(A)],$
(20) $\displaystyle\Xi$ $\displaystyle=$
$\displaystyle\mathrm{exp}\bigg{\\{}\frac{i}{2}m^{2}\tau\xi^{2}[\mathrm{sinc}^{2}(\mu\omega
q_{0})-1]\bigg{\\}},$ (21) $\displaystyle A$ $\displaystyle=$
$\displaystyle\frac{1}{2}m^{2}\tau\xi^{2}[\mathrm{sinc}(2\mu\omega
q_{0})-\mathrm{sinc}^{2}(\mu\omega q_{0})].$ (22)

Here $J_{n}$ are the Bessel functions of the first kind. We will assume
hereinafter $k^{\mu}=q^{\mu}$. We also note that the elements $\Pi^{00}$ and
$\Pi^{33}$ are equal, which preserves the gauge invariance and the Ward-
Takahashi identity blp . These components will not be evaluated in our study
as they do not affect the phenomena under consideration.

The birefringent and dichroic properties of the vacuum in the presence of
strong fields manifest themselves in the difference between $\Pi^{11}$ and
$\Pi^{22}$ elements: photon polarizations along the $x$ and $y$ axes
correspond to different refractive and absorption indexes. In what follows, we
will compute these elements. As was stated above, these quantities involve the
three parameters $\xi$, $\omega$, and $q_{0}$, but they depend, in fact, on
$\xi$ and $\eta=\omega q_{0}/m^{2}$ as becomes evident from Eqs. (17)–(22).

Figure 2: Feynman diagrams corresponding to the leading-order contribution
within the PT expansion in terms of the external field (the amplitudes are
proportional to $\xi^{2}$). The interaction with the classical external field
is denoted by the cross. Depending on the energy-momentum transfer at the
cross vertices, the process is either elastic (2-to-2 process) or corresponds
to $k\neq q$.

##  IV Perturbation theory

Here we will consider the leading-order term of Eq. (17) with respect to the
small-$\xi$ expansion. This contribution is proportional to $\xi^{2}$ and
corresponds to the three Feynman diagrams displayed in Fig. 2. Expanding the
function $\Xi$ and the Bessel functions in Taylor series, one obtains

$\displaystyle\tilde{b}_{\text{LO}}$ $\displaystyle=$ $\displaystyle
m^{2}\xi^{2}\bigg{\\{}\frac{1}{2}[\mathrm{sinc}^{2}(\mu\omega
q_{0})-1]+\frac{\tau}{\mu}\,\sin^{2}(\mu\omega q_{0})\bigg{\\}},$ (23)
$\displaystyle\Delta\tilde{b}_{\text{LO}}$ $\displaystyle=$ $\displaystyle
m^{2}\xi^{2}[-2\sin^{2}(\mu\omega q_{0})+\mathrm{sinc}^{2}(\mu\omega q_{0})$
(24) $\displaystyle-$ $\displaystyle 2\,\mathrm{sinc}(2\mu\omega q_{0})+1],$
$\displaystyle\tilde{c}_{\text{LO}}$ $\displaystyle=$
$\displaystyle-\frac{i}{2}q_{0}^{2}\mu
m^{2}\xi^{2}[\mathrm{sinc}^{2}(\mu\omega q_{0})-1].$ (25)

Here “LO” stands for “low order”. It turns out that one can replace $\mu$ with
Eq. (13) and perform the $\tau$ integration analytically. Let us first
introduce the following general representation:

$\Pi_{\text{elastic}}^{\mu\nu}(q,k)=-(2\pi)^{3}\alpha\delta(k-q)m^{2}\xi^{2}M^{\mu\nu}.$
(26)

Within PT we find

$\displaystyle M^{11}_{\text{LO}}$ $\displaystyle=$
$\displaystyle\int\limits_{-1}^{1}dv\bigg{[}\frac{2v^{2}}{1-v^{2}}\,I_{1}(v)+\frac{1}{2}\,I_{2}(v)+I_{3}(v)\bigg{]},$
(27) $\displaystyle M^{22}_{\text{LO}}$ $\displaystyle=$
$\displaystyle\int\limits_{-1}^{1}dv\bigg{[}\frac{2}{1-v^{2}}\,I_{1}(v)+\frac{1}{2}\,I_{2}(v)\bigg{]},$
(28)

where

$\displaystyle I_{1}(v)$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\infty}\frac{dt}{t}\,\sin^{2}(\gamma
t)\mathrm{e}^{-it}=\frac{1}{4}\ln\big{|}1-4\gamma^{2}\big{|}-\frac{i\pi}{4}\theta(\gamma-1/2),$
(29) $\displaystyle I_{2}(v)$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\infty}\frac{dt}{t}\,[\mathrm{sinc}^{2}(\gamma
t)-1]\mathrm{e}^{-it}=\frac{3}{2}-\frac{1}{2}\bigg{(}1+\frac{1}{4\gamma^{2}}\bigg{)}\ln\big{|}1-4\gamma^{2}\big{|}-\frac{1}{2\gamma}\,\ln\bigg{|}\frac{1+2\gamma}{1-2\gamma}\bigg{|}$
(30) $\displaystyle+$
$\displaystyle\frac{i\pi}{2}\bigg{(}1-\frac{1}{\gamma}+\frac{1}{4\gamma^{2}}\bigg{)}\theta(\gamma-1/2),$
$\displaystyle I_{3}(v)$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\infty}\frac{dt}{t}\,[1+\mathrm{sinc}^{2}(\gamma
t)-2\,\mathrm{sinc}(2\gamma
t)]\mathrm{e}^{-it}=-\frac{1}{2}+\frac{1}{2}\bigg{(}1-\frac{1}{4\gamma^{2}}\bigg{)}\Big{[}\ln\big{|}1-4\gamma^{2}\big{|}-i\pi\theta(\gamma-1/2)\Big{]},$
(31) $\displaystyle\gamma$ $\displaystyle=$
$\displaystyle\gamma(v)=\frac{\omega
q_{0}}{2m^{2}}(1-v^{2})=\frac{1}{2}\eta(1-v^{2}).$ (32)

The expressions (27) and (28) depend only on $\eta=\omega q_{0}/m^{2}$, while
the nonperturbative values of $M^{\mu\nu}$ [see Eq. (26)] also involve $\xi$.
Below we will compare the leading-order terms with the nonperturbative
results. Let us now present the low- and high-energy asymptotic expressions
for $M^{11}_{\text{LO}}$ and $M^{22}_{\text{LO}}$. In the low-energy case
$\varepsilon\equiv 2\eta=2\omega q_{0}/m^{2}\ll 1$,

$\displaystyle M^{11}_{\text{LO}}$ $\displaystyle=$
$\displaystyle-\frac{4}{45}\varepsilon^{2}-\frac{17}{3150}\varepsilon^{4}+\mathcal{O}(\varepsilon^{6}),$
(33) $\displaystyle M^{22}_{\text{LO}}$ $\displaystyle=$
$\displaystyle-\frac{7}{45}\varepsilon^{2}-\frac{131}{9450}\varepsilon^{4}+\mathcal{O}(\varepsilon^{6}).$
(34)

In the high-energy limit, we obtain [$\varepsilon\equiv
1/(2\eta)=m^{2}/(2\omega q_{0})\ll 1$]

$\displaystyle M^{11}_{\text{LO}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\ln^{2}\varepsilon+\bigg{(}1-\ln
2+\frac{i\pi}{2}\bigg{)}\ln\varepsilon+\bigg{[}\frac{5}{2}-\ln
2+\frac{1}{2}\ln^{2}2-\frac{\pi^{2}}{4}+\frac{i\pi}{2}(1-\ln 2)\bigg{]}$ (35)
$\displaystyle+$ $\displaystyle
i\pi\varepsilon\ln\varepsilon+\bigg{[}-\frac{\pi^{2}}{2}+\frac{i\pi}{2}(3-2\ln
2)\bigg{]}\varepsilon+\mathcal{O}(\varepsilon^{2}\ln^{2}\varepsilon),$
$\displaystyle M^{22}_{\text{LO}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\ln^{2}\varepsilon+\bigg{(}1-\ln
2+\frac{i\pi}{2}\bigg{)}\ln\varepsilon+\bigg{[}\frac{7}{2}-\ln
2+\frac{1}{2}\ln^{2}2-\frac{\pi^{2}}{4}+\frac{i\pi}{2}(1-\ln 2)\bigg{]}$ (36)
$\displaystyle+$ $\displaystyle
i\pi\varepsilon\ln\varepsilon+\bigg{[}-\frac{\pi^{2}}{2}+\frac{i\pi}{2}(1-2\ln
2)\bigg{]}\varepsilon+\mathcal{O}(\varepsilon^{2}\ln^{2}\varepsilon).$

While the low-energy result (33), (34) is real, the expressions (35) and (36)
possess imaginary parts, which describe the process of photon decay. The
imaginary part of the difference $\delta M_{\text{LO}}\equiv
M^{11}_{\text{LO}}-M^{22}_{\text{LO}}\approx-1+i\pi\varepsilon$ governs the
dichroic properties of the vacuum and appears once $\eta>1$. In Sec. VI we
will discuss how the imaginary part appears in a direct evaluation of the
Feynman diagrams in Fig. 2.

##  V Locally-constant field approximation

Here we will employ relatively simple closed-form expressions treating the
external background as locally constant. There are basically two different
approaches. The first one is based on calculating the polarization tensor in
constant crossed fields and the using the actual spatiotemporal dependence of
the plane-wave field (1) when integrating over $\varphi$. The second method
employs the Heisenberg-Euler effective Lagrangian computed in a constant
electromagnetic field and takes into account the leading-order quantum
correction with respect to the field amplitude $E_{0}$. The first approach is
generally more accurate as it incorporates the higher-order terms in $E_{0}$
and involves the expression for the polarization operator which is derived for
arbitrary photon energies $q_{0}$. The second technique based on the
Heisenberg-Euler Lagrangian is only valid for sufficiently low photon
energies, when there is only a small momentum transfer into the $e^{+}e^{-}$
loop in the diagram in Fig. 1. Besides, the applicability of this method is
limited since it involves the PT expansion with respect to the field
amplitude. In what follows, we will describe the both approaches and then
thoroughly analyze their validity.

###  V.1 Polarization operator in constant crossed fields

In the setup under consideration, the vector potential (1) is assumed to be a
monochromatic plane wave (2). If one replaces $\sin\leavevmode\nobreak\
\varphi$ in Eq. (2) with $\varphi$, the external background will obviously
become a combination of constant crossed electric and magnetic fields,
$E_{x}=B_{y}=-E_{0}$. In this case, one can also perform nonperturbative
calculations of the polarization tensor narozhnyi_28_371_1969 ;
batalin_shabad_1968 ; ritus_ann_phys_1972  and then locally approximate a
generic external background by constant crossed fields meuren_2013 . Applying
this technique to the field configuration (2), one obtains

$\displaystyle M^{11}_{\text{LCFA}}$ $\displaystyle=$
$\displaystyle\frac{1}{3\pi\xi^{2}}\int\limits_{-1}^{1}dv\,\bigg{(}\frac{\chi}{w}\bigg{)}^{2/3}\big{(}w-1\big{)}g(v),$
(37) $\displaystyle M^{22}_{\text{LCFA}}$ $\displaystyle=$
$\displaystyle\frac{1}{3\pi\xi^{2}}\int\limits_{-1}^{1}dv\,\bigg{(}\frac{\chi}{w}\bigg{)}^{2/3}\big{(}w+2\big{)}g(v),$
(38)

where $\chi=2\xi\eta$, $w=4/(1-v^{2})$, and

$\displaystyle g(v)$ $\displaystyle=$
$\displaystyle\int\limits_{-\pi}^{\pi}d\varphi
f^{\prime}(u)(\cos\varphi)^{2/3},$ (39) $\displaystyle u$ $\displaystyle=$
$\displaystyle\bigg{(}\frac{w}{\chi\cos\varphi}\bigg{)}^{2/3},$ (40)
$\displaystyle f(u)$ $\displaystyle=$ $\displaystyle
i\int\limits_{0}^{\infty}d\tau\mathrm{e}^{-i(u\tau+\tau^{3}/3)}=\pi\mathrm{Gi}(u)+i\pi\mathrm{Ai}(u).$
(41)

Here $\mathrm{Gi}$ and $\mathrm{Ai}$ are the Scorer and Airy functions,
respectively.

Note that the integrals in Eqs. (37) and (38) depend only on $\chi$, i.e. the
product $\xi\eta$, which simplifies the further analysis. This fact is a well-
known property of the LCFA ritus_1985 . This approximation is well justified
if the parameter $\xi$ is sufficiently large, so one can expect that the
predictions (37) and (38) significantly differ from the exact nonperturbative
result given in Eq. (17) once $\xi\lesssim 1$. This issue will be discussed in
detail in Sec. VII.

Finally, we present the asymptotic forms of Eqs. (37) and (38) in the case
$\chi\ll 1$. One obtains

$\displaystyle\mathrm{Re}\,M^{11}_{\text{LCFA}}$ $\displaystyle=$
$\displaystyle-\frac{4\chi^{2}}{45\xi^{2}}\bigg{[}1+\frac{1}{4}\,\chi^{2}+\mathcal{O}(\chi^{4})\bigg{]},$
(42) $\displaystyle\mathrm{Re}\,M^{22}_{\text{LCFA}}$ $\displaystyle=$
$\displaystyle-\frac{7\chi^{2}}{45\xi^{2}}\bigg{[}1+\frac{13}{49}\,\chi^{2}+\mathcal{O}(\chi^{4})\bigg{]},$
(43) $\displaystyle\mathrm{Im}\,M^{11}_{\text{LCFA}}$ $\displaystyle=$
$\displaystyle-\frac{3\chi^{3/2}}{8\xi^{2}}\sqrt{\frac{\pi}{2}}\,\mathrm{e}^{-8/(3\chi)}\big{[}1+\mathcal{O}(\chi)\big{]},$
(44) $\displaystyle\mathrm{Im}\,M^{22}_{\text{LCFA}}$ $\displaystyle=$
$\displaystyle-\frac{3\chi^{3/2}}{4\xi^{2}}\sqrt{\frac{\pi}{2}}\,\mathrm{e}^{-8/(3\chi)}\big{[}1+\mathcal{O}(\chi)\big{]}.$
(45)

For small $\chi$ the imaginary part is exponentially suppressed corresponding
to tiny probabilities of the photon decay. Note that the ratio $\chi/\xi$
coincides with $\varepsilon=2\eta$ in Eqs. (33) and (34), so the leading-order
contribution is reproduced by the LCFA. Nevertheless, the validity of the LCFA
and that of the PT expansion correspond to substantially different domains of
parameters. Whereas for given $\xi$ they both are accurate for sufficiently
small $\eta<\eta_{\text{max}}(\xi)$, with increasing $\xi$ the bound
$\eta_{\text{max}}(\xi)$ increases in the case of the LCFA and decreases in
the case of PT. This will be quantitatively demonstrated in Sec. VII. Finally,
we note that both the LCFA and PT capture the imaginary part of the
polarization tensor.

###  V.2 Heisenberg-Euler approximation

Another approach is based on the PT expansion of the polarization operator
derived from the one-loop effective Lagrangian in the presence of a constant
electromagnetic background felix_rashid_2015 . The approximate formula for the
$\xi^{2}$ contribution to the polarization tensor has the following form:

$\Pi^{\mu\nu}_{\text{LCFA-
HE}}(q,k)=\frac{\alpha}{45\pi}\frac{e^{2}}{m^{4}}\int\\!d^{4}x\,\mathrm{e}^{i(k-q)x}\Big{[}4(qF)^{\mu}(kF)^{\nu}+7(qG)^{\mu}(kG)^{\nu}\Big{]}.$
(46)

Here $(kF)^{\mu}\equiv k_{\rho}F^{\rho\mu}$. The electromagnetic tensor
$F_{\mu\nu}=\partial_{\mu}\mathcal{A}_{\nu}-\partial_{\nu}\mathcal{A}_{\mu}$
and the dual tensor
$G^{\mu\nu}=(1/2)\varepsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}$ are evaluated
at the spacetime point $x$ according to the local treatment of the external
field. In the case of the plane-wave background (1), the integrals in Eq. (46)
lead to the conservation laws which may change the photon momentum by $\pm
2\omega$ or keep it the same. We are interested in the latter contribution
governing the elastic process. The explicit form of Eq. (46) then reads

$\Pi_{\text{LCFA-HE,
elastic}}^{\mu\nu}(q,k)=\frac{32\pi^{3}\alpha}{45}m^{2}\xi^{2}\bigg{(}\frac{\omega
q_{0}}{m^{2}}\bigg{)}^{2}\delta(k-q)\begin{pmatrix}0&0&0&0\\\ 0&4&0&0\\\
0&0&7&0\\\ 0&0&0&0\end{pmatrix}.$ (47)

This exactly corresponds to the leading low-energy terms in Eqs. (33) and (34)
and to the leading-order terms in Eqs. (42) and (43). In what follows, they
will be denoted by $M^{11}_{\text{LCFA-LO}}$ and $M^{22}_{\text{LCFA-LO}}$,
respectively. Note that the leading-order LCFA expressions completely
disregard the imaginary part of the polarization tensor, i.e., fail to
describe the process of dichroism.

Figure 3: Branch cuts (red) of the electron propagators in the case $q_{0}<m$
before the $z$ integration in Eq. (50) and a possible integration contour
(blue).

Figure 4: Real and imaginary parts of the difference $\delta M\equiv
M^{11}-M^{22}$ calculated within the leading-order of perturbation theory
[Eqs. (27) and (28)], by means of the Heisenberg-Euler approximation (47) and
according to the low-energy expansions (33) and (34). The latter two
approaches yield zero imaginary part.

##  VI Direct evaluation of the Feynman diagrams

Here we will directly compute the Feynman diagrams depicted in Fig. 2. The
corresponding amplitudes and accordingly the contributions to the polarization
tensor are proportional to $E_{0}^{2}$, i.e. $\xi^{2}$ [cf. Eq. (26)]. Each
interaction vertex involves the energy-momentum transfer with the four-vector
$\pm K$, where $K^{\mu}\equiv\omega n^{\mu}$ is the four-momentum of the
photons that constitute the external plane wave. As we are interested in
studying the elastic contributions, the two vertices in each diagram should
correspond to one emission and one absorption, so the diagram represents
essentially a two-to-two scattering process. Since one has to evaluate three
diagrams, the leading-order matrix $M_{\text{LO}}^{\mu\nu}$ is a sum of three
terms, $M_{\text{LO}}^{\mu\nu}=M_{1}^{\mu\nu}+M_{2}^{\mu\nu}+M_{3}^{\mu\nu}$.
Considering, for instance, the first diagram and using the Feynman rules, we
obtain the following expression for $M_{1}^{\mu\nu}$:

$M_{1}^{\mu\nu}=-\frac{i}{8\pi^{2}}\sum_{s=\pm 1}\int
d^{4}p\mathrm{Tr}\,\big{[}\gamma^{\nu}S(p+q/2-sK/2)\gamma^{1}S(p+q/2+sK/2)\gamma^{\mu}S(p-q/2+sK/2)\gamma^{1}S(p-q/2-sK/2)\big{]}.$
(48)

Here $s$ indicates at which of the two vertices the external-field photon is
emitted (absorbed). The integration variables $p^{\mu}$ are shifted, so that
the integrand has a more symmetric form (cf. Ref. OS ). The electron
propagator is given by

$S(p)=\frac{\gamma^{\mu}p_{\mu}+m}{m^{2}-p^{2}-i\varepsilon},$ (49)

where $\varepsilon\to 0$.

One can explicitly verify that the total expression for
$M_{\text{LO}}^{\mu\nu}$ depends only on the product $\omega q_{0}$, i.e.
$\eta=\omega q_{0}/m^{2}$, in accordance with Eqs. (27) and (28). Therefore,
we will assume that $q_{0}=\omega=\sqrt{\eta}m$, so
$\boldsymbol{K}=-\boldsymbol{q}$. Then Eq. (48) takes the form

$\displaystyle M_{1}^{\mu\nu}$ $\displaystyle=$
$\displaystyle-\frac{i}{8\pi^{2}}\int\limits_{-\infty}^{\infty}dz\int
d^{3}\boldsymbol{p}\,\mathrm{Tr}\,\big{[}\gamma^{\nu}S(z,\boldsymbol{p}+\boldsymbol{q})\gamma^{1}S(z+q_{0},\boldsymbol{p})\gamma^{\mu}S(z,\boldsymbol{p}-\boldsymbol{q})\gamma^{1}S(z-q_{0},\boldsymbol{p})$
(50) $\displaystyle+$
$\displaystyle\gamma^{\nu}S(z+q_{0},\boldsymbol{p})\gamma^{1}S(z,\boldsymbol{p}+\boldsymbol{q})\gamma^{\mu}S(z-q_{0},\boldsymbol{p})\gamma^{1}S(z,\boldsymbol{p}-\boldsymbol{q})\big{]}.$

The trace contains denominators that for each $\boldsymbol{p}$ turn to zero at
complex points $z$ with small nonzero imaginary parts for nonzero values of
$\varepsilon$. After the $\boldsymbol{p}$ integration, the trace as a function
of $z$ possesses six branch cuts depicted in Fig. 3 for $q_{0}<m$. The $z$
integration over the real axis in Eq. (50) can be, in fact, performed over any
contour like that displayed in Fig. 3, provided it does not intersect any of
the branch cuts. In the case $q_{0}<m$ ($\eta<1$), one can, for instance,
rotate the contour, so that it coincides with the imaginary axis. Substituting
then $z=iw$, where $w\in\mathbb{R}$, one can explicitly demonstrate that the
total contribution
$M_{\text{LO}}^{\mu\nu}=M_{1}^{\mu\nu}+M_{2}^{\mu\nu}+M_{3}^{\mu\nu}$ is real
in accordance with Eqs. (27) and (28).

In order to address the high-energy case $\eta>1$, we employ the following
numerical procedure. We change the order of the $z$ and $\boldsymbol{p}$
integrations and first integrate over $z\in\mathbb{R}$. Accordingly, the $z$
integrand has a number of isolated poles $\xi_{j}-i\sigma\varepsilon$ where
$\sigma=\pm 1$ and the real parts $\xi_{j}$ depend on $\boldsymbol{p}$. In
each vicinity $(\xi_{j}-\delta,\xi_{j}+\delta)$ we perform the integration
semi-analytically by means of the Sokhotski–Plemelj identity. This allows us
to set $\varepsilon=0$ while performing the rest integrations numerically and
avoid computational singularities.

Our procedure was also generalized to compute the diagrams for arbitrary
independent $q_{0}$ and $\omega$. The main steps here are generally the same.
After that, we confirmed the results obtained by means of the technique
described above. Finally, we note that the expression (50) has a similar form
to the amplitude of photon emission via the so-called tadpole diagram (see
Ref. aleksandrov_prd_2019_2 , where it was evaluated in the regime $\eta<1$).

Figure 5: Real and imaginary parts of the difference $\delta M\equiv
M^{11}-M^{22}$ evaluated within the leading-order of perturbation theory (LO),
by means of the LCFA [Eqs. (37) and (38)] and exact nonperturbative expression
(17) for $\xi=0.1$ (top), $\xi=1.0$ (middle), $\xi=10.0$ (bottom). For
$\xi=0.1$ the “LO” and exact curves coincide.

##  VII Numerical results

We will now perform numerical calculations of the difference $\delta M\equiv
M^{11}-M^{22}$, whose real and imaginary parts govern the effects of vacuum
birefringence and dichroism, respectively. First, we will evaluate $\delta M$
within the leading order with respect to the field amplitude. In this case,
the results do not depend on $\xi$. In Fig. 4 we present $\delta M$ as a
function of $\eta$. First, one observes that the Heisenberg-Euler
approximation within the leading order of perturbation theory can be accurate
only in the low-energy regime. If one takes into account the $1/\eta^{4}$
terms according to Eqs. (33) and (34), the results become slightly more
accurate although they completely fail to reproduce the full PT results for
$\eta>1$. Second, the more general expressions (27) and (28) yield a nonzero
imaginary part for $\eta>1$, so the PT approach may allow one to describe the
effects of dichroism. Finally, we note that our approach based on direct
computations of the Feynman diagrams as described in Sec. VI provides exactly
the same results as Eqs. (27) and (28), which benchmarks the corresponding
numerical procedures. To judge whether the leading-order approximation is
justified, one has to perform nonperturbative calculations for various values
of $\xi$, which will be done next.

Figure 6: Domains of the validity of the various approximate methods: LCFA
based on using the polarization tensor in constant crossed fields (LCFA),
Heisenberg-Euler approximation (LCFA-HE), and PT calculations within the
leading order in terms of the external-field amplitude (LO).

In Fig. 5 we display the real and imaginary parts of $\delta M$ as a function
of $\eta$ for three different values of $\xi$: 0.1, 1.0, and 10.0. We refer to
Eq. (17) as the exact result. First, we observe that the $\eta$ dependence
very nontrivially changes as a function of $\xi$, which cannot be taken into
account by means of the PT approach. Whereas for $\xi\ll 1$, this
approximation provides indeed very accurate results within a broad range of
$\eta$, for $\xi\gtrsim 1$, it fails to reproduce the exact values unless
$\eta\ll 1$. Second, as was mentioned above, the LCFA predictions have the
form $\delta M_{\text{LCFA}}(\xi,\eta)=(1/\xi^{2})\delta
M_{\text{LCFA}}(1,\xi\eta)$, so the different LCFA curves can be obtained by
simply rescaling the plot axes. This approach does not allow one to describe
the nontrivial structure that takes place for $\xi\lesssim 1$, although it is
accurate for very small $\eta$, where the expansions (42) and (43) are valid.

Let us now quantitatively identify the domains of validity of various
approximations for describing the vacuum birefringence effects. In Fig. 6 we
identify the values of $\xi$ and $\eta$ for which the approximate predictions
match the exact results with a relative uncertainty on the level of $10\%$.
First, let us discuss the PT approach, which yields the leading-order
estimates (27) and (28). In the regime $\xi\gg 1$, it is only valid for
$\eta\ll 1$. It turns out that in the corresponding domain of parameters
$\chi\lesssim 0.5$. Since for large values of $\xi$ one can employ the LCFA,
it is possible to estimate the exact result for the real part of $M^{\mu\nu}$
by means of Eqs. (42) and (43). Comparing these with the low-energy asymptotic
expansions (33) and (34), one can obtain the threshold value of $\chi$. For
instance, requiring that the relative uncertainty of PT be less than $10\%$,
one obtains $\chi<\sqrt{(7/2)0.1}\approx 0.59$. According to our numerical
analysis, this condition, in fact, reads $\chi<0.55$. In the regime
$\xi\lesssim 1$, the validity of the LCFA (37), (38) is very limited, so one
has to directly compare the leading-order PT results with the nonperturbative
predictions. In this domain, the applicability of perturbation theory is not
solely governed by $\chi$ as can be seen in Fig. 6, where the domain of the PT
applicability is no longer bounded by a straight line. Finally, we note that
in the region $\xi\lesssim 1$, even if the PT approach fails to reproduce the
exact results for $\eta\sim 1$, it may provide quite accurate predictions for
sufficiently large values of $\eta$, where $\mathrm{Re}\leavevmode\nobreak\
\delta M$ becomes close to $-1$ [see Fig. 5 (middle)]. Moreover, in this
region the nonzero imaginary part of the polarization operator can also be
obtained by means of perturbation theory.

In order to identify the validity domain of the leading-order Heisenberg-Euler
approximation (47), it is sufficient to compare its predictions with the
leading-order PT result (27), (28). Since within these approaches the matrix
$M^{\mu\nu}$ is independent of $\xi$, one should only determine the threshold
value of $\eta$. For the $10\%$ uncertainty level, it amounts to
$\eta_{\text{max}}\approx 0.44$. The validity domain of the Heisenberg-Euler
approximation is then the intersection of the region $\eta<0.44$ and the
validity domain of the PT approach.

The applicability of the LCFA (37), (38) corresponds to a much larger region
than that where the Heisenberg-Euler approximation is justified. It not only
describes the effect of birefringence in the low-energy domain but is also
valid in the case of high-energy probe photons ($\eta\gtrsim 1$), provided
$\xi\gg 1$.

As was indicated above, the imaginary part of the polarization tensor, which
is responsible for dichroic properties of the vacuum, cannot be estimated by
means of the leading-order Heisenberg-Euler approximation (47). Nevertheless,
both the PT approach and the LCFA (37), (38) are very useful here — they can
be employed within the corresponding regions indicated in Fig. 6.

According to our results, the validity domain of the Heisenberg-Euler
approximation is the smallest. The corresponding results can always be
additionally confirmed by either perturbation theory or the LCFA based on the
calculation of the polarization operator in constant crossed fields. The
advantage of the latter approach is the possibility to consider larger values
of $\eta$ once $\xi\gtrsim 1$. Note also that a considerable part of the plot
in Fig. 6 relates to large values of the parameter $\chi$, which are not
realistic at present. Nevertheless, given the logarithmic scale in the graph,
the LCFA covers a domain of parameters which is substantially broader than the
validity region of the Heisenberg-Euler approximation. The PT approach is
always accurate once the LCFA-HE technique is justified. In addition, the
leading-order predictions coincide with the exact results for any values of
$\eta$ if $\xi$ is sufficiently small.

##  VIII Conclusion

In the present investigation, we examined the effects of vacuum birefringence
and dichroism in strong plane-wave backgrounds by means of several theoretical
methods allowing one to evaluate the leading one-loop contribution to the
polarization operator. First, we employed closed-form expressions exactly
incorporating the interaction between the electron-positron field and
classical external background depending on the spatiotemporal coordinates.
Second, we performed calculations within the leading order with respect to the
field amplitude, i.e., by means of perturbation theory. This was done by
expanding the nonperturbative result and by means of our numerical method
based on a direct evaluation of the leading-order Feynman diagrams. It was
found that these two approaches yield identical quantitative predictions both
for real and imaginary parts of the polarization tensor. Varying the field
parameters and the probe-photon energy, we examined the validity of the
perturbative methods. Third, we utilized the locally-constant field
approximation (LCFA) in two different forms: Heisenberg-Euler approximation
and the technique involving exact expressions for the polarization operator in
constant crossed fields. By comparing the approximate predictions with the
exact results, we evidently identified the field and probe-photon parameters
for which each of the approximate techniques is justified.

An important prospect for future studies is the analogous analysis beyond the
plane-wave scenario, where the exact analytical expressions are unknown. In
this case, for instance, the applicability of the LCFA may be additionally
limited if the external electric and magnetic fields are not crossed in
contrast to the field configuration examined in the present investigation.

###### Acknowledgements.

The study was funded by RFBR and ROSATOM, project No. 20-21-00098. I.A.A. also
acknowledges the support from the Foundation for the advancement of
theoretical physics and mathematics “BASIS”.

## References

  * (1) H. Euler and B. Kockel, Naturwiss. 23, 246 (1935).
  * (2) W. Heisenberg and H. Euler, Z. Phys. 98, 714 (1936).
  * (3) V. Weisskopf, Kong. Dans. Vid. Selsk., Mat.-fys. Medd. XIV, 6 (1936).
  * (4) R. Karplus and M. Neuman, Phys. Rev. 80, 380 (1950); 83, 776 (1951).
  * (5) F. Sauter, Z. Phys. 69, 742 (1931).
  * (6) J. Schwinger, Phys. Rev. 82, 664 (1951).
  * (7) A. Di Piazza, C. Müller, K. Z. Hatsagortsyan, and C. H. Keitel, Rev. Mod. Phys. 84, 1177 (2012).
  * (8) B. S. Xie, Z. L. Li, and S. Tang, Matter Radiat. Extremes 2, 225 (2017).
  * (9) A. Fedotov, A. Ilderton, F. Karbstein, B. King, D. Seipt, H. Taya, and G. Torgrimsson, Phys. Rep. 1010, 1 (2023).
  * (10) J. S. Toll, Ph.D. thesis, Princeton Univ., 1952 (unpublished).
  * (11) R. Baier and P. Breitenlohner, Acta Phys. Austriaca 25, 212 (1967).
  * (12) R. Baier and P. Breitenlohner, Nuovo Cimento B 47, 117 (1967).
  * (13) V. N. Baier, A. I. Milstein, and V. M. Strakhovenko, Zh. Eksp. Teor. Fiz. 69, 1893 (1975) [Sov. Phys. JETP 42, 961 (1976)].
  * (14) W. Becker and H. Mitter, J. Phys. A 8, 1638 (1975).
  * (15) I. A. Aleksandrov, A. Di Piazza, G. Plunien, and V. M. Shabaev, Phys. Rev. D 105, 116005 (2022).
  * (16) A. Di Piazza, K. Z. Hatsagortsyan, and C. H. Keitel, Phys. Rev. Lett. 97, 083603 (2006).
  * (17) T. Heinzl, B. Liesfeld, K. U. Amthor, H. Schwoerer, R. Sauerbrey, and A. Wipf, Opt. Commun. 267, 318 (2006).
  * (18) V. Dinu, T. Heinzl, A. Ilderton, M. Marklund, and G. Torgrimsson, Phys. Rev. D. 89, 125003 (2014).
  * (19) F. Karbstein, H. Gies, M. Reuter, and M. Zepf, Phys. Rev. D. 92, 071301 (2015).
  * (20) F. Karbstein and E. A. Mosman, Phys. Rev. D 101, 113002 (2020).
  * (21) F. Karbstein, D. Ullmann, E. A. Mosman, and M. Zepf, Phys. Rev. Lett. 129, 061802 (2022).
  * (22) F. Karbstein and R. Shaisultanov, Phys. Rev. D 91, 085027 (2015).
  * (23) S. Meuren, C. H. Keitel, and A. Di Piazza, Phys. Rev. D 88, 013007 (2013).
  * (24) S. Bragin, S. Meuren, C. H. Keitel, and A. Di Piazza, Phys. Rev. Lett. 119, 250403 (2017).
  * (25) B. King and N. Elkina, Phys. Rev. A 94, 062102 (2016).
  * (26) Y. Nakamiya and K. Homma, Phys. Rev. D 96, 053002 (2017).
  * (27) I. A. Aleksandrov, G. Plunien, and V. M. Shabaev, Phys. Rev. D 99, 016020 (2019).
  * (28) D. G. Sevostyanov, I. A. Aleksandrov, G. Plunien, and V. M. Shabaev, Phys. Rev. D 104, 076014 (2021).
  * (29) I. A. Aleksandrov, D. G. Sevostyanov, and V. M. Shabaev, Symmetry 14, 2444 (2022).
  * (30) I. A. Aleksandrov, D. G. Sevostyanov, and V. M. Shabaev, arXiv:2210.15626.
  * (31) V. B. Berestetskii, E. M. Lifshitz, and L. P. Pitaevskii, Quantum Electrodynamics (Elsevier Butterworth-Heinemann, Oxford, 1982).
  * (32) I. A. Batalin and A. E. Shabad. FIAN Preprint 166, (1968).
  * (33) N. B. Narozhnyi, Sov. Phys. JETP 28, 371 (1969).
  * (34) V. I. Ritus, Ann. Phys. 69, 555 (1972).
  * (35) V. I. Ritus, J. Sov. Laser Res. 6, 497 (1985).
  * (36) I. A. Aleksandrov and V. M. Shabaev, Optics and Spectroscopy 129, 890 (2021).
  * (37) I. A. Aleksandrov, G. Plunien, and V. M. Shabaev, Phys. Rev. D 100, 116003 (2019).


0 Research Research algorithm/technique Max Franke and Steffen Koch are with
the University of Stuttgart, Germany. E-mail<EMAIL_ADDRESS>stuttgart.de

# Compact Phase Histograms for Guided Exploration of Periodicity

Max Franke  Steffen Koch

###### Abstract

Periodically occurring accumulations of events or measured values are present
in many time-dependent datasets and can be of interest for analyses. The
frequency of such periodic behavior is often not known in advance, making it
difficult to detect and tedious to explore. Automated analysis methods exist,
but can be too costly for smooth, interactive analysis. We propose a compact
visual representation that reveals periodicity by showing a phase histogram
for a given period length that can be used standalone or in combination with
other linked visualizations. Our approach supports guided, interactive
analyses by suggesting other period lengths to explore, which are ranked based
on two quality measures. We further describe how the phase can be mapped to
visual representations in other views to reveal periodicity there.

Our approach is centered around the phase histogram widget (a), which
indicates the von Mises distribution [36] for given period lengths as a row-
wise heat map (b). The duration of the period can be adjusted by vertically
scrolling the heat map. The current period length is indicated by a frame (c)
and visualized additionally as a bar chart (d). The quality measures for
period lengths are visualized as a vertical bar chart (e). Bigger changes of
the period length or picking a dedicated one is possible with the time slider
(f). Other interesting period lengths that are fractional multiples of the
current one are suggested above (g) the slider.  Human-centered
computingVisualizationVisualization techniques Introduction

Many datasets contain temporal information. For many years, visualization
techniques have been used to understand developments, trends, and other
temporal patterns in such time data series [1]. Such patterns include
periodically occurring accumulations of events, measured values, or
frequencies. In time series data, one or more attributes depend on time:
$a(t)$. Accordingly, periodic behavior of period length $\tau$ can be seen as
a similarity of the characteristics of $a(t)$ in the intervals
$\left[t_{0}+k\tau,~{}t_{0}+k\tau+\tau\right],~{}k\in\mathbb{Z}$. A special
case are event data, where the temporal distribution of data points is of
interest, rather than the distribution of values over time. Periodic behavior
can be interesting for analyses; for instance, to predict future values of a
data attribute, or to find hidden dependencies in the data. Several well-known
automated procedures for detecting periodic effects exist, such as Fourier
transforms [6], seasonal-trend decomposition based on Loess (STL) [10], and
dynamic mode decomposition (DMD) [27, 20]. However, some of these methods
require the period length of interest as an input parameter, or can be
expensive regarding the compute time.

In certain situations, the _interactive_ visual identification and exploration
of periodic effects has benefits over automatic procedures and static visual
approaches; for example, if the interplay of periodicity in events with their
geographical location or other data attributes is of interest. Different
periodicities can be explored without making presumptions on interesting
periods. In addition, important facets of the time-dependent data become
visible immediately. Various approaches [37, 35, 5, 19, 30] map phase and time
to positions in a grid- or spiral-based layout, revealing periodicity in an
intuitive manner. We extend previous work with a compact, aggregated, and
interactive representation (LABEL:fig:teaser:widget) of time series data that
reveals periodic behavior even for datasets with larger temporal extent. We
offer guidance [9] in the form of suggestions (LABEL:fig:teaser:suggestions)
and quality measures [28, 36]. Since exploring the interplay with other data
attributes can be interesting, our approach can be integrated within a larger,
multiple-view visualization.

The contributions of the approach presented in this work comprise the
introduction of a new composite widget that  (i) helps users detect and
explore periodic occurrences of aggregated quantitative information in long
time series data,  (ii) offers guidance to point users to potentially
interesting periods exhibiting seasonal effects, and  (iii) lets users pick
suitable mappings including glyphs and color scales to understand periodicity
in other views.

##  1 Related Work

We briefly discuss related automatic methods for analyzing periodic behavior,
but focus on interactive visual approaches. Various possibilities for
representing cyclical temporal data exist [1, 5]. Hence, we limit our
discussion here to approaches closely related to ours.

#### Automatic Periodicity Analysis.

Fourier analysis [6] converts an input signal to the frequency domain. This
requires a high sampling rate on the input signal, which can be expensive for
event data, where it must be generated as a long, sparse histogram of the
events. Fourier analysis also returns many false positives for non-sinusoidal
periodicity. Cleveland et al. [10] introduced STL, which splits a time-
dependent data signal into a seasonally recurring component, a linear trend,
and a remainder. Constrained DMD (cDMD) [20] produces similar results using a
different method. Cycle plots [11, 7] produce visualizations that also reveal
such seasonal characteristics. While these methods produce good results even
for noisy data, they are usually expensive to compute, and require prior
knowledge on the frequency of the periodic behavior. STL also requires the
period length to be an integer multiple of the sampling rate. Unconstrained
DMD [27] does not require the period length as input, but has similar
computation costs as cDMD [20], and does not always output the expected period
length. Fourier analysis [6], DMD [27], and cDMD [20] find sinusoidal signal
components; but struggle with other signal characteristics that could occur
with event data. In contrast to these methods, we focus on interactive and
exploratory analysis of periodic behavior in event data, often in the context
of other data attributes.

#### Visualization of Periodicity.

Various works have explored non-linear layouts of data to reveal periodicity.
Carlis et al. [8] and Weber et al. [37] proposed Archimedean spirals, where
one turn along the spiral represents one period. Periodic behavior of that
period length would then appear as lines or cones going radially outward from
the center (LABEL:fig:patterns:spiral). The same effect can be reached with a
line-wise, rectangular representation (LABEL:fig:patterns:rect,
LABEL:fig:patterns:rect-multiple and LABEL:fig:patterns:rect-fraction), where
each line represents one period [19]. In both representations, spiral and
rectangular, the period length can be adjusted interactively, by tightening
the spiral or changing the aspect ratio of the rectangle, to find
periodicities. Our approach focuses particularly on scalability, as well as
guidance of users to interesting period lengths. We show these representations
as detail views on demand. A special case of the rectangular representations
are calendar-based layouts [35, 21, 29], which reveal periodic behavior in
human-made time concepts on multiple scales, but are limited to fixed time
concepts. Several visualizations utilizing concentric circles [2, 3, 16, 23]
or stacks [22] also fit this category. Frey et al. [14] present a matrix
representation which reveals self-similar patterns in time series data, but
requires considerable screen estate.

_Pinus_ view [30] shows an aggregated visual summary of a time series in a
triangular fashion, where each position in the triangle corresponds with an
interval and an offset in the temporal domain of the dataset. The residuality
model by Van de Weghe et al. [34] describes a similar concept. These concepts,
like ours, show a visual summary of the temporal domain. However, they are
orthogonal to ours in that they consider aggregated intervals, rather than
looking at repeating patterns. Closer to our approach is the work by Suschnigg
et al. [31], who compare characteristics between periods via glyphs, or a
matrix of anomaly scores. The approach by Ishii and Misue [17] is also closely
related to our quality metrics, as they essentially visualize the von Mises
[36] vector strength and direction for different data attributes. Recurrence
quantification analysis [24, 25] also share similarities to our work. In
contrast to these works, we focus on interactive discovery of interesting
period lengths, and offer a more compact visual representation.

##  2 Approach

Our approach is centered around an aggregated view on all periods at once for
a given period length, which we call the _phase histogram._ Periodic behavior
for that period length is visible by non-uniform distribution of values (Figs.
LABEL:fig:teaser:heatmap, LABEL:fig:teaser:histogram, LABEL:fig:patterns:ours-
sharp, and LABEL:fig:patterns:ours-fuzzy; Table 1). We calculate quality
measures for different period lengths to guide towards interesting ones,
focusing on fractional multiples of the current period length to balance out
the challenges discussed in Section 2.2.

###  2.1 Phase Histogram

Our approach considers a time series of events $S$, $\left|S\right|=n$. The
time series
$S=\left\\{t_{i}~{}\middle|~{}i\in\left\\{1,\dots,n\right\\}\subset\mathbb{N}\right\\}$
consists of points in time $t_{i}$. It has an extent from $t_{0}$ to $t_{1}$:
$t_{i}\in\left[t_{0},t_{1}\right]$, where $t_{0},t_{1}\in\mathbb{R}$. For a
user-selectable _period length_ $\tau$, each point in time $t$ has a _phase_
$\varphi(t,\tau)\in\left[0,2\pi\right)$. The phase is the offset from the
start of the period relative to the period length, where the first period
starts at $t_{0}$:

$\varphi(t,\tau)=2\pi\cdot\frac{\left(t-t_{0}\right)~{}\mathrm{mod}~{}\tau}{\tau}$

For a given $\tau$, we then calculate a histogram with $N$ bins over the
phases of all time series events (LABEL:fig:patterns:schematic). The choice of
$N$ depends mainly on available space. It affects the granularity of the
resulting phase histogram, as well as the scaling of the entropy estimation
values. In the subsequent examples, we use $N=25$. The shape of the histogram
then reveals potential periodic behavior if one or some of the bins contain
considerably more items than the others. Figures 2b, c, e to k, and Table 1
(rows 2–5) show some examples of histograms which represent periodic behavior.
Multiples and fractions of the actual periodic behavior’s period length will
also produce interesting-looking patterns. We consider such multiples and
fractions to evaluate whether they show clearer periodic behavior than the
current period length $\tau$ (see Section 2.2). While our approach focuses on
event data, it can also be employed to show periodic behavior of time-
dependent data attributes (Section 2.3). Our representation can be understood
as a vertical aggregation of the rectangular binned representations [19, 35,
21, 29], which reveals overarching patterns in all rows (see
LABEL:fig:patterns:schematic).

Figure 1:  A schematic explanation (a) of our approach: Event data gets binned
over the temporal domain. For the Cartesian representations [19, 35, 21, 29]
the bins then get placed so that each row represents one period. Our approach
shows an aggregated view on this, and varies the period durations slightly
from row to row. The current period length is framed in red in the figure.
Example patterns for periodic behavior in our approach for a sharp (b) and a
less sharp (c) periodic pattern, and for uniform noise (d). As a comparison,
periodic behavior (e) in the Cartesian binned representation is shown for the
actual signal period length, and for multiples (f), integer fractions (g), and
non-integer fractions (h, here $\nicefrac{{6}}{{5}}$) thereof. Here, the
periodicity manifests as vertical lines. Patterns appear as straight lines
emanating from the center in the Archimedean binned representation (i) [8,
37], but are also visible if the period length is nearly right (k).

###  2.2 Pre-calculation and Guidance

Table 1:  Example phase histograms alongside the two quality measures used:
Shannon entropy [28] and von Mises vector strength [36]. For the vector
strength, the distribution of the data points on the unit circle (blue) and
their barycenter (red) are also visualized.

Interesting period lengths are those that match periodically reoccurring
events in the data. In these cases, the phases of a larger-than-average number
of events are similar, and the phase histogram is not uniform. Our approach
calculates two quality measures that help decide which period lengths show
promising patterns. The first measure is the Shannon entropy [28] of a phase
histogram. Shannon entropy is high for more uniformly distributed histograms,
and low for histograms with clear peak outliers (Table 1). The second measure
is von Mises vector strength [36], which is a measure between 0 and 1, and is
highest when the phases of all data points are the same. Von Mises vector
strength is approximated by the distance from the center of a unit circle to
the barycenter of all data points after these have been projected onto the
unit circle at the angle of their respective phase. Both measures are suited
for the task, since they generally indicate when the phase histograms are not
uniformly distributed. They are also fast to compute, and so can be used as
part of a highly interactive exploration approach. They have individual
downsides, which is why we calculate both for a more nuanced view on the data:
For Shannon entropy, the values quickly rise with only little noise. Von Mises
vector strength cannot recognize whole multiples of the period length of
periodic behavior, as the barycenter then moves to the center of the unit
circle (Table 1, last two rows).

For a dataset, our approach initially pre-calculates the histograms and
quality measures for a set of period lengths between a lower bound (e.g., 1
min) and an upper bound (the temporal extent of the dataset). Within this
range, multiples of time units are sampled ($1\dots 59\,\mathrm{min}$, $1\dots
23\,\mathrm{h}$, etc.). In addition, exponentially increasing period lengths
are included. Data for additional period lengths can then be requested during
analysis, and are calculated ad-hoc.

We also offer suggestions based on the current period length as a guidance. A
key challenge here is to determine the actual frequency or period length of
periodic behavior: A repeating pattern of period length $\tau$ looks repeating
also for multiples of $\tau$ (Tables 1 and LABEL:fig:patterns:rect-multiple),
but also for fractions thereof, if only an aggregation of all periods is
considered (LABEL:fig:patterns:rect-fraction). Furthermore, fractions of
multiples (such as $\nicefrac{{4}}{{3^{\mathrm{rds}}}}$ of the actual period
length) can look interesting as well in an aggregated view on the data
(LABEL:fig:patterns:rect-real-fraction). To facilitate discovery of better
period lengths that are located at such factors, our approach samples a small
set of fractions $\nicefrac{{k}}{{n}}$, where $n$ is a small natural number
$\geq 2$ and $k\in\left\\{1,\dots,2n-1\right\\}$. Small natural-numbered
multiples of the current period length are sampled as well, and a selection of
the best-ranking suggestions are returned based on available space. This
sampling happens on the fly for the current period length $\tau$.

###  2.3 Visual Representation

We propose a compact widget (LABEL:fig:teaser:widget) that can be included as
part of a larger, coordinated-views visualization. The widget shows the
aggregated phase distribution for the currently selected period
(LABEL:fig:teaser:histogram) length $\tau$, as well as for the immediate
neighborhood period lengths. This context (LABEL:fig:teaser:heatmap) is shown
as a pixel-based [18] visualization or heat map, where each row shows a period
length. Each row visualizes the histogram as a line of colored rectangles,
where the bin’s respective value is mapped to color. The period lengths are
shown in ascending order from the pre-calculated data (Sections 2.2 and
LABEL:fig:patterns:schematic), with a context of $n$ rows above and below the
current row, which is located at the center of the heat map
(LABEL:fig:teaser:current-row). By scrolling on the heat map with the mouse
wheel, or clicking on rows, users can browse the rows, changing the current
period length accordingly. By default, the color of the heat map cells
represents the number of data items in that histogram bin. It is also possible
to map this to other measures calculated on these subsets of data, such as the
mean value or variance of another data attribute (shown in
LABEL:fig:teaser:suggestions for the number of sun spots per day).

To the right of the heat map, an interest measure is indicated for each row by
a vertical bar chart (LABEL:fig:teaser:quality-histogram). For the vector
strength (Section 2.2), which is shown by default, the measure is mapped
directly. For the entropy, the highest possible entropy (least interesting) is
mapped to an empty bar, and the lowest possible entropy (most interesting) is
mapped to a full bar. At the top of the widget, the current period length’s
histogram is visualized again, as a bar chart (LABEL:fig:teaser:histogram). In
applications where the phase is mapped to visual attributes in other views
(see Section 2.4), this visual mapping is displayed in a legend
(LABEL:fig:mapping:legend-color and LABEL:fig:mapping:legend-glyph) below the
bar chart.

Our approach also contains a time slider (LABEL:fig:teaser:slider), which
represents the domain of possible period lengths on a logarithmic scale and
can be used for larger adjustments to the current period length $\tau$. Period
lengths with the highest-ranked quality measure values are indicated as longer
tick marks to guide users. During exploration, multiples and fractions for the
current period length are sampled. The sampled period lengths are ranked by
the selected quality measure, and the most promising ones are suggested as
small thumbnails (LABEL:fig:teaser:suggestions) above the time slider.
Clicking on one such thumbnail changes the current period length to that of
the thumbnail.

###  2.4 Visually Mapping the Phase

The presence of periodicity in the context of other data attributes can be of
interest as well. For example, users might be interested in understanding
whether certain events displayed on a map occur at recurring times (see
Section 3). We have tested this mapping with position-based visualizations in
our approach, but it could be extended to other visualization types. Our
prototype implementation contains a scatter plot (LABEL:fig:mapping:color),
which represents the spatial attribute of the data. To reveal spatio-temporal
periodicity patterns, we propose two different types of mapping that can be
applied to the markers of the scatter plot: color and shape (Fig. 2).

The first variant maps the phase of the data item to a position on a
continuous color scale. Cyclical color scales, such as a rainbow color scale,
can be a reasonable choice in this case. However, some datasets might favor
color scales with a clear cut between the end and the beginning. The second
variant maps to a visual mark with a parameterized shape. For our prototype,
we offer a mapping to a moon phase-like mapping, or the angular rotation
within $\left[0,\pi\right)$ of a rectangle, both cyclical. As acyclic
mappings, we offer star shapes that morph into circles. Previous works [12,
26, 4] on glyph design offer additional options here. Future work might
explore suitable shapes for representing cyclical data attributes.

Within our central widget, we show a legend (LABEL:fig:mapping:legend-color
and LABEL:fig:mapping:legend-glyph) of this mapping. Especially for the
acyclic mappings, it can be beneficial to change the offset of the mapping.
Hence, the offset can be adjusted from within the legend.

##  3 Case Study

We tested our approach with various synthetic and real-world datasets. Here,
we demonstrate its applicability on tidal sea level data obtained from NOAA
[32], with two specific datasets. Such data contains various periodicities
pertaining to the rotation of the earth in relation to the moon, the orbit of
the moon around the earth and the earth around the sun, but also to annual or
perennial patterns stemming from meteorological phenomena. The interplay of
these factors is not always obvious, demonstrating the need for exploratory
analysis of the resulting patterns in the data.

Figure 2:  Example mappings of phase to color (a, c) or shape of visual marks
(b) in a scatter plot. Periodic behavior related to the spatial aspect of the
data is revealed by uniform areas. The mapping is shown in our widget
(LABEL:fig:teaser:widget) as a legend (d, e) that can be interactively
adjusted to change the mapping of glyph or color to phase.

The first dataset was generated from hourly measurements of the mean sea level
(MSL) at one station located in Honolulu, Hawai’i [33]. The station
measurements cover 118 years, from 1905 to the present day. We reduced this
data to the measurements that exceeded a threshold value ($+0.5\,\mathrm{m}$),
leaving us with 5060 events. Our approach guided us to find the spring and
neap tide periodicity at half a sidereal month ($13.66\,\mathrm{d}$,
LABEL:fig:teaser:widget). We also could see and verify the $18.61$-year nodal
cycle discussed by Haigh et al. [15]. By mapping time of year to the $x$ and
year to the $y$ coordinate in a scatter plot (LABEL:fig:mapping:calendar-
honolulu, days increase to the right, years to the top), we can also see an
increase of events over the years, and a seasonal yearly component (denser and
less dense vertical areas).

The second dataset contains the dates with highest extreme water levels for
all observation stations, as provided by NOAA [32]. Over a record interval of
122.5 years, this data only contains 571 events, distributed over 74 unique
stations. Hence, the observed patterns are not as expressive; additional
preprocessing of raw data with more domain knowledge might yield better
results here. Still, we can see some periodicity at the $18.61$-year period
length [15]. Furthermore, seasonal differences in events between the US east
and west coast are clearly visible in the mapping of phase to color in a
geographical scatter plot (LABEL:fig:mapping:color) for a one-year period
length.

##  4 Discussion

We implemented a web-based prototype for our approach. Its backend calculates
the required data, as well as additional detail data on demand. We discuss
initial results, and benefits and disadvantages compared to other visual
representations and automated methods.

#### Suitable Datasets and Use Cases.

Our approach can be applied to use cases where the periodic re-occurrence of
events is of interest. Such events can also be generated from general time
series data by determining points in time where time-dependent data values
match some criteria. We demonstrated the viability of this strategy for
finding periodic behavior in real-world data in our case study (Section 3). We
offer a few sample datasets in our prototype, but also allow for arbitrary
datasets to be loaded by the user.

#### Implementation Details and Scalability.

Our web-based prototype visualizes pre-calculated data, and additional data is
calculated on demand asynchronously in a backend. With this strategy, even for
larger datasets with tens of thousands of events, interaction with the user
interface is very responsive and smooth, with no or very few frame drops.
Suggestions are calculated ad-hoc by the backend and typically appear within
half a second after the last interaction. The thumbnails at the top
(LABEL:fig:teaser:suggestions) work well to guide users to interesting period
lengths. We have found that the two measures supplement each others’ drawbacks
(Section 2.2) quite well. The mapping of phase to color or shape in the
scatter plot has proven useful to detect periodic behavior that is local in
other data attributes (Fig. 2). For performance, redrawing the scatter plot is
the bottleneck, but we have found datasets with thousands or tens of thousands
of data points to still render interactively. The implementation is available
on GitHub and Zenodo [13].

#### Comparison to automated analyses.

STL [10] and related methods [20] allow for the automated analysis of periodic
behavior, and handle noisy data better. These methods require prior knowledge
of the period length for which the periodic behavior should be analyzed. We
compared mean computation time over ten runs of our method, STL, and DMD for
two datasets: For a smaller synthetic dataset, STL took $227.2\,\mathrm{ms}$
to calculate for _one_ period length of interest, and $3.45\,\mathrm{s}$ for
the Hawai’i tide dataset [33]. DMD took $2.53\,\mathrm{s}$ and
$1.63\,\mathrm{s}$, respectively, to produce decompositions with six
components and a sensible delay parameter, but did not find the expected
periodic behavior. In contrast, our method took $470\,\mathrm{ms}$ and
$589.5\,\mathrm{ms}$, respectively, to calculate the phase histograms and
quality measures for over $1800$ sampled period lengths. The fairly high
computation costs of methods like STL or cDMD prohibit extensive pre-
calculations, or on-the-fly computation during smooth interaction. We also see
that DMD struggles to produce good results for non-sinusoidal periodic
behavior. However, these methods could still be utilized as a second step to
verify results found interactively with our approach on demand. Fourier
analysis [6] is another option, but not always suitable to real-world
periodicity patterns. We argue that automated analysis is most useful for
targeted use cases, but that exploratory analysis, especially in the context
of other data attributes (Section 2.4), still requires more interactive
methods.

#### Comparison to other visual representations.

Cycle plots [11, 7] are a powerful way to visualize periodicity, but again
presume knowledge of the periodic behavior. The calendar-based representations
[35, 21] presume specific period lengths as well, but are often suitable for
data based on from human-specified time concepts, such as months. We think
other rectangle-[19] and spiral-based [8, 37] visualizations show periodic
behavior in an intuitive manner, especially when period length can be
adjusted. However; these representations do not scale well for larger time
spans; and interactive use can often lead to large shifts, bad aspect ratios,
and flickering effects. Our approach condenses these two-dimensional
representations down to a one-dimensional aggregation. This opens up space for
visualizing different data attributes to understand relations between them,
and allows showing close-by period lengths as a context to the current one.
This supports users in finding the best local period for periodic events. The
aggregation hides periodicity that only appears for a part of the temporal
extent, and multiples or fractions of the actual periodicity are not clearly
identifiable as such (Section 2.2). Hence, we offer guidance, as well as
detailed views as a tooltip to closer examine interesting period lengths (Fig.
1).

#### Open challenges.

We demonstrated the general applicability of our approach with a prototype and
different datasets, and plan to explore its integration into larger
visualization systems. Understanding the best choices with respect to phase
mapping (Section 2.4) in different scenarios, and determining the most
appropriate visual mapping in terms of color scales and glyph designs to
represent periodicity, is another interesting research direction we plan to
pursue. The integration of other quality measures such as variance of the
phase, domain-specific measures, or the combination of multiple measures,
could be an interesting extension. Incorporating more complex automatic
analysis methods, such as STL, could be another good addition to verify
promising findings on demand.

##  5 Conclusion

We have presented a novel, aggregated visual representation of temporal data
to explore periodicity. The time domain is mapped to the phase or remainder
for a given and adjustable period length and visualized in a binned manner.
Patterns that repeat with the same period length appear clearly in our
representation. Our approach offers a more compact, interactive alternative to
existing representations. It can be employed as a part of visualization
systems with many views to let users configure the mapping of phases for
depicting periodicity in other visual contexts.

## Acknowledgments

This work has been partially funded by the German Research Foundation (DFG)
project # 314647693. We thank Tim Krake in particular; as well as Frank Heyen,
Daniel Klötzl, and Kuno Kurzhals; for their valuable feedback and assistance.

## References

  * [1] W. Aigner, S. Miksch, H. Schumann, and C. Tominski.  Visualization of Time-oriented Data.  Springer, London, UK, 2011.
  * [2] E. N. Argyriou and A. Symvonis.  Detecting periodicity in serial data through visualization.  In Advances in Visual Computing, pp. 295–304. Springer Berlin Heidelberg, Jul 2012. doi: 10 . 1007/978-3-642-33191-6_29
  * [3] K. Bale, P. Chapman, N. Barraclough, J. Purdy, N. Aydin, and P. Dark.  Kaleidomaps: A new technique for the visualization of multivariate time-series data.  Information Visualization, 6(2):155–167, Jan 2007. doi: 10 . 1057/palgrave . ivs . 9500154
  * [4] R. Borgo, J. Kehrer, D. H. S. Chung, E. Maguire, R. S. Laramee, H. Hauser, M. Ward, and M. Chen.  Glyph-based visualization: Foundations, design guidelines, techniques and applications.  In Eurographics STARs, pp. 39–63. The Eurographics Association, Aug 2013. doi: 10 . 2312/conf/eg2013/stars/039-063
  * [5] M. Brehmer, B. Lee, B. Bach, N. H. Riche, and T. Munzner.  Timelines revisited: A design space and considerations for expressive storytelling.  IEEE Transactions on Visualization and Computer Graphics, 23(9):2151–2164, Sep 2017. doi: 10 . 1109/tvcg . 2016 . 2614803
  * [6] E. O. Brigham and R. E. Morrow.  The fast Fourier transform.  IEEE Spectrum, 4(12):63–70, Dec 1967. doi: 10 . 1109/mspec . 1967 . 5217220
  * [7] M. Bögl, P. Filzmoser, T. Gschwandtner, T. Lammarsch, R. A. Leite, S. Miksch, and A. Rind.  Cycle plot revisited: Multivariate outlier detection using a distance-based abstraction.  Computer Graphics Forum, 36(3):227–238, Jun 2017. doi: 10 . 1111/cgf . 13182
  * [8] J. V. Carlis and J. A. Konstan.  Interactive visualization of serial periodic data.  In Proc. UIST, pp. 29–38. ACM, Nov 1998.
  * [9] D. Ceneda, T. Gschwandtner, and S. Miksch.  A review of guidance approaches in visual data analysis: A multifocal perspective.  Computer Graphics Forum, 38(3):861–879, Jun 2019. doi: 10 . 1111/cgf . 13730
  * [10] R. B. Cleveland, W. S. Cleveland, J. E. McRae, and I. Terpenning.  STL: A seasonal-trend decomposition procedure based on Loess.  Journal of Official Statistics, 7(1):3–73, 1990.
  * [11] W. S. Cleveland.  Visualizing data.  Hobart Press, Sep 1993.
  * [12] D. S. Ebert, R. M. Rohrer, C. D. Shaw, P. Panda, J. M. Kukla, and D. A. Roberts.  Procedural shape generation for multi-dimensional data visualization.  Computers & Graphics, 24(3):375–384, Jun 2000. doi: 10 . 1016/s0097-8493(00)00033-9
  * [13] M. Franke.  UniStuttgart-VISUS/periodic-time-vis.  Zenodo, Jul 2023.  v1.0.0. doi: 10 . 5281/zenodo . 8164734
  * [14] S. Frey, F. Sadlo, and T. Ertl.  Visualization of temporal similarity in field data.  IEEE Transactions on Visualization and Computer Graphics, 18(12):2023–2032, Aug 2012. doi: 10 . 1109/tvcg . 2012 . 284
  * [15] I. D. Haigh, M. Eliot, and C. Pattiaratchi.  Global influences of the 18.61 year nodal cycle and 8.85 year cycle of lunar perigee on high tidal levels.  Journal of Geophysical Research, 116(C6):C06025, Jun 2011. doi: 10 . 1029/2010jc006645
  * [16] M. Hao, M. Marwah, S. Mittelstädt, H. Janetzko, D. Keim, U. Dayal, C. Bash, C. Felix, C. Patel, M. Hsu, Y. Chen, and M. Hund.  Visual analytics of cyber physical data streams using spatio-temporal radial pixel visualization.  In P. C. Wong, D. L. Kao, M. C. Hao, C. Chen, and C. G. Healey, eds., SPIE Proceedings. SPIE, Feb 2013. doi: 10 . 1117/12 . 2002948
  * [17] T. Ishii and K. Misue.  2.5D extension of ChronoView for exploring periodic features of temporal data.  In Proc. IV. IEEE, Jul 2018. doi: 10 . 1109/iv . 2018 . 00014
  * [18] D. Keim.  Designing pixel-oriented visualization techniques: Theory and applications.  IEEE Transactions on Visualization and Computer Graphics, 6(1):59–78, 2000. doi: 10 . 1109/2945 . 841121
  * [19] R. Kosara.  Spirals for periodic data.  https://eagereyes.org/techniques/spirals, Aug 2011.  Accessed: 2023-02-22.
  * [20] T. Krake, D. Klötzl, B. Eberhardt, and D. Weiskopf.  Constrained dynamic mode decomposition.  IEEE Transactions on Visualization and Computer Graphics, pp. 182–192, Sep 2022. doi: 10 . 1109/tvcg . 2022 . 3209437
  * [21] T. Lammarsch, W. Aigner, A. Bertone, J. Gärtner, E. Mayr, S. Miksch, and M. Smuc.  Hierarchical temporal patterns and interactive aggregated views for pixel-based visualizations.  In Proc. IV. IEEE, Jul 2009. doi: 10 . 1109/iv . 2009 . 52
  * [22] T.-Y. Lee, A. Chaudhuri, F. Porikli, and H.-W. Shen.  CycleStack: Inferring periodic behavior via temporal sequence visualization in ultrasound video.  In Proc. PacificVis. IEEE, Mar 2010. doi: 10 . 1109/pacificvis . 2010 . 5429602
  * [23] G. C. Mariano, V. G. Staggemeier, L. P. Cerdeira Morellato, and R. da Silva Torres.  Multivariate cyclical data visualization using radial visual rhythms: A case study in phenology analysis.  Ecological Informatics, 46:19–35, Jul 2018. doi: 10 . 1016/j . ecoinf . 2018 . 05 . 003
  * [24] N. Marwan, M. C. Romano, M. Thiel, and J. Kurths.  Recurrence plots for the analysis of complex systems.  Physics Reports, 438(5–6):237–329, Jan 2007. doi: 10 . 1016/j . physrep . 2006 . 11 . 001
  * [25] T. Rawald, M. Sips, and N. Marwan.  PyRQA—Conducting recurrence quantification analysis on very long time series efficiently.  Computers & Geosciences, 104:101–108, Jul 2017. doi: 10 . 1016/j . cageo . 2016 . 11 . 016
  * [26] T. Ropinski, S. Oeltze, and B. Preim.  Survey of glyph-based visualization techniques for spatial multivariate medical data.  Computers & Graphics, 35(2):392–401, Apr 2011. doi: 10 . 1016/j . cag . 2011 . 01 . 011
  * [27] P. J. Schmid.  Dynamic mode decomposition of numerical and experimental data.  Journal of Fluid Mechanics, 656:5–28, Jul 2010. doi: 10 . 1017/s0022112010001217
  * [28] C. E. Shannon.  A mathematical theory of communication.  Bell System Technical Journal, 27(4):623–656, Oct 1948. doi: 10 . 1002/j . 1538-7305 . 1948 . tb00917 . x
  * [29] P. Silva, C. Macas, E. Polisciuc, and P. Machado.  Visualisation tool to support fraud detection.  In Proc. IV. IEEE, Jul 2021. doi: 10 . 1109/iv53921 . 2021 . 00022
  * [30] M. Sips, P. Kothur, A. Unger, H.-C. Hege, and D. Dransch.  A visual analytics approach to multiscale exploration of environmental time series.  IEEE Transactions on Visualization and Computer Graphics, 18(12):2899–2907, Dec 2012. doi: 10 . 1109/tvcg . 2012 . 191
  * [31] J. Suschnigg, B. Mutlu, G. Koutroulis, V. Sabol, S. Thalmann, and T. Schreck.  Visual exploration of anomalies in cyclic time series data with matrix and glyph representations.  Big Data Research, 26:100251, Nov 2021. doi: 10 . 1016/j . bdr . 2021 . 100251
  * [32] Tides & Currents, National Oceanic and Atmospheric Administration, United States of America.  Extreme water levels.  https://tidesandcurrents.noaa.gov/est/.  Accessed: 2023-04-25.
  * [33] Tides & Currents, National Oceanic and Atmospheric Administration, United States of America.  Station home page of Honolulu, HI.  https://tidesandcurrents.noaa.gov/stationhome.html?id=1612340.  Accessed: 2023-04-26.
  * [34] N. Van de Weghe, R. Docter, P. De Maeyer, B. Bechtold, and K. Ryckbosch.  The triangular model as an instrument for visualising and analysing residuality.  Journal of Archaeological Science, 34(4):649–655, Apr 2007. doi: 10 . 1016/j . jas . 2006 . 07 . 007
  * [35] J. J. van Wijk and E. R. van Selow.  Cluster and calendar based visualization of time series data.  In Proc. InfoVis. IEEE, Oct 1999. doi: 10 . 1109/infvis . 1999 . 801851
  * [36] R. von Mises.  Über die “Ganzzahligkeit” der Atomgewichte und verwandte Fragen.  Physikalische Zeitschrift, 19:490–500, 1918.
  * [37] M. Weber, M. Alexa, and W. Müller.  Visualizing time-series on spirals.  In Proc. InfoVis. IEEE, Oct 2001. doi: 10 . 1109/infvis . 2001 . 963273