# Energy momentum tensor on and off the light cone

Xianghui Cao  Department of Modern Physics, University of Science and
Technology of China, Hefei 230026, China  Siqi Xu  Institute of Modern
Physics, Chinese Academy of Sciences, Lanzhou 730000, China  School of Nuclear
Science and Technology, University of Chinese Academy of Sciences, Beijing
100049, China  Yang Li  Department of Modern Physics, University of Science
and Technology of China, Hefei 230026, China  Department of Physics and
Astronomy, Iowa State University, Ames, Iowa 50011  Guangyao Chen  Department
of Physics, Jacksonville University, Jacksonville, FL, 32211  Xingbo Zhao
Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou 730000,
China  School of Nuclear Science and Technology, University of Chinese Academy
of Sciences, Beijing 100049, China  CAS Key Laboratory of High Precision
Nuclear Spectroscopy, Institute of Modern Physics, Chinese Academy of
Sciences, Lanzhou 730000, China  V. A. Karmanov  Lebedev Physical Institute,
Leninsky Prospekt 53, 119991 Moscow, Russia  J. P. Vary  Department of Physics
and Astronomy, Iowa State University, Ames, Iowa 50011

###### Abstract

We compute the gravitational form factors $A_{i}$, $D_{i}$ and $\bar{c}_{i}$
of a scalar theory using both the light-cone and covariant perturbation theory
at the one-loop level. The light-cone formalism provides a potential approach
to access these form factors beyond the perturbative regime. However, unlike
the covariant formulation, the Poincaré symmetry on the light cone is not
manifest. In this work, we use perturbation theory as a benchmark to extract
the gravitational form factors from the light-front energy-momentum tensor. By
comparing results on and off the light cone, we identify
$T^{++},T^{+a},T^{+-},T^{12}$ as the “good currents” that are properly
renormalized and can be used to extract the gravitational form factors.

##  I Introduction

The energy-momentum tensor (EMT) $T^{\mu\nu}$ is an observable of fundamental
significance. It encodes the internal energy, momentum, angular momentum and
stress distributions of a system [1]. Therefore, this quantity is particularly
interesting for unraveling the structures of hadrons, quarks and gluons bound
together by the strong force [2, 3]. The hadronic matrix element (HME) of the
EMT can be parametrized in terms of Lorentz tensors [4, 5]. The scalar
functions that are specific to each hadron are known as the gravitational form
factors (GFFs). These quantities then probe the internal dynamics of the
hadron [6, 7].

Since QCD is strongly-coupled at the hadron scale, non-perturbative methods
are needed to access the GFFs. Among various approaches, the light-front
Hamiltonian formulation is poised as a natural framework to describe hadrons
as relativistic self-bound states [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18].
This approach is demonstrated in one of our recent works for a strongly-
coupled scalar Yukawa theory [19]. In that work, we extracted GFFs $A(Q^{2})$
and $D(Q^{2})$ from $T^{++}$ and $T^{+-}$, two tensor components of the EMT.
We showed that these two components can be renormalized using the same
counterterms obtained from solving the light-front Schrödinger equation, i.e.
the renormalization of the light-front Hamiltonian [20, 21, 22, 23].

Further quantities of interest are the GFFs associated with each constituent,
which allow us to decompose a hadronic density into contributions from each
constituent, e.g. quarks and gluons [24, 25, 26, 27]. For the scalar theory,
apart from GFFs $A_{i}(Q^{2})$ and $D_{i}(Q^{2})$, where $i$ enumerates over
the constituents, there is an additional GFF, $\bar{c}_{i}(Q^{2})$ [2]. It
represents the force between constituent $i$ and the rest of the system.
Current conservation $\partial_{\mu}T^{\mu\nu}=0$ implies that the total
$\bar{c}$ vanishes,

$\bar{c}(Q^{2})\equiv\sum_{i}\bar{c}_{i}(Q^{2})=0\,.$ (1)

Namely, all inter-particle forces balance out for an isolated system – the
quantum version of Newton’s third law. To extract $\bar{c}_{i}(Q^{2})$, using
the above two currents $T_{i}^{++}$ and $T^{+-}_{i}$ are not sufficient, and
we need additional currents such as $T^{12}_{i},T^{11}_{i},T^{22}_{i}$.
However, these operators contain uncanceled divergences and lead to a
violation of the current conservation. These difficulties are caused by the
breaking of the Poincaré symmetry due to truncation within the light-front
formalism. This effect is better encapsulated in the explicitly covariant
formulation of light-front dynamics (CLFD) (see Ref. [12] for a review). In
CLFD, $T^{12},T^{11},T^{22}$ are contaminated by the spurious contributions
that depend on the orientation of the light front [19]. The most general
Lorentz structures of the EMT taking into account CLFD are given by Lorcé
[28].

Facing these difficulties in the non-perturbative calculation, we turn to
perturbation theory for guidance. Recently, More et. al. evaluated the EMT of
a quark dressed by gluons using leading order (LO) light-cone perturbation
theory (LCPT) [29, 30]. Their results can be compared with the recent LO QED
calculations in covariant perturbation theory (CPT) [31, 32, 33, 34]. More et.
al. suggested to use spin-flip HMEs to extract GFFs $D_{i}(Q^{2})$ and
$\bar{c}_{i}(Q^{2})$, which are unfortunately not available in the scalar
case. Furthermore, they obtained a non-vanishing $\bar{c}(Q^{2})$ except at
$Q=0$. They attributed the violation of current conservation to zero modes
[30]. In this work, we compute the EMT of the scalar Yukawa theory using both
covariant and LCPT at the one-loop level. Scalar theory is among the simplest
yet non-trivial quantum field theories in 3+1 dimensions [35, 36, 37, 38, 39,
40, 41, 42, 43, 44, 45, 46, 19, 47]. The zero modes within the scalar theory
can be removed using Pauli-Villars regularization [48, 16]. Therefore, it is
an ideal testbed for developing computational techniques in light-front field
theory. By comparing the covariant and light-cone perturbation theory for the
scalar Yukawa model, we identify the sources of the divergences and explore
how to properly extract GFFs $A_{i}(Q^{2})$, $D_{i}(Q^{2})$ and
$\bar{c}_{i}(Q^{2})$ from the light-front EMT.

The remaining part of the article is organized as follows. We introduce the
scalar Yukawa model including the one-loop renormalization of this theory in
Sect. II. We then compute the energy-momentum tensor up to one-loop in both
the covariant perturbation theory and in light-cone perturbation theory in
Sect. III. We analyze these results using covariant light-front dynamics in
Sect. IV. Finally, we summarize in Sect. V.

##  II Scalar Yukawa theory

In this work, we consider the scalar Yukawa theory. The classical Lagrangian
density of this model reads [20, 21, 22, 23],

$\mathscr{L}=\partial_{\mu}N^{\dagger}\partial^{\mu}N-m^{2}N^{\dagger}N\\\
+\frac{1}{2}\partial_{\mu}\pi\partial^{\mu}\pi-\frac{1}{2}\mu^{2}\pi^{2}+gN^{\dagger}N\pi\,.$
(2)

This Lagrangian describes a complex scalar field $N$ interacting with a real
scalar field $\pi$ via a Yukawa type coupling. It can be used as a rudimentary
model for the nucleon-pion interaction. Therefore, we will refer to the
complex scalar field $N$ as the (mock) nucleon field and tentatively assign
the mass of the nucleon $m=0.94\,\mathrm{GeV}$ to its physical mass.
Similarly, the real scalar field $\pi$ will be referred to as the (mock) pion
field, and tentatively assigned the pion mass $\mu=0.14\,\mathrm{GeV}$. The
coupling $g$ is dimensionful. Hence, this model is super-renormalizable.
However, there are still divergences within the one-loop radiative correction.
It is useful to introduce a dimensionless coupling $\alpha=g^{2}/(16\pi
m^{2})$. This coupling is the Yukawa coupling appearing in the tree-level
Yukawa potential,

$V(r)=-\frac{\alpha}{r}e^{-\mu r}\,.$ (3)

Specifically, we consider the nucleon dressed by the pion up to order
$O(\alpha)$. We adopt the bare renormalization scheme [49], i.e. the field
operator is not normalized, $\langle
0|N^{\dagger}(x)|p\rangle=\sqrt{Z}e^{ip\cdot x}$. Instead, the state vector
will be normalized to unity, $\langle
p^{\prime}|p\rangle=2p^{+}(2\pi)^{3}\delta^{3}(p-p^{\prime})$. The
normalization of the state vector is easier to be generalized to the non-
perturbative regime in the light-front Hamiltonian formulation [23]. Within
this scheme, the Lagrangian can be written as,

$\mathscr{L}=\partial_{\mu}N^{\dagger}\partial^{\mu}N-m^{2}N^{\dagger}N\\\
+\frac{1}{2}\partial_{\mu}\pi\partial^{\mu}\pi-\frac{1}{2}\mu^{2}\pi^{2}+g_{0}N^{\dagger}N\pi+\delta
m^{2}N^{\dagger}N+\frac{1}{2}\delta\mu^{2}\pi^{2}\,,$ (4)

where, $\delta m^{2}$ and $\delta\mu^{2}$ are mass counterterms and $g_{0}$ is
the bare coupling. These parameters can be obtained from the mass and coupling
constant renormalization.

Figure 1: Leading order diagram for the self-energy of the scalar Yukawa model

We adopt two types of regularization schemes, the dimensional regularization
(DR) and Pauli-Villars regularization (PV). The former is the standard choice
in covariant perturbation theory and the latter is useful for both the
covariant and light-cone perturbation theory [50, 51, 52, 53, 54, 48]. Then
the nucleon self-energy up to the LO order (see Fig. 1) is,

$\Sigma(p^{2})=-\frac{g_{0}^{2}}{16\pi^{2}}\begin{dcases*}\frac{1}{\epsilon}+\ln\frac{4\pi}{e^{\gamma_{E}}}-\int_{0}^{1}{\mathrm{d}}x\ln\big{[}(1-x)\mu^{2}+xm^{2}-x(1-x)p^{2}\big{]},&\text{DR}\\\
\int_{0}^{1}{\mathrm{d}}x\ln\frac{(1-x)\mu^{2}_{\textsc{pv}}+xm^{2}-x(1-x)p^{2}}{(1-x)\mu^{2}+xm^{2}-x(1-x)p^{2}},&\text{PV}\end{dcases*}$
(5)

where, $\epsilon=2-d/2$, and $d$ is the spacetime dimension.
$\gamma_{E}\approx 0.577216$ is the Euler’s constant. With the on-shell
scheme, the mass counterterm is then $\delta m^{2}=\Sigma(m^{2})$. The
$Z$-factor is,

$Z=\Big{[}1-\frac{{\mathrm{d}}\Sigma(p^{2})}{{\mathrm{d}}p^{2}}\Big{]}^{-1}\Big{|}_{p^{2}\to
m^{2}}.$ (6)

Note that $Z$ is not divergent in the scalar Yukawa model. The bare coupling
is obtained from imposing the physical coupling $g$ to the nucleon-pion
scattering amplitude at the on-shell point. And the result, according to the
Lehmann-Symanzik-Zimmermann reduction formula [49], is,

$g_{0}=\frac{g}{Z\sqrt{Z^{\prime}}}\,,$ (7)

where $Z^{\prime}$ is the field renormalization constant of the pion,

$Z^{\prime}=\Big{[}1+\frac{g_{0}^{2}}{(4\pi)^{2}}\int_{0}^{1}dx\frac{x(1-x)}{m^{2}-x(1-x)\mu^{2}}\Big{]}^{-1}\,.$
(8)

##  III Energy-momentum tensor

The EMT of the scalar Yukawa theory is,

$T^{\mu\nu}=\partial^{\\{\mu}N^{\dagger}\partial^{\nu\\}}N-g^{\mu\nu}(\partial_{\rho}N^{\dagger}\partial^{\rho}N-m^{2}N^{\dagger}N)-g^{\mu\nu}g_{0}N^{\dagger}N\pi-g^{\mu\nu}\delta
m^{2}N^{\dagger}N\\\
+\partial^{\mu}\pi\partial^{\nu}\pi-\frac{1}{2}g^{\mu\nu}(\partial_{\sigma}\pi\partial^{\sigma}\pi-\mu^{2}\pi^{2})-\frac{1}{2}g^{\mu\nu}\delta\mu^{2}\pi^{2}\,.$
(9)

Here, $a^{\\{\mu}b^{\nu\\}}\equiv a^{\mu}b^{\nu}+a^{\nu}b^{\mu}$. This EMT is
symmetric, and there is no distinction between the canonical EMT and the
Belinfante–Rosenfeld-Hilbert-Einstein EMT. In analogy to mass decomposition in
QCD [24, 25], we split the EMT into two parts,

$\displaystyle T^{\mu\nu}_{N}=\,$
$\displaystyle\partial^{\\{\mu}N^{\dagger}\partial^{\nu\\}}N-g^{\mu\nu}(\partial_{\rho}N^{\dagger}\partial^{\rho}N-m^{2}N^{\dagger}N)-g^{\mu\nu}gN^{\dagger}N\pi-g^{\mu\nu}\delta
m^{2}N^{\dagger}N\,,$ $\displaystyle T^{\mu\nu}_{\pi}=\,$
$\displaystyle\partial^{\mu}\pi\partial^{\nu}\pi-\frac{1}{2}g^{\mu\nu}(\partial_{\sigma}\pi\partial^{\sigma}\pi-\mu^{2}\pi^{2})-\frac{1}{2}g^{\mu\nu}\delta\mu^{2}\pi^{2}\,,$
(10)

i.e. the nucleon EMT contains the interaction. The HMEs of the EMT for a
spin-0 particle can be parametrized by the following Lorentz tensors [2],

$t^{\alpha\beta}_{i}=\langle
p^{\prime}|T^{\alpha\beta}_{i}(0)|p\rangle=2P^{\alpha}P^{\beta}A_{i}(-q^{2})+\frac{1}{2}\big{(}q^{\alpha}q^{\beta}-g^{\alpha\beta}q^{2}\big{)}D_{i}(-q^{2})+2M^{2}g^{\alpha\beta}\bar{c}_{i}(-q^{2})\,.$
(11)

Here, $P=(p^{\prime}+p)/2$, $q=p^{\prime}-p$. $M$ is the mass of the scalar,
i.e. $p^{2}=p^{\prime 2}=M^{2}$. The Lorentz scalar functions $A_{i}(-q^{2})$,
$D_{i}(-q^{2})$ and $\bar{c}_{i}(-q^{2})$ are the GFFs that encode the
internal dynamics of the physical particle. The subscript $i=N,\pi$ enumerates
over the constituents. Summing over all constituents gives the total EMT.
Accordingly,

