11institutetext:  Centre for Nonlinear Science & Engineering, School of
Electrical & Electronics Engineering,
SASTRA Deemed University, Thanjavur -613 401, Tamil Nadu, India.
Centre for Nonlinear Dynamics, School of Physics, Bharathidasan University,
Tiruchirappalli - 620 024,
Tamil Nadu, India.

Synchronization; coupled oscillators Self-organized systems Nonlinear dynamics
and chaos

# Conjugate coupling induced symmetry breaking and quenched oscillations

K. Ponrasu 11 K. Sathiyadevi 11 V. K. Chandrasekar 11 M. Lakshmanan 221122

###### Abstract

Spontaneous symmetry breaking (SSB) is essential and plays a vital role many
natural phenomena, including formation of Turing pattern in organisms and
complex patterns in brain dynamics. In this work, we investigate whether a set
of coupled Stuart-Landau oscillators can exhibit spontaneous symmetry breaking
when the oscillators are interacting through dissimilar variables or conjugate
coupling. We find the emergence of SSB state with coexisting distinct
dynamical states in the parametric space and show how the system transits from
symmetry breaking state to out-of-phase synchronized (OPS) state while
admitting multistabilities among the dynamical states. Further, we also
investigate the effect of feedback factor on SSB as well as oscillation
quenching states and we point out that the decreasing feedback factor
completely suppresses SSB and oscillation death states. Interestingly, we also
find the feedback factor completely diminishes only symmetry breaking
oscillation and oscillation death (OD) states but it does not affect the
nontrivial amplitude death (NAD) state. Finally, we have deduced the
analytical stability conditions for in-phase and out-of-phase oscillations, as
well as amplitude and oscillation death states.

###### pacs:

05.45.Xt

###### pacs:

05.65.+b

###### pacs:

05.45.−a

##  1 Introduction

Coupled nonlinear oscillatory systems exhibit distinct collective dynamical
behaviors including various synchronization patterns [1]-[4], and oscillation
quenching states [5]-[8], which have been observed in physical, chemical and
biological systems [9]-[11]. Among them, one of the intriguing phenomena
observed recently is the symmetry breaking states which play a crucial role in
many natural systems. Importantly, in biological systems the spontaneous
symmetry breaking (SSB) is essential for cell movement, developmental
patterning of vertebrates and also for Turing pattern formation in
organism[12]-[13]. In the neural network, SSB leads to the formation of
diverse complex patterns [14]. Due to the above features, the mechanism
underlying such SSB and the situation in which such SSB occurs are of great
interest. Recently, the trade-off between the attractive and repulsive
coupling induced SSB has been reported [15].

Oscillation quenching is another emergent intriguing phenomenon found in
various chemical and biological systems [16]-[19]. The quenching of
oscillation is further subdivided into two categories, i.e. (i) amplitude
death (AD) and (ii) oscillation death (OD). The AD exists when the coupled
oscillators interact in such a way as to suppress each other’s oscillations
and collectively attain a homogeneous steady state [20]. In the earlier
reports, it has been reported that three main factors can cause AD, namely (i)
coupling through dissimilar variables [21], (ii) time delay [22]-[25] and
(iii) frequency mismatch of the oscillators [26]. Later, such dynamical
behavior was also observed in dynamic coupling [27], direct-indirect coupling
[28] and even in mean field diffusive coupling [7]. In many practical
situations the nontrivial AD is a desirable control mechanism to suppress
harmful oscillations in laser systems and neuronal activity, where a constant
output is needed and fluctuations should be suppressed [29]-[31]. Further, the
coupled oscillators can exhibit another type of oscillation quenching
phenomenon, namely oscillation death (OD) where the strong interaction among
the oscillators leads the system to attain an inhomogeneous steady state
(IHSS) [17],[32]. The OD can occur under distinct coupling schemes such as
time-delayed coupling [33]-[34] conjugate coupling [35], mean field
diffusive[7] and repulsive interactions [36]. The OD state exhibits
significant implications in many biological systems such as synthetic genetic
networks[37]-[38], cellular differentiation [39] and cardiovascular
phenomena[40]. Even though the OD state has many practical applications, in
some situations the oscillation death is detrimental and should be avoided.
Later, many investigations have been carried out to control oscillation death
or revival of oscillation from death and several authors have reported the
ubiquitous phenomenon of feedback as a mechanism to control death state
[41]-[43].

Moreover, many investigations have been accomplished to study the dynamical
behavior in coupled oscillators where the coupling among the oscillations is
established through similar variables. On the other hand, in many real
situations coupling via dissimilar variable is also desirable and such a
coupling is known as conjugate coupling. In a variety of experimental
situations the coupling via a conjugate variable is essential where the
subsystems are coupled by feeding the output of one into the other. For
example, in laser experiments, [44] the light emitted from each laser diode is
monitored by a photo-diode detector (whose AC signal is amplified or reduced)
by feed back to the other. Further, in the ecological context, the dynamics of
two predator-prey systems that are coupled via cross predation is also
investigated, in which each predator consumes prey in the other system which
leads to stabilizing the food web to a new equilibrium [45]. In addition,
theoretical predictions reveal that time delay is not necessary, and coupling
with dissimilar variables can lead to AD even when the coupling is
bidirectional and diffusive [21]. Further, the emergence of a large number of
states such as homogeneous and inhomogeneous oscillation death states and
homogeneous oscillatory states have been found in a ring of identical Stuart-
Landau oscillators where each oscillator interacts with its nearest neighbors
via a diffusive conjugate coupling [35]. Recently, shrinking of AD regions
were reported by adjusting the asymmetry factors in conjugately coupled
systems [46] and the emergence of explosive death has also been found in a
network of coupled van der Pol oscillators while coupling via conjugate
variables [47].

Motivated by the above studies, in this letter, we study the dynamical
behavior in a system of conjugately coupled Stuart-Landau oscillators. We
investigate, whether the conjugate coupling can be able to induce spontaneous
symmetry breaking state. Interestingly, we find that the conjugate coupling
among the oscillators breaks the system symmetry spontaneously and gives rise
to SSB state along with distinct multistable regions. Through a study of the
associated basin of attraction [48], we show how the transition occurs from
SSB state to OPS state. In particular, the transition to OPS state takes place
via two distinct types of transitions: (i) sudden transition from SSB state to
OPS state and (ii) continuous transition (i.e. increasing parameter increases
the OPS basin while decreasing the SSB basin) to OPS state. In addition, we
analyze the effect of feedback by introducing the feedback factor in the
system and show that decreasing feedback factor suppresses the symmetry
breaking oscillations. Moreover, we find the mean field conjugate coupling can
induce the oscillation death (OD) state which coexists with nontrivial
amplitude death (NAD) state. The impact of feedback factor is also
investigated on both OD and NAD states. The observed oscillation death state
shrinks while decreasing the feedback factor. The feedback factor only
suppresses oscillation death but it does not affect the nontrivial amplitude
death state. This type of dynamical behavior can be useful and have potential
applications in laser systems where we need a constant output and the
fluctuation should be shorted out.

##  2 Model

In order to exemplify our results, we consider a general, paradigmatic model
of Stuart-Landau limit cycle oscillators [49], which are coupled through mean-
field of the dissimilar variables. The governing equations can be given by

$\displaystyle\dot{z_{k}}=f(z_{k})+\epsilon[{\mathrm{Im}}{(\mathbf{Z})}-a\,{\mathrm{Re}}{(z_{k})}],\quad{k=1,2}.$
(1)

where $f(z_{k})=(\lambda+i\omega-|z_{k}|^{2})z_{k}$ and
$z_{k}=re^{i\phi_{k}}=x_{k}+iy_{k}\in C$, (${k=1,2}$). $x_{k}$ and $y_{k}$ are
the state variables of the system. $\mathbf{Z}=\frac{z_{1}+z_{2}}{2}$, is the
mean field of the coupled system. The intrinsic parameter of the system is
denoted by $\lambda$ which corresponds to the distance from the Hopf
bifurcation point [50]. $\omega$ is the natural frequency, $\epsilon$ be the
coupling strength and $a$ $(0\leq a\leq 1)$ is the feedback control parameter.
To solve Eq. (1) we use Runge-Kutta fourth order method with time step value
of 0.01.

Figure 1: Temporal dynamics of the collective system (1) as a function of
coupling strength $\epsilon$. (a) In-phase synchronized oscillation (IPS) for
$\epsilon=0.2$, (b) nontrivial amplitude death (NAD) state for $\epsilon=0.4$,
(c) symmetry breaking oscillation (SSB) for $\epsilon=0.5$, (d) out-phase
synchronized oscillation (OPS) state for $\epsilon=0.6$ and (e) oscillation
death (OD) state for $\epsilon=0.9$. Other parameters: $\omega=0.5$, $a=1.0$
and $\lambda=1.0$.

##  3 Dynamical behavior and symmetry breaking dynamics

At first, the dynamical behavior of the collective system is inspected through
time series for a fixed value of natural frequency $\omega=0.5$ (for fixed
values of $\lambda=1.0$ and $a=1.0$) and by varying the coupling strength
$\epsilon$ in Fig. 1. Emergence of in-phase synchronized (IPS) oscillations is
observed at $\epsilon=0.2$ as shown in Fig. 1(a). In this state, both the
oscillators oscillate with the same amplitude and phase. On increasing the
coupling strength to $\epsilon=0.4$, both the oscillators attain a nontrivial
homogeneous steady state or amplitude death (NAD) state (see Fig. 1(b)). Fig.
1(c) depicts that at a critical value of the coupling strength $\epsilon=0.5$,
the system exhibits a spontaneous symmetry breaking (SSB) state, where both
the oscillators oscillate with different amplitudes. It is due to spontaneous
breaking of the permutational or translational symmetry $(z_{1}\pm z_{2})$ of
the system. Further, increasing of the coupling strength to $\epsilon=0.6$,
one can note that the oscillators adjust their amplitudes so as to give rise
out-of-phase synchronized (OPS) oscillations as shown in Fig. 1(d). Increasing
the coupling strength to a higher value $\epsilon=0.9$ leads to oscillation
death (OD) state where both the oscillators attain inhomogeneous steady states
(IHSS) ( see Fig. 1(e)).

