# Theory of transient chimeras in finite Sakaguchi-Kuramoto networks

Roberto C. Budzinski  James W. C. Graham  Ján Mináč  Lyle E. Muller
Department of Mathematics, Western University, London, ON, Canada  Western
Institute for Neuroscience, Western University, London, ON, Canada  Western
Academy for Advanced Research, Western University, London, ON, Canada

###### Abstract

Chimera states are a phenomenon in which order and disorder can co-exist
within a network that is fully homogeneous. Precisely how transient chimeras
emerge in finite networks of Kuramoto oscillators with phase-lag remains
unclear. Utilizing an operator-based framework to study nonlinear oscillator
networks at finite scale, we reveal the spatiotemporal impact of the adjacency
matrix eigenvectors on the Sakaguchi-Kuramoto dynamics. We identify a specific
condition for the emergence of transient chimeras in these finite networks:
the eigenvectors of the network adjacency matrix create a combination of a
zero phase-offset mode and low spatial frequency waves traveling in opposite
directions. This combination of eigenvectors leads directly to the coherent
and incoherent clusters in the chimera. This approach provides two specific
analytical predictions: (1) a precise formula predicting the combination of
connectivity and phase-lag that creates transient chimeras, (2) a mathematical
procedure for rewiring arbitrary networks to produce transient chimeras.

Chimera states are a central phenomenon in the study of nonlinear dynamics.
Early computer simulations found that pockets of synchronous and asynchronous
nodes could co-exist in a network where all nodes and connections were
otherwise identical [1]. This example prompted extensive theoretical work to
understand how phase synchrony could spontaneously appear in a subset of nodes
without systematic differences in either their intrinsic properties or their
connectivity [2]. Since their initial discovery, chimera states have been
observed in networks of phase oscillators [3, 4, 5, 6], neural networks [7, 8,
9, 10], and other dynamical systems [11, 12, 13, 14]. Chimeras have also been
found in experiments [15, 16, 17, 18, 19, 20, 21], demonstrating their
existence in real-world systems and prompting further theoretical study of
their underlying mechanisms.

The mathematical approaches that have been developed to study chimeras have
been well suited to describe these states in continuous systems [5, 22, 23,
24]. The assumption of the continuum or asymptotic limit has made analytical
insights possible in previous work, such as the identification of bifurcations
[2, 5, 25, 23, 26]. Further, numerical results have studied bifurcations
leading to the emergence of chimera states [27, 28]. Several studies have
analyzed the existence of weak chimeras, which are characterized by frequency
synchronization [29, 30, 31], and also explored these states in minimal
networks [32, 33]. Recent work has developed analytical approaches for strong
chimera states in finite networks, where the synchronized and asynchronous
clusters mutually interact to permanently stabilize the chimera in networks
with two or more modules [34, 35]. Transient states, however, are increasingly
understood to be important in neural dynamics and computation [36], and
chimeras are known to be transient states in finite networks of Sakaguchi-
Kuramoto oscillators [37]. While techniques such as the master stability
function or analysis of eigenvalues of the graph Laplacian have provided
fundamental insight into the dynamics of oscillator networks, it is
increasingly recognized that these techniques have limitations in studying
these transient dynamics related to the brain [38]. At present, the underlying
network mechanisms and specific mathematical conditions for the emergence of
transient chimeras states in finite networks of Sakaguchi-Kuramoto oscillators
remain unknown.

Here, we report a analytical approach for chimera states in finite networks of
Sakaguchi-Kuramoto oscillators. This mathematical approach can be understood
in geometric terms, where combinations of eigenvectors of the adjacency matrix
play a fundamental role. Our results demonstrate that transient chimera states
can be understood as a combination of a zero phase-difference mode and low
spatial frequency wave patterns traveling in opposite directions. This is
valid for finite networks where different chimera patterns are observed. Our
approach also identifies the specific roles played by connectivity and phase-
lag between oscillators in generating transient chimera states.

We use a recently introduced complex-valued, operator-based approach to the
Kuramoto model [39], which is able to reproduce a diversity of synchronization
patterns observed in Kuramoto networks [40]. This approach allows a direct
connection between the eigenvectors describing the system and the shape of the
resulting nonlinear dynamics [40, 41]. We start with the standard equation for
the Sakaguchi-Kuramoto network dynamics

$\dot{\theta}_{i}(t)=\omega+\epsilon\sum_{j=1}^{N}A_{ij}\sin\big{(}\theta_{j}(t)-\theta_{i}(t)-\phi\big{)}\,,$
(1)

where $\theta_{i}(t)$ is the phase of the $i^{\mathrm{th}}$ oscillator at time
$t$, $\omega$ is the natural frequency, $\epsilon$ scales the coupling
strength, $\bm{A}$ is the adjacency matrix, and $\phi$ is the phase-lag. Here,
we consider a distance-dependent power-law connectivity, where the weights
$A_{ij}$ decay with the distance between nodes $i$ and $j$ following a power-
rule with periodic boundary conditions [42]. Without loss of generality, we
consider oscillators with natural frequency of $\omega=10\pi$. We then
introduce a complex-valued system that is related to the original, nonlinear
Kuramoto model and allows connecting the eigenspectrum of the adjacency matrix
to the resulting nonlinear dynamics [39, 40, 41]. Specifically, starting with
unit-modulus initial conditions $|x_{i}(0)|=1\forall i$, with arguments
$\mathrm{Arg}[x_{i}(0)]$ equal to the initial conditions of the Kuramoto
system $\theta_{i}(0)$, we can use two operators to evolve a complex-valued
system:

$\bm{x}(t+\varsigma)=\Lambda\big{[}e^{\mathrm{i}\omega\varsigma}e^{\varsigma\bm{K}}\bm{x}(t)\big{]}\,,$
(2)

where $\varsigma$ is small but finite, $\Lambda$ represents an elementwise
operator mapping the modulus of each state vector element $x_{i}$ to unity,
and $\bm{K}=\epsilon e^{-\mathrm{i}\phi}\bm{A}$. We have shown in recent work
that Eq. (2) precisely reproduces the trajectories of the original Kuramoto
system, specifically by comparing the variables $\theta_{i}(t)$ to
$\mathrm{Arg}[x_{i}(t)]$ [40]. Because all evolution of the arguments
$\mathrm{Arg}[x_{i}(t)]$ is governed by the linear matrix exponential operator
(with the elementwise $\Lambda$ operator only changing the moduli), we can
analyze the nonlinear dynamics in this system in terms of the arguments of the
eigenvectors of $\bm{K}$ [41, 43]. Thus, while connecting the eigenspectrum of
a network’s adjacency matrix to the resulting nonlinear dynamics is a
difficult problem in general, here we demonstrate that this is possible in the
Sakaguchi-Kuramoto oscillator system.

We then study the emergence of transient chimera states using this approach.
When each node in a network is connected by the same connectivity rule, such
as with the complete graph, ring graphs, and distance-dependent networks, we
can obtain the eigenvalues and eigenvectors of the system analytically via the
circulant diagonalization theorem (CDT) [44] (see Supplemental Sec. II). The
CDT states that the eigenvectors are given by
$(\bm{v}_{k})_{s}=(\nicefrac{{1}}{{\sqrt{N}}})\exp\left[(\nicefrac{{-2\pi\mathrm{i}}}{{N}})(k-1)(s-1)\right]$,
which are the columns of the discrete Fourier transform (DFT). Specifically,
the first eigenvector is $\bm{v}_{1}=[1,1,\cdots,1]$, and the arguments of its
elements are $\mathrm{Arg}[\bm{v}_{1}]=[0,0,\cdots,0]$. In the context of our
analysis, this eigenvector represents the zero phase-difference state. In a
network arranged on a one-dimensional ring, the other eigenvectors have a
phase configuration (or, argument of each element) representing waves with
increasing spatial frequency. These eigenvectors appear in pairs that travel
in opposite directions (clockwise and counter-clockwise; cf. inset of Fig.
1c). Here, we quantify the contributions of each eigenvector in terms of the
“eigenmode contribution” $\mu_{i}(t)=\langle\bm{x}(t),\bm{v}_{i}\rangle$,
where $\langle.\rangle$ represents the complex inner product, and $i\in[1,N]$.

Figure 1: A substrate for transient chimera states. (a) We observe the
emergence of chimera states in the network. (b) The eigenmode contribution, at
first, appears to be uniform, however, a close look reveals a different
pattern when the chimera state emerges (see inset). (c) When the network is
asynchronous, at $t=500$, the modes contribution $|\mu|$ is uniform (blue line
with dots); at $t=4500$, when a chimera state is observed, there are specific
modes with higher contribution (red line with squares). These modes are
related to eigenvectors with specific phase configuration (see inset): zero
phase-difference and low spatial frequency waves traveling in opposite
directions. Here, we consider $\epsilon=20$ and $\phi=1.40$.