$\displaystyle\sum_{i}A_{i}(-q^{2})=\,$ $\displaystyle A(-q^{2}),$ (12)
$\displaystyle\sum_{i}D_{i}(-q^{2})=\,$ $\displaystyle D(-q^{2}),$ (13)
$\displaystyle\sum_{i}\bar{c}_{i}(-q^{2})=\,$
$\displaystyle\bar{c}(-q^{2})\,.$ (14)

Current conservation

$\partial_{\nu}T^{\mu\nu}=0\,,$ (15)

implies

$\bar{c}(-q^{2})=0\,.$ (16)

Physically, $\bar{c}_{i}$ represents the force between the $i$-th constituent
and the rest of the system. Eq. (16) is nothing but Newton’s third law applied
to hadrons: the total force vanishes for a self-bound system. Current
conservation also protects the EMT $T^{\alpha\beta}$ against the quantum
fluctuation. Consequently, the GFFs $A(-q^{2})$ and $D(-q^{2})$ do not depend
on the renormalization scheme or scale. By contrast, the GFFs $A_{i}(-q^{2})$,
$D_{i}(-q^{2})$ and $\bar{c}_{i}(-q^{2})$ for each constituent are
renormalization scheme and scale dependent. Unless otherwise stated, these
dependences are suppressed for simplicity. Despite this arbitrariness, the
decomposition of the EMT provides insights into the system. For example, in
QCD, mass decomposition based on the EMT provides answers to the origin of the
proton mass in terms of quark mass, quark kinetic energy, gluon kinetic energy
and the anomalous mass [24, 25, 26, 27].

Figure 2: Leading-order diagrams for the energy-momentum tensor. The solid
lines with an arrow represent the nucleon. The wavy lines represent the pion.
The cross represents the mass counterterm. The circled crosses represent the
operator insertion.

Figure 2 shows the LO diagrams for the EMT. Diagrams ($a$)–($e$) contribute to
$T^{\mu\nu}_{N}$ while diagram ($f$) is the contribution to
$T^{\mu\nu}_{\pi}$. We first evaluate these diagrams in covariant perturbation
theory and extract the relevant GFFs. These results will serve as a benchmark
for the LCPT.

###  III.1 Covariant perturbation theory

The covariant calculation is straightforward. The results are,

$\displaystyle t^{\mu\nu}_{a}=\,$ $\displaystyle
Z\Big{[}2P^{\mu}P^{\nu}-\frac{1}{2}(q^{\mu}q^{\nu}-q^{2}g^{\mu\nu})\Big{]},$
(17) $\displaystyle t^{\mu\nu}_{b}=\,$
$\displaystyle-t^{\mu\nu}_{c}=-t^{\mu\nu}_{d}=-Zg^{\mu\nu}\Sigma(m^{2}),$ (18)
$\displaystyle t^{\mu\nu}_{e}=\,$ $\displaystyle
Zg_{0}^{2}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\int\frac{{\mathrm{d}}^{4}\ell_{E}}{(2\pi)^{4}}\frac{\mathcal{O}^{\mu\nu}_{e}}{(\ell_{E}^{2}+\Delta_{e})^{3}},$
(19) $\displaystyle t^{\mu\nu}_{f}=\,$ $\displaystyle
Zg_{0}^{2}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\int\frac{{\mathrm{d}}^{4}\ell_{E}}{(2\pi)^{4}}\frac{\mathcal{O}^{\mu\nu}_{f}}{(\ell_{E}^{2}+\Delta_{f})^{3}}\,,$
(20)

where, $\Delta_{e}=(x+y)^{2}m^{2}+(1-x-y)\mu^{2}-xyq^{2}$,
$\Delta_{f}=(1-x-y)^{2}m^{2}+(x+y)\mu^{2}-xyq^{2}$ and
$\mathcal{O}^{\mu\nu}_{e}=g^{\mu\nu}\ell^{2}_{E}+2(1-x-y)^{2}2P^{\mu}P^{\nu}+((x-y)^{2}-1)(q^{\mu}q^{\nu}-q^{2}g^{\mu\nu})+2g^{\mu\nu}\Big{\\{}\big{[}1-(1-x-y)^{2}\big{]}m^{2}+\frac{1}{2}(x^{2}-x+y^{2}-y)q^{2}\Big{\\}}$,
$\mathcal{O}^{\mu\nu}_{f}=\ell^{\mu\nu}_{E}+2(1-x-y)^{2}2P^{\mu}P^{\nu}+\big{[}(x-y)^{2}-1\big{]}(q^{\mu}q^{\nu}-g^{\mu\nu}q^{2})+2g^{\mu\nu}\big{\\{}\mu^{2}-(1-x-y)^{2}m^{2}+\frac{1}{2}(x^{2}-x+y^{2}-y)q^{2}\big{\\}}$.

From these HMEs, we can extract the GFFs. It turned out, GFFs $A_{i}$ and
$D_{i}$ are divergence free,

$\displaystyle A_{N}(Q^{2})=\,$ $\displaystyle
Z+Z\frac{g_{0}^{2}}{(4\pi)^{2}}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\frac{(1-x-y)^{2}}{\Delta_{e}},$
(21) $\displaystyle A_{\pi}(Q^{2})=\,$ $\displaystyle
Z\frac{g_{0}^{2}}{(4\pi)^{2}}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\frac{(1-x-y)^{2}}{\Delta_{f}},$
(22) $\displaystyle D_{N}(Q^{2})=\,$ $\displaystyle-
Z-Z\frac{g_{0}^{2}}{(4\pi)^{2}}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\frac{1-(x-y)^{2}}{\Delta_{e}},$
(23) $\displaystyle D_{\pi}(Q^{2})=\,$
$\displaystyle-Z\frac{g_{0}^{2}}{(4\pi)^{2}}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\frac{1-(x-y)^{2}}{\Delta_{f}}\,.$
(24)

On the other hand, GFFs $\bar{c}_{i}$ are UV divergent. The results with
dimensional regularization are,

$\bar{c}_{N}^{(\overline{\textsc{ms}})}(Q^{2},\mu^{2}_{R})=-\frac{Zg_{0}^{2}}{64\pi^{2}m^{2}}\Big{(}\frac{1}{\epsilon}+1\Big{)}+\frac{Zg_{0}^{2}}{32\pi^{2}m^{2}}\int_{0}^{1}{\mathrm{d}}x\ln\frac{(1-x)\mu^{2}+x^{2}m^{2}}{\mu^{2}_{R}}-\frac{Zg_{0}^{2}}{32\pi^{2}m^{2}}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\ln\frac{\Delta_{e}}{\mu^{2}_{R}}\\\
+\frac{Zg_{0}^{2}}{64\pi^{2}m^{2}}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\frac{2(x+y)(2-x-y)m^{2}+(x^{2}-x+y^{2}-y)q^{2}}{\Delta_{e}},$
(25)

and,

$\bar{c}_{\pi}^{(\overline{\textsc{ms}})}(Q^{2},\mu_{R}^{2})=\frac{Zg_{0}^{2}}{64\pi^{2}m^{2}}\Big{(}\frac{1}{\epsilon}-1\Big{)}-\frac{Zg_{0}^{2}}{32\pi^{2}m^{2}}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\ln\frac{\Delta_{f}}{{\mu}^{2}_{R}}\\\
+\frac{Zg_{0}^{2}}{64\pi^{2}m^{2}}\int_{0}^{1}{\mathrm{d}}x\int\limits_{0}^{1-x}{\mathrm{d}}y\frac{2\mu^{2}-2(1-x-y)^{2}m^{2}+(x^{2}-x+y^{2}-y)q^{2}}{\Delta_{f}}.$
(26)

Here, $\mu_{R}$ is the renormalization scale in the $\overline{\mathrm{MS}}$
scheme. The GFFs $\bar{c}_{i}$ can also be obtained with the Pauli-Villars
regularization. The results are,

$\bar{c}_{N}^{(\textsc{pv})}(Q^{2},\mu_{\textsc{pv}}^{2})=-\frac{Zg_{0}^{2}}{32\pi^{2}m^{2}}\int_{0}^{1}{\mathrm{d}}x\ln\frac{(1-x)\mu^{2}_{\textsc{pv}}+x^{2}m^{2}}{(1-x)\mu^{2}+x^{2}m^{2}}+\frac{Zg_{0}^{2}}{32\pi^{2}m^{2}}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\ln\frac{\Delta_{e}^{(1)}}{\Delta_{e}}\\\
+\frac{Zg_{0}^{2}}{64\pi^{2}m^{2}}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\sum_{j}(-1)^{j}\frac{2(x+y)(2-x-y)m^{2}+(x^{2}-x+y^{2}-y)q^{2}}{\Delta_{e}^{(j)}},$
(27)

Here, the index “$j$” runs over 0 and 1 for the original field and the Pauli-
Villars field respectively. And,

$\bar{c}_{\pi}^{(\textsc{pv})}(Q^{2},\mu_{\textsc{pv}}^{2})=\frac{Zg_{0}^{2}}{32\pi^{2}m^{2}}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\ln\frac{\Delta_{f}^{(1)}}{\Delta_{f}}\\\
+\frac{Zg_{0}^{2}}{64\pi^{2}m^{2}}\int_{0}^{1}{\mathrm{d}}x\int_{0}^{1-x}{\mathrm{d}}y\sum_{j}(-1)^{j}\frac{2\mu_{j}^{2}-2(1-x-y)^{2}m^{2}+(x^{2}-x+y^{2}-y)q^{2}}{\Delta_{f}^{(j)}}.$
(28)

In these expressions, $\Delta_{e}^{(0)}=\Delta_{e}$, and $\Delta_{e}^{(1)}$ is
obtained from $\Delta_{e}$ by replacing $\mu$ with
$\mu_{1}=\mu_{\textsc{pv}}$. Similarly, $\Delta_{f}^{(0)}=\Delta_{f}$, and
$\Delta_{f}^{(1)}$ is obtained from $\Delta_{f}$ by replacing $\mu$ with
$\mu_{1}=\mu_{\textsc{pv}}$.

The obtained $\bar{c}_{i}(Q^{2})$ for both the nucleon and the pion is
constant, i.e. independent of $Q^{2}$, consistent with the ultra-local
character of the scalar Yukawa theory. The total
$\bar{c}(Q^{2})=\bar{c}_{N}(Q^{2})+\bar{c}_{\pi}(Q^{2})$, including the
divergences, vanishes with both the dimensional (DR) and Pauli-Villars (PV)
regularizations, as required by current conservation.

Comparing the results from the dimensional regularization and the Pauli-
Villars regularization, it is readily established that the Pauli-Villars
results can be obtained from the dimensional regularization with the following
replacement,

$\frac{1}{\epsilon}\;\leftrightarrow\;\ln\frac{\mu^{2}_{\textsc{pv}}}{\mu_{R}^{2}}-\frac{3}{2},$
(29)

in the physical limit ($\mu_{\textsc{pv}}\to\infty$ and $\epsilon\to 0$).

###  III.2 Light-cone perturbation theory

We adopt the Pauli-Villars regularization for calculations in LCPT [50, 51,
52]. This regularization scheme can be extended to the non-perturbative light-
front Hamiltonian formalism including gauge theories [53, 54, 48, 55]. In LO
LCPT, the obtained counterterms are the same as the counterterms in covariant
perturbation theory with Pauli-Villars regularization, i.e. Eqs. (5–8).

For evaluating the HMEs of the EMT, we adopt the Drell-Yan frame, in which the
longitudinal momentum of the graviton vanishes $q^{+}=0$ [56]. The vacuum
creation and annihilation are suppressed in this frame. Mathematically,
$q^{+}=0$ is equivalent to integrating out the longitudinal $x^{-}$ direction
on the light cone. In the Drell-Yan frame, the diagrams needed for the EMT are
the same as Fig. 2, except all diagrams are now ordered in light-front time
$x^{+}$. Applying the diagrammatic rules [11], we obtain the HMEs as,

$\displaystyle t_{a}^{\mu\nu}=\,$ $\displaystyle
Z\left[2P^{\mu}P^{\nu}-\frac{1}{2}(q^{\mu}q^{\nu}-q^{2}g^{\mu\nu})\right],$
(30) $\displaystyle t_{b}^{\mu\nu}=\,$ $\displaystyle-
t_{c}^{\mu\nu}=-t_{d}^{\mu\nu}=-Zg^{\mu\nu}\Sigma(m^{2}),$ (31) $\displaystyle
t_{e}^{\mu\nu}=\,$ $\displaystyle
Zg_{0}^{2}\int_{0}^{1}\frac{{\mathrm{d}}x}{2x(1-x)}\int\frac{{\mathrm{d}}^{2}k_{\perp}}{(2\pi)^{3}}\sum_{j}(-1)^{j}\frac{\frac{{1}}{2}q^{2}g^{\mu\nu}+p_{2}^{\\{\mu}p_{2}^{\prime\nu\\}}}{(1-x)(s_{2}^{j}-m^{2})(s_{2a}^{j}-m^{2})},$
(32) $\displaystyle t_{f}^{\mu\nu}=\,$ $\displaystyle
Zg_{0}^{2}\int_{0}^{1}\frac{{\mathrm{d}}x}{2x(1-x)}\int\frac{{\mathrm{d}}^{2}k_{\perp}}{(2\pi)^{3}}\sum_{j}(-1)^{j}\frac{\frac{1}{2}q^{2}g^{\mu\nu}+p_{1j}^{\\{\mu}p_{1j}^{\prime\nu\\}}}{x(s^{j}_{2}-m^{2})(s^{j}_{2b}-m^{2})},$
(33)

where,

$\displaystyle p^{\mu}_{1j}=\,$
$\displaystyle\big{(}xP^{+},\frac{\vec{p}_{1j\perp}^{2}+\mu^{2}_{j}}{p_{1j}^{+}},\vec{k}_{\perp}+x\vec{P}_{\perp}-\frac{1}{2}x\vec{q}_{\perp}\big{)},$
(34) $\displaystyle p^{\prime\mu}_{1j}=\,$
$\displaystyle\big{(}xP^{+},\frac{\vec{p}_{1j\perp}^{\prime
2}+\mu^{2}_{j}}{p_{1j}^{\prime+}},\vec{k}_{\perp}+x\vec{P}_{\perp}+(1-\frac{1}{2}x)\vec{q}_{\perp}\big{)},$
(35) $\displaystyle p^{\mu}_{2}=\,$
$\displaystyle\big{(}(1-x)P^{+},\frac{\vec{p}_{2\perp}^{2}+m^{2}}{p_{2}^{+}},-\vec{k}_{\perp}+(1-x)\vec{P}_{\perp}-\frac{1}{2}(1-x)\vec{q}_{\perp}\big{)},$
(36) $\displaystyle p^{\prime\mu}_{2}=\,$
$\displaystyle\big{(}(1-x)P^{+},\frac{\vec{p}_{2\perp}^{\prime
2}+m^{2}}{p_{2}^{\prime+}},-\vec{k}_{\perp}+(1-x)\vec{P}_{\perp}+\frac{1}{2}(1+x)\vec{q}_{\perp}\big{)}.$
(37)

