# Anti-phase synchronization in a population of swarmalators

Samali Ghosh<EMAIL_ADDRESS>Physics and Applied Mathematics Unit,
Indian Statistical Institute, 203 B. T. Road, Kolkata 700108, India  Gourab
Kumar Sar<EMAIL_ADDRESS>Physics and Applied Mathematics Unit, Indian
Statistical Institute, 203 B. T. Road, Kolkata 700108, India  Soumen Majhi
<EMAIL_ADDRESS>Department of Ecology and Evolution, University of
Chicago, Chicago, Illinois 60637, USA  Dibakar Ghosh<EMAIL_ADDRESS>Physics and Applied Mathematics Unit, Indian Statistical Institute, 203 B. T.
Road, Kolkata 700108, India

###### Abstract

(Received XX MONTH XX; accepted XX MONTH XX; published XX MONTH XX)
Swarmalators are oscillatory systems endowed with a spatial component, whose
spatial and phase dynamics affect each other. Such systems can demonstrate
fascinating collective dynamics resembling many real-world processes. Through
this work, we study a population of swarmalators where they are divided into
different communities. The strengths of spatial attraction, repulsion as well
as phase interaction differ from one group to another. Also, they vary from
inter-community to intra-community. We encounter, as a result of variation in
the phase coupling strength, different routes to achieve the static
synchronization state by choosing several parameter combinations. We observe
that when the inter-community phase coupling strength is sufficiently large,
swarmalators settle in the static synchronization state. On the other hand,
with a significant small phase coupling strength the state of anti-phase
synchronization as well as chimera-like coexistence of sync and async are
realized. Apart from rigorous numerical results, we have been successful to
provide semi-analytical treatment for the existence and stability of global
static sync and the anti-phase sync states.

DOI: XXXXXXX

##  I Introduction

Synchronization Pecora and Carroll (1990); Boccaletti et al. (2002); Arenas et
al. (2008) refers to the phenomenon in which interacting dynamical systems
adjust their rhythm in time, and the study of diverse aspects of this process
has been at the forefront of research in nonlinear dynamics of networked
systems. This is mainly due to its applicability in social and physical
methods to biological and technological systems Rosenblum and Kurths (2003);
Arenas et al. (2008). From opinion formation Jusup et al. (2022), flashing of
ﬁreﬂies Buck (1988) to Josephson junction Wiesenfeld et al. (1996) or firing
neurons Montbrió et al. (2015), synchronization takes place in a variety of
natural and man-made systems. Alternative to the synchrony in the states of
the systems, self-organization in space takes place in flocking birds, school
of fish, swarming insects, a herd of sheep Buhl et al. (2006); Couzin (2007);
Ballerini et al. (2008); Bialek et al. (2012); Leonard et al. (2012);
Garcimartín et al. (2015) and even in micro-organisms Zhang et al. (2010);
Levy and Requeijo (2008); Tian et al. (2021). This phenomenon of spatial self-
organization without explicit alteration in the internal states is referred to
as the swarming Fetecau et al. (2011); Darnton et al. (2010); Mogilner and
Edelstein-Keshet (1999). On the contrary, as mentioned earlier, the internal
state dynamics play the primary role in synchronization in which spatial
movement does not necessarily participate. In the last few decades, these two
fields, synchrony, and swarming, have been studied independently and parallel
to each other. The study of mobile oscillators or moving agents brought these
two fields into contact by considering the effect of oscillators’ motion on
the internal dynamics Frasca et al. (2008); Majhi et al. (2019a). In nature,
there are many cases of this behavior such as vinegar eels Peshkov et al.
(2022), Japanese tree frogs Aihara et al. (2014), starfish embryos Tan et al.
(2022) where the spatial and internal dynamics are dependent on each other.
Recently, the combined effect of sync and swarming is expressed by a coupled
system, namely swarmalators. The swarmalators O’Keeffe et al. (2017); Sar and
Ghosh (2022); Yoon et al. (2022) are systems with concurrent existence of
synchronization and swarming dynamics, particularly oscillatory systems having
spatial and phase dynamics coupled. The study of such a fascinating interplay
between the internal states and the positions in space was initiated by Tanaka
et al. Tanaka (2007); Iwasa et al. (2010) who presented a swarm-oscillator
model, followed by the recent swarmalator model proposed by O’Keeffe et al.
O’Keeffe et al. (2017). The latter model built upon a space-dependent
generalized version of globally coupled Kuramoto oscillators predicts few
novel collective dynamics.

Significant efforts have been made since then in the last few years in order
to understand the dynamics of swarmalators under different system-interaction
setups. Phase similarity arising through both spatial attraction and spatial
repulsion can result in ring phase wave states O’Keeffe et al. (2018).
External periodic forcing affecting the phases leads to phase transitions from
the states of the non-forced model through partial to full synchrony Lizarraga
and de Aguiar (2020). Besides the impact of thermal noise on the swarmalator
system Hong et al. (2023), a number of other aspects, such as distributed
coupling O’Keeffe and Hong (2022), time-delayed interactions Blum et al.
(2022), finite-cutoff interaction distance Lee et al. (2021) are examined. The
outcomes are a plethora of collective patterns, some observable in as
disparate as Japanese tree frogs and electroporated Quincke rollers. Sar et
al. Sar et al. (2022) has come up with a swarmalator model subject to time-
varying competitive phase interaction in which the competition between the
attractive and repulsive interaction takes place depending on the sensing
radius of the units. Swarmalators on a ring subject to random pinning are
investigated Sar et al. (2023) and found to result in low-dimensional chaos,
an abrupt transition to synchronous state, along with phase wave and split
phase wave. Lately, Ceron et al. Ceron et al. (2023) demonstrates that the
edition of non-identical frequencies of the oscillators, local coupling, and
chirality lead to new dynamics including beating clusters and lattices of
vortices. Nevertheless, research in this fascinating world of swarmalators is
still in its infancy and there are adequate scopes of further investigation
leading to the possible revelation of new emerging collective dynamics due to
the bidirectional reciprocity between the phase and the spatial dynamics.

One of the most pivotal characteristics of many real-world networked systems
is that of community structures or clustering Girvan and Newman (2002);
Fortunato (2010) referring to the compartmental subdivisions of networked
systems. This, precisely, corresponds to the organization of the units of the
system in strongly intra-connected communities or groups while possessing
weaker inter-group connections. From numerous social systems including
collaboration networks to biological networks, such as metabolic networks,
regulatory networks, and food webs, are naturally found to exhibit community
structures Girvan and Newman (2002); Fortunato (2010); Newman (2006). The
problem of detection and characterization of these communities Fortunato
(2010); Malliaros and Vazirgiannis (2013); Clauset et al. (2004); Newman and
Girvan (2004) is one of the preeminent issues in the study of structural
network theory. Through this article, we assume a community-structured
framework of the underlying network and demonstrate the genesis of multiple
variants of collective patterns in interacting communities of swarmalators.
We, specifically, study a population of swarmalators where they are
distributed in two communities. We analyze how the trade-off between the
intra- and inter-community interactions affects the generic interplay between
the phase and spatial dynamics of swarmalators. The phase interactions along
with the spatial attraction and repulsion differ in each community. Under such
a network setup, we encounter diverse routes to the static synchronization
state as the inter-community phase coupling strength increases, for different
choices of parameter values. Besides the states like active and static async
or active phase wave, we detect anti-phase synchrony and the chimera state in
the process towards the emergence of in-phase synchronization. We must here
emphasize the remarkable fact that the anti-phase synchronization state arises
in the sole presence of repulsive coupling even when the network size is
considerably large, which we do not experience in the case of simple phase
oscillator models without any spatial dynamics Tsimring et al. (2005). We have
also provided semi-analytical treatment concerning the stability analysis of
the global static synchronization and the anti-phase synchronization state. We
should mention here that, in a stereotypical community network structure, the
strength of inter-community interaction is usually considered to be smaller
than that of the intra-community interaction strength Girvan and Newman
(2002). In this work, however, we have not necessarily followed this
convention. We have varied the inter-community phase coupling strength over a
feasible range where a number of diverse collective states are observed.

##  II Proposed mathematical model

We consider $N$ number of swarmalators moving in a two-dimensional region. We
randomly distribute them in $p$ groups. Let $C_{i}$ denote the set of indices
of swarmalators belonging to the $i$-th group. Then trivially we have,
$\sum_{i=1}^{p}|C_{i}|=N$, where $|C_{i}|$ denotes the cardinality of the set
$C_{i}$ and $\cup_{i=1}^{p}C_{i}=\\{1,2,\cdots,N\\}$. Suppose, without loss of
generality, that the $i$-th swarmalator belongs to the $n$-th group. Then we
can write the governing equations as,

$\displaystyle\dot{\mathbf{x}}_{i}=\mathbf{v}_{i}$
$\displaystyle+\sum_{m=1}^{p}\dfrac{1}{|C_{m}|}\sum_{j\in
C_{m}\setminus\\{i\\}}\Bigg{[}\dfrac{\mathbf{x}_{j}-\mathbf{x}_{i}}{|\mathbf{x}_{j}-\mathbf{x}_{i}|}\big{(}1+$
$\displaystyle
J_{n,m}\cos(\theta_{j}-\theta_{i})\big{)}-\dfrac{\mathbf{x}_{j}-\mathbf{x}_{i}}{|\mathbf{x}_{j}-\mathbf{x}_{i}|^{2}}\Bigg{]},$
(1)
$\displaystyle\dot{\theta}_{i}=\omega_{i}+\sum_{m=1}^{p}\dfrac{K_{n,m}}{|C_{m}|}\sum_{j\in
C_{m}\setminus\\{i\\}}\dfrac{\sin(\theta_{j}-\theta_{i})}{|\mathbf{x}_{j}-\mathbf{x}_{i}|},$
(2)

where $i=1,2,\cdots,N$. $\textbf{x}_{i}\equiv(x_{i},y_{i})$ is the spatial
position in two-dimensional plane and $\theta_{i}$ is the internal phase of
the $i$-th swarmalator. $\omega_{i}$ and $\mathbf{v}_{i}$ are the natural
frequency, and self-propulsion velocity of the $i$-th swarmalator,
respectively. The spatial attraction, repulsion as well as phase interaction
functions are chosen the same as in Ref. O’Keeffe et al., 2017 where all the
swarmalators belong to a single group, i.e., $p=1$. The spatial attraction
term ensures that the swarmalators remain close to each other without
dispersing indefinitely, whereas spatial repulsion among them is necessary to
avoid collision. They can be perceived as long-range attraction and short-
range repulsion. $J_{n,m}$ highlights how phases of those two swarmalators
affect their spatial attraction. We assume $J_{n,m}>0$ so that swarmalators
which are in nearby phases attract each other spatially due to the presence of
the term $\cos(\theta_{j}-\theta_{i})$. Similarly, $K_{n,m}$ indicates the
phase coupling strength between the two groups $C_{n}$ and $C_{m}$ (note that,
here by group $C_{n}$, we mean the swarmalators belonging to the $n$-th group,
without ambiguity). When $K_{n,m}>0$, swarmalators’ phases are attractively
coupled and the phase coupling is repulsive when $K_{n,m}<0$. For symmetry,
$J_{n,m}=J_{m,n}$, and $K_{n,m}=K_{m,n}$. Then, for $p$ groups, the number of
distinct parameters related to $J$ and $K$ are $(p^{2}+p)/2$ each. We work
with $p=2$ groups in this article which leaves us with $J_{1,1}\equiv J_{1}$,
$J_{2,2}\equiv J_{2}$, $J_{1,2}=J_{2,1}\equiv J_{3}$, $K_{1,1}\equiv K_{1}$,
$K_{2,2}\equiv K_{2}$, and $K_{1,2}=K_{2,1}\equiv K_{3}$, say. Effectively we
have these six parameters in hand which we vary to obtain different collective
behaviors. Also note that, the model defined by Eqs. (1)-(2) is a
generalization of the model proposed by O’Keeffe et al. O’Keeffe et al.
(2017). We work with swarmalators having identical natural frequencies and
velocities, i.e., $\omega_{i}=\omega$ and $\mathbf{v}_{i}=\mathbf{v}$ for all
$i$. By moving to a proper reference frame, we set $\omega=|\mathbf{v}|=0$.

##  III Results

First, we assume that the swarmalators are distributed in equal numbers in two
populations (we remove this assumption later in Appendix A to show that the
results do not change if they are distributed unequally as long as $N$ is
sufficiently large). For simplicity, let $C_{i}$ denote both the $i$-th
population and the set of indices of swarmalators belonging to that
population, whenever appropriate, for $i=1,2$. Then $J_{1}$ measures the
extent to which the phases of swarmalators belonging to $C_{1}$ affect their
spatial attraction and similarly $J_{2}$ for the group $C_{2}$. $J_{3}$ gauges
the phase-dependent spatial attraction when swarmalators belong to different
groups. $K_{1}$, and $K_{2}$ are the phase coupling strengths between
swarmalators in $C_{1}$ and $C_{2}$, respectively, whereas $K_{3}$ is the
strength of phase interaction when they are in different groups. The values of
these control parameters decide the fate of the swarmalator system where we
observe various emerging collective states by changing these values. Before
moving forward to describe these states, first, we define some order
parameters that are useful to measure several properties of the emerging
states.

