# Prolonged hysteresis in the Kuramoto model with inertia and higher-order
interactions

Narayan G Sabhahit1<EMAIL_ADDRESS>Akanksha S. Khurd2 Sarika
Jalan3<EMAIL_ADDRESS>1\. Department of Physical Sciences, Indian Institute
of Science Education and Research, Kolkata  2\. Department of Physics, Indian
Institute of Science Education and Research, Tirupati  3.Complex Systems Lab,
Department of Physics, Indian Institute of Technology Indore, Khandwa Road,
Simrol, Indore-453552, India

###### Abstract

The inclusion of inertia in the Kuramoto model has been long reported to
change the nature of phase transition, providing a fertile ground to model the
dynamical behaviors of interacting units. More recently, higher-order
interactions have been realized as essential for the functioning of real-world
complex systems ranging from the brain to disease spreading. Yet, analytical
insights to decipher the role of inertia with higher-order interactions remain
challenging. Here, we study Kuramoto model with inertia on simplicial
complexes, merging two research domains. We develop an analytical framework in
a mean-field setting using self-consistent equations to describe the steady-
state behavior, which reveals a prolonged hysteresis in the synchronization
profile. Inertia and triadic interaction strength exhibit isolated influence
on system dynamics by predominantly governing, respectively, the forward and
backward transition points. This work sets a paradigm to deepen our
understanding of real-world complex systems such as power grids modeled as the
Kuramoto model with inertia.

##  I Introduction

The emergence of collective behavior in complex real-world systems has been a
long-standing research interest [1]. It was initially in the landmark paper
[2] that Kuramoto modeled the phenomenon of synchronization using a system of
network-coupled oscillators in an analytically tractable setting, illustrating
that the system underwent a second-order phase transition from incoherent to a
coherent state. Since then, numerous works on various extensions of the
Kuramoto Model have been done, revealing several phenomena [4, 3, 5, 6, 7]. Of
particular interest to us is the Kuramoto Model with inertia (also known as
the second-order Kuramoto model). Inspired by the modeling of synchronized
flashing in Pteroptix malaccae by Ermentrout [8], a second-order extension of
the Kuramoto model was first proposed by Tanaka et al. [9, 10]. They showed
that the system experienced a first-order phase transition upon introducing
inertia rather than the smooth second-order phase transition observed in the
Kuramoto model. They put forth a self-consistent method akin to the one
proposed by Kuramoto to study the steady-state behavior of the coupled
oscillator system. Since then, the second-order Kuramoto model has been
extensively explored in diluted networks [11] and various real-world complex
systems like Josephson junctions [12] and power grids [13, 14, 15, 16, 17]. In
[18], Filatrella et al. explained how the second-order Kuramoto model
originates in power grids by simply accounting for power conservation at each
node of the grid, rendering it more than just a mathematical convenience.

However, all these results were obtained by focusing on the interactions to be
purely dyadic in nature. Recent research highlights that such a reductionist
view might not reveal the complete picture of the underlying mechanism of
exotic phenomena observed in some real-world complex systems where the
interactions between agents are inherently higher-order in nature [19, 20, 21,
22]. Using the Ott-Antonsen (OA) dimensionality reduction method [23] Skardal
& Arenas [24] showed that incorporating higher-order interactions into the
Kuramoto model resulted in abrupt (de)synchronization transitions. It was
remarkable to observe that adding an inertia term to the Kuramoto model or
incorporating higher-order interactions independently gave rise to a first-
order phase transition in the system. Here, we are interested in understanding
how the interplay of inertia and higher-order interactions manifests itself in
the system and affects the synchronization profile, which has not been
explored before.

This article unifies these two disparate fields by providing a generalized
analytical framework to predict the steady-state behavior of coupled
oscillator systems with inertia interacting via higher-order interactions. The
first challenge lies in the fact that using dimensionality reduction methods
like the OA does not carry over easily towards analyzing second-order Kuramoto
models [26] as the density function, which is usually Fourier expanded in
terms of the phase now also depends on the velocity of the oscillators. Hence,
in this study, we develop novel analytics using the self-consistency
methodology originally proposed by Kuramoto and later extended to the second-
order Kuramoto model by Tanaka et al. [9]. We show that the effect of inertia
and higher-order interactions manifest independently, resulting in a prolonged
hysteresis-driven synchronization transition within the system. The forward
and backward transition points are found to be predominantly dependent on
inertia and higher-order interaction, respectively Fig. 1.

##  II Model

Figure 1: (Color online) Prolonged Hysteresis. a) Schematic depiction of
emerging collective behavior in the Kuramoto Model (KM). ($a^{\prime}$),
($c^{\prime}$), and ($d^{\prime}$) plot the usual behavior of KM [2] in the
sole impression of higher-order [24] or inertia [9], whereas ($b^{\prime}$)
illustrates a simultaneous forward and backward shift in the transition points
upon introduction of $m$ and $K_{2}$ in KM (Eq. 1), revealing a prolonged
hysteresis. The green arrow indicates the direction of the shift in the
transition points. b) $r_{1}$ versus $K_{1}$ plot for $K_{2}=1$ and $m=1$
(blue-circles) and $K_{2}=7$ and $m=3$ (red-squares). Filled circles and
squares represent the simulation results for the forward, and hollow circles
and squares represent the backward processes. The dashed and continuous curves
represent the forward and backward analytical predictions, respectively.

We study the inertial effects in a globally coupled oscillator network
considering the simultaneous presence of dyadic and triadic interactions.
Phases of $N$-coupled oscillators, each with mass $m$, evolve based on the
following coupled nonlinear equations,

$\begin{split}m\ddot{\theta_{i}}=-\dot{\theta_{i}}+\omega_{i}+\frac{K_{1}}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i})\\\
+\frac{K_{2}}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}\sin(2\theta_{j}-\theta_{k}-\theta_{i}).\end{split}$
(1)

In Eq. 1, $\theta_{i}$ and $\dot{\theta}_{i}$ refer to the instantaneous phase
and angular velocity of the $i^{th}$ oscillator, respectively. $\omega_{i}$ is
the intrinsic frequency of the $i^{th}$ oscillator derived from a unimodal
symmetric probability distribution $g(\omega)$ with mean $\Omega$. The
coupling constants $K_{1}\geq 0$ and $K_{2}$ are the dyadic and triadic
coupling strengths, respectively. The form of the higher order coupling term
directly falls from the phase reduction of mean-field complex Ginzburg-Landau
equation [25].

We decouple the differential equations in Eq. 1 and write them in terms of
mean-field quantities by introducing the following general order parameter for
$p\in\\{1,2\\}$,

$r_{p}e^{i\psi_{p}}=\frac{1}{N}\sum_{j=1}^{N}e^{ip\theta_{j}}.$ (2)

From the above definition, $r_{1}$ measures the global phase coherence and can
be interpreted as the centroid of phases of oscillators on a unit circle in
the complex plane, and $\psi_{1}$ measures the average phase of the
oscillators. $r_{2}$, referred to as the Daido order parameter [27] captures
cluster synchronization. As we are interested in the steady state behavior of
the system, we omit the time dependence in the definition of the general order
parameter. In the incoherent state, the phases of the oscillators are
scattered uniformly on the unit circle and hence $r_{1}\approx r_{2}\approx
0$. Meanwhile, in the coherent state, a single group of oscillators is formed
locked to the mean phase $\psi_{1}$ rotating uniformly at angular velocity
$\Omega$, hence $r_{1}\approx r_{2}\approx 1$. Using Eq. 2, Eq. 1 can be
written in terms of mean-field quantities as,

$\begin{split}m\ddot{\theta}_{i}=-\dot{\theta}_{i}+\omega_{i}+K_{1}r_{1}\sin(\psi_{1}-\theta_{i})\\\
+K_{2}r_{1}r_{2}\sin(\psi_{2}-\psi_{1}-\theta_{i}).\end{split}$ (3)

Because of the rotational symmetry in the model, the mean of the $g(\omega)$
distribution can be set to zero by moving into the rotating frame at the
frequency $\Omega$. This can be facilitated by making the transformation
$\theta_{i}\rightarrow\theta_{i}+\Omega t$ in Eq. 1. Once in the rotating
frame, $\psi_{1}$ and $\psi_{2}$ can be set to zero by appropriately shifting
the origin, i.e., $\theta_{i}(0)\rightarrow\theta_{i}(0)+\psi_{1}(0)$. Eq. 3,
now takes the following form,

$m\ddot{\theta}_{i}=-\dot{\theta}_{i}+\omega_{i}-q\sin(\theta_{i}),$ (4)

where, for the ease of notation, $q=r_{1}(K_{1}+K_{2}r_{2})$.

##  III Analytical Results

Note that for fixed parameter values $K_{1}$ and $K_{2}$, Eq. 4 has two
variables $r_{1}$ and $r_{2}$. Hence, to chalk out the steady state behavior
of Eq. 4, we develop a system of self-consistent equations and seek the values
of $(K_{1},r_{1},r_{2})$ which simultaneously satisfy them. We start by taking
the thermodynamic limit ($N\rightarrow\infty$); the coupled oscillator system
in the steady state is then described by a probability density
$\rho(\theta,\omega)$ where for a given intrinsic frequency $\omega$,
$\rho(\theta,\omega)d\theta$ represents the fraction of oscillators with their
phase between $\theta$ and $\theta+d\theta$. The general order parameter in Eq
2 takes the following form in the continuum limit,

$r_{p}e^{i\psi_{p}}=\int_{-\infty}^{\infty}\int_{-\pi}^{\pi}e^{ip\theta}\rho(\theta,\omega)g(\omega)d\omega
d\theta.$ (5)

In the steady state, the oscillator population splits up into two groups
depending on their intrinsic frequency. One group of oscillators is locked to
the mean phase; meanwhile, the other oscillators drift over the locked
oscillators. Hence the overall phase coherence ($r_{p}$) can be split into
contributions from the locked ($r_{p}^{l}$) and drifting ($r_{p}^{d}$)
oscillators, i.e, $r_{p}=r_{p}^{l}+r_{p}^{d}$. Before calculating $r_{p}^{l}$
and $r_{p}^{d}$, we point out that systems whose motion is governed by Eq. 4
are known to depict hysteresis and have been well studied in [9, 10, 28]. For
completeness, we briefly summarise the reason for the hysteretic behavior
here. Dropping the subscript $i$ and by introducing a new timescale
$\tau=\sqrt{\frac{q}{m}}t$, Eq. 4 is transformed to a second order
differential equation with just two parameters as,

$\ddot{\theta}=-\alpha\dot{\theta}+\beta-\sin(\theta),$ (6)