Figure 1a illustrates a chimera state in a network of Sakaguchi-Kuramoto
oscillators with distance-dependent power-law connectivity [42]. Due to random
initial conditions, the network begins in an asynchronous state. The network
then evolves to exhibit a chimera state, with clusters of synchronous and
asynchronous nodes. Figure 1b shows $|\mu_{i}(t)|$ across time. When the
network is an asynchronous state ($t=500$, blue dotted line, Fig. 1b), the
eigenmode contribution is uniform across modes (blue line and dots, Fig. 1c).
When the network exhibits a chimera ($t=4500$, red dotted line, Fig. 1b),
however, the eigenmode contribution is no longer uniform: the $\mu_{i}(t)$ are
increased for the lowest and highest modes relative to the asynchronous state
(red line and squares, Fig. 1c). This result indicates that, during transient
chimera states, the network exhibits a combination of low-spatial-frequency
modes, traveling in both the clockwise and counter-clockwise directions around
the ring, in addition to the mode representing zero phase offset $\mu_{1}$
(inset of Fig. 1c).

Figure 2: Eigenmode contribution for finite networks. (a) We analyze the
eigenmode contribution $|\mu|$ for the final state of the network ($t=5000$).
For $\phi=1.40$ (red line) chimera states are possible, and the contribution
of the first and last modes increase in the same manner as observed in Fig.
1c. As a control case, we consider $\phi=\nicefrac{{\pi}}{{2}}$ (blue line)
where the network remains in an asynchronous state and the eigenmode
contribution is uniform across modes. These results consider $1000$
realizations with different initial conditions, where we plot the average with
shaded area representing the standard mean error. We also consider different
examples of chimera patterns for $N=501$ and (b) $\phi=1.35$, (c) $\phi=1.40$,
(d) $\phi=1.45$, (e) and $\phi=1.50$, where we observe that, in all cases, the
modes contribution follows the similar fashion as shown in Fig. 1c. Here, we
consider $\epsilon=20$.

This analysis suggests a potential substrate for transient chimera states in
finite Sakaguchi-Kuramoto oscillator networks as a specific combination of
modes. To verify that this is the case, we next studied a set of networks with
increasing numbers of nodes across many different initial conditions. We
consider three distance-dependent networks [42] with size $N=101$, $N=501$,
and $N=1001$ and plot $|\mu|$ for the final state of the system at $t=5000$.
Figure 2a illustrates the case with phase-lag ($\phi=1.40$, red line), where
the network displays chimera states, and with neutral coupling
($\phi=\nicefrac{{\pi}}{{2}}$, blue line), where the network remains in an
asynchronous state. When the network is asynchronous, the $|\mu|$ are uniform,
demonstrating equal contribution across all eigenmodes. When the network
exhibits a chimera state, however, the contribution of the first and last
eigenmodes increases in the same manner as in Fig. 1. This result demonstrates
that the condition identified in Fig. 1 occurs in the average over one
thousand simulations with different, random initial conditions. Further
simulations across different phase-lag values demonstrate that the identified
pattern of eigenmode contributions holds across the range of parameters that
lead to transient chimera states in these networks, and for conditions where
the synchronized cluster appears at different positions in the network (Figs.
2b, 2c, 2d, and 2e). It is important to note that our approach generalizes to
a diversity of networks, where the eigenspectrum of matrices that are not
circulant can be obtained numerically (see Fig. S2 in the Supplement). Taken
together, these results provide further evidence that a specific combination
of modes occurs generally, across initial conditions and networks, when
transient chimera states appear.

What network mechanisms lead to this specific combination of modes, and thus
to chimera states? To answer this question, we look to the eigenvalues of
$\bm{K}$ to understand how these specific combinations of modes arise. Figure
3a shows the eigenvalues of $\bm{K}=\epsilon e^{-\mathrm{i}\phi}\bm{A}$ for
increasing values of $\phi$. Our analytical approach allows us to understand
the effect of $\phi$ as a rotation of these eigenvalues in the complex plane.
When $\phi=0$, the eigenvalues of $\bm{K}$ fall along the real axis (blue
dots, Fig. 3a). When $\phi=\nicefrac{{\pi}}{{2}}$, the eigenvalues fall along
the imaginary axis (red dots, Fig. 3a). This rotation reduces the real part of
the $\lambda_{i}$ that were initially positive, reduces the real part of
$\lambda_{1}$ relative to the other eigenvalues, and - for the range of $\phi$
where chimera states are produced - the rotation leaves a small set of modes
making a significant contribution to the dynamics (Fig. 3b). For all cases
with $\phi<\nicefrac{{\pi}}{{2}}$, however, the leading eigenvalue (in real
part), is given by $\lambda_{1}$ which is associated with the phase
synchronizing mode. Our approach thus analytically explains previous numerical
observations that chimeras are transient states that are followed by phase
synchronization in finite networks of Sakaguchi-Kuramoto oscillators [37].

Figure 3: Phase-lag and the emergence of transient chimera states. (a) As the
phase-lag increases, the eigenvalues of $\bm{K}$ rotate with an angle $-\phi$
in the complex plane going from being purely real (for $\phi=0$) to purely
imaginary (for $\phi=\nicefrac{{\pi}}{{2}}$). (b) This rotation decreases the
difference in the real part of the eigenvalues, which allows the network to
access the modes relevant for the appearance of transient chimera states. (c)
We use a chimera index $\overline{\sigma}$ [45] for the numerical simulations
(gray bars), where we observe that chimeras appear for a specific range of
phase-lag $\phi$. These results are given by an average over $100$ simulations
of $t=10^{6}$ timesteps with errorbars indicating the standard deviation. We
then use the spectral properties of $\bm{K}$ through Eq. (3) to analytically
study this (black line). (d) This approach generalizes for networks with
different connectivity (given by different values of $\alpha$), where our
theory (black line) is able to capture the value of phase-lag where the
chimera index is maximum (gray dots). (e), Lastly, we show examples of the
synchronization level (blue line) and the chimera index (yellow line) as a
function of time for different phase-lag values based on panel (c) – each
represented by a different symbol.

Our theory offers further analytical insights into the role of phase-lag and
connectivity on the emergence of transient chimeras. First, we use a chimera
index $\overline{\sigma}$ [45] introduced in [46] to measure the existence of
chimera states in numerical simulations (see Section III in the Supplement).
We consider a distance-dependent network with $\alpha=1.0$ and vary the phase-
lag $\phi\in[0,\nicefrac{{\pi}}{{2}}]$. Consistent with previous observations
[46], chimera states appear for a specific range of phase-lag (gray bars, Fig.
3c). We then consider an analytical quantifier, based only on the spectral
properties of $\bm{K}$, to study this phenomenon:

${\color[rgb]{0,0,0}\zeta(\alpha,\phi)={\Big{|}\lambda_{1}^{\ast}-\frac{1}{{\Delta_{1,2}}^{\ast}}\Big{|}}^{-1},}$
(3)

where
$\lambda_{1}^{\ast}=\eta(\alpha)\frac{\mathrm{Real}[\lambda_{1}(\alpha,\phi)]}{\mathrm{Real}[\lambda_{1}(\alpha,0)]}$
and
${\Delta_{1,2}}^{\ast}=\frac{\mathrm{Real}[\lambda_{1}(\alpha,\phi)-\lambda_{2}(\alpha,\phi)]}{\mathrm{Real}[\lambda_{1}(\alpha,0)-\lambda_{2}(\alpha,0)]}$
(see Supplemental Sec. IV), with $\eta(\alpha)$ being the normalization factor
for the distance-dependent network [42]. This quantifier captures the balance
between the first mode and the low spatial frequency modes, and thus expresses
the emergence of chimera states. By using only two eigenvalues of $\bm{K}$, we
are able to capture the non-monotonic behavior as a function of $\phi$ (black
line, Fig. 3c). This behavior has been numerically observed in [46], and our
theory now offers an analytical explanation for it. This result generalizes
for networks with different connectivity (different values of $\alpha$), where
we can use our theory to analytically predict the precise phase-lag value
where the chimera index is maximum (Fig. 3d). Lastly, to illustrate the
complexity of this phenomenon, we plot the order parameter and the chimera
index as a function of time (blue lines and yellow lines, Fig. 3e,
respectively) for different values of $\phi$ based on Fig. 3c (represented by
different symbols). We observe that if the phase-lag is too small, the network
quickly reaches phase synchronization (pink hexagon, Fig. 3e). As the phase-
lag increases, the transient becomes longer and chimera states appear (green
square and brown star, Fig. 3e). If the phase-lag is too close to
$\nicefrac{{\pi}}{{2}}$, however, the network remains in an asynchronous state
and no chimera state is observed (purple circle, Fig. 3e), and the order
parameter and chimera index fluctuate due to the finite-size effects.

Figure 4: Designing networks to exhibit chimeras. (a) We start with a global
network, which adjacency matrix $\bm{A}$ is given by the complete graph. (b)
We use the CDT to obtain the eigenvalues of the matrix $\bm{K}=\epsilon
e^{-\mathrm{i}\phi}\bm{A}$ for this network (blue dots). We then manually
modify these eigenvalues in order to obtain a similar fashion as revealed by
our approach in Fig. 3 with the modes relevant for chimera states having a
higher real part in $\lambda$ (insets b). (c) We build the rewired network,
where the adjacency matrix now has weighted connections. We then numerically
integrate Eq. (1) considering: (d) the global network, where the system
transitions to phase synchronization; (e) the rewired network with the same
initial conditions, where the system exhibits chimera states. Here,
$\epsilon=0.08$ and $\phi=0$.