The $s$ factors in the energy denominators are,

$\displaystyle s_{2}^{j}$
$\displaystyle=\frac{\vec{k}_{\perp}^{2}+\mu_{j}^{2}}{x}+\frac{\vec{k}_{\perp}^{2}+m^{2}}{1-x},$
(38) $\displaystyle s_{2a}^{j}$
$\displaystyle=\frac{(\vec{k}_{\perp}-x\vec{q}_{\perp})^{2}+\mu_{j}^{2}}{x}+\frac{(\vec{k}_{\perp}-x\vec{q}_{\perp})^{2}+m^{2}}{1-x},$
(39) $\displaystyle s^{j}_{2b}$
$\displaystyle=\frac{[\vec{k}_{\perp}+(1-x)\vec{q}_{\perp}]^{2}+\mu_{j}^{2}}{x}+\frac{[\vec{k}_{\perp}+(1-x)\vec{q}_{\perp}]^{2}+m^{2}}{1-x}.$
(40)

As mentioned, in Pauli-Villars regularization, we introduce $\mu_{0}=\mu$ and
$\mu_{1}=\mu_{\textsc{pv}}$.

From these expressions, it is clear that the light-front results are not
manifestly covariant and different components have to be evaluated separately.
In our previous work, we proposed to use $t^{++}_{i}$ to extract GFF $A_{i}$
and to use $t^{++}_{i},t^{+-}_{i}$ to extract GFF $D_{i}$ [19, 57]. Using the
perturbative results obtained above, one readily shows that $t^{++}_{i}$ and
$t^{+-}_{i}$ from LCPT and from the covariant perturbation theory are
identical diagram by diagram, which validates our previous proposal in the
context of perturbation theory. Unfortunately, this equivalence could not be
extended to all components of the EMT. Indeed, one can check that current is
not conserved, viz.

$q_{\nu}t^{1\nu}\neq 0\,.$ (41)

Since individually,

$q_{\nu}t^{1\nu}_{i}=2m^{2}q_{\perp}^{1}\bar{c}_{i}(-q^{2})\,,$ (42)

the violation of the current conservation leads to a non-vanishing total
$\bar{c}$:

$\bar{c}(-q^{2})\equiv\sum_{i}\bar{c}_{i}(-q^{2})=\frac{Z}{m^{2}}\Sigma(m^{2})\,.$
(43)

Note that $\bar{c}$ is also divergent. Fortunately, we can compare with the
covariant results diagram by diagram to pinpoint the sources of the violation.
We find that the violation of the current conservation originates from
$t^{11}_{i},t^{22}_{i}$ with contributions from diagrams ($e$) and ($f$). The
origin of the discrepancy is that the counterterms derived from renormalizing
the Hamiltonian operator are not sufficient to remove the divergences in
$t^{11},t^{22}$. This is in contrast to $t^{+-}$ which in fact involves the
same operator as the Hamiltonian. On the other hand,
$t^{+1}_{i},t^{+2}_{i},t^{12}_{i}$ are identical to the respective covariant
results. The minus HMEs, $t^{-i}_{i},t^{--}_{i}$ are very singular and cannot
be regularized with a single Pauli-Villars pion. More copies of the Pauli-
Villars ghosts are needed to tame these components [47].

In summary, $t^{++}_{i},t^{+-}_{i},t^{+1}_{i},t^{+2}_{i},t^{12}_{i}$ on the
light front are identical to their counterparts in covariant perturbation
theory, while $t^{11}_{i},t^{22}_{i}$ are not. $t^{11}_{i},t^{22}_{i}$ also
contain uncanceled divergences. We can use these “good currents”
$t^{++}_{i},t^{+-}_{i},t^{+1}_{i},t^{+2}_{i},t^{12}_{i}$ to extract the GFFs,

$\displaystyle t_{i}^{++}$ $\displaystyle=2(P^{+})^{2}A_{i}(Q^{2}),$ (44)
$\displaystyle t_{i}^{+a}$
$\displaystyle=2P^{+}P^{a}_{\perp}A_{i}(Q^{2})+\frac{1}{2}q^{+}q_{\perp}^{a}D_{i}(Q^{2}),$
(45) $\displaystyle t_{i}^{12}$
$\displaystyle=2P^{1}_{\perp}P^{2}_{\perp}A_{i}(Q^{2})+\frac{1}{2}q^{1}_{\perp}q_{\perp}^{2}D_{i}(Q^{2}),$
(46) $\displaystyle t_{i}^{+-}$
$\displaystyle=2\big{(}P^{2}_{\perp}+M^{2}+\frac{1}{4}q^{2}_{\perp}\big{)}A_{i}(Q^{2})+{q}_{\perp}^{2}D_{i}(Q^{2})+4M^{2}\bar{c}_{i}(Q^{2}),$
(47)

where, $a=1,2$, $i=N,\pi$ and $Q^{2}=-q^{2}$. We can use the combinations of
these HMEs to extract the GFFs $A_{i},D_{i},\bar{c}_{i}$. The above
expressions are further simplified in the transverse Breit frame with
$\vec{P}_{\perp}=0$. Note that to use $t^{+a}$ to extract $D_{i}$, we need to
compute the EMT in a general frame $q^{+}\neq 0$ and take the Drell-Yan limit
$q^{+}\to 0$ in the end. Since $\bar{c}_{i}$ is only contained in $t^{+-}_{i}$
for the above set of currents, this HME is indispensable to obtain the forces
between the nucleon and the pion.

##  IV Covariant light-front dynamics

The conventional wisdom regarding currents in light-front physics is that one
should use the “good currents” while avoiding the “bad” and “terrible
currents” if possible [10, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77]. The good currents are those associated with the
plus component of the current, e.g. $J^{+}$. The bad currents are currents
with transverse components, i.e. $\vec{J}_{\perp}$. And the terrible (or
worst) currents are those associated with the minus component of the current,
e.g. $J^{-}$ [68]. Quantitatively, one can define $\gamma$ the “goodness” of
an operator $O$ as its scaling dimension under the longitudinal boost, viz
[78].

$\big{[}K_{3},O\big{]}=-i\gamma O,$ (48)

where, $K_{3}$ is the generator of the longitudinal boost. From the infinite-
momentum perspective, good currents are those that increase as $P_{z}$
increases, i.e. relevant. Bad and terrible currents are those that decrease as
$P_{z}$ increases [79, 68]. Our LCPT calculation certainly follows the
conventional wisdom. On the other hand, the traditional good currents $t^{++}$
(goodness 2) and $t^{+a}$ (goodness 1) are not sufficient to extract the full
set of GFFs. And the goodness of the currents does not distinguish the UV
divergent $t^{11},t^{22}$ from the UV finite $t^{+-}$, $t^{12}$, all of which
have goodness 0.

A systematic way to analyze the current components based on symmetry is the
covariant light-front dynamics proposed by Karmanov et. al [80, 81, 82, 83,
84, 12, 85, 86, 87, 88]. The rationale of this method is as follows. The
Lorentz decomposition of the HME is based on the Poincaré symmetry. In light-
front dynamics, some of the Poincaré symmetries are not manifest [90]. In
practical calculations, these symmetries are likely broken. As a result, the
HMEs of the EMT have to be re-parametrized with the reduced symmetries to
account for the possible breaking of the Poincaré symmetries. And the full
Poincaré symmetry is restored in the continuum limit. An effective way to
parametrize the HME with the reduced symmetries is to include a null vector
$\omega^{\mu}$ ($\omega_{\mu}\omega^{\mu}=0$), which is perpendicular to the
light front, i.e. $\omega\cdot x=0$ [83, 84, 12, 85, 86, 87, 88]. Standard
light-front coordinates correspond to the choice
$\omega^{\mu}=(\omega^{+},\omega^{-},\vec{\omega}_{\perp})=(0,2,0)$. There is
a conformal symmetry related to the scaling of $\omega$. Namely, physics does
not change if we multiply $\omega$ with an arbitrary non-zero real number.
This symmetry requires that the $\omega^{\mu}$ dependence always appears in
the form $\omega^{\mu}/(\omega\cdot P)$ [12].

The presence of $\omega$ also adds one more scalar $\zeta=(\omega\cdot
q)/(\omega\cdot P)$ into the calculation. Now, the form factors depend on two
Lorentz scalars, $q^{2}$ and $\zeta$. It can be shown that Hermiticity of the
EMT operator leads to $F_{i}(\zeta,q^{2})=F^{*}_{i}(-\zeta,q^{2})$, where
$F_{i}$ is a form factor [28]. In the Drell-Yan frame $\zeta=0$ and the form
factors become real functions that depend on $q^{2}$ only. This choice
dramatically simplifies the Lorentz decomposition in covariant light-front
dynamics [85]. We will adopt the Drell-Yan frame unless otherwise stated.
Then, one of the most general Lorentz decompositions of the HMEs of the EMT
operator for a spin-0 hadron on the light front in the Drell-Yan frame
$q^{+}=0$ is111This decomposition differs from our previous work [19] by
replacing tensor $V^{\alpha}V^{\beta}$ with
$V^{\alpha}V^{\beta}+q^{\alpha}q^{\beta}$. Another minor difference is to
replace ${q^{4}\omega^{\alpha}\omega^{\beta}}/{(\omega\cdot P)^{2}}$ with
${M^{4}\omega^{\alpha}\omega^{\beta}}/{(\omega\cdot P)^{2}}$. The physical
difference of these two choices will be analyzed below.,

$\langle
p^{\prime}|T_{i}^{\alpha\beta}(0)|p\rangle=2P^{\alpha}P^{\beta}A_{i}(-q^{2})+\frac{1}{2}(q^{\alpha}q^{\beta}-q^{2}g^{\alpha\beta})D_{i}(-q^{2})+2M^{2}g^{\alpha\beta}\bar{c}_{i}(-q^{2})\\\
+\frac{M^{4}\omega^{\alpha}\omega^{\beta}}{(\omega\cdot
P)^{2}}S_{1i}(-q^{2})+(V^{\alpha}V^{\beta}+q^{\alpha}q^{\beta})S_{2i}(-q^{2}),$
(49)

where, $P=(p^{\prime}+p)/2$, $q=p^{\prime}-p$. Vector $V^{\alpha}$ is defined
as
$V^{\alpha}=\varepsilon^{\alpha\beta\rho\sigma}P_{\beta}q_{\rho}\omega_{\sigma}/(\omega\cdot
P)$. This vector is perpendicular to $q^{\alpha}$, $P^{\alpha}$ as well as
$\omega^{\alpha}$. Hermiticity excludes terms like
$P^{\\{\alpha}q^{\beta\\}}$, $\omega^{\\{\alpha}q^{\beta\\}}$ or
$\omega^{\\{\alpha}V^{\beta\\}}$. Form factors $S_{1i}$ and $S_{2i}$ are
allowed by the reduced kinematic symmetries but forbidden by the full Poincaré
symmetry. Therefore, they are called the spurious form factors. They are
expected to vanish in the continuum limit. Using this decomposition Eq. (49),
we can express the HMEs of the currents in terms of the GFFs as,

$\displaystyle A_{i}=\,$ $\displaystyle\frac{1}{2(\omega\cdot
P)^{2}}\omega_{\alpha}\omega_{\beta}t^{\alpha\beta}_{i},$ (50) $\displaystyle
D_{i}=\,$
$\displaystyle\frac{2}{q^{4}}\big{(}q_{\alpha}q_{\beta}-V_{\alpha}V_{\beta})t^{\alpha\beta}_{i},$
(51) $\displaystyle\bar{c}_{i}=\,$
$\displaystyle\frac{1}{4M^{2}q^{2}}\Big{\\{}(q_{\alpha}q_{\beta}-3V_{\alpha}V_{\beta}+q^{2}g_{\alpha\beta})t^{\alpha\beta}_{i}-(M^{2}-\frac{1}{4}q^{2})\frac{q^{2}}{(\omega\cdot
P)^{2}}\omega_{\alpha}\omega_{\beta}t^{\alpha\beta}_{i}\Big{\\}},$ (52)
$\displaystyle S_{1}=\,$
$\displaystyle\frac{1}{2M^{4}q^{2}}\Big{[}(M^{2}-\frac{1}{4}q^{2})(V_{\alpha}V_{\beta}+q_{\alpha}q_{\beta}-q^{2}g_{\alpha\beta})t_{i}^{\alpha\beta}+2q^{2}P_{\alpha}P_{\beta}t_{i}^{\alpha\beta}-(M^{2}-\frac{1}{4}q^{2})^{2}\frac{q^{2}}{(\omega\cdot
P)^{2}}\omega_{\alpha}\omega_{\beta}t_{i}^{\alpha\beta}\Big{]},$ (53)
$\displaystyle S_{2}=\,$
$\displaystyle\frac{1}{2q^{4}}\Big{[}(q_{\alpha}q_{\beta}-q^{2}g_{\alpha\beta}+3V_{\alpha}V_{\beta})t_{i}^{\alpha\beta}+(M^{2}-\frac{1}{4}q^{2})\frac{q^{2}}{(\omega\cdot
P)^{2}}\omega_{\alpha}\omega_{\beta}t_{i}^{\alpha\beta}\Big{]}.$ (54)

It is more transparent to write down the relations between the GFFs and the
current components in standard light-front dynamics. These relations are,

$\displaystyle t^{++}_{i}=2(P^{+})^{2}A_{i}(q^{2}_{\perp}),$ (55)
$\displaystyle t^{+a}=2P^{+}P^{a}_{\perp}A_{i}(q^{2}_{\perp}),$ (56)
$\displaystyle
t^{ab}_{i}=2P^{a}_{\perp}P^{b}_{\perp}A_{i}(q^{2}_{\perp})+\frac{1}{2}(q^{a}_{\perp}q^{b}_{\perp}-\delta^{ab}q^{2}_{\perp})D_{i}(q^{2}_{\perp})-2M^{2}\delta^{ab}\bar{c}_{i}(q^{2}_{\perp})+\delta^{ab}q_{\perp}^{2}S_{2i}(q_{\perp}^{2}),$
(57) $\displaystyle
t^{+-}_{i}=2(M^{2}+P_{\perp}^{2}+\frac{1}{4}q_{\perp}^{2})A_{i}(q^{2}_{\perp})+q^{2}_{\perp}D_{i}(q^{2}_{\perp})+4M^{2}\bar{c}_{i}(q_{\perp}^{2}),$
(58) $\displaystyle
t^{--}_{i}=2\Big{(}\frac{M^{2}+P^{2}_{\perp}+\frac{1}{4}q_{\perp}^{2}}{P^{+}}\Big{)}^{2}A_{i}(q_{\perp}^{2})+2\Big{(}\frac{\vec{q}_{\perp}\cdot\vec{P}_{\perp}}{P^{+}}\Big{)}^{2}D_{i}(q^{2}_{\perp})+\frac{4M^{4}}{(P^{+})^{2}}S_{1i}(q_{\perp}^{2})+\frac{4P^{2}_{\perp}q_{\perp}^{2}}{(P^{+})^{2}}S_{2i}(q_{\perp}^{2}),$
(59) $\displaystyle
t^{-a}_{i}=\frac{M^{2}+P^{2}_{\perp}+\frac{1}{4}q_{\perp}^{2}}{P^{+}}2P^{a}_{\perp}A_{i}(q_{\perp}^{2})+\frac{\vec{q}_{\perp}\cdot\vec{P}_{\perp}}{P^{+}}q^{a}_{\perp}D_{i}(q^{2}_{\perp})+\frac{2q^{2}_{\perp}}{P^{+}}P^{a}_{\perp}S_{2i}(q^{2}_{\perp}).$
(60)

