# Collective steady-state patterns of swarmalators with finite-cutoff
interaction distance

Hyun Keun Lee  Department of Physics, Sungkyunkwan University, Suwon 16419,
Korea  Kangmo Yeo  Department of Physics, Jeonbuk National University, Jeonju
54896, Korea  Hyunsuk Hong<EMAIL_ADDRESS>Department of Physics, Jeonbuk
National University, Jeonju 54896, Korea  Research Institute of Physics and
Chemistry, Jeonbuk National University, Jeonju 54896, Korea

###### Abstract

We study the steady-state patterns of population of the coupled oscillators
that sync and swarm, where the interaction distances among oscillators have
finite-cutoff in interaction distance. We examine how the static patterns
known in the infinite-cutoff are reproduced or deformed, and explore a new
static pattern that does not appear until a finite-cutoff is considered. All
steady-state patterns of the infinite-cutoff, static sync, static async, and
static phase wave are respectively repeated in space for proper finite-cutoff
ranges. Their deformation in shape and density takes place for the other
finite-cutoff ranges. Bar-like phase wave states are observed, which has not
been the case for the infinite-cutoff. All the patterns are investigated via
numerical and theoretical analysis.

> Swarming [1] is the collective behavior of biological/artificial entities in
> the absence of centralized coordination. Interaction between agents is the
> reason for the behavior, and thus its modeling is a core task to understand
> or control the emergent phenomena. Recently, a steady state pattern was
> studied [2] for the overdamped limit model with potential while motivated by
> biological aggregations in viscous liquid, and the static feature later
> becomes various [3] as self-propelling factor is implemented through the
> phase dynamics of synchronization [4]. The interaction therein is mediated
> by a potential function that may change in time. We here consider the cutoff
> of interaction distance to model the fact that the distance any
> living/engineering entities can manage is necessarily finite. The previously
> understood static feature lasts till the cutoff is larger than the size of
> the pattern formed in the absence of cutoff. For the smaller cutoff, the
> pattern is deformed in an explainable manner, and then it becomes amorphous
> or non-stationary for further decrease. A new class of static structures not
> appearing in the no-cutoff model is observed and specified. Phase diagram
> for the steady-state patterns according to the cutoff distance is
> illustrated. Finite-cutoff of interaction distance makes the swarm behavior
> prolific as combined with underlying potential.

##  I introduction

Synchronization in the population of coupled oscillators has been widely
explored. In particular, not only mathematical models but also various
experimental systems have been considered to investigate the synchronization
behavior [5, 6, 7, 8, 9, 10, 11, 4, 12, 13, 14, 15]. In one recent study [3],
the oscillators that also swarm [2, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]
was introduced while called swarmalator, and the self-organization in space
and time was reported. Three steady states and two nonstationary ones have
been found in the swarmalators model with long-range attraction and repulsion.
The swarmalator model is followed by various researches [26, 27, 28, 29, 30,
31, 32, 33] in its novel modeling of collective behaviors and practical
applicability as well.

According to Refs. [2, 3], the spatial sizes of the steady-state patterns are
finite when the interaction distance is long-ranged of no cutoff or infinite-
cutoff. Here, one may expect there can be a finite-cutoff in interaction
range, for which the result is qualitatively same in the sense that all the
steady-state patterns by the long range interaction are reproduced. If a
specific pattern known for the infinite-cutoff is the target of biological or
engineering swarms, a proper interaction-range cutoff larger than the
pattern’s size could be enough. If an interaction range is smaller than this
cutoff but still comparable, a describable change of the pattern is expected.
Moreover, since the realization of infinite range interaction is by itself
impractical, the system with finite-cutoff has actual implications. The
various spatiotemoral constraints [29, 30] are indeed the factual barriers in
realizing the notion of swarmalators. All of these motivate the current work.

In this paper, we consider the local interaction among the swarmalators, where
the interaction range is restricted by a control parameter, and explore the
collective steady states of the system. We examine how the steady-state
patterns of the system with long-range interaction are deformed in the
understandable way and also pay special attention to the possibility of new
steady states. This paper consists of five sections. Section II introduces the
model, and Sec. III revisits the system with infinite-cutoff (global
interaction). In Sec. IV, the finite-cutoff interaction is considered, and the
various steady-state patterns are demonstrated. Section V gives a brief
summary.

##  II model

We consider a population of oscillators governed by

$\displaystyle\dot{\bf r}_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{i}(r_{\rm c})}\sum_{j\in\Lambda_{i}(r_{\rm
c})}\left[\Big{(}1+J\cos(\theta_{j}-\theta_{i})\Big{)}-\frac{1}{r^{2}_{ij}}\right]$
(1)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\Big{(}{\bf
r}_{j}-{\bf r}_{i}\Big{)},$ $\displaystyle\dot{\theta}_{i}$ $\displaystyle=$
$\displaystyle\omega_{i}+\frac{K}{N}\sum_{j\neq
i}\frac{\sin(\theta_{j}-\theta_{i})}{r_{ij}},$ (2)

where ${\bf r}_{i}={\bf r}_{i}(t)=(x_{i}(t),y_{i}(t))$ is the position vector
of agent $i$ at time $t$ in 2-dimensional space with $r_{ij}\equiv|{\bf
r}_{i}-{\bf r}_{j}|$ and $\theta_{i}=\theta_{i}(t)$ is the $i$’s phase angle.
Note that the dynamics governed by Eqs. (1) and (2) makes the oscillators
swarm and sync. Considering such characteristic we here call the oscillators
swarmalators, following the Ref. [3]. $N_{i}(r_{\rm c})$ is the number of
members in the collection $\Lambda_{i}(r_{\rm c})$ composed of such
swarmalators whose distances to $i$ are not greater than $r_{c}$ (see Fig. 1).
This way, the interaction range of spatial dynamics is bounded by $r_{\rm c}$.

Figure 1: (Color Online) Schematic diagram for the one swarmalator’s neighbors
controlled by the range parameter $r_{c}$. Here the dots represent the
swarmalators and the colors their phases $\theta(\in[-\pi,\pi))$.

$J$ is the parameter for attraction with $|J|\leq 1$, for which $J=1$ is used
in the numerical work below. $\omega_{i}$ is the natural frequency of $i$, $K$
is a coupling strength of phase and $N$ is the total number of swarmalators in
system. We consider identical swarmalators with $\omega_{i}=0$ for all $i$ in
this work.

In the $r_{\rm c}\to\infty$ limit, when $J=0$, Eq. (1) recovers the model
suggested in Ref. [2] of special interest in the formation of uniform density
and identical radius disc of swarms. The motion of ${\bf r}_{i}$ for $J=0$ is
basically a gradient flow as each summand can be rewritten by the gradient of
a scalar function for $r_{ij}$. As Eq. (1) with $J=0$ is motivated by the
biological aggregations in viscous liquid [21, 25], the inertia effect is
ignored to set up an overdamped limit equation [34]. Next when $J\neq 0$, Eqs.
(1) and (2) corresponds to one of the models in Ref. [3], which share common
feature of various active and static phases of the swarmalators. Due to the
phase variables in Eq. (1), which change in time by Eq. (2), the motion of
${\bf r}_{i}$ now has self-propelling property [17, 18, 19]. Equation (2) is a
variant of the Kuramoto model [6, 7], i.e., the coupling strength $K$ in this
prototype model is generalized to $K_{ij}=K/r_{ij}$ in Eq. (2) in order to
reflect the spatial configuration in phase dynamics. This was suggested in
Ref. [3] to generate pairing current of the swarmalators in a proper range of
negative $K$ with $r_{\rm c}=\infty$.

##  III steady states for infinite cutoff

According to Refs. [2, 3], there are three steady-state patterns when $r_{\rm
c}=\infty$. As we are interested in their reproduction or moderate change for
finite $r_{\rm c}$, we begin with explaining the three patterns, in a
comprehensive way as possible. In the steady states, $\dot{r}_{i}=0$ and
$\dot{\theta}_{i}=0$ for all $i$. Therefore, as its continuum counterpart such
as velocity field ${\bf v}({\bf r})$ at position ${\bf r}$, it is reasonable
to set ${\bf v}({\bf r})=0$. When one considers a proper coarse-grain of Eq.
(1) in the steady states, the velocity field reads

${\bf v}({\bf r})=\int d{\bf r}^{\prime}\left(-\nabla_{{\bf r}}U({\bf r},{\bf
r}^{\prime})\right)\rho({\bf r}^{\prime})=0,$ (3)

where