Figure 2:  Bifurcation diagram for a pair of conjugately coupled Stuart-Landau
oscillators for $\omega=0.5$ and the enlarged image of dashed lines (grey) is
shown in the inset. In the figure dotted (black) lines and solid (red) line
represent the unstable and stable oscillation death states. Filled circles and
triangles denote stable in-phase synchronized and out-of-phase synchronized
states, respectively. Spontaneous symmetry breaking states are denoted by
diamond points (pink) and the unfilled circles represent unstable limit cycle
oscillations. $SN$, $PB_{P}$ and $PB_{s}$ are the saddle-node, pitchfork
bifurcation of periodic orbits and subcritical pitchfork bifurcation points,
respectively. Fig. 2 is also plotted along the line $L_{1}$ in Fig. 3. Other
parameters are same as in Fig. 1.

In addition, to understand the dynamical transitions and multistability of the
considered system, a bifurcation analysis has been carried out by fixing the
natural frequency at $\omega=0.5$. The corresponding bifurcation diagrams are
illustrated in Fig 2 (obtained using XPPAUT software). At this frequency
$\omega=0.5$, the systems exhibit IPS state for lower values of coupling
strength $\epsilon$ (see Fig. 2) and it is evident that increasing $\epsilon$
destabilizes the IPS state and gives birth to nontrivial amplitude death (NAD)
state. The observed NAD state is stable for further entire range of coupling
strength. In addition, the unstable out-of-phase synchronized (OPS) state is
coexistant with stable IPS state and gets stabilized via pitchfork bifurcation
of periodic orbit $(PB_{P})$ which is clearly depicted in the inset of Fig. 2.
Before stabilization of OPS state, we find that there arises spontaneous
symmetry breaking (SSB) state due to the breaking up of the permutational or
translational symmetry $(z_{1}\pm z_{2})$ of the system. The OPS state is
further destabilized while increasing $\epsilon$, which gives birth to
oscillation death (OD) state at strong values of $\epsilon$ via subcritical
pitchfork bifurcation ($PB_{s}$). Additionally, the coexistence of NAD state
is observed in the OD region (see Fig. 2). Besides, the bifurcation analysis
indicates that both the oscillators in the IPS state (i.e. both the
oscillators oscillate with same amplitude and phase) transit to nontrivial AD
state (i.e. the oscillators attain the homogeneous steady state with
nontrivial values) on increasing the strength of coupling. On the other hand,
the OPS state (oscillations with same amplitude and anti-phase) emerges from
the SSB (oscillators oscillate with different amplitudes) state and transit to
oscillation death state on increasing the coupling strength. The above
detailed bifurcation analysis was carried out along the line $L_{1}$ in Fig.
3(a). The inset in Fig. 2 illustrates the enlarged image of the dashed lines
(grey), which delineates the dynamical transition in the symmetry breaking
regions. Before the stabilization of the OPS state, the emergence of stable
SSB state is evident from the inset of Fig. 2, which is further destabilized
via saddle-node (SN) bifurcation. In addition, as noted before we find that
the OPS state is stabilized via pitchfork bifurcation of periodic orbits
($PB_{P}$). The shaded region in the inset denotes the coexistence of OPS and
SSB states. Comparing Fig. 2 and the inset, one can note that the shaded
region denotes the tristabilty region where SSB, OPS and NAD states are
stable. The global dynamical behavior of Eq. (1) is further detailed in the
following section.

##  4 Global behavior of the coupled system

Figure 3: (a) Two parameter phase diagram in ($\epsilon,\omega$) space. (b)
Enlarged image of spontaneous symmetry breaking region. IPS and OPS represent
stable in-phase synchronized and out-of-phase synchronized regions. NAD and OD
are the nontrivial amplitude death and oscillation death states. Spontaneous
symmetry breaking state is denoted by SSB. $R_{1}$, $R_{2}$, $R_{3}$ and
$R_{4}$ represent the bistability regions of IPS-OPS, OPS-NAD, SSB-NAD and
NAD-OD, respectively. The tristability region of SSB-NAD-OPS is represented by
$R_{5}$. The boundary lines of each state denote the analytically obtained
stability curves. Other parameters are $\lambda=1.0$, $a=1.0$.

In order to understand the global dynamical behavior of the system (1), a two
parameter phase diagram is constructed in $(\epsilon,\omega)$ space as shown
in Fig. 3. From Fig. 3(a), at lower values of $\omega$, for increasing
coupling we have observed stable IPS states at very low ranges of $\epsilon$,
and increment of $\epsilon$ further gives rise to NAD state for all values of
coupling strength. In this case, the IPS state loses its stability and
ultimately the system transits to stable NAD state. Increase of coupling
strength also exhibits bistability between OD and NAD states. At lower values
of coupling strength, on increasing the natural frequency the system transits
from monostable nature of IPS state to bistable state. In other words,
coexistence of IPS and OPS states can be observed. It is also noticed that at
critical values of $\epsilon$ and $\omega$, the system exhibits SSB. In this
state the transition to OPS state takes place through symmetry breaking
oscillations and the NAD and OPS states also coexist with SSB. At strong
values of $\epsilon$ and $\omega$ the system admits only NAD state. The
regions $R_{1}$, $R_{2}$, $R_{3}$ and $R_{4}$ are the bistability regions
between IPS-OPS, OPS-AD, SSB-NAD and NAD-OD, respectively. Further, to inspect
the dynamical transitions in the symmetry breaking region, we enlarge the
symmetry breaking region in Fig. 3(b) where we find two distinct types of
transitions along the lines $L_{1}$ and $L_{2}$ which are illustrated through
a study of basin of attraction in Fig. 4. Interestingly, in the SSB region we
find the coexistence of three stable states (i.e. tristability), namely SSB,
OPS and NAD, which is denoted as $R_{5}$ in Fig. 3(b). Additionally, we have
investigated how the transition takes place from the SSB to OPS state in the
following section.

##  5 Dynamical transition through basin of attraction in SSB region

Figure 4:  Basin of attraction obtained by fixing $y_{1}(0)=-0.45$ and
$y_{2}(0)=0.5$ while varying $x_{1}(0)$ and $x_{2}(0)$. Figs. (A) are plotted
along the line $L_{1}$ in Fig. 3(b) at $\omega=0.5$, for coupling strengths
(a) $\epsilon=0.50$, (b) $\epsilon=0.507$, (c) $\epsilon=0.509$ and (d)
$\epsilon=0.52$. Figs. (B) are pictured along the line $L_{2}$ at
$\omega=0.43$ for (a) $\epsilon=0.576$ and (b) $\epsilon=0.577$. Other
parameters: $\lambda=1.0$ and $a=1.0$.

In order to understand the multistabilities among the dynamical states and
transition route to OPS in the SSB region, we portray the basin of attraction
in Fig. 4 by fixing $y_{1}(0)=-0.45$ and $y_{2}(0)=0.5$ and by varying
$x_{1}(0)$ and $x_{2}(0)$. Fig. 4(A): (a)-(d) are plotted for the dashed line
along $L_{1}$ in Fig. 3(b) at $\omega=0.5$. At coupling strength
$\epsilon=0.5$, we find the coexistence of SSB and NAD states in the basin of
attraction, that is SSB state for some initial conditions and NAD state for
other initial conditions (see Fig 4A(a)). Here, we notice that the symmetric
initial conditions (i.e. $x_{1}(0)$ and $x_{2}(0)$ both are either positive or
negative) lead to NAD whereas asymmetric initial states give rise to SSB
states (i.e. in the first and third quadrants where either $x_{1}(0)$ or
$x_{2}(0)$ is positive while the other one is negative). On increasing the
coupling strength $\epsilon=0.507$, the emergence of OPS branch inside the SSB
region is observed in the basin of attraction. Here, the coexistence of NAD,
SSB and OPS states are noticed with respect to distinct initial conditions. On
further increment of $\epsilon$ to 0.509, the basin of attraction for OPS
state is increased and the corresponding SSB basin is decreased. Increasing of
coupling strength further, it is seen that the SSB basin is completely
occupied by the OPS basin. In this way the SSB state transits to OPS state
through increasing of OPS basin smoothly as a function of coupling strength.
The transition is also analyzed along the line $L_{2}$ in Fig. 3(b) at
$\omega=0.43$ in the SSB region, and the corresponding basin of attraction is
depicted in Figs. 4 B(a)-B(b). The SSB state in the first and third quadrant
is directly replaced by OPS basin during a small increment of coupling
strength from $\epsilon=0.576$ to $\epsilon=0.577$. From the basin stability
analysis, we find that the coexistence of distinct dynamical states in the
basin of attraction causes multistability in the parametric space and we
confirm that two distinct types of transition to OPS from the SSB exist and
also coexistence of SSB-NAD, SSB-OPS-NAD and OPS-NAD states. From the basin of
attraction, it is also clear that the symmetric initial states lead only to
NAD and IPS states whereas the asymmetric initial states give rise to OPS and
SSB states. Further, the effect of feedback factor on stable SSB, NAD and OD
states are demonstrated in the following section.

##  6 Impact of feedback factor on spontaneous symmetry breaking (SSB) state
and oscillation quenching states