Here, $a=1,2$ is the transverse index and $i=N,\pi$. Further simplification
can be achieved by taking the Breit frame,
$\vec{P}_{\perp}=\frac{1}{2}(\vec{p}_{\perp}+\vec{p}^{\prime}_{\perp})=0$,

$\displaystyle t^{++}_{i}=2(P^{+})^{2}A_{i}(q^{2}_{\perp}),$ (61)
$\displaystyle t^{+a}_{i}=0,$ (62) $\displaystyle
t^{ab}_{i}=\frac{1}{2}(q^{a}_{\perp}q^{b}_{\perp}-\delta^{ab}q^{2}_{\perp})D_{i}(q^{2}_{\perp})-2M^{2}\delta^{ab}\bar{c}_{i}(q^{2}_{\perp})+\delta^{ab}q_{\perp}^{2}S_{2i}(q_{\perp}^{2}),$
(63) $\displaystyle
t^{+-}_{i}=2(M^{2}+\frac{1}{4}q_{\perp}^{2})A_{i}(q^{2}_{\perp})+q^{2}_{\perp}D_{i}(q^{2}_{\perp})+4M^{2}\bar{c}_{i}(q^{2}_{\perp}),$
(64) $\displaystyle
t^{--}_{i}=2\Big{(}\frac{M^{2}+\frac{1}{4}q_{\perp}^{2}}{P^{+}}\Big{)}^{2}A_{i}(q_{\perp}^{2})+\frac{4M^{4}}{(P^{+})^{2}}S_{1i}(q_{\perp}^{2}),$
(65) $\displaystyle t^{-a}_{i}=0,$ (66) $\displaystyle\Rightarrow\quad$
$\displaystyle\mathrm{tr}_{\perp}\tensor{t}_{i}=t^{11}_{i}+t^{22}_{i}=-\frac{1}{2}q^{2}_{\perp}D_{i}(q^{2}_{\perp})-4M^{2}\bar{c}_{i}(q^{2}_{\perp})+2q_{\perp}^{2}S_{2i}(q_{\perp}^{2}),$
(67) $\displaystyle
t^{12}_{i}=\frac{1}{2}q^{1}_{\perp}q^{2}_{\perp}D_{i}(q^{2}_{\perp}).$ (68)

Note that there are 5 independent HMEs and 5 GFFs. The physical GFFs
$A_{i},D_{i},\bar{c}_{i}$ can only be extracted from combinations of the three
“good currents” $t^{++}_{i}$, $t^{+-}_{i}$ and $t^{12}_{i}$, without the
contamination of the spurious form factors $S_{1,2}$. This observation is
consistent with the perturbative results [19]. For the total GFFs, in order
for the Newton’s third law $\sum_{i}\bar{c}_{i}(Q^{2})=0$ to hold, a non-
trivial relation between these three good HMEs,

$\frac{M^{2}+\frac{1}{4}\vec{q}^{2}_{\perp}}{(P^{+})^{2}}t^{++}+2\frac{\vec{q}^{2}_{\perp}}{q^{1}_{\perp}q^{2}_{\perp}}t^{12}=t^{+-},$
(69)

has to be satisfied. We note that this relation is not automatically
satisfied, except for $q^{2}_{\perp}=0$. Therefore, Newton’s third law (16)
holds when integrated over space on the light front $x^{+}=0$. For the scalar
Yukawa model, $\bar{c}_{i}(Q^{2})$ are constants in one-loop order. Therefore,
Newton’s third law (16) is satisfied everywhere in this case.

Decomposition (49) is not unique, and different forms are linearly related to
each other with redefinitions of the form factors. There is still
arbitrariness in the Lorentz decomposition even if the physical part of the
Lorentz structures is fixed [28, 89, 19]. For example, in our previous work
[19], we adopted the Lorentz tensor $V^{\alpha}V^{\beta}$ instead of
$V^{\alpha}V^{\beta}+q^{\alpha}q^{\beta}$ in Eq. (49). Yet another possible
tensor is $P^{\\{\alpha}\omega^{\beta\\}}/(\omega\cdot P)$ [28], which can be
expressed in terms of the present tensors as,

$\frac{P^{\\{\alpha}\omega^{\beta\\}}}{\omega\cdot
P}=g^{\alpha\beta}+(M^{2}-\frac{1}{4}q^{2})\frac{\omega^{\alpha}\omega^{\beta}}{(\omega\cdot
P)^{2}}-\frac{V^{\alpha}V^{\beta}+q^{\alpha}q^{\beta}}{q^{2}}.$ (70)

Obviously, different decompositions are related to each other by redefinition
of the spurious form factors $S_{1i}(-q^{2})$ and $S_{2i}(-q^{2})$ as well as
the physical form factors $A_{i}(-q^{2})$, $D_{i}(-q^{2})$,
$\bar{c}_{i}(-q^{2})$. Note that the physical GFFs may also get redefined,
which is in contrast to the vector current $J^{\mu}$ where the physical and
spurious form factors do not mix. Since the spurious form factors in general
carry uncanceled divergences, this redefinition also shuffles around the
divergences.

One of the pressing questions is how to choose a proper decomposition scheme
such that the GFFs associated with the physical tensors, i.e.
$P^{\alpha}P^{\beta}$,$q^{\alpha}q^{\beta}-q^{2}g^{\alpha\beta}$ and
$g^{\alpha\beta}$, correspond to the genuine physical GFFs. In a non-
perturbative case, there is no easy answer to this problem, since even the
physical form factors are only approximations of the true physical ones,
unless a perturbative matching can be achieved. Fortunately, in the
perturbative case, we can use the covariant results as a benchmark. It is
evident that (49) is the only decomposition scheme that is consistent with the
LO perturbation theory shown in Sect. III.2, cf. Eqs. (55–58). The matching to
the perturbative results also guarantees the exact cancellation of divergences
in the physical GFFs. And, all divergences are confined within the spurious
form factors $S_{1i}(-q^{2})$ and $S_{2i}(-q^{2})$. In particular, adoption of
the spurious tensor $V^{\alpha}V^{\beta}+q^{\alpha}q^{\beta}$ ensures that the
light-front energy density $T^{+-}$ is divergence free, which is consistent
with our previous analysis using light-front energy conservation [19].

Another interesting observation is that this tensor structure implies a hidden
symmetry. The spurious terms within $t_{i}^{\mu\nu}$ depend only on
$q^{2}_{\perp}$,

$\displaystyle t^{+-}_{\text{spur}}(q^{2}_{\perp})=\,$ $\displaystyle
t^{12}_{\text{spur}}(q^{2}_{\perp})=0,$ (71) $\displaystyle
t^{11}_{\text{spur}}(q_{\perp}^{2})=\,$ $\displaystyle
t^{22}_{\text{spur}}(q_{\perp}^{2})=q_{\perp}^{2}S_{2i}(q^{2}_{\perp}),$ (72)

i.e. they are isotropic in the transverse plane.When the theory is quantized
on the light cone, a special direction, $z$ direction in standard light front
coordinates, is introduced. This choice underlines the rise of the spurious
terms. However, the rotational symmetry in the transverse plane is kept and
there is no potential source, kinematical or dynamical, to break this
symmetry. Therefore, the above observation is physically plausible. By
contrast, the spurious tensor $V^{\alpha}V^{\beta}$ is anisotropic in the
transverse direction. This symmetry suggests that our decomposition (49) works
beyond the leading-order perturbation theory.

##  V Summary and outlooks

In this work, we compare the gravitational form factors of the scalar Yukawa
theory obtained from covariant perturbation theory and from light-cone
perturbation theory at the one-loop level. Our intention is to find a reliable
way to extract the GFFs in light-front dynamics beyond the perturbative
regime. In the light-front formalism, the full Poincaré symmetry is not
manifest. As such, GFFs extracted from some current components of the EMT
contain uncanceled divergences, which further lead to violation of
conservation laws and known sum rules. Traditionally, the “good currents” are
used to extract physical form factors. In the case of the GFFs, the good
currents $T^{++}$ (and $T^{+a}$) are not sufficient to obtain the $D$-term as
well as the $\bar{c}$-term, i.e. the stress within hadrons. Using perturbation
theory as a benchmark, we found that hadronic matrix elements $t^{++}_{i}$,
$t^{+a}_{i}$, $t^{12}_{i}$ and $t^{+-}_{i}$ in light-cone perturbation theory
reproduce those obtained from covariant perturbation theory whereas
$t^{11}_{i}$, $t^{22}_{i}$ contains uncanceled divergences. Based on these
analysis, we proposed a new covariant Lorentz decomposition of the energy-
momentum tensor on the light front (49), which takes into account the presence
of a special direction $\omega^{\mu}$ when quantizing the theory.

The new CLFD decomposition avoids the mixing of the uncanceled spurious
divergence within the physical GFFs. In particular, the light-front energy
density $T^{+-}$ does not contain spurious contributions, which is consistent
with our previous derivation that this operator can be properly renormalized
using the same set of counterterms that renormalize the light-front
Hamiltonian [19]. Furthermore, the new decomposition also introduces a new
symmetry for the spurious terms on the transverse plane, viz. the breaking of
the Poincaré symmetry is isotropic in the transverse direction.

One immediate application is the same system, the scalar Yukawa theory, in the
non-perturbative regime. The present work overcomes the divergence issues
encountered in the mass decomposition of the strongly coupled scalar nucleon
[19]. While the decomposition Eq. (49) differs from Eq. (69) of our previous
work, it can be shown that for the total GFF $D(Q^{2})$ extracted from
$T^{+-}$, these two decompositions give the same result.

The method presented in this work is very general and can also be applied to
other systems [91, 29, 30, 92, 31, 32, 33, 34, 93, 94] and to hadrons with
spin. One of the prime systems of interest is the nucleon. True nucleons are
spin-1/2 hadrons. And they have one more set of GFFs $J_{i}$ (or $B_{i}$),
representing the spin distribution inside the system [2]. On the other hand,
there are also spin-flip hadron matrix elements. Some of the spurious Lorentz
structures contain zero-model contributions [30]. Therefore, it is interesting
to see whether covariant light-front dynamics can capture at least some of the
zero-mode effects.

## Acknowledgements

The authors acknowledge fruitful discussions with J. More, C. Mondal, S. Nair,
and K. Serafin.

This work was supported in part by the National Natural Science Foundation of
China (NSFC) under Grant No. 12375081, by the Chinese Academy of Sciences
under Grant No. YSBR-101, and by the US Department of Energy (DOE) under Grant
No. DE-SC0023692.

X. Z. is supported by new faculty startup funding by the Institute of Modern
Physics, Chinese Academy of Sciences, by Key Research Program of Frontier
Sciences, Chinese Academy of Sciences, Grant No. ZDBS-LY-7020, by the
Foundation for Key Talents of Gansu Province, by the Central Funds Guiding the
Local Science and Technology Development of Gansu Province, Grant No.
22ZY1QA006, by Gansu International Collaboration and Talents Recruitment Base
of Particle Physics (2023-2027), by International Partnership Program of the
Chinese Academy of Sciences, Grant No. 016GJHZ2022103FN, by National Natural
Science Foundation of China, Grant No. 12375143, by National Key R&D Program
of China, Grant No. 2023YFA1606903 and by the Strategic Priority Research
Program of the Chinese Academy of Sciences, Grant No. XDB34000000. V.A.K. is
supported by the Chinese Academy of Sciences President’s International
Fellowship Initiative Grant No. 2023VMA0010.