We now show that we can “rewire” a network to exhibit the eigenvalue
configuration in Fig. 3 and then express chimera states. To do this, we start
with a globally connected network (the complete graph on $N$ nodes, $A_{ij}=1$
for all $i\in[1,N],i\neq j$; Fig. 4a) and with no phase-lag ($\phi=0$). We
obtain the eigenvalues and eigenvectors of $\bm{K}=\epsilon
e^{-\mathrm{i}\phi}\bm{A}$ using the CDT, as above (blue dots, Fig. 4b). In
the complete graph, only $\lambda_{1}$ has positive real part (Fig. 4b). Our
analytical approach predicts this will lead the network to phase synchrony
because $\lambda_{1}$ is associated with the zero phase-difference mode
(inset, Fig. 4b). We then manually adjusted the eigenvalues to exhibit a
configuration that our approach reveals to give rise to chimera states, as
shown in Fig. 3, and then obtained a “rewired” matrix
$\bm{A}^{\prime}=\bm{V}\bm{D}^{\prime}\bm{V}^{-1}$ (Fig. 4c) with the new
eigenvalues in $\bm{D}^{\prime}$. In this case, while a numerical simulation
of the original Kuramoto network connected by matrix $\bm{A}$ (complete graph)
transitions to synchrony (Fig. 4d), the network with the new matrix
$\bm{A}^{\prime}$ now generates a chimera (Fig. 4e). This result demonstrates
a specific analytical procedure for designing networks to produce chimera
states in Kuramoto networks. This topic has received increasing attention in
the past years with development of new algorithms to create specific networks
to produce a desired pattern [47, 48, 49], and our analytical approach could
have important impact for future machine learning algorithms for training
neural networks to produce chimeras and other patterns.

There has been much recent interest in the possibility that neural systems
exhibit chimera states, from the dynamics of spiking neural networks [9], to
the oscillations associated with memory storage and retrieval [50], and to the
rhythms of sleep [12]. As a specific example, the 0.1-1 Hz slow oscillation
that occurs in the deepest stages of non-rapid-eye-movement (non-REM) sleep
has been observed to occur only in one hemisphere of the brain at a time,
while the other half of the brain remains de-synchronized [51, 52, 53, 54].
Explaining how one half of the brain can remain desynchronized while the other
half exhibits highly synchronized neuronal activity has suggested connections
to chimera states in oscillator networks [49]. In uni-hemispheric sleep,
synchronization occurs across networks of 50-100 cortical regions. While many
mathematical approaches have provided insight into the origin of chimera
states and provided bifurcation analyses [2, 26], our mathematical approach
provides specific insight into the spatiotemporal impact of the eigenvectors
of the network adjacency matrix in shaping these dynamics. This analytical
insight links the eigenspectrum of the graph adjacency matrix to the resulting
nonlinear dynamics in Sakaguchi-Kuramoto networks. This is, in general, a very
difficult problem in applied mathematics [55], but our analysis allows us to
make this connection in this case. A key feature of our mathematical approach
is to analyze the equations defining the network dynamics at the finite scale,
which in turn allows us to understand the emergence of transient chimeras at
the scale of tens to hundreds of nodes that is relevant for understanding
neural dynamics during sleep.

Our results demonstrate how this operator-based approach to networked systems
at the finite scale can provide deep insight into the interplay of
connectivity and nonlinearity at individual network nodes. Here, we have
provided analytical insight into the interaction of local connectivity -
including distance-dependent power-law graphs, ring graphs, clustered
networks, and Watts-Strogatz graphs - and a phase-lagged sine interaction
function to produce transient chimera states. Our theory provides a simple and
precise mathematical formulation of how eigenvectors of the adjacency matrix
combine to produce transient chimera states in finite Sakaguchi-Kuramoto
dynamics (Figs. 1 and 2), the interplay between phase-lag and network
connectivity that leads to the emergence of these states (Fig. 3), and even
how to design patterns of connectivity to produce chimera states (Fig. 4).
Previous work has observed connections between the eigenvalues of the
Laplacian matrix and phase synchronization in oscillator networks [56, 57,
58]. Our theoretical approach provides a precise connection between network
spectra and dynamics in two ways. First, our approach focuses on the network
adjacency matrix, which is more directly related to the network dynamics,
instead of the Laplacian. Second, our approach connects both the eigenvectors
of the adjacency matrix and their corresponding eigenvalues to the resulting
dynamics (specifically, when analyzing the matrix $\bm{K}$).

Lastly, this work opens new avenues in the study of sophisticated
spatiotemporal dynamics in oscillator networks. Future work can leverage our
analytical predictions to control chimera states, a topic which has been
studied with interest via numerical simulations until now [59, 60], or to
study and control transient dynamics, which has also received much interest
[61, 62], because of its empirical relevance [63, 64, 56]. Taking all these
points together, it is possible to utilize chimera states and the analytical
framework introduced here, to design controllable and interpretable dynamics
that can advance our understanding of computation with spatiotemporal dynamics
in nonlinear systems [65].

###### Acknowledgements.