To exemplify the effect of feedback factor on SSB we have plotted the regions
of SSB for three different values of feedback factors, namely $a=1.0$, $a=0.4$
and $a=0.2$ in the ($\epsilon,\omega$) space in Fig. 5(a). It is evident that
there is a decreasing of SSB regions with decreasing feedback factor $a$. From
Fig. 5(a), it is observed that the SSB region is relatively large for $a=1$
compared to the SSB regions with $a=0.4$ and $a=0.2$. Decreasing the feedback
factor decreases the SSB region and finally there will be no symmetry breaking
oscillations for $a=0$ which completely suppresses the SSB states.
Additionally the impact of feedback factor is also analyzed in NAD and OD
regions. The shrinking of oscillation death is evident from Fig. 5(b). The
shaded region in Fig. 5(b) denotes the OD region for $a=1$ and distinct lines
represent different boundaries of NAD and OD for different values of feedback
factor. Decreasing the feedback factor only shrinks the oscillation death
state but not the nontrivial AD state which sustains even in the absence of
feedback factor. Fig. 5(c) is plotted for $\omega$ as a function of feedback
factor $a$. With respect to the feedback factor the boundary of the OD region
attains null values of $\omega$ but the nontrivial AD states take nonzero
values of $\omega$ which confirms the retaining of AD regions. Furthermore to
confirm the observed results, the bifurcation diagram is plotted at $a=0$ in
Fig. 5(d) which clearly depicts the emergence of NAD even in the absence of
the feedback factor. It is observed that the oscillatory region is increased
compared to the presence of feedback factor. Finally the system shows stable
NAD at large coupling strengths and it confirms that the OD is state entirely
suppressed. From the above results we conclude that the the feedback
completely suppresses the SSB and OD regions only but it sustains the NAD
region.

Figure 5: (a) Regions of symmetry breaking oscillations in ($\epsilon,\omega$)
space with respect to three distinct values of feedback factor $a$. (b) Stable
boundaries of NAD and OD region in ($\epsilon,\omega$) space as a function of
feedback factor $a$. (c) maximum stability region of NAD and OD at
$\epsilon=0.5$ with respect to $a$. (d) Bifurcation diagram for $a=0$,
$\omega=0.5$. Here $\lambda$ is fixed as $1$.

The stability of the observed in-phase and out-phase oscillations, NAD and OD
states are detailed in the next section.

##  7  Stability of distinct dynamical states

In order to study the stability of the observed dynamical states, we first
rewrite Eq. (1) by separating the symmetric ($z_{+}$) and anti-symmetric
($z_{-}$) variables, where $z_{+}={(z_{1}+z_{2})}/{2}$ and
$z_{-}={(z_{1}-z_{2})}/{2}$. Eq. (1) in terms of these new variables is given
by

$\displaystyle\dot{z}_{+}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(f(z_{+}+z_{-})+f(z_{+}-z_{-})\right)+i\epsilon\,\mathrm{Im}(z_{+})$
(2a) $\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad-\epsilon
a\,\mathrm{Re}(z_{+}),$ $\displaystyle\dot{z}_{-}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(f(z_{+}+z_{-})-f(z_{+}-z_{-})\right)-\epsilon
a\,\mathrm{Re}(z_{-}).$ (2b)

In the in-phase subspace, $Z_{+}=\\{(z_{+},z_{-})|z_{-}\equiv 0\\}$, and in
the anti-symmetric subspace, $Z_{-}=\\{(z_{+},z_{-})|z_{+}\equiv 0\\}$. Thus
in the symmetric and anti-symmetric subspaces, the dynamical equations can be
reduced, respectively, to

$\displaystyle\dot{z}_{+}=f(z_{+})+i\epsilon\,\mathrm{Im}(z_{+})-\epsilon
a\,\mathrm{Re}(z_{+}),\;\;\dot{z}_{-}=0,$ (3)

and

$\displaystyle\dot{z}_{-}=f(z_{-})-\epsilon{a}\,\mathrm{Re}(z_{-}),\;\;\dot{z}_{+}=0.$
(4)

We find the explicit expressions for the different oscillatory states, steady
states and their stabilities in the following. (a) Dynamical states in
symmetric subspace: The solution of the corresponding dynamical equation in
the symmetric subspace Eq. (3) is found to be expressed as

$\displaystyle x(t)$ $\displaystyle=$
$\displaystyle\frac{{e^{\frac{\bar{\lambda}}{2}t}}\cos[\delta-\frac{1}{2}\psi
t]}{(C-\frac{1}{2\bar{\omega}^{2}}e^{\bar{\lambda}t}[U_{0}+U_{1}\cos\psi
t+U_{2}\sin\psi t])^{\frac{1}{2}}},$ $\displaystyle y(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2\bar{\omega}}\left(-{\epsilon
a}+\psi\,\mathrm{\tan}(\frac{1}{2}\psi\,t-\delta)\right)x(t),$ (5)

where $\psi=\sqrt{4\omega\bar{\omega}-{\epsilon}^{2}a^{2}}$,
$\bar{\lambda}=2\lambda-\epsilon a$,
$\beta_{1}=-\lambda((2+a^{2})\epsilon-2\omega)+a((1+a^{2})\epsilon^{2}+\epsilon\omega-2\omega^{2})$,
$\beta_{2}=(-a\lambda+\epsilon(1+a^{2})-\omega)\psi$ and
$\bar{\omega}=\omega-\epsilon$. $C$ and $\delta$ are integration constants.
The other constants $U_{0}$, $U_{1}$ and $U_{2}$ are

$\displaystyle
U_{0}=\frac{2(\epsilon-2\omega)\bar{\omega}}{\bar{\lambda}},\quad
U_{1}=\frac{\epsilon(\beta_{1}\cos
2\delta+\beta_{2}sin2\delta)}{2(\lambda^{2}-\lambda\epsilon
a+\omega\bar{\omega})},$ $\displaystyle U_{2}=\frac{\epsilon(\beta_{1}\sin
2\delta-\beta_{2}cos2\delta)}{2(\lambda^{2}-\lambda\epsilon
a+\omega\bar{\omega})}.$ (6)

Further, we can find that the solution of Eq. (5) will be periodic when
$4\omega\bar{\omega}>{\epsilon}^{2}a^{2}$. In this case, we can write the
state variables of $x_{k}(t)$ and $y_{k}(t)$, $k=1,2$, in the asymptotic limit
($t\rightarrow\infty$) as

$\displaystyle
x_{1}(t)=\frac{\cos(\delta-\frac{1}{2}\psi\,t)}{\left(U_{0}+U_{1}\cos\psi
t+U_{2}\sin\psi t\right)^{\frac{1}{2}}},$ $\displaystyle
y_{1}(t)=\frac{-1}{2\bar{\omega}}\frac{({\epsilon
a}\cos(\delta-\frac{1}{2}\psi\,t)\,+\psi\,sin(\delta-\frac{1}{2}\psi\,t))}{\left(U_{0}+U_{1}\cos\psi
t+U_{2}\sin\psi t\right)^{\frac{1}{2}}},$ (7)

with $x_{2}=x_{1}$ and $y_{2}=y_{1}$. Using the Floquet theory and by
integrating the above equations from $t=0$ to
$t=T=\frac{2\pi}{\sqrt{4\omega\bar{\omega}-{\epsilon}^{2}a^{2}}}$, we
determine the Floquet multipliers $\rho_{i}\,(i=1,...,4)$ and we conclude that
when all the four eigen values lie within the unit circle then the
corresponding periodic orbit is stable. Using the Floquet multipliers we have
obtained the boundary for in-phase synchronized oscillations.

Whenever $4\omega\bar{\omega}<{\epsilon}^{2}a^{2}$, the solution in Eq. (5)
tends to a homogeneous steady state which is given by

$\displaystyle x_{1}^{*}$
$\displaystyle=\mp\frac{1}{\sqrt{2}}\bigg{(}\frac{1}{\hat{a}\epsilon}(-\epsilon
a^{2}\beta_{3}+\bar{\omega}(2\lambda+\psi_{1})-a(\beta_{4}+\lambda\psi_{1}))\bigg{)}^{\frac{1}{2}},$
$\displaystyle y_{1}^{*}$ $\displaystyle=\mp\frac{1}{2\bar{\omega}}(\epsilon
a+\psi_{1})x_{1}^{*},$ (8)

with $x_{2}^{*}=x_{1}^{*}$ and $y_{2}^{*}=y_{1}^{*}$. Other constants are
$\psi_{1}=\sqrt{\epsilon^{2}a^{2}-4\bar{\omega}\omega}$, $\hat{a}=(1+a^{2})$,
$\beta_{3}=(\epsilon a+\lambda+\psi)$,
$\beta_{4}=(\epsilon^{2}+\epsilon\omega-2\omega^{2})$. The eigen values
corresponding to Eq. (8) can be expressed as

$\displaystyle\mu_{1,2}=\frac{1}{2}\,(-\bar{\lambda}-2\psi_{1}\mp\sqrt{\bar{\lambda}^{2}}),\,\,\mu_{3,4}=\frac{1}{2}\bigg{(}-\bar{\lambda}-2\psi_{1}$
$\displaystyle\mp\quad\sqrt{\frac{1}{\hat{a}}(\beta_{5}-4a\omega(2\lambda+\psi_{1})+4(\beta_{6}+\lambda\psi_{1}))}\bigg{)},$
(9)

where $\beta_{5}=a^{3}\epsilon(\epsilon
a-4\lambda)+a^{2}(\epsilon^{2}+4\lambda^{2})$ and
$\beta_{6}=\lambda^{2}+(\epsilon-2\omega)\omega$. The stable region of NAD
emerges when
$\omega<\frac{1}{2}(\epsilon+\sqrt{\epsilon^{2}+a^{2}\epsilon^{2}})$. It is
also clear that when varying the coupling strength $\epsilon$, there exists a
direct transition from stable IPS to NAD state. The analytical boundaries for
IPS and NAD states get well fitted with the numerically obtained boundaries
(see Fig. 3).

(b) Dynamical states in the anti-symmetric subspace: In this section, we have
also found the solution for the dynamical equation in the anti-symmetric
subspace, Eq. (4), which can be written as

$\displaystyle
x(t)=\frac{e^{\frac{\bar{\lambda}}{2}t}\cos(\delta-\frac{\theta}{2}t)}{\left(C+e^{\bar{\lambda}t}(V_{0}+V_{1}\cos(2\delta-\theta
t)+V_{2}\sin(2\delta-\theta t))\right)^{\frac{1}{2}}},$ $\displaystyle
y(t)=\frac{1}{2\omega}\left(-{\epsilon
a}+\theta\,\mathrm{tan}(\frac{1}{2}\theta t-\delta)\right)x(t),$ (10)