###  III.1 Order parameters

To measure the amount of synchrony in swarmalators’ phases throughout the
population, we define

$re^{i\psi}=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}}.$ (3)

Here $r$ lies between $0$ to $1$ by definition and gives an indication of the
overall synchrony in swarmalators’ phases. Phases are completely synchronized
when $r=1$, or else asynchronous behavior is present. $\psi$ is the mean phase
of the overall population. We measure the phase coherence among swarmalators
belonging to the $p$-th group by

$r_{p}e^{i\psi_{p}}=\frac{1}{|C_{p}|}\sum_{j\in C_{p}}e^{i\theta_{j}},$ (4)

where $r_{p}$ again lies between $0$ to $1$ and $\psi_{p}$ is the average
phase of $p$-th group. We also define

$Re^{i\Psi}=\frac{1}{N}\sum_{j=1}^{N}e^{2i\theta_{j}},$ (5)

which is useful to examine anti-phase synchrony where a phase difference of
$\pi$ is observed among swarmalators’ phases. In the anti-phase synchrony
state, $R=1$ but $r\neq 1$. In some of the collective states (discussed later
in Sec. III.2) we observe a correlation between swarmalators’ phases
$\theta_{j}$ and their spatial angle $\phi_{j}=\tan^{-1}(y_{j}/x_{j})$. For
this, we define the following order parameters,

$S_{\pm}e^{i\Psi\pm}=\frac{1}{N}\sum_{j=1}^{N}e^{i(\phi_{j}\pm\theta_{j})},$
(6)

which quantifies the correlation between phases and spatial angles. We take
the maximum of $S_{\pm}$ and define $S=\max\\{S_{+},S_{-}\\}$. A nonzero value
of $S$ indicates the presence of a correlation between swarmalators’ spatial
angles and phases. In one of the collective states, swarmalators arrange
themselves inside an annular-like structure and they rotate around this
annulus. Their phases also keep changing from $0$ to $2\pi$. To distinguish
this state from others, $\gamma$ is defined as

$\gamma=\frac{N_{rot}}{N},$ (7)

where $N_{rot}$ is the number of swarmalators executing at least one full
circle rotation in both spacial location and phase. $\gamma$ gives the
fraction of such swarmalators and subsequently lies between $0$ to $1$. We
find both stationary and non-stationary states in our model where swarmalators
become static in position and phase in the stationary states but keep moving
in the non-stationary ones. To separate these, we measure the mean velocity
denoted by $V$, and is defined as

$V=\Bigg{\langle}\frac{1}{N}\sum_{i=1}^{N}\sqrt{\dot{x}_{i}^{2}+\dot{y}_{i}^{2}+\dot{\theta}_{i}^{2}}\Bigg{\rangle}_{t},$
(8)

where $\langle\cdots\rangle_{t}$ stands for the time average, which is taken
after discarding the initial transients. With the knowledge of these order
parameters, we proceed to study the emerging collective states of our model.

Figure 1: Order parameters along with the snapshots of the emerging states. (a) Variation of different order parameters with $K_{3}$. (b) Anti-phase sync for $K_{3}=-0.6$, (c) chimera for $K_{3}=-0.4$, (d) active async for $K_{3}=-0.25$ & static async for $K_{3}=0.0$, (e) active phase wave for $K_{3}=0.1$ and (f) static sync for $K_{3}=0.3$. Simulations are performed for $N=200$ swarmalators for $T=5000$ time units and step-size $dt=0.01$ using Heun’s method. In all cases, swarmalators are initially placed inside the box $[-1,1]\times[-1,1]$ uniformly at random, while their phases are drawn randomly from $[0,2\pi]$. We fix $J_{1}=J_{2}=J_{3}=0.1$ and $K_{1}=-0.1,K_{2}=-0.2$. Note that there is a long transient until the states are achieved. The order parameters are calculated with the last $10\%$ data after discarding the initial transients. Table 1: This table shows how the emerging states of the population of swarmalators are identified with the order parameters $r$, $R$, $S$, $\gamma$, and $V$.  $r$ | $R$ | $S$ | $\gamma$ | $V$ | Emerging state
---|---|---|---|---|---
$\approx 1$ | $\approx 1$ | $0<S\ll 1$ | $\approx 0$ | $\approx 0$ | Static sync
$\approx 0$ | $0<R\ll 1$ | $\approx 0$ | $\approx 0$ | $\approx 0$ | Static async
$\approx 0$ | $0<R\ll 1$ | $0<S\ll 1$ | $\ll 1$ | $\neq 0$ | Active async
$\approx 0$ | $0<R\ll 1$ | $0\ll S<1$ | $0\ll\gamma<1$ | $\neq 0$ | Active phase wave
$0<r\ll 1$ | $0\ll R<1$ | $0\ll S<1$ | $\neq 0$ | $\neq 0$ | Chimera
$\approx 0$ | $\approx 1$ | $0<S<R$ | $\approx 0$ | $\approx 0$ | Anti-phase sync

###  III.2 Emerging collective states

We investigate the twin activities of synchronization and swarming in our
model. For simplistic purpose, we take the values of $J_{1}$, $J_{2}$, and
$J_{3}$ to be equal to $0.1$ and fix $K_{1}=-0.1$, $K_{2}=-0.2$. These choices
of parameter values are arbitrary and solely made for a case study of our
model. We relax this choice in the subsequent sections. However, the natural
indication after performing numerical simulations is that $K_{3}$ is the most
crucial parameter which determines the inter-group phase coupling. That is why
we keep it as a free parameter and study our model’s collective states while
varying it. The model exhibits six long-term collective states: anti-phase
sync, chimera state, active async, static async, active phase wave, and static
sync when we vary $K_{3}$ inside an interval $[-0.75,0.5]$. Figure 1(b)-(f)
display the states by scatter plots in the $(x,y)$ plane where the
swarmalators are represented by dots and they are colored according to their
phases $\theta$. Figure 1(a) reveals the variation of order parameters as a
function of $K_{3}$. The order parameters $r,S,\gamma,R,V$ are plotted by
blue, red, magenta, purple, and green-colored dotted lines, respectively.
Table 1 provides information regarding the values of the order parameters in
these states. Next, we discuss these collective states and their structural
properties in detail.

We start from the left end point of the interval. Here, $K_{3}\ll
K_{1},K_{2}$. The population breaks into two disjoint clusters formed by the
two groups of swarmalators. Both clusters are stationary in spatial position
and phase. Swarmalators inside each cluster are completely synchronized. But,
one cluster is synchronized at a common phase which is at $\pi$ difference
from the common phase of the other cluster. We call this as anti-phase sync
(see Fig. 1(b)). Look at the white region of the parameter space of Fig. 1(a),
whereby the definition of $R$ in Eq. 5, $R\approx 1$ in this state (purple
curve). Since the overall population’s phases are distributed in $\pi$
difference in two sub-populations and they are equal in size, we get $r\approx
0$ (blue curve). Being a static state, it also gives $\gamma\approx 0$
(magenta curve) and $V\approx 0$ (green curve). The other order parameter $S$
holds a nonzero value that is less than $R$ in this state (red curve). For a
compact view of the order parameters, we refer the reader to Table 1. Section
V.2 presents this state in a more detailed way. See Movie 1 of the
Supplementary Material for the time evolution of this state.

When we gradually increase $K_{3}$, swarmalators in one community remain fully
synchronized in phase but in the other community, asynchrony starts to appear.
The cluster formation in the anti-phase sync state remains intact here (Fig.
1(c)). However, asynchrony in one cluster brings some activity inside that
cluster in the sense that swarmalators now move. This state is best visualized
when studied in terms of $r_{1}$ and $r_{2}$. The synchronized cluster gives
$r_{1}=1$, and the desynchronized one gives $r_{2}<1$. This co-existence of
synchronized and desynchronized swarmalator communities is reminiscent of the
chimera state found in the study of coupled oscillators Abrams et al. (2008);
Abrams and Strogatz (2004); Majhi et al. (2019b) and we simply name this state
as chimera state. All the five order parameters $r,S,\gamma,R,V$ show nonzero
values here (pink region in Fig. 1). See Table 1 for more details. We also
discuss this state in detail in Sec. VI. Movie 2 of the Supplementary Material
demonstrates the time evolution of the chimera state.

On further increment of $K_{3}$ from the chimera state, we encounter the
swarmalators moving and arranging themselves within a circular disc and their
phases are totally incoherent, i.e., $r\approx 0$. The activity never dies and
they keep moving in the two-dimensional plane which gives $V\neq 0$. This
state is named as active async as the swarmalators maintain movement in the
$(x,y)$ plane, and their phases are desynchronized. Find Table 1 for the
description of order parameters in this state. Also, observe the cyan region
in Fig. 1. The activity dies keeping the disc structure with the incoherent
phase nature when $K_{3}$ is increased beyond this state. This is the static
async state. The only difference between this state and the active async state
is that $V=0$ is in this state. Static async state prevails over the yellow
region in Fig. 1. Figure 1(d) represents a snapshot at a particular time
instant for both these states. See Movies 3 & 4 of the Supplementary Material
for the time evolution to this state.

Moving to the right with increasing $K_{3}$ from the static async state, we
observe another collective state where the swarmalators arrange themselves
inside an annular ring and oscillate to achieve regular cycles in both phase
and space. This state was termed as active phase wave in previous studies
O’Keeffe et al. (2017). A snapshot of this motion is best illustrated in Fig.
1(e). By our definition, $\gamma$ is nonzero in this state. Find the green
region of Fig. 1 for the occurrence of this state (Movie 5 of the
Supplementary Material best describes this state). Finally, to the extreme
right of this parameter region where $K_{3}$ is sufficiently large and
positive, phases of the swarmalators throughout the population get
synchronized and they form a disc structure in the plane. This previously
reported state is known as the static sync O’Keeffe et al. (2017). The value
of $r$ is the maximum here which is observed by the blue curve in the purple
region. $R$ is also close to $1$ here by definition. Figure 1(f) illustrates a
snapshot of this state (also see Movie 6 of the Supplementary Material).

Till now we have only varied $K_{3}$ and studied the emerging six collective
states. Now, we simultaneously vary $J_{3}$ along with $K_{3}$ and observe the
dynamical behaviors. The resulting parameter space is shown in Fig. 2. In this
figure, the $J_{3}$-$K_{3}$ parameter plane is divided into $100\times 100$
mesh points, and at each point, we simulate our model for $T=5000$ time units.
The value of order parameter $R$ is calculated over the last $10\%$ data and
the mesh point is colored according to this value. We observe from Fig. 2 that
the emerging states are robust concerning variation in $J_{3}$. The top yellow
region corresponds to the static sync state where $K_{3}$ is positive and
$r,R\approx 1$. The yellow region towards the bottom corresponds to the anti-
phase sync state where $R\approx 1$, but $r\approx 0$ (not shown here). The
red and black curves are the analytical predictions for achieving the static
sync and anti-phase sync state, respectively. Find Sec. V.1 and V.2 for the
derivation of these curves. So far we have always considered $J_{1}=J_{2}$.
Appendix B demonstrates the picture when we work with $J_{1}\neq J_{2}$. The
emerging states remain the same which can be seen from Fig. 9.

Figure 2: $J_{3}$-$K_{3}$ Parameter space for $J_{1}=J_{2}$. (a)
$J_{1}=J_{2}=0.1$. (b) $J_{1}=J_{2}=0.5$. The model is integrated with $N=200$
swarmalators using Heun’s method with step-size $dt=0.01$ for $T=5000$ time
units. Order parameter $R$ is calculated with the last $10\%$ data after
discarding the transients. Colorbar stands for the value of $R$. Red and black
curves are analytical predictions Eqs. (21) and (25), respectively. Here
$K_{1}=-0.1,K_{2}=-0.2$.

##  IV Emerging collective states from identical communities: Dynamical routes

We know from Ref. O’Keeffe et al. (2017) that with a single community
structure, our model exhibits five long-term states depending on the parameter
values. These states are static sync, static async, static phase wave,
splinter phase wave, and active phase wave of which the last two are non-
stationary states. Here we assume that both communities are in the same state
which belongs to one of these five states. This means the communities are
identical with $J_{1}=J_{2}$ and $K_{1}=K_{2}$. Furthermore, we fix
$J_{3}=0.1$ and analyze the routes from static anti-phase sync to static sync
by varying the parameter $K_{3}$ over a range to perceive the collective
states.

Figure 3: Behavior of the order parameters for identical swarmalator
communities. Order parameters as a function of $K_{3}$ where initially both
the communities are in (a) static sync ($J_{1}=J_{2}=0.1$, $K_{1}=K_{2}=1.0$),
(b) static async ($J_{1}=J_{2}=0.1$, $K_{1}=K_{2}=-1.0$), (c) static phase
wave ($J_{1}=J_{2}=1.0$, $K_{1}=K_{2}=0.0$), (d) splintered phase wave
($J_{1}=J_{2}=1.0$, $K_{1}=K_{2}=-0.1$) and (e) active phase wave
($J_{1}=J_{2}=1.0$, $K_{1}=K_{2}=-0.75$). Simulation parameters
$(dt,T,N)=(0.01,5000,200)$. The order parameters are calculated with the last
$10\%$ data. Here, we fix $J_{3}=0.1$.