$\displaystyle U({\bf r},{\bf r}^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{2}|{\bf r}-{\bf r}^{\prime}|^{2}(1+J\cos(\theta({\bf
r})-\theta({\bf r}^{\prime})))-\ln|{\bf r}-{\bf r}^{\prime}|$ (4)
$\displaystyle\equiv$ $\displaystyle U_{\rm att}({\bf r},{\bf
r}^{\prime})+U_{\rm rep}({\bf r},{\bf r}^{\prime})$

with phase field $\theta({\bf r})$, and $\rho({\bf r})$ is the steady-state
density of swarmalators at position ${\bf r}$ with normalization $\int d{\bf
r}\rho({\bf r})=1$. As Eq. (3) holds for all ${\bf r}$, it also holds that
$\nabla_{\bf r}\cdot{\bf v}({\bf r})=0$. Then, in the fact that $\nabla_{\bf
r}^{2}U_{\rm rep}({\bf r},{\bf r}^{\prime})=2\pi\delta({\bf r}-{\bf
r}^{\prime})$ [35], one finds a self-consistent equation for $\rho({\bf r})$,
which is given by

$\displaystyle\rho({\bf r})$ $\displaystyle=$ $\displaystyle\frac{1}{2\pi}\int
d{\bf r}^{\prime}\nabla_{\bf r}^{2}U_{\rm att}({\bf r},{\bf
r}^{\prime})\rho({\bf r}^{\prime}).$ (5)

This self-consistent equation plays a central role in finding the steady-state
patterns in the help of other information like symmetry and/or numerical
observations.

One of the three steady-state patterns appears when $K>0$. For positive $K$,
Eq. (2) derives the system to phase synchronization state of a common phase
$\theta_{\rm s}$. Then, Eq. (1) becomes free from phase as
$\cos(\theta_{j}-\theta_{i})\to 1$ as evolution goes on. In the steady state,
therefore, it reads

$U_{\rm att}({\bf r},{\bf r}^{\prime})=\frac{1}{2}|{\bf r}-{\bf
r}^{\prime}|^{2}(1+J).$ (6)

When this is considered in Eq. (5) in the two-dimensional $(x,y)$ coordinate,
since $\nabla_{\bf r}^{2}U_{\rm att}({\bf r},{\bf r}^{\prime})=2(1+J)$ is
constant, using $\int d{\bf r}^{\prime}\rho({\bf r}^{\prime})=1$, $\rho({\bf
r})=(1+J)/\pi$ is immediate. Next, it should be answered where to reside.

The center of positions is conserved by the symmetric pairs of changes in Eq.
(1). Thus when the center is used as the origin, it is natural to expect the
rotational symmetry of density, i.e., $\rho({\bf r})=\rho(r)$ for $r=|{\bf
r}|$. The numerical results always give such distributions that strongly
suggest a circular $\rho({\bf r})$ as shown in Fig. 2(a).

Figure 2: (Color Online) Three steady-state patterns known in the $r_{\rm
c}\to\infty$ limit [3]: (a) is the static sync (SS) when $K>0$. (b) is the
static async (SA) when $K<K_{\rm c}$ for a negative critical value $K_{\rm
c}$. (c) is the static phase wave (SPW) when $K=0$. For numerical realization,
$N=200$ system with $J=1$ is considered. $K=1$ and $K=-2$ are used for (a) and
(b), respectively. The black circles therein guide the associated radius
analytically understood (see text).

Thus when this is considered in Eq. (5), it follows that

$\rho({\bf r})=\frac{1}{\pi R_{\rm s}^{2}}~{}~{}{\rm for}~{}~{}r<R_{\rm s},$
(7)

where $R_{\rm s}$ is given by

$R_{\rm s}=1/\sqrt{1+J}.$ (8)

Interestingly, the density is uniform and the radius depends only on $J$. The
pattern by Eqs. (7) and (8) was named the static sync (SS) [3].

The steady-state condition [Eq. (3)] holds for Eq. (7) as follows. When

$V({\bf r})=\int d{\bf r}^{\prime}U({\bf r},{\bf r}^{\prime})\rho({\bf
r}^{\prime})$ (9)

is introduced, $\nabla_{\bf r}\cdot{\bf v}({\bf r})=0$ corresponds to
Laplace’s equation $\nabla_{\bf r}^{2}V({\bf r})=0$. Thus, $\rho({\bf r})$ of
Eq. (5) gives $V({\bf r})$ of the Laplace’s equation, and this is mediated by
Eq. (9). As $V({\bf r})=V(r)$ by $U_{\rm att}({\bf r},{\bf r}^{\prime})=U_{\rm
att}(|{\bf r}-{\bf r}^{\prime}|)$ in Eq. (6) and $\rho({\bf r})=\rho(r)$ in
Eq. (7), the solution of the Laplace’s equation is given by such cylindrical
harmonics of no longitudinal neither angular dependence. Then, the choice
should be [35] $V({\bf r})\propto 1$ exclusively or $\ln r$. Here, one
excludes the latter as $V({\bf r})$ remains finite at ${\bf r}=0$ in a simple
calculation of Eq. (9) with Eqs. (6) and (7). Thus, ${\bf v}({\bf
r})=-\nabla_{\bf r}\cdot V({\bf r})=0$ is immediate for constant $V({\bf r})$.

The second steady-state pattern appears [3] when $K<K_{\rm c}<0$ holds for a
negative critical value $K_{\rm c}\approx-1.2J$. For negative $K$, uniform
$\theta_{i}$ is unstable in Eq. (2), and an erratic profile (asynchrony) would
appear. Below $K_{\rm c}$, if the configuration of $\theta_{i}$ is erratic
enough, the summation of the sinusoidal part of Eq. (1) may vanish. When this
is the case, $U_{\rm att}({\bf r},{\bf r}^{\prime})=(1/2)|{\bf r}-{\bf
r}^{\prime}|^{2}$ reads. This is no more than Eq. (6) with $J=0$, from which
it follows that

$\rho({\bf r})=\frac{1}{\pi}~{}~{}{\rm for}~{}~{}r<R_{\rm a}=1.$ (10)

This result of uniform density and unit radius successively explains the
numerical data shown in Fig. 2(b). This pattern was named the static async
(SA) [3].

The third steady-state pattern appears [3] when $K=0$. For $K=0$, one
trivially reads $\dot{\theta}_{i}=0$ in Eq. (2). The numerical data shown in
Fig. 2(c) indicates $\theta_{i}$ can be replaced with the spatial angle
$\phi_{i}$ of ${\bf r}_{i}$; $\theta_{i}=\phi_{i}+c$ for a constant $c$. Then,
in polar coordinate ${\bf r}=(r,\phi)$, it reads

$\displaystyle U_{\rm att}({\bf r},{\bf r}^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(r^{2}+r^{\prime
2}-2rr^{\prime}\cos(\phi-\phi^{\prime})\right)$ (11)
$\displaystyle~{}~{}\times\left(1+J\cos(\phi-\phi^{\prime})\right).$

Using Eq. (11), we find its Laplacian is given by

$\displaystyle\nabla_{\bf r}^{2}U_{\rm att}({\bf r},{\bf r}^{\prime})$
$\displaystyle=$ $\displaystyle
2-\frac{Jr^{\prime}}{2r}+f_{1}(r,r^{\prime})\cos(\phi-\phi^{\prime})$ (12)
$\displaystyle~{}~{}~{}~{}~{}~{}+f_{2}(r,r^{\prime})\cos
2(\phi-\phi^{\prime})$

for proper functions $f_{1}(r,r^{\prime})$ and $f_{2}(r,r^{\prime})$. Also,
the numerically observed circular strip suggests that $\rho({\bf r})=\rho(r)$.
Plugging this and Eq. (12) into Eq. (5), one finds

$\rho({\bf r})=\frac{1}{\pi}-\frac{C}{r}~{}~{}{\rm for}~{}~{}R_{1}<r<R_{2}$
(13)

for constant $C$. One again considers Eq. (13) in Eq. (5) to know

$C=\frac{2J(R_{2}^{3}-R_{1}^{3})}{3\pi(4+J(R_{2}^{2}-R_{1}^{2}))}.$ (14)

Normalization of $\rho({\bf r})$ gives

$(2\pi C-R_{1}-R_{2})(R_{1}-R_{2})=1.$ (15)

To fix the solution, one may refer to Eq. (3) of ${\bf v}({\bf r})=0$. This
requires vanishing radial component because ${\bf v}({\bf r})=(v_{\rm att}(r)\
+v_{\rm rep}(r))\hat{\bf r}$ in $U_{\rm att/rep}({\bf r},{\bf
r}^{\prime})=U_{\rm att/rep}(|{\bf r}-{\bf r}^{\prime}|)$ and $\rho({\bf
r})=\rho(r)$, where $\hat{\bf r}\equiv{\bf r}/r$ is the unit vector in radial
direction. Attraction part is

$\displaystyle v_{\rm att}(r)$ $\displaystyle=$ $\displaystyle\int d{\bf
r}^{\prime}(-\partial_{r}U_{\rm att}({\bf r},{\bf
r}^{\prime}))\rho(r^{\prime})$ $\displaystyle=$
$\displaystyle(R_{2}^{2}-R_{1}^{2}-2\pi C(R_{2}-R_{1}))r-2\pi C.$

For repulsion part, one can use the divergence theorem [35] and $\nabla_{{\bf
r}^{\prime\prime}}^{2}U_{\rm rep}({\bf r}^{\prime\prime},{\bf
r}^{\prime})=2\pi\delta({\bf r}^{\prime\prime}-{\bf r}^{\prime})$ to write

$\displaystyle 2\pi rv_{\rm rep}(r)$ $\displaystyle=$
$\displaystyle-\int_{R_{1}\leq r^{\prime\prime}\leq r}d{\bf
r}^{\prime\prime}\int d{\bf r}^{\prime}\nabla_{{\bf
r}^{\prime\prime}}^{2}U_{\rm rep}({\bf r}^{\prime\prime},{\bf
r}^{\prime})\rho(r^{\prime})$ (17) $\displaystyle=$ $\displaystyle
2\pi\left(2\pi C(r-R_{1})-(r^{2}-R_{1}^{2})\right).$

With Eqs. (III) and (17), one knows $v_{\rm att}(r)+v_{\rm rep}(r)=0$ is
possible only when

$2\pi C=R_{1}.$ (18)

Using Eq. (18), one obtains in Eqs. (14) and (15) that

$\displaystyle R_{1}$ $\displaystyle=$
$\displaystyle\Lambda_{J}((12-J)\sqrt{3}-3\sqrt{(4+J)(12-5J)})/12J$
$\displaystyle R_{2}$ $\displaystyle=$ $\displaystyle\Lambda_{J}/2\sqrt{3}$
(19)

for $\Lambda_{J}\equiv((12-3J+\sqrt{3(4+J)(12-5J)})/(2-J))^{1/2}$. Then, Eqs.
(18) and (III) fixes the solution in Eq. (13). This was named the static phase
wave (SPW) [3].

##  IV Finite-cutoff of interaction distance

In the numerical study with finite $r_{\rm c}$, we have observed various
steady-state patterns. Some of them are same as those in the $r_{\rm
c}\to\infty$ limit, another ones are their proper deformations, and the others
are classified as a new class. The three classes appear in order as $r_{\rm
c}$ decreases, and the shape becomes amorphous or non-stationary for further
decrease. In the following, the observed patterns are understood and/or
described by exploiting Eq. (5), the self-consistent equation for density.

###  IV.1 SS for $r_{\rm c}>2R_{\rm s}$

For finite $r_{\rm c}>2R_{\rm s}$, the system with positive coupling strength
$K$ is found to be in SS. It is interesting to note that the steady state with
multiple discs in SS is discovered (see Fig. 3). Multiple discs of SS have not
been found in the system with $r_{\rm c}=\infty$. Figure 3 is the numerical
results of $N=800$ system for $r_{\rm c}=1.5$ when $K=0.1$ and $J=1$ (we below
fix $J=1$), displaying the process how the multiple discs of SS are formed.

Figure 3: (Color Online) Formation of sync state pattern for $r_{\rm c}=1.5$
in $N=800$ system with $K=0.1$ and $J=1$ in the long time limit. Each panel is
a snapshot at the specified time $T$ at left upper corner. Multiple SS discs
are finally formed as shown in the last panel. For initial condition,
uniformly random positions sampled out of $L\times L$ square at origin was
used for $L=2$ (this type of initialization is applied below in this work).

The radius of the discs is, in common, $\approx 1/\sqrt{2}$ that reminds us of
$R_{\rm s}$ in Eq. (8) with $J=1$. A simple explanation of this numerical
observation is following; the swarmalators are divided into three groups apart
from others farther than $r_{\rm c}$ during transient period, and each group
evolves into its own steady state. The snapshots are shown in Fig. 3, from
initial $T=0$ via transient $T=50,~{}100$ to final $T=1000$.

Once a number of swarms is included in a circular region of radius $r_{\rm c}$
and this region is separated from the others farther than $r_{\rm c}$, the
dynamics taking place thereafter in the region is identical to that by $r_{\rm
c}=\infty$. If there are several regions separated that way (the number of
regions depends on initial condition), the same number of patterns finally
appears. We note, when $r_{\rm c}>2R_{\rm s}$, the self-consistent equation
for density [Eq. (5)] does not change from that for $r_{\rm c}=\infty$. Then,
each final pattern becomes no more than the SS of Eqs. (7) and (8). As
expected, the discs in the last snapshot are, in common, of radius
$1/\sqrt{2}$ with their own uniform densities. The minimal separation between
the discs is $2R_{\rm s}$ as a consequence of the condition $r_{\rm c}>2R_{\rm
s}$.

With no loss of generality, the condition $r_{\rm c}>2R_{\rm s}$ considered
for SS case is generalized to

$r_{\rm c}>D_{\infty}~{},$ (20)

where $D_{\infty}$ is the diameter of a pattern appearing in the $r_{\rm
c}\to\infty$ limit. It can be SA or SPW beside SS depending on $K$. That is,
any of them is expected in the same reason when Eq. (20) holds and, when
repeated, the minimal separation is given by the diameter.

###  IV.2 SA for $r_{\rm c}>2R_{\rm a}$

In case of static async, $D_{\infty}=2R_{\rm a}$. Thus by Eq. (20), SA should
appear for finite $r_{\rm c}>2R_{\rm a}$ when the negatively strong coupling
$(K<K_{c}<0)$ is provided. As expected, the system is found to be in SA where
the phases of all the swarmalators are random in the full range of
$[-\pi,\pi)$. Also, similar to the previous sync case, the steady state
composed of multiple discs of SA is discovered, which has not been found in
the system with $r_{\rm{c}}=\infty$. Figure 4 shows multiple SA discs for
$r_{\rm c}=2.1$ (note $D_{\infty}=2R_{\rm a}=2$ for SA).

Figure 4: (Color Online) Multiple SA discs in $N=400$ system with $r_{\rm
c}=2.1$ when $K=-8$ and $J=1$. For initial condition, that for Fig. 3 was used
with $L=10$.

###  IV.3 Anomalous SS/SA with nonuniform radial density

Next, we find there appear anomalous steady-state patterns for $r_{\rm
c}<2R_{\rm s}$ when $K>0$, as shown in Fig. 5(a) and (b).

Figure 5: (Color Online) Anomalous static states of $N=200$ system and
supporting schematics: Anomalous static sync (aSS) in (a) and (b) are,
respectively, the case when $r_{\rm c}=1.3$ and $r_{\rm c}=1.2$ for $K=0.1$
and $J=1$. Anomalous static async (aSA) in (c) is when $r_{\rm c}=1.7$ for
$K=-2$ and $J=1$. The geometric condition required for these patterns are
depicted in (d). The red arrow therein guides to the eyes, showing the
relation given by Eq. (23). Schematic of radial density is illustrated in (e).
In (d) and (e), $R_{\rm s}$ and $R_{\rm s}^{\prime}$ are replaced with $R_{\rm
a}$ and $R_{\rm a}^{\prime}$, respectively, in the async case.

Because Eq. (20) does not hold this time for $D_{\infty}=2R_{\rm s}$, what a
different pattern may form was in question, and the shown anomalous disc is
the result to be understood. The radius becomes larger to $R_{\rm s}+\Delta$.
The density looks uniform inside the region of radius $R_{\rm
s}^{\prime}(<R_{\rm s})$ while decreases thereafter. These observations can be
explained with Eq. (5) as follows.

For finite $r_{\rm c}$, it is given from Eq. (6) that

$\nabla_{\bf r}^{2}U_{\rm att}({\bf r},{\bf r}^{\prime})=2(1+J)H(r_{\rm
c}-|{\bf r}-{\bf r}^{\prime}|),$ (21)

where $H(x)$ is the Heaviside step-function assigned with $1$ for $x\geq 0$,
or $0$ otherwise. Plugging Eq. (21) into Eq. (5) with the numerical
observation $\rho({\bf r})=\rho(r)$, one may consider the geometric situations
by $R_{\rm s}+\Delta$, $R_{\rm s}^{\prime}$, and $r_{\rm c}$. Figure 5(d)
shows a situation when the interaction range of radius $r_{\rm c}$ marginally
covers the whole distribution of radius $R_{\rm s}+\Delta$; see the blue arc
and the dashed circle tangentially meet. Then, since $\nabla_{\bf r}^{2}U_{\rm
att}({\bf r},{\bf r}^{\prime})$ is constant in the interaction range, Eq. (5)
gives $\rho(r=R_{\rm s}^{\prime})=(2\pi)^{-1}2(1+J)\int d{\bf
r}^{\prime}H(r_{\rm c}-|{\bf r}-{\bf
r}^{\prime}|)\rho(r^{\prime})=((1+J)/\pi)\int d{\bf
r}^{\prime}\rho(r^{\prime})=(1+J)/\pi$. This value does not change as long as
the blue circle covers the dashed one, and this is the case when $r\leq R_{\rm
s}^{\prime}$. Note $(1+J)/\pi$ is the very $(\pi R_{\rm s}^{2})^{-1}$ of Eq.
(7).

When $r>R_{\rm s}^{\prime}$, the blue circle cannot include the dashed one,
and the covered region decreases as $r$ increases. So that, in Eq. (5), the
integration range becomes smaller for $r$, and this leads to decreasing
$\rho(r)$ in the left hand side. Then, one may write

$\rho({\bf r})=\begin{cases}\rho_{0}&{\rm for}~{}r\leq R_{\rm s}^{\prime}\\\
g(r)&{\rm for}~{}R_{\rm s}^{\prime}<r<R_{\rm s}+\Delta~{},\end{cases}$ (22)

where $\rho_{0}\equiv(\pi R_{\rm s}^{2})^{-1}$ and $g(r)$ is a proper function
that decreases for $r$ from $g(R_{\rm s}^{\prime})=\rho_{0}$. Figure 5(e)
shows the schematic of Eq. (22). The density begins to decrease from
$\rho_{0}$ following $O(\epsilon^{3/2})$ for small $\epsilon=r-R_{\rm
s}^{\prime}>0$ as the area $\sim\epsilon^{3/2}$ is excluded from the
integration in Eq. (5). A lower bound of $g(R_{\rm s}+\Delta)$ is simply
arguable in the $R_{\rm s}^{\prime}=0^{+}$ case, which is accompanied with
$r_{\rm c}=R_{\rm s}+\Delta$. In this case, when the blue circle is moved to
locate its center at $r=R_{\rm s}+\Delta$, Eq. (5) gives $g(R_{\rm
s}+\Delta)>\rho_{0}/4$ as the integration range covers, at least, a quarter
sector of the disc in the dashed circle. Note the coverage is minimized when
$R_{\rm s}^{\prime}=0^{+}$, so that $\rho_{0}/4$ can be a lower bound of
$g(R_{\rm s}+\Delta)$. As the truncated $U_{\rm att}({\bf r},{\bf
r}^{\prime})$ by $r_{\rm c}$ still shows $U_{\rm att}({\bf r},{\bf
r}^{\prime})=U_{\rm att}(|{\bf r}-{\bf r}^{\prime}|)$ and the density in Eqs.
(21) and (22) also does $\rho({\bf r})=\rho(r)$, $|V(0)|<\infty$ follows.
Then, a consequent constant $V({\bf r})$ gives ${\bf v}({\bf r})=-\nabla_{\bf
r}V({\bf r})=0$. Hence, the density in Eq. (22) can be a steady-state
configuration. We call this anomalous static sync (aSS).

The aSS explained so far is conditioned in $r_{\rm c}<2R_{\rm s}$. Meanwhile,
the numerical data manifest there should be a lower bound on $r_{\rm c}$. For
a survey on lower bound, we revisit Fig. 5(d). It states

$R_{\rm s}+R_{\rm s}^{\prime}+\Delta=r_{\rm c}$ (23)

with $R_{\rm s}^{\prime}=R_{\rm s}^{\prime}(r_{\rm c})$ and
$\Delta=\Delta(r_{\rm c})$. As $R_{\rm s}$ is independent of $r_{\rm c}$, the
lower bound of $r_{\rm c}$ is reduced to that of $R_{\rm s}^{\prime}(r_{\rm
c})+\Delta(r_{\rm c})$. Also in the independence, Eq. (23) shows $R_{\rm
s}^{\prime}(r_{\rm c})+\Delta(r_{\rm c})$ is decreasing for $r_{\rm c}$. Since
$R_{\rm s}^{\prime}$ approaches $0$ as $r_{\rm c}$ decreases, the lower bound
$l$ of $r_{\rm c}$ is $l=R_{\rm s}+\Delta(l)$. Thus when $r_{\rm c}=l$, it
reads $\rho(r)=g(r)$ in the absence of the central constant region. For a
simple estimation of $l$, we regard $g(r)$ is linear and
$g(l)\approx\rho_{0}/4$. From normalization $\int d{\bf r}\rho(r)=1$, it
follows $l\approx\sqrt{2}R_{\rm s}$. Thus considering the above-mentioned
upper bound, one finds an interval

$\sqrt{2}R_{\rm s}\lesssim r_{\rm c}<2R_{\rm s}$ (24)

for the formation of aSS, and this is, interestingly, consistent with our
numerical data. In the other side in $r_{\rm c}\lesssim\sqrt{2}R_{\rm s}=1$
(note we set $J=1$ in numerics), we numerically observe the pattern becomes
amorphous. $l\approx\sqrt{2}R_{\rm s}$ is by itself the upper bound of $R_{\rm
s}+\Delta$, the radius of aSS. Depending on initial condition, there can
appear multiple aSS. After transient period, each region separated farther
than $r_{\rm c}$ from the others has its own disc (not shown here).

Our argument so far is also applicable for the change of the asynchronized-
swam disc for $K<K_{\rm c}<0$ while $R_{\rm s}$ is replaced with $R_{\rm a}$.
We call the pattern anomalous static async (aSA). An instance is shown in Fig.
5(c). We remark the async-counterpart of Eq. (24), $\sqrt{2}R_{\rm a}\lesssim
r_{\rm c}<2R_{\rm a}$ is valid for the smaller criterion than $K_{\rm c}$. In
the numerical data obtained for $N=400$, the lower bound shows deviation till
$K\gtrsim-2$ (see the upper-left part of the aSA region in Fig. 8). There is
more deviation for larger $K$ and smaller $r_{\rm c}$ in the corner. This
observation is consistent with the reported stable-region of SA in Ref. [3]
for $r_{\rm c}=\infty$, $K\lesssim K_{\rm c}\approx-1.2J$, in that the
criterion $K_{\rm c}$ is a consequence of linear stability and a smaller
$r_{\rm c}$ means a stronger perturbation. Different from aSS case, aSA does
not become amorphous but a non-stationary pattern when it disappears by the
decrease of $r_{\rm c}$.

###  IV.4 Variations in SPW

The diameter of SPW that appears when $K=0$ for $r_{\rm c}=\infty$ is
$D_{\infty}=2R_{2}$. Thus by Eq. (20), for $r_{\rm c}>2R_{2}$, a formation of
the same SPW(s) is expected, and this is numerically observed. Next, for
$r_{\rm c}<2R_{2}$, the results become different, also as expected. However,
this time, the difference is not that apparent in shape; the inner and outer
radii change a little though the cutoff value is decreased a lot to $r_{\rm
c}\approx 1.6$ from $2.5\approx 2R_{2}$. Figure 6(a) is a pattern at $r_{\rm
c}=1.6$, whose inner and outer radii are not discriminated so clearly from
$R_{1}$ and $R_{2}$, respectively, known for $r_{\rm c}=\infty$ in Eq. (III).

Figure 6: (Color Online) Variations of SPW by finite $r_{\rm c}$ in $N=800$
system with $K=0$ and $J=1$: (a) shows the formed SPW for $r_{\rm c}=1.6$
slightly exceeds the region between the radius $R_{1}$ and $R_{2}$ of Eq.
(III) known for $r_{\rm c}=\infty$. In (b), radial densities for $r_{\rm
c}=2.5$ (red) and $r_{\rm c}=1.6$ (blue) are shown. The red pluses and blue
crosses represent the numerical data, and the bold curves show the overall
profile in each cases.

Instead, we have found that a substantial change in density takes place while
$r_{\rm c}$ decreases a lot. In Fig. 6(b), the red pluses come from the
numerical data obtained at $r_{\rm c}=2.5$ that is slightly smaller than
$2R_{2}=2.53..\,$. Each of them is the count of swarms in a narrow circular
region between radii $r-\delta$ and $r+\delta$, where $r$ is discrete and
small $\delta>0$ is a few times of the resolution of $r$. Their overall
profile is guided by the bold red curve from Eqs. (13), (18), and (III). This
guidance is plausible as $r_{\rm c}=2.5\approx 2R_{2}$. The blue crosses in
Fig. 6(b) come from the numerical data at $r_{\rm c}=1.6$, and the blue bold
curve is their overall shape as a guide to eyes. The substantial change in
density is supposed to compensate for the small change in shape despite the
considerable decrease of $r_{\rm c}$. We also observed the small change of
shape mainly occurs below $r_{\rm c}\approx 1.8$. The inner and outer radii
remain apparently same for $1.8\lesssim r_{\rm c}<2R_{2}$ while the density
change takes place.

A heuristic argument for this observation is following. Let $r_{\rm c}$ be
smaller that $2R_{2}$, and consider a circle of radius $r_{\rm c}$ whose
center is in somewhere in SPW. If the center is near the outer boundary of
SPW, the circle cannot cover the SPW. When this incomplete coverage is
considered in Eq. (5), the reduced integral region would lead to density
decrease at the center. Next, consider the center is near the inner boundary.
In this case, the $r_{\rm c}$-radius circle covers the SPW when $r_{\rm
c}-R_{2}>R_{1}$ is provided. For the density of such center, the integral
region of Eq. (5) is conserved, and there is no explicit reason to decrease
the density. Instead, an opposite reason was already due because the expected
density decrease near the outer region requires mass conservation in the
system. We consider the cross between the density profiles shown in Fig. 6(b)
is the consequence.

If $r_{\rm c}-R_{2}>R_{1}$ does not hold, the integral region of Eq. (5) is
always reduced whatever points of SPW is used for the center of $r_{\rm
c}$-radius circle. This may indicate a termination of the shape-preserving
density change. One then considers an interval

$R_{1}+R_{2}\lesssim r_{\rm c}<2R_{2}\,,$ (25)

where such SPWs of negligible/considerable change in radius/density for
$r_{\rm c}$ are respectively expected. The lower bound suggested this way is
rather comparable with the numerical lower bound $r_{\rm c}\approx 1.8$ in
that $R_{1}+R_{2}=1.74..\approx 1.8$ (see Eq. (III) with $J=1$ for the
$R_{1}+R_{2}$ value). Thus when SPW is specified by the shape only, the
condition for its formation is relaxed from $r_{\rm c}>2R_{2}$ to $r_{\rm
c}\gtrsim R_{1}+R_{2}$.

We remark that some cases do not give SPW when $r_{\rm c}$ begins to be
smaller than $1.8$, depending on initial conditions. No appearance of SPW
below $r_{\rm c}\approx 1.8$ becomes more frequent for the smaller $r_{\rm
c}$. Furthermore, SPW does not appear when $r_{\rm c}\lesssim 1.5$, at least,
numerically. In brief, the numerical data suggests that the basin of attractor
[34] for SPW begins to decrease from $r_{\rm c}\approx 1.8$, and then
disappears below $r_{\rm c}\approx 1.5$.

###  IV.5 Bar-like patterns at $K=0$

It is worthy of noting that, when the pattern’s shape is not circular below
$r_{\rm c}\approx 1.8$, it looks linear. According to the numerical data,
linear shapes begin to appear for $r_{\rm c}\approx 1.8$, and disappear when
$r_{\rm c}\lesssim 1$. Figure 7 shows a few linear bar-like patterns (Bar).

Figure 7: (Color Online) Various bar-like patterns (Bar) for $1.2\lesssim
r_{\rm c}\lesssim 1.8$ in $N=400$ system with $K=0$ and $J=1$: (a)-(d) are
obtained for $r_{\rm c}=1.8,1.6,1.4$, and $1.2$, respectively. The black
rectangle in the pattern of (a) is of length $\pi$ and width $1$.

It is inappropriate we consider to view such patterns as the deformed APW by
finite $r_{\rm c}$. As a first step to the understanding of such bar-like
patterns, we below argue that each of them could be a deformed shape, by its
own finite $r_{\rm c}$, of the bar given as the solution of Eq. (5) in the
$r_{\rm c}\to\infty$ limit. An ironical situation here is no bar-like pattern
has been observed yet in our numerical study neither reported in literature as
far as we know. We thus conjecture that a solution bar exists for Eq. (5), of
no stability, but it becomes stable for finite $1\lesssim r_{\rm c}\lesssim
1.8$ accompanied with proper deformation. If our scenario is valid, the
pattern of Fig. 7(a) is most close to the solution bar among the figures
therein because its $r_{\rm c}=1.8$ is larger than the others’.

We choose $x$-axis as the direction of bar with no loss of generality. In the
numerical data, phase smoothly increases along the $x$-axis. The positive
direction of the $x$-axis is selected in line with the phase increase. We thus
introduce phase field $\theta({\bf r})=\theta(x)$ that is differentiable for
$x$, and then write

$\displaystyle U_{\rm att}({\bf r},{\bf r}^{\prime})$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left((x-x^{\prime})^{2}+(y-y^{\prime})^{2}\right)$
(26)
$\displaystyle~{}~{}~{}~{}\times\left(1+J\cos(\theta(x)-\theta(x^{\prime}))\right).$

in rectangular coordinate. No dependence on $r_{\rm c}$ shows Eq. (26) is that
for the infinite cutoff.

Consider a bar of length $2l$ and of width $2w$, whose center is the origin.
By the symmetry of $\rho(x,y)=\rho(|x|,|y|)$, we choose $\theta(0)=0$ for
simplicity. We then give a relation between $\rho(x,y)$ and $\theta(x)$ as

$\int dy\,\rho({\bf r})=\rho(x)=\frac{1}{2\pi}\frac{d\theta}{dx}$ (27)

in the consideration of i) $\theta_{i}$s are initialized with uniform random
variables in our numerics and ii) any of them does not change in time in Eq.
(2) with $K=0$. By Eq. (27), $\theta(x)$ increases from $\theta(0)=0$ to
$\theta(l)=\pi$, and then remains same thereafter, i.e., $\theta(x)=\pi$ for
$x>l$. By differentiabilty, $d\theta(l)/dx=0$ is given. By the symmetry
mentioned above, it reads that $\theta(-x)=-\theta(x)$.