where $\alpha=\frac{1}{\sqrt{qm}}$ and $\beta=\frac{\omega}{q}$. This equation
has two fixed point solutions, a saddle point and a stable node for $\beta<1$,
obtained by setting $\dot{\theta}=0$ and $\ddot{\theta}=0$. At $\beta=1$, the
system undergoes a saddle-node bifurcation, annihilating the two fixed point
solutions and admitting a unique stable limit cycle solution for all $\beta>1$
[29]. However, it so happens that as we decrease the value of $\beta$ to be
less than one, the limit cycle persists for some small values of $\alpha$.
Hence, bistability exists in the system, where a stable limit cycle and a
stable node coexist. A further decrease in $\beta$ will disintegrate the limit
cycle via a homoclinic bifurcation. Fig. 2a displays these three dynamical
regimes in the $\alpha-\beta$ parameter space. For small values of the damping
term $\alpha$, ensured by keeping finite inertia, the homoclinic bifurcation
curve is seen to be approximated by a straight line Fig. 2a. Upon implementing
Melnikov’s method, [28, 30] the equation of the straight line comes out to be
$\beta=\frac{4}{\pi}\alpha$. In conclusion, we see the presence of three
different dynamical regimes, namely a limit-cycle regime ($\beta>1$), a bi-
stable regime which can be approximated by ($\frac{4}{\pi}\alpha<\beta\leq
1$), and a fixed point regime approximated by ($\beta\leq\frac{4}{\pi}\alpha$)
[28].

The bi-stable region is responsible for hysteresis in systems governed by
equations like Eq. 6. Hence, following [9], instead of studying the system in
its full generality, we break down the self-consistency analysis for our model
into forward $\left(f\right)$ and backward $\left(b\right)$ processes. In the
forward process, we start from a small $K_{1}$ value, and therefore, the
system is in an incoherent state $\left(r_{1}\approx 0\right)$. This leads to
high $\alpha$ and $\beta$ values, indicating that the oscillators are in the
limit cycle regime. As we adiabatically increase $K_{1}$, the oscillators stay
in the basin of attraction of the stable limit cycle even after crossing
$\beta=1(\omega=q)$ and fall into the locked cluster only after
$\beta\approx\frac{4}{\pi}$
$\alpha(\omega\approx\frac{4}{\pi}\sqrt{\frac{q}{m}})$, below which the limit
cycle vanishes. For the backward process, we start from a high $K_{1}$ value,
and hence the oscillators exist in the fixed point regime, i.e., the
oscillators are locked in a cluster and therefore, the system is in a coherent
state $\left(0<<r_{1}<1\right)$. As we adiabatically decrease $K_{1}$, the
oscillators remain in the basin of attraction of the stable node until
$\beta=1$, when the stable node vanishes via a saddle-node bifurcation. Thus,
in the backward process, oscillators having $|\omega|\leq q=\omega_{b}$
contribute to the locked oscillators, while in the forward process, only those
with $|\omega|\leq\frac{4}{\pi}\sqrt{\frac{q}{m}}=\omega_{f}$ are in a locked
state and all the oscillators with $\omega>\omega_{f,b}$ drift around the
locked cluster. We point out that $K_{2}$ is concealed in $q$ and directly
affects the fraction of oscillators in a locked or drifting state.

The contribution of the locked oscillator($r_{p}^{l}$) to overall coherence
for the forward/backward process can now be calculated as
$r_{p}^{l}=\int_{-\omega_{f,b}}^{\omega_{f,b}}e^{ip\sin^{-1}(\frac{\omega}{q})}g(\omega)d\omega$.
The imaginary part of $r_{p}^{l}$ is zero as $g(-\omega)=g(\omega)$. Hence
taking only the real part and noting that
$\theta_{f,b}=\sin^{-1}(\omega_{f,b}/q)$, we arrive at the expression for
$r_{p}^{l}$ as follows,

$\begin{split}r_{p}^{l}=q\int_{-\theta_{f,b}}^{\theta_{f,b}}\cos(\theta)\cos(p\theta)g(q\sin(\theta))d\theta.\end{split}$
(7)

The contribution to overall coherence from the drifting oscillators can be
accounted for by calculating
$r_{p}^{d}=\int_{|\omega|>\omega_{f,b}}\int_{-\pi}^{\pi}e^{ip\theta}\rho_{d}(\theta,\omega)g(\omega)d\omega
d\theta$ where $\rho_{d}(\theta,\omega)$ is the density of drifting oscillator
which satisfies $\rho_{d}(\theta,\omega)\propto 1/|\dot{\theta}|$ [9]. The
normalization condition for $\rho_{d}(\theta,\omega)$ gives,
$\int_{-\pi}^{\pi}\rho_{d}(\theta,\omega)d\theta=\int_{0}^{T}\rho_{d}(\theta,\omega)\dot{\theta}dt=1$
(for a given $\omega$), where $T$ is the time period of the limit cycle
solution. Hence we end up with the relation
$\rho_{d}(\theta,\omega)=\frac{1}{|\dot{\theta}|T}$, which when plugged into
the form of $r_{p}^{d}$ gives us,

$\begin{split}r_{p}^{d}=\int_{|\omega|>\omega_{f,b}}\left[\frac{1}{T}\int_{0}^{T}e^{ip\theta}dt\right]g(\omega)d\omega.\end{split}$
(8)

To calculate $r_{p}^{d}$, we first need to obtain an approximate analytic
expression for the limit cycle solution of Eq. 6. We follow the method
specified in [31] of writing $\dot{\theta}$ as a Fourier series in $\theta$ by
only considering the first harmonics
$(\dot{\theta}(\theta)=A_{0}+A_{1}\cos(\theta)+B_{1}\sin(\theta))$. On
substituting this in Eq. 6, we find the expression of the coefficients in
terms of $\alpha$($=\frac{1}{\sqrt{qm}}$) and $\beta$($=\frac{\omega}{q}$)
such that the first harmonic vanishes, giving us,

$\dot{\theta}(\theta)=\frac{\beta}{\alpha}+\frac{\alpha^{2}}{\alpha^{4}+\beta^{2}}\left[\frac{\beta}{\alpha}\cos(\theta)-\alpha\sin(\theta)\right],$
(9)

and $\theta(t,\omega)$ by integrating Eq. 9 with time [31]. As it turns out
that $\theta(t,-\omega)=-\theta(t,\omega)$, and $g(-\omega)=g(\omega)$, the
imaginary part in Eq. 8 is zero. Thus,

$r_{p}^{d}=\int_{|\omega|>\omega_{f,b}}\left<\cos(p\theta)\right>g(\omega)d\omega.$
(10)

The expression for $\left<\cos(p\theta)\right>$ (for $p\in\\{1,2\\}$) can now
be readily calculated as
$\left<\cos(p\theta)\right>=\frac{1}{T}\int_{0}^{T}\cos(p\theta)dt=\int_{0}^{2\pi}\frac{\cos(p\theta)}{\dot{\theta}}d\theta\text{\LARGE$/$}\int_{0}^{2\pi}\frac{1}{\dot{\theta}}d\theta$
to obtain,

$\displaystyle\left<\cos(\theta)\right>$
$\displaystyle=\frac{\beta}{\alpha}\left[\sqrt{\frac{\beta^{2}}{\alpha^{2}}-\frac{\alpha^{2}}{\beta^{2}+\alpha^{4}}}-\frac{\beta}{\alpha}\right],$
(11) $\displaystyle\left<\cos(2\theta)\right>$
$\displaystyle=\left[\frac{\beta^{2}-\alpha^{4}}{\beta^{2}+\alpha^{4}}\right]\times$
$\displaystyle\left[\frac{2\beta(\beta^{2}+\alpha^{4})}{\alpha^{3}}\left(\frac{\beta}{\alpha}-\sqrt{\frac{\beta^{2}}{\alpha^{2}}-\frac{\alpha^{2}}{\beta^{2}+\alpha^{4}}}\right)-1\right].$

We are now finally ready to write down the self-consistent equations that let
us describe the steady state of the coupled oscillator system governed by Eq.
1. For the remainder of the work, we consider the intrinsic frequency to be
derived from Lorentz distribution,
$g(\omega)=\frac{1}{\pi}\frac{1}{1+\omega^{2}}$ centered around zero. Noting
that the integrands in Eqs. 7 and 10 for $p\in\\{1,2\\}$ are even functions,
we arrive at,

$\begin{split}r_{p}=2q\int_{0}^{\theta_{f,b}}\cos(\theta)\cos(p\theta)g(q\sin(\theta))d\theta\\\
+2\int_{\omega_{f,b}}^{\infty}\left<\cos(p\theta)\right>g(\omega)d\omega.\end{split}$
(12)

These two equations together describe the steady-state behavior of the system.

##  IV Numerical Results

We numerically solve the set of self-consistent equations, Eq. 12 to find the
nontrivial branch of solutions for the order parameter (both forward and
backward processes). Fig. 1a provides a schematic representation of the
synchronization profiles of our result in comparison to previously explored
models [2, 9, 24]. Fig. 1b presents analytical and simulation results for the
$r_{1}$ vs. $K_{1}$ curves for $(m,K_{2})=(1,1)$ and $(m,K_{2})=(3,7)$. As for
the simulation protocol, we simulate Eq. 3 on a network of $N=10^{4}$ nodes by
splitting it into a pair of first-order differential equations and integrating
them using the Runge-Kutta 4 algorithm (time-step 0.1). For a chosen value of
$m$ and $K_{2}$, we start with random initial conditions for
$\theta(\in[0,2\pi))$ and $\dot{\theta}(\in[-1,1])$ and $K_{1}=0$. We
adiabatically increase $K_{1}$ in steps of $\Delta K_{1}$ ($=0.1$, unless
specified otherwise) till $K_{1}=12$ is reached (forward), followed by an
adiabatic decrease till $K_{1}=0$ (backward). By adiabatic increase/decrease,
we imply that for every $K_{1}$ except the first ($K_{1}=0$), the initial
conditions are taken as the final state obtained for the previous $K_{1}$
value. At all coupling strengths $K_{1}$, the order parameter values are
calculated after discarding transients by averaging over the steady state.

Fig. 1b displays a good agreement between the simulation and analytical
results. For the forward process, as $K_{1}$ is increased from zero, the
system undergoes a first-order phase transition from incoherent to coherent
state at a finite critical coupling value ($K_{1}^{f}$). However, for the
backward process, the system undergoes abrupt desynchronization at a value
($K_{1}^{b}$), which is less than $K_{1}^{f}$. Hence, hysteresis is observed
where the system stays in two different states depending on the initial
configuration. The derived self-consistency equations can also be used with
other extended-tailed distributions like the Gaussian distribution. In the
backward process, there exists a discrepancy between analytical and numerical
values (Fig. 1b) as also reported in [9]. This happens because the maximum
value of the $K_{1}$ chosen for the simulations is 12, when the system starts
from a partially coherent state. Whereas, the analytical solutions are for
when the system starts from a fully coherent state. Hence, a better fit
between analytical and numerical values for the backward process can be
obtained by increasing the maximum value of $K_{1}$ in the simulation
protocol. We point out that when $m$ and $K_{2}$ values are both increased,
$K_{1}^{f}$ shifts to the right while $K_{1}^{b}$ shifts to the left,
revealing a prolonged hysteresis region as illustrated in Fig. 1b.