###  IV.1 Static sync

We start from the static sync state for both the communities
($J_{1}=J_{2}=0.1$, $K_{1}=K_{2}=1.0$) and change $K_{3}$. Firstly, in the
negative $K_{3}$ region, we notice the population forms two clusters that are
static in both phase and position. They are separated by a phase difference of
$\pi$ from each other, which is the anti-phase sync state. An increment of
$K_{3}$ shows the persistence of cluster structure but with a lower phase
difference (see Movie 7 of the Supplementary Material). This cluster
synchronization state (which is not the anti-phase sync state) exists over a
small interval of $K_{3}$ before finally yielding the static sync state. With
increasing $K_{3}$, we find

anti-phase sync $\rightarrow$ cluster sync $\rightarrow$ static sync.

See Fig. 3(a) for the behavior of order parameters here.

###  IV.2 Static async

Here, we take the two groups initially at the static async state by choosing
the parameter values as $J_{1}=J_{2}=0.1$ and $K_{1}=K_{2}=-1.0$ and vary
$K_{3}$. Starting from a relatively lower value of $K_{3}$ in comparison to
$K_{1}$ and $K_{2}$, we notice the presence of an anti-phase sync state where
two synchronized, stable clusters maintain a phase difference of $\pi$. When
we increase $K_{3}$ in very small magnitude, we observe the chimera state
where one group of swarmalators are fully phase coherent and in the other
group they are out of synchrony. We encounter the active async state as we
further increase $K_{3}$. After this, activity dies and the swarmalators
arrange themselves in a static async by adjusting their spatial position with
further increments of $K_{3}$. From this, we spot the emergence of an active
phase wave state by increasing $K_{3}$. As $K_{3}$ is further increased, the
whole community accomplishes themselves in static sync finally. Figure 3(b)
portrays the phase transitions in this case. The route is

anti-phase sync $\rightarrow$ chimera $\rightarrow$ active async $\rightarrow$
static async $\rightarrow$ active phase wave $\rightarrow$ static sync.

Compared to the earlier case where both the communities were in a static sync
state, here we observe that the intermediate dynamics are relatively richer
when the communities are in static async.

###  IV.3 Static phase wave

Here both communities are in a static phase wave state. Primarily, here we
deal with a phase-dependent aggregation model as $J_{1}=J_{2}=1.0$ and
$K_{1}=K_{2}=0$. So, intra-community phase coupling is absent here. Phase
interaction only takes place through the inter-community structure via
$K_{3}$. We can trace the anti-phase sync in the negative $K_{3}$ region as
before. The swarmalators follow a path from anti-phase sync to static sync
through a static phase wave state which is deformed in nature, i.e., they are
distributed in a non-uniform pattern in the $2$D plane (in the existing static
phase wave state, they are distributed uniformly in an annular ring). We can
trace this deformed state very close to the $K_{3}=0$ region. Movie 8 of the
Supplementary Material best describes this state. Swarmalators arrange
themselves into static sync for $K_{3}>0$. Here, the route can be noted down
as:

anti-phase sync $\rightarrow$ deformed static phase wave $\rightarrow$ static
sync.

The order parameters can be found in Fig. 3(c).

###  IV.4 Splintered phase wave

So far we deal with the scenario where both the communities are in static
states initially. Here we start with two identical non-stationary states
namely splinter phase wave. We keep the parameter values $J_{1}=J_{2}=1.0$ and
$K_{1}=K_{2}=-0.1$ and vary $K_{3}$. We observe anti-phase sync where two
static, synchronized clusters persist with a phase difference $\pi$ for a
relatively smaller value of $K_{3}$ compared to $K_{1}$ and $K_{2}$.
Increasing the value of $K_{3}$, we mark the splintered phase wave state
$-0.4<K_{3}<0.22$. Here the swarmalators split into two clusters where the
mean phases of the clusters differ from each other by approximately $\pi$
(Movie 9 of the Supplementary Material). Moving to the right with an
increasing value of $K_{3}$, we notice some of the swarmalators start to
execute a full cycle rotation spatially but their phases do not change from
$0$ to $2\pi$ as in the active phase wave state. This peculiar state can be
deciphered as the simultaneous coexistence of splintered phase wave and active
phase wave states (see Movie 10 of the Supplementary Material for an
illustration of the state). We observe the mixed activity of splintered and
active phase wave states when $-0.22<K_{3}<0.22$. Further increasing $K_{3}$,
the swarmalators are again divided into two clusters but this time they
maintain a difference of mean phases around $0$ (see the time evolution of
this state in Movie 11 of the supplementary Material). This state exists for
$0.22<K_{3}<0.4$. Finally, the whole population reaches static synchrony after
a certain value of $K_{3}$ ($\approx 0.4$). The overall route is depicted as:

anti-phase sync $\rightarrow$ splintered phase wave (mean phase difference
close to $\pi$) $\rightarrow$ mixed (splintered and active phase waves)
$\rightarrow$ splintered phase wave (mean phase difference close to $0$)
$\rightarrow$ static sync.

###  IV.5 Active phase wave

Here the story starts with two identical active phase wave states. The
parameter values are $J_{1}=J_{2}=1.0$ and $K_{1}=K_{2}=-0.75$. To analyze the
route from anti-phase synchrony to static synchrony, we vary $K_{3}$ over a
broad range. In this case, the anti-phase sync state is found for a relatively
larger negative value of $K_{3}$ compared to the previous cases
($K_{3}<-2.12$). Increasing the value of $K_{3}$, swarmalators start to
segregate into two clusters and we observe activity emerging in the system.
Some of the swarmalators undergo a full circle rotation in space and phase and
consequently, $\gamma$ exhibits a small non-zero value around
$-2.16<K_{3}<-1.28$. Movie 12 of the Supplementary Material demonstrates the
state best. On further increment of $K_{3}$, we notice their oscillations
increase in amplitude until all of them start to execute regular cycles in
both phase and spatial angle, i.e., the swarmalators settle in the active
phase wave state. The value of $\gamma$ is close to $1$ and $S$ is very small.
With a further increment of $K_{3}$, their activity begins to diminish
gradually and they are again separated into two clusters. The phase difference
also decreases between the two clusters and $\gamma$ is very small compared to
$1$ near $1.12<K_{3}<2.0$ (see Movie 13 of the Supplementary Material).
Ultimately we find static sync for $K_{3}>2.0$. The route in this case becomes

anti-phase sync $\rightarrow$ mixed (splintered and active phase waves)
$\rightarrow$ active phase wave $\rightarrow$ mixed (splintered and active
phase waves) $\rightarrow$ static sync.

##  V Analytical findings

In the previous section, we explored the dynamic states of our model with
various parameter values. The most striking result that we encountered is the
occurrence of anti-phase sync with a reasonably small value of $K_{3}$ and on
the other hand, a sufficiently large and positive value of $K_{3}$ results in
the whole population in the static sync state. These two static states exist
at the opposite extremes of $K_{3}$ values. Now, we try to establish the
criteria for achieving these states.

###  V.1 Static sync state

Before going into the study of the static sync state, we first analyze the
phase dynamics of our model where spatial positions do not affect the phases.
In that case, the phase equation becomes

$\displaystyle\dot{\theta}_{i}=\omega_{i}+\sum_{m=1}^{2}\dfrac{K_{n,m}}{|C_{m}|}\sum_{j\in
C_{m}\setminus\\{i\\}}\sin(\theta_{j}-\theta_{i}),$ (9)

where $i\in C_{n}$. We move to the continuum limit where
$|C_{p}|\rightarrow\infty$, $p=1,2$. Considering the probability density
function $\rho_{n}(\theta,t)$ of oscillators belonging to the $n$-th group,
the Fokker-Planck equation can be written as

$\dfrac{\partial\rho_{n}}{\partial
t}+\dfrac{\partial}{\partial\theta}(\rho_{n}v_{n})=0,$ (10)

where the velocity $v_{n}(\theta^{n},t)$ is given by

$v_{n}(\theta^{n},t)=\omega+\sum_{m=1}^{2}K_{n,m}\int\sin(\theta^{m}-\theta^{n})\rho_{m}(\theta^{m},t)d\theta^{m}.$
(11)

We define the complex order parameter

$z_{n}(t)=\sum_{m=1}^{2}K_{n,m}\int
e^{i\theta^{m}}\rho_{m}(\theta^{m},t)d\theta^{m}.$ (12)

Using this, Eq. (11) is re-written as

$v_{n}(\theta^{n},t)=\omega+\frac{1}{2i}(z_{n}e^{-i\theta^{n}}-z_{n}^{*}e^{i\theta^{n}}),$
(13)

where $*$ denotes complex conjugate. Following Ott-Antonsen ansatz Ott and
Antonsen (2008), we choose a special class of density functions that has an
invariant manifold of Poisson kernels,

$\rho_{n}(\theta^{n},t)=\frac{1}{2\pi}\left\\{1+\left[\sum_{k=1}^{\infty}[a_{n}(t)e^{i\theta^{n}}]^{k}+c.c.\right]\right\\}.$
(14)

where the unknown function $a_{n}(t)$ must be found self-consistently.
Inserting this form of $\rho_{n}$ given by Eq. (14) into Eq. (10), we find
that $\rho_{n}$ satisfies the Fokker-Planck equation for all harmonics $k$ if
$a_{n}$ satisfies

$\dot{a}_{n}+i\omega
a_{n}+\frac{1}{2}\left[a_{n}^{2}z_{n}-z_{n}^{*}\right]=0.$ (15)

Further inserting Eq. (14) into Eq. (12) and after performing the integration,
the complex order parameter $z_{n}$ is expressed in terms of $a_{n}$ as

$z_{n}(t)=\sum_{m=1}^{2}K_{n,m}a_{m}^{*}(t).$ (16)

Then the amplitude equation for $a_{1}$ becomes

$\displaystyle\dot{a}_{1}=$ $\displaystyle-i\omega
a_{1}-\frac{1}{2}(K_{1,1}a_{1}^{*}+K_{1,2}a_{2}^{*})$
$\displaystyle+\frac{1}{2}(K_{1,1}a_{1}+K_{1,2}a_{2}).$ (17)

Similarly, we find the equation for $\dot{a}_{2}$ by interchanging $1$’s and
$2$’s in Eq. (17). We move to the polar coordinates to rewrite the amplitude
equations by defining $a_{n}=r_{n}e^{-i\phi_{n}}$, $n=1,2$. We further define
$\Phi=\phi_{1}-\phi_{2}$. Substituting these into the amplitude equations and
after simplifying, we get

$\displaystyle\dot{r_{1}}$
$\displaystyle=\frac{1-r_{1}^{2}}{2}(K_{1}r_{1}+K_{3}r_{2}\cos\Phi),$ (18)
$\displaystyle\dot{r_{2}}$
$\displaystyle=\frac{1-r_{2}^{2}}{2}(K_{2}r_{2}+K_{3}r_{1}\cos\Phi),$ (19)
$\displaystyle\dot{\Phi}$
$\displaystyle=-K_{3}\left(\frac{r_{1}^{2}+r_{2}^{2}+2r_{1}^{2}r_{2}^{2}}{2r_{1}r_{2}}\right)\sin\Phi$
(20)

(note that, $K_{1,1}\equiv K_{1}$, $K_{2,2}\equiv K_{2}$, and
$K_{1,2}=K_{2,1}\equiv K_{3}$).  We integrate Eqs. (18)-(20) with initial
conditions $(r_{1}(0),r_{2}(0),\Phi(0))=(0.9,0.9,\pi-0.1)$ and demonstrate
their variation as functions of $K_{3}$ in Fig. 4 where $K_{1}$ and $K_{2}$
are fixed to $-0.1$ and $-0.2$, respectively. When $K_{3}<-0.2$,
$r_{1}=r_{2}=1$ and $\Phi=\pi$ which represents the anti-phase sync state (we
study this state in detail in the next section). On the opposite side, for
$K_{3}>0.2$, we get $r_{1}=r_{2}=1$ and $\Phi=0$ which stand for the sync
state. In the middle region $-0.2<K_{3}<0.2$, chimera-like states appear.