Plugging Eqs. (26) and (27) into Eq. (5), after a little algebra, one can
obtain

$\rho(x,y)=\frac{1}{\pi}-\frac{J}{2\pi^{2}}\frac{d^{2}}{dx^{2}}\Big{(}c(l)x\sin\theta+s(l)\cos\theta\Big{)}$
(28)

for $c(l)=l+\int_{0}^{l}dx\cos\theta$ and $s(l)=\int_{0}^{l}dxx\sin\theta$.
Obviously, $\rho(x,y)$ of Eq. (28) is considered for $|x|<l$ and $|y|<w$ while
$\rho(x,y)=0$ elsewhere.

One integrates Eq. (28) to know the area $S$ of the solution bar is

$S=4wl=\pi.$ (29)

Interestingly, the pattern in Fig. 7(a) covers a bar of area $\pi$ by length
$\pi$ and width $1$, with a rather tolerable margin (see the guiding black
rectangle). This numerical observation suggests $w\approx 1/2$ and
$l\approx\pi/2$. We also observed the pattern becomes elongated and broadening
as $r_{\rm c}$ decreases. Thus the patterns in Fig. 7(b)-(d) for the smaller
$r_{\rm c}$ also cover the rectangle with the more margins (one can imagine
the rectangle in each ones without difficulty). This is consistent with our
scenario for the birth of bar-like patterns with finite $r_{\rm c}$, proposed
above. We also observe in our numerical data that the bar-like pattern is so
deformed for $r_{\rm c}\lesssim 1$ that it becomes amorphous.