Figure 2: (Color online) a) $\alpha(=1/\sqrt{qm})$-$\beta(=\omega/q)$
parameter space. Different dynamical regimes are present in the $\dot{\theta}$
vs. $\theta$ phase space of Eq. 6. b) Synchronization profile $r_{1}$ versus
$K_{1}$ for $m=1$ and different values of $K_{2}=$ 0 (orange squares), 3 (blue
triangles), and 5 (red circles). c) Synchronization profile for a fixed value
of $K_{2}=2$ and different values of $m=$ 3 (orange squares), 1 (blue
triangles), and 0 (red circles). In both b) and c), the filled and hollow
symbols indicate the simulation results for the forward and backward cases,
respectively. The dashed and continuous curves represent the analytically
calculated values for the forward and backward processes, respectively. d)
Backward transition points. The dashed curve represents the analytical
predictions of $K_{1}^{b}$ for $m=0$ and different $K_{2}$
($K_{1}^{b}=2\sqrt{2K_{2}}-K_{2}$ as derived in Section VI.1 for $K_{2}\geq
2$). The scatter plots are the $K_{1}^{b}$ vs. $K_{2}$ for different $m$ = 0,
1, 5, 10 obtained via numerical simulation.

A natural question would then be to address the dependency of the forward and
backward transition points on $m$ and $K_{2}$. To analytically obtain the
expression for ($K_{1}^{f}$), we evaluate Eq. 12 in the limit
$r_{1}\rightarrow 0^{+}$ ($q\rightarrow 0^{+}$). As we take this limit, we see
that $\beta/\alpha$($=\omega\sqrt{\frac{m}{q}}$) tends to very high value as
compared to
$\alpha^{2}/(\beta^{2}+\alpha^{4})$($=\frac{qm}{1+\omega^{2}m^{2}}$). This
allows us to perform a Taylor series expansion of $\left<\cos(\theta)\right>$
for $\epsilon=\alpha^{2}/(\beta^{2}+\alpha^{4})<<1$ which gives,
$\left<\cos(\theta)\right>=\frac{-\alpha^{2}}{2(\beta^{2}+\alpha^{4})}+\mathcal{O}(\epsilon^{4})\approx\frac{-mq}{2(1+m^{2}\omega^{2})}$.
However, in the limit $r_{1}\rightarrow 0^{+}$, $r_{2}\rightarrow 0^{+}$ and
the parameter $\alpha\rightarrow\infty$ implying that the limit of the
integrals for the forward and backward processes become the same as there
exists no bistability region in the parameter space. Taking
$\theta_{f,b}=\frac{\pi}{2}$, dividing both sides of Eq. 12 by $q$, and
evaluating the limit (at which the two equations in Eq. 12 decouple) we have,

$\frac{1}{K_{1}^{f}}=\frac{\pi}{2}g(0)-m\int_{0}^{\infty}\frac{1}{1+m^{2}\omega^{2}}g(\omega)d\omega.$
(13)

After evaluating the integral and rearranging the terms, we end up with
$K_{1}^{f}=2(m+1)$. We see that the forward transition point is independent of
$K_{2}$ and purely depends on $m$ and, hence, matches the previously derived
value of the forward transition point in [33, 7]. Fig. 2b illustrates the
effect of varying $K_{2}(0.0,3.0,5.0)$ for the case of fixed $m(=1)$. As
expected, $K_{1}^{f}$ remains the same for all three cases, validating our
analytical result. At this $K_{1}^{f}(=4)$, the magnitude of the first-order
jump for fixed $m$ increases with the value of $K_{2}$. In Fig. 2c, we study
the effect of varying mass ($0.0,1.0,3.0$) for the case of fixed
$K_{2}(=2.0)$. As inertia increases, $K_{1}^{f}$ shifts to higher values as
predicted analytically. However, we note that the analytically calculated
values of $K_{1}^{f}$ do not match exactly with numerical simulations owing to
the finite size effects. A detailed study has been done in [32].

A fairly good analytical approximation for $K_{1}^{b}$, as also pointed out in
[32], would be to obtain the minimum value of $K_{1}$ along the non-trivial
branch of the backward self-consistent curve. The simulation results in Fig.
1b and Fig. 2b and 2c are seen to back up this observation for our model.
However, obtaining a clean analytical expression for the same by calculating
$\frac{dK_{1}}{dr_{1}}=0$ is not possible because of the complexity of the
integrand of the drift oscillator contribution in $r_{2}$. Alternatively, we
resort to simulation results to decipher the dependency of $K_{1}^{b}$ on $m$
and $K_{2}$. From Fig. 2b, it can be seen that for the backward process, the
coherent branch persists till increasingly smaller values of $K_{1}$ with an
increase in the $K_{2}$ value, after which the system undergoes an abrupt
transition to the incoherent state. Hence, it is clear that an increase in
$K_{2}$ leads to a decrease in $K_{1}^{b}$. To study the effect of mass on
$K_{1}^{b}$, we fix $K_{2}$ and vary $m$ as in Fig. 2c. It is observed that
the backward branches for fixed $K_{2}$(=2) and different $m$ values merge for
high $K_{1}$ values and get separated for low $K_{1}$ values. As there is an
influence of $m$ on the nature of the curve for low $K_{1}$ values, this
indicates the possibility of dependency of $K_{1}^{b}$ on $m$.

It was shown in [32] that for the pure dyadic case ($K_{2}=0$), $K_{1}^{b}$
decreases with an increase in $m$ and plateaus out for high $m$ values. In
Fig. 2d, we address how this changes with the introduction of finite $K_{2}$.
The $K_{1}^{b}$ obtained via simulation (performed for $N=10^{3}$ number of
nodes) for values of $K_{2}$ ranging from 0 to 10 and different values of
$m$(0,1,5,10) are plotted. We see that for small values of $K_{2}$ and finite
inertia, an increase in the values in $m$ leads to a decrease in $K_{1}^{b}$.
However, we point out that for higher values of $K_{2}$, the effect of $m$ on
$K_{1}^{b}$ becomes less pronounced, and desynchronization happens at the same
value irrespective of mass. An analytical prediction of $K_{1}^{b}$ becomes
possible following this observation by considering the $m=0$ case. We derive
self-consistent equations for this case in Section VI.1 and obtain $K_{1}^{b}$
as a function of $K_{2}$ ($K_{1}^{b}=2\sqrt{2K_{2}}-K_{2}$ for $K_{2}\geq 2$)
by finding the minimum value of $K_{1}$ in the self-consistency curve. These
analytically calculated $K_{1}^{b}$ values for the $m=0$ case are represented
by the dashed line in Fig. 2d. It can be clearly observed that for higher
values of $K_{2}(\geq 2)$, the analytical predictions of $K_{1}^{b}$ match
closely with the ones obtained via simulation for different masses.

##  V Conclusion and Discussion

In this article, we have put forward a generalized analytical framework to
study the steady-state behavior of coupled oscillator systems with inertia
interacting via higher-order interactions. The analytical predictions, backed
up by numerical simulation, show a prolonged hysteretic first-order phase
transition to an (in)coherent state. We show that the forward transition point
increases linearly with $m$ and is independent of $K_{2}$. Meanwhile, the
backward transition point decreases with $K_{2}$ and is independent of $m$ for
high $K_{2}$ values.

We also highlight here the analytical challenges emanating from combining both
the inertia and higher-order interaction terms in the Kuramoto model. For the
coupled Kuramoto oscillator model, two widespread analytical techniques are
the self-consistency method and the Ott-Antonsen (OA) approach to studying the
synchronization profile of the entire system. OA dimensionality reduction
method provides a very easy way to obtain the time evolution form of $r_{1}$,
solving which yields a complete description of bifurcations leading to
synchronization. While the OA approach was successfully extended for the
Kuramoto model with higher-order interactions, for the second-order Kuramoto
models, the density function containing phase term also depends on velocity,
posing restrictions on practicing the OA method. Hence, the Kuramoto model
with inertia has been analytically solved using the self-consistency method.
The challenge in this method is first finding the limit of integration in the
self-consistent equation by finding an approximate frequency bound for the
oscillators, which are in the locked state as a function of external
parameters and the order parameter. Further, the drifting oscillators also
contribute to the overall coherence, which needs to be accounted for, unlike
in the first-order Kuramoto model, where it is zero. Here, we have employed
the self-consistency method to predict the steady-state behavior of the
Kuramoto oscillators having both inertia and higher-order terms.

We have presented the results for triadic interactions; however, extending our
analysis to other powers of higher-order interactions is easy, as long as the
sinusoidal coupling function contains $\theta_{i}$ term only. For example, the
detailed analysis of quartic interactions is presented in Section VI.2.
Further, developing the self-consistent method for other choices of higher-
order coupling functions, such as $\sin(\theta_{j}+\theta_{k}-2\theta_{i})$
[34, 35] along with pairwise coupling proves to be complicated because of the
existence of higher order harmonics in the mean-field equation; however, the
self-consistency approach can work for this form of the coupling in the
absence of pairwise coupling (pure triadic case) which we have explored in
Section VI.3. An immediate future direction of our work would be to extend our
analysis to diluted simplicial complexes, which can provide fundamental
insights into the dynamics of various real-world complex systems such as power
grids.

###### Acknowledgements.

SJ gratefully acknowledges SERB Power grant SPF/2021/000136, and thanks
Mehrnaz Anvari and Baruch Barzel for useful comments on the manuscript.

##  VI Appendix

###  VI.1 Derivation for $m=0$

This section elaborates on the self-consistency analysis and the derivation
for the forward and backward transition points for $m=0$. In this case, a
change from a smooth (second-order) transition to synchronization to an abrupt
(first-order) one, along with hysteresis, is observed with an increase in the
$K_{2}$ value. The occurrence of hysteresis can be accounted for by the shift
of the backward transition point to lower $K_{1}$ values with an increase in
$K_{2}$. Note that the results for the first-order Kuramoto model with higher-
order interactions were reported in [24] using the Ott-Antonsen dimensionality
reduction method. Here, we derive the closed-form solution for the bifurcation
curve using the self-consistency method.

The dynamical equation for our model takes the following form.

$\dot{\theta}_{i}=\omega_{i}+\frac{K_{1}}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i})+\frac{K_{2}}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}\sin(2\theta_{j}-\theta_{k}-\theta_{i})$
(A-1)