Figure 4: Variation of order parameters $r_{1}$, $r_{2}$ and phase difference
$\Phi$ as functions of $K_{3}$. We integrate Eqs. (18)-(20) starting from
initial conditions $(r_{1}(0),r_{2}(0),\Phi(0))=(0.9,0.9,\pi-0.1)$ for
$T=5000$ time units. Then they are time averaged over the last $10\%$ data and
plotted as functions of $K_{3}$. We fix $K_{1}=-0.1,K_{2}=-0.2$ like in Fig.
1. (a) The variations of $r_{1}$ (blue) and $r_{2}$ (red) are plotted. (b)
Represents the change of $\Phi$ (magenta) with varying $K_{3}$. Figure 5:
Anti-phase sync state. The entire swarmalator population forms two disjoint
clusters in space where the clusters belong to the two communities. Simulation
parameters: $J_{1}=J_{2}=0.5$. $K_{1}=-0.1,K_{2}=-0.2,K_{3}=-1.5$.
$(dt,T,N)=(0.01,5000,200)$. $J_{3}=0.5$ in (a) and (b). (a) Snapshot at
$T=5000$ time units showing the spatial structure of the swarmalators in the
anti-sync state where they are colored according to their phases. (b) The
phases of the swarmalators are plotted against their respective indices at
$T=5000$ time units where red and blue dots correspond to swarmalators
belonging to the first and second communities, respectively. The phase
difference is $\pi$ between the communities. (c) The distance between the
center of masses of the two clusters is plotted as a function of $J_{3}$. Red
dots are simulation results and the black curve indicates the analytical
prediction, Eq. (24).

In the global sync state, the phases of the swarmalators throughout the entire
population become identical, which yields $r_{1}=r_{2}=1$ and $\Phi=0$. This
is the trivial solution to Eqs. (18)-(20) and the Jacobian matrix at this
steady state gives eigenvalues $-2K_{3},-K_{1}-K_{3},-K_{2}-K_{3}$. From this,
the sync state stability condition is achieved as

$K_{3}>\max\\{0,-K_{1},-K_{2}\\}.$ (21)

When we consider the phase dynamics of swarmalators, Eq. (2), the spatial
effect is to be dealt with. However, from numerical simulations, we observe
that when $J_{1}=J_{2}$ and $J_{3}>J_{1}$, the stability condition Eq. (21)
holds for achieving the static sync state. This is demonstrated by the red
lines in Fig. 2. For non-identical values of $J_{1}$ and $J_{2}$, the spatial
distributions of swarmalators in the two communities do not remain the same.
The spatial positions having an impact on the phase dynamics, in turn, affect
the critical $K_{3}$ in Eq. (21). The small deviation from the condition given
by Eq. (21) (plotted by the red curve) is observed in Fig. 9 (Appendix B)
where we present our results with $J_{1}=0.1$ and $J_{2}=0.5$.

###  V.2 Anti-phase sync state

In the anti-phase sync state, the two groups get separated from each other in
the phase component. Their phases are fully synchronized within each group,
but there is a phase difference of $\pi$ between these two groups (see Fig.
5(b)). When $J_{3}$ is absent, i.e., $J_{3}=0$, these two groups arrange
themselves spatially in the form of a disc where these discs overlap. The
radius of these discs depends on the choices of $J_{1}$ and $J_{2}$. But when
the value of $J_{3}$ is nonzero, swarmalators belonging to different groups
start to reduce the attraction between them (since the strength of attraction
between these two groups is $1-J_{3}$ as phase difference is exactly $\pi$).
As a result, these two groups form disjoint clusters in the plane. See Fig.
5(a).

In the anti-phase sync state, $r_{1}=r_{2}=1$ and $\Phi=\pm\pi$. These also
satisfy Eqs. (18)-(20) in the steady state. Linearizing these equations around
this steady state and calculating the Jacobian matrix, yields the eigenvalues
$2K_{3}$, $K_{3}-K_{1}$, and $K_{3}-K_{2}$. This gives the stability condition
of the anti-phase sync state as

$K_{3}<\min\\{0,K_{1},K_{2}\\}.$ (22)

We use Eq. (22) to find the stability condition of the anti-phase sync state
found in our systems defined by Eqs (1)-(2). Since in our model, the phase
dynamics of the swarmalators are influenced by the spatial dynamics, we first
study this effect in the anti-phase sync state. From simulation results, we
find that in the anti-phase sync state with nonzero $J_{3}$, swarmalators form
disjoint clusters in a two-dimensional plane. Swarmalators belonging to group
$C_{1}$ make a cluster among them where their phases are synchronized and the
other cluster is formed by the swarmalators similarly belonging to $C_{2}$.
This can be considered as a two-particle system where swarmalators belonging
to the same group are represented by their center of positions and
synchronized phase Sar et al. (2022). Let $\mathbf{x}_{C_{1}}$ and
$\mathbf{x}_{C_{2}}$ be the center of positions of $C_{1}$ and $C_{2}$ and
$\theta_{C_{1}}$, $\theta_{C_{2}}$ be their synchronized phase angles,
respectively. Then from Eq. (1), we can write

$\displaystyle
0=\Bigg{[}\dfrac{\mathbf{x}_{C_{2}}-\mathbf{x}_{C_{1}}}{|\mathbf{x}_{C_{2}}-\mathbf{x}_{C_{1}}|}\big{(}1+J_{3}\cos(\theta_{C_{2}}-\theta_{C_{1}})\big{)}-\dfrac{\mathbf{x}_{C_{2}}-\mathbf{x}_{C_{1}}}{|\mathbf{x}_{C_{2}}-\mathbf{x}_{C_{1}}|^{2}}\Bigg{]}.$
(23)

This gives us the distance between the center of positions of $C_{1}$ and
$C_{2}$ as

$|\mathbf{x}_{C_{2}}-\mathbf{x}_{C_{1}}|=\frac{1}{1-J_{3}},$ (24)

since $|\mathbf{x}_{C_{2}}-\mathbf{x}_{C_{1}}|\neq 0$ and
$\theta_{C_{2}}-\theta_{C_{1}}=\pm\pi$. This is plotted in the black line in
Fig. 5(c) where the red dots are simulation results. When the swarmalators
form separate groups in spatial positions, their effective phase coupling
strength changes since it depends on the distance between the swarmalators.
This is why Eq. (22) does not stand valid for our model. To find the stability
condition of the anti-phase sync state, we need to investigate the effect of
spatial position carefully. The average distance $R_{1}$ between two particles
in $C_{1}$ can be considered as the half of its diameter (maximum distance
between particles in $C_{1}$) which is a function of $J_{1}$,$J_{2}$, and
$J_{3}$, i.e., $R_{1}(J_{1},J_{2},J_{3})$. Similarly for $C_{2}$ it is
$R_{2}(J_{1},J_{2},J_{3})$. On the other hand, the average distance of the
particle throughout the whole population then becomes
$R_{1}+R_{2}+1/(1-J_{3})=R_{3}$, say. The effective ratio of $K_{3}$ to
$K_{1}$ can be written down as $R_{1}K_{3}/R_{3}K_{1}$ and that of $K_{3}$ to
$K_{2}$ is $R_{2}K_{3}/R_{3}K_{1}$. Then from Eq. (22), we write down the
stability condition of the anti-phase sync state as

$K_{3}<\min\\{0,\frac{R_{3}K_{1}}{R_{1}},\frac{R_{3}K_{2}}{R_{2}}\\}.$ (25)

Due to the complexity of the model, we are unable to find explicit expressions
for $R_{1}$ and $R_{2}$. But from numerical simulations, we observe that these
quantities depend majorly on the values of $J_{1}$ and $J_{2}$ for respective
groups and not on $J_{3}$. We can approximately write $R_{1}\approx
R_{1}(J_{1})$ and $R_{2}\approx R_{2}(J_{2})$. To verify our results, we take
$J_{1}=J_{2}=0.1$ and $K_{1}=-0.1,K_{2}=-0.2$. Numerical simulations suggest
$R_{1}\approx 0.98\approx R_{2}$. The curve defined by Eq. (25) is drawn with
these values and is plotted in black in Fig. 2(a). With $J_{1}=J_{2}=0.5$ and
same $K_{1}$ and $K_{2}$ we find $R_{1}\approx 0.7\approx R_{2}$. The
separatrix curve is again calculated and plotted in black in Fig. 2(b). Both
curves match very well with our numerical results. We also verify our findings
with unequal $J_{1}$ and $J_{2}$ in Appendix B.

##  VI Chimera state

The co-existence of coherence and incoherence is known as chimera state Abrams
et al. (2008); Abrams and Strogatz (2004); Majhi et al. (2019b). We find that
for certain parameter values, there is complete synchrony among one group of
swarmalators but the other group is desynchronized. This means one of $r_{1}$
and $r_{2}$ is $1$ and the other one is strictly less than $1$. We display one
such occurrence of chimera state in Fig. 6. Snapshot at $t=300$ time units
with $N=200$ swarmalators is shown in Fig. 6(a) where the two groups arrange
themselves in the $x$-$y$ plane in the shape of non-overlapping half discs.
Here $r_{1}=1$ and $r_{2}<1$. This is evident when we look at the phases of
the swarmalators in Fig. 6(b). We plot the phases of the swarmalators against
their indices where red and blue dots stand for groups one and two,
respectively.

Figure 6: Chimera state. One of the communities is completely phase
synchronized and the presence of asynchrony is found in the other community.
Simulation parameters: $J_{1}=J_{2}=J_{3}=0.1$. $K_{1}=-0.1$, $K_{2}=-0.2$,
$K_{3}=-0.4$, $(dt,T,N)=(0.01,5000,200)$. (a) Snapshot at $T=5000$ time units
demonstrating the chimera state. (b) Snapshots of the swarmalators’ phases are
plotted against their indices. The red and blue dots refer to the first and
second communities, respectively. Swarmalators are synchronized in the first
community (red dots) but desynchronized in the second one (blue dots).

We further study the nature of the chimera state. For a case study, we fix
$J_{1}=J_{2}=J_{3}=0.1$ and set $K_{1}$ and $K_{2}$ to $-0.1$ and $-0.2$,
respectively. By careful investigation, we find that a chimera state exists
for these parameter values when $-0.56<K_{3}<-0.28$. In the chimera state,
$r_{1}$ stays fixed to $1$ but $r_{2}$ is always less than $1$. Moreover, we
observe oscillation in $r_{2}$, which means it varies with time. So, the
chimera we report in this work is breathing chimera. We establish this by
drawing Fig. 7 where $r_{2}$ is plotted as a function of time for various
values of $K_{3}$. It is to be noted that, with decreasing $K_{3}$ the
magnitude of $r_{2}$ keeps increasing. Eventually around $K_{3}\approx-0.57$,
$r_{2}$ goes to $1$, which is the anti-phase sync state.

Figure 7: Breathing chimera. We delineate the breathing nature of the chimera
state. Parameter values used: $J_{1}=J_{2}=J_{3}=0.1$.
$K_{1}=-0.1,K_{2}=-0.2$. Simulation parameters $(dt,T,N)=(0.01,2000,200)$. The
order parameter $r_{1}$ measuring the phase coherence among the swarmalators
in the first community acquires the value $1$. But $r_{2}$ is less than $1$
due to the presence of asynchrony in the second community. The temporal
evolution of $r_{2}$ is plotted for several $K_{3}$ values, $K_{3}=-0.3$
(blue), $-0.35$ (red), $-0.4$ (yellow), and $-0.57$ (magenta). We observe
oscillatory behavior of $r_{2}$ which reveals the breathing nature of the
chimera state. The magnitude of $r_{2}$ increases and the oscillation decays
with decreasing $K_{3}$ until it reaches the maximum value $1$ where the
oscillation completely dies.

##  VII Conclusion

The phase-dependent spatial aggregation and position-dependent phase
synchronization are at the core of swarmalator dynamics. Swarmalators endowed
with spatial and phase interactions are competent to exhibit complex
collective behaviors. These states can be found in real-world systems like
Japanese tree frogs Aihara et al. (2014), magnetic domain walls Hrabec et al.
(2018), Janus matchsticks Chaudhary et al. (2014), robotic swarms Schilcher et
al. (2021); Barciś and Bettstetter (2020) etc. To this end, studies are being
carried out on swarmalator models by defining suitable interaction functions,
network structures, coupling schemes, etc. (We refer the reader to Sar and
Ghosh (2022) for a recent review on swarmalator systems.)

In this article, we have studied a population of swarmalators where they are
distributed in two communities. The intra and inter-community coupling
strengths have been carefully varied to observe different emerging states. Two
of them, viz., the anti-phase sync and the chimera state are not commonly
observed in swarmalator systems and to the best of our knowledge, have not
been studied rigorously (the anti-phase state has been reported previously in
Sar et al. (2022); Ceron et al. (2023) and chimera like states were observed
in Hong et al. (2021)). The novelty of our work lies in the fact that we have
found an anti-phase sync state with all the intra and inter-phase coupling
strengths being negative. It can be inferred that the imposed community
structure is responsible for this. The chimera state encountered is also due
to the interplay between swarmalators belonging to different communities.
Although we were not able to provide any mathematical formulation for the
chimera state, our model still can be used as a testbed for future works on
chimera states in swarmalator systems.

We have also conspicuously illustrated the phase transitions by varying the
inter-community phase coupling strength $K_{3}$. The emerging states are
characterized in terms of some order parameters. Anti-phase sync state is
perceived for a sufficiently small (negative) value of $K_{3}$ and the sync
state is detected for a positive large value of it. We study these two states
in detail and provide semi-analytical conditions for achieving these states.
We also study the different routes from the anti-phase sync state to the sync
state by assuming that the two communities are identical to start with.
Moreover, we have established our results when the parameters $J_{3}$ and
$K_{3}$ are varied simultaneously.