## References

  * [1] O. V. Teryaev, Front. Phys. (Beijing) 11, no.5, 111207 (2016) doi:10.1007/s11467-016-0573-6
  * [2] M. V. Polyakov and P. Schweitzer, Int. J. Mod. Phys. A 33, no.26, 1830025 (2018) doi:10.1142/S0217751X18300259 [arXiv:1805.06596 [hep-ph]].
  * [3] V. D. Burkert, L. Elouadrhiri, F. X. Girod, C. Lorcé, P. Schweitzer and P. E. Shanahan, Rev. Mod. Phys. 95, no.4, 041002 (2023) doi:10.1103/RevModPhys.95.041002 [arXiv:2303.08347 [hep-ph]].
  * [4] I. Y. Kobzarev and L. B. Okun, Zh. Eksp. Teor. Fiz. 43, 1904-1909 (1962)
  * [5] H. Pagels, Phys. Rev. 144, 1250-1260 (1966) doi:10.1103/PhysRev.144.1250
  * [6] X. D. Ji, J. Tang and P. Hoodbhoy, Phys. Rev. Lett. 76, 740-743 (1996) doi:10.1103/PhysRevLett.76.740 [arXiv:hep-ph/9510304 [hep-ph]].
  * [7] X. D. Ji, Phys. Rev. Lett. 78, 610-613 (1997) doi:10.1103/PhysRevLett.78.610 [arXiv:hep-ph/9603249 [hep-ph]].
  * [8] J. M. Namyslowski, Prog. Part. Nucl. Phys. 14, 49-104 (1985) doi:10.1016/0146-6410(85)90051-1
  * [9] W. M. Zhang, Chin. J. Phys. 32, 717-808 (1994) [arXiv:hep-ph/9412244 [hep-ph]].
  * [10] M. Burkardt, Adv. Nucl. Phys. 23, 1-74 (1996) doi:10.1007/0-306-47067-5_1 [arXiv:hep-ph/9505259 [hep-ph]].
  * [11] S. J. Brodsky, H. C. Pauli and S. S. Pinsky, Phys. Rept. 301, 299-486 (1998) doi:10.1016/S0370-1573(97)00089-6 [arXiv:hep-ph/9705477 [hep-ph]].
  * [12] J. Carbonell, B. Desplanques, V. A. Karmanov and J. F. Mathiot, Phys. Rept. 300, 215-347 (1998) doi:10.1016/S0370-1573(97)00090-2 [arXiv:nucl-th/9804029 [nucl-th]].
  * [13] T. Heinzl, Lect. Notes Phys. 572, 55-142 (2001) doi:10.1007/3-540-45114-5_2 [arXiv:hep-th/0008096 [hep-th]].
  * [14] G. A. Miller, Prog. Part. Nucl. Phys. 45, 83-155 (2000) doi:10.1016/S0146-6410(00)00103-4 [arXiv:nucl-th/0002059 [nucl-th]].
  * [15] B. L. G. Bakker, A. Bassetto, S. J. Brodsky, W. Broniowski, S. Dalley, T. Frederico, S. D. Glazek, J. R. Hiller, C. R. Ji and V. Karmanov, et al. Nucl. Phys. B Proc. Suppl. 251-252, 165-174 (2014) doi:10.1016/j.nuclphysbps.2014.05.004 [arXiv:1309.6333 [hep-ph]].
  * [16] J. R. Hiller, Prog. Part. Nucl. Phys. 90, 75-124 (2016) doi:10.1016/j.ppnp.2016.06.002 [arXiv:1606.08348 [hep-ph]].
  * [17] S. J. Brodsky, A. Deur and C. D. Roberts, Nature Rev. Phys. 4, no.7, 489-495 (2022) doi:10.1038/s42254-022-00453-3 [arXiv:2202.06051 [hep-ph]].
  * [18] F. Gross, E. Klempt, S. J. Brodsky, A. J. Buras, V. D. Burkert, G. Heinrich, K. Jakobs, C. A. Meyer, K. Orginos and M. Strickland, et al. Eur. Phys. J. C 83, 1125 (2023) doi:10.1140/epjc/s10052-023-11949-2 [arXiv:2212.11107 [hep-ph]].
  * [19] X. Cao, Y. Li and J. P. Vary, Phys. Rev. D 108, no.5, 056026 (2023) doi:10.1103/PhysRevD.108.056026 [arXiv:2308.06812 [hep-ph]].
  * [20] Y. Li, V. A. Karmanov, P. Maris and J. P. Vary, Few Body Syst. 56, no.6-9, 495-501 (2015) doi:10.1007/s00601-015-0965-0 [arXiv:1411.1707 [nucl-th]].
  * [21] Y. Li, V. A. Karmanov, P. Maris and J. P. Vary, Phys. Lett. B 748, 278-283 (2015) doi:10.1016/j.physletb.2015.07.014 [arXiv:1504.05233 [nucl-th]].
  * [22] Y. Li,
  * [23] V. A. Karmanov, Y. Li, A. V. Smirnov and J. P. Vary, Phys. Rev. D 94, no.9, 096008 (2016) doi:10.1103/PhysRevD.94.096008 [arXiv:1610.03559 [hep-th]].
  * [24] X. D. Ji, Phys. Rev. Lett. 74, 1071-1074 (1995) doi:10.1103/PhysRevLett.74.1071 [arXiv:hep-ph/9410274 [hep-ph]].
  * [25] X. D. Ji, Phys. Rev. D 52, 271-281 (1995) doi:10.1103/PhysRevD.52.271 [arXiv:hep-ph/9502213 [hep-ph]].
  * [26] C. Lorcé, H. Moutarde and A. P. Trawiński, Eur. Phys. J. C 79, no.1, 89 (2019) doi:10.1140/epjc/s10052-019-6572-3 [arXiv:1810.09837 [hep-ph]].
  * [27] C. Lorcé, A. Metz, B. Pasquini and S. Rodini, JHEP 11, 121 (2021) doi:10.1007/JHEP11(2021)121 [arXiv:2109.11785 [hep-ph]].
  * [28] C. Lorcé, JHEP 08, 045 (2015) doi:10.1007/JHEP08(2015)045 [arXiv:1502.06656 [hep-ph]].
  * [29] J. More, A. Mukherjee, S. Nair and S. Saha, Phys. Rev. D 105, no.5, 056017 (2022) doi:10.1103/PhysRevD.105.056017 [arXiv:2112.06550 [hep-ph]].
  * [30] J. More, A. Mukherjee, S. Nair and S. Saha, Phys. Rev. D 107, no.11, 116005 (2023) doi:10.1103/PhysRevD.107.116005 [arXiv:2302.11906 [hep-ph]].
  * [31] A. Metz, B. Pasquini and S. Rodini, Phys. Lett. B 820, 136501 (2021) doi:10.1016/j.physletb.2021.136501 [arXiv:2104.04207 [hep-ph]].
  * [32] A. Freese and W. Cosyn, Phys. Rev. D 106, no.11, 114014 (2022) doi:10.1103/PhysRevD.106.114014 [arXiv:2207.10788 [hep-ph]].
  * [33] A. Freese, A. Metz, B. Pasquini and S. Rodini, Phys. Lett. B 839, 137768 (2023) doi:10.1016/j.physletb.2023.137768 [arXiv:2212.12197 [hep-ph]].
  * [34] M. I. Eides, Eur. Phys. J. C 83, no.5, 356 (2023) doi:10.1140/epjc/s10052-023-11535-6 [arXiv:2301.02883 [hep-ph]].
  * [35] G. Baym, Phys. Rev. 117, 886-888 (1960) doi:10.1103/PhysRev.117.886
  * [36] T. Nieuwenhuis and J. A. Tjon, Phys. Rev. Lett. 77, 814-817 (1996) doi:10.1103/PhysRevLett.77.814 [arXiv:hep-ph/9606403 [hep-ph]].
  * [37] F. Gross, C. Savkli and J. Tjon, Phys. Rev. D 64, 076008 (2001) doi:10.1103/PhysRevD.64.076008 [arXiv:nucl-th/0102041 [nucl-th]].
  * [38] D. S. Hwang and V. A. Karmanov, Nucl. Phys. B 696, 413-444 (2004) doi:10.1016/j.nuclphysb.2004.06.049 [arXiv:hep-th/0405035 [hep-th]].
  * [39] M. Emami-Razavi, N. Bergeron and J. W. Darewych, J. Phys. G 37, 025007 (2010) doi:10.1088/0954-3899/37/2/025007
  * [40] M. Emami-Razavi, N. Bergeron and J. W. Darewych, J. Phys. G 38, 065004 (2011) doi:10.1088/0954-3899/38/6/065004
  * [41] C. R. Ji and Y. Tokunaga, Phys. Rev. D 86, 054011 (2012) doi:10.1103/PhysRevD.86.054011 [arXiv:1205.3812 [nucl-th]].
  * [42] J. Carbonell, F. de Soto and V. A. Karmanov, Few Body Syst. 54, 2255-2269 (2013) doi:10.1007/s00601-012-0502-3 [arXiv:1211.5474 [nucl-th]].
  * [43] V. E. Rochev, J. Phys. A 46, 185401 (2013) doi:10.1088/1751-8113/46/18/185401 [arXiv:1302.1311 [hep-th]].
  * [44] V. Gigante, J. H. Alvarenga Nogueira, E. Ydrefors, C. Gutierrez, V. A. Karmanov and T. Frederico, Phys. Rev. D 95, no.5, 056012 (2017) doi:10.1103/PhysRevD.95.056012 [arXiv:1611.03773 [hep-ph]].
  * [45] S. S. Chabysheva and J. R. Hiller, Phys. Rev. D 105, no.5, 056027 (2022) doi:10.1103/PhysRevD.105.056027 [arXiv:2112.14950 [hep-ph]].
  * [46] S. S. Chabysheva and J. R. Hiller, Phys. Rev. D 105, no.11, 116006 (2022) doi:10.1103/PhysRevD.105.116006 [arXiv:2201.00123 [hep-th]].
  * [47] A. P. Bray, S. S. Chabysheva and J. R. Hiller, [arXiv:2401.00336 [hep-ph]].
  * [48] S. S. Chabysheva and J. R. Hiller, Phys. Rev. D 82, 034004 (2010) doi:10.1103/PhysRevD.82.034004 [arXiv:1006.1077 [hep-ph]].
  * [49] M. E. Peskin and D. V. Schroeder, Addison-Wesley, 1995, ISBN 978-0-201-50397-5
  * [50] S. J. Brodsky, J. R. Hiller and G. McCartor, Phys. Rev. D 64, 114023 (2001) doi:10.1103/PhysRevD.64.114023 [arXiv:hep-ph/0107038 [hep-ph]].
  * [51] S. J. Brodsky, J. R. Hiller and G. McCartor, Annals Phys. 305, 266-285 (2003) doi:10.1016/S0003-4916(03)00063-0 [arXiv:hep-th/0209028 [hep-th]].
  * [52] S. J. Brodsky, J. R. Hiller and G. McCartor, Annals Phys. 321, 1240-1264 (2006) doi:10.1016/j.aop.2005.09.005 [arXiv:hep-ph/0508295 [hep-ph]].
  * [53] S. S. Chabysheva and J. R. Hiller, Phys. Rev. D 79, 114017 (2009) doi:10.1103/PhysRevD.79.114017 [arXiv:0902.4443 [hep-ph]].
  * [54] S. S. Chabysheva and J. R. Hiller, Phys. Rev. D 81, 074030 (2010) doi:10.1103/PhysRevD.81.074030 [arXiv:0911.4455 [hep-ph]].
  * [55] S. S. Chabysheva and J. R. Hiller, [arXiv:1506.05429 [hep-th]].
  * [56] S. D. Drell and T. M. Yan, Phys. Rev. Lett. 24, 181-185 (1970) doi:10.1103/PhysRevLett.24.181
  * [57] S. Xu, X. Cao, T. Hu, Y. Li, X. Zhao and J. P. Vary, [arXiv:2404.06259 [hep-ph]].
  * [58] B. L. G. Bakker and C. R. Ji, Phys. Rev. D 65, 073002 (2002) doi:10.1103/PhysRevD.65.073002 [arXiv:hep-ph/0109005 [hep-ph]].
  * [59] S. J. Brodsky and D. S. Hwang, Nucl. Phys. B 543, 239-252 (1999) doi:10.1016/S0550-3213(98)00807-4 [arXiv:hep-ph/9806358 [hep-ph]].
  * [60] H. M. Choi and C. R. Ji, Phys. Rev. D 70, 053015 (2004) doi:10.1103/PhysRevD.70.053015 [arXiv:hep-ph/0402114 [hep-ph]].
  * [61] Y. Li, P. Maris and J. Vary, Phys. Rev. D 97, no.5, 054034 (2018) doi:10.1103/PhysRevD.97.054034 [arXiv:1712.03467 [hep-ph]].
  * [62] W. Jaus, Phys. Rev. D 41, 3394 (1990) doi:10.1103/PhysRevD.41.3394
  * [63] W. Jaus, Phys. Rev. D 60, 054026 (1999) doi:10.1103/PhysRevD.60.054026
  * [64] A. J. Arifi, H. M. Choi and C. R. Ji, Phys. Rev. D 108, no.1, 013006 (2023) doi:10.1103/PhysRevD.108.013006 [arXiv:2306.08536 [hep-ph]].
  * [65] S. Tang, S. Jia, P. Maris and J. P. Vary, Phys. Rev. D 104, no.1, 016002 (2021) doi:10.1103/PhysRevD.104.016002 [arXiv:2011.05454 [hep-ph]].
  * [66] M. Li, Y. Li, P. Maris and J. P. Vary, Phys. Rev. D 98, no.3, 034024 (2018) doi:10.1103/PhysRevD.98.034024 [arXiv:1803.11519 [hep-ph]].
  * [67] A. F. Krutov, R. G. Polezhaev and V. E. Troitsky, Phys. Rev. D 97, no.3, 033007 (2018) doi:10.1103/PhysRevD.97.033007 [arXiv:1801.01458 [hep-ph]].
  * [68] M. Strikman and C. Weiss, Phys. Rev. C 97, no.3, 035209 (2018) doi:10.1103/PhysRevC.97.035209 [arXiv:1706.02244 [hep-ph]].
  * [69] H. M. Choi and C. R. Ji, Phys. Rev. D 72, 013004 (2005) doi:10.1103/PhysRevD.72.013004 [arXiv:hep-ph/0504219 [hep-ph]].
  * [70] J. R. Cooke and G. A. Miller, Phys. Rev. C 66, 034002 (2002) doi:10.1103/PhysRevC.66.034002 [arXiv:nucl-th/0112037 [nucl-th]].
  * [71] A. F. Krutov and V. E. Troitsky, Phys. Rev. C 65, 045501 (2002) doi:10.1103/PhysRevC.65.045501 [arXiv:hep-ph/0101327 [hep-ph]].
  * [72] H. M. Choi and C. R. Ji, Phys. Rev. D 59, 034001 (1999) doi:10.1103/PhysRevD.59.034001 [arXiv:hep-ph/9807500 [hep-ph]].
  * [73] W. Melnitchouk, M. Sargsian and M. I. Strikman, Z. Phys. A 359, 99-109 (1997) doi:10.1007/s002180050372 [arXiv:nucl-th/9609048 [nucl-th]].
  * [74] J. Bolz and P. Kroll, Z. Phys. A 356, 327 (1996) doi:10.1007/s002180050186 [arXiv:hep-ph/9603289 [hep-ph]].
  * [75] M. Burkardt, Phys. Rev. D 54, 2913-2920 (1996) doi:10.1103/PhysRevD.54.2913 [arXiv:hep-ph/9601289 [hep-ph]].
  * [76] M. Burkardt, Phys. Rev. D 62, 071503 (2000) [erratum: Phys. Rev. D 66, 119903 (2002)] doi:10.1103/PhysRevD.62.071503 [arXiv:hep-ph/0005108 [hep-ph]].
  * [77] M. Vanderhaeghen, P. A. M. Guichon and M. Guidal, Phys. Rev. D 60, 094017 (1999) doi:10.1103/PhysRevD.60.094017 [arXiv:hep-ph/9905372 [hep-ph]].
  * [78] H. Leutwyler and J. Stern, Annals Phys. 112, 94 (1978) doi:10.1016/0003-4916(78)90082-9
  * [79] L. L. Frankfurt and M. I. Strikman, Phys. Rept. 160, 235-427 (1988) doi:10.1016/0370-1573(88)90179-2
  * [80] F. E. Close and H. Osborn, Phys. Lett. B 34, 400-404 (1971) doi:10.1016/0370-2693(71)90637-X
  * [81] H. Osborn, Nucl. Phys. B 38, 429-465 (1972) doi:10.1016/0550-3213(72)90319-7
  * [82] V. B. Berestetsky and M. V. Terentev, Sov. J. Nucl. Phys. 25, 347-354 (1977)
  * [83] V. A. Karmanov and A. V. Smirnov, Nucl. Phys. A 546, 691-717 (1992) doi:10.1016/0375-9474(92)90004-4
  * [84] V. A. Karmanov and J. F. Mathiot, Nucl. Phys. A 602, 388-404 (1996) doi:10.1016/0375-9474(96)00101-7
  * [85] J. Carbonell and V. A. Karmanov, Eur. Phys. J. A 6, 9-19 (1999) doi:10.1007/s100500050311 [arXiv:nucl-th/9902053 [nucl-th]].
  * [86] J. Carbonell and V. A. Karmanov, Nucl. Phys. A 663, 361-364 (2000) doi:10.1016/S0375-9474(99)00617-X
  * [87] V. A. Karmanov and J. Carbonell, [arXiv:nucl-th/0207075 [nucl-th]].
  * [88] S. J. Brodsky, J. R. Hiller, D. S. Hwang and V. A. Karmanov, Phys. Rev. D 69, 076001 (2004) doi:10.1103/PhysRevD.69.076001 [arXiv:hep-ph/0311218 [hep-ph]].
  * [89] A. Freese and W. Cosyn, Phys. Rev. D 106, no.11, 114013 (2022) doi:10.1103/PhysRevD.106.114013 [arXiv:2207.10787 [hep-ph]].
  * [90] P. A. M. Dirac, Rev. Mod. Phys. 21, 392-399 (1949) doi:10.1103/RevModPhys.21.392
  * [91] S. J. Brodsky, D. S. Hwang, B. Q. Ma and I. Schmidt, Nucl. Phys. B 593, 311-335 (2001) doi:10.1016/S0550-3213(00)00626-X [arXiv:hep-th/0003082 [hep-th]].
  * [92] S. Rodini, A. Metz and B. Pasquini, JHEP 09, 067 (2020) doi:10.1007/JHEP09(2020)067 [arXiv:2004.03704 [hep-ph]].
  * [93] B. Pasquini and S. Boffi, Phys. Lett. B 653, 23-28 (2007) doi:10.1016/j.physletb.2007.07.037 [arXiv:0705.4345 [hep-ph]].
  * [94] A. Amor-Quiroz, W. Focillon, C. Lorcé and S. Rodini, Eur. Phys. J. C 83, no.11, 1012 (2023) doi:10.1140/epjc/s10052-023-12190-7 [arXiv:2304.10339 [hep-ph]].