The rectangle we have specified is a consequence of numerically motivated
speculation encapsulated in Eq. (28). This equation is again converted with
Eqs. (27) and (29) to

$\frac{d\theta}{dx}=\frac{\frac{2\pi}{J}\left(x-l\theta/\pi\right)-c(l)\sin\theta}{c(l)x\cos\theta-s(l)\sin\theta}\,.$
(30)

Its solution will lead to the self-consistent equations for $s(l)$ and $c(l)$,
which can fix the value of $l$ and next that of $w$ in use of Eq. (29). This
way, finding the size of the bar is dual to knowing its density profile. The
shooting method looks one of the practical approaches for finding $l$. We
leave related interesting task for future work as remarking $l\approx\pi/2$
and $w\approx 1/2$ is expected with the numerical data up to now.

We add that an 1-dimension result could be relevant to the size of the bar. In
1-dimension, $U_{\rm rep}(r,r^{\prime})=-|r-r^{\prime}|$ giving
$\nabla_{r}^{2}U_{\rm rep}(r,r^{\prime})=2\delta(r-r^{\prime})$ is responsible
for the repulsion part. Thus when $U_{\rm
att}(r,r^{\prime})=(1/2)(r-r^{\prime})^{2}(1+J)$ is considered, uniform
distribution of width $W(J)=2/(1+J)$ follows. We here consider the
swarmalators on $y$-direction segment in the bar are effectively in
1-dimension in that their phases are same so that $U_{\rm att}({\bf r},{\bf
r^{\prime}})=(1/2)(y-y^{\prime})^{2}(1+J)$, the distribution is uniform in the
segment as written in Eq. (28), and its width in numerics with $J=1$ is
$2w\approx 1=W(1)$. We leave expected explicit connection between the
1-dimension result and the bar for future work.