We can highlight the limitation of our work by pointing at the inability to
explicitly incorporate the spatial dynamics in the analysis of the anti-phase
sync and static sync state. This might be wiped out if some simpler type of
spatial interaction functions is used other than the power laws used in our
model. It also remains to be seen what happens when more than two communities
are considered. The model can be simplified by reducing the spatial dimension
placing the swarmalators on a ring and then imposing the community
interactions.  Future works can also be carried out with our model by
considering nonidentical swarmalators by drawing frequencies from Gaussian or
Lorentzian distributions. Through preliminary inspection, we observed that
some of the emerging states that we reported here (static async, static sync)
will have their analogous counterparts for nonidentical swarmalators. But for
the existence of other states like anti-phase, chimera, etc., a deep and
systematic investigation is required.

## Data availability statement

The data that support the findings of this study are openly available in the
GitHub repository web .

##  Appendix A Unequal community sizes

Figure 8: Phase transition with unequal community sizes. Simulations are
performed for a total population of $N=200$ swarmalators. (a) $p_{1}=0.6$ and
(b) $p_{1}=0.7$. Other parameter values used are $J_{1}=J_{2}=J_{3}=0.1$ and
$K_{1}=-0.1,K_{2}=-0.2$. We observe the same states as in Fig. 3 with equal
community sizes. It establishes the fact that our reported results are robust
and independent of the initial distribution of swarmalators in the
communities.

In the main text of our paper, we discussed the case where the communities are
of equal size and studied different states. Here, we cover the scenario where
the two communities have unequal sizes. The total population size is $N$.
These swarmalators are distributed in two communities. Let $p_{1}$ and $p_{2}$
denote the probabilities that the $i$-th swarmalator belongs to the first, and
second communities, respectively. Clearly, $p_{1}+p_{2}=1$. For equal
community sizes, $p_{1}$ is essentially equal to $p_{2}$. Here, we take
$p_{1}\neq p_{2}$ so that the communities are unequal in size. We study two
cases, one where $p_{1}=0.6$ and the other one $p_{1}=0.7$. The parameter
values are $J_{1}=J_{2}=J_{3}=0.1$, $K_{1}=-0.1$, and $K_{3}=-0.2$, the same
values which were used in Sec. III.2. In both cases, what we observe that the
same qualitative behavior of all the order parameters. As a result, the
emerging states remain unaltered. In Fig. 8, we have shown the phase
transition. In the case of equal population sizes, the order parameter $r$ is
approximately zero in anti-phase sync. Due to an equal number of swarmalators
in each group, the terms within the summation in Eq. (3) nullify each other.
But if we choose unequal sub-populations, $r$ has a non-zero value depending
on the ratio of swarmalators.

##  Appendix B $J_{1}\neq J_{2}$

We study the case where the $J$’s (phase-dependent spatial coupling strengths
among communities) are not equal i.e., $J_{1}\neq J_{2}$. For instance, we
take $J_{1}=0.1,J_{2}=0.5$. $K_{1}$ and $K_{2}$ are kept fixed at $-0.1$ and
$-0.2$, respectively. The resulting behavior is demonstrated through Fig. 9.
The overall collective states remain the same. It can be observed if we
compare Fig. 9 with Fig. 2 (where $J_{1}=J_{2}$).

Figure 9: $J_{3}$-$K_{3}$ parameter space for $J_{1}\neq J_{2}$. Here
$K_{1}=-0.1$, $K_{2}=-0.2$, $J_{1}=0.1$, $J_{2}=0.5$. The system is integrated
with $N=200$ swarmalators using Heun’s method with step-size $dt=0.01$ for
$T=5000$ time units. Order parameter $R$ is calculated with the last $10\%$
data after discarding the transients. Colorbar stands for the value of $R$.
Red and black curves are drawn using Eqs. (21) and (25), respectively.

##  Appendix C Nonidentical swarmalators

For our study, we have considered swarmalators with identical frequencies in
both the communities, i.e., $\omega_{i}=\omega$ for $i=1,2,\ldots,N$, and it
is further set to zero by moving to a proper reference frame. Here, we draw
the frequencies from the Gaussian distribution with zero mean and unit
standard deviation to make the nonidentical swarmalators. We observe that the
sync state takes place for a larger inter-community phase coupling strength
$K_{3}$ compared to identical swarmalators. The phases never become static and
keep evolving with time which is seen via Fig. 10(a)-(c) where snapshots are
taken at different time units. For small coupling strength $K_{3}$, the async
state is realized. Here also, the phases are non-stationary. In Fig.
10(d)-(f), snapshots of the async state are shown at $T=2000,3500$, and $5000$
time units, respectively. However, we were unable to detect the emergence of
anti-phase and chimera states. A rigorous study through minute exploration of
the parameters is needed when one considers nonidentical swarmalators.

Figure 10: Sync and async states for nonidentical swarmalators. We choose
$\omega_{i}$ randomly from Gaussian distribution centered at zero with
standard deviation $1$. Simulation parameters: $J_{1}=J_{2}=J_{3}=0.5$.
$K_{1}=-0.1,K_{2}=-0.2$. $(dt,N)=(0.01,200)$. Snapshots of the sync state are
shown at (a) $T=2000$, (b) $T=3500$, and (c) $T=5000$ time units where
$K_{3}=4.0$. The async state is shown for $K_{3}=0.5$ at (d) $T=2000$, (e)
$T=3500$, and (f) $T=5000$ time units. Swarmalators are colored according to
their phases. The numbers of swarmalators in both communities are same here.

## References

  * Pecora and Carroll (1990) L. M. Pecora and T. L. Carroll, Physical Review Letters 64, 821 (1990).
  * Boccaletti et al. (2002) S. Boccaletti, J. Kurths, G. Osipov, D. Valladares, and C. Zhou, Physics Reports 366, 1 (2002).
  * Arenas et al. (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Physics Reports 469, 93 (2008).
  * Rosenblum and Kurths (2003) M. Rosenblum and J. Kurths, _Synchronization: a universal concept in nonlinear science_ (Cambridge University Press, 2003).
  * Jusup et al. (2022) M. Jusup, P. Holme, K. Kanazawa, M. Takayasu, I. Romić, Z. Wang, S. Geček, T. Lipić, B. Podobnik, L. Wang, et al., Physics Reports 948, 1 (2022).
  * Buck (1988) J. Buck, The Quarterly Review of Biology 63, 265 (1988).
  * Wiesenfeld et al. (1996) K. Wiesenfeld, P. Colet, and S. H. Strogatz, Physical Review Letters 76, 404 (1996).
  * Montbrió et al. (2015) E. Montbrió, D. Pazó, and A. Roxin, Physical Review X 5, 021028 (2015).
  * Buhl et al. (2006) J. Buhl, D. J. Sumpter, I. D. Couzin, J. J. Hale, E. Despland, E. R. Miller, and S. J. Simpson, Science 312, 1402 (2006).
  * Couzin (2007) I. Couzin, Nature 445, 715 (2007).
  * Ballerini et al. (2008) M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, et al., Proceedings of the National Academy of Sciences 105, 1232 (2008).
  * Bialek et al. (2012) W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. M. Walczak, Proceedings of the National Academy of Sciences 109, 4786 (2012).
  * Leonard et al. (2012) N. E. Leonard, T. Shen, B. Nabet, L. Scardovi, I. D. Couzin, and S. A. Levin, Proceedings of the National Academy of Sciences 109, 227 (2012).
  * Garcimartín et al. (2015) A. Garcimartín, J. Pastor, L. Ferrer, J. Ramos, C. Martín-Gómez, and I. Zuriguel, Physical Review E 91, 022808 (2015).
  * Zhang et al. (2010) H.-P. Zhang, A. Be’er, E.-L. Florin, and H. L. Swinney, Proceedings of the National Academy of Sciences 107, 13626 (2010).
  * Levy and Requeijo (2008) D. Levy and T. Requeijo, Bulletin of Mathematical Biology 70, 1684 (2008).
  * Tian et al. (2021) M. Tian, C. Zhang, R. Zhang, and J. Yuan, Biophysical Journal 120, 1615 (2021).
  * Fetecau et al. (2011) R. C. Fetecau, Y. Huang, and T. Kolokolnikov, Nonlinearity 24, 2681 (2011).
  * Darnton et al. (2010) N. C. Darnton, L. Turner, S. Rojevsky, and H. C. Berg, Biophysical Journal 98, 2082 (2010).
  * Mogilner and Edelstein-Keshet (1999) A. Mogilner and L. Edelstein-Keshet, Journal of Mathematical Biology 38, 534 (1999).
  * Frasca et al. (2008) M. Frasca, A. Buscarino, A. Rizzo, L. Fortuna, and S. Boccaletti, Physical Review Letters 100, 044102 (2008).
  * Majhi et al. (2019a) S. Majhi, D. Ghosh, and J. Kurths, Physical Review E 99, 012308 (2019a).
  * Peshkov et al. (2022) A. Peshkov, S. McGaffigan, and A. C. Quillen, Soft Matter 18, 1174 (2022).
  * Aihara et al. (2014) I. Aihara, T. Mizumoto, T. Otsuka, H. Awano, K. Nagira, H. G. Okuno, and K. Aihara, Scientific Reports 4, 3891 (2014).
  * Tan et al. (2022) T. H. Tan, A. Mietke, J. Li, Y. Chen, H. Higinbotham, P. J. Foster, S. Gokhale, J. Dunkel, and N. Fakhri, Nature 607, 287 (2022).
  * O’Keeffe et al. (2017) K. P. O’Keeffe, H. Hong, and S. H. Strogatz, Nature Communications 8, 1504 (2017).
  * Sar and Ghosh (2022) G. K. K. Sar and D. Ghosh, Europhysics Letters 139, 53001 (2022).
  * Yoon et al. (2022) S. Yoon, K. O’Keeffe, J. Mendes, and A. Goltsev, Physical Review Letters 129, 208002 (2022).
  * Tanaka (2007) D. Tanaka, Physical Review Letters 99, 134103 (2007).
  * Iwasa et al. (2010) M. Iwasa, K. Iida, and D. Tanaka, Physical Review E 81, 046220 (2010).
  * O’Keeffe et al. (2018) K. P. O’Keeffe, J. H. Evers, and T. Kolokolnikov, Physical Review E 98, 022203 (2018).
  * Lizarraga and de Aguiar (2020) J. U. Lizarraga and M. A. de Aguiar, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 053112 (2020).
  * Hong et al. (2023) H. Hong, K. P. O’Keeffe, J. S. Lee, and H. Park, Physical Review Research 5, 023105 (2023).
  * O’Keeffe and Hong (2022) K. O’Keeffe and H. Hong, Physical Review E 105, 064208 (2022).
  * Blum et al. (2022) N. Blum, A. Li, K. O’Keeffe, and O. Kogan, arXiv preprint arXiv:2210.11417 (2022).
  * Lee et al. (2021) H. K. Lee, K. Yeo, and H. Hong, Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 033134 (2021).
  * Sar et al. (2022) G. K. Sar, S. N. Chowdhury, M. Perc, and D. Ghosh, New Journal of Physics 24, 043004 (2022).
  * Sar et al. (2023) G. K. Sar, D. Ghosh, and K. O’Keeffe, Physical Review E 107, 024215 (2023).
  * Ceron et al. (2023) S. Ceron, K. O’Keeffe, and K. Petersen, Nature Communications 14, 940 (2023).
  * Girvan and Newman (2002) M. Girvan and M. E. Newman, Proceedings of the National Academy of Sciences 99, 7821 (2002).
  * Fortunato (2010) S. Fortunato, Physics Reports 486, 75 (2010).
  * Newman (2006) M. E. Newman, Proceedings of the National Academy of Sciences 103, 8577 (2006).
  * Malliaros and Vazirgiannis (2013) F. D. Malliaros and M. Vazirgiannis, Physics Reports 533, 95 (2013).
  * Clauset et al. (2004) A. Clauset, M. E. Newman, and C. Moore, Physical Review E 70, 066111 (2004).
  * Newman and Girvan (2004) M. E. Newman and M. Girvan, Physical Review E 69, 026113 (2004).
  * Tsimring et al. (2005) L. Tsimring, N. Rulkov, M. Larsen, and M. Gabbay, Physical Review Letters 95, 014101 (2005).
  * Abrams et al. (2008) D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A. Wiley, Physical Review Letters 101, 084103 (2008).
  * Abrams and Strogatz (2004) D. M. Abrams and S. H. Strogatz, Physical Review Letters 93, 174102 (2004).
  * Majhi et al. (2019b) S. Majhi, B. K. Bera, D. Ghosh, and M. Perc, Physics of Life Reviews 28, 100 (2019b).
  * Ott and Antonsen (2008) E. Ott and T. M. Antonsen, Chaos: An Interdisciplinary Journal of Nonlinear Science 18, 037113 (2008).
  * Hrabec et al. (2018) A. Hrabec, V. Křižáková, S. Pizzini, J. Sampaio, A. Thiaville, S. Rohart, and J. Vogel, Physical Review Letters 120, 227204 (2018).
  * Chaudhary et al. (2014) K. Chaudhary, J. J. Juárez, Q. Chen, S. Granick, and J. A. Lewis, Soft Matter 10, 1320 (2014).
  * Schilcher et al. (2021) U. Schilcher, J. F. Schmidt, A. Vogell, and C. Bettstetter, in _2021 IEEE International Conference on Autonomic Computing and Self-Organizing Systems (ACSOS)_ (IEEE, 2021), pp. 90–99.
  * Barciś and Bettstetter (2020) A. Barciś and C. Bettstetter, IEEE Access 8, 218752 (2020).
  * Hong et al. (2021) H. Hong, K. Yeo, and H. K. Lee, Physical Review E 104, 044214 (2021).
  * (56) https://github.com/gourab-sar/Anti-phase-synchronization-in-a-population-of-swarmalators.


# Compact extra dimensions as the source of
primordial black holes

Valery V. Nikulin1
<EMAIL_ADDRESS>& Maxim A. Krasnov1
<EMAIL_ADDRESS>& Sergey G. Rubin1,2
<EMAIL_ADDRESS>
( 1 National Research Nuclear University MEPhI (Moscow Engineering Physics
Institute),
115409, Kashirskoe shosse 31, Moscow, Russia
2 N. I. Lobachevsky Institute of Mathematics and Mechanics, Kazan Federal
University,
420008, Kremlevskaya street 18, Kazan, Russia)

###### Abstract

We discuss the model of the primordial black holes formation at the reheating
stage. These small massive black holes appear due to specific properties of
the compact extra dimensions. The latter gives rise the low energy model
containing the effective scalar field potential capable for the domain walls
production. Formed during inflation, these walls are quite dense so they
collapse soon after inflation ends. The discussion is performed within the
scope of multidimensional $f(R)$-gravity.

_Keywords_ extra dimensions, modified gravity, $f(R)$-gravity, primordial
black holes, cosmology, domain wall.

##  1 Introduction

Extra dimensions are usually studied within the framework of elementary
particle physics [1], for example in the context of the unification of
interactions [2, 3], explaining the nature of the Standard Model fields [4],
searching for their manifestations in collider experiments [5, 6]. This paper
explores another possible cosmological consequence — shows that compact extra
dimensions may be the cause of the primordial black holes (PBHs) formation
immediately after the end of the inflation.

It is known that one of the central tasks of theories with compact extra
dimensions is to provide their compactification and stabilization [7] during
cosmological evolution. This can be done, for example, by introducing
additional scalar fields [8] or $f(R)$-modification of gravity [9, 10]. The
latter approach is particularly promising because the Starobinsky quadratic
$f(R)$-gravity [11, 12] gives the best fit to observational constraints on the
parameters of the inflation [13]. Moreover, in multidimensional
$f(R)$-gravity, the processes of cosmological inflation and compactification
are manifestations of general gravitational dynamics in different subspaces
[14].

The possibilities of $f(R)$-gravity are widely studied [15, 16], they offer
solutions to many cosmological problems [17, 18, 19, 20]. One of the problems
that $f(R)$-gravity can solve is the existence of primordial black holes.
Today, the primordial origin of some discovered black holes (quasars at small
$z$ [21, 22], BHs of intermediate masses detected by gravitational-wave
observatories [23]) is hotly discussed [24, 25, 26]. In this paper we
demonstrate how the primordial black holes can appear as a result of
inflationary dynamics in the framework of $f(R)$-gravity model.

The idea of our proposed mechanism is based on the known possibility of domain
walls formation during cosmological inflation followed by their collapse into
primordial black holes [27, 28]. The formation of such domain walls requires a
scalar field with a nontrivial potential containing several vacuums. Such kind
of scalar field effectively arises in the multidimensional $f(R)$-models in
the Einstein frame [29, 10, 14]. This field controls the size of the compact
extra space, and its different vacuums correspond to different universes. In
this paper we calculate the parameters of the domain walls formed by the field
and conclude that appearing at inflationary stage they will immediately
collapse into PBHs during reheating. For a remote observer in the Jordan frame
the appearance of such PBH is interpreted as a manifestation of non-trivial
$f(R)$-gravitational dynamics of multidimensional space. These formed PBHs
will grow rapidly in the process of further cosmological evolution due to
accretion and are capable to turn into observable supermassive quasars at
small $z$.

##  2 Model

The study is based on the modified $f(R)$ gravity acting in $D=4+n$
dimensions. It is described by the action111 In this paper we use the
following conventions for the Riemann curvature tensor
$R_{\mu\nu\alpha}^{\beta}=\partial_{\alpha}\Gamma_{\mu\nu}^{\beta}-\partial_{\nu}\Gamma_{\mu\alpha}^{\beta}+\Gamma_{\sigma\alpha}^{\beta}\Gamma_{\nu\mu}^{\sigma}-\Gamma_{\sigma\nu}^{\beta}\Gamma_{\mu\alpha}^{\sigma}$,
and the Ricci tensor is defined as follows:
$R_{\mu\nu}=R^{\alpha}_{\mu\alpha\nu}\,.$

$S=\frac{m_{D}^{D-2}}{2}\int
d^{4+n}x\sqrt{|g_{D}|}\,\left[f(R)+c_{1}R_{AB}R^{AB}+c_{2}R_{ABCD}R^{ABCD}\right]\,,$
$f(R)=a_{2}R^{2}+R-2\Lambda_{D}\,,$ (1)

where $m_{D}$ is the multidimensional Planck mass. The parameters
$a_{2},c_{1},c_{2}$ of the Lagrangian (2) have dimensionality $[m_{D}^{-2}]$,
the cosmological constant $\Lambda$ has dimensionality $[m_{D}^{2}]$.
Hereafter we will work everywhere in units $m_{D}\equiv 1$, unless it is
explicitly specified. In this work we will represent multidimensional space as
the direct product of $\mathbb{M}=\mathbb{M}_{4}\times\mathbb{M}_{n}$, where
$\mathbb{M}_{4}$ is– four-dimensional space, $\mathbb{M}_{n}$ – an $n$-dim
compact extra space with metric

$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}-e^{2\beta(t,x)}d\Omega^{2}_{n}\,.$ (2)