#  $\mathcal{PT}-$symmetry and chaos control via dissipative optomechanical
coupling

S. R. Mbokop Tchounda<EMAIL_ADDRESS>Department of
Physics, Faculty of Science, University of Yaounde I, P.O. Box 812, Yaounde,
Cameroon  P. Djorwé<EMAIL_ADDRESS>Department of Physics, Faculty of
Science, University of Ngaoundere, P.O. Box 454, Ngaoundere, Cameroon
Stellenbosch Institute for Advanced Study (STIAS), Wallenberg Research Centre
at Stellenbosch University, Stellenbosch 7600, South Africa  M. V. Tchakui
<EMAIL_ADDRESS>Department of Electrical and Electronic Engineering,
National Higher Polytechnic Institute, University of Bamenda, PO Box 39,
Bamenda, Cameroon  S. G. Nana Engo<EMAIL_ADDRESS>Department of Physics, Faculty of Science, University of Yaounde I, P.O. Box
812, Yaounde, Cameroon

###### Abstract

We study a dissipative mechanically coupled optomechanical system that hosts
gain and loss. The gain (loss) is engineered by driving a purely dispersive
optomechanical cavity with a blue-detuned (red-detuned) electromagnetic field.
By taking into account the dissipative coupling, the Exceptional Point (EP),
which is the $\mathcal{PT}-$symmetry phase transition, occurs at low threshold
driving strength compared to what happens in a solely dispersive system. In
the linear regime, the $\mathcal{PT}-$symmetry is unbroken and the dissipative
term induces strong coupling between the mechanical resonators, leading to an
increase of energy exchange. For strong enough driving, the system enters into
a nonlinear regime where the $\mathcal{PT}-$symmetry is broken. In this
regime, the mechanical resonators exhibit chaotic beats like-behaviour in the
purely dispersive system. By switching on the dissipative coupling, the
complex dynamics is switched off, and this restores regular dynamics into the
system. This work suggests a way to probe quantum phenomena in dissipative
$\mathcal{PT}-$symmetric systems at low-threshold driving strength. It also
provides a new control scheme of complex dynamics in optomechanics and related
fields.

Optomechanics, exceptional point, frequency locking, chaos

###### pacs:

42.50.Wk, 42.50.Lc, 05.45.Xt, 05.45.Gg

##  I Introduction

Optomechanical coupling, an interaction between an electromagnetic field and
mechanical objects, is a nice platform to foster quantum breakthroughs for
both fundamental and technological purposes Aspelmeyer et al. (2014). Two main
optomechanical coupling are investigated in the literature, the dispersive
couplings Aspelmeyer et al. (2014) and the dissipative/reactive coupling Li et
al. (2009); Huang and Agarwal (2010); Fu et al. (2012); Weiss et al. (2013);
ju Gu et al. (2013); Qu and Agarwal (2015); Wu et al. (2014); Tagantsev et al.
(2018). The former is well known and widely used, while the latter, which also
leads to nice features as well seems to attract less attention. The process
where the cavity frequency is modulated by the mechanical vibrations leads to
the dispersive coupling, while the dissipative coupling results from the
modulation of both the resonance frequency and the linewidth of the cavity
mode by the mechanical vibrations. In the dispersive case, optomechanical
cavities (COM) have led to several investigations including quantum ground
state cooling Clark et al. (2017); Djorwé et al. (2012), quantum entanglement
Kotler et al. (2021); Tchodimou et al. (2017); Djorwé et al. (2014), squeezing
Wollman et al. (2015); Djorwé et al. (2013a) , nonlinear dynamics Djorwé et
al. (2013a); Alphonse et al. (2022); Djorwe et al. (2019a); Djorwé et al.
(2020) and chaos Monifi et al. (2016); Djorwe et al. (2018); Kingni et al.
(2020). In terms of array systems, dispersive COM have been used to engineer
Parity-Time ($\mathcal{PT}$) symmetry and exceptional points (EPs), which are
Non-Hermitian degeneracies. These EPs have triggered interesting and
intriguing effects in physics EPs Miri and Alù (2019). Among them are stopping
light Goldzak et al. (2018), loss-induced suppression and revival of lasing,
pump-induced lasing death, unidirectional invisibility Peng et al. (2014),
collective phenomena Djorwé et al. (2020); Djorwe et al. (2018) , sensors
Djorwe et al. (2019b); Chen et al. (2017); Hodaei et al. (2017) and
topological features Xu et al. (2016); Ren et al. (2022); Peano et al. (2015).
Similarly, interesting achievements have been made in the dissipative’s COM
framework. These achievements include normal mode splitting Huang and Agarwal
(2010), strong-coupling effect Weiss et al. (2013), squeezing ju Gu et al.
(2013); Qu and Agarwal (2015); Tagantsev et al. (2018) and detection Wu et al.
(2014). To the best of our knowledge, there are no $\mathcal{PT}-$symmetry nor
EPs studies in the literature involving dissipative COM, although they share
common promising features with their dispersive counterparts. Addressing and
bridging this gap may provide new opportunities for optomechanical platforms
to uncover interesting features.

Our work aims to investigate the effect of dissipative coupling on
$\mathcal{PT}-$ symmetry and EPs engineering in coupled optomechanical
cavities. The benchmark system consists of two optomechanical systems where
the involved mechanical resonators are mechanically coupled. The gain (loss)
in the system is controlled by driving the resonator with a blue (red) detuned
laser. More interestingly, the dissipative coupling rate in our proposal has
been revealed as an additional parameter to adjust the gain (loss) rate even
at a fixed driving field. This flexibility to control gain and loss is a good
prerequisite for $\mathcal{PT}-$symmetry engineering, which is an advantage
over the dispersive COM structures. In the linearized regime, we carried out
the linear stability study and confirmed that this stability fails at the EP
when the system reaches the phonon lasing or the parametric gain threshold
Foulla et al. (2017). For purely dissipative coupling, we found that the EP is
shifted towards low values of the driving field as the dissipation rate
increases. This leads to strong coupling and energy exchange between the
hybridized supermodes in the system, even for weak driving strength. In the
nonlinear regime, the dissipative coupling quenches the complex dynamics by
switching them to regular behavior, providing a chaos control scheme. Due to
these interesting features, the strong coupling effect induced here can be
deeply explored for quantum investigations, while the chaos control may offer
great opportunities for collective phenomena Djorwé et al. (2020) and
dynamical entanglement Wang et al. (2014).

This work is organized as follows. The section II describes the model and its
dynamical equations up to the EP features. The section III is devoted to the
dissipatively induced strong coupling while the chaos control through
dissipative coupling is presented in section IV. The work is concluded in
section V.

##  II Modelling and dynamical equations

###  II.1 Modelling

Our benchmark system consists of two optomechanical cavities in which the
mechanical resonators are mechanically coupled as sketched in FIG. 1. One of
the COM is driven by the red-detuned laser to engineer losses in the cavity,
while the other is driven by a blue-detuned field to provide gain to the
cavity with that offsets the losses. This system is assumed to contain both
dispersive and dissipative couplings.

Figure 1: Optomechanical cavities containing both dispersive and dissipative
couplings. The mechanical resonators are mechanically coupled.

In the rotating frame of the driving fields ($\omega_{p}^{j}$ with $j=1,2$),
the Hamiltonian ($\hslash=1$) describing this system is

$H=H_{\rm OM}+H_{\rm int}+H_{\rm drive}+H_{\kappa}+H_{\gamma_{m}},$ (1)

where

$\displaystyle H_{\rm OM}$
$\displaystyle=\sum_{j=1,2}\left[-\Delta_{j}a^{{\dagger}}_{j}a_{j}+\omega_{j}b^{{\dagger}}_{j}b_{j}+\frac{L}{c}\bar{n}_{g}\omega_{p}^{j}(a^{in}_{j})^{2}\right],$
(2) $\displaystyle H_{\rm int}$
$\displaystyle=-\sum_{j=1,2}\left[g_{j}a^{{\dagger}}_{j}a_{j}(b^{{\dagger}}_{j}+b_{j})\right]-J_{m}(b^{{\dagger}}_{1}b_{2}+b_{1}b^{{\dagger}}_{2}),$
(3) $\displaystyle H_{drive}$ $\displaystyle=i\sqrt{\kappa_{x}}\alpha^{\rm
in}(a^{{\dagger}}_{j}-a_{j}).$ (4)

where $H_{\rm OM}$, $H_{\rm int}$ and $H_{\rm drive}$ are the optomechanical,
coupling interactions and the driving Hamiltonians respectively. The cavity
and the mechanical decay are captured by the Hamiltonians $H_{\kappa}$ and
$H_{\gamma_{m}}$. The annihilation bosonic field operators describing the
optical and mechanical resonators are $a_{j}$ and $b_{j}$. The mechanical
displacements $x_{j}$ are related to the operators $b_{j}$ as $x_{j}=x_{\rm
ZPF}(b_{j}+b^{{\dagger}}_{j})$, where $x_{\rm
ZPF}=\sqrt{\frac{\hslash}{2m\omega_{j}}}$ is the zero-point fluctuation
amplitude of the mechanical resonator, where $\omega_{j}$ being the mechanical
frequency of the $j^{th}$ resonator as shown in FIG. 1. The frequency
detunings between the driving ($\omega_{p}^{j}$) and the cavity
($\omega_{cav}^{j}$) are defined by
$\Delta_{j}=\omega_{p}^{j}-\omega_{cav}^{j}$, and the phonon-hopping rate
between the two mechanical resonators is $J_{m}$. $L$ is the length of the
free-standing waveguide interacting with the optical cavity, $c$ the speed of
light in vacuum, and $\bar{n}_{g}$ the group index of the waveguide optical
mode. The amplitude of the driving fields is $\alpha^{in}$, which is assumed
to be the same for both cavities, and is related to the input power $P_{in}$
as $a^{in}=\sqrt{\frac{P_{in}}{\hbar\omega_{p}}}$. The dispersive
optomechanical rates are $g_{j}$ while the dissipative coupling $g_{\kappa}$
results from the modulation of the cavity dissipation $\kappa_{x}$ by the
mechanical displacement as,

$\kappa_{x}\approx\kappa+\kappa_{\rm
om}x=\kappa+g_{\kappa}(b^{{\dagger}}+b)=\kappa(1+\eta(b^{{\dagger}}+b)),$ (5)

where $g_{\kappa}=\kappa_{\rm om}x_{\rm ZPF}$ and
$\eta=\frac{g_{\kappa}}{\kappa}$.

###  II.2 Dynamical equations

From the Heisenberg’s equation
$\dot{\mathcal{O}}=i[H,\mathcal{O}]+\mathcal{N}$, where
$\mathcal{O}\equiv(a_{j},b_{j})$ and
$\mathcal{N}\equiv(a^{in}_{j},b^{in}_{j})$ the related noise operator, the
quantum Langevin equations (QLEs) of our system read,

$\displaystyle\begin{split}\dot{a_{j}}&=\left(i\left(\Delta_{j}+(g_{j}+i\frac{\eta_{j}\kappa}{2})(b^{{\dagger}}_{j}+b_{j})\right)-\frac{\kappa}{2}\right)a_{j}\\\
&+\sqrt{\kappa}\left(1+\frac{\eta_{j}}{2}(b^{{\dagger}}_{j}+b_{j})\right)\alpha^{in}+\sqrt{\kappa}a^{in},\end{split}$
(6a)
$\displaystyle\begin{split}\dot{b_{j}}&=-(i\omega_{j}+\frac{\gamma_{m}}{2})b_{j}+iJ_{m}b_{3-j}\\\
&+\eta_{j}\sqrt{\kappa}(a^{{\dagger}}_{j}-a_{j})\alpha^{in}_{j}+ig_{j}a^{{\dagger}}a+\sqrt{\gamma_{m}}b^{in}_{j}.\end{split}$
(6b)

These noise operators $\mathcal{N}\equiv(a^{in}_{j},b^{in}_{j})$ have zero
mean and are characterized by the following autocorrelation functions
Tchodimou et al. (2017),

$\displaystyle\langle\mathcal{N}(t)\mathcal{N}^{{\dagger}}(t^{\prime})\rangle$
$\displaystyle=(n_{\nu}+1)\delta(t-t^{\prime}),$ (7)
$\displaystyle\langle\mathcal{N}^{{\dagger}}(t)\mathcal{N}(t^{\prime})\rangle$
$\displaystyle=n_{\nu}\delta(t-t^{\prime}),$ (8)

with $n_{\nu}\equiv(n_{th},n_{a})$, where
$n_{th}=\rm{exp(\tfrac{\hslash\omega_{j}}{k_{B}T}-1)^{-1}}$ and $n_{a}=\langle
a^{in{\dagger}}_{j}a^{in}_{j}\rangle$. For a seek of simplicity and without
loss of generality, we assume in the following that
$\omega_{j}\equiv\omega_{m}$, $g_{j}\equiv g_{m}$ and $\eta_{j}\equiv\eta$.