Following the same procedure as in the article, the mean-field equation in the
rotating frame is obtained as $\dot{\theta_{i}}=\omega_{i}-q\sin(\theta_{i})$,
where $q=r_{1}(K_{1}+K_{2}r_{2})$. All the oscillators with intrinsic
frequency $|\omega_{i}|\leq q$ go to the fixed point state, while others are
in a drift state. Unlike in the finite-inertia case, here the limit of the
integrals for forward and backward cases become the same
($\theta_{f}=\theta_{b}=\frac{\pi}{2}$), enabling us to study the system in
full generality. In the thermodynamic limit, the locked oscillator
contribution has the same form as derived in the article,

$\begin{split}r_{p}^{l}=q\int_{-\pi/2}^{\pi/2}\cos(\theta)\cos(p\theta)g(q\sin(\theta))d\theta,\end{split}$
(A-2)

which upon integrating for $p\in\\{1,2\\}$ gives us,

$\displaystyle r_{1}^{l}$ $\displaystyle=\dfrac{\sqrt{q^{2}+1}-1}{q},$ (A-3)
$\displaystyle r_{2}^{l}$ $\displaystyle=\frac{2}{\pi
q^{2}}\left[(q^{2}+2)\tan^{-1}(q)-2q\right].$

The contribution of drifting oscillators to the global order parameter,
however, has the following expression,

Figure A-1: Synchronization profile of Eq. A-1, for and $K_{2}=0$ (cyan), $3$
(green), $6$ (red), and $9$ (yellow). The circles represent simulation results
for the forward and backward cases, while the continuous line represents the
analytically predicted values.

$\displaystyle r_{p}^{d}$
$\displaystyle=\int_{|\omega|>q}\int_{-\pi}^{\pi}e^{ip\theta}\rho_{d}(\theta,\omega)g(\omega)d\omega
d\theta,$ $\displaystyle r_{p}^{d}$
$\displaystyle=\int_{|\omega|>q}\int_{-\pi}^{\pi}e^{ip\theta}\frac{\sqrt{\omega^{2}-q^{2}}}{2\pi|\omega-q\sin(\theta)|}g(\omega)d\omega
d\theta.$

The expression for the density of the drifting oscillators is derived by
noting that $\rho_{d}(\theta,\omega)\propto
1/|\dot{\theta}|=c/|\omega-q\sin(\theta)|$, where $c$ is the normalization
constant such that $\int_{-\pi}^{\pi}\rho_{d}(\theta,\omega)d\theta=1$. The
value of $r_{1}^{d}$ is zero since
$\rho_{d}(\theta+\pi,-\omega)=\rho_{d}(\theta,\omega)$ and
$g(-\omega)=g(\omega)$. Meanwhile, in $r_{2}^{d}$, only the imaginary term is
zero, since $\rho_{d}(-\theta,-\omega)=\rho_{d}(\theta,\omega)$ and
$g(-\omega)=g(\omega)$ while the real part is not zero. Since the integrand is
even, we can evaluate the integral as follows:

$\displaystyle r_{2}^{d}$
$\displaystyle=2\int_{q}^{\infty}\left[\int_{-\pi}^{\pi}\cos(2\theta)\frac{\sqrt{\omega^{2}-q^{2}}}{2\pi(\omega-q\sin(\theta))}d\theta\right]g(\omega)d\omega,$
(A-4)
$\displaystyle=2\int_{q}^{\infty}\left[\dfrac{2\omega}{q^{2}}\left(\sqrt{\omega^{2}-q^{2}}-\omega\right)+1\right]g(\omega)d\omega,$
$\displaystyle=-\dfrac{2}{\pi
q^{2}}\left((q^{2}+2)\tan^{-1}(q)-2q\right)+\dfrac{2+q^{2}-2\sqrt{q^{2}+1}}{q^{2}}.$

Using the self-consistency condition $(r_{p}=r_{p}^{l}+r_{p}^{d})$ and Eq.
A-3, Eq. A-4 we arrive at a set of self-consistent equations which when solved
give us the analytical predictions of the steady state behavior,

$\displaystyle r_{1}$ $\displaystyle=\dfrac{\sqrt{q^{2}+1}-1}{q},$ (A-5a)
$\displaystyle r_{2}$ $\displaystyle=\dfrac{2+q^{2}-2\sqrt{q^{2}+1}}{q^{2}},$
(A-5b)

where $q=r_{1}(K_{1}+r_{2}K_{2})$. Solving Eq. A-5 gives us the relation
$r_{2}=r_{1}^{2}$, which, when plugged back into Eq. A-5a leads to the closed
form solution for the bifurcation curve,

$K_{1}=\dfrac{2}{1-r_{1}^{2}}-r_{1}^{2}K_{2}.$ (A-6)

We point out that this equation matches with the one obtained by [24] using
the OA method.

Fig. A-1 presents analytical and simulation results for multiple values of
triadic coupling strength $K_{2}$. As predicted by our analysis, the forward
synchronization transition happens at $K_{1}^{f}=2$, while the backward
transition point shifts to lower values of $K_{1}$ as $K_{2}$ increases. From
Fig. A-1, we notice that the minima of the self-consistency curve fairly well
approximates the backward de-synchronization point. To this end, setting
$\frac{dK_{1}}{dr_{1}}=0$, we obtain
$r_{1}^{b}=\sqrt{1-\sqrt{\frac{2}{K_{2}}}}$ for $K_{2}\geq 2$. Plugging the
form of $r_{1}^{b}$ back into Eq. A-6 gives $K_{1}^{b}=2\sqrt{2K_{2}}-K_{2}$.
This is used in Fig. 2d in the article.

###  VI.2 Extension to quartic interactions

Figure B-1: Synchronization profiles for the Kuramoto model with inertia
involving quartic interactions with $m=1$ and different combinations of
$K_{2}$ and $K_{3}$ values. The red circles and blue squares represent the
simulation results for the forward and backward cases of Eq. B-1 for $N=1000$
nodes, respectively. Whereas the dashed and continuous curves are the
analytical predictions from B-2 for the forward and backward cases.

This section explains how the self-consistency analysis presented in the
article can be easily extended to include quartic interactions. We add a
quartic interaction term to our model as proposed in [24]. In our study, we
find a good agreement between the analytical predictions of our self-
consistency analysis and the simulation results. The detailed analysis of the
impact of quartic interactions on the dynamics is out of the scope of this
discussion. The dynamical equation for our model is given as follows:

$\displaystyle m\ddot{\theta_{i}}=$
$\displaystyle-\dot{\theta_{i}}+\omega_{i}+\frac{K_{1}}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i})$
(B-1)
$\displaystyle+\frac{K_{2}}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}\sin(2\theta_{j}-\theta_{k}-\theta_{i})$
$\displaystyle+\frac{K_{3}}{N^{3}}\sum_{j=1}^{N}\sum_{k=1}^{N}\sum_{l=1}^{N}\sin(\theta_{j}+\theta_{k}-\theta_{l}-\theta_{i}).$

Upon using the definition of the generalized order parameter as defined in the
article, we can convert the above equation into a mean-field form as,

$\displaystyle m\ddot{\theta}_{i}$
$\displaystyle=-\dot{\theta}_{i}+\omega_{i}+K_{1}r_{1}\sin(\psi_{1}-\theta_{i})$
$\displaystyle+K_{2}r_{1}r_{2}\sin(\psi_{2}-\psi_{1}-\theta_{i})+K_{3}r_{1}^{3}\sin(\psi_{1}-\theta_{i}).$

By moving to the rotating frame, we set $\psi_{1}$ and $\psi_{2}$ as zero,
which gives
$m\ddot{\theta}_{i}=-\dot{\theta}_{i}+\omega_{i}-q\sin(\theta_{i})$, where
$q=r_{1}(K_{1}+K_{2}r_{2}+K_{3}r_{1}^{2})$ is the overall coupling constant.
We proceed exactly, as in the article, to arrive at the self-consistency
equations defining the steady-state behavior of the model.

$\begin{split}r_{1}=2q\int_{0}^{\theta_{f,b}}\cos^{2}(\theta)g(q\sin(\theta))d\theta+2\int_{\omega_{f,b}}^{\infty}\omega\sqrt{\frac{m}{q}}\left[\sqrt{\omega^{2}\frac{m}{q}-\frac{qm}{1+\omega^{2}m^{2}}}-\omega\sqrt{\frac{m}{q}}\right]g(\omega)d\omega,\end{split}$
(B-2a)
$\begin{split}r_{2}&=2q\int_{0}^{\theta_{f,b}}\cos(\theta)\cos(2\theta)g(q\sin(\theta))d\theta+\\\
&2\int_{\omega_{f,b}}^{\infty}\left[\frac{\omega^{2}m^{2}-1}{\omega^{2}m^{2}+1}\right]\left[\frac{2\omega}{q}\left(\frac{\omega^{2}m^{2}+1}{\sqrt{qm}}\right)\left(\sqrt{\omega^{2}\frac{m}{q}-\frac{qm}{1+\omega^{2}m^{2}}}-\omega\sqrt{\frac{m}{q}}\right)-1\right]g(\omega)d\omega.\end{split}$
(B-2b)

Fig. B-1 plots $r_{1}$ as a function of $K_{1}$ for different combinations of
$K_{2}$ and $K_{3}$, and we see that the analytical predictions match the
numerical simulation results.

###  VI.3 Derivation for the $\sin(\theta_{j}+\theta_{k}-2\theta_{i})$ model

In this section, we present the self-consistency analysis for the Kuramoto
model with inertia involving purely triadic interactions via the
$\sin(\theta_{j}+\theta_{k}-2\theta_{i})$ coupling function. Research on the
dynamics of Kuramoto oscillators coupled via the said sinusoidal function
($\dot{\theta}_{i}=\omega_{i}+\frac{K_{2}}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}\sin(\theta_{j}+\theta_{k}-2\theta_{i})$)
[35] has shown the presence of cluster formation and a continuum of abrupt de-
synchronization transitions based on initial conditions. Crucially, no
synchronization transition has been reported for this model. Surprisingly, our
studies report that inertia has no effect on the synchronization profile of
this system. We infer that considering purely triadic interactions removes the
distinction between inertia-less and finite-inertia cases. The dynamical
equations for this model are given as:

$m\ddot{\theta_{i}}=\omega_{i}-\dot{\theta_{i}}+\frac{K_{2}}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}\sin(\theta_{j}+\theta_{k}-2\theta_{i}).$
(C-1)

Using the definition of the generalized order parameter, we write Eq. C-1 in
the mean-field format as,

$m\ddot{\theta_{i}}=\omega_{i}-\dot{\theta_{i}}+K_{2}r_{1}^{2}\sin(2\psi_{1}-2\theta_{i}).$
(C-2)