Here $g_{\mu\nu}$ is a 4-dimensional metric $\mathbb{M}_{4}$, $\beta$ is a
metric function and $d\Omega^{2}_{n}$ is the linear element square of a
maximally symmetric compact extra space $\mathbb{M}_{n}$ with positive
curvature.

Cosmological scenarios based on this theory [29, 14, 30]. are finished by an
effective 4-dim theory at the low energy. Its properties are defined by the
Lagrangian parameters of of action (2). The procedure for obtaining such a
theory is also described in [10].

Following this procedure, we assume maximal symmetry of the extra space
$\mathbb{M}_{n}$. Its radius $\rho\equiv e^{\beta}$ does not depend on the
extra coordinates. In addition, the useful approximation of this effective
theory works in a region where the extra curvature $R_{n}=n(n-1)/e^{2\beta}$
is large compared to the 4-dimensional curvature and changes slowly:

$R=R_{4}+R_{n}+P_{k}\,,\quad
P_{k}=2n\,\partial^{2}\beta+n(n+1)(\partial\beta)^{2}\,,$ $R_{4},P_{k}\ll
R_{n}\,.$ (3)

where $R_{4},R_{n}$ are Ricci scalars for $\mathbb{M}_{4},\mathbb{M}_{n}$. The
approximation (3) applied to action (2) leads to the effective 4-dim model
with a non-minimal coupling between the observed 4-dimensional gravity $R_{4}$
and the scalar field [10]. Therefore the resulting theory will be written in
the Jordan frame and we will use this frame as the physically observed one
[31].

##  3 Effective scalar field action

The minimal coupling between the observed 4-dimensional gravity and the scalar
field is achieved by the transition to the Einstein frame, which strongly
facilitate the analysis. This transition has been performed in detail in [10],
[32]:

$S=\frac{m_{4}^{2}}{2}\int d^{4}x\sqrt{-g_{4}}\
\text{sign}(f^{\prime})\left[R_{4}+K(\phi)(\partial\phi)^{2}-2V(\phi)\right]\,,$
(4)

where the Planck mass in the Einstein frame is
$m_{4}=\sqrt{2\pi^{\frac{n+1}{2}}/\Gamma(\frac{n+1}{2})}$, and
$g_{4}^{\mu\nu}$ is the observed 4-dimensional metric. The Ricci scalar of the
extra space is perceived as the scalar field, $\phi\equiv R_{n}$.

The action (4) contains a potential and a nontrivial kinetic term, which are
expressed through the initial parameters of the Lagrangian (2) (see [14]):

$K(\phi)=\frac{1}{4\phi^{2}}\left[6\phi^{2}\left(\frac{f^{\prime\prime}}{f^{\prime}}\right)^{2}-2n\phi\left(\frac{f^{\prime\prime}}{f^{\prime}}\right)+\frac{n(n+2)}{2}\right]+\frac{c_{1}+c_{2}}{f^{\prime}\phi}\,,$
(5)
$V(\phi)=-\frac{\operatorname{sign}\left(f^{\prime}\right)}{2(f^{\prime})^{2}}\left[\frac{|\phi|}{n(n-1)}\right]^{n/2}\left[f(\phi)+\frac{c_{1}+2c_{2}/(n-1)}{n}\phi^{2}\right]\,.$
(6)

The transition validity from (2) to (4) dictated by ineaqualities (3) will be
analyzed later.

For a wide range of parameters, the potential term (6) has two minima
corresponding to different vacuums of the Universe (Fig. 1). The value of the
potential at these minima must be almost zero according to almost zero value
of the cosmological constant $\Lambda_{4}$. It leads to a relation of the
parameter $c=1/4(a+c_{v}/n)$, where $c_{v}=c_{1}+2c_{2}/(n-1)$.

Figure 1: Graphs of the potential term and kinetic factor (5), (6) for
parameters: $n=6$, $c_{1}=-8000$, $c_{2}=-5000$, $a_{2}=-500$. The left
minimum of the potential is at $\phi=0$, but it is impossible to roll down to
it in a finite time because of the increasing friction effect of the kinetic
factor, the right minimum is at $\phi_{min}=2.3\cdot 10^{-4}.$

The rolling of the field $\phi\equiv R_{n}$ into the right minimum of the
potential (6) corresponds to the stabilization of the compact extra space
$R_{n}\neq 0$ (the extra space is compactified and has some radius
$\rho_{0}$), and leads to the observed cosmology [14]. The presence of the
left minimum $R_{n}\approx 0$ suggests the possibility of another scenario
[30], in which the extra space is unstable and expands to macroscopic scale.

The non-trivial kinetic factor (5) significantly modifies the character of the
field evolution $\phi(t)$ in comparison to the standard scalar-field theory,
providing an increasing friction when rolling to the left minimum, see Fig. 1.
One can simplify the Lagrangian by substituting

$\psi=m_{4}\int\limits_{\phi_{0}}^{\phi}\sqrt{K(\phi^{\prime})}\,d\phi^{\prime}\,,\quad\tilde{V}(\psi)=m_{4}^{2}\,V(\phi(\psi)),\quad
K(\phi^{\prime})>0\,.$ (7)

In this case, $\psi(\phi)$ is monotonic and invertible (which is required to
find the potential in the expression (8)). After simple algebra, the
Lagrangian is reduced to the standard form

$S=\frac{m_{4}^{2}}{2}\int d^{4}x\sqrt{-g_{4}}\,R_{4}+\int
d^{4}x\sqrt{-g_{4}}\,\left[\frac{1}{2}(\partial\psi)^{2}-\tilde{V}(\psi)\right]\,,$
(8)

where $m_{4}^{2}=16\pi^{3}/{15}$ for chosen extra space dimensionality $n=6$.

Figure 2: Graph of the potential term $\tilde{V}(\psi)$ of the effective
theory (8) for the parameters chosen in Fig. 1. The right minimum is at
$\psi_{min}=84.3$. Strictly speaking, the left minimum is at $\psi=-\infty$,
but for our purposes a very accurate approximation $\psi=0$ is sufficient
(corresponding to $\phi=\phi_{0}$), up to which the effective theory (4) still
works. The maximum value of the potential in the transition region between
vacuums $\tilde{V}_{max}\approx 4.1\cdot 10^{-20}$ and the steepness of the
slope $\sqrt{\tilde{V}^{\prime\prime}(\psi_{max})}=2.7\cdot 10^{-11}$.

##  4 Domain walls

It is well known that potentials like (2), containing several minima
(vacuums), can lead to formation of non-trivial field configurations [33] —
"bubbles" of one vacuum inside another, surrounded by a domain wall.

Let us investigate such a configuration by deriving the field equation for
$\psi$ from the effective action (8). For simplicity, we consider it
spherically symmetric and static, which gives the equation:

$\psi_{uu}+\frac{2\psi_{u}}{u}-\tilde{V}^{\prime}(\psi)=0\,,$ (9)

where $u$ is the radial coordinate. When considering a sufficiently large
"bubble" (such that its radius is much larger than the characteristic
thickness of the domain wall $\psi_{uu}\gg 2\psi_{u}/u$) it can easily be
reduced to a first-order equation

$\partial_{u}\psi=\pm\sqrt{2\tilde{V}(\psi)}\,.$ (10)