The steady-state condition Eq. (3) is a sufficient condition for Eq. (5). That
is, the solution of Eq. (5) is not necessarily a steady state. We note
$\rho(x)$ of Eq. (28) is not a steady state. For this, one can check the
$y$-component of $v_{\rm att}({\bf r})$ given by $\int d{\bf
r}^{\prime}(-\partial_{y}U_{\rm att}({\bf r},{\bf r}^{\prime}))\rho({\bf
r}^{\prime})=-y$ cannot balance with its repulsion-counterpart for such
$\rho(x)$ that vanishes at $x=l$ (note $\rho(l)\propto d\theta(l)/dx=0$).
However, with $r_{\rm c}=1.8$, we observe a steady-state pattern that is
rather close to the rectangle expected with Eq. (28), as depicted in Fig.
7(a). We interpret that the bar we speculated with $r_{\rm c}=\infty$ is
properly deformed for $1\lesssim r_{\rm c}\lesssim 1.8$, and then it becomes
stable while still keeping linear shape more or less. The patterns in Fig. 7
are the instances.

##  V Summary

In this paper, we considered the population of coupled oscillators that sync
and swarm under finite-ranged interaction. We explored how the interaction
cutoff affects the steady states of the system, with particular attention to
the patterns of the steady states. We found that new patterns are induced by
the finite-ranged interaction. When the positive/negative phase coupling is
present (in the negative case, the coupling below a critical value is
considered), multiple clusters forming uniform sync/async discs are found in
the static sync/async states. Anomalous discs with nonuniform radial density
were also found. In the absence of the phase coupling, APW with abnormal
radial density and bar-like patterns of the phase wave state are discovered.
We analyzed the new patterns theoretically and numerically to understand their
birth and character. The results are summarized in Fig. 8.

Figure 8: (Color Online) Diagram for steady-state patterns for $J=1$ in
$(K,r_{\rm c})$ plane: Boundaries are drawn with solid lines while
approximated ones with grey blocks. Diagram for $K=0$ is depicted in the $K=0$
band. The boundary by $K_{c}$ is the consequence of the stability analysis in
Ref. [3]. The numerical values for the lower bounds of Bar/SPW and that below
$K_{c}$ are specific to $J=1$ only. The data points around the left-upper part
of aSA show the numerical boundary of aSA near there (aSA patterns appear at
blue crosses while not at red pluses for $N=400$). The other data points below
indicate the argued region of aSA is numerically realized for $K\lesssim-2$.

The patterns found for $K\geq 0$ are static in the fact the dynamics becomes a
gradient flow in the long time limit. But, it is not trivial whether all the
patterns for $K<0$ are stationary because the phase dynamics is frustrated for
negative $K$. In this situation, we gave an argument for aSA based on the
stability of SA in $K\lesssim-1.2J$, reported in Ref. [3]. The lower bound of
$r_{\rm c}$ for aSA, assessed by $\sqrt{2}R_{\rm a}$ in our argument, does not
reach $\sqrt{2}R_{\rm a}$ when $K\gtrsim-2$, as shown in Fig. 8. The numerical
data obtained at the red-cross points near there but still in $r_{\rm
c}\gtrsim\sqrt{2}R_{\rm a}$ give completely different patterns that hardly
look like instances of aSA. Those are rather comparable with bar-like ones
(not shown here). This implies that that part does not belong to aSA but to
non-stationary region, on which understanding is out of the scope of this
work. A systematic investigation on the left and upper boundary of aSA could
be a starting point in studying the non-stationary states by finite $r_{\rm
c}$. We believe that the steady-state features reported in this work will
provide their own starting points in extending the understanding to non-
stationary region. Related works will appear elsewhere.

###### Acknowledgements.

We thank to Agata Barcis and Prof. Christian Bettstetter for useful
discussions in the early stage of the study. This study was supported by
NRF-2018R1A2B6001790 and 2021R1A2B5B01001951 (H.H), and 2018R1D1A1B07049254
(H.K.L).

DATA AVAILABILITY

The data that support the findings of this study are available within the
article.

REFERNCES

## References

  * Bouffanais [2016] R. Bouffanais,  _Design and Control of Swarm Dynamics_ (Springer, Singapore, 2016).
  * Fetecau, Huang, and Kolokolnikov [2011] R. C. Fetecau, Y. Huang, and T. Kolokolnikov, “Swarm dynamics and equilibria for a nonlocal aggregation model,” Nonlinearity 24, 2681–2716 (2011).
  * O’Keeffe, Hong, and Strogatz [2017] K. P. O’Keeffe, H. Hong, and S. H. Strogatz, “Oscillators that sync and swarm,” Nature Comm. 8, 1504 (2017).
  * Pikovsky, Rosenblum, and Kurths [2003] A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: A universal Concept in Nonlinear Sciences_ (Cambridge University Press, 2003).
  * Winfree [1967] A. T. Winfree, “Biological rhythms and the behavior of populations of coupled oscillators,” J. Theor. Biol. 16, 15–42 (1967).
  * Kuramoto [1984] Y. Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_ (Springer, Berlin, 1984).
  * Kuramoto and Nishikawa [1987] Y. Kuramoto and I. Nishikawa, “Statistical macrodynamics of large dynamical systems. case of a phase transition in oscillator communities,” J. Stat. Phys. 49, 569–605 (1987).
  * Crawford [1994] J. D. Crawford, “Amplitude expansions for instabilities in populations of globally-coupled oscillators,” J. Stat. Phys. 74, 1047–1084 (1994).
  * Crawford [1995] J. D. Crawford, “Scaling and singularities in the entrainment of globally coupled oscillators,” Phys. Rev. Lett. 74, 4341–4344 (1995).
  * Crawford and Davies [1999] J. D. Crawford and K. T. R. Davies, “Synchronization of globally coupled phase oscillators: singularities and scaling for general couplings,” Physica D 125, 1–46 (1999).
  * Strogatz [2000] S. H. Strogatz, “From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators,” Physica D 143, 1–20 (2000).
  * Acebron  _et al._ [2005] J. A. Acebron, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, “The kuramoto model: A simple paradigm for synchronization phenomena,” Rev. Mod. Phys. 77, 137–185 (2005).
  * Hong, Park, and Choi [2005] H. Hong, H. Park, and M. Y. Choi, “Collective synchronization in spatially extended systems of coupled oscillators with random frequencies,” Phys. Rev. E 72, 036217 (2005).
  * Hong  _et al._ [2007] H. Hong, H. Chate, H. Park, and L.-H. Tang, “Entrainment transition in populations of random frequency oscillators,” Phys. Rev. Lett. 99, 184101 (2007).
  * Hong and Strogatz [2011] H. Hong and S. H. Strogatz, “Kuramoto model of coupled oscillators with positive and negative coupling parameters: An example of conformist and contrarian oscillators,” Phys. Rev. Lett. 106, 054102 (2011).
  * Partridge [1982] B. L. Partridge, “The structure and function of fish schools,” Scientific American 246, 114–123 (1982).
  * Reynolds [1987] C. W. Reynolds, “Flocks, herds, and schools: A distributed behavioral model,” Computer Graphics 21, 25–34 (1987).
  * Vicsek  _et al._ [1995] T. Vicsek, A. Czirok, E. Ben-Jacob, and I. Cohen, “Novel type of phase transition in a system of self-driven particles,” Phys. Rev. Lett. 75, 1226–1229 (1995).
  * O’Loan and Evans [1998] O. J. O’Loan and M. R. Evans, “Alternating steady state in one-dimensional flocking,” J. Phys. A: Math. Gen. 32, L99–L105 (1998).
  * Czirok, Barabasi, and Vicsek [1999] A. Czirok, A.-L. Barabasi, and T. Vicsek, “Collective motion of self-propelled particles: kinetic phase transition in one dimension,” Phys. Rev. Lett. 82, 209–212 (1999).
  * Mogilner and Edelstein-Keshet [1999] A. Mogilner and L. Edelstein-Keshet, “A non-local model for a swarm,” J. Math. Biol. 38, 534–570 (1999).
  * Rauch, Millonas, and Chialvo [1995] E. Rauch, M. Millonas, and D. Chialvo, “Pattern formation and functionality in swarm models,” Phys. Lett. A 207, 185–193 (1995).
  * Hubbard  _et al._ [2004] S. Hubbard, P. Babak, S. T. Sigurdsson, and K. G. Magnússon, “A model of the formation of fish schools and migrations of fish,” Ecol. Model. 174, 359–374 (2004).
  * Lee  _et al._ [2004] H. K. Lee, R. Barlovic, M. Schreckenberg, and D. Kim, “Mechanical restriction versus human overreaction triggering congested traffic states,” Phys. Rev. Lett. 92, 238702 (2004).
  * Topaz, Bertozzi, and Lewis [2006] C. M. Topaz, A. L. Bertozzi, and M. A. Lewis, “A nonlocal continuum model for biological aggregation,” Bull. Math. Bio. 68, 1601–1623 (2006).
  * O’Keeffe, Evers, and Kolokolnikov [2018] K. P. O’Keeffe, J. H. M. Evers, and T. Kolokolnikov, “Ring states in swarmalator systems,” Phys. Rev. E 98, 022203 (2018).
  * O’Keeffe and Bettstetter [2019] K. O’Keeffe and C. Bettstetter, “A review of swarmalators and their potential in bio-inspired computing,” in  _Proc. SPIE 10982, Micro- and Nanotechnology Sensors, Systems, and Applications XI_ (2019) p. 10982.
  * Ha  _et al._ [2019] S.-Y. Ha, J. Jung, J. Kim, J. Park, and X. Zhang, “Emergent behaviors of the swarmalator model for position-phase aggregation,” Mathematical Models and Methods in Applied Sciences 29, 2271–2320 (2019).
  * Barciś, Barciś, and Bettstetter [2019] A. Barciś, M. Barciś, and C. Bettstetter, “Robots that sync and swarm: A proof of concept in ros 2,” (IEEE, 2019).
  * Baris and Bettstetter [2020] A. Baris and C. Bettstetter, “Sandsbots: Robots that sync and swarm,” IEEE Access 8, 218752–218764 (2020).
  * McLennan-Smith, Roberts, and Sidhu [2020] T. A. McLennan-Smith, D. O. Roberts, and H. S. Sidhu, “Oscillatory behavior in a system of swarmalators with a short-range repulsive interaction,” Phys. Rev. E 102, 032607 (2020).
  * Lizarraga and de Aguiar [2020] J. U. F. Lizarraga and M. A. M. de Aguiar, “Synchronization and spatial patterns in forced swarmalators,” Chaos 30, 053112 (2020).
  * Jiménez-Morales [2020] F. Jiménez-Morales, “Oscillatory behavior in a system of swarmalators with a short-range repulsive interaction,” Phys. Rev. E 101, 062202 (2020).
  * Strogatz [1994] S. H. Strogatz,  _Nonlinear Dynamics and Chaos_ (Addison-Wesley, MA, 1994).
  * Arfken and Weber [2005] G. B. Arfken and H. J. Weber,  _Mathematical Methods for Physicists_ , 6th ed. (Elsevier, Amsterdam, 2005).


#  Entopy release in Electroweak Phase Transition in 2HDM

Arnab Chaudhuri
Department of Physics and Astronomy, Novosibirsk State University,
email<EMAIL_ADDRESS>
Maxim Yu. Khlopov
Institute of Physics, Southern Federal University,
Université de Paris, CNRS, Astroparticule et Cosmologie and
National Research Nuclear University "MEPHI"
email<EMAIL_ADDRESS>
Shiladitya Porey
Department of Physics and Astronomy, Novosibirsk State University,
email<EMAIL_ADDRESS>
###### Abstract