To gain insight into our system and to understand its dynamics, we linearize
the set of nonlinear equations given in Eq. (6), which consists of splitting
the field operators as $\mathcal{O}=\langle O\rangle+\delta\mathcal{O}$, where
$\langle O\rangle\equiv\beta_{j},\alpha_{j}$ are the coherent complex parts of
the operator and $\delta\mathcal{O}\equiv\delta\beta_{j},\delta\alpha_{j}$
their associated fluctuations. Therefore, the coherent dynamics reads

$\displaystyle\begin{split}&\dot{\alpha_{j}}=(i\tilde{\Delta_{j}}-\frac{\kappa}{2})\alpha_{j}+\sqrt{\kappa}(1+\eta\Re{(\beta_{j})})\alpha^{in},\\\
&\begin{aligned}
\dot{\beta_{j}}=&-(i\omega_{m}+\frac{\gamma_{m}}{2})\beta_{j}+iJ_{m}\beta_{3-j}\\\
&+\eta\sqrt{\kappa}(\alpha^{{\dagger}}_{j}-\alpha_{j})\alpha^{in}+ig\left|\alpha_{j}\right|^{2},\end{aligned}\end{split}$
(9)

and the fluctuation dynamics yields,

$\begin{split}\dot{\delta\alpha_{j}}&=(i\tilde{\Delta_{j}}-\frac{\kappa}{2})\delta\alpha_{j}+\sqrt{\kappa}\left(1+\eta\Re{(\beta_{j})}\right)\delta\alpha^{in}\\\
&+\left(i\alpha_{j}(g+i\frac{\eta\kappa}{2})+\frac{\eta\sqrt{\eta}\alpha^{in}}{2}\right)(\delta\beta_{j}+\delta\beta^{in}_{j}),\\\
\dot{\delta\beta_{j}}&=-(i\omega_{m}+\frac{\gamma_{m}}{2})\delta\beta_{j}+iJ_{m}\delta\beta_{3-j}\\\
&+\eta\sqrt{\kappa}(\delta\alpha^{{\dagger}}_{j}-\delta\alpha_{j})\alpha^{in}+\eta\sqrt{\kappa}(\alpha^{{\dagger}}_{j}-\alpha_{j})\delta\alpha^{in}\\\
&+ig(\alpha^{{\dagger}}_{j}\delta\alpha_{j}+\alpha_{j}\delta\alpha^{{\dagger}}_{j}),\end{split}$
(10)

with
$\tilde{\Delta_{j}}=\Delta_{j}+2\Re{(\beta_{j})}(g+i\frac{\eta\kappa}{2})$.

In the following we will work out the steady state equation of our system from
Eq. (9) and use it to study both the stability and to derive the EP. The
steady state condition is that the averaged variables in Eq. (9) are no longer
time dependent. This means that $\dot{\alpha_{j}}=\dot{\beta_{j}}=0$, which
leads to,

$\begin{split}&0=(i\tilde{\Delta_{j}}-\frac{\kappa}{2})\bar{\alpha_{j}}+\sqrt{\kappa}(1+\eta\Re{(\bar{\beta_{j})}})\alpha^{in},\\\
&\begin{split}0&=(i\omega_{m}+\frac{\gamma_{m}}{2})\bar{\beta_{j}}-iJ_{m}\bar{\beta}_{3-j}\\\
&-\eta\sqrt{\kappa}(\bar{\alpha}^{{\dagger}}_{j}-\bar{\alpha}_{j})\alpha^{in}-ig\left|\bar{\alpha}_{j}\right|^{2}.\end{split}\end{split}$
(11)

This set of equations can be solved to obtain the steady state equation for
either $\bar{\alpha_{j}}$ or $\bar{\beta_{j}}$. After tedious computation, and
because the expressions are cumbersome, we present only the equation for
$\bar{\alpha_{j}}$ which reads,

$\left|\bar{\alpha_{j}}\right|^{6}+a_{0}\left|\bar{\alpha_{j}}\right|^{5}+a_{1}\left|\bar{\alpha_{j}}\right|^{4}+a_{2}\left|\bar{\alpha_{j}}\right|^{3}+a_{3}\left|\bar{\alpha_{j}}\right|^{2}+a_{4}\left|\bar{\alpha_{j}}\right|+a_{5}=0,$
(12)

where the involved coefficients are, $a_{0}=\frac{b_{0}}{c}$,
$a_{1}=\frac{b_{1}}{c}$, $a_{2}=\frac{b_{2}}{c}$, $a_{3}=\frac{b_{3}}{c}$,
$a_{4}=\frac{b_{4}}{c}$ and $a_{5}=\frac{b_{5}}{c}$ with

$\displaystyle\begin{split}&b_{0}=-4g\Omega^{2}\eta\sqrt{\kappa}\alpha^{in}(\eta^{2}\kappa^{2}+4),\\\
\hfil\displaystyle\begin{split}&b_{1}=4\Omega^{2}\eta^{2}\kappa\alpha^{in}(\eta^{2}\kappa^{2}+4)-\kappa(g\Omega\eta\alpha^{in})^{2}\\\
&+g\Omega(4\Delta_{j}+\eta\kappa^{2}),\end{split}\\\
&b_{2}=4\Omega^{2}g(\eta\alpha^{in})^{3}\kappa\sqrt{\kappa}-2\Omega\sqrt{\kappa}\alpha^{in}(4\Delta_{j}+\eta\kappa^{2}),\\\
&b_{3}=2\Omega\eta\kappa
g(\alpha^{in})^{2}-(2\Omega\kappa)^{2}(\eta\alpha^{in})^{4}+\Delta_{j}^{2}+\frac{\kappa^{2}}{4},\\\
&b_{4}=-4\eta\kappa(\alpha^{in})^{3},\\\
&b_{5}=-\kappa\alpha^{in},\end{split}$ (13)

and the other of the parameters are, $c=(g\Omega)^{2}((\eta\kappa)^{2}+4)$ and
$\Omega=\tfrac{(\omega_{m}^{2}-J_{m}^{2})(\omega_{m}+J_{m})}{(J_{m}^{2}-\omega_{m}^{2})^{2}+(\omega_{m}\gamma_{m})}$.
The solution of Eq. (12) can then be used in the fluctuation dynamics Eq. (10)
for the EP analysis. For this purpose we change to another interaction picture
by introducing the following slowly varying operators with tildes, $\delta
a=\delta\tilde{a}e^{i\tilde{\Delta}t}$ and $\delta
b_{j}=\delta\tilde{b_{j}}e^{-i\omega_{j}t}$. In the limit of
$\omega_{j}\gg(G_{j},\kappa)$ the rotating wave approximation can be invoked
depending on the blue or red sideband resonances, where the fast oscillating
terms are neglected. Moreover, under the condition $\kappa\gg(G,\gamma_{j})$,
which is satisfied in our analysis, the cavity field can be eliminated
adiabatically, and this leads us to the effective dynamical system for the two
mechanical resonators Jiang et al. (2021):

$\begin{split}\dot{\delta\beta_{1}}&=(\frac{\gamma_{eff}^{1}}{2}-i\omega_{eff}^{1})\delta\beta_{1}+iJ_{m}\delta\beta_{2}+\sqrt{\gamma_{m}}\delta\beta^{in}_{1}\\\
&+(i\sqrt{\Gamma}+2\eta\alpha^{in})(1+\eta\Re{(\bar{\beta})})\delta\alpha^{in{\dagger}},\\\
\dot{\delta\beta_{2}}&=(\frac{\gamma_{eff}^{2}}{2}-i\omega_{eff}^{2})\delta\beta_{2}+iJ_{m}\delta\beta_{1}+\sqrt{\gamma_{m}}\delta\beta^{in}_{2}\\\
&+(i\sqrt{\Gamma}-2\eta\alpha^{in})(1+\eta\Re{(\bar{\beta})})\delta\alpha^{in}.\\\
\end{split}$ (14)

$\Gamma=\tfrac{4|G|^{2}}{\kappa}$ is the optical damping with
$G=g\bar{\alpha}$ the optomechanical induced coupling rate. We have assumed
that the two cavities are driven with the same driving strength
$\alpha^{in}_{j}\equiv\alpha^{in}$ so that
$\bar{\alpha_{j}}\equiv\bar{\alpha}$ and $\bar{\beta_{j}}\equiv\bar{\beta}$.
We have also defined both the effective dissipation and the frequency as,

$\displaystyle\begin{split}\gamma_{eff}^{1}&=\Gamma-\gamma_{m}-2\eta^{2}\alpha^{in}(\alpha^{\ast}\sqrt{\kappa}-\alpha^{in}),\\\
\gamma_{eff}^{2}&=-(\Gamma+\gamma_{m})+2\eta^{2}\alpha^{in}(\alpha\sqrt{\kappa}-\alpha^{in}),\\\
\omega_{eff}^{1}&=\omega_{m}+\frac{\eta\sqrt{\Gamma}}{2}(\alpha^{\ast}\sqrt{\kappa}+\alpha^{in}),\\\
\omega_{eff}^{2}&=\omega_{m}+\frac{\eta\sqrt{\Gamma}}{2}(\alpha\sqrt{\kappa}+\alpha^{in}).\end{split}$
(15)

###  II.3 EP analysis

Figure 2: (a) and (b) are the photon number and the mechanical displacement
obtained from numerical simulation of the full nonlinear equation in Eq. (6),
where the noise terms are dropped and the operators are averaged, i.e.,
$\langle a_{j}\rangle\equiv\alpha_{j}$, and $\langle
b_{j}\rangle\equiv\beta_{j}$. (c) and (d) are the eigenfrequencies and
eigendampings of our system. It can be seen that they both coalesce at the EP,
which agrees well with the parametric instability threshold (see (a) and (b))
where the system lost its stability. The parameters used are
$\Delta_{1}=\omega_{m}$, $\Delta_{2}=-\omega_{m}$, $g_{m}=1.076\times
10^{-4}\omega_{m}$, $\kappa=7.3\times 10^{-2}\omega_{m}$, $\eta=0.1g_{m}$,
$\gamma_{m}=1.076\times 10^{-5}\omega_{m}$, and $J_{m}=4\times
10^{-4}\omega_{m}$.

In order to process with the EP analysis, drop the noise term in Eq. (14) and
write the noiseless effective equations in the compact form,

$i\partial_{t}\psi=\rm{H_{eff}}\psi,$ (16)

with $\psi=(\delta\beta_{1},\delta\beta_{2})$ and the effective Hamiltonian
given by,

$\rm{H_{eff}}=\begin{pmatrix}\omega_{eff}^{1}+i\frac{\gamma_{eff}^{1}}{2}&-J_{m}\\\
-J_{m}&\omega_{eff}^{2}+i\frac{\gamma_{eff}^{2}}{2}\end{pmatrix}.$ (17)

The eigenvalues of the Hamiltonian given in Eq. (17) are,

$\lambda_{\pm}=\frac{1}{2}(\omega_{eff}^{1}+\omega_{eff}^{2})+\frac{i}{4}(\gamma_{eff}^{1}+\gamma_{eff}^{2})\pm\frac{\sigma}{4},$
(18)

where,

$\sigma=\sqrt{(2\Delta\omega_{eff}+i\Delta\gamma_{eff})^{2}+16J_{m}^{2}},$
(19)

with $\Delta\omega_{eff}=\omega_{eff}^{1}-\omega_{eff}^{2}$ and
$\Delta\gamma=\gamma_{eff}^{1}-\gamma_{eff}^{2}$. From these eigenvalues, the
eigenfrequencies and eigendampings of the system are defined as the real
($\omega_{\pm}=\Re{(\lambda_{\pm})}$) and imaginary
($\gamma_{\pm}=\Im{(\lambda_{\pm})}$) parts of $\lambda_{\pm}$ respectively.
The EP occurs when these two pairs of frequencies and dampings coalesce, i.e.
$\omega_{-}=\omega_{+}$ and $\gamma_{-}=\gamma_{+}$, which means that
$\sigma=0$. This condition on $\sigma$ can be fulfilled in our proposal by
tuning the driving field instead of $\eta$, which is not easily tunable as a
driving field. Fig. 2 shows the agreement between the numerical stability and
the analytical investigation by the EP. In fact, Fig. 2(a-b) are the numerical
simulations of the photon number and the mechanical displacement from Eq. (6),
where the noise terms are dropped and the operators are averaged, i.e.,
$\langle a_{j}\rangle\equiv\alpha_{j}$ and $\langle
b_{j}\rangle\equiv\beta_{j}$. These bosonic quantities are extracted from the
steady state time series for different values of the driving strength. The
photon number is given by $\langle\alpha^{{\dagger}}_{j}\alpha_{j}\rangle$
while the mechanical position is quantified by
$x_{j}=\beta_{j}+\beta_{j}^{{\dagger}}$. It can be seen that these quantities
are exponentially amplified around $\alpha^{in}\sim 130\omega^{1/2}_{m}$,
known as the parametric instability threshold where the system loses its
stability Foulla et al. (2017). Furthermore, Figs. 2(c-d) show the
eigenfrequencies and eigendampings of our system, which both coalesce at the
EP, and this also happens around the threshold mentioned in Figs. 2(a-b). This
agreement confirms the validity of our analytical investigation. From this
analysis, we devote the next section to study the effect of the dissipative
coupling $\eta$ on the dynamical behavior of our system.

##  III Dissipation induces strong coupling

Figure 3: (a) and (b) are the 3D representation of the eigenfrequencies in
purely dispersive ($\eta=0$) and with a fixed dissipative coupling rate at
$\eta=g_{m}$. (c) and (d) highlight the EP evolution, which is depicted with
black line in (a) and (b). The other parameters are the same as in Fig. 2.

Due to the relationship between the effective parameters and the dissipative
coupling in Eq. (15), it is obvious that the gain-loss engineering will depend
on $\eta$, and this will induce a modification on the EP behavior in our
system. Fig. 3a shows the structure of the eigenfrequencies for the purely
dispersive case ($\eta=0$), while Fig. 3b shows the eigenfrequencies when the
dissipative coupling is fixed at $\eta=g_{m}$. It can be seen that in the
purely dispersive case the EPs appear for relatively high values of the
driving strength compared to the case where the dissipative coupling is taken
into account. This observation is highlighted in Fig. 3c and Fig. 3d, where
the 2D projections of Fig. 3a and Fig. 3b are extracted. Indeed, Fig. 3d shows
that for strong dispersive coupling, EP can be engineered for low values of
the driving field when $\eta$ is involved. It follows that dissipative
coupling induces a low power $\mathcal{PT}-$symmetry phase transition compared
to the purely dispersive coupling case. From an experimental point of view,
this may be a crucial aspect to consider and a good motivation for more
dissipative $\mathcal{PT}-$symmetry features. To confirm the above analysis,
we have further performed numerical simulations in both the linear and
nonlinear regime. From Fig. 2, it is clear that the linear regime corresponds
to the unbroken $\mathcal{PT}-$symmetry regime, where the system exhibits Rabi
oscillation-like dynamics due to the two frequencies involved, as shown in
Fig. 2c. In the nonlinear regime, however, Fig. 2 shows that the system is in
the broken $\mathcal{PT}-$symmetry phase, where the mechanical resonators
dissipate differently, as shown in Fig. 2d.