The characteristic solution of the equation (10) connecting the left vacuum of
the potential $\tilde{V}(\psi)$ to the right one (Fig. 2) is shown in Fig. 3
(blue line).

Figure 3: Numerical solution of (9) under boundary conditions: the left
minimum $\tilde{V}(\psi)$ (Fig. 2) is reached inside the bubble $\psi(0)=0$,
and the right minimum is reached for the remote observer
$\psi(u\rightarrow\infty)=\psi_{min}$ (and forms our Universe). In the graphs:
the blue line is the direct solution of $\psi(u)$ (domain wall) and the red
dashed line is its corresponding function $\phi(u)$. In this example the
bubble radius is $u_{0}\approx 8\cdot 10^{13}$.

We can calculate the domain wall energy density $\varepsilon_{\psi}$ as a
component $T^{00}$ of energy-momentum tensor for the scalar field Lagrangian
$\psi$:

$\varepsilon_{\psi}(u)=T^{00}(u)=\cfrac{\partial\mathcal{L}_{\psi}}{\partial\left(\partial_{0}\psi\right)}\,\partial^{0}\psi-\mathcal{L}_{\psi}g^{00}=\frac{1}{2}(\partial_{u}\psi)^{2}+V(\psi)=2\tilde{V}(\psi(u)).$
(11)

Integration over the radial coordinate (11) gives the surface energy density
of the domain wall in multidimensional Planck units ($m_{D}=1$):

$\sigma=\int\limits_{0}^{\infty}\varepsilon_{\psi}(u)\,du=\int\limits_{0}^{\psi_{min}}\frac{2\tilde{V}(\psi)}{\psi_{u}}\,d\psi=\int\limits_{0}^{\psi_{min}}\sqrt{2\tilde{V}(\psi)}\,d\psi\,.$
(12)

We can also formally estimate the characteristic wall thickness $\delta$ as

$\displaystyle\sigma=$
$\displaystyle\int\limits_{0}^{\infty}2\tilde{V}(\psi)\,du\approx
2\frac{\tilde{V}_{max}}{2}\,\delta\quad\implies\quad\delta\approx\frac{\sigma}{\tilde{V}_{max}}\,.$
(13)

Figure 4: Energy density (11) profile along the radial coordinate $u$ for the
solution shown in Fig. 3. In the case of spherical symmetry, the domain wall
forms a bubble with surface energy density $\sigma\approx 5\cdot 10^{-9}$ and
wall thickness $\delta\approx 1.2\cdot 10^{11}$, according to (12), (13).

Domain walls considered in this paper appears to be very massive (Fig. 4). As
a result, they could be a source of primordial black holes. The mechanism of
formation and collapse of such domain walls is well studied in [34, 35, 27]
and leads to the formation of primordial black holes in a wide range of
masses.

##  5 Generation of domain walls during inflation

As stated above, the rolling of the $\psi$ field to the right minimum of Fig.
2, creates the observable Universe. For its formation, a mechanism of the
cosmological inflation is required. In order not to complicate our
consideration, we will consider the inflation as an external process with a
characteristic Hubble parameter $H$.

The mechanism of inflationary production of the bubbles of alternative vacuum
is well known [27]. As a result of repeated quantum fluctuations, the field
$\psi$ can be flipped from the region of rolling down to the right minimum to
the region of rolling down to the left minimum in some area of the
inflationary Universe. Thy size of this region is growing during inflation
while the $\psi$ field is frozen near the potential maximum. After the end of
inflation, the field tends to one minimum inside the bubble and another
minimum outside it. Increasing energy density gradually forms the domain wall
(Fig. 4) around the bubble.

Several constraints must be imposed on the working model of domain wall
production in the considered $f(R)$-gravity:

  1. 1.

During cosmological inflation, 4-dimensional space is described by a de Sitter
metric with curvature $R_{4}\simeq 12H^{2}$, where $H$ is the Hubble
parameter. Therefore, the approximation (3) is applicable only for field
values $\phi\equiv R_{n}\gg R_{4}\simeq 12H^{2}$. From Fig. 1 we see that
$\phi\sim 10^{-6}-10^{-4}$, so to satisfy the constraint we need $H\ll
10^{-4}$.

  2. 2.

To generate walls via quantum fluctuations of field $\psi$ (near the maximum
of potential 2) its slow rolling is required:
$\sqrt{\tilde{V}^{\prime\prime}(\psi_{max})}\ll H$. According to Fig. 2 we see
that $\sqrt{\tilde{V}^{\prime\prime}(\psi_{max})}\sim 10^{-11}$, so to satisfy
the constraint we need $H\gg 10^{-11}$.

  3. 3.

Domain walls should not be too dense so as not to dominate the inflaton:
$\varepsilon_{\psi}\ll\varepsilon_{inf}\sim H^{2}m_{4}^{2}$. From Fig. 4 we
see that $\varepsilon_{\psi}\sim 10^{-20}$, so to satisfy the constraint we
need $H\gg 10^{-11}$.

  4. 4.

The fluctuations of the $\psi$ field during inflation should not be too large
to prevent overproduction of the domain walls: $\delta\psi=H/2\pi\ll\psi$. As
one can see from Fig. 2, $\psi\sim 10$, so to satisfy the constraint we need
$H\ll 10$.

All constraints above are satisfied if the inflation has the characteristic
scale $10^{-11}\ll H\ll 10^{-4}$. Here $m_{4}=5.75$ for $n=6$ according to
(4). Remind that all estimations are given in the Einstein frame and if the
dimensionality is not explicitly specified, all calculations are carried out
in $m_{D}\equiv 1$ units.

In the previous sections, all calculations were performed in the Einstein
frame. As a physically observable one it is often accepted to consider the
Jordan frame, in which it is possible to manipulate the 4-dimensional Planck
scale $M_{4}=\Omega(\phi_{\text{min}})m_{4}$ by adjusting the size of extra
space. So we will transform final physical results to the Jordan frame.
Cosmological inflation should be described in a physically observable frame —
that is the Jordan frame (let’s denote it by the index $J$):

$S^{J}_{\text{inf}}=\int
d^{4}x\sqrt{-g_{4}^{J}}\,\left[\frac{1}{2}\left(\partial\chi^{J}\right)^{2}-U^{J}\left(\chi^{J}\right)\right]\,,$
(14)

where $\chi$ is the inflaton whose potential $U(\chi)$ determines the Hubble
parameter $H\sim\sqrt{G\,U(\chi)}$ during the inflation. The constraint above
on $H$ have been checked for the inflation in the Einstein frame (denote it by
the index $E$). The transition to it from (14) is known [36, 14]:

$g_{\mu\nu}^{J}=\Omega^{-2}g_{\mu\nu}^{E}\,,\quad\text{where}\quad\Omega^{2}=e^{n\beta(\phi)}|f^{\prime}(\phi)|.$
(15)

Substituting (15) into (14), we obtain the action written in Einstein frame:

$\displaystyle S^{J}_{\text{inf}}$ $\displaystyle=\int
d^{4}x\sqrt{-g_{4}^{E}}\,\left[\frac{1}{2}\Omega^{-2}\left(\partial\chi^{J}\right)^{2}-\Omega^{-4}U^{J}\left(\chi^{J}\right)\right]\approx$
$\displaystyle\approx\int
d^{4}x\sqrt{-g_{4}^{E}}\,\left[\frac{1}{2}\left(\partial\chi^{E}\right)^{2}-U^{E}\left(\chi^{E}\right)\right]=S^{E}_{\text{inf}}\,,$
(16)

where transformations of the inflaton $\chi^{E}=\Omega^{-1}\chi^{J}$ and its
potential $U^{E}=\Omega^{-4}U^{J}$ to the Einstein frame are made. Here we use
the fact that during inflation the field $\phi$ is practically frozen
(constraint 2) and is in the region of maximum of potential, so the factor
$\Omega(\phi)$ can be considered constant and equal
$\Omega\approx\Omega(\phi_{\text{max}})\sim\Omega(\phi_{\text{min}})\sim
10^{8}$ for the parameters chosen in Fig. 1.

The Hubble parameter for the inflaton $\chi^{E}$ observed in the Einstein
frame is related to the Hubble parameter for the $\chi^{J}$ inflaton observed
in the Jordan frame as follows:

$H^{E}\sim\sqrt{G\,U^{E}(\chi)}=\sqrt{G\
\Omega^{-4}U^{J}(\chi)}\sim\Omega^{-2}H^{J}\,.$ (17)

From this, in accordance with the constraint for $H$ obtained above for the
Einstein frame, it follows that $H^{J}>10^{5}\sim 10^{-4}\ [M_{4}]\sim
10^{14}\ [\text{GeV}]$ at the chosen parameters. It is consistent with the
observed data [13].

##  6 Formation of primordial black holes

During the inflationary stage the classical motion of the scalar field $\psi$
is frozen — this is determined by the constraint 2 in Section 5). If there a
seed for the future domain wall appear during quantum fluctuations (the field
jumps to the left slope of the potential in Figure 2), then in this place the
field values $\psi$ should lie near the potential maximum since the
fluctuations are small (constraint 4 in Section 5). At the end of inflation,
the inflaton the Hubble parameter vanishes and the $\psi$ field begins to
rapidly roll from the maximum to the potential minima (to the left minimum,
for the inner region of the bubble and to the right minimum for the outer
region). During this roll-off, the energy density $\sigma(t)$ in the
transition region gradually grows up to the value $\sigma$ calculated in the
previous sections. In addition, due to postinflationary expansion, the radius
of the region is growing $u_{w}(t)=u_{0}a(t)$ also, where $u_{0}$ is radius at
the end of inflation and $a(t)$ is scale factor. Because of this, the mass of
the formed domain wall increases $m_{w}(t)=4\pi u_{w}(t)^{2}\sigma(t)$
together with its gravitational radius $u_{g}(t)=2Gm_{w}(t)$. At a certain
time, the gravitational radius will cross the entire domain wall
$u_{g}(t_{*})=u_{w}(t_{*})$ and a primordial black hole with mass
$m_{w}(t_{*})$ will form for a remote observer.

The moment $t_{*}$ of crossing of the wall by the gravitational radius in our
model will come long before it reaches the final energy density $\sigma$. This
can be seen from the fact that the ratio of the gravitational radius to the
size of the wall (assuming that the wall has a final energy density $\sigma$)
is grater than unity

$\frac{u_{g}}{u_{w}}=8\pi G\sigma u_{w}>m_{4}^{-2}\sigma\delta\approx 16\,,$
(18)

where the fact that the wall radius is always larger than its thickness
$u_{w}(t)>\delta$ is used. The ratio (18) is calculated for the wall
parameters obtained in Fig. 4 and holds for any model parameters (1)
satisfying the set of constraints discussed in Section 5. Therefore, the wall
is formed long before the field rolls to the potential minima. The black hole
is formed at the time $t_{*}$ when the horizon radius $u_{g}$ is approximately
equal to the wall radius $u_{w}$. Hence, the mass of the PBH can be estimated
as

$M_{PBH}\simeq m_{w}(t_{*})\simeq u_{g}(t_{*})/2G\simeq u_{w}(t_{*})/2G\equiv
u_{0}a(t_{*})/2G\,.$ (19)

All processes described above occur in a very short time interval $t_{*}$,
which is significantly less than the characteristic time of the field $\psi$
roll-off to the minimum of potential,
$t_{*}<\tau_{\psi}\sim\sqrt{V^{\prime\prime}(\psi_{\text{max}})}^{-1}\sim
10^{6}$. This time in our model is less than the inflation time
$\tau_{\text{inf}}\approx 60\ H^{-1}\sim 10^{7}$. Therefore, the roll-off of
field $\psi$ process occurs at the reheating stage (as for the inflaton), the
scale factor at which we can approximately [37] consider
$a(t_{*})=\left((\tau_{\text{inf}}+t_{*})/\tau_{\text{inf}}\right)^{\frac{\nu+2}{3\nu}}<\left((\tau_{\text{inf}}+\tau_{\psi})/\tau_{\text{inf}}\right)^{\frac{\nu+2}{3\nu}}\sim\left(1+10^{-1}\right)^{\frac{\nu+2}{3\nu}}\sim
1$, where $\nu$ is the index of degree in the inflaton potential. It means
that the time interval $t_{*}$ is very short and one can neglect the post-
inflationary expansion.

The radius $u_{0}$ is determined by the moment of inflationary generation of
the fluctuation, leading to the subsequent formation of the wall. If the
fluctuation is formed at the e-fold $N$ of total number $N_{\text{inf}}$, then
its size by the end of inflation will be
$u_{0}=H_{\text{inf}}^{-1}e^{N_{\text{inf}}-N}$. For example, suppose that the
fluctuation of the $\psi$ field leading to the wall formation occurred at
$N=25$ e-fold of $N_{\text{inf}}=60$ and the inflation scale is $H\sim
10^{14}$ GeV (which are characteristic parameters [13]) Then, immediately
after the inflation ends, a PBH of mass $M_{PBH}\sim u_{w}/2G\sim
u_{0}/2G=4\pi M_{4}^{2}H_{\text{inf}}^{-1}e^{35}\sim 10^{41}\ [GeV]\sim
10^{-16}\ [M_{\odot}]$ is formed. During further cosmological evolution in the
reheating stage, such a PBH will gain mass as a result of accretion. The
dynamics of this process is rather complicated and depends on the properties
of the reheating stage [38]. Since such black holes are formed at the super-
dense Universe, immediately after the end of inflation, their masses can reach
many solar masses due to the accretion, as our estimates show.