Upon simulating the dynamics of this equation, we observe that analogous to
the studies reported in [35], there exists no forward synchronization in the
system rather, a sequence of de-synchronization transitions is observed based
on the asymmetry in the initial conditions. Thus, we focus the analytical
derivation only on the de-synchronization profiles.

Dropping the subscript $i$, by going to a suitable rotating frame, we set
$\psi=0$, which gives us
$m\ddot{\theta}=\omega-\dot{\theta}-K_{2}r_{1}^{2}\sin(2\theta)$, where the
probability distribution g($\omega$) is unimodal and symmetric about the mean
zero. To study the steady-state behavior, we invoke a time transformation as
$\tau=\sqrt{\frac{K_{2}r_{1}^{2}}{m}}t$, which gives:

$\ddot{\theta}=\beta-\alpha\dot{\theta}-\sin(2\theta),$ (C-3)

where $\beta=\frac{\omega}{K_{2}r_{1}^{2}}$ and
$\alpha=\frac{1}{\sqrt{K_{2}r_{1}^{2}m}}$. The parameter space of Eq. C-3 is
qualitatively similar to that of the model considered in the article. The
quantitative differences between the two are as follows: (i) For each $\beta$
such that $|\beta|\leq 1$, we now have two stable fixed points given by
$\theta_{1}^{*}=\frac{1}{2}\sin^{-1}{\beta},\,\,\,\text{and}\,\,\,\theta_{2}^{*}=\frac{1}{2}\sin^{-1}{\beta}+\pi$
separated by two saddles. (ii) The separatrix equation is given by
$\theta(t)=\sin^{-1}{\tanh(\sqrt{2}t)},\,\,\,\dot{\theta}(t)=\sqrt{2}\frac{1}{\cosh(\sqrt{2}t)}$.
Using this equation of separatrix and implementing Melnikov’s method [30], the
equation for homoclinic bifurcation can be obtained as
$\beta=\frac{2\sqrt{2}}{\pi}\alpha$. Thus, as seen in Fig. C-1

Figure C-1: The parameter space for Eq. C-3 representing the limit-cycle,
tristable, and fixed-point regimes, respectively. The line at $\beta=1$
represents the saddle-node bifurcation in the system. On the boundary of the
tri-stable and fixed-point regimes, the scatter plot indicates the actual
homoclinic bifurcation curve while the dotted line is the approximation
obtained using Melnikov’s method ($\beta=\frac{2\sqrt{2}}{\pi}\alpha$).

  1. 1.

if $|\beta|>1$, the system goes to a limit cycle;

  2. 2.

if $|\beta|\leq\frac{2\sqrt{2}}{\pi}\alpha$, the system goes to a stable fixed
point state, where the choice of the fixed point depends upon the initial
conditions.

  3. 3.

Finally, if $\frac{2\sqrt{2}}{\pi}\alpha<|\beta|<1$, we now have a tri-stable
state with the simultaneous presence of one stable limit cycle and two stable
fixed points.

Figure C-2:  The top panel displays the $r_{1}$ versus $K_{2}$ behaviour for
Eq. C-2 for a) $m=0$, b) $m=1$, and c) $m=3$ respectively. The violet, yellow,
and green curves represent the de-synchronization profiles for $\eta$ equal to
1, 0.9, and 0.8, respectively. For each case, the continuous curve represents
the analytically obtained values. The bottom panel displays the simulation-
obtained values for the Daido order parameter ($r_{2}$) for the corresponding
cases (d) $m=0$, e) $m=1$, and f) $m=3$.

In terms of $\omega$, $K_{2}$, and $r_{1}$, it means that for a particular
$K_{2}$ and corresponding steady-state value of $r_{1}$, the oscillators with
$|\omega|>K_{2}r_{1}^{2}$ become drifting oscillators. The oscillators with
$|\omega|\leq\frac{2\sqrt{2}}{\pi}\sqrt{\frac{K_{2}r_{1}^{2}}{m}}$ contribute
to the formation of two diametrically opposite clusters of locked oscillators.
Finally the oscillators with
$\frac{2\sqrt{2}}{\pi}\sqrt{\frac{K_{2}r_{1}^{2}}{m}}<|\omega|\leq 1$ are in
the tri-stable stable region. However, as we will deal only with the system
de-synchronization profile, the oscillators in the tri-stable region will also
go to their respective fixed points and contribute to cluster formation.

The next step is to calculate the locked and drifting oscillator contribution
to the order parameter. First, let us calculate the contribution from the
locked oscillators. To account for the two clusters of locked oscillators, we
introduce a variable $\eta(\omega)\in[0,1]$. The values of $\eta(\omega)$ and
$1-\eta(\omega)$ are the probabilities that an oscillator with intrinsic
frequency $\omega$ will go to the fixed points
$\theta_{1}^{*}=\frac{1}{2}\sin^{-1}{\beta}$ and
$\theta_{2}^{*}=\frac{1}{2}\sin^{-1}{\beta}+\pi$, respectively. For
simplicity, we consider only symmetric cases of the function $\eta(\omega)$,
such that $\eta(\omega)=\eta(-\omega)$. The contribution of locked oscillators
is then given by
$r_{1}^{l}=\int_{-K_{2}r_{1}^{2}}^{K_{2}r_{1}^{2}}[(1-\eta(\omega))e^{i\theta(\omega)+\pi}+\eta(\omega)e^{i\theta(\omega)}]g(\omega)d\omega$.
As $e^{i\theta(\omega)+\pi}=-e^{i\theta(\omega)}$ and $\sin(-x)=-\sin(x)$, we
have:

$r_{1}^{l}=\int_{-K_{2}r_{1}^{2}}^{K_{2}r_{1}^{2}}(2\eta(\omega)-1)\cos\theta(\omega)g(\omega)d\omega.$
(C-4)

Next, let us consider the contribution of the drifting oscillators to the
order parameter. Let $\rho_{d}(\theta,\omega)$ be the density of drifting
oscillators which satisfies
$\int_{-\pi}^{\pi}\rho_{d}(\theta,\omega)d\theta=1$. The continuity equation
for the conservation of the number of oscillators gives
$\rho_{d}(\theta,\omega)=c/\dot{\theta}(\theta,\omega)$. An expression for
$\dot{\theta}$ can be obtained by following an analogous method as in the
article by considering
$\dot{\theta}=A_{0}+A_{1}\cos(2\theta)+B_{1}\sin(2\theta)$. By substituting
this in Eq. C-3, the expression for $\dot{\theta}$ can be obtained as:

$\dot{\theta}(\omega,\theta)=\frac{\beta}{\alpha}+\frac{2\beta\alpha}{\alpha^{4}+4\beta^{2}}\cos(2\theta)-\frac{\beta^{3}}{\alpha^{4}+4\beta^{2}}\sin(2\theta).$
(C-5)

Eq. C-5 can be simplified by considering $h_{0}=\beta/\alpha$ and
$\frac{e^{ih_{2}}}{h_{1}}=2h_{0}+i\alpha$, which gives
$\dot{\theta}=h_{0}+h_{1}\cos(2\theta+h_{2})$. By integrating over the
normalization condition of $\rho_{d}(\theta,\omega)$, we get
$c=\sqrt{h_{0}^{2}-h_{1}^{2}}/2\pi$, which gives the following equation for
the density of drifting oscillators:

$\rho_{d}(\theta,\omega)=\left|\frac{1}{2\pi}\sqrt{\frac{\omega^{2}m}{K_{2}r_{1}^{2}}-\frac{K_{2}r_{1}^{2}m}{1+4m^{2}\omega^{2}}}\times\frac{1}{\frac{\omega\sqrt{m}}{\sqrt{K_{2}r_{1}^{2}}}+\sqrt{\frac{K_{2}r_{1}^{2}m}{1+4m^{2}\omega^{2}}}\cos(2\theta(\omega)+\tan^{-1}(\frac{1}{2m\omega}))}\right|,$
(C-6)
$r_{1}^{d}=\int_{|\omega|>K_{2}r_{1}^{2}}\int_{-\pi}^{0}[e^{i\theta(\omega)}\rho_{d}(\theta,\omega)+e^{i(\theta(\omega)+\pi)}\rho_{d}(\theta+\pi,\omega)]g(\omega)d\theta
d\omega.$ (C-7)

From Eq. C-6, we notice that
$\rho_{d}(\theta+\pi,\omega)=\rho_{d}(\theta,\omega)$. Thus, the contribution
of drifting oscillators to order parameter
$r_{1}^{d}=\int_{|\omega|>K_{2}r_{1}^{2}}\int_{-\pi}^{\pi}e^{i\theta(\omega)}\rho_{d}(\theta,\omega)g(\omega)d\theta
d\omega$ can be simplified as Eq. C-7

As $e^{i(\theta(\omega)+\pi)}=-e^{i\theta(\omega)}$, we get $r_{1}^{d}=0$.
Therefore, from Eq. C-4 and the self consistency condition
($r_{1}=r_{1}^{l}+r_{1}^{d}$), we get

$r_{1}=\int_{-K_{2}r_{1}^{2}}^{K_{2}r_{1}^{2}}(2\eta(\omega)-1)\cos\theta(\omega)g(\omega)d\omega.$
(C-8)

We now compare the results obtained by simulating Eq. C-2 with the analytical
values predicted by the Eq. C-8 in Fig. C-2. For simplicity, we consider
$\eta(\omega)$ to be a constant function with respect to $\omega$. While
simulating for a given value of $\eta$, we consider the initial phase of the
oscillators as $0$ with probability $\eta$ and $\pi$ with probability
$(1-\eta)$. With these initial conditions, we start from $K_{2}=K_{max}=16$,
and then adiabatically decrease $K_{2}$ till $K_{min}=0$. At each $K_{2}$, we
integrate Eq. C-2 using the RK4 algorithm and calculate the value of $r_{1}$
and $r_{2}$ by averaging over all the iterations after removing the
transients.

We note that due to the vanishing of the drift term, the form of the
bifurcation curve Eq. C-8 is the same as well for the $m=0$ case. Fig. C-2
also shows that the $r_{1}$ and $r_{2}$ graphs are almost identical for the
finite mass and the mass-less case, with identical backward transition points.
Hence, considering purely triadic interactions eliminates the distinction
between the finite inertia and inertia-less case. For $\eta=0.9$ (Yellow) and
$\eta=0.8$ (Green), we see that the $r_{2}$ values are higher than the $r_{1}$
values, indicating that a two-cluster state is present in the system, which is
the expected outcome.

A more generalized approach to further this study would be to incorporate
pair-wise interactions along with the triadic ones to give a model like:

$\displaystyle m\ddot{\theta_{i}}$
$\displaystyle=\omega_{i}-\dot{\theta_{i}}+\frac{K_{1}}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i})$
(C-9)
$\displaystyle+\frac{K_{2}}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}\sin(\theta_{j}+\theta_{k}-2\theta_{i})$