Electroweak phase transition in the simplest extension of the standard model
namely two Higgs doublet model and entropy production within this framework is
studied. We have considered several benchmark points which were called using
BSMPT, a C++ package, within the limit of $vev/T_{C}>0.2$ are studied, and
corresponding entropy productions are shown in this paper.

Keywords: Electroweak symmetry breaking, Entropy.

##  1 Introduction

For a successful explanation about the origin of excess baryons over
antibaryons in the universe through electroweak baryogenesis (EWBG), a strong
first-order electroweak phase transition (EWPT) in the early universe is
necessary. Cosmic EWPT happened when the hot universe cooled down enough in
the primeval time so that the potential of the Higgs field got and settled at
a non-zero minimum and in consequence, the symmetry of the theory
$SU(2)_{L}\times U(1)_{Y}$ broke to $U(1)_{\rm em}$. At the time of first-
order EWPT, bubbles of the broken phase originate and baryon-antibaryon
asymmetry generates outside the wall of the bubbles of the broken phase.
However, after the discovery of the standard model (SM) Higgs boson, it is
widely known that EWPT in SM with a single Higgs field is just a smooth cross-
over. Therefore, for a successful EWBG, a theory of EWPT in beyond SM (BSM) is
needed [1].

On the other side, $\sim 26.5\%$ of the total energy density of the universe
is contributed by the dark matter (DM) whose mysterious nature has not been
unveiled till now. Although, primordial black holes and MACHOs which are
considered as one of the viable baryonic DM candidates, it is now clear that
they are unable to contribute completely to the DM energy density of the
universe. There are theories about multicharged extension of the standard
model like dark atoms which can be viable dark matter candidates, [2]. But
there are no experimental evidences as of now.

Not only about the baryogenesis, but there is also no irrefutable theory in SM
about nonbaryonic DM particle which can successfully explain all the
observations.

Similar to the above mentioned facts, there are many limitations of the SM.
Thus scientists are desperately searching for experimental evidence of BSM.
For them the recent result from Fermilab about $g_{\mu}-2$ for muon may be a
ray of hope. $g_{\mu}$ is the gyromagnetic ratio of muon which is defined as
the ratio of magnetic moment to the angular moment of muon and whose value is
$2$ from tree-level calculation. If we define $a_{\mu}=(g_{\mu}-2)/2$, then
higher order loop corrections from SM gives
$a_{\mu}=116,591,810(43)~{}\times~{}10^{-11}$ where the value measured from
Fermilab is $16,592,061(41)~{}\times~{}10^{-11}$ which differs from SM at
$3.3\sigma$ level [3]. This contradiction is actually buttressed the
previously claimed result from the E821 experiment at Brookhaven National Lab
(BNL). There are numerous explanation for this anomalous result including the
existence of BSM.

Among all the BSM theories, the two Higgs doublet model (2HDM) is one of the
most popular theories which not only exhibits strong first-order EWPT for the
proper choice of parameter space but also provides the minimal
phenomenological description of some effects of the supersymmetric model
predicting two Higgs boson doublet. In addition to that, this model can
produce dark matter particles [4] and gives a satisfactory explanation for the
$g_{\mu}$ anomaly of muon [5] if the parameter space is properly chosen.

At or around the epoch of EWPT the energy density of the universe was
dominated by relativistic species with negligible chemical potential. In
addition to that, the universe was almost always in thermal equilibrium except
some special epochs. Thus entropy density per comoving volume of the
relativistic plasma was conserved. However, EWPT is a strongly thermally non-
equilibrium process and thus there is a possibility that entropy might have
been generated during this cosmic process.

In this work, we explored the increase in entropy during the epoch of EWPT in
the real type-I 2HDM framework. We have shown that entropy density per
comoving volume increases if EWPT happens as a first-order phase transition in
2HDM model.

The article is arranged as follows: In the next section the Lagriangian of the
model along with the results are given. A generic conclusion follows and in
the appendix the detailed potential is mentioned.

##  2 Lagrangian density of the model

The Lagrangian density of EWPT theory in real type-I 2HDM is given by

$\displaystyle\mathcal{L}=\mathcal{L}_{\rm
gauge,kin}+\mathcal{L}_{f}+\mathcal{L}_{\rm Yuk}+\mathcal{L}_{\rm
Higgs}-V(\Phi_{1},\Phi_{2},T)$ (1)

where $\mathcal{L}_{\rm gauge,kin}$, $\mathcal{L}_{f}$ and $\mathcal{L}_{\rm
Yuk}$ are the kinetic energy term of gauge bosons ($W_{\alpha}$ and
$B_{\alpha}$ with $\alpha=0,1,2,3$) , kinetic energy of fermions and Yukawa
interaction term of fermions with Higgs bosons. These terms are defined in [6,
7] and also discussed in Appendix Appendix A. Throughout this article, all the
Greek indices used in super or sub-script run from $0$ to $3$ and Latin
indices from $1$ to $3$ if not mentioned otherwise.

$\mathcal{L}_{\rm Higgs}$ incorporates the kinetic term of the Higgs field and
their interaction with the gauge bosons. Thus

$\displaystyle\mathcal{L}_{\rm Higgs}=\left\\{(\partial^{\mu}+i{\cal
W}^{\mu})\Phi_{a}\right\\}^{\dagger}\left\\{(\partial_{\mu}+i{\cal
W}_{\mu})\Phi_{a}\right\\}$ (2)

where $a=1,2$ for two Higgs field, $i=\sqrt{(-1)}$, $i{\cal
W}^{\mu}\equiv+igT^{k}{W^{k}}^{\mu}+ig^{\prime}YB^{\mu}$, and $g$ and
$g^{\prime}$ are coupling constants, $T^{i}$ is the generator of $SU(2)_{L}$
(left-Chiral), which is also a form of Pauli matrices, and $Y$ is the hyper-
charge generator of the $U(1)$.

The total CP-conserving potential for our 2HDM model considered is

$\displaystyle
V(\Phi_{1},\Phi_{2},T)=V_{tree}(\Phi_{1},\Phi_{2})+V_{CW}(\Phi_{1},\Phi_{2})+V_{T}(T){+V_{\rm
daisy}(T)}$ (3)

The tree-level potential can be written as

$\displaystyle V_{\rm tree}(\Phi_{1},\Phi_{2})=$ $\displaystyle
m_{11}^{2}\Phi_{1}^{\dagger}\Phi_{1}+m_{22}^{2}\Phi_{2}^{\dagger}\Phi_{2}-\left[m_{12}^{2}\Phi_{1}^{\dagger}\Phi_{2}+m_{12}^{*}\Phi_{2}^{\dagger}\Phi_{1}\right]+\frac{1}{2}\lambda_{1}\left(\Phi_{1}^{\dagger}\Phi_{1}\right)^{2}$
$\displaystyle+\frac{1}{2}\lambda_{2}\left(\Phi_{2}^{\dagger}\Phi_{2}\right)^{2}+\lambda_{3}\left(\Phi_{1}^{\dagger}\Phi_{1}\right)\left(\Phi_{2}^{\dagger}\Phi_{2}\right)+\lambda_{4}\left(\Phi_{1}^{\dagger}\Phi_{2}\right)\left(\Phi_{2}^{\dagger}\Phi_{1}\right)$
(4)
$\displaystyle+\left[\frac{1}{2}\lambda_{5}\left(\Phi_{1}^{\dagger}\Phi_{2}\right)^{2}+\frac{1}{2}\lambda_{5}^{*}\left(\Phi_{2}^{\dagger}\Phi_{1}\right)^{2}\right].$

$m_{12}^{2}$, $m_{11}^{2}$, and $m_{22}^{2}$ can be estimated from the
following formula

$\displaystyle m_{12}^{2}$ $\displaystyle=$ $\displaystyle 100^{2}\mathrm{\
GeV}^{2}$ (5) $\displaystyle m_{11}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{4v_{1}}\left(-2\lambda_{1}v_{1}^{3}+4m_{12}^{2}v_{2}-2\lambda_{3}v_{1}v_{2}^{2}-2\lambda_{4}v_{1}v_{2}^{2}-\lambda_{5}v_{1}v_{2}^{2}-v_{1}v_{2}^{2}\lambda_{5}\right)$
(6) $\displaystyle m_{22}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{4v_{2}}\left(4m_{12}^{2}v_{1}-2\lambda_{3}v_{1}^{2}v|_{2}-2\lambda_{4}v_{1}^{2}v_{2}-\lambda_{5}v_{1}^{2}v_{2}-2\lambda_{2}v_{2}^{3}-v_{1}^{2}v_{2}\lambda_{5}\right)$
(7)

The value of $m_{12}^{2}$ can alter for different parameter space. These
formulas are valid since $\lambda_{5}$ is real and
$\lambda_{6}=\lambda_{7}=0$. The $\lambda_{1-5}$ can be calculated from the
parameter space as

$\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{m_{H}^{2}\cos{\alpha}^{2}+m_{h}^{2}\sin\alpha^{2}-m_{12}^{2}\tan\beta}{v^{2}\cos\beta^{2}},$
(8) $\displaystyle\lambda_{2}$ $\displaystyle=$
$\displaystyle\frac{m_{H}^{2}\sin\alpha^{2}+m_{h}^{2}\cos{\alpha}^{2}-m_{12}^{2}\tan\beta^{-1}}{v^{2}\sin\beta^{2}},$
(9) $\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle\frac{(m_{H}^{2}-m_{h}^{2})\sin\alpha\cos{\alpha}+2m_{H^{\pm}}^{2}\cos\beta\sin\beta-
m_{12}^{2}}{v^{2}\sin\beta\cos\beta},$ (10) $\displaystyle\lambda_{4}$
$\displaystyle=$
$\displaystyle\frac{(m_{A}^{2}-2m_{H^{\pm}}^{2})\sin\beta\cos\beta+m_{12}^{2}}{v^{2}\sin\beta\cos\beta},$
(11) $\displaystyle\lambda_{5}$ $\displaystyle=$
$\displaystyle\frac{m_{12}^{2}-m_{A}^{2}\sin\beta\cos\beta}{v^{2}\sin\beta\cos\beta}.$
(12)

where $v$ is the standard model expectation value,
$v^{2}=v_{1}^{2}+v_{2}^{2}$, $\tan\beta=v_{2}/v_{1}$ and
$\cos(\beta-\alpha)\to 0$ leads to SM result. The details about the parameter
space of $m_{H}$, $m_{h}$, $m_{H\pm}$ can be found in the recent works [8, 9,
10].

The Coleman–Weinberg correction to the potential -

$\displaystyle V_{\rm
CW}\left(v_{1}+v_{2}\right)=\sum_{j}\frac{n_{j}}{64\pi^{2}}(-1)^{2s_{j}}m_{j}^{4}\left(v_{1},v_{2}\right)\left[\log\left(\frac{m_{j}^{2}\left(v_{1},v_{2}\right)}{\mu^{2}}\right)-c_{j}\right]$
(13)

