# Qubits on programmable geometries with a trapped-ion quantum processor

Qiming Wu<EMAIL_ADDRESS>Department of Physics, University of California,
Berkeley, CA, USA  Yue Shi  Department of Physics, Princeton University,
Princeton, NJ, USA  Jiehang Zhang<EMAIL_ADDRESS>1\. School of
Physical Sciences, University of Science and Technology of China, Hefei
230026, China  2\. Shanghai Research Center for Quantum Science and CAS Center
for Excellence in Quantum Information and Quantum Physics, University of
Science and Technology of China, Shanghai 201315, China  3\. Hefei National
Laboratory, University of Science and Technology of China, Hefei 230088, China

###### Abstract

Geometry and dimensionality have played crucial roles in our understanding of
the fundamental laws of nature, with examples ranging from curved space-time
in general relativity to modern theories of quantum gravity. In quantum many-
body systems, the entanglement structure can change if the constituents are
connected differently, leading to altered bounds for correlation growth and
difficulties for classical computers to simulate large systems. While a
universal quantum computer can perform digital simulations, an analog-digital
hybrid quantum processor offers advantages such as parallelism. Here, we
engineer a class of high-dimensional Ising interactions using a linear one-
dimensional (1D) ion chain with up to 8 qubits through stroboscopic sequences
of commuting Hamiltonians. In addition, we extend this method to non-commuting
circuits and demonstrate the quantum XY and Heisenberg models using Floquet
periodic drives with tunable symmetries. The realization of higher dimensional
spin models offers new opportunities ranging from studying topological phases
of matter or quantum spin glasses to future fault-tolerant quantum
computation.

††preprint: APS/123-QED

## Introduction

Quantum computers are expected to solve certain problems intractable in their
classical counterparts. A prominent example given by Feynman is to simulate
quantum many-body problems, ranging from strongly correlated materials [1, 2]
to quantum chemistry [3, 4] and the dynamical evolution of quantum many-body
systems [5, 6]. Quantum simulations can be performed in a well-controlled
platform with interactions native to the physical system: neutral atoms in
optical lattices or optical cavities are used for models with nearest-neighbor
hopping [7, 8] or all-to-all couplings [9, 10]; whereas platforms such as
trapped ions, Rydberg atoms, nitrogen-vacancy spin defects in diamonds, and
molecules are well-suited for simulating long-range interacting physics [11,
12, 13, 14]. The quest for higher flexibility and programmability has
attracted much attention recently, where the interactions are encoded beyond
these regimes. Qubits interacting on non-Euclidean geometries [15, 16] can
have eigenenergies and correlations different from their Euclidean
counterparts, with a prominent motivation to simulate physics in models of
quantum gravity [17]. While in solid-state systems such as superconducting
qubits, the interactions are hard-wired by design [15, 16], in atomic systems
such as cavity quantum electrodynamics (c-QED), new modulation schemes have
shown powerful tunability across different types of graphs without the need to
physically re-wire [18, 19].

Trapped ion systems afford such programmability: collective phonon modes
mediate the interactions, and modulated laser drive fields can generate
tunable non-local couplings. Quantum simulation experiments performed so far
have realized long-range approximate power-law decaying interactions generated
in the “dispersive coupling“ regime, where global laser fields couple to all
phonon modes[20, 21]. Though successful in demonstrating non-equilibrium
dynamics hard to tract by classical computers [20], the limited tunability
(finite interaction range) and the residual spin-phonon coupling limits future
scalability [22]. In recent years, new theory proposals have emerged to
engineer coupling graphs such as nearest-neighbor interactions in 1D and 2D
systems [23, 24, 25, 26], followed by experimental demonstrations of time-
reversal-breaking dynamics in triangular ladders [27].

Figure 1: Illustration of the trapped-ion hybrid quantum information
processor. (a) A chain of $\mathrm{{}^{9}Be^{+}}$ ions are confined in a RF
Paul trap in a linear configuration. A pair of Raman beams provides beatnote
frequencies symmetrically detuned from the qubit transition
$\omega_{0}\pm\mu$, mediating the spin-spin interactions. The blue (red)
curves denote near equal (zero) coupling strength between the connected ion
pairs. (b) With 6 spins, the interaction illustrated in (a) can be mapped to
nearest-neighbor interaction on a sphere. (c) These interactions are generated
through a weighted combination of coupling to modes 1,3,5 (blue solid lines).
The blue arrows visualize the corresponding normal mode eigenvectors.

In this work, we experimentally demonstrate, for the first time to our
knowledge, that global laser fields can be modulated to encode qubits in high-
dimensional non-Euclidean spaces, exemplified by up to 8 qubits interacting on
a four-dimensional hyper-sphere, an example of positively-curved elliptical
space. We study the temporal dynamics of many-body quantum systems out of
equilibrium without intrinsic spin-phonon errors under programmable
geometries. We experimentally tune between commuting Ising-type and non-
commuting Heisenberg-type Hamiltonian models with tunable coupling strengths
and durations and study the non-equilibrium dynamics with extensions to
Floquet models [6, 28]. Our hybrid analog-digital quantum simulation method
offers unique parallelism, thus complementing fully digital quantum computers
[29, 30].

## Tunable coupling graphs by frequency-domain modulation

The non-local Coulomb coupling in principle allows arbitrary connectivity in a
trapped ion system, but requires individual addressing and classical control
techniques to entangle the desired spin states without excess phonon
occupation [29]. With global laser fields, previous experiments have focused
on either collectively coupling to all modes [20, 21] or tuning the frequency
to near the center-of-mass mode [31, 12], while novel theories propose the
combination of selected phonon modes for engineering tunable graphs [25].
Here, we modulate the global laser fields in the frequency domain, where the
spin-dependent displacements of the ions can be decomposed into couplings with
different normal motional modes, forming different correlation patterns. By
predominantly coupling to a single mode, the phonon-mediated interaction
inherits a structure from the normal mode eigenvectors (Fig. 1). The spin-spin
interactions can be represented by a non-local weighted graph, where the
weights are generated by a linear combination of the coupling strength to the
individual modes, yielding tunable couplings with $(N-1)$ degrees of freedom
in an $N$-ion chain [25].

Here, we focus on engineering qubit interconnections beyond the physical
dimension. A key to our scheme is to reverse the sign of unwanted interaction
between certain ion pairs when coupling to certain modes, and by linking
together the global interactions, the average interactions between these pairs
vanish, forming connectivities for the desired geometries. For example,
Fig.1(c) illustrates that six qubits on a sphere can be encoded through
coupling to the first, third, and fifth axial modes with appropriate weights.
Moreover, we can change the coupling patterns and explore a more general class
of unique interactions by modifying the motional spectrum of the ion chain,
such as generating tree-like geometries and qubits on a torus under an
anharmonic trapping potential [32].

Figure 2: Illustration of qubit interaction on a plaquette with four ions. (a)
The interaction is generated with equal participation of the first and third
modes. (b) Calibration of the BSB time evolution on the modes used in
Hamiltonian engineering. Blue dots denote the measurements of the average
population on the $\ket{\downarrow}$ state and blue solid lines are numerical
simulations of BSB floppings, assuming $\bar{n}=0.1$. (c) Time evolution of 4
qubit interaction on a plaquette. The measurement of site-resolved spin
magnetizations is depicted as solid lines. The numerical simulations are
extrapolated through dashed lines, serving as a visual comparison to the
experimental results. The error bars denote one standard deviation.

## Hybrid quantum information processor

Our hybrid quantum information processor contains a linear chain of
$\mathrm{{}^{9}Be^{+}}$ ions confined in a radio-frequency Paul trap (Fig.
1(b)), held in their equilibrium positions with the Coloumb interactions
balancing the external harmonic confinement. The coupled ion motions in the
axial direction can be decomposed to $n$ eigenmodes, where $n$ is the ion
number. We encode the qubits in the ground state hyperfine levels:
$\ket{\downarrow}_{z}\equiv~{}^{2}S_{1/2}\ket{F=2,m_{F}=2}$ and
$\ket{\uparrow}_{z}\equiv~{}^{2}S_{1/2}\ket{F=1,m_{F}=1}$, with a frequency
difference of $\sim$1.018 GHz under a magnetic field of $\sim$120 G. The
experiment begins by cooling all the motional modes of the ions to near the
ground state using continuous Raman sideband cooling [33]. We calibrate the
spin-motion coupling by driving the blue sideband (BSB) transitions on the
desired modes used for mediating the Ising interaction. Fig. 2(b)) shows the
BSB calibration in a 4-ion chain. We apply a pair of global Raman laser beams,
one of which is bi-chromatic such that the beat-note frequencies are
$\omega_{0}\pm\mu$, symmetrically detuned from the qubit transition
$\omega_{0}$. One of the beams is linear polarized ($\pi$), the other has
linear polarization (equal $\sigma+$/$\sigma-$) plus a small non-zero term
circular polarized term to suppress the Stark shift to less than 500 Hz [34],
limited by the Zeeman-sensitive qubit frequency drifts under fluctuating
magnetic fields. The Raman beatnote generates a spin-dependent optical dipole
force via coupling to the axial motions of the ions, which mediates the spin-
spin interaction [35]. We choose axial modes because of their near-even
spacings well-separated in the frequency domain, easing frequency modulation
schemes [25]. We vary the detunings over the mode spectrum so that positive
and negative coupling terms could exist simultaneously [32]. When the laser
detuning $\mu$ is close to one specific motional mode frequency $\omega_{M}$
with $|\mu-\omega_{M}|\ll|\mu-\omega_{m}|_{M\neq m}$, we generate the mode-
specific Ising interaction with a predominant coupling to one mode:

$H_{M}=\sum_{i<j}J_{ij}(\mu_{M})\sigma_{i}^{\phi}\sigma_{j}^{\phi},$ (1)

where $J_{ij}(\mu_{M})$ is the spin-spin interaction term when coupling to
mode $\omega_{M}$ and $J_{ij}(\mu_{M})\propto
b_{i,M}b_{j,M}/(\mu_{M}-\omega_{M})$. $b_{i,M}$ is the eigenvector component
between the $i$-th ion and mode M,
$\sigma_{i}^{\phi}=\sigma_{i}^{x}\cos{\phi}-\sigma_{i}^{y}\sin{\phi}$ is the
Pauli operator for the $i$-th ion with a tunable phase $\phi$, and the reduced
Planck’s constant $\hbar$ is set to 1. In previous quantum simulation
experiments [22], the spin-phonon coupling term is neglected in the dispersive
regime when the detuning $|\mu-\omega_{m}|\gg\eta_{i,m}\Omega_{i}$ at the cost
of residual errors [11]. However, by making use of the axial mode structure,
we mitigate this effect by choosing the detunings to our desired coupling mode
$\delta_{M}=\mu-\omega_{M}$ to be a common divisor to detuning of all other
modes $\delta_{m}=\mu-\omega_{m}=k\delta_{M}(k\in\mathbb{Z})$. An additional
benefit is a speed-up of entangling gates compared to the adiabatic condition
[36].

We first choose the appropriate laser phase to generate the
$\sigma_{i}^{x}\sigma_{j}^{x}$ type Hamiltonian. We sequentially apply Ising
interactions $H_{a},H_{b},...,H_{N}$ with tunable time duration of
$\tau_{a},\tau_{b},...,\tau_{N}$, where every $H_{m}$ generates a
$\sigma_{x}\sigma_{x}$ interaction with different weights by coupling to the
$m$-th mode (Eq.1). In the case of all commuting terms, we create the
following effective Hamiltonian:

$\displaystyle{H_{eff}}$
$\displaystyle=\frac{\sum_{m}H_{m}\tau_{m}}{\sum_{m}\tau_{m}}=\frac{\sum_{i<j}\sum_{m}J_{ij}(\mu_{m})\tau_{m}}{\sum_{m}\tau_{m}}\sigma_{i}^{x}\sigma_{j}^{x}$
(2) $\displaystyle=\sum_{i<j}\bar{J_{ij}}\sigma_{i}^{x}\sigma_{j}^{x},$

where we denote $\bar{J_{ij}}$ to be the effective coupling strength and
$\bar{J_{ij}}\propto\sum_{m}b_{i,m}b_{j,m}\tau_{M}/(\mu_{m}-\omega_{m})$.
Through concatenating continuously tunable (analog) blocks with different
laser detunings, the system is subject to a programmable quantum evolution
characterized by $\bar{J_{ij}}$. Such couplings correspond to the nearest-
neighbor interactions in plaquette, sphere, and hypersphere for 4, 6, and 8
ions, respectively. At the end of each global interaction layer, we perform
single-site readouts of the qubit magnetizations by illuminating the ion chain
with resonant detection light and collecting the spin-dependent fluorescence.

Figure 3: Realization of Ising interaction in a higher dimension sphere.
(a)(c) Illustration of qubit interaction of nearest-neighbor interaction on a
sphere and hypersphere, with 6 and 8 qubits, respectively. (b)(d) Measurement
of the average spin magnetization $\overline{\langle\sigma_{z}\rangle}$. Solid
lines extrapolate the spin magnetizations measured in the experiment, and
solid rhombuses are the experimental data. Dash lines extrapolate the
numerical calculations as a side-by-side comparison to the experimental
result. Error bars denote one standard deviation.

## Interaction in higher-dimension with 1D chain

We first show that a 4-ion chain can be encoded in the smallest setting of 2D
nearest-neighbor interactions in a square plaquette (Fig. 2(a)). We prepare
the ions on the $\ket{\downarrow}_{z,N}$ state and apply successive $H_{1}$
and $H_{3}$ interactions iteratively. $H_{1}$ couples to the first mode
generating all-to-all couplings $J_{ij}=\mathcal{J}$ for all $i,j$ and $H_{3}$
couples to the third mode with reversed coupling in two of the pairs
$J_{14},J_{23}=-\mathcal{J}$. The laser detunings of these two interactions
are $\delta_{1}/2\pi=107.6$ kHz from the first mode and
$\delta_{3}/2\pi=-71.14$ kHz from the third mode, with interaction time
$\tau_{1}=2\pi/\delta_{1}=9.29~{}\mu s$ and
$\tau_{3}=2\pi/\delta_{3}=14.06~{}\mu s$. After the concatenated unitary
evolution, the effective Hamiltonian generates a 2D nearest-neighbor
interaction, which we measure via site-resolve spin magnetization
$\langle\sigma_{z}\rangle$ (Fig. 2(c)). The understand the deviations of the
experimental results from the implemented Hamiltonian, we perform numerical
simulations of the implemented Hamiltonian taking into account systematic
experimental offsets, including qubit frequency and laser intensity
gradient[32], and find reasonable agreements between the simulation and the
experimental data.

We next generalize the method to higher dimensional spaces. With six ions, a
minimum graph corresponding to a sphere with evenly distributed qubits (two on
the north and south poles, and four evenly distributed on the equator) can be
generated by removing selected vertices (red lines in Fig. 1(a)), which is
similar to removing the “cross” vertices in our previous case of 4 qubits. As
a primary example of non-Euclidean space with positive curvature, we can
extend this to 8 qubits and encode them on a 4-sphere. The protocols are
periodic drives applied to the system with coupling to the odd number modes,
namely modes 1,3,5 in the 6-ion chain and modes 1,3,5,7 in the 8-ion chain.
The detunings from the modes and layer durations can be found in [32].
Similarly, we capture the temporal dynamics of the system by measuring spin
magnetizations. Fig. 3(b)(c) show the average spin dynamics as a function of
the gate layers applied. Our results are consistent with numerical
calculation, considering the actual Hamiltonian with the experimental offset
(see S.I. [32]).

## Floquet XY and Heisenberg models

Figure 4: Realization of Floquet XY and XYZ Hamiltonians. Illustration of
(a)XY (c) Heisenberg XYZ model in a 2D nearest neighbor interaction with four
ions. The (b) navy and (d) green rhombuses denote the evolution of the
$\langle\sigma_{z}\rangle$ as a function of the applied analog blocks under an
isotropic (b) XY (d) XYZ 2D nearest neighbor interaction with four ions. The
error bars denote one standard deviation. The dashed lines extrapolate the
numerical simulation of the ideal quantum evolution, considering the
experimental offsets.

The flexible control of layers of global interaction makes our schemes
favorable to extend to other spin degrees of freedom and even non-commuting
spin-spin interactions, which native Hamiltonians do not exist in our
platform. As an example, we implement the same 2D interaction on a plaquette
on a Floquet XY model by applying a periodic drive with $T=2(t_{1}+t_{3})$ to
the system, which consists of four successive global interaction layers:

$H=\left\\{\begin{array}[]{l}H_{1}^{XX}=\sum_{i<j}J_{ij}^{(1)}\sigma_{i}^{x}\sigma_{j}^{x}\qquad
t=t_{1}\\\ H_{3}^{XX}=\sum_{i<j}J_{ij}^{(3)}\sigma_{i}^{x}\sigma_{j}^{x}\qquad
t=t_{3}\\\ H_{1}^{YY}=\sum_{i<j}J_{ij}^{(1)}\sigma_{i}^{y}\sigma_{j}^{y}\qquad
t=t_{1}\\\ H_{3}^{YY}=\sum_{i<j}J_{ij}^{(3)}\sigma_{i}^{y}\sigma_{j}^{y}\qquad
t=t_{3}\end{array}\right..$ (3)

Here, $H_{1}$ and $H_{3}$ are the two types of interaction patterns when
coupled to the first and third modes. A static quantum XY model
$H_{XY}=J_{ij}(\sigma_{i}^{x}\sigma_{j}^{x}+\sigma_{i}^{y}\sigma_{j}^{y})$ has
no dynamics for the ferromagnetic chain when the spins are initialized in
$\ket{\downarrow\downarrow\downarrow\downarrow}_{z}$ state (paramagnetic spin
chain). However, the system exhibits strong time-dependent dynamics under a
Floquet periodic drive with suitable interaction strength of each layer. Fig.
4(a) shows the experimental result of the emerging dynamics of the collective
spin magnetization. We find decent agreement with the numerical model
considering experimental offsets.

Furthermore, we can engineer Heisenberg Hamiltonians with tunable anisotropy
with periodic modulation. We present an experimental implementation of the
Heisenberg model to the same qubit interaction graph with four ions. To
generate the $\sigma_{z}\sigma_{z}$-type interaction, we drive the system into
a spin-rotating frame using a global $\pi/2$ y-rotation, followed by
$\sigma_{x}\sigma_{x}$ global interaction, and then a global $-\pi/2$
y-rotation. By that, a trotter step consists of 8 global operations[37]:

$\displaystyle U_{XYZ}(t)=$ $\displaystyle
e^{-iH_{1}^{XX}t_{1}}e^{-iH_{3}^{XX}t_{3}}e^{-iH_{1}^{YY}t_{1}}e^{-iH_{3}^{YY}t_{3}}$
(4) $\displaystyle
R_{y}(\pi/4)e^{-iH_{1}^{XX}t_{1}}e^{-iH_{3}^{XX}t_{3}}R_{y}^{\dagger}(\pi/4).$

Fig. 4 (b) shows the experimental result of the emerging dynamics of the
collective spin magnetization under the Floquet Heisenberg model. We compare
it with numerical simulations, validating that our protocol implements the
target Hamiltonian after up to 24 analog blocks.

## Conclusion and outlook

We have experimentally realized qubit interactions in a higher dimensional
space on a 1D ion chain of up to 8 qubits using the analog-digital hybrid
approach with global laser beams. Furthermore, we apply this protocol to
generate Floquet XY and Heisenberg models of the 2D nearest neighbor
interaction using up to 24 global interaction blocks with continuously tunable
parameters. Our results contribute to the engineering of multi-dimension
Hamiltonians using global interactions with analog-digital hybrid protocols,
which opens up the full potential of the trapped ion quantum simulator beyond
the long-range interaction regime, laying the foundations for more complex
geometries [26]. Combining our method with single-spin addressing and optical
tweezers [38], generating arbitrary Ising graphs at a lower engineering cost
than a universal quantum computer will be possible. This new platform will
allow us to perform quantum simulation of quantum systems with complex
geometries such as MaxCut problems[39], quantum spin glasses[40] and the 2D
Kagome lattice[41]. Furthermore, our scheme might inspire the direct encoding
of a quantum error correction circuit [42], with potential advantages in using
non-Euclidean geometries [43].

###### Acknowledgements.

We thank Luis Orozco, Dries Sels, Wenchao Ge, Ye Wang, and Zijian Ding for
helpful discussions and critical reading of the manuscript. This work was
performed at New York University, with financial support from the US
Department of Energy.

## References

  * Kennes et al. [2021] D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. Basov, A. N. Pasupathy, and A. Rubio, Nature Physics 17, 155 (2021).
  * Jafferis et al. [2022] D. Jafferis, A. Zlokapa, J. D. Lykken, D. K. Kolchmeyer, S. I. Davis, N. Lauk, H. Neven, and M. Spiropulu, Nature 612, 51 (2022).
  * Lanyon et al. [2010] B. P. Lanyon, J. D. Whitfield, G. G. Gillett, M. E. Goggin, M. P. Almeida, I. Kassal, J. D. Biamonte, M. Mohseni, B. J. Powell, M. Barbieri, et al., Nature chemistry 2, 106 (2010).
  * Quantum et al. [2020] G. A. Quantum, Collaborators*†, F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, et al., Science 369, 1084 (2020).
  * Zhang et al. [2017a] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, et al., Nature 543, 217 (2017a).
  * Geier et al. [2021] S. Geier, N. Thaicharoen, C. Hainaut, T. Franz, A. Salzinger, A. Tebben, D. Grimshandl, G. Zürn, and M. Weidemüller, Science 374, 1149 (2021).
  * Bloch et al. [2012] I. Bloch, J. Dalibard, and S. Nascimbene, Nature Physics 8, 267 (2012).
  * Gross and Bakr [2021] C. Gross and W. S. Bakr, Nature Physics 17, 1316 (2021).
  * Barontini et al. [2015] G. Barontini, L. Hohmann, F. Haas, J. Estève, and J. Reichel, Science 349, 1317 (2015).
  * Pedrozo-Peñafiel et al. [2020] E. Pedrozo-Peñafiel, S. Colombo, C. Shu, A. F. Adiyatullin, Z. Li, E. Mendez, B. Braverman, A. Kawasaki, D. Akamatsu, Y. Xiao, et al., Nature 588, 414 (2020).
  * Kim et al. [2010] K. Kim, M.-S. Chang, S. Korenblit, R. Islam, E. E. Edwards, J. K. Freericks, G.-D. Lin, L.-M. Duan, and C. Monroe, Nature 465, 590 (2010).
  * Britton et al. [2012] J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk, and J. J. Bollinger, Nature 484, 489 (2012).
  * Choi et al. [2017] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, et al., Nature 543, 221 (2017).
  * Kandala et al. [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Kollár et al. [2019] A. J. Kollár, M. Fitzpatrick, and A. A. Houck, Nature 571, 45 (2019).
  * Chen et al. [2023] A. Chen, H. Brand, T. Helbig, T. Hofmann, S. Imhof, A. Fritzsche, T. Kießling, A. Stegmaier, L. K. Upreti, T. Neupert, et al., Nature Communications 14, 622 (2023).
  * Ambjørn et al. [2004] J. Ambjørn, J. Jurkiewicz, and R. Loll, Physical review letters 93, 131301 (2004).
  * Bentsen et al. [2019] G. Bentsen, T. Hashizume, A. S. Buyskikh, E. J. Davis, A. J. Daley, S. S. Gubser, and M. Schleier-Smith, Physical review letters 123, 130601 (2019).
  * Periwal et al. [2021] A. Periwal, E. S. Cooper, P. Kunkel, J. F. Wienand, E. J. Davis, and M. Schleier-Smith, Nature 600, 630 (2021).
  * Zhang et al. [2017b] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, Nature 551, 601 (2017b).
  * Joshi et al. [2022] M. K. Joshi, F. Kranzl, A. Schuckert, I. Lovas, C. Maier, R. Blatt, M. Knap, and C. F. Roos, Science 376, 720 (2022).
  * Monroe et al. [2021] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, et al., Reviews of Modern Physics 93, 025001 (2021).
  * Rajabi et al. [2019] F. Rajabi, S. Motlakunta, C.-Y. Shih, N. Kotibhaskar, Q. Quraishi, A. Ajoy, and R. Islam, npj Quantum Information 5, 1 (2019).
  * Teoh et al. [2020] Y. H. Teoh, M. Drygala, R. G. Melko, and R. Islam, Quantum Science and Technology 5, 024001 (2020).
  * Shapira et al. [2020] Y. Shapira, R. Shaniv, T. Manovitz, N. Akerman, L. Peleg, L. Gazit, R. Ozeri, and A. Stern, Physical Review A 101, 032330 (2020).
  * Manovitz et al. [2020] T. Manovitz, Y. Shapira, N. Akerman, A. Stern, and R. Ozeri, PRX quantum 1, 020303 (2020).
  * Shapira et al. [2023] Y. Shapira, T. Manovitz, N. Akerman, A. Stern, and R. Ozeri, Physical Review X 13, 021021 (2023).
  * Scholl et al. [2022] P. Scholl, H. J. Williams, G. Bornet, F. Wallner, D. Barredo, L. Henriet, A. Signoles, C. Hainaut, T. Franz, S. Geier, et al., PRX Quantum 3, 020303 (2022).
  * Debnath et al. [2016] S. Debnath, N. M. Linke, C. Figgatt, K. A. Landsman, K. Wright, and C. Monroe, Nature 536, 63 (2016), ISSN 0028-0836.
  * Pino et al. [2021] J. M. Pino, J. M. Dreiling, C. Figgatt, J. P. Gaebler, S. A. Moses, M. Allman, C. Baldwin, M. Foss-Feig, D. Hayes, K. Mayer, et al., Nature 592, 209 (2021).
  * Martinez et al. [2016] E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl, P. Hauke, M. Dalmonte, T. Monz, P. Zoller, et al., Nature 534, 516 (2016).
  * wu [2] See supplemental materials.
  * Wu et al. [2023] Q. Wu, Y. Shi, and J. Zhang, Physical Review Research 5, 023022 (2023).
  * Wu [2023] Q. Wu, Ph.D. thesis, New York University (2023).
  * Porras and Cirac [2004] D. Porras and J. I. Cirac, Physical review letters 92, 207901 (2004).
  * Schäfer et al. [2018] V. Schäfer, C. Ballance, K. Thirumalai, L. Stephenson, T. Ballance, A. Steane, and D. Lucas, Nature 555, 75 (2018).
  * Arrazola et al. [2016] I. Arrazola, J. S. Pedernales, L. Lamata, and E. Solano, Scientific reports 6, 1 (2016).
  * Olsacher et al. [2020] T. Olsacher, L. Postler, P. Schindler, T. Monz, P. Zoller, and L. M. Sieberer, PRX Quantum 1, 020316 (2020).
  * Ebadi et al. [2022] S. Ebadi, A. Keesling, M. Cain, T. T. Wang, H. Levine, D. Bluvstein, G. Semeghini, A. Omran, J.-G. Liu, R. Samajdar, et al., Science (2022).
  * Gopalakrishnan et al. [2011] S. Gopalakrishnan, B. L. Lev, and P. M. Goldbart, Physical review letters 107, 277201 (2011).
  * Korenblit et al. [2012] S. Korenblit, D. Kafri, W. C. Campbell, R. Islam, E. E. Edwards, Z.-X. Gong, G.-D. Lin, L.-M. Duan, J. Kim, K. Kim, et al., New Journal of Physics 14, 095024 (2012).
  * Ryan-Anderson et al. [2022] C. Ryan-Anderson, N. Brown, M. Allman, B. Arkin, G. Asa-Attuah, C. Baldwin, J. Berg, J. Bohnet, S. Braxton, N. Burdick, et al., arXiv preprint arXiv:2208.01863 (2022).
  * Kitaev [2003] A. Y. Kitaev, Annals of physics 303, 2 (2003).


# Synchronization in Networks with Nonlinearly Delayed Couplings on Example of
Neural Mass Model

Sergei A. Plotnikov  Institute for Problems of Mechanical Engineering, Saint-
Petersburg, Russia, (e-mail<EMAIL_ADDRESS>
###### Abstract

The problem of synchronization in heterogeneous networks of linear systems
with nonlinear delayed diffusive coupling is considered. The network is
presented in new coordinates mean-field dynamics and synchronization errors.
Thus the problem of network synchronization is reduced to the studying of
synchronization-error system stability. The circle criterion for time-delay
systems is used to derive the stability conditions of synchronization-error
system. Obtained results are applied to a network of neural mass model
populations, and the synchronization conditions are established. Simulation
results are provided to illustrate the obtained analytical results.

###### keywords:

Synchronization, time-delay systems, circle criterion, oscillation, neural
mass model.

††thanks: This work was performed in IPME RAS and supported by Russian Science
Foundation (project no. 21-72-00107).

##  1 Introduction

Synchronization in networks of coupled oscillators is attractive phenomenon to
study for specialists in various fields of science Arenas et al. (2008).
Synchronization underlies many natural phenomena and is the cornerstone of
many technical concepts and engineering approaches. Examples of
synchronization, include, among others, numerous forms of collective behavior
in complex biological and artificial systems, such as flocks of birds,
swarming, and rendezvous Herbert-Read (2016); Sumpter (2010); ensembles of
oscillators Hong and Strogatz (2011) and a group of mobile robots Ren and
Beard (2008). Special attention is paid to synchronization in neural network
dynamics. Synchronization depends on various network parameters, and a time
delay in a signal propagation between the nodes plays a crucial role in this
phenomenon.

Time delays are always present in real physical systems, therefore, in order
to develop adequate realistic models of dynamic networks, one should to take
into account delays in signal propagation in order to properly analyze the
design of their dynamics. In neural networks, time delays can induce various
rhythmic spatiotemporal patterns Coombes and Laing (2009); Song et al. (2009),
change the stability of existing patterns Ermentrout and Ko (2009), and play a
crucial role in synchronization behavior Dahlem et al. (2009); Schnitzler et
al. (2009). Various works are devoted to study synchronization in delay
coupled networks, just to mention a few, Steur et al. (2012); Proskurnikov
(2013); Selivanov et al. (2015); Plotnikov and Fradkov (2018).

This paper continuous the work started in Plotnikov and Fradkov (2021). Here
synchronization in heterogeneous networks of linear systems with nonlinear
delayed diffusive coupling is considered. The network of neural mass model
(NMM) populations Jansen and Rit (1995) is an example of the networks of this
type. The problem of network synchronization is reduced to studying the
stability of synchronization error system which can be obtained using a
coordinate transformation proposed in Panteley and Loría (2017). The analogue
of the circle criterion for time-delay systems (TDS) proposed in Churilova
(1995) and generalized in Bryntseva and Fradkov (2019) can be used to study
the system stability.

The rest of the paper is organized as follows. Section 2 reminds some
important concepts related to synchronization and the circle criterion for
TDS. In Sec. 3 synchronization conditions of linear network with nonlinear
delayed diffusive couplings are obtained, and NMM network is considered as an
example Section 4 provides numerical results on synchronization. Finally,
conclusions are given in Sec. 5.

Notation. Throughout the paper the superscript T (∗) stands for matrix
transposition (complex conjugate); $\mathbb{R}^{n}$ denotes the $n$
dimensional real Euclidean space with vector norm $|\cdot|$; $j=\sqrt{-1}$ is
the imaginary unit; the notation $z=\operatorname{col}(x,y)$ means that $z$ is
a vector of two components $x,~{}y$; the notation
$D=\operatorname{diag}\\{d_{1},\dots,d_{n}\\}$ means that $D$ is a $n\times n$
diagonal matrix, where $d_{i}$ is its $i$th diagonal element; $A\otimes B$
means the Kronecker product of matrices $A$ and $B$; $I_{n}$ is an identity
$n\times n$ matrix, while $0_{n}$ is a $n\times n$ matrix of zeros.

##  2 Preliminaries

###  2.1 Circle Criterion for Time-Delay Systems

In this section the details about studying the stability of linear system with
multiple nonlinearities with time-varying delays are given. The analogue of
the circle criterion proposed in Churilova (1995) and its generalization for
multi input – multi output (MIMO) systems proposed in Bryntseva and Fradkov
(2019) will be used with the purpose of studying the global stability of
nonlinear delay-coupled networks.

Consider the system, which is described by the following equations

$\displaystyle\dot{x}$ $\displaystyle(t)=Ax(t)+B\varphi(t,\sigma(t-\tau(t))),$
(1) $\displaystyle\sigma$ $\displaystyle(t-\tau(t))=C^{\rm\scriptscriptstyle
T}x(t-\tau(t)),$

where $A$ is a constant $n\times n$ matrix, $B$, $C$ are constant $n\times m$
matrices and $x\in\mathbb{R}^{n}$ is a state vector,
$\sigma=\operatorname{col}\\{\sigma_{1},\dots,\sigma_{m}\\}\in\mathbb{R}^{m}$
is an input of the system,
$\tau_{i}=\operatorname{col}\\{\tau_{1},\dots,\tau_{m}\\}\in\mathbb{R}^{m}$,
$\tau_{i}(t)\in[0,T]$, $i=1,\dots,m$, $\forall t$ is a bounded time-varying
delay, $\varphi=\operatorname{col}\\{\varphi_{1},\dots,\varphi_{m}\\}$ is a
vector function of sector-bounded nonlinearities, i.e. the following
inequalities are fulfilled:

$\begin{gathered}\mu_{11}\leq\varphi_{1}(t,\sigma_{1})/\sigma_{1}\leq\mu_{21},\\\
\vdots\\\
\mu_{1m}\leq\varphi_{m}(t,\sigma_{m})/\sigma_{m}\leq\mu_{2m},\end{gathered}$
(2)

for $\sigma_{i}\neq 0$, $i=1,\dots,m$.

###### Theorem 1

For the system (1) denote the transfer function of its linear part
$W(p)=C^{\rm\scriptscriptstyle T}(A-pI_{n})^{-1}B$ and the characteristic
polynomial $\Delta(p)=\det(pI_{n}-A)$. Introduce the following diagonal
matrices $\mu_{1}=\operatorname{diag}\\{\mu_{11},\dots,\mu_{1m}\\}$,
$\mu_{2}=\operatorname{diag}\\{\mu_{21},\dots\mu_{2m}\\}$. Let the following
assumptions be fulfilled:

  1. 1.

nonlinearity $\varphi(\sigma)$ in system (1) satisfies the inequalities (2)
for $\sigma_{i}\neq 0$, $i=1,\dots,m$;

  2. 2.

There exists diagonal $m\times m$ matrix $\mu_{0}$ such that each element
$\mu_{0i}$ lies between $\mu_{1i}$ and $\mu_{2i}$, $i=1,\dots,m$, and matrix
$A+B\mu_{0}C^{\rm\scriptscriptstyle T}$ is Hurwitz;

  3. 3.

For some diagonal $m\times m$ matrix $\nu$ with positive diagonal elements
$\nu_{i}$ such that $1-4\nu_{i}\mu_{1i}\mu_{2i}>0$, $i=1,\dots,m$ the function

$\pi(\omega)=W(j\omega)^{*}(\mu_{1}\mu_{2}\nu-T^{2}\omega^{2}I_{m}/4)(I_{m}-4\nu\mu_{1}\mu_{2})\\\
\times
W(j\omega)+\operatorname{Re}[W(j\omega)(\mu_{1}+\mu_{2})(I_{m}-4\nu\mu_{1}\mu_{2})\nu]\\\
+\nu(I_{m}-\nu(\mu_{1}+\mu_{2})^{2}),$ (3)

satisfies the following conditions

$\displaystyle\lim\limits_{\omega\to\infty}\pi(\omega)$ $\displaystyle>$
$\displaystyle 0,$ (4a) $\displaystyle|\Delta(j\omega)|^{2}\pi(\omega)$
$\displaystyle>$ $\displaystyle 0,\quad\forall\omega\geq 0,$ (4b)

(for $\omega$ such that $|\Delta(j\omega)|=0$ holds the inequality (4b) is
understood as a limiting one).

Then there exist positive constants $C_{1}$, $C_{2}$, $\epsilon$ depending
only on the coefficients of the linear part of the system (1) $A$, $B$, $C$
and matrices $\mu_{1}$, $\mu_{2}$, $TI_{n}$, such that for all solutions
$x(t)$ of the system (1) with a continuous initial function $x_{0}(t)$ defined
for $t\in[-\tau_{\max},0]$, where $\tau_{\max}=\max_{i=1,\dots,m}\tau_{i}(0)$,
the following inequality holds

$\|x(t)\|\leq(C_{1}\|x_{0}(0)\|+C_{2}\max\limits_{-\tau_{\max}\leq t\leq
0}\|C^{\rm\scriptscriptstyle T}x_{0}(t)\|)\e^{-\epsilon t},$

for all $t\geq 0$.

###  2.2 Synchronization

Here the mathematical notion of synchronization will be introduced. The
networks of linear systems with heterogeneous delayed nonlinearities are
considered throughout this paper. Each system of this type can be presented in
normal form as the following:

$\displaystyle\dot{y}_{i}(t)$
$\displaystyle=A^{y}x_{i}(t)+\varphi_{i}(\sigma_{i}(t-\tau_{i}(t))),$ (5)
$\displaystyle\dot{z}_{i}(t)$ $\displaystyle=A^{z}x_{i}(t),\quad i=1,\dots,N,$

where $y_{i}\in\mathbb{R}^{m}$, $z_{i}\in\mathbb{R}^{n-m}$,
$x_{i}=\operatorname{col}\\{y_{i},z_{i}\\}\in\mathbb{R}^{n}$ and
$\sigma_{i}\in\mathbb{R}^{m}$ denote the input, the zero-dynamics, the state
and the output of the $i$th system, respectively. $A^{y}\in\mathbb{R}^{m\times
n}$ and $A^{z}\in\mathbb{R}^{(n-m)\times n}$ are constant matrices.
$\varphi_{i}=\operatorname{col}\\{\varphi_{i1}\dots\varphi_{im}\\}:\mathbb{R}^{m}\to\mathbb{R}^{m}$
is a vector function, while
$\tau_{i}=\operatorname{col}\\{\tau_{i1},\dots,\tau_{iN}\\}$ is a time-varying
delay. As input diffusive coupling is considered, which is described by

$\sigma_{i}(t-\tau_{i}(t))=\sum\limits_{k=1}^{N}\gamma_{ik}[y_{i}(t-\tau_{ik}(t))-y_{k}(t-\tau_{ik}(t))].$
(6)

Suppose that the graph of the network under consideration is connected and
undirected, therefore $\gamma_{ik}=\gamma_{ki}$ $\forall i\neq k$,
$i,k=1,\dots,N$. For this type of coupling both signals are time-delayed. Such
type of coupling may be observed, for instance, when the systems are
interconnected by a centralized control law.

Synchronization phenomenon is often defined as the asymptotically identical
evolution of the systems. One can easily introduce the notion of the
asymptotic coordinate synchronization Fradkov (2007):

$\lim\limits_{t\to\infty}(x_{i}(t)-x_{k}(t))=0,\quad i,k=1,\dots,N.$ (7)

The fulfilment of (7) means that the asymptotic behavior of all nodes of the
network (5), (6) is identical and can be described by the function
$x_{s}\in\mathbb{R}^{n}$. Thus by defining the synchronization errors as
$e_{i}=x_{i}-x_{s}$ one can study the problem of network synchronization as
the problem of stability of synchronization-error system (SES).

To obtain the equations of SES the approach proposed in Panteley and Loría
(2017) can be used. The idea behind this approach is the following: the system
state space is decomposed in two orthogonal subspaces, one on which is
projected the behavior of the mean-field state and one in which lay the
synchronization errors. To obtain the network equations in new coordinates
suppose that the graph of considered network is connected and undirected.
Following the steps described in Plotnikov and Fradkov (2021) up to a time-
delay one can present the equations of the network (5), (6) in the following
form:

$\displaystyle\dot{\tilde{x}}_{1}(t)$ $\displaystyle=$ $\displaystyle
A\tilde{x}_{1}(t)+(\mathbf{1}_{N}\otimes E_{m})$ (8a) $\displaystyle\times$
$\displaystyle\Phi[(LU_{1}\otimes E_{m}^{\rm\scriptscriptstyle
T})\tilde{x}_{2}(t-\tau(t))],$ $\displaystyle\dot{\tilde{x}}_{2}(t)$
$\displaystyle=$ $\displaystyle(I_{N-1}\otimes
A)\tilde{x}_{2}(t)+(U_{1}^{\dagger}\otimes E_{m})$ (8b) $\displaystyle\times$
$\displaystyle\Phi[(LU_{1}\otimes E_{m}^{\rm\scriptscriptstyle
T})\tilde{x}_{2}(t-\tau(t))],$

where

$\displaystyle A=\begin{bmatrix}A^{y}\\\
A^{z}\end{bmatrix}\in\mathbb{R}^{n\times
n},\quad\tau(t)=\begin{bmatrix}\tau_{1}(t)\\\ \vdots\\\
\tau_{N}(t),\end{bmatrix}\in\mathbb{R}^{N^{2}}$ (9)
$\displaystyle\Phi(\sigma)=\begin{bmatrix}\varphi_{1}(\sigma_{1})\\\ \vdots\\\
\varphi_{N}(\sigma_{N})\end{bmatrix}:\mathbb{R}^{mN}\to\mathbb{R}^{mN}$

are the matrix of the linear part of the individual system, the vector of
delays and the nonlinear vector function, respectively;

$L=\begin{bmatrix}\sum_{k=2}^{N}\gamma_{1k}&-\gamma_{12}&\cdots&-\gamma_{1N}\\\
-\gamma_{21}&\sum_{k=1,k\neq 2}^{N}\gamma_{2k}&\cdots&-\gamma_{2N}\\\
\vdots&\vdots&\ddots&\vdots\\\
-\gamma_{N1}&-\gamma_{N2}&\cdots&\sum_{k=1}^{N-1}\gamma_{Nk}\end{bmatrix}\in\mathbb{R}^{N\times
N},$ (10)

is the Laplace matrix defining the coupling links in the network, where
$\gamma_{ik}=\gamma_{ki}$ $\forall i\neq k,~{}i,k=1,\dots,N$ by the
assumption;

$\displaystyle\mathbf{1}_{N}=\begin{bmatrix}1\\\ \vdots\\\
1\end{bmatrix}\in\mathbb{R}^{N},\quad E_{m}=\begin{bmatrix}I_{m}\\\
0_{(m-n)\times n}\end{bmatrix}\in\mathbb{R}^{n\times m},$ (11) $\displaystyle
U=\begin{bmatrix}\mathbf{1}_{N}&U_{1}\end{bmatrix}\in\mathbb{R}^{N\times
N},\quad U^{-1}=\begin{bmatrix}\mathbf{1}_{N}^{\rm\scriptscriptstyle T}\\\
U^{\dagger}_{1}\end{bmatrix}\in\mathbb{R}^{N\times N}$
$\displaystyle\mbox{such
that~{}}U^{-1}LU=\operatorname{diag}\\{\lambda_{1},\dots,\lambda_{N}\\}$

are the auxiliary matrices and vectors;

$x=\begin{bmatrix}x_{1}\\\ \vdots\\\
x_{N}\end{bmatrix}\in\mathbb{R}^{nN},~{}{\tilde{x}}_{1}=\mathbf{1}_{N}^{\rm\scriptscriptstyle
T}\otimes I_{n}x,~{}{\tilde{x}}_{2}=U_{1}^{\dagger}\otimes I_{n}x$

are the state space vector of the whole network, the vector proportional to
the mean-field dynamics (MFD), and the vector proportional to the
synchronization error.

###  2.3 Neural Mass Model

As an example of linear system with nonlinear delayed couplings a time-delay
NMM will be considered. This model is based on the standard NMM proposed in
Jansen and Rit (1995) and incorporates a time delay Geng et al. (2014). The
presence of a time delay in neuronal signal transmission could cause seizure-
like activity in the brain. The NMM simulates the average firing activity of a
population of pyramidal neurons that interacts with two populations of
intercalary neurons and integrates inhibitory and excitatory signals from them
Jansen and Rit (1995); Grimbert and Faugeras (2006). The dynamics of each
neuronal population is described by two first-order nonlinear delay
differential equations as follows:

$\displaystyle\dot{y}_{1}(t)$ $\displaystyle=-2\alpha
y_{1}(t)-\alpha^{2}y_{2}(t)+\alpha\beta\varphi[\sigma(t-\tau(t))],$
$\displaystyle\dot{y}_{2}(t)$ $\displaystyle=y_{1}(t).$

where $y_{2}$ refers to the post-synaptic potential, i.e. the deviation of the
membrane from the resting potential, while $y_{1}$ is its derivative, $\sigma$
is a delayed input. The parameters $\alpha$ and $\beta$ are different in the
excitatory and inhibitory cases. $\alpha$ is the reciprocal of the
synaptic/membrane time constant; $\beta$ is the gain for the post-synaptic
response kernel; $\tau(t)$ is a time-varying delay. The function
$\varphi(\sigma)$ is a nonlinear centered sigmoidal function relating the
neuronal states

$\varphi(\sigma)=\frac{g}{1+\e^{k_{0}-\sigma}}-\frac{g}{1+\e^{k_{0}}},$ (12)

where $k_{0}$ represents the ratio of average inhibitory synaptic gain and
$g>0$ is the average excitatory synaptic gain.

##  3 Main Result

###  3.1 Synchronization Conditions. General Case

This paper considers a network of $N$ linear dynamical systems with
heterogeneous delayed nonlinear couplings in normal form (5). The graph of
considered network is supposed to be connected and undirected, and the
connections between the nodes of the network are diffusive ones (6). As
described in the Subsec. 2.2, such a network can be represented in the form of
MFD (8a) and SES (8b). To study the synchronization problem of such network
one can apply the circle criterion for TDS (see Subsec. 2.1) to the SES (8b)
to obtain the conditions of its stability. These conditions in turn will
guarantee synchronization in the original network.

The transfer function of the linear part of the system (8b) can be calculated
just like in Plotnikov and Fradkov (2021):

$W(p)=L\otimes[E_{m}^{\rm\scriptscriptstyle T}(A-pI_{n})^{-1}E_{m}],$ (13)

where $L$ is the Laplace matrix (10), $A$ is the matrix of linear part (9) and
$E_{m}$ is the supplementary matrix (11).

The nonlinear part of the system (8b) are described by the functions
$\varphi_{i}=\operatorname{col}\\{\varphi_{i1}\dots\varphi_{im}\\}$,
$i=1,\dots,N$. Suppose that they belong to the two-cavity sector between two
straight lines, i.e. the following inequalities are fulfilled:

$\begin{gathered}\mu_{1i1}\leq\varphi_{i1}(\sigma_{i1})/\sigma_{i1}\leq\mu_{2i1},\\\
\vdots\\\ \mu_{1im}\leq\varphi_{im}(\sigma_{im})/\sigma_{im}\leq\mu_{2im},\\\
i=1,\dots,N.\end{gathered}$ (14)

Introduce diagonal matrices matrices by the following way

$\begin{gathered}\mu_{1}=\operatorname{diag}\\{\mu_{111},\dots,\mu_{1m1},\dots,\mu_{11N},\dots,\mu_{1mN}\\},\\\
\mu_{2}=\operatorname{diag}\\{\mu_{211},\dots,\mu_{2m1},\dots,\mu_{21N},\dots,\mu_{2mN}\\}.\end{gathered}$

Suppose that all delays
$\tau_{i}=\operatorname{col}\\{\tau_{i1},\dots,\tau_{iN}\\}$, $i=1,\dots,N$ in
system (8b) are bounded functions, i.e. $\tau_{i}(t)\in[0,T]$, $i=1,\dots,m$,
$\forall t$.

Thus the theorem about synchronization of linear networks with nonlinear
delayed couplings can be formulated.

###### Theorem 2

If the following conditions are fulfilled

  1. 1.

The network (5), (6) nonlinearities lie in the sector, i.e. inequalities (14)
hold;

  2. 2.

The graph of the network (5), (6) is connected and undirected.

  3. 3.

There exists matrix $\mu_{0}\in\mathbb{R}^{mN}\times\mathbb{R}^{mN}$ such that
matrices $\mu_{0}-\mu_{1}$ and $\mu_{2}-\mu_{0}$ have only nonnegative
elements, and matrix

$\Psi=I_{N-1}\otimes A+[U_{1}^{\dagger}\otimes E_{m}]\mu_{0}[LU_{1}\otimes
E_{m}^{\rm\scriptscriptstyle T}]$

is Hurwitz;

  4. 4.

For some diagonal $mN\times mN$ matrix $\nu$ with positive diagonal elements
such that matrix $I_{mN}-4\nu\mu_{1}\mu_{2}$ has positive diagonal elements,
the function

$\pi(\omega)=W(j\omega)^{*}(\mu_{1}\mu_{2}\nu-T^{2}\omega^{2}I_{mN}/4)\\\
\times(I_{mN}-4\nu\mu_{1}\mu_{2})W(j\omega)\\\
+\operatorname{Re}[W(j\omega)(\mu_{1}+\mu_{2})(I_{mN}-4\nu\mu_{1}\mu_{2})\nu]\\\
+\nu(I_{mN}-\nu(\mu_{1}+\mu_{2})^{2}),$ (15)

satisfies the following conditions

$\displaystyle\lim\limits_{\omega\to\infty}\pi(\omega)$ $\displaystyle>$
$\displaystyle 0,$ (16a) $\displaystyle|\Delta(j\omega)|^{2}\pi(\omega)$
$\displaystyle>$ $\displaystyle 0,\quad\forall\omega\geq 0,$ (16b)

(for $\omega$ such that $|\Delta(j\omega)|=0$ holds the inequality (16b) is
understood as a limiting one), where $W(p)$ is the transfer function (13) and

$\Delta(p)=\det(pI_{n(N-1)}-I_{N-1}\otimes A)\\\
=\det(I_{N-1}\otimes(pI_{n}-A))=(\det(pI_{n}-A))^{N-1}$ (17)

is the characteristic polynomial of the matrix of the linear part (8b).

Then the systems in the network (5), (6) are asymptotically synchronized.

###  3.2 Synchronization Conditions of Neural Mass Model Populations

This section considers the heterogeneous network of nonlinearly delayed
coupled NMM populations

$\displaystyle\dot{y}_{1i}(t)$ $\displaystyle=-2\alpha
y_{1i}(t)-\alpha^{2}y_{2i}(t)$ (18)
$\displaystyle+\alpha\beta\varphi_{i}\left\\{\sum\limits_{k=1}^{N}\gamma_{ik}[y_{1i}(t-\tau_{ik}(t))-y_{1k}(t-\tau_{ik}(t))]\right\\},$
$\displaystyle\dot{y}_{2i}(t)$ $\displaystyle=y_{1i}(t),\quad i=1,\dots,N,$

where $x=\operatorname{col}\\{y_{11},y_{21},\dots,y_{1N},y_{2N}\\}$ is a state
vector; $\alpha$, $\beta$ are system parameters; $\gamma_{ik}$ are coupling
coefficients; $\tau_{ik}(t)\in[0,T]$, $\forall t$, $i,k=1,\dots,N$ are bounded
time-varying delays (all delay functions have the same upper bound $T$).
Functions $\varphi_{i}$ are sigmoidal ones, which are descirbed by (12) with
parameters $g_{i}$ and $k_{0i}$, $i=1,\dots,N$.

The network of NMMs can be presented in new coordinates (8a), (8b) with
matrices

$A=\begin{bmatrix}-2\alpha&-\alpha^{2}\\\ 1&0\end{bmatrix},\quad
E_{m}=\begin{bmatrix}1\\\ 0\end{bmatrix}.$

To study the network (18) synchronization one should check the conditions of
Theorem 2:

Sigmoidal functions (12) are sector ones, which lie between two straight lines
$0$ and $0.5g_{i}\sigma$ (see explanation in Gorshkov et al. (2017)),
$i=1,\dots,N$. One can find $g_{\max}=\max_{i=1,\dots,N}g_{i}$. Suppose that
$\alpha\beta>0$, then the matrices from the condition (1) of the Theorem 2 can
be expressed as $\mu_{1}=0_{N}$ and $\mu_{2}=0.25g_{\max}\alpha\beta I_{N}$.

Supposing that the graph of considered network (18) is connected and
undirected guarantees the fulfillment of the condition (2) of the Theorem 2.

The condition (3) of Theorem 2 is the same as for the case without delays:
this fact follows from the conditions of Circle criterion. This condition was
previously checked in Plotnikov and Fradkov (2021): If $\alpha>0$, then the
matrix $\Psi$ in the condition (3) of the Theorem 2 is Hurwitz.

To calculate the function $\pi(\omega)$ from the condition (4) of the Theorem
2 the frequency transfer function $W(j\omega)$ should be found using formula
(13):

$W(j\omega)=L\otimes\left[\begin{bmatrix}1&0\end{bmatrix}\begin{bmatrix}-2\alpha-j\omega&-\alpha^{2}\\\
1&-j\omega\end{bmatrix}^{-1}\begin{bmatrix}1\\\ 0\end{bmatrix}\right]\\\
=\frac{-j\omega}{\alpha^{2}-\omega^{2}+2j\alpha\omega}L=\frac{-2\alpha\omega^{2}-j\omega(\alpha^{2}-\omega^{2})}{(\alpha^{2}+\omega^{2})^{2}}L.$
(19)

Meaning that $\mu_{1}=0_{N}$ and $\mu_{2}=0.25g_{\max}\alpha\beta I_{N}$ and
choosing matrix $\nu$ as $\nu_{0}I_{N}$ one obtains:

$\pi(\omega)=\frac{-T^{2}\omega^{2}}{4}W(j\omega)^{*}W(j\omega)+\frac{g_{\max}\alpha\beta\nu_{0}}{4}\operatorname{Re}[W(j\omega)]\\\
+\nu_{0}\left(1-\frac{\nu_{0}g_{\max}^{2}\alpha^{2}\beta^{2}}{16}\right)I_{N}.$
(20)

Since the graph of considered network is undirected the corresponding Laplace
matrix $L$ is symmetric. Using this fact and (19) the function $\pi(\omega)$
equals:

$\pi(\omega)=\frac{-T^{2}\omega^{4}}{4(\alpha^{2}+\omega^{2})^{2}}L^{2}-\frac{g_{\max}\alpha^{2}\beta\nu_{0}\omega^{2}}{2(\alpha^{2}+\omega^{2})^{2}}L\\\
+\nu_{0}\left(1-\frac{\nu_{0}g_{\max}^{2}\alpha^{2}\beta^{2}}{16}\right)I_{N}.$
(21)

Now check the matrix inequality (16a):

$\lim\limits_{\omega\to\infty}\pi(\omega)=\frac{-T^{2}}{4}L^{2}+\nu_{0}\left(1-\frac{\nu_{0}g_{\max}^{2}\alpha^{2}\beta^{2}}{16}\right)I_{N}>0$

Consider some symmetric matrix $P$: it is positive definite if the
corresponding quadratic form $x^{\rm\scriptscriptstyle T}Px$ is positive
$\forall x\neq 0$. $\exists$ matrix $U$: $P=UDU^{\rm\scriptscriptstyle T}$,
where $D$ is a diagonal matrix. Then for $z=U^{\rm\scriptscriptstyle T}x$ one
obtains $z^{\rm\scriptscriptstyle T}Dz>0$. Therefore the obtained inequality
is fulfilled if and only if the following inequalities are fulfilled:

$\frac{-T^{2}\lambda_{i}^{2}}{4}+\nu_{0}\left(1-\frac{\nu_{0}g_{\max}^{2}\alpha^{2}\beta^{2}}{16}\right)>0,~{}i=1,\dots,N,$
(22)

where $\lambda_{i}$ are the eigenvalues of Laplace matrix $L$. Then the
maximal value of the delay $T$ can be estimated:

$T^{2}<\frac{16\nu_{0}-\nu_{0}^{2}g_{\max}^{2}\alpha^{2}\beta^{2}}{4\lambda_{i}^{2}},\quad
i=1,\dots,N.$ (23)

There is a quadratic equation depending on $\nu_{0}$ in the numerator of
resulting fraction, which has the maximal value for
$\nu_{0}=8/(g_{\max}^{2}\alpha^{2}\beta^{2})$. All equations (23) hold if

$T<\frac{4}{g_{\max}\lambda_{\max}\alpha\beta},$ (24)

where $\lambda_{\max}$ is the maximal eigenvalue of Laplace matrix $L$.

From (16b) one obtains that inequality

$\frac{-T^{2}\omega^{4}}{4}L^{2}-\frac{g_{\max}\alpha^{2}\beta\nu_{0}\omega^{2}}{2}L\\\
+\nu_{0}\left(1-\frac{\nu_{0}g_{\max}^{2}\alpha^{2}\beta^{2}}{16}\right)(\alpha^{2}+\omega^{2})^{2}I_{N}>0$
(25)

should be fulfilled $\forall\omega\geq 0$. As before the set of the following
inequalities can be consider instead of the obtained matrix inequality:

$\left[\nu_{0}\left(1-\frac{\nu_{0}g_{\max}^{2}\alpha^{2}\beta^{2}}{16}\right)-\frac{T^{2}\lambda_{i}^{2}}{4}\right]\omega^{4}\\\
+\alpha^{2}\nu_{0}\left[2\left(1-\frac{\nu_{0}g_{\max}^{2}\alpha^{2}\beta^{2}}{16}\right)-\frac{g_{\max}\beta\lambda_{i}}{2}\right]\omega^{2}\\\
+\alpha^{4}\nu_{0}\left(1-\frac{\nu_{0}g_{\max}^{2}\alpha^{2}\beta^{2}}{16}\right)>0.$
(26)

The coefficient before $\omega^{4}$ is the same as (22), and it is positive if
$\nu_{0}=8/(g_{\max}^{2}\alpha^{2}\beta^{2})$ and the inequality (24) is
fulfilled. Zero order term is also positive for chosen value of $\nu_{0}$. The
coefficient before $\omega^{2}$ is positive if
$g_{\max}\beta\lambda_{\max}<2$. Thus the inequality (26) is fulfilled for
chosen parameters.

All conditions of the Theorem 2 are fulfilled, therefore the network (18) is
synchronized. The following theorem holds.

###### Theorem 3

If the network (18) systems parameters $\alpha>0$ and $\beta>0$, $g_{i}>0$,
$i=1,\dots,N$, the graph of the network is connected and undirected, the
maximum eigenvalue of the Laplace matrix $L$ is less than $2/(\beta
g_{\max})$, and the delays in the signal propagation are bounded (24), then
the network of NMMs is asymptotically synchronized.

Note that while $T\to 0$ and choosing $\nu_{0}\to 0$ such that
$T^{2}/\nu_{0}\to 0$, we get the similar conditions of network synchronization
as in Plotnikov and Fradkov (2021).

##  4 Simulation

This section presents the results of simulation. The network of NMM
populations with $N=10$ is considered. The system parameters $\alpha$ and
$\beta$ are equal to $1$ and $0.8$, respectively. The parameters defining the
shape of sigmoidal function have uniform distribution: $g_{i}$, $i=1,\dots,N$
are distributed on the interval $[0;1]$, while $k_{0i}$, $i=1,\dots,N$ are
distributed on the interval $[-1;1]$. The graph of considered network is
weighted, connected and undirected, meaning that its adjacency matrix is a
symmetric sparse matrix with density $0.7$, which means that it has
approximately $0.7N^{2}$ nonzero entries. Let the delays
$\tau_{ik}(t)\in[0;T]$, $i,k=1,\dots,N$ be time-varying functions

$\tau_{ik}(t)=h_{1ik}+h_{2ik}\sin(h_{3ik}t+h_{4ik}),$

which are uniformly distributed on the interval $[0;T/2]$ such that
$h_{1ik}>h_{2ik}$, $i,k=1,\dots,N$. The maximum value of the delay $T$ will be
defined later. The initial functions $y_{1i}(t)$, $y_{2i}(t)$, $t\in[-T;0]$,
$i=1,\dots,N$ are constants, which are uniformly distributed on the interval
$[-1;1]$.

Figure 1: Synchronization of neural mass model (NMM) population network (18)
with $N=10$ nodes. (a) and (b) dynamics of post-synaptic potential derivatives
$y_{1}$ and post-synaptic potentials $y_{2}$ of all nodes, respectively.
System parameters: $N=10$, $\alpha=1$, $\beta=0.8$, $T=2$, $g_{i}$,
$i=1,\dots,N$ are uniformly distributed on the interval $[0;1]$, $k_{0i}$,
$i=1,\dots,N$ are uniformly distributed on the interval $[-1;1]$,
$\lambda_{\max}=2.5080$. Initial functions $y_{1i}(t)$, $y_{2i}(t)$,
$t\in[-T;0]$, $i=1,\dots,N$ are constants, which are uniformly distributed on
the interval $[-1;1]$.

First of all consider the case, when NMM network has parameters satisfying the
conditions of the Theorem 3. The maximum eigenvalue $\lambda_{\max}$ of the
Laplace matrix $L$ is equal to $2.5080$ is this case, while $g_{\max}=0.9724$.
For these parameters of the network the inequality $\lambda_{\max}<2/(\beta
g_{\max})\approx 2.5709$ is fulfilled. Choosing $T=2$ one can ensure the
fulfillment of the inequality (24), thereby guarantee the fulfillment of all
condition of the Theorem 3. This means that for these parameters the network
of NMMs will synchronize. Figure 1 presents the results of simulation. As one
can see, for the chosen parameter values, there is synchronization among the
state variables of the network, and all system trajectories tend to
equilibrium point, which confirms Theorem 3.

Figure 2: Desynchronization of neural mass model (NMM) population network (18)
with $N=10$ nodes. (a) and (b) dynamics of post-synaptic potential derivatives
$y_{1}$ and post-synaptic potentials $y_{2}$ of all nodes, respectively.
System parameters: $T=10$, $\lambda_{\max}=1.7983$. Other parameters and
initial functions are the same as in Fig. 1.

Now consider the case, when the delays in signal propagation between the nodes
are too large, that prevents the network synchronization. In this case
$\lambda_{\max}=1.7983$ and $g_{\max}=0.9630$, which means that the inequality
$\lambda_{\max}<2/(\beta g_{\max})\approx 2.5960$ is also fulfilled. Choosing
$T=10$ one violates the Theorem 3 condition (24). Therefore for these
parameters of the network, the Theorem 3 doesn’t guarantee the network
synchronization. One can see the results of simulation in Fig 2: there is no
synchronization among the state variables of the network.

##  5 Conclusion

In this paper, the problem of heterogeneous network synchronization of linear
systems with delayed nonlinear diffusive couplings. The delays are supposed to
be bounded time-varying functions. The nonlinear coupling functions can be
different but their graphs should lie in two-cavity sector between two
straight lines. As in Plotnikov and Fradkov (2021) the synchronization problem
is reduced to study the stability of the SES. The coordinate transformation
approach proposed in Panteley and Loría (2017) is used to present the network
in coordinates ”MFD - SES”. To find the conditions for SES stability the
circle criterion for TDS was applied. The theorem about network
synchronization of this type was formulated and proven.

As an example, the dynamics of NMM populations connected via delayed nonlinear
diffusive coupling was considered. Using the obtained theorem, the simple
condition for network synchronization was derived. In the case of the delay
absence this condition coincide with the result obtained in Plotnikov and
Fradkov (2021). Also, the simulation of NMM network dynamics was performed. In
the case, when the Theorem 3 conditions are fulfilled, one can observe the
synchronization between the network states. In the other case, when the delays
are large enough, there is no synchronization between the network nodes.

## References

  * Arenas et al. (2008) Arenas, A., Díaz-Guilera, A., Kurths, J., Moreno, Y., and Zhou, C. (2008).  Synchronization in complex networks.  _Phys. Rep._ , 469(3), 93–153.  10.1016/j.physrep.2008.09.002.
  * Bryntseva and Fradkov (2019) Bryntseva, T.A. and Fradkov, A.L. (2019).  Frequency-domain estimates of the sampling interval in multirate nonlinear systems by time-delay approach.  _Intern. J. Control_ , 92, 1985–1992.  10.1080/00207179.2017.1423394.
  * Churilova (1995) Churilova, M.Y. (1995).  Analog of the cyclic criterion of absolute stability for systems with variable delays.  _Autom. Remote Control_ , 56, 195–198.
  * Coombes and Laing (2009) Coombes, S. and Laing, C. (2009).  Delays in activity-based neural networks.  _Phil. Trans. R. Soc. A_ , 367(1891), 1117–1129.  10.1098/rsta.2008.0256.
  * Dahlem et al. (2009) Dahlem, M.A., Hiller, G., Panchuk, A., and Schöll, E. (2009).  Dynamics of delay-coupled excitable neural systems.  _Int. J. Bifur. Chaos_ , 19(2), 745–753.  10.1142/S0218127409023111.
  * Ermentrout and Ko (2009) Ermentrout, B. and Ko, T.W. (2009).  Delays and weakly coupled neuronal oscillators.  _Phil. Trans. R. Soc. A_ , 367(1891), 1097–1115.  10.1098/rsta.2008.0259.
  * Fradkov (2007) Fradkov, A. (2007).  _Cybernetical physics: From control of chaos to quantum control_.  Springer-Verlag, Berlin Heidelberg.  10.1007/978-3-540-46277-4.
  * Geng et al. (2014) Geng, S., Zhou, W., Zhao, X., Yuan, Q., Ma, Z., and Wang, J. (2014).  Bifurcation and oscillation in a time-delay neural mass model.  _Biol. Cybern._ , 108(6), 747–756.  10.1007/s00422-014-0616-4.
  * Gorshkov et al. (2017) Gorshkov, A.A., Plotnikov, S.A., and Fradkov, A. (2017).  Bifurcation and synchronization analysis of neural mass model subpopulations.  _IFAC-PapersOnLine_ , 50(1), 14741–14745.  10.1016/j.ifacol.2017.08.2577.
  * Grimbert and Faugeras (2006) Grimbert, F. and Faugeras, O. (2006).  Bifurcation analysis of Jansen’s neural mass model.  _Neur. Comput._ , 18(12), 3052–3068.  10.1162/neco.2006.18.12.3052.
  * Herbert-Read (2016) Herbert-Read, J.E. (2016).  Understanding how animal groups achieve coordinated movement.  _J. Exp. Biol._ , 219(19), 2971––2983.  10.1242/jeb.129411.
  * Hong and Strogatz (2011) Hong, H. and Strogatz, S.H. (2011).  Kuramoto model of coupled oscillators with positive and negative coupling parameters: An example of conformist and contrarian oscillators.  _Phys. Rev. Lett._ , 106, 054102.  10.1103/PhysRevLett.106.054102.
  * Jansen and Rit (1995) Jansen, B. and Rit, V. (1995).  Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns.  _Biol. Cybern._ , 73(4), 357–366.  10.1007/BF00199471.
  * Panteley and Loría (2017) Panteley, E. and Loría, A. (2017).  Synchronization and dynamic consensus of heterogeneous networked systems.  _IEEE Trans. Automat. Control_ , 62(8), 3758–3773.  10.1109/TAC.2017.2649382.
  * Plotnikov and Fradkov (2018) Plotnikov, S.A. and Fradkov, A.L. (2018).  On synchronization in fitzhugh-nagumo networks with small delays.  In _2018 European Control Conference (ECC)_ , 2052–2056.  10.23919/ECC.2018.8550552.
  * Plotnikov and Fradkov (2021) Plotnikov, S.A. and Fradkov, A.L. (2021).  Synchronization of nonlinearly coupled networks based on circle criterion.  _Chaos_ , 31(10), 103110.  10.1063/5.0055814.
  * Proskurnikov (2013) Proskurnikov, A.V. (2013).  Average consensus in networks with nonlinearly delayed couplings and switching topology.  _Automatica_ , 49(9), 2928–2932.  10.1016/j.automatica.2013.06.007.
  * Ren and Beard (2008) Ren, W. and Beard, R.W. (2008).  _Distributed consensus in multi-vehicle cooperative control: theory and applications_.  Springer-Verlag, London.  10.1007/978-1-84800-015-5.
  * Schnitzler et al. (2009) Schnitzler, A., Munks, C., Butz, M., Timmermann, L., and Gross, J. (2009).  Synchronized brain network associated with essential tremor as revealed by magnetoencephalography.  _Mov. Disorders_ , 24(11), 1629–1635.  10.1002/mds.22633.
  * Selivanov et al. (2015) Selivanov, A., Fradkov, A., and Fridman, E. (2015).  Passification-based decentralized adaptive synchronization of dynamical networks with time-varying delays.  _J. Franklin Inst._ , 352(1), 52–72.  10.1016/j.jfranklin.2014.10.007.
  * Song et al. (2009) Song, Y., Makarov, V.A., and Velarde, M.G. (2009).  Stability switches, oscillatory multistability, and spatio-temporal patterns of nonlinear oscillations in recurrently delay coupled neural networks.  _Biol. Cybern._ , 101(2), 147–167.  10.1007/s00422-009-0326-5.
  * Steur et al. (2012) Steur, E., Oguchi, T., van Leeuwen, C., and Nijmeijer, H. (2012).  Partial synchronization in diffusively time-delay coupled oscillator networks.  _Chaos_ , 22, 043144.  10.1063/1.4771665.
  * Sumpter (2010) Sumpter, D.J. (2010).  _Collective Animal Behavior_.  Princeton University Press, Princeton.