where $\theta=\sqrt{4\omega^{2}-{\epsilon}^{2}a^{2}}$. $C$ and $\delta$ are
integration constants. The other constants $V_{0}$, $V_{1}$ and $V_{2}$ are

$\displaystyle V_{0}=\frac{2}{2\lambda-a\epsilon},\quad V_{1}=\frac{\epsilon
a(\lambda\epsilon
a-\epsilon^{2}a^{2}+2\omega^{2})}{4\omega^{2}(\lambda^{2}-\lambda\epsilon
a+\omega^{2})},$ $\displaystyle V_{2}=\frac{\epsilon a(\lambda-\epsilon
a)\theta}{4\omega^{2}(\lambda^{2}-\lambda\epsilon a+\omega^{2})}.$ (11)

The solution given by Eq. (10) is periodic when
$4\omega^{2}>{\epsilon}^{2}a^{2}$. In the asymptotic limit
($t\rightarrow\infty$), the state variables $x_{k}$ and $y_{k}$, $k=1,2,$ can
be written as

$\displaystyle x_{1}(t)=\frac{\cos(\delta-\frac{1}{2})\theta
t}{\left(V_{0}+V_{1}\cos(2\delta-\theta t)-V_{2}\sin(2\delta-\theta
t)\right)^{\frac{1}{2}}},$ $\displaystyle
y_{1}(t)=\frac{-1}{\omega}\frac{({\epsilon a}\cos(\delta-\frac{1}{2}\theta
t)+\theta\sin(\delta-\frac{1}{2}\theta
t))}{\left(V_{0}+V_{1}\cos(2\delta-\theta\,t)+V_{2}\sin(2\delta-\theta\,t)\right)^{\frac{1}{2}}},$
(12)

with $x_{2}=-x_{1}$ and $y_{2}=-y_{1}$. Note that the solution given in Eq.
(12) is periodic with respect to the period t=T=
$\frac{2\pi}{\sqrt{4\omega^{2}-{\epsilon}^{2}a^{2}}}$. Then using the Floquet
multipliers, we have figured out the boundary of stable OPS state.

Further, the solution in Eq. (10) tends toward a steady state when
$4\omega^{2}<{\epsilon}^{2}a^{2}$, which can be expressed as

$\displaystyle x_{1}^{*}$
$\displaystyle=\mp\frac{1}{\sqrt{2}}(\lambda-\epsilon a+\frac{1}{\epsilon
a}(2\omega^{2}-\lambda\theta_{1})+\theta_{1})^{\frac{1}{2}},$ $\displaystyle
y_{1}^{*}$ $\displaystyle=\pm\frac{1}{2\omega}(\epsilon
a+\theta_{1})x_{1}^{*}$ (13)

with $x_{2}^{*}=-x_{1}^{*}$ and $y_{2}^{*}=-y_{1}^{*}$.
$\theta_{1}=\sqrt{\epsilon^{2}a^{2}-4\omega^{2}}$, we have studied the
stability of such anti-symmetric steady states Eq. (13) too and found that the
eigenvalues corresponding to the states are given by