The values of $n_{j}$, $s_{j}$, $c_{j}$ and different mass-values
$m_{j}^{2}\left(v_{1},v_{2}\right)$ are mention in Appendix 4 and
$\mu=246\mathrm{\ GeV}$.

Temperature correction of potential and its series expansion in Landau gauge
are

$\displaystyle V_{T}=\frac{T^{4}}{2\pi^{2}}\left(\sum_{j={\rm
bosons}}n_{j}J_{B}\left[\frac{m_{j}^{2}(v_{1},v_{2})}{T^{2}}\right]+\sum_{j={\rm
fermions}}n_{j}J_{F}\left[\frac{m_{j}^{2}(v_{1},v_{2})}{T^{2}}\right]\right)$
(14) $\displaystyle
T^{4}J_{B}\left[\frac{m^{2}}{T}\right]=-\frac{\pi^{4}T^{4}}{45}+\frac{\pi^{2}}{12}T^{2}m^{2}-\frac{\pi}{6}T(m^{2})^{3/2}-\frac{1}{32}m^{4}\ln\frac{m^{2}}{a_{b}T^{2}}+\cdots,$
(15) $\displaystyle
T^{4}J_{F}\left[\frac{m^{2}}{T}\right]=\frac{7\pi^{4}T^{4}}{360}-\frac{\pi^{2}}{24}T^{2}m^{2}-\frac{1}{32}m^{4}\ln\frac{m^{2}}{a_{f}T^{2}}+\cdots,$
(16)

where $a_{b}=16a_{f}=16\pi^{2}\exp(3/2-2\gamma_{E})$ with $\gamma_{E}$ being
the Euler-Mascheroni constant.

The daisy term is defined as

$\displaystyle V_{\text{daisy}}(T)$
$\displaystyle=-\frac{T}{12\pi}\sum\limits_{i=1}\left[\left(M_{i}^{2}\left(v_{1},v_{2},T\right)\right)^{3/2}-\left(m_{i}^{2}\left(v_{1},v_{2}\right)\right)^{3/2}\right]$
(17)

Details about the $M_{i}^{2}\left(v_{1},v_{2},T\right)$ term can be found in
[12, 13]. Actually, we will see later that all these terms will be taken care
of the software package we have used for this work.

At sufficiently high temperature, the total potential of eq.(3) has only one
minimum at $\left<\Phi_{1}\right>=\left<\Phi_{2}\right>=0$ and there is no
symmetry breaking. The critical temperature ($T_{c}$) is defined as the
temperature at which if the temperature drops down, the total potential gets a
second minimum at $(\Phi_{a,{\rm
min}})\equiv\left\\{\left<\Phi_{1}\right>=v_{1},\left<\Phi_{2}\right>=v_{2}\right\\}$.
For simplicity, we are assuming in this work that both of the Higgs field
$\Phi_{1}$ and $\Phi_{2}$ get the second minimumat the same temperature
$T_{c}$ at the same time. Thus

$\displaystyle
V\left(\Phi_{1}=0,\Phi_{2}=0,T_{c}\right)=V\left(\Phi_{1}=v_{1},\Phi_{2}=v_{2},T_{c}\right).$
(18)

As soon as the Higgs potential gets a non-zero minimum, the other relativistic
particles starts to gain mass and becomes non-relativistic. The reaction rate
among them and also with photon becomes comparable with the Hubble parameter
and thus decouples from relativistic plasma. The mass of the particle and
coupling constant determine the decoupling temperature. For instance, top
quark decouples earlier than electron or other quarks.

Now, at the time of EWPT the universe can be assumed as perfectly homogeneous
and isotropic and thus we can neglect the spatial partial derivatives of the
Higgs fields. Therefore, when the Higgs fields start to oscillate around their
minima $(\Phi_{a,{\rm min}})$ then energy density $\rho$ and pressure $P$ are

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle\dot{\Phi}_{a,{\rm
min}}^{2}+V_{\rm tot}(\Phi_{1},\Phi_{2},T)+\frac{g_{*}\pi^{2}}{30}T^{4}.$ (19)
$\displaystyle P$ $\displaystyle=$ $\displaystyle\dot{\Phi}_{a,{\rm
min}}^{2}-V_{\rm
tot}(\Phi_{1},\Phi_{2},T)+\frac{1}{3}\frac{g_{*}\pi^{2}}{30}T^{4}$ (20)

The last terms in eq.(19) and eq.(20) arise from the Yukawa interaction
between fermions and Higgs bosons and from the energy density of the fermions,
the gauge bosons, and the interaction between the Higgs and gauge bosons.
$g_{*}$ depends on the effective number of particles present in the
relativistic soup at or near the EWPT. It’s value in our model is greater than
the value in SM.

Since the oscillation of $\Phi_{a}$ around $\Phi_{a,{\rm min}}$ is small
compared to Hubble expansion, we can neglect the time derivative of
$\dot{\Phi}_{a,{\rm min}}$ [11] for simplicity in this work.

Again, entropy density per comoving volume is defined as

$\displaystyle s=\frac{\rho+P}{T}a^{3}$ (21)

which is conserved for relativistic species with negligible chemical
potential. From eq.(19) and eq.(20) we get

$\displaystyle\rho+P=2\dot{\Phi}_{a,{\rm
min}}^{2}+\frac{4}{3}\frac{g_{*}\pi^{2}}{30}T^{4}$ (22)

As discussed earlier, $g_{*}$ will change with the decoupling process and thus
$s$ for relativistic plasma will increase for our considered scenario. Then
the increase in entropy can be calculated using conservation of energy
momentum tensor

$\displaystyle\dot{\rho}=-3H\left(\rho+P\right)$ (23)

To solve eq.(23), we have used BSMPT [13, 14], a C++ package to calculate the
vacuum expectation value (VEV) of the total potential, value of the total
potential at VEV for different temperatures including $T_{c}$. We have chosen
the parameter in such a way so that $VEV/T_{c}>0.02$. We have considered five
sets of benchmark values and the corresponding figures are shown in Fig.1

Table 1: 2HDM Benchmark points for entropy production

|  $m_{h}$ [GeV] |  $m_{H}$ [GeV] |  $m_{H^{\pm}}$ [GeV] |  $m_{A}$ [GeV] | $\tan\beta$ | $\cos\left(\beta-\alpha\right)$ | $m_{12}^{2}~{}\text{GeV}^{2}$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$ | $\lambda_{5}$ | $T_{c}$ | $vev/T_{c}$ | $\delta s/s[\%]$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
BM 1 | $125$ | $500$ | $500$ | $500$ | $2$ | $0$ | $10^{5}$ | $0.258$ | $0.258$ | $0.258$ | $0$ | $0$ | $161.36$ | 1.4 | $57$
BM 2 | " | " | $485$ | $500$ | $2$ | $0.00$ | $10^{5}$ | $0.258$ | $.258$ | $-0.23$ | $0.49$ | $0$ | $153.27$ | 1.25 | $53$
BM 3 | " | " | $485$ | $485$ | $2$ | $0.07$ | $10^{5}$ | $1.28$ | $0.002$ | $0.21$ | $0.244$ | $0.244$ | ${168.61}$ | 1.7 | $59$
BM 4 | " | $485$ | $485$ | $485$ | $10$ | $0.1$ | 23,289.6 | $3.9$ | $0.22$ | $3.9$ | $0$ | $0$ | ${230.18}$ | 1.86 | $70$
BM 5 | " | $90$ | $200$ | $300$ | $10$ | $0$ | $801.98$ | $0.258$ | $0.258$ | $1.31$ | $0.3$ | $-1.35$ | ${135.38}$ | 1.06 | $37$

Figure 1: The figures show the entropy production for five different benchmark
(BM) points: Pink line (BM 5, $T_{c}=135.38~{}\text{GeV}$ and $\delta
s/s=37\%$), Red line (BM 2, $T_{c}=153.27~{}\text{GeV}$ and $\delta
s/s=53\%$), Green line (BM 1, $T_{c}=161.36~{}\text{GeV}$ and $\delta
s/s=30\%$), Black line (BM 3, $T_{c}=168.61~{}\text{GeV}$ and $\delta
s/s=59\%$), and Blue line (BM 4 $T_{c}=230.18~{}\text{GeV}$ and $\delta
s/s=70\%$).

##  3 Conclusion

As seen from Fig.1, the entropy productions for some benchmark points are
shown here. A proper difference can be noticed from the standard model
scenario. As seen in [11], the entropy released is around $~{}13\%$ and in the
present scenario, we see that the production is considerably higher. This is
because first-order phase transition as seen in 2HDM can release more entropy
compared to smooth crossover or second-order in the case of the standard
model. The massive scalar particles in 2HDM contribute considerably to this
production as well.

## Acknowledgements

The work of S.P. and A.C. is funded by RSF Grant 19-42-02004. The research by
M.K. was supported by the Ministry of Science and Higher Education of the
Russian Federation under Project "Fundamental problems of cosmic rays and dark
matter", No. 0723-2020-0040.

## References

  * [1] D. E. Morrissey and M. J. Ramsey-Musolf, New J. Phys. 14 (2012), 125003 doi:10.1088/1367-2630/14/12/125003 [arXiv:1206.2942 [hep-ph]].
  * [2] A. Chaudhuri and M. Yu. Khlopov, Universe 7, (2021) 8, 275 doi:10.3390/universe7080275 [arXiv: 2106.11646 [hep-ph]].
  * [3] B. Abi et al. [Muon g-2], Phys. Rev. Lett. 126 (2021) no.14, 141801 doi:10.1103/PhysRevLett.126.141801 [arXiv:2104.03281 [hep-ex]].
  * [4] L. Wang and X. F. Han, Phys. Lett. B 739 (2014), 416-420 doi:10.1016/j.physletb.2014.11.016 [arXiv:1406.3598 [hep-ph]].
  * [5] L. Wang and X. F. Han, JHEP 05 (2015), 039 doi:10.1007/JHEP05(2015)039 [arXiv:1412.4874 [hep-ph]].
  * [6] A. Chaudhuri and M. Y. Khlopov, MDPI Physics 3 (2021) no.2, 275-289 doi:10.3390/physics3020020 [arXiv:2103.03477 [hep-ph]].
  * [7] A. Chaudhuri, M. Y. Khlopov and S. Porey, Galaxies 9 (2021) no.2, 45 doi:10.3390/galaxies9020045 [arXiv:2105.10728 [hep-ph]].
  * [8] O. Eberhardt, A. P. Martínez and A. Pich, JHEP 05 (2021), 005 doi:10.1007/JHEP05(2021)005 [arXiv:2012.09200 [hep-ph]].
  * [9] S. Karmakar, Springer Proc. Phys. 248 (2020), 193-198 doi:10.1007/978-981-15-6292-1_23
  * [10] G. C. Dorsch, S. J. Huber, K. Mimasu and J. M. No, Phys. Rev. D 93 (2016) no.11, 115033 doi:10.1103/PhysRevD.93.115033 [arXiv:1601.04545 [hep-ph]].
  * [11] A. Chaudhuri and A. Dolgov, JCAP 01 (2018), 032 doi:10.1088/1475-7516/2018/01/032 [arXiv:1711.01801 [hep-ph]].
  * [12] J. Bernon, L. Bian and Y. Jiang, JHEP 05 (2018), 151 doi:10.1007/JHEP05(2018)151 [arXiv:1712.08430 [hep-ph]].
  * [13] P. Basler and M. Mühlleitner, Comput. Phys. Commun. 237 (2019), 62-85 doi:10.1016/j.cpc.2018.11.006
  * [14] P. Basler, M. Mühlleitner and J. Müller, Comput. Phys. Commun. 269 (2021), 108124 doi:10.1016/j.cpc.2021.108124 [arXiv:2007.01725 [hep-ph]].
  * [15] P. Basler, M. Krause, M. Muhlleitner, J. Wittbrodt and A. Wlotzka, JHEP 02 (2017), 121 doi:10.1007/JHEP02(2017)121 [arXiv:1612.04086 [hep-ph]].