which in the mean-field format and rotating frame, with $\psi_{1}$ set to
zero, reads as
$m\ddot{\theta_{i}}=\omega_{i}-\dot{\theta_{i}}-K_{1}r_{1}\sin(\theta_{i})-K_{2}r_{1}^{2}\sin(2\theta_{i})$.
However, the simultaneous presence of the first and second-order harmonics of
the sinusoidal term makes it complicated to analytically study the parameter
space of Eq. C-9.

## References

  * [1] Boccaletti, S., Latora, V., Moreno, Y., Chavez, M., & Hwang, D. U. Complex networks: Structure and dynamics. Physics reports 424 (2006).
  * [2] Kuramoto, Y. Self-entrainment of a population of coupled non-linear oscillators.( Springer. Berlin Heidelberg. 1975). International Symposium on Mathematical Problems in Theoretical Physics: January 23–29, 1975, Kyoto University (pp. 420-422)
  * [3] Sethia, Gautam C., Abhijit Sen, and Fatihcan M. Atay. Clustered chimera states in delay-coupled oscillator systems. Physical Review Letters 100, 144102 (2008).
  * [4] Childs, L. M., & Strogatz, S. H. Stability diagram for the forced Kuramoto model. Chaos: An Interdisciplinary Journal of Nonlinear Science, 18 (4), 043128 (2008).
  * [5] Omel’Chenko, E., & Wolfrum, M. Nonuniversal transitions to synchrony in the Sakaguchi-Kuramoto model. Physical Review Letters, 109 (16), 164101 (2012).
  * [6] Olmi, S. Chimera states in coupled Kuramoto oscillators with inertia. Chaos: An Interdisciplinary Journal of Nonlinear Science, 25 (12), 123125 (2015).
  * [7] Rodrigues, F. A., Peron, T. K. D., Ji, P., & Kurths, J. (2016). The Kuramoto model in complex networks. Physics Reports, 610, 1-98.
  * [8] Ermentrout, B. An adaptive model for synchrony in the firefly Pteroptyx malaccae. Journal of Mathematical Biology, 29 (6) (1991).
  * [9] Tanaka, H. A., Lichtenberg, A. J., & Oishi, S. I. First order phase transition resulting from finite inertia in coupled oscillator systems. Physical Review Letters, 78 (11), 2104 (1997).
  * [10] Tanaka, H. A., Lichtenberg, A. J., & Oishi, S. I. Self-synchronization of coupled oscillators with hysteretic responses. Physica D: Nonlinear Phenomena, 100 (3-4) (1997).
  * [11] Ji, P., Peron, T. K. D., Menck, P. J., Rodrigues, F. A., & Kurths, J. Cluster explosive synchronization in complex networks. Physical Review Letters, 110 (21), 218701 (2013).
  * [12] Trees, B. R., Saranathan, V., & Stroud, D. Synchronization in disordered Josephson junction arrays: Small-world connections and the Kuramoto model. Physical Review E, 71(1), 016215 (2005).
  * [13] Ji, P., & Kurths, J. Basin stability of the Kuramoto-like model in small networks. The European Physical Journal Special Topics, 223(12) (2014).
  * [14] Grzybowski, J. M. V., Macau, E. E. N., & Yoneyama, T. On synchronization in power-grids modelled as networks of second-order Kuramoto oscillators. Chaos: An Interdisciplinary Journal of Nonlinear Science, 26(11), 113113 (2016).
  * [15] Rohden, M., Sorge, A., Timme, M., & Witthaut, D. Self-organized synchronization in decentralized power grids. Physical Review Letters, 109 (6), 064101 (2012).
  * [16] Dörfler, F., Chertkov, M., & Bullo, F. Synchronization in complex oscillator networks and smart grids. Proceedings of the National Academy of Sciences, 110 (6) (2013).
  * [17] Witthaut, D., Hellmann, F., Kurths, J., Kettemann, S., Meyer-Ortmanns, H., & Timme, M. Collective nonlinear dynamics and self-organization in decentralized power grids.Reviews of Modern Physics, 94(1), 015005 (2022).
  * [18] Filatrella, G., Nielsen, A. H., & Pedersen, N. F. Analysis of a power grid using a Kuramoto-like model. The European Physical Journal B, 61 (2008).
  * [19] Benson, A. R., Abebe, R., Schaub, M. T., Jadbabaie, A., & Kleinberg, J. Simplicial closure and higher-order link prediction. Proceedings of the National Academy of Sciences, 115 (48) (2018).
  * [20] Boccaletti, S. et. al. The structure and dynamics of networks with higher order interactions. Physics Reports, 1018 1 (2023).
  * [21] Battiston, F., Cencetti, G., Iacopini, I., Latora, V., Lucas, M., Patania, A., & Petri, G. Networks beyond pairwise interactions: structure and dynamics. Physics Reports, 874, (2020).
  * [22] Rajwani, Suman, A., & P., Jalan, S. Tiered synchronization in Kuramoto oscillators with adaptive higher-order interactions. Chaos (Fast Track), 33 , (2023).
  * [23] Ott, E., & Antonsen, Thomas M. Low dimensional behavior of large systems of globally coupled oscillators. Chaos. 18, 3 (2008).
  * [24] Skardal, P. S., & Arenas, A. Higher order interactions in complex networks of phase oscillators promote abrupt synchronization switching. Communications Physics, 3 (1), 218 (2020).
  * [25] León, I., & Pazó, D. Phase reduction beyond the first order: The case of the mean-field complex Ginzburg-Landau equation. Physical Review E, 100(1), (2019).
  * [26] Ji, P., Peron, Thomas K. D. M., Rodrigues, F. A., & Kurths, J. Low-dimensional behavior of Kuramoto model with inertia in complex networks. Scientific Reports, 4 4783 (2014).
  * [27] Daido, H. Multibranch entrainment and scaling in large populations of coupled oscillators. Physical Review Letters, 77(7), 1406 (1996).
  * [28] Strogatz, S. H. Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering.(CRC press. 2018).
  * [29] Levi, M., Hoppensteadt, F. C., & Miranker, W. L. Dynamics of the Josephson junction. Quarterly of Applied Mathematics, 36 (2) (1978).
  * [30] Guckenheimer, J., & Holmes, P. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. (Springer Science & Business Media, 2013) Vol. 42.
  * [31] Gao, J., & Efstathiou, K. Self-consistent method and steady states of second-order oscillators. Physical Review E, 98 (4), 042201 (2018).
  * [32] Olmi, S., Navas, A., Boccaletti, S., & Torcini, A. Hysteretic transitions in the Kuramoto model with inertia. Physical Review E, 90 (4), 042905 (2014)
  * [33] Gupta, S., Campa, A.,& Ruffo, S. Nonequilibrium first-order phase transition in coupled oscillator systems with inertia and noise. Physical Review E, 89(2), 022123 (2014).
  * [34] Tanaka, T., & Aoyagi, T. Multistable attractors in a network of phase oscillators with three-body interactions. Physical Review Letters, 106(22), 224101 (2011).
  * [35] Skardal, P. S., & Arenas, A. Abrupt desynchronization and extensive multistability in globally coupled oscillator simplexes. Physical Review Letters, bf 122(24), 248301 (2019)


#  Higgs to diphoton decay rate and the antisymmetric tensor unparticle
mediation

E. O. Iltan
Physics Department, Middle East Technical University
Ankara, Turkey
E-mail address<EMAIL_ADDRESS>
###### Abstract

We study the contribution of the antisymmetric tensor unparticle mediation to
the diphoton production rate of the Higgs boson and try to explain the
discrepancy between the measured value of the decay width of the discovered
new resonance and that of the standard model Higgs boson. We observe that tree
level contribution of the antisymmetric unparticle mediation is a possible
candidate to explain the measured value of the diphoton decay rate.

The standard model (SM) electroweak symmetry breaking mechanism is based on
the existence of a scalar particle, the Higgs boson $H_{0}$, which is crucial
for productions of the masses of fundamental particles. Recently, the ATLAS
and CMS collaborations [1, 2] discovered a resonance with the invariant mass
$125-126\,(GeV)$. At this stage one needs a conformation that the properties
of the discovered resonance coincide with that of the SM Higgs boson. Current
data shows that there is no significant deviation in the decay widths of the
processes $H_{0}\rightarrow W\,W^{*}$ and $H_{0}\rightarrow Z\,Z^{*}$,
however, in the $H_{0}\rightarrow\gamma\gamma$ channel, there is a deviation
from the SM result, namely, the diphoton production rate reaches $1.5$ to $2$
times that of the SM prediction [1, 2, 3, 4]. Even if there needs more data in
order to check whether the excess is based on the statistical fluctuations or
not, a possible attempt to explain this excess from the theoretical side would
be worthwhile and it has been studied in various models beyond the SM
[5]-[38].

In the present work, we consider the antisymmetric tensor unparticle mediation
in order to explain the excess in the diphoton production and we restrict the
free parameters existing in the scenario. Unparticles, being massless, having
non integral scaling dimension $d_{U}$, around the scale $\Lambda_{U}\sim
1.0\,TeV$, are proposed by [39, 40]. They are new degrees of freedom arising
from a hypothetical scale invariant high energy ultraviolet sector with non-
trivial infrared fixed point. In the low energy level the effective
interaction of the SM-unparticle sector reads (see for example [41])

${\cal{L}}_{eff}=\frac{\eta}{\Lambda_{U}^{d_{U}+d_{SM}-n}}\,O_{SM}\,O_{U}\,,$
(1)

where $O_{U}$ ($O_{SM}$) is the unparticle (the SM) operator, $\Lambda_{U}$ is
the energy scale, $n$ is the space-time dimension and $\eta$ is the effective
coefficient [39, 40, 42]. The antisymmetric tensor unparticle propagator which
drives one of the outgoing photon in diphoton production is obtained by the
two point function arising from the scale invariance and it becomes

$\displaystyle\int\,d^{4}x\,e^{ipx}\,<0|T\Big{(}O^{\mu\nu}_{U}(x)\,O^{\alpha\beta}_{U}(0)\Big{)}0>=i\,\frac{A_{d_{U}}}{2\,sin\,(d_{U}\pi)}\,\Pi^{\mu\nu\alpha\beta}(-p^{2}-i\epsilon)^{d_{U}-2}\,,$
(2)

with

$\displaystyle
A_{d_{U}}=\frac{16\,\pi^{5/2}}{(2\,\pi)^{2\,d_{U}}}\,\frac{\Gamma(d_{U}+\frac{1}{2})}{\Gamma(d_{U}-1)\,\Gamma(2\,d_{U})}\,,$
(3)

and the projection operator