Calculation of the initial mass spectrum $N(M)$ of the above described PBHs is
reduced to the calculation of the size spectrum $N(u)$ for the $\psi$ field
fluctuations generated at the inflationary stage. It can be done by
substituting the connection $u\equiv u_{w}=2GM$. The spectrum $N(u)$ is
calculated in such works as [28, 39, 40]. They investigated the dependence of
this spectrum on inflation parameters and the initial field value $\psi_{0}$,
with which the modern horizon appears. In these works it was shown that the
width of the spectrum, the characteristic masses, and the total number of PBHs
strongly depend on the choice of $\psi_{0}$. We will not investigate the
spectrum in the present paper, since its goal was to demonstrate a new
mechanism for PBH formation, arising in $f(R)$-gravity with the extra
dimensions.

##  7 Conclusion

In this paper we discuss the model of multidimensional $f(R)$-gravity that
gives rise the primordial black holes formation at the early stage of
cosmological evolution. For a remote observer in the Jordan frame, the
formation of such black holes can be interpreted as a manifestation of non-
trivial gravitational dynamics of multidimensional space after the
cosmological inflation.

In this paper the dynamics of the extra dimensional metric is reduced to the
dynamics of the effective scalar field. The domain walls formation occurs just
after the completion of cosmological inflation. We show that PBH nucleation is
inevitable in our model for that range of parameters of $f(R)$-gravity which
both satisfy cosmological constraints and lead to a non-trivial effective
scalar field potential. The mechanism of collapse into primordial black holes
has been studied in [28, 39, 40].

Primordial black holes arising in the developed mechanism appear at the
reheating stage and then will actively grow due to accretion. In addition, the
formation of several domain walls under one gravitational radius is possible,
which modifies our mechanism. These rather complex processes will affect the
formation of the final mass spectrum of primordial black holes and will be
considered by us in future works.

Cosmological consequences of multidimensional $f(R)$-gravity are quite rich
and lead to a variety of exotic observational manifestations. As our study
shows, the widely discussed primordial black holes are added to such
manifestations. Therefore, the confirmation of the primordial origin of some
classes of black holes may be the evidence in favor of existence of extra
dimensions.

##  8 Acknowledgements

The work of V.V.N. and M.A.K. was supported by the MEPhI Program Priority
2030. The work of S.G.R. has been supported by the Kazan Federal University
Strategic Academic Leadership Program.

## References

  * [1] Ferruccio Feruglio  “Extra dimensions in particle physics”  In _The European Physical Journal C_ 33.S1 Springer ScienceBusiness Media LLC, 2004, pp. s114–s128  DOI: 10.1140/epjcd/s2004-03-1699-8
  * [2] Keith R Dienes, Emilian Dudas and Tony Gherghetta  “Grand unification at intermediate mass scales through extra dimensions”  In _Nuclear Physics B_ 537.1-3 Elsevier, 1999, pp. 47–108  DOI: 10.1016/S0550-3213(98)00669-5
  * [3] Lawrence Hall and Yasunori Nomura  “Gauge unification in higher dimensions”  In _Physical Review D_ 64.5 APS, 2001, pp. 055003  DOI: 10.1103/PhysRevD.64.055003
  * [4] A.. Grobov and S.. Rubin  “Higgs-Like Field and Extra Dimensions”  In _International Journal of Theoretical Physics_ 52.12 Springer ScienceBusiness Media LLC, 2013, pp. 4283–4292  DOI: 10.1007/s10773-013-1742-9
  * [5] JoAnne Hewett and Maria Spiropulu  “Particle physics probes of extra spacetime dimensions”  In _Annual Review of Nuclear and Particle Science_ 52.1, 2002, pp. 397–424  DOI: 10.1146/annurev.nucl.52.050102.090706
  * [6] Nicolas Deutschmann, Thomas Flacke and Jong Soo Kim  “Current LHC constraints on minimal universal extra dimensions”  In _Physics Letters B_ 771 North-Holland, 2017, pp. 515–520  DOI: 10.1016/J.PHYSLETB.2017.06.004
  * [7] Edward Witten  “Instability of the Kaluza-Klein vacuum”  In _Nuclear Physics B_ 195.3 Elsevier, 1982, pp. 481–492  DOI: 10.1016/0550-3213(82)90007-4
  * [8] Sean M. Carroll, James Geddes, Mark B. Hoffman and Robert M. Wald  “Classical stabilization of homogeneous extra dimensions”  In _Physical Review D_ 66.2 American Physical Society (APS), 2002  DOI: 10.1103/physrevd.66.024036
  * [9] Tonguç Rador  “Acceleration of the Universe via f(R) gravities and the stability of extra dimensions”  In _Physical Review D_ 75.6 American Physical Society (APS), 2007  DOI: 10.1103/physrevd.75.064033
  * [10] K.. Bronnikov and S.. Rubin  “Self-stabilization of extra dimensions”  In _Phys. Rev. D_ 73.12, 2006, pp. 124019  DOI: 10.1103/PhysRevD.73.124019
  * [11] Alexei A. Starobinsky  “A New Type of Isotropic Cosmological Models Without Singularity”  In _Phys. Lett._ B91, 1980, pp. 99–102  DOI: 10.1016/0370-2693(80)90670-X
  * [12] Alexander Vilenkin  “Classical and quantum cosmology of the Starobinsky inflationary model”  In _Physical Review D_ 32.10 APS, 1985, pp. 2511  DOI: 10.1103/PhysRevD.32.2511
  * [13] Y. Akrami et al.  “Planck2018 results”  In _Astronomy & Astrophysics_ 641 EDP Sciences, 2020, pp. A10  DOI: 10.1051/0004-6361/201833887
  * [14] Júlio C Fabris, Arkady A Popov and Sergey G Rubin  “Multidimensional gravity with higher derivatives and inflation”  In _Physics Letters B_ 806 Elsevier, 2020, pp. 135458  DOI: 10.1016/j.physletb.2020.135458
  * [15] Antonio De Felice and Shinji Tsujikawa  “f(R) Theories”  In _Living Reviews in Relativity_ 13.1 Springer ScienceBusiness Media LLC, 2010  DOI: 10.12942/lrr-2010-3
  * [16] Salvatore Capozziello and Mariafelicia De Laurentis  “Extended Theories of Gravity”  In _Physics Reports_ 509.4, 2011, pp. 167–321  DOI: 10.1016/j.physrep.2011.09.003
  * [17] S. Capozziello and M. Laurentis  “The Dark Matter problem from f(R) gravity viewpoint”  In _Annalen der Physik_ 524, 2012, pp. 545–578  DOI: 10.1002/andp.201200109
  * [18] K A Bronnikov, R V Konoplich and S G Rubin  “The diversity of universes created by pure gravity”  In _Classical and Quantum Gravity_ 24.5 IOP Publishing, 2007, pp. 1261–1277  DOI: 10.1088/0264-9381/24/5/011
  * [19] K.A. Bronnikov et al.  “Inhomogeneous compact extra dimensions”  In _Journal of Cosmology and Astroparticle Physics_ 2017.10 IOP Publishing, 2017, pp. 001–001  DOI: 10.1088/1475-7516/2017/10/001
  * [20] Kirill A. Bronnikov, Arkady A. Popov and Sergey G. Rubin  “Inhomogeneous compact extra dimensions and de Sitter cosmology”  In _The European Physical Journal C_ 80.10 Springer ScienceBusiness Media LLC, 2020  DOI: 10.1140/epjc/s10052-020-08547-x
  * [21] R. Falomo et al.  “Low-redshift quasars in the Sloan Digital Sky Survey Stripe 82. The host galaxies”  In _Monthly Notices of the Royal Astronomical Society_ 440.1, 2014, pp. 476–493  DOI: 10.1093/mnras/stu283
  * [22] VI Dokuchaev, Yu N Eroshenko and SG Rubin  “Origin of supermassive black holes”  In _arXiv preprint arXiv:0709.0070_ , 2007  DOI: 10.48550/arXiv.0709.0070
  * [23] The Virgo Collaboration. The LIGO Scientific Collaboration  “Search for intermediate mass black hole binaries in the first and second observing runs of the Advanced LIGO and Virgo network”  In _Physical Review D_ 100.6 American Physical Society (APS), 2019  DOI: 10.1103/physrevd.100.064064
  * [24] Bernard Carr and Florian Kuhnel  “Primordial Black Holes as Dark Matter Candidates”  In _arXiv preprint arXiv:2110.02821_ , 2021  DOI: 10.48550/arXiv.2110.02821
  * [25] V. Luca, G. Franciolini, P. Pani and A. Riotto  “Primordial black holes confront LIGO/Virgo data: current situation”  In _Journal of Cosmology and Astroparticle Physics_ 2020.06 IOP Publishing, 2020, pp. 044–044  DOI: 10.1088/1475-7516/2020/06/044
  * [26] Alexander S. Sakharov, Yury N. Eroshenko and Sergey G. Rubin  “Looking at the NANOGrav signal through the anthropic window of axionlike particles”  In _Phys. Rev. D_ 104.4, 2021, pp. 043005  DOI: 10.1103/PhysRevD.104.043005
  * [27] Sergey Rubin, M. Khlopov and A. Sakharov  “Primordial Black Holes from Non-Equilibrium Second Order Phase Transition”  In _Grav. Cosmol._ 6, 2000  DOI: 10.48550/arXiv.hep-ph/0005271
  * [28] Konstantin M. Belotsky et al.  “Clusters of Primordial Black Holes”  In _The European Physical Journal C_ 79.3 Springer ScienceBusiness Media LLC, 2019  DOI: 10.1140/epjc/s10052-019-6741-4
  * [29] Yana Lyakhova, Arkady A Popov and Sergey G Rubin  “Classical evolution of subspaces”  In _The European Physical Journal C_ 78.9 Springer, 2018, pp. 1–13  DOI: 10.1140/epjc/s10052-018-6251-9
  * [30] Kirill A Bronnikov and Sergey G Rubin  “Local regions with expanding extra dimensions”  In _Physics_ 3.3 Multidisciplinary Digital Publishing Institute, 2021, pp. 781–789  DOI: 10.3390/physics3030048
  * [31] S. Capozziello, P. Martin-Moruno and C. Rubano  “Physical non-equivalence of the Jordan and Einstein frames”  In _Physics Letters B_ 689.4, 2010, pp. 117–121  DOI: 10.1016/j.physletb.2010.04.058
  * [32] “Black Holes, Cosmology and Extra Dimensions”  In _Black Holes, Cosmology and Extra Dimensions. Edited by Bronnikov Kirill A et al. Published by World Scientific Publishing Co. Pte. Ltd., 2013. ISBN #9789814374217_ , 2013  DOI: 10.1142/9789814374217
  * [33] Alexander Vilenkin and E Paul S Shellard  “Cosmic strings and other topological defects”  Cambridge University Press, 1994
  * [34] Jing Liu, Zong-Kuan Guo and Rong-Gen Cai  “Primordial black holes from cosmic domain walls”  In _Physical Review D_ 101.2 American Physical Society (APS), 2020  DOI: 10.1103/physrevd.101.023513
  * [35] Jaume Garriga, Alexander Vilenkin and Jun Zhang  “Black holes and the multiverse”  In _Journal of Cosmology and Astroparticle Physics_ 2016.02 IOP Publishing, 2016, pp. 064–064  DOI: 10.1088/1475-7516/2016/02/064
  * [36] Guillem Domènech and Misao Sasaki  “Conformal frame dependence of inflation”  In _Journal of Cosmology and Astroparticle Physics_ 2015.04 IOP Publishing, 2015, pp. 022–022  DOI: 10.1088/1475-7516/2015/04/022
  * [37] George Lazarides  “Inflationary cosmology”  In _Cosmological Crossroads_ Springer, 2002, pp. 351–391  DOI: 10.48550/arXiv.hep-ph/0111328
  * [38] Bernard Carr, Konstantinos Dimopoulos, Charlotte Owen and Tommi Tenkanen  “Primordial black hole formation during slow reheating after inflation”  In _Physical Review D_ 97.12 APS, 2018, pp. 123535  DOI: 10.1103/PhysRevD.97.123535
  * [39] VV Nikulin, AV Grobov and SG Rubin  “A mechanism for protogalaxies nuclei formation from primordial black holes clusters”  In _Journal of Physics: Conference Series_ 934, 2017, pp. 012040  IOP Publishing  DOI: 10.1088/1742-6596/934/1/012040
  * [40] V. Nikulin, Sergey Rubin and L. Khromykh  “Formation of Primordial Black Hole Clusters from Phase Transitions in the Early Universe”  In _Bulletin of the Lebedev Physics Institute_ 46, 2019, pp. 97–99  DOI: 10.3103/S1068335619030060