## Appendix A

The kinetic energy term of gauge bosons, kinetic energy of fermions and Yukawa
interaction term of fermions with Higgs bosons are

$\displaystyle\mathcal{L}_{f}=\sum_{\Psi=Q,L,u,d,l}i\left(\bar{\Psi}_{L}\not{D}\Psi_{L}+\bar{\Psi}_{R}\not{D}\Psi_{R}\right)$
(24) $\displaystyle\mathcal{L}_{\rm
Yuk}=-\left[y_{e}\bar{e_{R}}\Phi_{a}^{\dagger}L_{L}+y_{e}^{*}\bar{L_{L}}\Phi_{a}^{\dagger}e_{R}+\cdots\right]$
(25) $\displaystyle\mathcal{L}_{\rm
gauge,kin}=-\frac{1}{4}G^{j}_{\mu\nu}{G^{j}}^{\mu\nu}-\frac{1}{4}F^{B}_{\mu\nu}{F^{B}}^{\mu\nu}$
(26)

where $\Psi$ is the fermionic field, subscript $L$ ($R$) is for the left
(right) chiral field. The sum in eq.(25) is also over quarks. $y_{e}$ is the
complex constant and

$\displaystyle
G^{j}_{\mu\nu}=\partial_{\mu}W^{j}_{\nu}-\partial_{\nu}W^{j}_{\mu}-g\epsilon^{jkl}W_{\mu}^{k}W_{\nu}^{l}$
(27) $\displaystyle
F^{B}_{\mu\nu}=\partial_{\mu}B_{\nu}-\partial_{\nu}B_{\mu}$ (28)
$\displaystyle\not{D}\Psi^{(j)}_{L,R}\equiv\gamma^{\mu}(\partial_{\mu}+igW_{\mu}+ig^{\prime}Y_{L,R}B_{\mu})\Psi^{(j)}_{L,R}$
(29)

##  4 Appendix B: Masses of new Scalars

$\displaystyle
c_{i}=\begin{cases}\frac{5}{6},&\left(i=W^{\pm},Z,\gamma\right)\\\
\frac{3}{2},&\text{otherwise}\end{cases}$ (30)

Bosons | $n_{i}$ | $s_{i}$ | $m(v)^{2}$ |
---|---|---|---|---
$h$ | $1$ | $1$ | eigenvalues of 42 | Higgs
$H$ | $1$ | $1$ | eigenvalues of 42 | Higgs
$A$ | $1$ | $1$ | eigenvalues of 42 | Higgs
$G^{0}$ | $1$ | $1$ | eigenvalues of 42 | Goldstone
$H^{\pm}$ | $2$ | $1$ | Eq.34 | Charged Higgs
$G^{\pm}$ | $2$ | $1$ | Eq.35 | Charged Goldstone
$Z_{L}$ | $1$ | $1$ | Eq.32 | Higgs
$Z_{T}$ | $2$ | $2$ | Eq.32 | Higgs
$W_{L}$ | $2$ | $1$ | Eq.31 | Higgs
$W_{T}$ | $4$ | $2$ | Eq.31 | Higgs
$\gamma_{L}$ | $1$ | $2$ | Eq.33 |
$\gamma_{T}$ | $2$ | $2$ | Eq.33 |

$\displaystyle m_{W}^{2}=\frac{g^{2}}{4}v^{2}.$ (31) $\displaystyle
m_{Z}^{2}=\frac{g^{2}+g^{\prime 2}}{4}v^{2}.$ (32) $\displaystyle
m_{\gamma}^{2}=0.$ (33)

$\displaystyle\bar{m}_{H^{\pm}}^{2}$
$\displaystyle=\frac{1}{2}\left(\mathcal{M}_{11}^{C}+\mathcal{M}_{22}^{C}\right)+\frac{1}{2}\sqrt{4\left(\left(\mathcal{M}_{12}^{C}\right)^{2}+\left(\mathcal{M}_{13}^{C}\right)^{2}\right)+\left(\mathcal{M}_{11}^{C}-\mathcal{M}_{22}^{C}\right)^{2}}.$
(34) $\displaystyle\bar{m}_{G^{\pm}}^{2}$
$\displaystyle=\frac{1}{2}\left(\mathcal{M}_{11}^{C}+\mathcal{M}_{22}^{C}+\right)-\frac{1}{2}\sqrt{4\left(\left(\mathcal{M}_{12}^{C}\right)^{2}+\left(\mathcal{M}_{13}^{C}\right)^{2}\right)+\left(\mathcal{M}_{11}^{C}-\mathcal{M}_{22}^{C}\right)^{2}}.$
(35)

where

$\displaystyle c_{1}$
$\displaystyle=\frac{1}{48}\left(12\lambda_{1}+8\lambda_{3}+4\lambda_{4}+3\left(3g^{2}+g^{\prime
2}\right)\right)$ (36) $\displaystyle c_{2}$
$\displaystyle=\frac{1}{48}\left(12\lambda_{2}+8\lambda_{3}+4\lambda_{4}+3\left(3g^{2}+g^{\prime
2}\right)+\frac{24}{v_{2}^{2}}m_{t}^{2}(T=0)\right)$ (37)
$\displaystyle+\frac{1}{2v_{2}^{2}}m_{b}^{2}(T=0)$

where $m_{t}(T=0)=172.5{\rm Gev}$ and $m_{b}(T=0)=4.92{\rm GeV}$. For our case
$(v_{3}=0)$,

$\displaystyle\mathcal{M}_{11}^{C}$
$\displaystyle=m_{11}^{2}+\lambda_{1}\frac{v_{1}^{2}}{2}+\lambda_{3}\frac{v_{2}^{2}}{2}$
(38) $\displaystyle\mathcal{M}_{22}^{C}$
$\displaystyle=m_{22}^{2}+\lambda_{2}\frac{v_{2}^{2}}{2}+\lambda_{3}\frac{v_{1}^{2}}{2}$
(39) $\displaystyle\mathcal{M}_{12}^{C}$
$\displaystyle=\frac{v_{1}v_{2}}{2}\left(\lambda_{4}+\lambda_{5}\right)-m_{12}^{2}$
(40) $\displaystyle\mathcal{M}_{13}^{C}$ $\displaystyle=0$ (41)

Masses of $h$, $H$ and $A$ are the eigen values of the matrix

$\displaystyle\bar{\mathcal{M}}^{N}=\left(\mathcal{M}^{N}\right)$ (42)

For our case $(v_{3}=0)$,

$\displaystyle\mathcal{M}_{11}^{N}$
$\displaystyle=m_{11}^{2}+\frac{3\lambda_{1}}{2}v_{1}^{2}+\frac{\lambda_{3}+\lambda_{4}}{2}v_{2}^{2}+\frac{1}{2}\lambda_{5}v_{2}^{2}$
(43) $\displaystyle\mathcal{M}_{22}^{N}$
$\displaystyle=m_{11}^{2}+\frac{\lambda_{1}}{2}v_{1}^{2}+\frac{\lambda_{3}+\lambda_{4}}{2}v_{2}^{2}-\frac{1}{2}\lambda_{5}v_{2}^{2}$
(44) $\displaystyle\mathcal{M}_{33}^{N}$
$\displaystyle=m_{22}^{2}+\frac{3\lambda_{2}}{2}v_{2}^{2}+\frac{1}{2}\left(\lambda_{3}+\lambda_{4}+\lambda_{5}\right)v_{1}^{2}$
(45) $\displaystyle\mathcal{M}_{44}^{N}$
$\displaystyle=m_{22}^{2}+\frac{\lambda_{2}}{2}v_{2}^{2}+\frac{1}{2}\left(\lambda_{3}+\lambda_{4}-\lambda_{5}\right)v_{1}^{2}$
(46) $\displaystyle\mathcal{M}_{12}^{N}$ $\displaystyle=0$ (47)
$\displaystyle\mathcal{M}_{13}^{N}$
$\displaystyle=-m_{12}^{2}+\left(\lambda_{3}+\lambda_{4}+\lambda_{5}\right)v_{1}v_{2}$
(48) $\displaystyle\mathcal{M}_{14}^{N}$ $\displaystyle=0$ (49)
$\displaystyle\mathcal{M}_{23}^{N}$ $\displaystyle=0$ (50)
$\displaystyle\mathcal{M}_{24}^{N}$
$\displaystyle=-m_{12}^{2}+\lambda_{5}v_{1}v_{2}$ (51)
$\displaystyle\mathcal{M}_{34}^{N}$ $\displaystyle=0$ (52)

Table 2: Field dependent mass of all fermions Fermions | $n_{i}$ | $s_{i}$ | $m_{f}(T=0)$ |
---|---|---|---|---
$e$ | $4$ | $\frac{1}{2}$ | $\frac{y_{e}}{\sqrt{2}}v_{k}$ | lepton
$\mu$ | $4$ | $\frac{1}{2}$ | $\frac{y_{\mu}}{\sqrt{2}}v_{k}$ | lepton
$\tau$ | $4$ | $\frac{1}{2}$ | $\frac{y_{\tau}}{\sqrt{2}}v_{k}$ | lepton
$u$ | $12$ | $\frac{1}{2}$ | $\frac{y_{u}}{\sqrt{2}}v_{k}$ | quark
$c$ | $12$ | $\frac{1}{2}$ | $\frac{y_{c}}{\sqrt{2}}v_{k}$ | quark
$t$ | $12$ | $\frac{1}{2}$ | $\frac{y_{t}}{\sqrt{2}}v_{k}$ | quark
$d$ | $12$ | $\frac{1}{2}$ | $\frac{y_{d}}{\sqrt{2}}v_{k}$ | quark
$s$ | $12$ | $\frac{1}{2}$ | $\frac{y_{s}}{\sqrt{2}}v_{k}$ | quark
$b$ | $12$ | $\frac{1}{2}$ | $\frac{y_{b}}{\sqrt{2}}v_{k}$ | quark