$\displaystyle\Pi_{\mu\nu\alpha\beta}=\frac{1}{2}(g_{\mu\alpha}\,g_{\nu\beta}-g_{\nu\alpha}\,g_{\mu\beta})\,.$
(4)

Notice that the projection operator has the transverse and the longitudinal
parts, namely,

$\displaystyle\Pi^{T}_{\mu\nu\alpha\beta}=\frac{1}{2}(P^{T}_{\mu\alpha}\,P^{T}_{\nu\beta}-P^{T}_{\nu\alpha}\,P^{T}_{\mu\beta})\,,\,\,\,\,\,\,\Pi^{L}_{\mu\nu\alpha\beta}=\Pi_{\mu\nu\alpha\beta}-\Pi^{T}_{\mu\nu\alpha\beta}\,,$
(5)

where $P^{T}_{\mu\nu}=g_{\mu\nu}-p_{\mu}\,p_{\nu}/{p^{2}}$ (see [43] and
references therein). At this stage we consider that the scale invariance is
broken at some scale $\mu_{U}$ and we take the antisymmetric tensor unparticle
propagator as

$\displaystyle\int\,d^{4}x\,e^{ipx}\,<0|T\Big{(}O^{\mu\nu}_{U}(x)\,O^{\alpha\beta}_{U}(0)\Big{)}0>=i\,\frac{A_{d_{U}}}{2\,sin\,(d_{U}\pi)}\,\Pi^{\mu\nu\alpha\beta}(-(p^{2}-\mu_{U}^{2})-i\epsilon)^{d_{U}-2}\,,$
(6)

by considering a simple model [44, 45] which provides a rough connection
between the unparticle sector and the particle sector.

Now, we are ready to present the low energy effective Lagrangian which drives
the new contribution to the diphoton production (see [43]):

$\displaystyle{\cal{L}}_{eff}$ $\displaystyle=$
$\displaystyle\frac{g^{\prime}\,\lambda_{B}}{\Lambda_{U}^{d_{U}-2}}\,B_{\mu\nu}\,O^{\mu\nu}_{U}+\frac{g\,\lambda_{W}}{\Lambda_{U}^{d_{U}}}\,(H^{\dagger}\,\tau_{a}\,H)\,W^{a}_{\mu\nu}\,O^{\mu\nu}_{U}\,,$
(7)

where $H$ is the Higgs doublet, $g$ and $g^{\prime}$ are weak couplings,
$\lambda_{B}$ and $\lambda_{W}$ are unparticle-field tensor couplings,
$B_{\mu\nu}$ is the field strength tensor of the $U(1)_{Y}$ gauge boson
$B_{\mu}=c_{W}\,A_{\mu}+s_{W}\,Z_{\mu}$ and $W^{a}_{\mu\nu}$, $a=1,2,3$, are
the field strength tensors of the $SU(2)_{L}$ gauge bosons with
$W^{3}_{\mu}=s_{W}\,A_{\mu}-c_{W}\,Z_{\mu}$ where $A_{\mu}$ and $Z_{\mu}$ are
photon and Z boson fields respectively. The gauge invariant amplitude of the
$H_{0}\rightarrow\gamma\gamma$ decay is

$M=C_{eff}\,(k_{1}.k_{2}\,g^{\mu\nu}-k_{1}^{\nu}\,k_{2}^{\mu})\,\epsilon_{1\mu}\,\epsilon_{2\nu}\,,$
(8)

with the effective coefficient $C_{eff}$ and $i^{th}$ photon polarization
(momentum) four vector $\epsilon_{i\alpha}$ ($k_{i\beta}$). In the framework
of the SM this decay appears at least at the one loop level [46, 47] (see
Appendix for details). On the other hand the antisymmetric tensor unparticle
mediation results in the contribution to the decay in the tree level, with the
transition $H_{0}\rightarrow\gamma O_{U}\rightarrow\gamma\,\gamma$ (see
Fig.1). Here $H_{0}\rightarrow\gamma O_{U}$ transition is carried by the
vertex

$\displaystyle
i\,\frac{e\,v\,\lambda_{W}}{\Lambda_{U}^{d_{U}}}\,k_{1\mu}\,\epsilon_{1\nu}\,O^{\mu\nu}_{U}\,H_{0},$

which arises from the second term in eq.(7). The $O_{U}\rightarrow\gamma$
transition appears with the vertex

$\displaystyle
2\,i\,e\,\Bigg{(}\frac{\lambda_{B}}{\Lambda_{U}^{d_{U}-2}}-\frac{v^{2}\,\lambda_{W}}{4\,\Lambda_{U}^{d_{U}}}\Bigg{)}\,k_{2\mu}\,\epsilon_{2\nu}\,O^{\mu\nu}_{U}\,\,,$

which is coming from the first term and the second term in eq.(7). Here $v$ is
the vacuum expectation value of the SM Higgs $H_{0}$ and $a=3$ is taken in
both vertices. Finally the effective coefficient $C_{eff}$ reads

$\displaystyle C_{eff}=C_{SM}+C_{U}\,,$ (9)

where

$\displaystyle
C_{U}=\frac{-i\,e^{2}\,\lambda_{W}\,v\,\mu_{U}^{2\,(d_{U}-2)}\,A_{d_{U}}}{2\,sin\,(d_{U}\pi)\,\Lambda_{U}^{2\,d_{U}}}\,\Bigg{(}\lambda_{B}\,\Lambda_{U}^{2}-\frac{v^{2}\,\lambda_{W}}{4}\Bigg{)}$
(10)

(see appendix for $C_{SM}$).