This work was supported by BrainsCAN at Western University through the Canada
First Research Excellence Fund (CFREF), the NSF through a NeuroNex award
(#2015276), the Natural Sciences and Engineering Research Council of Canada
(NSERC) grant R0370A01, Compute Ontario (computeontario.ca), Digital Research
Alliance of Canada (alliancecan.ca), and the Western Academy for Advanced
Research. R.C.B gratefully acknowledges the Western Institute for Neuroscience
Clinical Research Postdoctoral Fellowship.

## References

  * Kuramoto and Battogtokh [2002] Y. Kuramoto and D. Battogtokh, arXiv preprint cond-mat/0210694 (2002).
  * Abrams and Strogatz [2004] D. M. Abrams and S. H. Strogatz, Physical Review Letters 93, 174102 (2004).
  * Montbrió  _et al._ [2004] E. Montbrió, J. Kurths, and B. Blasius, Physical Review E 70, 056125 (2004).
  * Zhu  _et al._ [2014] Y. Zhu, Z. Zheng, and J. Yang, Physical Review E 89, 022914 (2014).
  * Abrams  _et al._ [2008] D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A. Wiley, Physical Review Letters 101, 084103 (2008).
  * Laing  _et al._ [2012] C. R. Laing, K. Rajendran, and I. G. Kevrekidis, Chaos: An Interdisciplinary Journal of Nonlinear Science 22, 013132 (2012).
  * Bansal  _et al._ [2019] K. Bansal, J. O. Garcia, S. H. Tompson, T. Verstynen, J. M. Vettel, and S. F. Muldoon, Science Advances 5, eaau8535 (2019).
  * Boaretto  _et al._ [2021] B. R. R. Boaretto, R. C. Budzinski, K. L. Rossi, C. Manchein, T. L. Prado, U. Feudel, and S. R. Lopes, Physical Review E 104, 024204 (2021).
  * Majhi  _et al._ [2019] S. Majhi, B. K. Bera, D. Ghosh, and M. Perc, Physics of Life Reviews 28, 100 (2019).
  * Santos  _et al._ [2017] M. S. Santos, J. D. Szezech, F. S. Borges, K. C. Iarosz, I. L. Caldas, A. M. Batista, R. L. Viana, and J. Kurths, Chaos, Solitons & Fractals 101, 86 (2017).
  * Banerjee  _et al._ [2016] T. Banerjee, P. S. Dutta, A. Zakharova, and E. Schöll, Physical Review E 94, 032206 (2016).
  * Parastesh  _et al._ [2021] F. Parastesh, S. Jafari, H. Azarnoush, Z. Shahriari, Z. Wang, S. Boccaletti, and M. Perc, Physics Reports 898, 1 (2021).
  * Dudkowski  _et al._ [2014] D. Dudkowski, Y. Maistrenko, and T. Kapitaniak, Physical Review E 90, 032920 (2014).
  * Viktorov  _et al._ [2014] E. A. Viktorov, T. Habruseva, S. P. Hegarty, G. Huyet, and B. Kelleher, Physical Review Letters 112, 224101 (2014).
  * Hagerstrom  _et al._ [2012] A. M. Hagerstrom, T. E. Murphy, R. Roy, P. Hövel, I. Omelchenko, and E. Schöll, Nature Physics 8, 658 (2012).
  * Totz  _et al._ [2018] J. F. Totz, J. Rode, M. R. Tinsley, K. Showalter, and H. Engel, Nature Physics 14, 282 (2018).
  * Gambuzza  _et al._ [2020] L. V. Gambuzza, L. Minati, and M. Frasca, Chaos, Solitons & Fractals 138, 109907 (2020).
  * Patzauer and Krischer [2021] M. Patzauer and K. Krischer, Physical Review Letters 126, 194101 (2021).
  * Hart  _et al._ [2016] J. D. Hart, K. Bansal, T. E. Murphy, and R. Roy, Chaos: an interdisciplinary journal of nonlinear science 26, 094801 (2016).
  * Tinsley  _et al._ [2012] M. R. Tinsley, S. Nkomo, and K. Showalter, Nature Physics 8, 662 (2012).
  * Wickramasinghe and Kiss [2013] M. Wickramasinghe and I. Z. Kiss, PloS One 8, e80586 (2013).
  * Omel’chenko  _et al._ [2008] E. Omel’chenko, Y. L. Maistrenko, and P. A. Tass, Physical Review Letters 100, 044105 (2008).
  * Kotwal  _et al._ [2017] T. Kotwal, X. Jiang, and D. M. Abrams, Physical Review Letters 119, 264101 (2017).
  * Nicolaou  _et al._ [2017] Z. G. Nicolaou, H. Riecke, and A. E. Motter, Physical Review Letters 119, 244101 (2017).
  * Laing [2009] C. R. Laing, Chaos: An Interdisciplinary Journal of Nonlinear Science 19, 013113 (2009).
  * Panaggio and Abrams [2015] M. J. Panaggio and D. M. Abrams, Nonlinearity 28, R67 (2015).
  * Bogomolov  _et al._ [2017] S. A. Bogomolov, A. V. Slepnev, G. I. Strelkova, E. Schöll, and V. S. Anishchenko, Communications in Nonlinear Science and Numerical Simulation 43, 25 (2017).
  * Schülen  _et al._ [2022] L. Schülen, A. Gerdes, M. Wolfrum, and A. Zakharova, Physical Review E 106, L042203 (2022).
  * Ashwin and Burylko [2015] P. Ashwin and O. Burylko, Chaos: An Interdisciplinary Journal of Nonlinear Science 25 (2015).
  * Bick and Ashwin [2016] C. Bick and P. Ashwin, Nonlinearity 29, 1468 (2016).
  * Bick  _et al._ [2017] C. Bick, M. Sebek, and I. Z. Kiss, Physical Review Letters 119, 168301 (2017).
  * Maistrenko  _et al._ [2017] Y. Maistrenko, S. Brezetsky, P. Jaros, R. Levchenko, and T. Kapitaniak, Physical Review E 95, 010203 (2017).
  * Burylko  _et al._ [2022] O. Burylko, E. A. Martens, and C. Bick, Chaos: An Interdisciplinary Journal of Nonlinear Science 32 (2022).
  * Zhang  _et al._ [2020] Y. Zhang, Z. G. Nicolaou, J. D. Hart, R. Roy, and A. E. Motter, Physical Review X 10, 011044 (2020).
  * Zhang and Motter [2021] Y. Zhang and A. E. Motter, Physical Review Letters 126, 094101 (2021).
  * Palmigiano  _et al._ [2017] A. Palmigiano, T. Geisel, F. Wolf, and D. Battaglia, Nature Neuroscience 20, 1014 (2017).
  * Wolfrum and Omel’chenko [2011] M. Wolfrum and E. Omel’chenko, Physical Review E 84, 015201 (2011).
  * Papo and Buldú [2019] D. Papo and J. M. Buldú, NeuroImage 196, 195 (2019).
  * Muller  _et al._ [2021] L. Muller, J. Mináč, and T. T. Nguyen, Physical Review E 104, L022201 (2021).
  * Budzinski  _et al._ [2022] R. C. Budzinski, T. T. Nguyen, J. Đoàn, J. Mináč, T. J. Sejnowski, and L. E. Muller, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 031104 (2022).
  * Budzinski  _et al._ [2023] R. C. Budzinski, T. T. Nguyen, G. B. Benigno, J. Đoàn, J. Mináč, T. J. Sejnowski, and L. E. Muller, Physical Review Research 5, 013159 (2023).
  * [42] The distance-dependent adjacency matrix is given by: $A_{ij}=\nicefrac{{1}}{{(\eta(d_{ij})^{\alpha})}}$, where $d_{ij}=\mathrm{min}(|i-j|,N-|i-j|)$ and $\eta=\sum_{j=1,j\neq i}^{N}(d_{ij})^{-\alpha}$. See Supplement Material (Sec. I) for more details.
  * Nguyen  _et al._ [2023] T. T. Nguyen, R. C. Budzinski, J. Ðoàn, F. W. Pasini, J. Mináč, and L. E. Muller, SIAM Journal on Applied Dynamical Systems 22, 802 (2023).
  * Davis [1979] P. Davis,  _Circulant matrices_ (John Wiley & Sons, 1979).
  * [45] The chimera index is based on the phase synchronization level of different parts of the network: we evaluate $\sigma(t)=\sqrt{\frac{\sum_{k=1}^{M}{R_{k}(t)}-\overline{R}(t)}{M}},$where $M$ is the number of groups, $R_{k}(t)$ is the order parameter of each group, and $\overline{R}(t)$ is the average order parameter over all groups. We then take the temporal average of $\sigma(t)$, represented by $\overline{\sigma}$. See Supplement Material (Sec. III) for more details.
  * Shanahan [2010] M. Shanahan, Chaos: An Interdisciplinary Journal of Nonlinear Science 20, 013108 (2010).
  * Forrow  _et al._ [2018] A. Forrow, F. G. Woodhouse, and J. Dunkel, Physical Review X 8, 041043 (2018).
  * Menara  _et al._ [2022] T. Menara, G. Baggio, D. Bassett, and F. Pasqualetti, Nature Communications 13, 4721 (2022).
  * Masoliver  _et al._ [2022] M. Masoliver, J. Davidsen, and W. Nicola, Communications Physics 5, 205 (2022).
  * Masoliver  _et al._ [2023] M. Masoliver, J. Davidsen, and W. Nicola, bioRxiv , 2023 (2023).
  * Rattenborg  _et al._ [1999] N. C. Rattenborg, S. L. Lima, and C. J. Amlaner, Nature 397, 397 (1999).
  * Mascetti [2016] G. G. Mascetti, Nature and Science of Sleep , 221 (2016).
  * Tamaki  _et al._ [2016] M. Tamaki, J. W. Bang, T. Watanabe, and Y. Sasaki, Current Biology 26, 1190 (2016).
  * Rattenborg  _et al._ [2016] N. C. Rattenborg, B. Voirin, S. M. Cruz, R. Tisdale, G. Dell’Omo, H.-P. Lipp, M. Wikelski, and A. L. Vyssotski, Nature Communications 7, 12468 (2016).
  * Curto and Morrison [2019] C. Curto and K. Morrison, Current Opinion in Neurobiology 58, 11 (2019).
  * Ocampo-Espindola  _et al._ [2021] J. L. Ocampo-Espindola, E. Omel’chenko, and I. Z. Kiss, Journal of Physics: Complexity 2, 015010 (2021).
  * Skardal  _et al._ [2014] P. S. Skardal, D. Taylor, and J. Sun, Physical Review Letters 113, 144101 (2014).
  * McGraw and Menzinger [2007] P. N. McGraw and M. Menzinger, Physical Review E 75, 027104 (2007).
  * Bick and Martens [2015] C. Bick and E. A. Martens, New Journal of Physics 17, 033030 (2015).
  * Sieber  _et al._ [2014] J. Sieber, E. Omel’chenko, and M. Wolfrum, Physical Review Letters 112, 054102 (2014).
  * Lilienkamp and Parlitz [2020] T. Lilienkamp and U. Parlitz, Physical Review E 102, 032219 (2020).
  * Omel’chenko and Tél [2022] E. Omel’chenko and T. Tél, Journal of Physics: Complexity 3, 010201 (2022).
  * Rabinovich  _et al._ [2008] M. Rabinovich, R. Huerta, and G. Laurent, Science 321, 48 (2008).
  * Brinkman  _et al._ [2022] B. A. W. Brinkman, H. Yan, A. Maffei, I. M. Park, A. Fontanini, J. Wang, and G. La Camera, Applied Physics Reviews 9, 011313 (2022).
  * Budzinski  _et al._ [tion] R. C. Budzinski, A. N. Busch, S. Mestern, E. Martin, L. H. B. Liboni, F. W. Pasini, J. Minac, T. Coleman, W. Inoue, and L. Muller, (in preparation).


# Asymptotic average solutions
to linear second order semi-elliptic PDEs:
a Pizzetti-type Theorem

Alessia E. Kogoj  Dipartimento di Scienze Pure e Applicate (DiSPeA)
Università degli Studi di Urbino Carlo Bo
Piazza della Repubblica 13, 61029 Urbino (PU), Italy<EMAIL_ADDRESS>and  Ermanno Lanconelli  Dipartimento di Matematica
Alma Mater Studiorum Università di Bologna
Piazza di Porta San Donato 5, 40126 Bologna, Italy.
<EMAIL_ADDRESS>
###### Abstract.

By exploiting an old idea first used by Pizzetti for the classical Laplacian,
we introduce a notion of asymptotic average solutions making pointwise
solvable every Poisson equation ${\mathcal{L}}u(x)=-f(x)$ with continuous data
$f$, where ${\mathcal{L}}$ is a hypoelliptic linear partial differential
operator with positive semidefinite characteristic form.

###### Key words and phrases:

asymptotic mean value formulas, semi-elliptic operators, hypoelliptic
operators, Poisson type-equations

###### 2010 Mathematics Subject Classification:

35D99, 35B05, 35J70, 35H10

##  1\. Introduction

The Poisson-type equations related to hypoelliptic linear second order PDE’s
with nonnegative characteristic form cannot be studied in $L^{p}$ spaces due
to the lack of a suitable Calderon-Zygmund theory for the relevant singular
integrals. Our paper presents a result allowing to satisfactory study such
equations in spaces of continuous functions. We follow a procedure introduced
by Pizzetti in his 1909’s paper [14] based on the asymptotic average solutions
for the classical Poisson-Laplace equation.

### 1.1.

Let $\Omega$ be a bounded open subset of ${{\mathbb{R}}^{{n}}}$, and let
$f:\Omega\longrightarrow{\mathbb{R}}$ be a continuous bounded function. Let us
denote by $u_{f}$ the Newtonian potential of $f$, i.e.,

$u_{f}:{{\mathbb{R}}^{{n}}}\longrightarrow{\mathbb{R}},\qquad
u_{f}(x):=\int_{\Omega}\Gamma(y-x)f(y)\ dy.$

Here $\Gamma$ denotes the fundamental solution of the Laplace equation, i.e.,

$\Gamma(x)=c_{n}|x|^{2-n},x\in{{\mathbb{R}}^{{n}}}\smallsetminus\\{0\\},$

$\omega_{n}$ being the volume of the unit ball in ${{\mathbb{R}}^{{n}}}$ and
$c_{n}:=\dfrac{1}{n(n-2)\omega_{n}}$.

It is well known that $u_{f}\in C^{1}({{\mathbb{R}}^{{n}}},{\mathbb{R}}),$
while, in general, ${{u_{f}}|}_{\Omega}\notin C^{2}(\Omega,{\mathbb{R}})$.
However, in the weak sense of distributions,

(1.1) $\varDelta u_{f}=-f\mbox{ in }\Omega.$

As a consequence, if the continuous function $f$ is such that

(1.2) $u_{f}\notin C^{2}(\Omega,{\mathbb{R}}),$

then the Poisson equation

(1.3) $\varDelta v=-f$

has no classical solutions, i.e., there does not exist a function $v\in
C^{2}(\Omega,{\mathbb{R}})$ satisfying

$\varDelta v(x)=-f(x)\mbox{ for every }x\in\Omega.$

Indeed, assume by contradiction that such a function exists. Then, by (1.1),

$\varDelta(u_{f}-v)=0\mbox{ in }\Omega$

in the weak sense of distributions, so that, by Caccioppoli–Weyl’s Lemma,
there exists a function $h$, harmonic in $\Omega$, such that

$u_{f}(x)-v(x)=h(x)$

a.e. in $\Omega$. Therefore, $u_{f}-v$ being continuous in $\Omega$,

$u_{f}=v+h\in C^{2}(\Omega,{\mathbb{R}}),$

in contradiction with (1.2). This proves the existence of continuous functions
$f$ such that the Poisson equation (1.3) is not pointwise solvable. In his
paper [14], Pizzetti introduced a notion of pointwise weak Laplacian, making
pointwise solvable every Poisson equation with continuous data. Pizzetti
started from the following remark. Given a function $u$ of class $C^{2}$ in
$\Omega$ one has

(1.4) $\lim_{r\longrightarrow
0}\frac{M_{r}(u)(x)-u(x)}{r^{2}}=\frac{1}{2(n+2)}\varDelta u(x)$

for every $x\in\Omega.$ Here $M_{r}$ denotes the Gauss average

$M_{r}(u)(x):=\frac{1}{|B(x,r)|}\int_{\partial B(x,r)}u(y)\ dy,$

$|B(x,r)|$ being the volume of $B(x,r)$, the Euclidean ball centered at $x$
with radius $r$. Then, if $u\in C(\Omega,{\mathbb{R}})$ is such that the limit
at the left hand side of (1.4) exists at a point $x\in\Omega$, Pizzetti
defines

$\varDelta_{a}u(x):=2(n+2)\lim_{r\longrightarrow
0}\frac{M_{r}(u)(x)-u(x)}{r^{2}}.$

We call $\varDelta_{a}u(x)$ the asymptotic average Laplacian of $u$ at $x$.
Keeping in mind (1.4), if $u\in C^{2}(\Omega,{\mathbb{R}})$, then

$\varDelta_{a}u(x)=\varDelta u(x)\mbox{ for every }x\in\Omega.$

We denote by

$\mathcal{A}(\Omega,\varDelta)$

the class of functions $u\in C(\Omega,{\mathbb{R}})$, such that
$\varDelta_{a}u(x)$ exists at any point $x\in\Omega.$ Obviously,
$\mathcal{A}(\Omega,\varDelta)$ is a (linear) sub-space of
$C(\Omega,{\mathbb{R}}).$ Moreover, by the previous remark,

$C^{2}(\Omega,{\mathbb{R}})\subseteq\mathcal{A}(\Omega,\varDelta).$

Pizzetti proved that the Newtonian potentials of continuous bounded functions
are contained in $\mathcal{A}(\Omega,\varDelta).$ Precisely he proved the
following theorem.

######  Theorem A (Pizzetti Theorem).

Let $\Omega\subseteq{{\mathbb{R}}^{{n}}},n\geq 3,$ be a bounded open subset of
${{\mathbb{R}}^{{n}}}$ and let $f:\Omega\longrightarrow{\mathbb{R}}$ be a
bounded continuous function. Then

$u_{f}\in\mathcal{A}(\Omega,\varDelta)$

and

$\varDelta_{a}u_{f}=-f\mbox{ in }\Omega.$

The aim of this paper is to extend the notion of asymptotic average solution
and Pizzetti’s Theorem to the class of linear second order semi-elliptic
partial differential operators that we will introduce in the next subsection.

### 1.2.

We will deal with partial differential operators of the type

(1.5)
${\mathcal{L}}=\sum_{i,j=1}^{n}\partial_{x_{i}}(\partial_{x_{j}}a_{ij}(x)),\
x\in{{\mathbb{R}}^{{n}}},$

where $A(x):=(a_{ij}=a_{ji})_{i,j=1,\ldots,n}$ is a symmetric nonnegative
definite matrix,

$x\longmapsto a_{ji}(x),\qquad i,j=1,\ldots,n$

are smooth functions in ${{\mathbb{R}}^{{n}}}$ and

$\sum_{i=1}^{n}a_{ii}(x)>0\mbox{ for every }x\in{{\mathbb{R}}^{{n}}}.$

Together with these qualitative properties we assume that ${\mathcal{L}}$ is
hypoelliptic in ${{\mathbb{R}}^{{n}}}$ and endowed with a smooth fundamental
solution

$\Gamma:\\{(x,y)\in{{\mathbb{R}}^{{n}}}\times{{\mathbb{R}}^{{n}}}\ |\ x\neq
y\\}\longrightarrow{\mathbb{R}},$

such that

  * (i)

$\Gamma(x,y)=\Gamma(y,x)>0,\mbox{ for every }x\neq y;$

  * (ii)

$\lim_{x\rightarrow y}\Gamma(x,y)=\infty,\mbox{ for every
}y\in{{\mathbb{R}}^{{n}}};$

  * (iii)

$\lim_{x\rightarrow\infty}\left(\sup_{y\in K}\Gamma(x,y)\right)=0$, for every
compact set $K\subseteq{{\mathbb{R}}^{{n}}};$

  * (iv)

$\Gamma(x,\cdot)$ belongs to $L^{1}_{\rm loc}({{\mathbb{R}}^{{n}}})$, for
every $x\in{{\mathbb{R}}^{{n}}}$.

We recall that when we say that $\Gamma$ is a fundamental solution of
${\mathcal{L}}$ we mean that, for every $\varphi\in
C_{0}^{\infty}({{\mathbb{R}}^{{n}}},{\mathbb{R}})$ and
$x\in{{\mathbb{R}}^{{n}}}$:

$\int_{{\mathbb{R}}^{{n}}}\Gamma(x,y)\,{\mathcal{L}}\varphi(y)\
dy=-\varphi(x).$

### 1.3.

Important examples of operators satisfying our assumptions are the “sum of
squares”of homogeneous Hörmander vector fields. Precisely: let

$X=\\{X_{1},\ldots,X_{m}\\}$

be a family of linearly independent smooth vector fields such that

  * (H1)

$X_{1},\ldots,X_{m}$ satisfy the Hörmander rank condition at $x=0$, that is,

$\dim\\{Y(0)\ |\ Y\in\mathrm{Lie}\\{X_{1},\ldots,X_{m}\\}\\}=n;$

  * (H2)

$X_{1},\ldots,X_{m}$ are homogeneous of degree $1$ with respect to a group of
dilations $(\delta_{\lambda})_{\lambda>0}$ of the following type

$\delta_{\lambda}:{{\mathbb{R}}^{{n}}}\longrightarrow{{\mathbb{R}}^{{n}}},$
$\delta_{\lambda}(x)=\delta_{\lambda}(x_{1},\ldots,x_{n})=(\lambda^{\sigma}x_{1},\ldots,\lambda^{\sigma_{n}}x_{n}),$

where the $\sigma_{j}$’s are natural numbers such that
$1\leq\sigma_{1}\leq\ldots\leq\sigma_{n}.$

Then,

(1.6) ${\mathcal{L}}=\sum_{j=1}^{m}X_{j}^{2}$

satisfies all the assumptions listed in subsection 1.2 (see [1], [2]).

We stress that the sub-Laplacians on stratified Lie groups in
${{\mathbb{R}}^{{n}}}$ are particular cases of the operator ${\mathcal{L}}$ in
(1.6).

### 1.4.

The extension of Pizzetti’s Theorem to the operator ${\mathcal{L}}$ in (1.5)
rests on some representation formulas on the superlevel set of $\Gamma$. If
$x\in{\mathbb{R}}$ and $r>0$, define

$\Omega_{r}(x):=\left\\{y\in{{\mathbb{R}}^{{n}}}\ :\
\Gamma(x,y)>\frac{1}{r}\right\\}.$

We will call $\Omega_{r}(x)$ the ${\mathcal{L}}$-ball centered at $x$ and with
radius $r$. It is easy to recognize that $\Omega_{r}(x)$ is a nonempty bounded
open set of ${{\mathbb{R}}^{{n}}}$. Moreover

(1.7) $\bigcap_{r>0}\Omega_{r}(x)=\\{x\\}$

and

$\frac{|\Omega_{r}(x)|}{r}\longrightarrow 0\mbox{ as }r\longrightarrow 0.$

######  Remark 1.1.

If ${\mathcal{L}}=\varDelta$, then

$\Omega_{r}(x)=B(x,\rho),\mbox{ \ with }\rho=(c_{n}r)^{\frac{1}{n-2}}.$

Let $\Omega\subseteq{{\mathbb{R}}^{{n}}}$ be open and let $u\in
C^{2}(\Omega,{\mathbb{R}})$. Then, for every ${\mathcal{L}}$-ball,
$\Omega_{r}(x)$ such that $\overline{\Omega_{r}(x)}\subseteq\Omega$ and for
every $\alpha>-1$ we have

(1.8) $u(x)=M_{r}(u)(x)-N_{r}({\mathcal{L}}u)(x),$

where $M_{r}$ and $N_{r}$ are the following average operators:

(1.9)
$M_{r}(u)(x):=\frac{\alpha+1}{r^{\alpha+1}}\int_{\Omega_{r}(x)}u(y)K(x,y)\
dy,$

where

$K(x,y):=\frac{\langle
A(y)\nabla_{y}\Gamma(x,y),\nabla_{y}\Gamma(x,y)\rangle}{(\Gamma(x,y))^{\alpha+2}};$

(1.10)
$N_{r}(w)(x):=\frac{\alpha+1}{r^{\alpha+1}}\int_{0}^{r}\rho^{\alpha}\left(\int_{\Omega_{r}(x)}\left(\Gamma(x,y)-\frac{1}{\rho}\right)w(y)\
dy\right)\ d\rho.$

The proof of the representation formula (1.8) can be found in [4].

######  Remark 1.2.

If ${\mathcal{L}}=\varDelta$ and $\alpha=\dfrac{2}{n-2}$, then the kernel $K$
is constant and $M_{r}$ becomes the Gauss average on the Euclidean ball
$B(x,\rho),$ with $\rho=(c_{n}r)^{\frac{1}{n-2}}.$

Letting

(1.11)
$Q_{r}(x):=N_{r}(1)=\frac{\alpha+1}{r^{\alpha+1}}\int_{0}^{r}\rho^{\alpha}\left(\int_{\Omega_{r}(x)}(\Gamma(x,y)-\frac{1}{\rho}\right)\
dy)\ d\rho,$

an easy computation shows that

$Q_{r}(x)=\int_{0}^{r}\frac{\Omega_{\rho}(x)}{\rho^{2}}\left(1-\left(\frac{\rho}{r}\right)^{\alpha+1}\right)\
d\rho.$

######  Remark 1.3.

If ${\mathcal{L}}=\varDelta$ and $\alpha=\dfrac{2}{n-2}$, then, letting
$\rho=(c_{n}r)^{\frac{1}{n-2}}$, we get

$\frac{M_{r}(u)(x)-u(x)}{Q_{r}(x)}=2(n+2)\frac{\frac{1}{|B(x,\rho)|}\int_{B(x,\rho)}u(y)\
dy-u(0)}{\rho^{2}},$

so that, by (1.4),

(1.12) $\lim_{r\longrightarrow 0}\frac{M_{r}(u)(x)-u(x)}{Q_{r}(x)}=\varDelta
u(x).$

The limit in (1.12) extends to all the operators ${\mathcal{L}}$ in (1.5).
Indeed, if $u$ is a $C^{2}$ function in an open set
$\Omega\subseteq{{\mathbb{R}}^{{n}}}$, from the representation formula (1.8)
and the identity (1.7), using Corollary 2.5 in Section 2, one immediately gets

$\lim_{r\longrightarrow
0}\frac{M_{r}(u)(x)-u(x)}{Q_{r}(x)}={\mathcal{L}}u(x).$

Then, in analogy with the case ${\mathcal{L}}=\varDelta$, we introduce the
following definition.

######  Definition 1.4.

Let ${\mathcal{L}}$ be a partial differential operator satisfying the
assumptions of subsection 1.2 and let $u$ be a continuous function in an open
set $\Omega\subseteq{{\mathbb{R}}^{{n}}}$. We say that

$u\in\mathcal{A}(\Omega,{\mathcal{L}}),$

if

$\lim_{r\longrightarrow 0}\frac{M_{r}(u)(x)-u(x)}{Q_{r}(x)}$

exists in ${\mathbb{R}}$ at every point $x\in\Omega.$ In this case we define

$({\mathcal{L}}_{a}(u))(x):=\lim_{r\longrightarrow
0}\frac{M_{r}(u)(x)-u(x)}{Q_{r}(x)}.$

Furthermore, if $f\in C(\Omega,{\mathbb{R}})$ and there exists
$u\in\mathcal{A}(\Omega,{\mathcal{L}})$ such that

$({\mathcal{L}}_{a}u)(x)=f(x)\mbox{ for every }x\in\Omega,$

we say that $u$ is an asymptotic average solution to

${\mathcal{L}}_{a}u=f\mbox{ in }\Omega.$

In the case $f=0$ this definition was first introduced in the paper [6].

The main result of our paper is the following theorem which extends Pizzetti’s
Theorem to the operators (1.5).

######  Theorem 1.5.

Let $f:{{\mathbb{R}}^{{n}}}\longrightarrow{\mathbb{R}}$ be a compactly
supported continuous function. Define

$u_{f}(x):=\int_{{\mathbb{R}}^{{n}}}\Gamma(x,y)f(y)\ dy,\qquad
x\in{{\mathbb{R}}^{{n}}}.$

Then, $u_{f}\in\mathcal{A}({{\mathbb{R}}^{{n}}},{\mathcal{L}})$ and

${\mathcal{L}}_{a}u_{f}=-f\mbox{ in }{{\mathbb{R}}^{{n}}}.$

We will prove this theorem in the next section. Here, by using a result in
[6], we show a consequence of Theorem 1.5.

######  Theorem 1.6.

Let $f,u:{{\mathbb{R}}^{{n}}}\longrightarrow{\mathbb{R}}$ be compactly
supported continuous functions. Then,

${\mathcal{L}}_{a}u=-f\mbox{ in }{{\mathbb{R}}^{{n}}}$

if and only if

${\mathcal{L}}u=-f\mbox{ in }\mathcal{D}^{\prime}({{\mathbb{R}}^{{n}}}).$

###### Proof.

By the previous Theorem 1.5,

${\mathcal{L}}_{a}u=-f\mbox{ in }{{\mathbb{R}}^{{n}}}$

if and only if

${\mathcal{L}}_{a}(u-u_{f})=0\mbox{ in }{{\mathbb{R}}^{{n}}}.$

Then, by Corollary 3.4 in [6], $u-u_{f}\in
C^{\infty}({{\mathbb{R}}^{{n}}},{\mathbb{R}})$ and

${\mathcal{L}}(u-u_{f})=0$

in the classical sense (and vice versa). Since ${\mathcal{L}}$ is
hypoelliptic, this is equivalent to say that

${\mathcal{L}}(u-u_{f})=0\mbox{ in
}\mathcal{D}^{\prime}({{\mathbb{R}}^{{n}}}),$

or that

(1.13) ${\mathcal{L}}(u)={\mathcal{L}}(u_{f})\mbox{ in
}\mathcal{D}^{\prime}({{\mathbb{R}}^{{n}}}).$

On the other hand, $\Gamma$ being a fundamental solution of ${\mathcal{L}}$,
${\mathcal{L}}(u_{f})=-f\mbox{ in
}\mathcal{D}^{\prime}({{\mathbb{R}}^{{n}}})$. Then, (1.13) can be written as
follows:

${\mathcal{L}}u=-f\mbox{ in }\mathcal{D}^{\prime}({{\mathbb{R}}^{{n}}}).$

This completes the proof. ∎

###  1.5. Bibliographical note

In recent years asymptotic mean value formulas characterizing classical or
viscosity solutions to linear and nonlinear second order Partial Differential
Equations have been proved by many authors; we refer to [6, 12, 11, 8, 5, 10,
7, 13, 3]. In those papers one can find quite exhaustive bibliography on this
subject.

We would also like to quote the papers [4] and [9] where the notion of
asymptotic sub-harmonic function is introduced in sub-Riemannian settings to
extend classical results by Blaschke, Privaloff, Reade and Saks.

##  2\. Proof of Theorem 1.5

For the readers’ convenience, we split this section in two subsections.

### 2.1.

Let $G$ be a compact subset of ${{\mathbb{R}}^{{n}}}$ and let $r>0$. Define

(2.1) $G_{r}:=\bigcup_{x\in G}\Omega_{r}(x).$

Then, we have the following lemma.

######  Lemma 2.1.

For every compact set $G\subseteq{{\mathbb{R}}^{{n}}}$ and for every $r>0$,
the set $\overline{G}_{r}$ is compact.

###### Proof.

It is enough to prove that $G_{r}$ is bounded. We argue by contradiction and
assume that $G_{r}$ is not bounded. Then, there exists a sequence $(z_{n})$ in
$G_{r}$ such that

$|z_{n}|\longrightarrow\infty.$

By the very definition of $G_{r}$, for every $n\in{\mathbb{N}}$, there exists
$x_{n}\in G$ such that $z_{n}\in\Omega_{r}(x_{n})$. This means that

$\Gamma(x_{n},z_{n})>\frac{1}{r}.$

As a consequence,

$\frac{1}{r}<\Gamma(x_{n},z_{n})\leq\sup_{x\in G}\Gamma(x,z_{n}),$

so that, by the assumption (iii) related to $\Gamma$

$0<\frac{1}{r}\leq\lim_{n\longrightarrow\infty}\left(\sup_{x\in
G}\Gamma(x,z_{n})\right)=0.$

This contradiction shows that $G_{r}$ is bounded. ∎

### 2.2.

In this subsection we prove the following lemma.

######  Lemma 2.2.

Let $G$ be a compact subset of ${{\mathbb{R}}^{{n}}}$ and let $r>0$. Then,
there exists a positive constant $C_{r}(G)$ such that

(2.2) $\sup_{x\in G}Q_{r}(x)\leq C_{r}(G).$

###### Proof.

Keeping in mind the definition of $Q_{r}(x)$ (see (1.11)) for every $x\in G$
we get

$\displaystyle Q_{r}(x)$ $\displaystyle\leq$
$\displaystyle\frac{\alpha+1}{r^{\alpha+1}}\int_{0}^{r}\rho^{\alpha}\left(\int_{\Omega_{r}(x)}\Gamma(x,y)\
dy\right)\ d\rho$ $\displaystyle\leq$ $\displaystyle(\mbox{by}\
\eqref{2.1})\quad\frac{\alpha+1}{r^{\alpha+1}}\int_{0}^{r}\rho^{\alpha}\left(\int_{G_{r}}\Gamma(x,y)\
dy\right)\ d\rho.$

On the other hand, if $\varphi\in
C_{0}^{\infty}({{\mathbb{R}}^{{n}}},{\mathbb{R}})$ is such that $\varphi=1$ on
$G_{r},$ $\varphi\geq 0$ (such a function exists thanks to Lemma 2.1), we have

$\displaystyle\int_{G_{r}}\Gamma(x,y)\ dy$ $\displaystyle\leq$
$\displaystyle\int_{{{\mathbb{R}}^{{n}}}}\varphi(y)\Gamma(x,y)\ dy$
$\displaystyle\leq$ $\displaystyle\sup_{x\in
G}\int_{{{\mathbb{R}}^{{n}}}}\varphi(y)\Gamma(x,y)\ dy$ $\displaystyle=$
$\displaystyle C_{\varphi}(G).$

Using this estimate in (2.2) we obtain

$\displaystyle\sup_{x\in G}Q_{r}(x)$ $\displaystyle\leq$ $\displaystyle
C_{\varphi}(G)\frac{\alpha+1}{r^{\alpha+1}}\int_{0}^{r}\rho^{\alpha}\ d\rho$
$\displaystyle=$ $\displaystyle C_{\varphi}(G):=C_{r}(G).$

∎

######  Remark 2.3.

Since $Q_{\rho}(x)\subseteq Q_{r}(x)$ for every $\rho\in]0,r[$, we can assume

$C_{\rho}(G)\leq C_{r}(G)$

for every $0<\rho<r$.

### 2.3.

Now, we show a kind of continuity property of the $\Omega_{r}(x)$ balls with
respect to the Euclidean topology. Precisely, we prove the following lemma.

######  Lemma 2.4.

For every $x\in{{\mathbb{R}}^{{n}}}$ and for every $R>0$ there exists $r>0$
such that

$\Omega_{r}(x)\subseteq B(x,R).$

###### Proof.

We still argue by contradiction and assume the existence of $R>0$ such that
$\Omega_{r}(x)\nsubseteq B(x,R)$ for every $r>0$. Then, if $(r_{n})$ is a
sequence of real positive numbers such that $r_{n}\searrow 0$, for every
$n\in{\mathbb{N}}$ there exists $y_{n}\in\Omega_{r_{n}}(x)$ such that

$y_{n}\notin B(x,R).$

This means

$y_{n}\notin B(x,R)\quad\mbox{ and }\quad\Gamma(x,y_{n})>\frac{1}{r_{n}}.$

Since $\Gamma(x,y)\longrightarrow 0$ as $y\longrightarrow\infty$ and
$\dfrac{1}{r_{n}}\longrightarrow\infty$, the sequence $(y_{n})$ is bounded. As
a consequence, we may assume

$\lim_{n\longrightarrow\infty}y_{n}=y^{\ast}$

for a suitable $y^{\ast}\in{{\mathbb{R}}^{{n}}}.$ Then $y^{\ast}\notin
B(x,R).$ In particular $y\neq x$ so that $\Gamma(x,y)<\infty.$ On the other
hand,

$\Gamma(x,y^{\ast})=\lim_{n\longrightarrow\infty}\Gamma(x,y_{n})\geq\lim_{n\longrightarrow\infty}\frac{1}{r_{n}}=\infty.$

This contradiction proves the lemma. ∎

From the previous lemma we obtain the following corollary.

######  Corollary 2.5.

Let $f:{{\mathbb{R}}^{{n}}}\longrightarrow{\mathbb{R}}$ be a continuous
function. Then, for every $x\in{{\mathbb{R}}^{{n}}},$

$\sup_{y\in\Omega_{r}(x)}|f(y)-f(x)|\longrightarrow 0\mbox{ as
}r\longrightarrow 0.$

###### Proof.

Since $f$ is continuous at $x$, for every $\varepsilon>0$ there exists $R>0$
such that

$\sup_{y\in B(x,r)}|f(y)-f(x)|<\varepsilon.$

By the previous lemma, there exists $r_{0}>0$ such that
$\Omega_{r_{0}}(x)\subseteq B(x,r).$ Then, for every $r<r_{0}$,

$\sup_{y\in\Omega_{r}(x)}|f(y)-f(x)|\leq\sup_{y\in\Omega_{r_{0}}(x)}|f(y)-f(x)|\leq\sup_{y\in
B(x,r)}|f(y)-f(x)|<\varepsilon.$

We have so proved that for every $\varepsilon>0$ there exists $r_{0}>0$ such
that

$\sup_{y\in\Omega_{r}(x)}|f(y)-f(x)|<\varepsilon$

for every $r<r_{0}.$ Hence,

$\lim_{r\longrightarrow 0}\left(\sup_{y\in\Omega_{r}(x)}|f(y)-f(x)|\right)=0.$

∎

### 2.4.

Let $f$ as in Theorem 1.5 and, to simplify the notation, let us denote $u_{f}$
by $u$. The aim of this subsection is to prove the following identity:

(2.4) $u(x)=M_{r}(u)(x)+N_{r}(f)(x)\quad\forall x\in{{\mathbb{R}}^{{n}}}.$

To this end we choose a sequence $(f_{p})$ in
$C_{0}^{\infty}({{\mathbb{R}}^{{n}}},{\mathbb{R}})$ with the following
properties:

  * (i)

there exists a compact set $K\subseteq{{\mathbb{R}}^{{n}}}$ such that
$\mathrm{supp\,}f\subseteq K$ and $\mathrm{supp\,}f_{p}\subseteq K$ for every
$p\in{\mathbb{N}}$;

  * (ii)

$\sup_{K}|f_{p}-f|\longrightarrow 0\mbox{ as }p\rightarrow\infty.$

For simplicity reasons, let us put $u_{p}=u_{f_{p}}$, i.e.,

$\displaystyle u_{p}(x)=\int_{{\mathbb{R}}^{{n}}}\Gamma(x,y)f_{p}(y)\
dy=\int_{K}\Gamma(x,y)f_{p}(y)\ dy.$

Then, by Lebesgue’s dominated convergence Theorem,

$\displaystyle
u(x)=\lim_{p\longrightarrow\infty}u_{p}(x)=\int_{K}\Gamma(x,y)\lim_{p\longrightarrow\infty}f_{p}(y)\
dy,$

for every $x\in{{\mathbb{R}}^{{n}}}.$ Actually, we have a stronger result. For
every compact set $G\subseteq{{\mathbb{R}}^{{n}}}$,

$\displaystyle\sup_{G}|u_{p}-u|$ $\displaystyle\leq$ $\displaystyle\sup_{x\in
G}\left|\int_{K}\Gamma(x,y)(f_{p}(y)-f(y))\ dy\right|$ $\displaystyle\leq$
$\displaystyle\sup_{K}|f_{p}-f|\sup_{x\in G}\int_{K}\Gamma(x,y)\ dy$
$\displaystyle=$ $\displaystyle C(G,K)\sup_{K}|f_{p}-f|.$

We explicitly observe that $C(G,K)$ is a strictly positive finite constant.

Hence,

(2.5) $\sup_{G}|u_{p}-u|\longrightarrow 0\mbox{ as }p\longrightarrow\infty.$

Moreover, for every $p\in{\mathbb{N}}$,

$u_{p}\in C^{\infty}({{\mathbb{R}}^{{n}}},{\mathbb{R}})\quad\mbox{ and
}\quad{\mathcal{L}}u_{p}=-f_{p}.$

Then, by identity (1.8),

$\displaystyle u_{p}(x)$ $\displaystyle=$ $\displaystyle
M_{r}(u_{p})(x)-N_{r}({\mathcal{L}}u_{p})(x)$ $\displaystyle=$ $\displaystyle
M_{r}(u_{p})(x)+N_{r}(f_{p})(x)$

for every $p\in{\mathbb{N}}.$

We have already noticed that $u_{p}(x)\longrightarrow u(x)$ as
$p\longrightarrow\infty.$

To prove (2.4) we now show that

(2.6) $\displaystyle\lim_{p\longrightarrow\infty}M_{r}(u_{p})(x)=M_{r}(u)(x)$

and

(2.7) $\displaystyle\lim_{p\longrightarrow\infty}N_{r}(f_{p})(x)=N_{r}(f)(x).$

For every $x\in{{\mathbb{R}}^{{n}}}$ we have

$\displaystyle|M_{r}(u_{p})(x)-M_{r}(u)(x)|$ $\displaystyle=$
$\displaystyle|M_{r}(u_{p}-u)(x)|$ $\displaystyle\leq$
$\displaystyle\sup_{\Omega_{r}(x)}|u_{p}-u|M_{1}(1)(x)$ $\displaystyle=$
$\displaystyle\sup_{\Omega_{r}(x)}|u_{p}-u|.$

Since $\overline{\Omega_{r}(x)}$ is compact (see Lemma 2.1), and keeping in
mind (2.5), the last right hand side goes to zero as $p\longrightarrow\infty.$
Then,

$\displaystyle|M_{r}(u_{p})(x)-M_{r}(u)(x)|\longrightarrow 0\mbox{ as
}p\longrightarrow\infty,$

proving (2.6).

Let us now prove (2.7). For every $x\in{{\mathbb{R}}^{{n}}}$, we have

$\displaystyle|N_{r}(f_{p})(x)-N_{r}(f)(x)|$ $\displaystyle\leq$
$\displaystyle|N_{r}(|f_{p}-f|)(x)|$ $\displaystyle\leq$
$\displaystyle\sup_{K}|f_{p}-f|Q_{r}(x).$

Then, for every compact set $G\subseteq{{\mathbb{R}}^{{n}}}$,

$\displaystyle\sup_{G}|N_{r}(f_{p})-N_{r}(f)|$ $\displaystyle\leq$
$\displaystyle\sup_{K}|f_{p}-f|\sup_{x\in G}|Q_{r}(x)|$ $\displaystyle\leq$
$\displaystyle(\mbox{by}\ \eqref{2.2})\quad C_{r}(G)\sup_{K}|f_{p}-f|.$

So we have proved that $(N_{r}(f_{p}))$ is uniformly convergent to $N_{r}(f)$
on every compact subset of ${{\mathbb{R}}^{{n}}}$. This, in particular,
implies (2.7).

### 2.5.

In this subsection we complete the proof of Theorem 1.5. To this end we first
remark that, thanks to (2.4), for every $x\in{{\mathbb{R}}^{{n}}},$ we have

$\displaystyle\frac{M_{r}(u)(x)-u(x)}{Q_{r}(x)}=-\frac{N_{r}(f)(x)}{Q_{r}(x),}$

so that, as $f(x)$ is constant with respect to $y\in\Omega_{r}(x)$,

$\displaystyle\left|\frac{M_{r}(u)(x)-u(x)}{Q_{r}(x)}+f(x)\right|$
$\displaystyle=$
$\displaystyle\frac{1}{Q_{r}(x)}\left|{N_{r}(f(x)-f})(x)\right|$
$\displaystyle\leq$
$\displaystyle\sup_{y\in\Omega_{r}(x)}|f(u)-f(y)|Q_{r}(x).$

By Corollary 2.5 and Remark 2.3, the left hand side of the previous inequality
goes to zero as $r\longrightarrow 0$. Hence,

$\displaystyle\lim_{r\longrightarrow
0}\frac{M_{r}(u)(x)-u(x)}{Q_{r}(x)}=-f(x)$

for every $x\in{{\mathbb{R}}^{{n}}}$. This completes the proof of Theorem 1.5.

## Declarations

  * -

Conflict of interest: The authors declare that they have no conflict of
interest.

  * -

Data availability: Data sharing not applicable to this article as no datasets
were generated or analysed during the current study.

  * -

Funding: The first author has been partially supported by the Gruppo Nazionale
per l’Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of
the Istituto Nazionale di Alta Matematica (INdAM).

## References

  * [1] S. Biagi and A. Bonfiglioli.  The existence of a global fundamental solution for homogeneous Hörmander operators via a global lifting method.  Proc. Lond. Math. Soc. (3), 114(5):855–889, 2017.
  * [2] S. Biagi, A. Bonfiglioli, and M. Bramanti.  Global estimates for the fundamental solution of homogeneous Hörmander operators.  Ann. Mat. Pura Appl. (4), 201(4):1875–1934, 2022.
  * [3] P. Blanc, F. Charro, J. J. Manfredi, and J. D. Rossi.  Asymptotic mean-value formulas for solutions of general second-order elliptic equations.  Adv. Nonlinear Stud., 22(1):118–142, 2022.
  * [4] A. Bonfiglioli and E. Lanconelli.  Subharmonic functions in sub-Riemannian settings.  J. Eur. Math. Soc. (JEMS), 15(2):387–441, 2013.
  * [5] F. Ferrari.  Mean value properties of fractional second order operators.  Commun. Pure Appl. Anal., 14(1):83–106, 2015.
  * [6] C. E. Gutiérrez and E. Lanconelli.  Classical viscosity and average solutions for PDE’s with nonnegative characteristic form.  Atti Accad. Naz. Lincei Cl. Sci. Fis. Mat. Natur. Rend. Lincei (9) Mat. Appl., 15(1):17–28, 2004.
  * [7] M. Ishiwata, R. Magnanini, and H. Wadade.  A natural approach to the asymptotic mean value property for the $p$-Laplacian.  Calc. Var. Partial Differential Equations, 56(4):Paper No. 97, 22, 2017.
  * [8] B. Kawohl, J. J. Manfredi, and M. Parviainen.  Solutions of nonlinear PDEs in the sense of averages.  J. Math. Pures Appl. (9), 97(2):173–188, 2012.
  * [9] A. E. Kogoj and G. Tralli.  Blaschke, Privaloff, Reade and Saks theorems for diffusion equations on Lie groups.  Potential Anal., 38(4):1103–1122, 2013.
  * [10] P. Lindqvist and J. J. Manfredi.  On the mean value property for the $p$-Laplace equation in the plane.  Proc. Amer. Math. Soc., 144(1):143–149, 2016.
  * [11] J. J. Manfredi, M. Parviainen, and J. D. Rossi.  An asymptotic mean value characterization for a class of nonlinear parabolic equations related to tug-of-war games.  SIAM J. Math. Anal., 42(5):2058–2081, 2010.
  * [12] J. J. Manfredi, M. Parviainen, and J. D. Rossi.  An asymptotic mean value characterization for $p$-harmonic functions.  Proc. Amer. Math. Soc., 138(3):881–889, 2010.
  * [13] B. Mebrate and A. Mohammed.  Harnack inequality and an asymptotic mean-value property for the Finsler infinity-Laplacian.  Adv. Calc. Var., 14(3):365–382, 2021.
  * [14] P. Pizzetti.  Sulla media dei valori che una funzione dei punti dello spazio assume alla superficie di una sfera.  Atti Accad. Naz. Lincei Rend. Cl. Sci. Fis. Mat. Nat. Serie V, 18:182–185, 1909.