$\displaystyle\mu_{1,2}=\frac{1}{2}\,(-\bar{\lambda}-2\theta_{1}\mp\sqrt{\bar{\lambda}^{2}}),\,\,\,\mu_{3,4}=\frac{1}{2}\,(-\bar{\lambda}-2\theta_{1}$
$\displaystyle\mp\sqrt{\frac{1}{a}(\epsilon a^{2}(\epsilon
a-4\lambda)+4(\lambda^{2}-2\epsilon\omega)+4\omega(2\lambda+\theta_{1})}),$
(14)

where $\bar{\lambda}=(\lambda-\epsilon a)$. From the above eigen values we
have found the stability condition for oscillation death state. When
$\omega<\frac{a\epsilon}{2}$, the stable OD region emerges. The obtained
analytical stability curves of OPS and OD regions exactly matches with the
numerically obtained boundaries (see Fig. 3). From this, we have verified the
stability curves of observed dynamical states of IPS, OPS, NAD and OD regions
analytically too.

##  8 Conclusion

We have investigated the dynamical behavior in a pair of conjugately coupled
Stuart-Landau oscillators. Interestingly, we found the emergence of
spontaneous symmetry breaking due to breaking of permutational symmetries. We
also noticed the coexistence of distinct dynamical states such as SSB-NAD,
NAD-OD, and IPS-OPS states in the parametric space. Through a basin of
attraction study, the dynamical transition to OPS state was illustrated in the
symmetry breaking region. In particular, we report that the transition to OPS
was found as smooth and sudden transition as a function of coupling strength.
Further we have analyzed the effect of feedback factor on the SSB, NAD and OD
states. We found the decreasing feedback factor completely shrinks the
symmetry breaking oscillation state and oscillation death state but it does
not completely suppress the nontrivial AD state. Finally, we have also deduced
the analytical expressions for the observed in-phase and out-of-phase
oscillations as well as nontrivial amplitude death and oscillation death
states. We also find that the feedback factor entirely shrinks oscillation
death states and nontrivial AD states sustained even in the absence of
feedback factor. The obtained results are having good agreement with laser
applications where we need to suppress unwanted fluctuations and attain
constant output. We believe that the present study will shed light on the
dynamics of symmetry breaking states and control of such dynamical behavior,
which will be helpful in many physical and biological systems.

###### Acknowledgements.

KP and VKC are supported by SERB-DST Fast Track scheme for young scientists
under Grant No. YSS/2014/000175. KS sincerely thanks the CSIR for fellowship
under SRF Scheme (09/1095(0037)/18-EMR-I). The work of M.L. is supported by a
DST- SERB Distinguished Fellowship program.

## References

  * [1] Pikovsky A. S., Rosenblum M. G., and Kurths J. Synchronization. A Universal Concept in Nonlinear Sciences Cambridge University Press, Cambridge, UK 2001
  * [2] Tinsley M. R., Nkomo S., and Showalter K. Nat. Phys.8 2012662.
  * [3] Williams C. R. S., Murphy T. E., Roy R., Sorrentino F., Dahms T., and Schöll E. Phys. Rev. Lett. 110 2013064104.
  * [4] Pecora L. M., Sorrentino F., Hagerstrom A. M., Murphy T. E., and Roy R. Nat. Commun. 52014 4079.
  * [5] Aronson D. G., Ermentrout G. B., and Kopell N. Physica D 411990403.
  * [6] Mirollo R. E. and Strogatz S. H. J. Stat. Phys.60 1990245.
  * [7] Banerjee T. and Ghosh D. Phys. Rev. E892014 062902; Phys. Rev. E 89 (2014) 052912.
  * [8] Zou W., Senthilkumar D. V., Duan J., and Kurths J. Phys. Rev. E902014032906.
  * [9] Matthews P. C. and Strogatz S. H. Phys. Rev. Lett.651990 1701.
  * [10] Dolnik M. and Epstein I. R. Phys. Rev. E 5419963361.
  * [11] Hohl A., Gavrielides A., Erneux T., and Kovanis V. Phys. Rev. Lett. 7819974745.
  * [12] Sawai S., Maeda Y., and Sawada Y. Phys. Rev. Lett.85 20002212.
  * [13] Pismen L. M. J. Chem. Phys. 10119943135.
  * [14] Singh R., Menon S. N., and Sinha S. Sci. Rep. 62016 22074.
  * [15] Sathiyadevi K., Karthiga S., Chandrasekar V. K, Senthilkumar D. V, Lakshmanan M. Phys. Rev. E952017 042301.
  * [16] Hynne F., and Graae Soerensen P. J. Phys. Chem.91 19876573.
  * [17] Koseska A., Volkov E., Kurths J. Phys. Rep. 5312013173.
  * [18] Zou W., Senthilkumar D. V.,Nagao R., Kiss I. Z., Tang Y., Koseska A., Duan J., and Kurths J. Nat. Commun.62015 7709.
  * [19] Chen J., , Liu W., Zhu Y., and Xiao J., Euro. Phys. Lett.115201620011
  * [20] Saxena G., Prasad A. and Ramaswamy R. Phys. Rep.5212012 205.
  * [21] Karnatak R., Ramaswamy R. , and Prasad A. Phys. Rev. E76 2007 035201.
  * [22] Ramana Reddy D. V., Sen A. and Johnston G. L. Phys. Rev. Lett.8019985109; Phys. Rev. Lett.85 2000 3381;
  * [23] Strogatz S. H. Nature, London 3941998316.
  * [24] Prasad A. Phys. Rev. E722005 056204; Vicente R., Tang S., Mulet J., Mirasso C. R., and Liu J. M. Phys. Rev. E732006047201.
  * [25] Bera B.K., Hens C. R., Ghosh D. Phys. Lett. A38020162366-2373.
  * [26] Mirollo R. E. and Strogatz S. H. J. Stat. Phys. 60 1990 245.
  * [27] Konishi K. Phys. Rev. E682003067202.
  * [28] Majhia S., Bera B. K., Bhowmick S. K. Ghosh. D. Phys. Lett. A 3802016 3617-3624.
  * [29] Kim M.Y., Roy R., Aron J.L., Carr T.W., Schwartz I.B. Phys. Rev. Lett.942005088101.
  * [30] Prasad A., Lai Y. C., Gavrielides A., and Kovanis V. Phys. Lett. A 318200371.
  * [31] Zhai Y., Kiss I.Z., Hudson J.L. Phys. Rev. E692004026208.
  * [32] Koseska A., Volkov E., Kurths J. Phys. Rev. Lett. 1112013 024103.
  * [33] Zou W., Senthilkumar D. V., Tang Y., and Kurths J. Phys. Rev. E 862012 036210.
  * [34] Zou W., Senthilkumar D. V., Duan J., and Kurths J. Phys. Rev. E 902014 032906.
  * [35] Han W., Cheng H., Dai Q., Li H., Ju P., Yang J. Commun Nonlinear Sci. Numer. Simulat.39 2016 73.
  * [36] Hens C. R., Olusola O. I., Pal P., and Dana S. K. Phys. Rev. E882013034902.
  * [37] Koseska A., Volkov E., Zaikin A., and Kurths J. Phys. Rev. E752007 031916.
  * [38] Ullner E., Zaikin A., Volkov E. I., and García-Ojalvo J. Phys. Rev. Lett.992007 148103.; Ullner E., Koseska A., Kurths J., Volkov E., Kantz H., and García-Ojalvo J., Phys. Rev. E78 2008 031904.
  * [39] Goto Y. and Kaneko K. Phys. Rev. E882013 032718.
  * [40] Suárez-Vargas J. J., González J. A., Stefanovska A., and McClintock P. V. E. Euro. Phys. Lett.85200938008.
  * [41] Deng T., Liu W., Xiao J., and Kurth J. Choas262016094813
  * [42] Chandrasekar V. K., Karthiga S., and Lakshmanan M. Phys. Rev. E 922015 012903.
  * [43] Majhi S., and Ghosh D. Euro. Phys. Lett.1182017 40002.
  * [44] Kim M.Y., Roy R., Aron J.L., Carr T.W., Schwartz I.B. Phys. Rev. Lett.942005088101.
  * [45] Karnatak R., Ramaswamy R., Feudel U. Chaos Solitons Fract.68201448.
  * [46] Zhao N., Sun Z., Yang X. and Xu. W. Euro. Phys. Lett.1182017 30005.
  * [47] Zhao N., Sun Z., Yang X., and Xu. W. Phys. Rev. E 972018 062203.
  * [48] Rakshit S., Bera B. K., Majhi S., Hens C.R., and Ghosh D. Scientific Reports7201745909.
  * [49] Schmidt L., Schönleber K., Krischer K., and García-Morales V. Chaos 242014 013102.
  * [50] Kundu S., Majhi S., Karmakar P., Ghosh D., and Rakshit B. Euro. Phys. Lett.1232018 30001.


††thanks: contributed equally††thanks: contributed equally

# Analysing four major air pollutants’ statistics in Europe - A
superstatistical approach with q-exponential distribution

Hankun He  School of Mathematical Sciences, Queen Mary University of London,
London E1 4NS, United Kingdom  Benjamin Schäfer  School of Mathematical
Sciences, Queen Mary University of London, London E1 4NS, United Kingdom
Institute for Automation and Applied Informatics, Karlsruhe Institute of
Technology, 76344 Eggenstein-Leopoldshafen, Germany  Faculty of Science and
Technology, Norwegian University of Life Sciences, 1432 Ås, Norway  Christian
Beck  School of Mathematical Sciences, Queen Mary University of London, London
E1 4NS, United Kingdom  The Alan Turing Institute, London NW1 2DB, United
Kingdom

# Systematic statistical analysis of air pollution

Hankun He  School of Mathematical Sciences, Queen Mary University of London,
London E1 4NS, United Kingdom  Benjamin Schäfer  School of Mathematical
Sciences, Queen Mary University of London, London E1 4NS, United Kingdom
Institute for Automation and Applied Informatics, Karlsruhe Institute of
Technology, 76344 Eggenstein-Leopoldshafen, Germany  Faculty of Science and
Technology, Norwegian University of Life Sciences, 1432 Ås, Norway  Christian
Beck  School of Mathematical Sciences, Queen Mary University of London, London
E1 4NS, United Kingdom  The Alan Turing Institute, London NW1 2DB, United
Kingdom

# Spatial heterogeneity of air pollution statistics

Hankun He  School of Mathematical Sciences, Queen Mary University of London,
London E1 4NS, United Kingdom  Benjamin Schäfer  School of Mathematical
Sciences, Queen Mary University of London, London E1 4NS, United Kingdom
Institute for Automation and Applied Informatics, Karlsruhe Institute of
Technology, 76344 Eggenstein-Leopoldshafen, Germany  Faculty of Science and
Technology, Norwegian University of Life Sciences, 1432 Ås, Norway  Christian
Beck  School of Mathematical Sciences, Queen Mary University of London, London
E1 4NS, United Kingdom  The Alan Turing Institute, London NW1 2DB, United
Kingdom

###### Abstract

Air pollution is one of the leading causes of death globally, and continues to
have a detrimental effect on our health. In light of these impacts, an
extensive range of statistical modelling approaches has been devised in order
to better understand air pollution statistics. However, the time-varying
statistics of different types of air pollutants are far from being fully
understood. The observed probability density functions (PDFs) of
concentrations depend very much on the spatial location and on the pollutant
substance. In this paper, we analyse a large variety of data from 3544
different European monitoring sites and show that the PDFs of nitric oxide
($NO$), nitrogen dioxide ($NO2$) and particulate matter ($PM10$ and $PM2.5$)
concentrations generically exhibit heavy tails and are asymptotically well
approximated by $q$-exponential distributions with a given width parameter
$\lambda$. We observe that the power-law parameter $q$ and the width parameter
$\lambda$ vary widely for the different spatial locations. For each substance,
we find different patterns of parameter clouds in the $(q,\lambda)$ plane.
These depend on the type of pollutants and on the environmental
characteristics (urban/suburban/rural/traffic/industrial/background). This
means the effective statistical physics description of air pollution exhibits
a strong degree of spatial heterogeneity.

## Introduction

Air pollution is among the highest contributors to premature death and disease
worldwide, causing a significant number of deaths from stroke, lung cancer and
heart diseases [1, 2]. Besides human health, air pollution affects vegetation,
natural ecosystems, climate change, the built environment and subsequently the
economy [3]. In Europe, air pollution is the single largest environmental
health risk [4, 5, 6] and its long-term effects are very serious. Although
emissions and ambient concentrations have fallen steadily in Europe over the
past few decades, as stated in [7, 8], many European countries still exceed
European Union’s (EU) standards [7] and World Health Organization’s (WHO)
guidelines, see [9] for the levels of air pollutants in 2018. Two key air
pollutants, namely, particulate matter ($PM$) and nitrogen oxides ($NOx$),
pose a considerable threat to the health of citizens. About $55000$ and
$417000$ premature deaths in 41 European countries in 2018 were attributed to
$NO2$ and $PM2.5$, respectively [3]. Pollutants such as ozone ($O3$), sulfur
dioxide ($SO2$) and carbon monoxide ($CO$) negatively affect human health as
well. Particularly, ground level $O3$ has been estimated to have caused
$20600$ premature deaths in Europe in 2018; this yearly number has risen by
$20\%$ since 2009 [3]. In this paper we focus on two of the most dangerous
pollutants, $NOx$ and $PM$, but the methodologies presented in our paper can
be similarly applied to other substances.

The impact of air pollution on health does not only depend on the pollutant
type but also on the type of surrounding environment, i.e. people living next
to traffic-heavy roads or industries face higher exposure to air pollution.
The EU [10] uses environmental surrounding types to classify air quality
monitoring sites into traffic, industrial, background, urban, suburban and
rural, based on predominant emission sources and building density. From a
policy perspective, this allows for evaluating the effectiveness of measures
targeting specific emissions sectors and assessing the impact of those
associated pollutants which dominate the area surrounding a given monitoring
station, such as for example traffic or industry (or their absence). Despite
progress made by EU and UK policies [7, 8] addressing all sectors to reduce
emissions and protect citizens from pollutants, meeting the emission reduction
commitments by 2030 remains a challenge [11]. Similar challenges to satisfy
given policies are also existing in other parts of the world.

Having a thorough understanding of the time-varying statistics, i.e. of the
entire probability density function (PDF) of air pollution, is crucial for
policymakers involved in defining thresholds or reducing overall exposure to
air pollution. It is also crucial for the construction of suitable statistical
physics models. PDFs such as gamma, log-normal and Weibull distributions [12]
have been widely used for fitting air pollutant concentration data. However,
these distributions decay approximately like exponential functions at large
values, while earlier investigations have found heavy tails in air pollution
statistics [13], which are not well-described by the above distributions. Some
recent studies [14, 15, 16] have explored the COVID-19 lockdown effects on air
quality (in Europe and in megacities such as Delhi), focusing on comparing the
PDFs or given moments of the PDFs before and during the lockdown.
Superstatistical methods, originating from turbulence modelling [17] and
applied to many fields [18, 19], offer a powerful effective approach to
describe the dynamics of air pollution assuming the existence of well-
separated time scales [13]. Air pollutants such as $NOx$ have been dealt with
success in the superstatistical approach, taking into account nonequilibrium
situations with fluctuating variance parameters [13]. However, this approach
has been verified for limited data sets only, chosen from the UK (London), and
also only for a limited set of pollutants, mainly $NO$ and $NO2$. On a
European scale, and for much larger data sets, it remains unclear whether
heavy tails are generically observed and whether an effective superstatistical
description is applicable. This is the topic of this paper.

The above consideration leads us to a problem that is of general interest for
statistical physics approaches to environmental science. Can we apply standard
methods from nonequilibrium statistical physics, such as the above-mentioned
superstatistical methods, to environmentally relevant time series of pollution
concentrations, and if yes, how can we extract the corresponding
superstatistical parameters from the time series? And how spatially
heterogeneous are the observed results? What are the values of the relevant
superstatistical parameters for different air pollutants and different
geographical environments? Furthermore, what are typical distributions of
observed PDF fitting parameters for the large number of sites distributed
across the European continent? These important types of questions relating to
large ensembles of different measuring stations will be dealt with in the
following.

The paper is organized as follows. First, we introduce our large data set
involving 3544 measuring sites. Next, we investigate the relation between mean
and standard deviation for the observed PDFs to clarify if the PDFs can be
approximated by simple exponential distributions or if more complicated
functions are needed. We then systematically investigate the PDFs of all
sites, in particular the tail behaviour, and show that the tails are generally
much better described by $q$-exponential functions with a given width
parameter $\lambda$ than by functions such as exponential and log-normal,
meaning there is generically power-law decay. Subsequently, we use the maximum
likelihood estimation method (MLE) to extract the $q$ and $\lambda$ parameters
for the best-fitting $q$-exponential distribution, and present plots of
scattered points in the $(q,\lambda)$ plane which exhibit interesting patterns
for our large number of spatial locations investigated. Our main result is
that air pollution statistics is extremely heterogeneous, with the local
variations of best-fitting parameters spanning many orders of magnitudes. Our
investigation is the first one that investigated this in a systematic way for
very large ensembles of different measuring stations. We show that there is a
complex pattern structure in the 2-dimensional ($q$, $\lambda$) parameter
space that depends both on the pollutant type as well as on the classification
type of the local surroundings.

## Results

### The data set considered

In this paper, we aim to conduct a large-scale statistical analysis of air
pollutants, typically $NOx$ and $PM$, on a European scale. Technically, we
access our air quality monitoring data from a large number of locations in
Europe through the interface "Saqgetr", which is an R package available on the
Comprehensive R Archive Network (CRAN) [20]. The vast majority of the
accessible data are openly available from the European Commission’s Airbase
and air quality e-reporting (AQER) repositories [21, 22]. To utilise the data
efficiently, they have been processed into a harmonised form with consistent
and careful treatment of the observations and metadata by Stuart K. Grange
[20]. Concentration level readings and site environment types are the two key
quantities we investigate. We import 9698 locations’ data throughout Europe
within the time span of January 2017 to December 2021, recorded at 1-hour
intervals. To minimise the influences of seasonal fluctuations in time series,
we eliminate sites whose data are too short, typically less than one year. We
also exclude sites where a high percentage of measurements falls below the
detection limit, since sites with clean air are not our primary analysis goal.
Furthermore, we filter out sites whose data are corrupted, the used code and
further details are described in the Methods section. We arrive at 3544 sites
with data that meet our criteria before we proceed with our statistical
analysis. Each data set contains at least $8760~{}(24\times 365)$ data points,
up to about $43800~{}(5\times 24\times 365)$ if the full 5 year period is
available.

To provide a general overview of our analysed data, we show all the data
sites’ locations in Fig. 1(a) as well as an example time series of a selected
site: Bahnhofstrasse, Weiz in Austria for illustration purposes. Measured
concentration time series and histograms are shown in Fig.1(b)-(e), for $NO$,
$NO2$, $PM10$, and $PM2.5$. $NOx$ and $PM$ show seasonal cycles, i.e. during
winter higher pollutant concentrations are more common. We also observe that
for this example site the probability density of $NO$ decays at a slower rate
to zero than those of the other three pollutants. Apparently, typical
distributions exhibit some heavy tails, which we will analyze in much more
detail in the following.

Figure 1:  (a) Illustration of the available data sites on a map of Europe,
with the red circle labeling our example site: Bahnhofstrasse, Weiz, Austria.
Measured time series are shown in (b) and (d), the corresponding probability
densities in (c) and (e). All pollutants display clear seasonality in their
trajectories.

Instead of considering the full distribution of each site’s pollutants, in the
following we will concentrate onto the tails. One reason for doing so is that
we are particularly interested in the statistics of high pollution states,
which are most damaging and described by the tails of the distribution.

With data sets from 3544 air pollution monitoring sites we require an
efficient and context-based automated approach of analysing the sites, details
are described in the Methods section. Refs [7, 10] give more details on macro-
and micro-scale sampling. In brief, stations are divided into three
categories: traffic, industrial, and background based on predominant emission
sources; the surrounding areas are classified as urban, suburban, or rural
based on the density/distribution of buildings. Station types are combined
with area types to provide an overall station classification, and we analyse
our data conditioned on this station classification. We use the following
definitions: "urban traffic": a site located in close proximity to a busy road
in a continuously built-up urban area; "suburban/rural industrial": a site
whose pollution level is influenced predominantly by emissions from an
industrial area or an industrial source in largely built-up or remote areas;
"rural background": a site whose pollution level is influenced by the combined
contribution from all sources upwind of the station and not in built-up areas.
Eventually, seven environmental area types "urban traffic", "suburban/rural
traffic", "urban background", "suburban background", "rural background",
"urban industrial" and "suburban/rural industrial" are used in our statistical
analysis.

### Checking Exponentiality

Let us first consider the simplest possible hypothesis, namely that pollution
concentrations follow an exponential distribution. In this case the PDF is
given by

$f_{\lambda}(x)=\begin{cases}\lambda e^{-\lambda x}&x\geq 0\\\
0&x<0\end{cases}.$ (1)

For exponential distributions one has the general fact that

$\mbox{mean}=\mbox{standard deviation}=\frac{1}{\lambda}.$ (2)

Thus, for each of our 3455 measuring stations we can easily test the
hypothesis of an exponential distribution by plotting mean versus standard
deviation for the measured data. If pollutants were to follow an exponential
distribution, we would expect a clustering along the diagonal in such a plot.
Stations with larger $\lambda$ (smaller mean and variance) would correspond to
cleaner air, and are expected to be found closer to the origin (0,0) as
compared to highly polluted locations. Our results are shown in Fig. 2.

Figure 2:  For each of the pollutants and measuring stations, we plot mean
versus standard deviation. The area type surrounding the measuring station is
color-coded. Data do not follow an exponential distribution, as evidenced by
the fact that the majority of dots do not fall onto the diagonal lines.
Different patterns are observed for the four different substances $NO$, $NO2$,
$PM2.5$, $PM10$. Green colors (rural stations with clean air) cluster near the
origin. However, clustering of the same color patches is observed to be
stronger for $NOx$ as compared to $PM$.

The majority of the points are clustered above (b-d) or below (a) the diagonal
lines, indicating deviations from an exponential distribution. These deviation
patterns are different for each of the four substances.

Apparently, the PDFs of $NO$ and $NO2$ are very different, as the points
scatter mainly below (NO) and mainly above (NO2) the diagonal. The scattering
plots for $PM2.5$ and $PM10$ are more centered around the diagonal, but there
are some unusual $PM2.5$ states with large standard deviation and low means.

Generally, connected clusters of the same colors are more pronounced for $NO$
and $NO2$ as compared to $PM$. This could be attributed to the long-range
transport of the $PM$-particles by moving air. Patterns and spatiotemporal
scales for $PM$ have been extensively discussed previously in [23, 24, 25,
26]. The transport depends on weather conditions and removes the memory to the
site where the particles were originally produced. Thus the colors are more
mixed in our plots. As the weather patterns and the local meteorological
features each contribute to the transport of $PM$-particles, the measuring
site type has limited impact on the observed PDFs of $PM$.

### Fitting power-law tails for the data

As exponential tails apparently do not fit the data well, as illustrated by
the deviations from the diagonal in Fig.2, we now propose a different fitting
function, motivated by many previous investigations in generalized versions of
statistical mechanics [27]. This is a fitting by a so-called $q$-exponential,
which asymptotically decays with a power-law exponent $-\frac{1}{q-1}$. The
normalized PDF is defined as follows:

$f_{q,\lambda}(x)=(2-q)\lambda[1-\lambda(1-q)x]^{\frac{1}{1-q}}\text{ for
}1-\lambda(1-q)x\geq 0,x>0,$ (3)

where $q$ is the entropic index [27, 28, 29], $\lambda$ is a positive width
parameter and $x$, in our case, denotes the air pollutant concentration. Eq.
(3) contains the exponential distribution as a special case, namely for $q=1$,
as the $q$-exponential function, defined as
$e_{q}(x)=[1+(1-q)x]^{\frac{1}{1-q}}$, converges to the exponential function
in the limit $q\rightarrow 1$. For $q<1$, $f_{q,\lambda}(x)$ lives on a finite
support and becomes exactly zero above a critical value x, since, by
definition, $e_{q}(x)=0$ for $1-\lambda(1-q)x<0$. In contrast, if $q>1$,
$1-\lambda(1-q)x>0$, then eq. (3) exhibits power-law asymptotic behavior.

The occurrence of $q$-exponentials with $q>1$ in PDFs of complex systems is
very well-motivated by superstatistical models [30, 31]. In these types of
models, one assumes a temporally fluctuating parameter $\lambda$ for local
exponential distributions as given in eq.(1). These fluctuations of $\lambda$
take place on a long time scale, much longer than local air pollution
concentration fluctuations. The marginal distribution, obtained by integration
over all possible values of $\lambda$, and describing the long-term behaviour
of the air pollution concentration dynamics, is then a $q$-exponential, with

$q=\frac{\langle\lambda^{2}\rangle}{\langle\lambda\rangle^{2}}.$ (4)

Here $\langle\cdots\rangle$ denotes the expectation with respect to the PDF of
$\lambda$, see [30] for more details. Strictly speaking, a $q$-exponential is
only obtained exactly if $\lambda$ is $\Gamma$-distributed, but the general
idea of superstatistics is that a parameter $q$ can be defined by eq.(4) for
more general distributions different than the $\Gamma$ distribution as well.
The concept that wind changes and other effects (such as traffic fluctuations)
can lead to a superstatistical dynamics for pollution concentrations was first
worked out in [13], where further details can be found, in that case for the
special example of $NO$ and $NO2$ concentrations as measured in London. Our
investigation here is much more general, as we include data of thousands of
measuring stations, and also investigate $PM2.5$ and $PM10$ concentrations.

Figure 3:  Best-fitting parameters of $q$-exponentials. We observe an
increasing trend of $q$ versus $\log\lambda$ for $NO$ (a) and $NO2$ (b),
whereas a more disk-shaped pattern is observed for $PM2.5$ (c) and $PM10$ (d).
The environmental characterizations of the measuring stations are again
encoded by colors. Again we observe correlated patches of a single given color
for $NO$ and $NO2$, where for $PM2.5$ and $PM10$ the patterns is more mixed.

For all our 3544 measuring stations we extract histograms of the pollution
concentration from the measured time series, and determine the best-fitting
parameters $q$ and $\lambda$ for the given data set. More details on the
numerical procedure are described in the Method section. Our results are shown
in Fig. 3.

A truly surprising result of our analysis is the fact that we observe an
immensely large range of values of the parameter $\lambda$ for the best-
fitting $q$-exponential as given in eq.(4) for the various measuring stations.
Note the logarithmic scale of the plots, the parameter $\lambda$ can take on
values as small as $10^{-2}$ up to values as large as $10^{2}$, which spans
four orders of magnitude. Typical $q$-values are in the range $0.8-1.4$, but
there are subtle differences between the various substances, with $NO$
reaching large $q$-values such as 1.6, and $NO2$ reaching small $q$-values
such as $0.6$ in the scattering plots. Also the shape of the scattering cloud
of points is different for the different substances. For example, the typical
range of $\lambda$ for $PM2.5$ and $PM10$ varies only by a factor 10, whereas
for $NO$ and $NO2$ it varies by a factor $10^{3}$. The scattering plot data
look more spherically symmetric for $PM2.5$ and $PM10$, as compared to $NO$
and $NO2$.

Fig. 3 panels a and b indicate a roughly linear approximate relationship
between $q$ and $\log\lambda$ for $NOx$, with most of the points corresponding
to traffic clustered at the left, while urban/suburban background points are
in the middle part, and rural background points are scattered widely at the
right. The PDF decay rate increases as $\lambda$ increases from highly
polluted urban traffic sites to less polluted rural background areas. For
$PM2.5$ and $PM10$, there is a different weak uphill relationship between $q$
values and $\lambda$, as can be seen in Fig.3 panels c and d. The urban
background points are reaching small $\lambda$ values such as $0.01$ for
$PM10$, and suburban/rural traffic, suburban/rural industrial and rural
background points cluster on the right hand side with large $\lambda$. The
attained range of $\lambda$ values is smaller as compared to the case of NOx,
and the shape of possible values $(q,\lambda)$ as displayed in the Figure is
more spherical. The colors appear to be more randomly mixed.

The stronger colour mixing for $PM2.5$ and $PM10$ can again be interpreted by
the fact that by air movement transport the distributions cannot be uniquely
identified with the original environmental types where the $PM$-particles were
produced. The large scattering of parameters $(q,\lambda)$ shows that for a
given substance at a given environmental type there is not just one possible
distribution, but a large range of possible distributions. These distributions
may also vary in time, according to the weather conditions. Still, our
scattering plots allow us, in principle at least, to select typical parameter
values for a given environmental type. The fact that there is a broad
distribution of parameters is very much in line with the basic modelling
assumption of superstatistics, in this case however applied to a spatial
ensemble of different locations. There is a strong heterogeneity in space,
meaning different spatial measuring locations have quite different PDFs. This
spatial heterogeneity is an second effect, adding to the temporal
heterogeneity of local exponentials, which effectively leads to
$q$-exponentials at individual locations as explained above.

### Spatial distribution plots of $\lambda$-values

Finally, we are interested in the PDFs of $\lambda$ values for our fits of the
various classified locations where the measurements are taken. We compare the
summary statistics (such as distribution, range and quartiles) of $\lambda$
for the four air pollutants with the aid of so-called violin plots, see Fig.4.
Within these, we visualise the distribution of $\lambda$ using density curves,
which corresponds to the approximate frequency of data. In addition to the
density curves, an overlaid box plot is shown with the rectangle depicting the
ends of the first (25%) and third quartiles (75%) and the central dot
indicating the median (50%).

Figure 4:  The violin plots show the distribution of $\lambda$ for seven
environment types as well as the median as a white dot, the interquartile
range as a thick black bar, and the 95% confidence interval as a thin black
bar within the colored violin. The environment types were ranked by medians of
$\lambda$ from lowest to highest. At each y-axis the number of sites evaluated
in each category is reported. In the case of $NO$ (a), a rescaling has been
applied to capture the different scale for rural background (green). Likewise,
there is also a rescaling for suburban/rural industrial (yellow) and rural
background for $NO2$.

An interesting result of the violin plots shown in Fig. 4 is that the
$\lambda$ distributions extend to very large values, as indicated by the long
extensions to the right for area types such as urban background,
suburban/rural industrial and rural background. Additionally, the probability
distributions of $\lambda$ (represented by the "shapes" of the violins)
exhibit nontrivial behaviour for some of the environmental site types. For
example, for $NO$ urban industrial sites there is a rather unusual pattern
with several local maxima and minima. A much more generic pattern, with just a
single broad maximum, is observed for sites which are suburban/rural
industrial, as well as for those with a rural background, and this structure
is there for all 4 different types of pollutants.

Another intriguing result is that the order of the median rankings (from low
to high) for $NO$ and $NO2$ are almost the same, with the exception of a swap
between urban industrial and urban background. For $PM2.5$ and $PM10$ there
are more swaps. Fig. 4 panels c and d show that for $PM2.5$ and $PM10$ the
typical values of $\lambda$ are smaller, below 0.9. A smaller $\lambda$
indicates a more heavily polluted site. Furthermore, we observe in the case of
$PM$ only minor differences in the medians of $\lambda$ for different
environment categories. The reason for this is that the type of environment
has a direct effect on $NOx$ concentrations while they have only a minor
effect on $PM$ since the particles travel and lose the memory of their
environmental category. Nevertheless, suburban/rural traffic and rural
background sites have the largest and second largest medians for $\lambda$,
respectively.

## Discussion and Conclusion

The use of $q$-exponentials for air pollution statistics has been previously
advocated by Williams et al. [13], however that study was special since it was
only looking at locations in Greater London, and the substances investigated
were just $NO$ and $NO2$. In this paper, we have extended the statistical
analysis to a much larger database, taking into account data from 3544
measuring stations, and analysing the statistics of $PM2.5$ and $PM10$, in
addition to $NO$ and $NO2$. Naturally, for this vast amount of data novel
methods needed to be developed for the fitting procedures, and novel graphical
representations (scattering plots of parameter tuples) were used to illustrate
the spatial heterogeneity of the results. Our main findings can be summarized
as follows:

Firstly, we have clear evidence that generically PDFs of pollution
concentrations do not decay in an exponential way. A much better fit is given
by $q$-exponentials, which asymtotically decay as a power-law if $q>1$, with
exponent $-1/(q-1)$. Our analysis complements previous work, in which often
different distributions have been used, including log-normal, gamma and
Weibull distributions. Overall, we find that $q$-exponentials yield a better
fit of the tails. The $q$-exponential is also a very plausible physical model,
since –in the spirit of superstatistics– it simply arises from the
agglomeration of many exponential distributions that have temporal
fluctuations of the effective decay rate. In our investigation we also tested
the other candidate distributions mentioned above and found that they
sometimes yield a good fit of the low-concentration behaviour, close to the
maximum, but not of the tail behaviour, see our code for details.

Secondly, $q$ and $\lambda$, as obtained from the optimum fittings of data
from 3544 measuring stations, exhibit interesting patterns in the
$(q,\lambda)$ plane. The shape of these regions is characteristic for each of
the 4 substances investigated, with big differences between $NO$, $NO2$ and
the $PM$ statistics.

Thirdly, environmental types, i.e. the surroundings of the measuring station,
play an important role. We color-coded these different environments into 7
categories. For $NO$ and $NO2$ each category occupies a typical sub-region in
the $(q,\lambda)$ plane, whereas for $PM2.5$ and $PM10$ the picture is more
mixed. This can be explained by the fact that there is transport by moving air
for the $PM$-particles, so that the memory to the environment of the station
where the actual measurements are done is lost, i.e. these particles may
travel quite a long way and many of them are not produced locally. The
patterns and spatio-temporal scales for $PM$-dynamics have been extensively
examined before in [23, 24, 25, 26], the transport is dependent on the weather
conditions and removes the memory of the original source site.

As a next step, it would be desirable to examine correlations between $PM$ and
$NOx$ concentrations and to include more environmental factors, including wind
speed, wind direction and surface temperature. In light of these more detailed
investigations, a statistical analysis could enable policymakers to produce
more precise rules and thresholds for individual types of environmental
conditions and meteorological conditions, taking into account fluctuations and
extreme events. Our analysis could also be extended to other substances, such
as sulfur oxide, carbon oxide and ozone, besides $NOx$ and $PM$. Taking all
this into account could help policymakers to develop tailored guidelines to
reduce health damage caused by air pollution. We stress again that the
detailed description of the entire pollution spectrum, including the exact
behaviour of the tails, seems critical here to better estimate the risks of
very high pollution situations. Concluding, our analysis shows that there is
strong heterogeneity in the data, and one needs to be careful because the PDFs
vary strongly from one location to another.

## Methods

### Data processing

We import European air pollution data from "Saqgetr"[20], which is an R
package developed by Stuart Grange to analyse air quality data — or more
generally atmospheric composition data. The package provides users with fast
access to thousands of sites’ data from air quality networks, which are
supported by Ricardo Energy & Environment. We import all 9698 sites’ hourly
data between 2017 and 2021 from the package and save them as individual .csv
files. The raw data contain detailed information about the air pollution
monitoring sites and their hourly measurements. From them, we select the
following information:

  1. 1.

Name and site code. Each site has a unique site code for identification and
simplification in coding. (E.g. "gb0050a" for the Rosia Road in Gibraltar)

  2. 2.

Longitude and latitude, which are used to show stations’ locations on maps
(such as Fig. 1)

  3. 3.

$NO$, $NO2$, $PM2.5$ and $PM10$ pollutants’ hourly concentration data

  4. 4.

Station and area types, which are combined for classifying sites

The data sets retrieved with the Saqgetr package contain several problems, for
which our solutions are:

  1. 1.

We filter out sites with less than one year of $NO$, $NO2$, $PM2.5$ and $PM10$
pollution concentration data. We need a minimum number of 8760 data points
($365$ days $\times$ $24$ hours per day) for one pollutant for a meaningful
statistical analysis and to avoid analysing a single season.

  2. 2.

Measurements below the detection limit, including some that are even below
zero, are usually replaced by the detection limit divided by two. We filter
out sites with more than $15\%$ of data below the detection limit. This is
aligned with the recommendations of the US Environmental Protection Agency (US
EPA)[32], which suggests that substitution may be a viable approach when up to
$15\%$ of the data cannot be detected. This steps is further justified as
sites with extremely low pollution are not the main focus of this study.

  3. 3.

We remove sites whose data is repeating in a single measurement at least
$15\%$ of the time. We assume that these sites either lack precision in
measurement or contain too much corrupted data.

### Fitting procedure

The filtered data sets were analysed in Python. First, we plot the probability
density function (PDF) of each site’s pollutant concentration level using a
log scale. Then, we consider the distributions of the higher concentration
tail by finding the maximum of the distribution and filtering out the smaller
concentration distribution (left of the maximum of the distribution). One
reason for doing so is that we are particularly interested in the high
concentration tail behaviour. As an example, let us discuss the probability
densities for low pollution concentrations for $NO$ in Amstetten, Austria
(Fig. 5(a)) and for $PM2.5$ in Riadok, Slovakia (Fig. 5(b)). We determine the
highest point of the kernel line estimate and its corresponding concentration
level (black vertical line). The underlying density distribution of the red
line created an increasing slope which we cut off. We only analyse the
underlying density distributions of the turquoise line and the higher
concentration tail.

Figure 5:  Determination of the lower cutoff at peak density (vertical line).
We display the empirical probability density functions (PDFs) of the $NO$ (a)
and $PM2.5$ (b) concentrations with cut-off options at
$1.5\times\mbox{mean}~{}\mu g/m^{3}$. The empirical smoothed PDF estimates the
densities of the data before (red) and after the peak (turquoise).

We fit the so-obtained data with exponential and $q$-exponential distributions
derived from MLE methods. C. Shalizi [33] and R. Gonzalez-Val [34] described
methods for estimating the parameters of the $q$-exponential and log-normal
distributions respectively using MLE. We carry out the fit, using the
_scipy.curve_fit_ and _scipy.stats.rv_continuous.fit_ functions of the scipy
module[35]. The $q$-exponential distribution using MLE methods generally show
robust fitting results. See code for further details and fits of
non-$q$-exponential functions.

### Data availability

Data sets analyzed in the present study can be obtained via the saqgetr
package from https://github.com/skgrange/saqgetr. They should be downloaded as
csv files and named after their site codes. The code to generate the figures
in the paper, as well as the implementation of the method for the data sets
used in the paper, are available at https://github.com/hurst0415/Spatial-
heterogeneity-of-air-pollution-statistics.

## References

  * [1] World Health Organization (WHO).  How air pollution is destroying our health (2020).  URL https://www.who.int/news-room/spotlight/how-air-pollution-is-destroying-our-health.
  * [2] Shah, A. S. _et al._ Short term exposure to air pollution and stroke: systematic review and meta-analysis.  _BMJ_ 350, h1295 (2015).
  * [3] Ortiz, A. & Guerreiro, C.  Air Quality in Europe - 2020 report (2020).
  * [4] Institute, H. E.  State of global air 2019.  _Data source: Global Burden of Disease Study_ (2019).
  * [5] Organization, W. H. _et al._ Who global ambient air quality database (update 2018).  _World Health Organization: Geneva, Switzerland_ (2018).
  * [6] Collaborators, G. . R. F. _et al._ Global, regional, and national comparative risk assessment of 84 behavioural, environmental and occupational, and metabolic risks or clusters of risks, 1990–2016: a systematic analysis for the global burden of disease study 2016.  _Lancet (London, England)_ 390, 1345 (2017).
  * [7] European Union.  Directive 2008/50/EC of the European Parliament and of the Council of 21 May 2008 on ambient air quality and cleaner air for Europe (2008).
  * [8] European Union.  Directive 2016/2284/EC of the European Parliament and of the Council of 14 December 2016 on the reduction of national emissions of certain atmospheric pollutants (2016).
  * [9] Organization, W. H.  _Air quality guidelines: global update 2005: particulate matter, ozone, nitrogen dioxide, and sulfur dioxide_ (World Health Organization, 2006).
  * [10] European Parliament.  Directive 2011/24/EU of the European Parliament and of the Council of 9 March 2011 on the application of patients’ rights in cross-border healthcare (2011).
  * [11] European Commission.  The clean air programme (2013).  URL http://ec.europa.eu/environment/air/clean_air/index.htm.
  * [12] Lu, H.-C.  Comparisons of statistical characteristic of air pollutants in Taiwan by frequency distribution, journal of the air & waste management.  _Association_ 53:5, 608–616 (2003).
  * [13] Williams, G., Schäfer, B. & Beck, C.  Superstatistical approach to air pollution statistics.  _Physical Review Research_ 2, 013019 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.013019.  Publisher: American Physical Society.
  * [14] Schäfer, B. _et al._ Covid-19 impact on air quality in megacities.  _arXiv preprint arXiv:2007.00755_ (2020).
  * [15] Baldasano, J. M.  COVID-19 lockdown effects on air quality by NO2 in the cities of Barcelona and Madrid (Spain).  _Science of The Total Environment_ 741, 140353 (2020).  URL https://www.sciencedirect.com/science/article/pii/S0048969720338754.
  * [16] Group, E. R. _et al._ The effect of COVID-19 lockdown measures on air quality in London in 2020 (2020).
  * [17] Beck, C.  Statistics of three-dimensional lagrangian turbulence.  _Phys. Rev. Lett._ 98, 064502 (2007).  URL https://link.aps.org/doi/10.1103/PhysRevLett.98.064502.
  * [18] Metzler, R.  Superstatistics and non-gaussian diffusion.  _The European Physical Journal Special Topics_ 229, 711–728 (2020).
  * [19] Beck, C. _et al._ Nonextensive statistical mechanics, superstatistics and beyond: theory and applications in astrophysical and other complex systems (2020).
  * [20] Grange, S. K.  Technical note: saqgetr r package (2019).  URL https://drive.google.com/open?id=1IgDODHqBHewCTKLdAAxRyR7ml8ht6Ods.
  * [21] European Parliament, Council of the European Union.  DIRECTIVE 2008/50/EC OF THE EUROPEAN PARLIAMENT AND OF THE COUNCIL of 21 May 2008 on ambient air quality and cleaner air for Europe (2008).  URL https://eur-lex.europa.eu/legal-content/EN/ALL/?uri=CELEX%3A32008L0050.
  * [22] European Commission.  2011/850/EU: Commission Implementing Decision of 12 December 2011 laying down rules for Directives 2004/107/EC and 2008/50/EC of the European Parliament and of the Council as regards the reciprocal exchange of information and reporting on ambient air quality (notified under document C(2011) 9068) (2011).  URL https://eur-lex.europa.eu/legal-content/EN/TXT/?uri=CELEX%3A32011D0850.
  * [23] Prospero, J. M., Olmez, I. & Ames, M.  Al and Fe in PM 2.5 and PM 10 suspended particles in South-Central Florida: The impact of the long range transport of African mineral dust. (2001).
  * [24] Uno, I. _et al._ Trans-Pacific yellow sand transport observed in April 1998: A numerical simulation. (2001).
  * [25] Bardouki, H. _et al._ Chemical composition of size resolved atmospheric aerosols in the eastern Mediterranean during summer and winter. (2003).
  * [26] Kallos, G., Papadopoulos, A., Katsafados, P. & Nickovic, S.  Trans-Atlantic Saharan dust transport: Model simulation and results. (2006).
  * [27] Tsallis, C.  Introduction to nonextensive statistical mechanics: approaching a complex world.  _Springer_ 1, 2–1 (2009).
  * [28] Hanel, R., Thurner, S. & Gell-Mann, M.  Generalized entropies and the transformation group of superstatistics.  _Proceedings of the National Academy of Sciences_ 108, 6390–6394 (2011).
  * [29] Jizba, P. & Korbel, J.  Maximum entropy principle in statistical inference: Case for non-shannonian entropies.  _Physical Review Letters_ 122, 120601 (2019).
  * [30] Beck, C. & Cohen, E. G. D.  Superstatistics.  _Physica A: Statistical mechanics and its applications_ 322, 267–275 (2003).
  * [31] Beck, C., Cohen, E. G. D. & Swinney, H. L.  From time series to superstatistics.  _Physical Review E_ 72, 056133 (2005).
  * [32] United States. Office of Environmental Information. Quality Staff.  _Guidance for Data Quality Assessment: Practical Methods for Data Analysis EPA QA/G-9, (QA00 Update)_ (U.S. Environmental Protection Agency, 2000).  URL https://books.google.com.hk/books?id=G3DmHAAACAAJ.
  * [33] Shalizi, C. R.  Maximum likelihood estimation for q-exponential (tsallis) distributions (2007).  eprint math/0701854.
  * [34] Gonzalez-Val, R., Ramos-Gutiarrez, A. & Sanz-Gracia, F.  Size Distributions for All Cities: Lognormal and q-exponential functions.  ERSA conference papers ersa11p554, European Regional Science Association (2011).  URL https://ideas.repec.org/p/wiw/wiwrsa/ersa11p554.html.
  * [35] Virtanen, P. _et al._ SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ 17, 261–272 (2020).

### Acknowledgements

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie-Sklodowska-Curie Grant
agreement No 840825. We also gratefully acknowledge funding from a 2022 QMUL
Research England Policy Impact grant awarded to C.B. and funding from the
Helmholtz Association under grant no. VH-NG-1727 to B.S.

### Author contributions

H.H., B.S. and C.B. conceived the project, H.H. and B.S. processed the data,
performed the data analysis and produced the plots. All authors interpreted
the results, wrote and reviewed the manuscript. Overall B.S. and C.B.
contributed equally.

### Competing interests

The authors declare no competing interests.