Discussion
In this section we study the discrepancy between the measured value of the
decay width of the discovered new resonance, interpreted as the Higgs boson,
and that of the SM one, i.e.,
$\frac{\Gamma(H_{0}\rightarrow\gamma\gamma)_{Measured}}{\Gamma(H_{0}\rightarrow\gamma\gamma)_{SM}}\sim
1.5$. We see that the intermediate antisymmetric tensor unparticle mediation
(see Fig.1) can explain the deviation of diphoton production rate from the SM
prediction. In the present scenario the couplings $\lambda_{B}$,
$\lambda_{W}$, the scale $\Lambda_{U}$, the scaling dimension $d_{U}$ of the
antisymmetric tensor unparticle operator and the scale $\mu_{U}$, which drives
the transition from the unparticle sector to the particle one, are the free
parameters. We take $\lambda_{B}$ and $\lambda_{W}$ as universal and choose
$\lambda_{B}=\lambda_{W}=1$. For the antisymmetric tensor unparticle scale
dimension $d_{U}$ one needs a restriction $d_{U}>2$ not to violate the
unitarity (see [48]). However, since we consider that the scale invariance is
broken at some scale $\mu_{U}$ we relax the restriction and we choose the
range $1<d_{U}<2$ for the scaling dimension $d_{U}$. Here we switched on the
scale invariance breaking by following the simple model [44, 45] which is
based on the redefinition of the unparticle propagator (see eq. (6)). Notice
that the unparticle sector flows to particle sector when $d_{U}$ converges to
one and the range of $d_{U}$ we consider above is appropriate to establish the
connection between these two sectors. For the scale $\mu_{U}$ where the scale
invariance is broken we choose different values $\mu_{U}=1-20\,GeV$ and we
take the scale $\Lambda_{U}$ at the order of the magnitude of $10^{4}\,GeV$.
In our numerical calculations we also consider constraints coming from the
Peskin-Takeuchi parameter, S, which is used to restrict the new physics
contribution to the gauge boson self energy ((see [43]). Here the operator
$(H^{\dagger}\,\tau^{a}\,H)\,W^{a}_{\mu\nu}\,B^{\mu\nu}$ which is induced by
the tensor unparticle exchange results in the S parameter

$\displaystyle S=\frac{A_{d_{U}}}{sin(d_{U}\pi)}\,\frac{g^{2}\,g^{\prime
2}\,\lambda_{B}\,\lambda_{W}\,c_{W}\,v^{2}\,\mu^{2\,(d_{U}-2)}}{s_{W}\,\Lambda_{U}^{2\,d_{U}-2}}\,.$
(11)

The parameter $S$ from the electroweak precision data reads
$S=0.00^{+0.11}_{-0.10}$ [49] which is due to new physics only. Since the
unparticle contribution is negative in our choice of parameters (we choose
$\lambda_{B}=\lambda_{W}=1$), we take the lower bound $S_{LB}=-0.10$ and we
plot the contour diagram of the S parameter with $1\,\sigma$ bound in
$\Lambda_{U},\mu_{U}$ plane for different values of $d_{U}$, as shown in
figure 2. In this figure the pair of parameters $\Lambda_{U},\mu_{U}$ under
the curves are excluded.

Now, we start to study the discrepancy between the measured value of the decay
width of the discovered new resonance and the SM Higgs boson by considering
the restriction region coming from the S parameter. Fig.3 represents $d_{U}$
dependence of the ratio
$r=\frac{\Gamma(H_{0}\rightarrow\gamma\gamma)_{SM+U}}{\Gamma(H_{0}\rightarrow\gamma\gamma)_{SM}}$
for different values of the scales $\Lambda_{U}$ and $\mu_{U}$. Here111The
solid straight line represents the ratio
$\frac{\Gamma(H_{0}\rightarrow\gamma\gamma)_{Measured}}{\Gamma(H_{0}\rightarrow\gamma\gamma)_{SM}}\sim
1.5$ in each figure. the upper-intermediate-lower solid (dashed) line
represents $r$ for $\mu_{U}=1-10-20\ (GeV)$,
$\Lambda_{U}=5000\,(10000)\,(GeV)$. The ratio is strongly sensitive to the
scale dimension $d_{U}$ for the values far from $1.9$ and the decrease in the
scale $d_{U}$ results in the increase in the unparticle contribution which
makes it possible to overcome the the discrepancy between the measured value
and the SM result. We observe that the measured value is reached if the
scaling dimension is in the range $1.48<d_{U}<1.68$ for the given numerical
values of $\Lambda_{U}$ and $\mu_{U}$. In the case of
$\Lambda_{U}=10000\,(5000)\,(GeV)$ and $\mu_{U}=1\,(GeV)$, the measured decay
rate is obtained for $d_{U}\sim 1.63\,(1.68)$. For $\mu_{U}=20\,(GeV)$ the
measured value is reached for $d_{U}\sim(1,49)\,1.54$.

Fig.4 is devoted to $\mu_{U}$ dependence of the ratio $r$ for
$\Lambda_{U}=10000\,(GeV)$ and different values $d_{U}$. Here the solid (long
dashed, dashed, dotted) line represents $r$ for
$d_{U}=1.4\,(1.5,\,1.6,\,1.7)$. The ratio is sensitive to the scale $\mu_{U}$
and increases with its decreasing value. One can reach the measured decay rate
for $2.3\,(GeV)<\mu_{U}<16\,(GeV)$ if $d_{U}$ is in the range $d_{U}\sim
1.50-1.60$.

In Fig.5, we present $\Lambda_{U}$ dependence of the ratio $r$ for different
values of $d_{U}$ and $\mu_{U}$. Here the solid (long dashed, dashed, dotted)
line represents $r$ for $d_{U}=1.5;\,\mu_{U}=1.0\,(GeV)$
$(d_{U}=1.5;\,\mu_{U}=10\,(GeV),\,d_{U}=1.6;\,\mu_{U}=1.0\,(GeV),\,d_{U}=1.6;\,\mu_{U}=10\,(GeV))$.
The measured decay rate is reached for $d_{U}=1.6;\,\mu_{U}=1.0\,(GeV)$ if the
energy scale reads $\Lambda_{U}\sim 17000\,(GeV)$.

As a summary, we show that the intermediate antisymmetric tensor unparticle
mediation is a possible candidate to overcome the deviation of diphoton
production rate from the SM prediction. We study the ratio
$r=\frac{\Gamma(H_{0}\rightarrow\gamma\gamma)_{SM+U}}{\Gamma(H_{0}\rightarrow\gamma\gamma)_{SM}}$
and see that $r\sim 1.5$ is reached if the scaling dimension is almost in the
range $1.48<d_{U}<1.68$ for the given numerical values of
$5000\,(GeV)<\Lambda_{U}<17000\,(GeV)$ and $1.0\,(GeV)<\mu_{U}<20\,(GeV)$.
This result makes it possible to explain the discrepancy between the measured
value of the decay width of the discovered new resonance and that of the SM
Higgs boson. In addition, it also gives an opportunity to understand the role
and the type of the unparticle scenario and to determine the existing free
parameters.

Appendix

In the framework of the SM, the $H_{0}\rightarrow\gamma\gamma$ decay appears
at least in the loop level with the internal W boson and fermions where the
top quark gives the main contribution. The gauge invariant amplitude reads

$\displaystyle
M=C_{SM}\,(k_{1}.k_{2}\,g^{\mu\nu}-k_{1}^{\nu}\,k_{2}^{\mu})\,\epsilon_{1\mu}\,\epsilon_{2\nu}\,,$
(12)

where $C_{SM}=\frac{\alpha_{EM}\,g}{4\,\pi\,m_{W}}\,F(x_{W},x_{f})$ and

$\displaystyle
F(x_{W},x_{f})=F_{1}(x_{W})+\sum_{f}N_{C}\,Q_{f}^{2}\,F_{2}(x_{f})\,,$ (13)

with

$\displaystyle F_{1}(x)=2+3\,x+3\,x\,(2-x)\,g(x)\,,$ $\displaystyle
F_{2}(x)=-2\,\Big{(}1+(1-x)\,g(x)\Big{)}\,.$ (14)

Here

$\displaystyle
g(x)=\left\\{\begin{array}[]{ll}arcsin^{2}(x^{-1/2}),&\hbox{$x\geq 1$;}\\\
\frac{-1}{4}\Bigg{(}ln\frac{1+\sqrt{1-x}}{1-\sqrt{1-x}}-i\pi\Bigg{)}^{2},&\hbox{$x<1$
,}\\\ \end{array}\right.$ (17)

and $Q_{f}$ is the charge of the fermion $f$, $N_{C}=1\,(3)$ for lepton
(quark), $x_{i}=\frac{4\,m_{i}^{2}}{m_{H_{0}}^{2}}$, $i=W,f$. Finally the
decay width $\Gamma(H_{0}\rightarrow\gamma\gamma)$ is obtained as

$\displaystyle\Gamma(H_{0}\rightarrow\gamma\gamma)=\frac{m^{3}_{H_{0}}\,|C_{SM}|^{2}}{64\,\pi}\,.$
(18)

## References

  * [1] ATLAS Collaboration, G. Aad et al., Phys. Lett. B716, 1 (2012).
  * [2] CMS Collaboration, S. Chatrchyan et al., Phys. Lett. B716, 30 (2012).
  * [3] Tevatron New Physics Higgs Working Group and CDF and D0 Collaborations, arXiv:1207.0449 [hep-ex].
  * [4] G. Aad et al. [ATLAS Collaboration], ATLAS-CONF-2012-091; S. Chatrchyan et al. [CMS Collaboration], CMS-PAS-HIG-12-015.
  * [5] V. Ilisie and A. Pich, Phys. Rev. D86, 033001 (2012).
  * [6] V. Barger, M. Ishida, W. Y. Keung, Phys. Rev. Lett. 108, 261801 (2012).
  * [7] M. Carena, S. Casagrande, F. Goertz, U. Haisch, M. Neubert, JHEP 1208, 156 (2012).
  * [8] P. Draper, D. McKeen, Phys. Rev. D85, 115023 (2012).
  * [9] K. Kumar, R. Vega-Morales, F. Yu, arXiv:1205.4244.
  * [10] S. Dawson, E. Furlan, Phys. Rev. D86, 015021 (2012).
  * [11] A. G. Akeroyd, S. Moretti, Phys. Rev. D86, 035015 (2012).
  * [12] L. Wang, X. -F. Han, Phys. Rev. D86, 095007 (2012).
  * [13] M. Carena, I. Low, C. E. Wagner, JHEP 1208, 060 (2012).
  * [14] W. F. Chang, J. N. Ng, J. M. S. Wu, Phys. Rev. D86, 033003 (2012).
  * [15] C. W. Chiang, K. Yagyu, arXiv:1207.1065.
  * [16] H. Cheon, S. K. Kang, arXiv:1207.1083.
  * [17] P. P. Giardinoa, K. Kannike, M.Raidal, A. Strumi, arXiv:1207.1347.
  * [18] M. R. Buckley, D. Hooper, Phys. Rev. D86, 075008 (2012).
  * [19] H. An, T. Liu, L. T. Wang, Phys. Rev. D86, 075030 (2012).
  * [20] A. Joglekar, P. Schwaller, C. E. Wagner, arXiv:1207.4235.
  * [21] N. A. Hamed, K. Blum, R. T. D. Agnolo, J. Fan, arXiv:1207.4482.
  * [22] N. Haba, K. Kaneta, Y. Mimura, R. Takahashi, arXiv:1207.5102.
  * [23] L. G. Almeida, E. Bertuzzo, P. A. Machado, R. Z. Funchal, JHEP 1211, 085 (2012).
  * [24] A. Delgado, G. Nardini, M. Quiros, arXiv:1207.6596.
  * [25] J. Kearney, A. Pierce, N. Weiner, arXiv:1207.7062.
  * [26] J. R. Espinosa, C. Grojean, V. Sanz, M. Trott, arXiv:1207.7355.
  * [27] M. Carena, S. Gori, N. R. Shah, C. E. Wagner, L. T. Wang, JHEP 1207, 175 (2012).
  * [28] M. A. Ajaib, I. Gogoladze, Q. Shafi, arXiv:1207.7068.
  * [29] N. D. Christensen, T. Han, S. Su, Phys. Rev. D85, 115018 (2012).
  * [30] K. Hagiwara, J. S. Lee, J. Nakamura, JHEP 1210, 002 (2012).
  * [31] R. Benbrik, M. G. Bock, S. Heinemeyer, O. Stal, G. Weiglein, L. Zeune, Eur. Phys. J. C72, 2171 (2012).
  * [32] I. Dorsner, S. Fajfer, A. Greljo, J. F. Kamenik, arXiv:1208.1266.
  * [33] K. S. Hoberg, F. Staub, JHEP 1210, 195 (2012).
  * [34] H. Davoudiasl, H.-S. Lee, W. J. Marciano, Phys. Rev. D86, 095009 (2012).
  * [35] M. Carena, S. Gori, N. R. Shah, C. E. Wagner, JHEP 1203, 014 (2012).
  * [36] J. J. Cao, Z.-X. Heng, J. M. Yang, Y. M. Zhang, J. Y. Zhu, JHEP 1203, 086 (2012).
  * [37] T. Kitahara, JHEP 1211, 21 (2012).
  * [38] M. Chala, arXiv:1210.6208.
  * [39] H. Georgi, Phys. Rev. Lett. 98, 221601 (2007).
  * [40] H. Georgi, Phys. Lett. B650, 275 (2007).
  * [41] S. L. Chen, X. G. He, Phys. Rev. D76, 091702 (2007).
  * [42] R. Zwicky, Phys. Rev. D77, 036004 (2008).
  * [43] T. Hur, P. Ko, X. H. Wu, Phys. Rev. D76, 096008 (2007).
  * [44] P. J. Fox, A. Rajaraman, Y. Shirman, Phys. Rev. D76, 075004 (2007).
  * [45] A. Rajaraman, Phys. Lett. B671, 411 (2009).
  * [46] J. R. Ellis, M. K. Gaillard and D. V. Nanopoulos, Nucl. Phys. B106, 292 (1976).
  * [47] M. A. Shifman, A. I. Vainshtein, M. B. Voloshin and V. I. Zakharov, Sov. J. Nucl. Phys. 30, 711 (1979) [Yad. Fiz. 30, 1368 (1979)].
  * [48] B. Grinstein, K. A. Intriligator, I. Z. Rothstein, Phys. Lett. B662, 367 (2008).
  * [49] J. Beringer et al. (Particle Data Group), Phys. Rev. D86, 010001 (2012).

Figure 1: Tree level diagram contributing to diphoton decay due to the
antisymmetric tensor unparticle mediation. Solid (wavy, double dashed) line
represents the Higgs (electromagnetic, antisymmetric tensor unparticle) field.
Figure 2:  $\Lambda_{U}$ with respect to $\mu_{U}$ within $1\,\sigma$ bound.
Here the solid (long dashed, dashed, dotted) line represents $\Lambda_{U}$ for
$d_{U}=1.4\,(1.5,1.6,1.7)$. Figure 3:  $r$ with respect to $d_{U}$. Here the
upper-intermediate-lower solid (dashed) line represents $r$ for
$\mu_{U}=1-10-20\ (GeV)$, $\Lambda_{U}=5000\,(10000)\,(GeV)$. Figure 4: $r$
with respect to $\mu_{U}$ for $\Lambda_{U}=10000\,(GeV)$. Here the solid (long
dashed, dashed, dotted) line represents $r$ for
$d_{U}=1.4\,(1.5,\,1.6,\,1.7)$.  Figure 5: $r$ with respect to $\Lambda_{U}$.
Here the solid (long dashed, dashed, dotted) line represents $r$ for
$d_{U}=1.5;\,\mu_{U}=1.0\,(GeV)$
$(d_{U}=1.5;\,\mu_{U}=10\,(GeV),\,d_{U}=1.6;\,\mu_{U}=1.0\,(GeV),\,d_{U}=1.6;\,\mu_{U}=10\,(GeV))$.