Our numerical simulations for the linear regime are displayed in Fig. 4.
Starting from $\eta=0$ in Fig. 4a for a fixed value of the driving strength
$\alpha^{in}=20[\omega_{m}^{1/2}]$, one can see that the system is driven
towards a strong coupling regime by tuning the value of $\eta$ (see Fig.
4(b-d)). In fact, Fig. 4a shows almost no Rabi oscillations, while these
oscillations are clearly visible in Fig. 4d, clearly showing the different
frequencies in the system. This strong coupling is reminiscent of a strong
energy exchange between the two mechanical resonators involved. It is
concluded that the dissipative coupling rate $\eta$ can be used to enhance PT-
symmetry features such as strong coupling Weiss et al. (2013); Huang and
Agarwal (2010) and to stimulate efficient energy transfer Xu et al. (2016).

##  IV Dissipative control of chaos

This section describes the dynamical behavior of our system above the
nonlinear regime, i.e., above the parametric threshold shown in Fig. 2. Above
this threshold and in the purely dispersive regime, the mechanical resonators
oscillate with a common frequency while each of them dissipates differently
according to the eigenvalues depicted in Fig. 2(c-d). In fact, these
mechanical resonators will first settle into a limit cycle regime (or in
regular oscillations) and they may exhibit collective phenomena Djorwé et al.
(2013b); Djorwé et al. (2020). Once that the driving force is strong enough,
multi-frequency and chaotic behavior can be induced Djorwe et al. (2018).

Figure 4: Dynamical evolution of the mechanical resonators in the unbroken
$\mathcal{PT}-$symmetric regime. (a)-(d) are the mechanical displacements for
different values of the dissipative coupling rate, i.e.
$\eta=0;0.1;0.5;g_{m}$. The driving force is fixed at
$\alpha^{in}=20[\omega_{m}^{1/2}]$, the mechanical dissipation at
$\gamma_{m}=1.076\times 10^{-5}\omega_{m}$ and the phonon hopping rate at
$J_{m}=4\times 10^{-4}\omega_{m}$. The other parameters are as in Fig. 2.

To illustrate the effect of the dissipative coupling in our proposal, we first
provide a general view of the dynamical behavior through a bifurcation diagram
as shown in Fig. 5, where the dynamics involved will be explained later. This
bifurcation diagram shows how our system exhibits chaotic dynamics for zero
and weak values of the dissipative coupling ($\eta\in[0-0.025]\omega_{m}$). As
the dissipative coupling increases, the system enters a quasi-periodic-like
behavior ($\eta\in[0.026-0.111]\omega_{m}$), where regular dynamics are
gradually recovered. For sufficiently strong $\eta$, the chaotic dynamics is
completely switched off and the regular dynamics is restored
($\eta>0.111\omega_{m}$). These three dynamical states are denoted as zones I,
II and III in Fig. 5, and are delineated by the dashed vertical lines to guide
the eye.

To illustrate the effect of dissipative coupling in our proposal, we first
provide a general view of the dynamical behavior through a bifurcation diagram
as shown in Fig. 5, where the dynamics involved will be explained later. This
bifurcation diagram shows how our system exhibits chaotic dynamics for weak
values of the dissipative coupling ($\eta\in[0-0.025]\omega_{m}$). As the
dissipative coupling increases, the system enters a quasi-periodic-like
behavior ($\eta\in[0.026-0.111]\omega_{m}$), where regular dynamics are
gradually recovered. For sufficiently large $\eta$, the chaotic dynamics is
completely switched off and the regular dynamics is restored
($\eta>0.111\omega_{m}$). These three dynamical states are denoted as zones I,
II and III in Fig. 5, and are delimited by the dotted vertical lines to guide
the eye.

Figure 5: Bifurcation diagram showing the vibrational states of the mechanical
resonators as a function of the dissipative coupling strength $\eta$. The
driving strength is fixed at $\alpha^{in}=10^{4}[\omega_{m}^{1/2}]$, and the
other parameters are as in Fig. 2.

To further highlight these dynamics, we have plotted the Poincaré section for
certain values of $\eta$ in Fig. 6. The dissipative coupling values used are
$\eta=0$, $\tfrac{\eta}{\omega_{m}}=2\times 10^{-2}$,
$\tfrac{\eta}{\omega_{m}}=9.7\times 10^{-2}$, and
$\tfrac{\eta}{\omega_{m}}=15.3$, respectively. One can see that these Poincaré
sections agree well with the dynamical state carried out in the bifurcation
diagram. In fact, Fig. 6(a,b) show a chaotic behavior in our system revealed
by an abundance of random points, Fig. 6c shows a quasi-periodic state, while
Fig. 6d reveals a regular dynamics in our system.

Figure 6: Poincaré section highlighting certain vibrational states based on
the bifurcation diagram shown in Fig. 5. The strengths of the dissipative
coupling are $\tfrac{\eta}{g_{m}}=0;2\times 10^{-2};9.7\times 10^{-2}$; and
$\tfrac{\eta}{g_{m}}=15.3\times 10^{-2}$, respectively. The driving force is
fixed at $\alpha^{in}=10^{4}[\omega_{m}^{1/2}]$, and the other parameters are
as in Fig. 2.

For more details about these dynamical states, we have plotted their
stationary time evolution in Fig. 7. It can be seen that Fig. 7(a-b) exhibit
complex dynamics as predicted both in the bifurcation diagram and in the
Poincaré sections shown in Fig. 5 and Fig. 6, respectively. Such complex
dynamics, resulting from a distortion of the oscillation beats induced by a
strong coupling in our system, are referred to as chaotic beats Ahamed et al.
(2011). Furthermore, the regular dynamical states carried out in Fig. 7(c-d)
are also in agreement with the above predictions. These analyses show how
dissipative coupling can be used to switch from a chaotic to a regular state.
This dynamical control shows how $\eta$ can be used to suppress unwanted
behavior in dynamical systems. Moreover, the oscillation beats shown in Fig.
7c suggest how the dissipative coupling tends to restore
$\mathcal{PT}-$symmetry in the system, which is similar to the investigation
carried out in Fig. 4.

Figure 7: Dynamical evolution of the mechanical resonators in the broken
$\mathcal{PT}-$symmetric regime. (a) and (b) show complex dynamics, and these
complex dynamics are swapped out in (c) and (d), where regular dynamics are
gradually restored. From (a) to (d), the dissipative coupling strengths
correspond to those in Fig. 6(a) to Fig. 6(d), respectively. The other
parameters are the same as in Fig. 2.

##  V Conclusion

We have investigated the effect of the dissipative coupling on a mechanically
coupled optomechanical system. In our proposal, one cavity is driven with
blue-detuned field to induce gain, while the other is driven with red-detuned
field to induce losses. These gain and loss rates are shown to depend on the
dissipative coupling parameter. We have shown that EP occurs at low driving
strengths when dissipative coupling is turned on, compared to a purely
dispersive system. In the linear regime, where the $\mathcal{PT}-$ symmetry is
unbroken, the dissipative coupling induces a strong coupling between the
mechanical resonators, leading to an enhancement of the energy exchange
between them. Above the parametric instability threshold, which occurs at the
EP, the $\mathcal{PT}-$symmetry is broken and the system enters into a
nonlinear regime. In this broken $\mathcal{PT}-$symmetry regime, the purely
dispersive system exhibits complex dynamics represented by chaotic beats.
These complex dynamics are controlled by the dissipative coupling, which
switches the chaotic beats into regular beat oscillations. These results show
how dissipative $\mathcal{PT}-$symmetry systems can be used to harness new
features of EP with less driving power compared to their dispersive
counterparts. Furthermore, this work extends the EP opportunities to
dissipative optomechanical systems, and paves the way to study
$\mathcal{PT}-$symmetry features at low-threshold power.

## Acknowledgments

This work was carried out under the Iso-Lomso Fellowship at the Stellenbosch
Institute for Advanced Study (STIAS), Wallenberg Research Centre at
Stellenbosch University, Stellenbosch 7600, South Africa.

P. Djorwe and S.G. Nana Engo thank the Ministry of Higher Education of
Cameroon (MINESUP) for the financial support within the framework of the
"Research Modernization" grants.

## References

  * Aspelmeyer et al. (2014) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Reviews of Modern Physics 86, 1391 (2014).
  * Li et al. (2009) M. Li, W. H. P. Pernice, and H. X. Tang, Physical Review Letters 103, 223901 (2009).
  * Huang and Agarwal (2010) S. Huang and G. S. Agarwal, Physical Review A 81, 053810 (2010).
  * Fu et al. (2012) C.-B. Fu, K.-H. Gu, X.-B. Yan, X. Yang, C.-L. Cui, and J.-H. Wu, Physics Letters A 377, 133 (2012).
  * Weiss et al. (2013) T. Weiss, C. Bruder, and A. Nunnenkamp, New Journal of Physics 15, 045017 (2013).
  * ju Gu et al. (2013) W. ju Gu, G. xiang Li, and Y. ping Yang, Physical Review A 88, 013835 (2013).
  * Qu and Agarwal (2015) K. Qu and G. S. Agarwal, Physical Review A 91, 063815 (2015).
  * Wu et al. (2014) M. Wu, A. C. Hryciw, C. Healey, D. P. Lake, H. Jayakumar, M. R. Freeman, J. P. Davis, and P. E. Barclay, Physical Review X 4, 021052 (2014).
  * Tagantsev et al. (2018) A. K. Tagantsev, I. V. Sokolov, and E. S. Polzik, Physical Review A 97, 063820 (2018).
  * Clark et al. (2017) J. B. Clark, F. Lecocq, R. W. Simmonds, J. Aumentado, and J. D. Teufel, Nature 541, 191 (2017).
  * Djorwé et al. (2012) P. Djorwé, J. H. T. Mbé, S. G. N. Engo, and P. Woafo, Physical Review A 86, 043816 (2012).
  * Kotler et al. (2021) S. Kotler, G. A. Peterson, E. Shojaee, F. Lecocq, K. Cicak, A. Kwiatkowski, S. Geller, S. Glancy, E. Knill, R. W. Simmonds, et al., Science 372, 622 (2021).
  * Tchodimou et al. (2017) C. Tchodimou, P. Djorwe, and S. G. N. Engo, Physical Review A 96, 033856 (2017).
  * Djorwé et al. (2014) P. Djorwé, S. G. N. Engo, and P. Woafo, Physical Review A 90, 024303 (2014).
  * Wollman et al. (2015) E. E. Wollman, C. U. Lei, A. J. Weinstein, J. Suh, A. Kronwald, F. Marquardt, A. A. Clerk, and K. C. Schwab, Science 349, 952 (2015).
  * Djorwé et al. (2013a) P. Djorwé, S. N. Engo, J. T. Mbé, and P. Woafo, Physica B: Condensed Matter 422, 72 (2013a).
  * Alphonse et al. (2022) H. Alphonse, P. Djorwe, S. Abbagari, S. Y. Doka, and S. N. Engo, Chaos, Solitons and Fractals 154, 111593 (2022).
  * Djorwe et al. (2019a) P. Djorwe, Y. Pennec, and B. Djafari-Rouhani, Scientific Reports 9, 1684 (2019a).
  * Djorwé et al. (2020) P. Djorwé, Y. Pennec, and B. Djafari-Rouhani, Physical Review B 102, 155410 (2020).
  * Monifi et al. (2016) F. Monifi, J. Zhang, Ş. K. Özdemir, B. Peng, Y. xi Liu, F. Bo, F. Nori, and L. Yang, Nature Photonics 10, 399 (2016).
  * Djorwe et al. (2018) P. Djorwe, Y. Pennec, and B. Djafari-Rouhani, Physical Review E 98, 032201 (2018).
  * Kingni et al. (2020) S. T. Kingni, C. Tchodimou, D. P. Foulla, P. Djorwe, and S. G. N. Engo, The European Physical Journal Special Topics 229, 1117 (2020).
  * Miri and Alù (2019) M.-A. Miri and A. Alù, Science 363, eaar7709 (2019), ISSN 0036-8075.
  * Goldzak et al. (2018) T. Goldzak, A. A. Mailybaev, and N. Moiseyev, Physical Review Letters 120, 013901 (2018).
  * Peng et al. (2014) B. Peng, Ş. K. Özdemir, S. Rotter, H. Yilmaz, M. Liertzer, F. Monifi, C. M. Bender, F. Nori, and L. Yang, Science 346, 328 (2014).
  * Djorwe et al. (2019b) P. Djorwe, Y. Pennec, and B. Djafari-Rouhani, Physical Review Applied 12, 024002 (2019b).
  * Chen et al. (2017) W. Chen, Ş. K. Özdemir, G. Zhao, J. Wiersig, and L. Yang, Nature 548, 192 (2017).
  * Hodaei et al. (2017) H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Nature 548, 187 (2017).
  * Xu et al. (2016) H. Xu, D. Mason, L. Jiang, and J. G. E. Harris, Nature 537, 80 (2016).
  * Ren et al. (2022) H. Ren, T. Shah, H. Pfeifer, C. Brendel, V. Peano, F. Marquardt, and O. Painter, Nature Communications 13, 3476 (2022), ISSN 2041-1723.
  * Peano et al. (2015) V. Peano, C. Brendel, M. Schmidt, and F. Marquardt, Physical Review X 5, 031011 (2015).
  * Foulla et al. (2017) D. P. Foulla, P. Djorwé, S. T. Kingni, and S. G. N. Engo, Physical Review A 95, 013831 (2017).
  * Wang et al. (2014) G. Wang, L. Huang, Y.-C. Lai, and C. Grebogi, Physical Review Letters 112, 110406 (2014).
  * Jiang et al. (2021) C. Jiang, Y.-L. Liu, and M. A. Sillanpää, Physical Review A 104, 013502 (2021).
  * Djorwé et al. (2013b) P. Djorwé, J. H. T. Mbé, S. G. N. Engo, and P. Woafo, The European Physical Journal D 67, 45 (2013b).
  * Ahamed et al. (2011) A. I. Ahamed, K. Srinivansan, K. Murali, and M. Lakshmanan, International Journal of Bifurcation and Chaos 21, 737 (2011).