# Quantum simulation of Dirac fermion mode, Majorana fermion mode and
Majorana-Weyl fermion mode in cavity QED lattice

Sujit Sarkar Poornaprajna Institute of Scientific Research, 4 Sadashivanagar,
Bangalore 5600 80, India.
e-mail<EMAIL_ADDRESS>

###### Abstract

Quantum simulation aims to simulate a quantum system using a controble
laboratory system that underline the same mathematical model. Cavity QED
lattice system is that prescribe system to simulate the relativistic quantum
effect. We quantum simulate the Dirac fermion mode, Majorana fermion mode and
Majorana-Weyl fermion mode and a crossover between them in cavity QED lattice.
We also present the different analytical relations between the field operators
for different mode excitations.
PACS: 42.50.Pq, 03.65.Vf, 42.50.-p

Introduction: The difficulties in observing the real quantum relativistic
effects have generated immense interest in the quantum simulation physics. In
recent years, there has been increased interest in the simulation of
relativistic quantum effects using different physical system in which
parameters tunibility allows access to different physical regimes gerri ; seo
; kane ; bardyn ; bkumar ; trif ; sola1 ; sola2 ; xzhang ; otterbach ; clark .
These difficulties of observing quantum relativistic effects stimulate us to
study the quantum simulation physics of Dirac fermion mode, Majorana fermion
mode and Majorana-Weyl fermion modes in one dimensional cavity QED lattice and
a crossover from one mode of excitations to the another using the tunibility
of the physical parameters of the system.
In quantum simulation, one aim is to simulate a quantum system using a
controllable laboratory system that underlines the same mathematical models.
Therefore it is possible to simulate a quantum system that can be neither
efficiently simulated on a classical computer nor easily accessed
experimentally.
The recent experimental success in engineering strong interaction between the
photons and atoms in high quality micro-cavities opens up the possibility to
use the light matter system as quantum simulators for many body physics
[9-29]. Many interesting results are coming out to understand the complicated
quantum many body system. A focus on the coupled cavities is one of the most
potential candidate for an efficient quantum simulator due to the control of
the microcavities parameters and success of fabrication of large scale cavity
arrays [25-26].
The concept and existence of Majorana fermion mode is one of the most advance
research area in quantum condensed matter system. Majorana introduced a
special kind of fermions which are their own antiparticle, i.e., the neutral
particle majo ; wil . He introduced this particle to describe neutrions. In
recent years, there are several candidates of Majorana fermions in quantum
condensed matter system like quantum Hall system with filling fraction $5/2$
read ; read2 . Kitaev first found the existence of Majorana fermion mode in
one dimensional model kitaev . Many research groups have already proposed the
physically existence of MFs at the edge state of 1D system like electrostatic
defects lines in superconductor, quasi-one dimensional superconductor and cold
atom trapped in one dimension wimmer ; kop . Majorana fermions obey the non-
Abelian statics both in 2D and 1D, allowing of certain gate operation required
in quantum computation nayak .
In the present study one of our goal is to predict the presence of Dirac
fermion mode, Majorana fermion mode and Majoran-Weyl fermion mode and
crossover in our model system. In the Majorana-Weyl fermion mode, the fermion
mode satisfy the condition for particle and antiparticle condition but the
excitation is massless.
There are few studies to simulate the Majorana fermion modes in cavity QED
system but the model Hamiltonian and the associate relevant physics are
different and at the same time there is no coossover study between the Dirac
and Majorana fermion modes in the literature of cavity QED system [9-11].
The authors of Ref.3 have found the process of Dirac to Majorana fermion
converator on the surface of a 3D topological insulator. A Dirac fermion
injected by the voltage source is split into a pair of Majorana fermion and
then finally recombine before going to drain. But in our present study there
is no such split and fusion process of Dirac and Majorana fermion. In our
study the system is in Dirac fermion mode when only a atom-photon coupling or
a single Rabi frequency oscillation present in the system and the existence of
Majorana fermions when two atom photon couplings and the two laser frequencies
are simultaneously present in the system. For a quantum simulated Hamiltonian
for Majorana fermion mode for a specific mathematical relation between the
Rabi frequencies oscillation, the atom-photon coupling strengths and the laser
field detuning.
To the best of our knowledge the quantum simulation physics for the different
kind of fermionic mode in a same cavity QED system is absent in the
literature.
In the present study, we simulate two model Hamiltonians through the proper
tuning of Rabi frequencies and the atom-photon coupling strengths in the
system which quantum simulate different modes of fermionic excitations.
The Model Hamiltonian:
The Hamiltonian of our present study consists of three parts:

$H~{}=~{}{H_{A}}~{}+~{}{H_{C}}~{}+~{}{H_{AC}}$ (1)

The Hamiltonians are the following

${H_{A}}~{}=~{}\sum_{j=1}^{N}{{\omega}_{e}}|e_{j}><e_{j}|~{}+~{}{\omega}_{ab}|b_{j}><b_{j}|$
(2)

where $j$ is the cavity index. ${\omega}_{ab}$ and ${\omega}_{e}$ are the
energies of the state $|b>$ and the excited state respectively. The energy
level of state $|a>$ is set as zero. $|a>$ and $|b>$ are the two stable state
of a atom in the cavity and $|e>$ is the excited state of that atom in the
same cavity. The following Hamiltonian describes the photons in the cavity,

${H_{C}}~{}=~{}{{\omega}_{C}}\sum_{j=1}^{N}{{a_{j}}}^{\dagger}{a_{j}}~{}+~{}{J_{C}}\sum_{j=1}^{N}({{a_{j}}}^{\dagger}{a_{j+1}}+h.c),$
(3)

where ${a_{j}}^{\dagger}({a_{j}})$ is the photon creation (annihilation)
operator for the photon field in the $j$’th cavity, ${\omega}_{C}$ is the
energy of photons and $J_{C}$ is the tunneling rate of photons between
neighboring cavities. The interaction between the atoms and photons and also
by the driving lasers are described by

${H_{AC}}~{}=~{}\sum_{j=1}^{N}[(\frac{{\Omega}_{a}}{2}e^{-i{{\omega}_{a}}t}+{g_{a}}{a_{j}})|e_{j}><a_{j}|+h.c]+[a\leftrightarrow
b].$ (4)

Here ${g_{a}}$ and ${g_{b}}$ are the couplings of the cavity mode for the
transition from the energy states $|a>$ and $|b>$ to the excited state.
${\Omega}_{a}$ and ${\Omega}_{b}$ are the Rabi frequencies of the lasers with
frequencies ${\omega}_{a}$ and ${\omega}_{b}$ respectively.
The authors of Ref. hart1 ; hart2 ; sujop  have derived an effective spin
model by considering the following physical processes: A virtual process
regarding the emission and absorption of photons between the two stable states
of neighboring cavity yields the resulting effective Hamiltonian as

${H_{xy}}=\sum_{j=1}^{N}B{{\sigma}_{j}}^{z}~{}+~{}\sum_{j=1}^{N}(\frac{J_{1}}{2}{{\sigma}_{j}}^{\dagger}{{\sigma}_{j+1}}^{-}~{}+~{}\frac{J_{2}}{2}{{\sigma}_{j}}^{-}{{\sigma}_{j+1}}^{-}+h.c)$
(5)

When $J_{2}$ is real then this Hamiltonian reduces to the XY model. Where
${{\sigma}_{j}}^{z}=|b_{j}><b_{j}|~{}-~{}|a_{j}><a_{j}|$,
${{\sigma}_{j}}^{+}=|b_{j}><a_{j}|$, ${{\sigma}_{j}}^{-}=|a_{j}><b_{j}|$ .

$H_{xy}=\sum_{i=1}^{N}B({{\sigma}_{i}}^{z}~{}+~{}{J_{x}}{{\sigma}_{i}}^{x}{{\sigma}_{i+1}}^{x}~{}+~{}{J_{y}}{{\sigma}_{i}}^{y}{{\sigma}_{i+1}}^{y}).$
(6)

With ${J_{x}}=(J_{1}+J_{2})$ and ${J_{y}}=(J_{1}-J_{2})$.
We follow the references james ; hart1 , to present the analytical expression
for the different physical parameters of the system.
$B=\frac{\delta_{1}}{2}-\beta$, $\beta$ is define in Ref. beta .

${J_{1}}=\frac{\gamma_{2}}{4}(\frac{{|{\Omega_{a}}|}^{2}{g_{b}}^{2}}{{{\Delta}_{a}}^{2}}+\frac{{|{\Omega_{b}}|}^{2}{g_{a}}^{2}}{{{\Delta}_{b}}^{2}}),{J_{2}}=\frac{\gamma_{2}}{2}(\frac{{\Omega_{a}}{\Omega_{b}}g_{a}g_{b}}{{\Delta}_{a}{\Delta_{b}}}).$
(7)

Where $\gamma_{a,b}=\frac{1}{N}\sum_{k}\frac{1}{{\omega}_{a,b}-{\omega}_{k}}$
$\gamma_{1}=\frac{1}{N}\sum_{k}\frac{1}{({\omega}_{a}+{\omega}_{b})/2-{\omega}_{k}}$
and
$\gamma_{2}=\frac{1}{N}\sum_{k}\frac{e^{ik}}{({\omega}_{a}+{\omega}_{b})/2-{\omega}_{k}}$
${\delta_{1}}={\omega}_{ab}-({\omega}_{a}-{\omega}_{b})/2$,
${\Delta}_{a}={\omega}_{e}-{\omega}_{a}$.
${\Delta}_{b}={\omega}_{e}-{\omega}_{a}-({\omega}_{ab}-{\delta_{1}})$.
${{\delta}_{a}}^{k}={\omega}_{e}-{\omega}_{k}$,
${{\delta}_{b}}^{k}={\omega}_{e}-{\omega}_{k}-({\omega}_{ab}-{\delta_{1}})$,
${\omega}_{k}={\omega}_{c}+J_{c}\sum_{k}cosk$. $g_{a}$ and $g_{b}$ are the
couplings of respective transition to the cavity mode, ${\Omega}_{a}$ and
${\Omega}_{b}$ are the Rabi frequency of laser with frequency $\omega_{a}$ and
$\omega_{b}$.

Quantum Simulation for Dirac Fermion mode and Mathematical Relation Between
the Fields :
Here we quantum simulate the Dirac fermion physics through the proper tuning
of cavity QED lattice parameters. This condition can be achieve when
$J_{x}=J_{y}$. This condition implies that $2J_{2}=0$. It is clear from the
analytical expression of $J_{1}$ and $J_{2}$ that to satisfy the condition one
of the atom-photon coupling strength, i.e., $g_{a}$ or $g_{b}$ be zero or one
of the Rabi frequency oscillation ${\Omega}_{a}$ or ${\Omega}_{b}$ should be
zero. The analytical expression for $J_{1}$ become
$J_{1}=\frac{\gamma_{2}}{4}\frac{{|{\Omega_{a}}|}^{2}{g_{b}}^{2}}{{\Delta_{a}}^{2}}$
or
$J_{1}=\frac{\gamma_{2}}{4}\frac{{|{\Omega_{b}}|}^{2}{g_{a}}^{2}}{{\Delta_{b}}^{2}}$.
The other condition is ${\delta_{1}}=2\beta$, i.e.,
$2{\omega}_{ab}-({\omega}_{a}-{\omega}_{b})=4\beta$. Therefore the condition
for quantum simulation relates upto the microscopic level.
In this limit the Hamiltonian reduced to

$H=J\sum_{i=1}^{N}({\sigma_{i}}^{x}{\sigma_{i+1}}^{x}+{\sigma_{i}}^{y}{\sigma_{i+1}}^{y})$
(8)

After the Jordan-Wigner transformation and Abelian Bosonization study one can
write the above Hamiltonian as gia ; frad

$H_{0}=\sum_{s}\int\frac{dk}{2\pi}\epsilon(k){\psi_{s}}^{\dagger}(k){\psi_{s}}(k)$
(9)
$H_{0}=\sum_{s}\int_{-\Lambda}^{\Lambda}\frac{dk}{2\pi}(kv_{F})({{\psi}_{s,R}}^{\dagger}(k){{\psi}_{s,R}}(k)-{{\psi}_{s,L}}^{\dagger}(k){{\psi}_{s,L}}(k)).$
(10)

Here we use ${\epsilon}_{k}=kv_{F}$ near the Fermi points and
${{\psi_{s}}}^{\dagger}(x)=({{\psi}_{s,R}}^{\dagger}(x),{{\psi}_{s,L}}^{\dagger}(x))$.
Following the discussions in addendum, we can write the above Hamiltonian in
the following form of Dirac equation without any mass term.

$H=2J\int dx\bar{\psi}(i{\gamma_{1}}){\partial_{x}}{\psi}$ (11)

and ${\psi}^{\dagger}(x)=({{\psi}_{R}}^{\dagger},{{\psi}_{L}}^{\dagger})$.
${\bar{\psi}}^{\dagger}(x)=({{\psi}_{L}}^{\dagger},{{\psi}_{R}}^{\dagger})$,
$\psi_{s}(x)=\int\frac{dk}{2\pi}{\psi_{s}}(k)e^{ikx}$. It is customery to
introduce the Dirac matrices in Dirac equation. The ${\gamma}$ matrices are
the following. ${\gamma}^{0}=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right)$ ,
${\gamma}^{1}=-i{\sigma}_{y}=\left(\begin{array}[]{cc}0&-1\\\
1&0\end{array}\right)$ ,
${\gamma}^{5}={\gamma_{0}}{\gamma_{1}}=\sigma_{z}=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)$ ,
${\psi}_{R}$ and ${\psi}_{L}$ are the fermionic field for the right and left
movers electron. The scalar field and its dual can be expressed as
$\phi(x)=\phi_{R}(x)+\phi_{L}(x)$, $\theta(x)=-{\phi}_{R}(x)+{\phi}_{L}(x)$.
The detail derivation of the analytical relation between the wave function and
the scalar field are relegated to the addendum. Here we only present the final
form Now we present the mathematical relation between the operators,
$e^{\pm\phi(x)}$ and $e^{\pm\theta(x)}$.

$e^{i\alpha\sqrt{\pi}\phi(x1)}e^{i\beta\sqrt{\pi}\theta(x2)}=e^{i\beta\sqrt{\pi}\theta(x2)}e^{i\alpha\sqrt{\pi}\phi(x1)}~{}sign(x2-x1)$
(12)
$e^{i\alpha\sqrt{\pi}\phi(x1)}{\psi_{R/L}}(x2)={\psi_{R/L}}(x2)e^{i\alpha\alpha\sqrt{\pi}\phi(x1)}~{}sign(x2-x1)$
(13)
$e^{i\beta\sqrt{\pi}\theta(x1)}{\psi_{R/L}}(x2)=-{\psi_{R/L}}(x2)e^{i\beta\sqrt{\pi}\theta(x1)}~{}sign(x2-x1).$
(14)

Therefore it is clear from our theoretical analysis that the cavity QED
lattice shows the Dirac fermion like mode when one of the atom-photon coupling
( $g_{a}$ or $g_{b}$ ) is zero or one of the applied Rabi frequencies (
${\Omega}_{a}$, ${\Omega}_{b}$) is zero.

Dirac Equation for Majorana Fermion mode and Condition of Majorana-Weyl
fermionic mode:
Here we quantum simulate the Majorana fermion mode and the Majorana-Weyl
fermion mode. Majorana thought whether it might be possible for spin-1/2
particle to be its own anti-particle. To get an equation alike to Dirac
equation but capable of governing a real field that requires the $\gamma$
matrices of that equation must satisfy the Clifford algebra are purely
imaginary. Here we show explicitly that the $\gamma$ matrices of Majorana
equation are purely imaginary and field satisfy the particle-antiparticle
equivalent condition.
We also derive the condition for the excitation of Majorana-Weyl fermion mode
where the simulated mode shows the gapless excitation.
We consider $J_{1}=J_{2}$, $J_{x}$ become $J_{1}+J_{2}$ and $J_{y}=0$. In the
microcavity array, the condition for $J_{1}=J_{2}$ achieve when

${{\Omega}_{a}}^{2}{g_{b}}^{2}{\Delta_{b}}^{2}+{{\Omega}_{b}}^{2}{g_{a}}^{2}{\Delta_{a}}^{2}=2{\Omega}_{a}{\Omega}_{b}g_{a}g_{b}{\Delta}_{a}{\Delta}_{b}.$
(15)

The above condition implies that
${\Omega}_{a}={\Omega}_{b}\frac{g_{a}\Delta_{a}}{g_{b}\Delta_{b}}$. The only
constraint is that ${\Delta}_{a}\neq{\Delta}_{b}$, the magnetic field diverge
when ${\Delta}_{a}={\Delta}_{b}$. At the same time,
${\Omega}_{a}={\Omega}_{b}$ and $g_{a}=g_{b}$ are also not possible because
this limit also leads to the condition ${\Delta}_{a}={\Delta_{b}}$. Suppose we
consider, ${\Omega}_{a}=\alpha_{1}{\Omega}_{b}$, $g_{a}=\alpha_{2}g_{b}$ and
${\Delta}_{a}=\alpha_{3}{\Delta}_{b}$. These relations implies that
${\alpha_{1}}^{2}+{\alpha_{2}}^{2}{\alpha_{3}}^{2}=2{\alpha}_{1}{\alpha}_{2}{\alpha}_{3}$.
${\alpha}_{1}={\alpha}_{2}{\alpha}_{3}$, ${\alpha}_{1},{\alpha_{2}}$ and
${\alpha_{3}}$ are the numbers. These analytical relations help to implement
the transverse Ising model Hamiltonian but $\alpha_{1}$, $\alpha_{2}$ and
$\alpha_{3}$ should not be equal to 1.
The quantum state engineering of cavity QED is in the state of art due to the
rapid progress of technological development of this field [1]. Therefore one
can achieve this limit to get the desire quantum state. One can write the
final Hamiltonian as,

$H_{T}=B\sum_{j=1}^{N}({{\sigma}_{z}}(j)+\lambda{{\sigma}_{x}}(j){{\sigma}_{x}}(j+1)),$
(16)

where $\lambda=\frac{J_{1}+J_{2}}{B}$. The effective Hamiltonian become the
transverse Ising model which studied in the previous literature ss ; mussardo
; druff .

One can also write the starting Hamiltonian as through the rotation in spin
basis

$H=\sum_{n}[\lambda\sigma_{z}(n)\sigma_{z}(n+1)+s_{x}(n)].$ (17)

We recast the Hamiltonian in the following form because it will help us to use
the order and disorder operator directly to the derivation of the equation of
motion and finally the Dirac equation for Majorana fermion.
Here our main motivation is to use some of important results of this model
Hamiltonian to discuss the relevant physics of array of cavity QED system. We
introduce the order and disorder operator in the addendum mussardo ; druff .
These operators are defining the sites of the dual lattice, i.e., we define
the operator between the nearest-neighbor site of the original lattice. Here
we define the Dirac spinor, $\chi_{1}(n)=\sigma_{z}(n)\mu_{z}(n+1/2)$ and
$\chi_{2}(n)=\sigma_{z}(n)\mu_{z}(n-1/2)$.
Now our main task is to find the equation of motion for the operators,
${\sigma}_{3}(n)$ and ${\mu}_{3}(n)$ which help us to build the Dirac
equation. The detail derivation are relagated to the appendix.
The equation of motion for the ${\sigma}_{z}(n)$ is the following:

$\frac{\partial\sigma_{z}(n)}{\partial\tau}=[H,\sigma_{z}(n)]={\sigma_{x}(n)\sigma_{z}(n)}$
(18)

The equation of motion for $\mu_{z}(n+1/2)$ is the following:

$\displaystyle\frac{\partial\mu_{z}(n+1/2)}{\partial\tau}$ $\displaystyle=$
$\displaystyle\lambda{\mu_{x}(n+1/2)\mu_{z}(n+1/2)}$ (19)
$\displaystyle=\lambda\sigma_{z}(n)\sigma_{z}(n+1/2)\mu_{z}(n+1/2)$
$\frac{\partial\chi_{1}(n)}{d\tau}=-\chi_{2}(n)+\lambda\chi_{2}(n+1).$ (21)
$\frac{\partial\chi_{2}(n)}{d\tau}=-\chi_{1}(n)+\lambda\chi_{1}(n-1).$ (22)

These two fields, $\chi_{1}(n)$ and $\chi_{2}(n)$ satisfy the following
relations, $\\{\chi_{1}(n1),\chi_{2}(n2)\\}=2\delta_{n1,n2}$. One can write
down the above equation in the following compact form,

$({\gamma}^{0}\frac{\partial}{\partial t}+{\gamma}^{3}\frac{\partial}{\partial
r}+m)\chi=0.$ (23)

where ${\chi}^{\dagger}=(\chi_{1},\chi_{2})$ and $m=\frac{1-\lambda}{\alpha}$,
${\gamma}^{0}=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right),{\gamma}^{3}=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)$ .
One can also write the above Majorana equation in a compact form:

$(i{\tilde{\gamma}}^{\mu}{\partial}_{\mu}-m){\chi}=0.$ (24)

Where ${\tilde{\gamma}}^{0}=\left(\begin{array}[]{cc}0&i\\\
i&0\end{array}\right),{\tilde{\gamma}}^{3}=\left(\begin{array}[]{cc}i&0\\\
0&-i\end{array}\right)$ .
Therefore we prove that the spinor field satisfy the Majorana condition of
spin-1/2 particle and also the $\tilde{\gamma}$ matrices are imaginary.
The condition for the massless Majorana fermion field is that
$m=0,i.e.,\lambda=1$. In this quantum simulation process, one can get through
this massless excitation through this analytical relation.
$\frac{\gamma_{2}}{2}(\frac{{|{\Omega_{a}}|}^{2}{g_{b}}^{2}}{{{\Delta}_{a}}^{2}}+\frac{{|{\Omega_{b}}|}^{2}{g_{a}}^{2}}{{{\Delta}_{b}}^{2}})=\frac{\omega_{ab}}{2}-\frac{1}{4}(\omega_{a}-\omega_{b})-\beta$.
We term this fermionic mode as Majorana-Weyl fermionic mode. In this study we
obtain the different kind of fermionic modes from the quantum simulation of
the system, here there is no splitting of Dirac fermion mode into the Majorana
fermion mode or Majorana-Weyl fermionic mode. There are no studies in the
previous literature of cavity QED where one has found a crossover from Dirac
fermion modes to Majorana fermion modes [9-11].

Here we present the analytical relations between the Majorana fermion
operators with the order and disorder operator with the free Dirac field in
Abelian bosonization theory. Then one obtains the following sets of
commutation relations.

${\sigma}_{z}(x1){\mu}_{z}(x2)={\mu}_{z}(x2){\sigma}_{z}(x1)sign(x1-x2)$ (25)
${\sigma}_{z}(x1){\chi}(x2)={\chi}(x2){\sigma}_{z}(x1)sign(x1-x2)$ (26)
${\mu}_{z}(x1){\chi}(x2)=-{\chi}(x2){\sigma}_{z}(x1)sign(x1-x2)$ (27)

It is very clear from the above analytical relations that
${{\chi}_{1}}^{\dagger}={\chi}_{1}$ and ${{\chi}_{2}}^{\dagger}={{\chi}_{2}}$.
The detail derivation is relegated in the appendix. The above relation has
similarity with the free Dirac field in Abalian bosonization theory, where
Dirac field operator is a local product of two phase exponential depending on
the scalar field and its dual gia ; gogo ; frad , as one study the Luttinger
liquid physics in Abelian bosonization theory.
Conclusions
We have presented the existence of Dirac fermion mode, Majorana fermion mode
and also Majorana-Weyl fermionic mode for the optical cavity array with the
relation between Rabi frequency oscillation and the atom photon coupling
strength. We have also presented the crossover between the Dirac fermion modes
to Majorana fermion modes and also the condition for appearence of Majorana-
Weyl massless excitation mode. We have also presented several analytical
relations between the Majorana field, order and disorder operators.

Acknowledgement: The author would like to acknowledge the discussions with
Prof. S. Girvin during the international workshop/school on Dirac Materials
and Chandrashekar lecture at ICTS and the library and LAMP group of Raman
Research Institute. The author would like to thank Dr. P. K. Mukherjee for
reading the manuscript carefully. Finally the author would like to acknowledge
the DST (SERC) project and HRI Library.

## References

  * (1) R. Gerritsma $et~{}al.$, Nature 463, 68 (2010).
  * (2) Kangjun Seo, Li Han and C. A. R. Sa de Melo, Phys. Rev. Lett 109, 105303 (2012).
  * (3) Liang Fu and C. L. Kane, Phys. Rev. Lett (2009).
  * (4) L. Lamta, J. Leo’n, T. Schatz and E. Solano, Phys. Rev. Lett (2007).
  * (5) A. Bermudez, M. A. Martin-Delgado and E. Solano, Phys. Rev. A 76 041801 (2007).
  * (6) , X. Zhang, Phys. Rev. Lett 100, 113903 (2008).
  * (7) J. Otterbach, R. G. Unanyan and M. Fleishhauer, Phys. Rev. Lett 102, 0636602 (2008).
  * (8) J. Y. Vaishanav and C. W. Clark, Phys. Rev. Lett 100, 153002 (2006).
  * (9) C. E. Bardyn and A. Imamoglu, arXiv(cond-mat)/1204.1238 .
  * (10) Brijesh Kumar and Somenath Jalal, arXiv(cond-mat)/1210.6922.
  * (11) Mircea Trif and Yaroslav Tserkovnyak, arXiv (cond-mat)/1202.2694.
  * (12) G. Agarwal, ”Quantum Optics”, (Cambridge University Press, Delhi 2013).
  * (13) Greentree. A D $et~{}al.$ 2006, Nature Phys. 466 856\.
  * (14) Hartmann Michael J, Fernando G S, Brando L and Plenio Martin B 2006 Nature Phys 462 849; Hartmann Michael J, Fernando G S, Brando L and Plenio Martin B 2008, Laser and Photonics Rev. 2 527\.
  * (15) Hartmann Michael J, Fernando G S, Brando L and Plenio Martin B 2007, Phys. Rev. Lett 99 160501\.
  * (16) S. Sarkar, Physica B 407, 44 (2012).
  * (17) A. C. Ji, X. C. Xie, and W. M. Liu , Phys. Rev.Lett. 99, 183602 (2007).
  * (18) T. Byrnes, N. Y. Kim, K. Kusudo, and Y. Yamamoto, Phys. Rev. B 78 075320 (2008).
  * (19) I. Carusotto I $et~{}al.$, arXiv:0812.4195 (2008).
  * (20) M. J. Bhaseen, M. Hohenadler, A. O. Silver, and B. D. Simons, Phys. Rev. Lett. 102 135301 (2009).
  * (21) J. Zhao, A. W. Sandvik and K. Ueda, arXiv:0806.3603 (2008).
  * (22) P. Pippan, H. G. Evertz, and M. Hohenadler , arXiv: 0904.1350 (2009).
  * (23) M. Aichhorn $et~{}al.$, Phys. Rev. Lett. 100 216401 (2008).
  * (24) S. Schmidt, and G. Blatter, arXiv:0905.3344 (2009).
  * (25) M. I. Hwang and M. S. Choi, arXiv:1207.0082 (2012).
  * (26) D. G. Angelakis, M. F. Santos and S. Bose, Phys. Rev. A 76 R031805 (2007).
  * (27) S. Horoche and J. M. Raimond 2006 in Exploring the Quantum Atoms, Cavities, and Photons, (Oxford University Press).
  * (28) E. Majorana, Nuovo Cimento 14, 171 (1937).
  * (29) F. Wilczek, Majorana returns, Nature Physics 5, 614
  * (30) G. Moore and N. Read, Nucl. Phys. B 360, 362 (1991).
  * (31) N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
  * (32) A. Y. Kitaev, Physics-Uspekhi 44, 131 (2001).
  * (33) M. Wimmer $et~{}al.$, Phys. Rev. Lett 105, 046803 (2010).
  * (34) N. B. Kopnin and M. M. Salomaa, Phys. Rev. B 44, 9667 (1991).
  * (35) C. Nayak $et~{}al.$, Rev. Mod. Phys 80, 1083 (2008).
  * (36) F. V. James and J. Jerke, Can. J. Phys 85, 625 (2007).
  * (37) $\displaystyle\beta$ $\displaystyle=$ $\displaystyle\frac{1}{2}[\frac{{|{\Omega_{b}}|}^{2}}{4{\Delta}_{b}}({\Delta}_{b}-\frac{{|{\Omega_{b}}|}^{2}}{4{\Delta}_{b}}-\frac{{|{\Omega_{b}}|}^{2}}{4({\Delta}_{a}-{\Delta}_{b})}-{\gamma_{b}}{g_{b}}^{2}-{\gamma_{1}}{g_{a}}^{2}$ (28) $\displaystyle+{\gamma_{1}}^{2}\frac{{g_{a}}^{4}}{{\Delta_{b}}}-(a\leftrightarrow b)].$
  * (38) Subir Sachdev, Quantum Phase Transitions, Cambridge University Press, 2001.
  * (39) G. Mussardo, ”Statistical Field Theory”, Oxford Graduate Texts (New Delhi, 2010).
  * (40) Itzyaskon and Druff, ”Statistical Field Theory”, (Cambridge University Press, New York 1989).
  * (41) T. Giamarchi, ” Quantum Physics in One Dimension”.
  * (42) A.O. Gogolin, A. A. Nersesyan and A. M. Tsvelik, ”Bosonization and Strongly Correlated Systems”, Cambridge University Press (New York, 1998).
  * (43) E. Fradkin, ””Field Theories in Condensed Matter Physics”, Cambridge University Press, Delhi 2013. .

Addendum
Dirac fermion mode and analytical relation between fields:

One dimensional quantum mechanical system both theoretically and
experimentally has studied extensively in last two decades with several
successful explanation for different physical systems. There are several
excellent research articles are available in the literature. Here we discuss
very briefly about the Luttinger liquid physics of one-dimensional quantum
many body system.
The low energy excitation of a non-interacting system. In one dimensional
quantum mechanical system consists of Fermi points at $\pm k_{F}$. The low
energy fermionic state thus have momentum, $k\sim\pm k_{F}$ and a single
particle energy close to Fermi energy , $E_{F}$.

$E(k)\simeq E_{F}+(|k|-k_{F})v_{F}+........$ (29)

Here we consider the one-dimensional infinite length system and the elementary
excitations across the Fermi points. Therefore we can write the Fermi field,

${\psi}_{s}(x)=\int\frac{dk}{2\pi}{\psi}_{s}(k)e^{ikx}$ (30)

Here we consider the modes of momentum expansion in a neighbourhood of $\pm
k_{F}$ of width $2\Lambda$. Therefore, we can write the fermionic field,

${\psi}_{s}(x)\simeq\int_{-\Lambda}^{\Lambda}\frac{dk}{2\pi}e^{i(k+k_{F})x}\psi_{s}(k+k_{F})+\int_{-\Lambda}^{\Lambda}\frac{dk}{2\pi}e^{i(k-k_{F})x}\psi_{s}(k-k_{F})$
(31)

One can write this analytical expression as,

${\psi}_{s}(x)\simeq e^{ik_{F}x}{\psi}_{s,R}(x)+e^{-ik_{F}x}{\psi}_{s,L},$
(32)

where ${\psi}_{s,R}(k)={\psi}_{s}(k+k_{F})$ and
${\psi}_{s,L}(k)={\psi}_{s}(k-k_{F})$. The free fermion Hamiltonian,

$H_{0}=\sum_{s}\int\frac{dk}{2\pi}\epsilon(k){\psi_{s}}^{\dagger}(k){\psi_{s}}(k)$
(33)
$H_{0}=\sum_{s}\int_{-\Lambda}^{\Lambda}\frac{dk}{2\pi}(kv_{F})({{\psi}_{s,R}}^{\dagger}(k){{\psi}_{s,R}}(k)-{{\psi}_{s,L}}^{\dagger}(k){{\psi}_{s,L}}(k)).$
(34)

Here we use ${\epsilon}_{k}=kv_{F}$ near the Fermi points and
${\psi_{s}}^{\dagger}(x)=({{\psi}_{s,R}}^{\dagger}(x),{{\psi}_{s,L}}^{\dagger}(x))$.

$H_{0}=\sum\int\frac{dk}{2\pi}{\psi_{s}}^{\dagger}(k)\sigma_{3}kv_{F}{\psi}_{s}(k).$
(35) $H_{0}=\sum_{s}\int
dx{\psi_{s}}^{\dagger}(x)\sigma_{3}kiv_{F}{\partial_{x}}{\psi}_{s}(k),$ (36)

where ${\sigma}^{3}=\left(\begin{array}[]{cc}1&0\\\ 0&-1\end{array}\right)$.
We can define scalar field ($\phi(x)$ ) and dual field ($\theta(x)$) by the
following relations
$\phi(x)=\phi_{R}(x)+\phi_{L}(x)$, $\theta(x)=-{\phi}_{R}(x)+{\phi}_{L}(x)$.
One canwrite the Dirac fermion density
$\rho(x)=\rho_{R}(x)+\rho_{L}(x)=-\frac{1}{\sqrt{\pi}}{\partial_{x}}{\phi}_{x}$
and the corresponding current density
$j(x)=v({\rho}_{R}-{\rho}_{L}=\frac{v_{F}}{\sqrt{\pi}}{\partial_{x}}{\theta}(x)$.
The commutattion relation between the scalar field (${\phi}$) and its dual
field ($\theta$).
$[\phi(x1),\phi(x2)]=0=[\theta(x1),\theta(x2)]$,
$[\phi(x1),\theta(x2)]=i~{}sign(x2-x1)$, $[\phi(x1),\Pi(x2)]=i\delta(x1-x2)$.
Where $\Pi(x)=\partial_{x}\theta(x)$. The corresponding Lagrangian debnsity
for massless Dirac fermions is

$L=i{\psi_{R}}^{\dagger}(\partial_{t}+v_{F}\partial_{x}){\psi_{R}}+i{\psi_{L}}^{\dagger}(\partial_{t}-v_{F}\partial_{x}){\psi_{L}}$
(37)

Now we present the algebra of the operators, $e^{\pm\phi(x)}$ and
$e^{\pm\theta(x)}$.

$e^{i\alpha\sqrt{\pi}\phi(x1)}e^{i\beta\sqrt{\pi}\theta(x2)}=e^{i\beta\sqrt{\pi}\theta(x2)}e^{i\alpha\sqrt{\pi}\phi(x1)}~{}sign(x2-x1)$
(38)
$e^{i\alpha\sqrt{\pi}\phi(x1)}{\psi_{R/L}}(x2)={\psi_{R/L}}(x2)e^{i\alpha\alpha\sqrt{\pi}\phi(x1)}~{}sign(x2-x1)$
(39)
$e^{i\beta\sqrt{\pi}\theta(x1)}{\psi_{R/L}}(x2)=-{\psi_{R/L}}(x2)e^{i\beta\sqrt{\pi}\theta(x1)}~{}sign(x2-x1)$
(40)

Analytical relations between order and disorder operators:
The analytical relation between the Pauli operators and $\mu$ operators are
relagated to the appendix B of this manuscript, please put the following
analytical expression in the appendix:

${{\mu}_{z}}^{2}=1={{\mu}_{x}}^{2},$ (41)
${{\mu}_{z}}(n-1/2){{\mu}_{z}}(n+1/2)={\sigma}_{x}(n).$ (42)
${{\mu}_{x}}(n+1/2)={\sigma}_{z}(n){\sigma}_{z}(n+1),$ (43)
${{\mu}_{z}}(n+1/2)=\Pi_{j=1}^{n}{\sigma_{x}}(j).$ (44)
${{\sigma_{z}}}(n)=\Pi_{j=0}^{n-1}{\mu_{x}}(j+1/2),$ (45)
$[{\mu}_{x}(n+1/2),{\mu}_{z}(n^{{}^{\prime}}+1/2)]=2\delta_{n,n^{{}^{\prime}}}$
(46) $[{\mu}_{z}(n+1/2),{\mu}_{z}(n^{{}^{\prime}}+1/2)]=0,$ (47)
$[{\mu}_{z}(n+1/2),{\sigma}_{x}(n^{{}^{\prime}})]=0$ (48)

The operator ${\mu}_{z}(n+1/2)$ acting on the original spin of the lattice
makes a spin flip of all those spin placed on the left hand side of spin at
the site n. Therefore ${\mu}_{z}(n+1/2)$ is a kink operator, it introduce the
disorder in the system. It is very clear from the above analytical relation of
the operators that Here we express the analytical relation between the
Majorana operators and the disorder and Pauli operators:

${\chi}_{1}(n)={\sigma}_{z}(n){\mu}_{z}(n+1/2)=-{\mu}_{z}(n+1/2){\sigma}_{z}(n)$
(49)
${\chi}_{2}(n)={\sigma}_{z}(n){\mu}_{z}(n-1/2)={\mu}_{z}(n-1/2){\sigma}_{z}(n)$
(50)
${\sigma}_{z}(n){\chi}_{2}(n)={\mu}_{z}(n-1/2)={\chi}_{2}(n){\sigma}_{z}(n)$
(51)
${\sigma}_{z}(n){\chi}_{1}(n)={\mu}_{z}(n+1/2)=-{\chi}_{1}(n){\sigma}_{z}(n)$
(52)
${\sigma}_{z}(n)={\mu}_{z}(n-1/2){\chi}_{2}(n)={\chi}_{2}(n){\mu}_{z}(n-1/2)$
(53)
${\mu}_{z}(n+1/2){\chi}_{1}(n)=-{\chi}_{1}(n){\mu}_{z}(n+1/2)={\sigma}_{z}(n)$
(54)

Derivation of Dirac equation for Majorana fermion field
The equation of motion for the ${\sigma}_{z}(n)$ is the following:

$\frac{\partial\sigma_{z}(n)}{\partial\tau}=[H,\sigma_{z}(n)]={\sigma_{x}(n)\sigma_{z}(n)}$
(55)

The equation of motion for $\mu_{z}(n+1/2)$ is the following:

$\displaystyle\frac{\partial\mu_{z}(n+1/2)}{\partial\tau}$ $\displaystyle=$
$\displaystyle\lambda{\mu_{x}(n+1/2)\mu_{z}(n+1/2)}$ (56)
$\displaystyle=\lambda\sigma_{z}(n)\sigma_{z}(n+1/2)\mu_{z}(n+1/2)$

Now we use the properties of the $\sigma$ and $\mu$ operators to derive the
equation of motion for the Majorana fields $\chi_{1}(n)$ and $\chi_{2}(n)$.

$\frac{\partial\chi_{1}(n)}{d\tau}=\frac{\partial\sigma_{z}(n)}{\partial\tau}\mu_{z}(n+1/2)+\sigma_{z}(n)\frac{\partial\mu_{z}(n)}{\partial\tau}.$
(58)
$\frac{\partial\chi_{1}(n)}{d\tau}=\sigma_{x}(n)\sigma_{z}(n)\mu_{z}(n+1/2)+\lambda\sigma_{z}(n)\sigma_{z}(n)\sigma_{z}(n+1)\mu_{z}(n+1/2).$
(59) $\displaystyle\frac{\partial\chi_{1}(n)}{d\tau}$ $\displaystyle=$
$\displaystyle-\sigma_{z}(n)\mu_{z}(n-1/2)\mu_{z}(n+1/2)\mu_{z}(n+1/2)$ (60)
$\displaystyle+\lambda\sigma_{z}(n)\sigma_{z}(n)\sigma_{z}(n+1)\mu_{z}(n+1/2).$
$\frac{\partial\chi_{1}(n)}{d\tau}=-\chi_{2}(n)+\lambda\chi_{2}(n+1).$ (61)

Now the equations of motion for $\chi_{2}(n)$ are

$\displaystyle\frac{\partial\chi_{2}(n)}{d\tau}$ $\displaystyle=$
$\displaystyle\frac{\partial\sigma_{z}(n)}{\partial\tau}\mu_{z}(n-1/2)$ (62)
$\displaystyle+\sigma_{z}(n)\frac{\partial\mu_{z}(n-1/2)}{\partial\tau}$
$\displaystyle\frac{\partial\chi_{2}(n)}{d\tau}$ $\displaystyle=$
$\displaystyle\sigma_{x}(n)\sigma_{z}(n)\mu_{z}(n-1/2)$ (63)
$\displaystyle+\lambda\sigma_{z}(n)\sigma_{z}(n-1)\sigma_{z}(n)\mu_{z}(n-1/2)$
$\displaystyle\frac{\partial\chi_{2}(n)}{d\tau}$ $\displaystyle=$
$\displaystyle\mu_{z}(n-1/2)\mu_{z}(n+1/2)\sigma_{z}(n)\mu_{z}(n-1/2)$ (64)
$\displaystyle+\lambda\sigma_{z}(n-1)\mu_{z}(n-1/2).$

After a little bit of calculations and using the relation between the disorder
operators (Eq. 23-30), we finally arrive the equation of motion of
$\chi_{2}(n)$ as,

$\frac{\partial\chi_{2}(n)}{d\tau}=-\chi_{1}(n)+\lambda\chi_{1}(n-1).$ (65)

These two fields, $\chi_{1}(n)$ and $\chi_{2}(n)$ satisfy the following
relations, $\\{\chi_{1}(n1),\chi_{2}(n2)\\}=2\delta_{n1,n2}$. One can write
down the above equation in the following compact form,

$({\gamma}^{0}\frac{\partial}{\partial t}+{\gamma}^{3}\frac{\partial}{\partial
r}+m)\chi=0.$ (66)

where ${\chi}^{\dagger}=(\chi_{1},\chi_{2})$ and $m=\frac{1-\lambda}{\alpha}$,
${\gamma}^{0}=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right),{\gamma}^{3}=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)$ .


# Viable massive gravity without nonlinear screening

Yusuke Manita<EMAIL_ADDRESS>Department of Physics, Kyoto
University, Kyoto 606-8502, Japan  Sirachak Panpanich
<EMAIL_ADDRESS>Department of Pure and Applied Physics, Graduate
School of Advanced Science and Engineering, Waseda University, Okubo 3-4-1,
Shinjuku, Tokyo 169-8555, Japan  Rampei Kimura<EMAIL_ADDRESS>Waseda
Institute for Advanced Study, Waseda University,1-6-1 Nishi-Waseda, Shinjuku,
Tokyo 169-8050, Japan

###### Abstract

We study nonlinear effects of perturbations around a cosmological background
in projected massive gravity, which admits self–accelerating solutions in an
open FLRW universe. Using the zero-curvature scaling limit, we derive
nonlinear equations containing all the relevant terms on subhorizon scales. We
find that the solution for a scalar graviton vanishes completely for all
scales, which agrees with the linear perturbation analysis in the previous
study. In addition, the effects on the gravitational potential due to the next
order perturbation are strongly suppressed within the horizon. Therefore, a
screening mechanism is no longer needed for consistency with solar-system
experiments in the projected massive gravity.

††preprint: KUNS-2951

##  I Introduction

In the last decade there has been considerable progress in the discovery and
development of ghost-free massive gravity, so called de Rham-Gabadadze-Tolley
(dRGT) massive gravity de Rham and Gabadadze (2010); de Rham  _et al._ (2011).
The theory is free from the Boulware-Deser (BD) ghost Boulware and Deser
(1972), typically appearing due to the breaking of general covariance, it thus
has only $5$ degrees of freedom ($2$ tensor $+$ $2$ vector $+$ $1$ scalar
modes). Although the dRGT massive gravity admits an open-Friedmann-Lemaître-
Robertson-Walker (FLRW) universe with the effective cosmological constant
driven by a mass term, it unfortunately suffers from nonlinear instabilities
De Felice  _et al._ (2012) associated with a strong coupling of the scalar and
vector graviton Gumrukcuoglu  _et al._ (2011). For this reason, the dRGT model
is not a viable theory of massive gravity, and various extension has been
intensively investigated in the context of new derivative interactions
Hinterbichler (2013); Kimura and Yamauchi (2013); de Rham  _et al._ (2014a);
Gao (2014); Folkerts  _et al._ (2011), adding an extra scalar degree of
freedom Huang  _et al._ (2013); D’Amico  _et al._ (2013), a minimal theory of
massive gravity De Felice and Mukohyama (2016); De Felice  _et al._ (2021,
2017a, 2017b) and bigravity Hassan and Rosen (2012).

One interesting and minimal way to extend the dRGT massive gravity without
introducing an extra degree of freedom is to break the translation invariance
of the Stückelberg fields while preserving the global Lorentz invariance, and
it allows to include new Lorentz-invariants in the Lagrangian. The simplest
extension where the constant parameters in the potential of the dRGT theory
are promoted to be arbitrary functions of the Lorentz-invariant is called
generalized massive gravity (GMG) theory de Rham  _et al._ (2014b, c).
Furthermore, with an idea of the degeneracy of the Stückelberg kinetic term,
the breaking of the translational invariance can introduce a nonminimal
coupling of the scalar curvature to the Lorentz-invariant quantity of the
Stückelberg field and two different classes of graviton mass terms without
invoking the BD ghost Gumrukcuoglu  _et al._ (2021). The first class is an
straightforward extension of the dRGT theory (as well as GMG) obtained by
applying the conformal and disformal deformations of the fiducial metric. In
spite of the failure of the dRGT theory, the GMG theory has self-accelerating
solutions which are free from any instabilities on an open-FLRW background
Kenna-Allison  _et al._ (2020a). The second class is called projected massive
gravity (PMG) theory. The mass term is composed of a projection tensor which
eliminates one of Stückelberg fields and manifestly prevents the appearance of
the BD ghost. The PMG theory similarly admits a perturbatively stable open-
FLRW solution Gumrukcuoglu  _et al._ (2020).

In massive gravity, the scalar mode typically induces a fifth force comparable
to the Newtonian force. In the case of the GMG theory, this can be seen, for
example, in the effective gravitational coupling in the evolution equation for
the matter perturbation Kenna-Allison  _et al._ (2020b). Although modified
gravity theories requires a screening mechanism to hide the fifth force in
order to satisfy various observations on small scales such as solar-system
Will (2014), massive gravity originally possess the Vainshtein mechanism
Vainshtein (1972) as a solution of the vDVZ discontinuity van Dam and Veltman
(1970); Zakharov (1970). This can be easily seen in the effective Lagrangian
obtained by taking the decoupling limit of the dRGT theory de Rham and
Gabadadze (2010), where the galileon interactions play an important role, and
the additional force due to the scalar graviton can be strongly suppressed
inside a characteristic scale called the Vainshtein radius Babichev  _et al._
(2009); Berezhiani  _et al._ (2013). In fact, the Vainshtein mechanism
operates in the case of the GMG theory Gumrukcuoglu  _et al._ (2021).

On the other hand, in the PMG theory, the linear perturbation analysis around
a cosmological background shows the fifth force associated with the scalar
graviton vanishes at the leading order in the subhorizon expansion, i.e., the
absence of the vDVZ discontinuity Manita and Kimura (2022). This intuitively
implies that the Vainshtein mechanism is not necessary in the PMG theory.
However, galileon-like interactions which potentially dominate at small scales
should be present in nonlinear equations as in the GMG case. Thus it is not
confirmed yet whether the fifth force comparable to the Newtonian one is
absent even at small scales. To this end, we, in the present paper,
investigate nonlinear solution of the scalar graviton by taking into account
the cosmological background in order to check the consistency with the local
experiments.

The paper is organized as follows. In Sec. II, we introduce action and basic
equations of the theory. We study background and linear perturbation equations
in the zero-curvature limit, and quasi-static limit in Sec. III. In Sec. IV,
we investigate nonlinear perturbations in a spherically symmetric setup.
Lastly, Sec. V is devoted to discussion and conclusion.

In the following, we use the notation $c=1$, and
$\eta_{\mu\nu}=\rm{diag.}(-1,1,1,1)$. With the Newton constant $G_{N}$, the
Planck mass is defined by $M_{\rm Pl}:=(8\pi G_{N})^{-2}$.

##  II Projected massive gravity

In this section, we briefly introduce the ghost-free theory of projected
massive gravity (PMG) Gumrukcuoglu  _et al._ (2020). The mass term in the PMG
theory is defined through the following fiducial metric:

$\displaystyle\bar{f}_{\mu\nu}:=P_{ab}\partial_{\mu}\phi^{a}\partial_{\nu}\phi^{b}\,.$
(1)

Here, the four scalar fields $\phi^{a}$ are the so-called Stückelberg fields,
and $P_{ab}$ is the projection operator on the field space defined by

$\displaystyle P_{ab}=\eta_{ab}-\frac{\phi_{a}\phi_{b}}{X}\,,$ (2)

with $X:=\eta_{ab}\phi^{a}\phi^{b}$, which projects onto surfaces in the field
space of the Stückelberg fields $\phi^{a}$. With the building block matrix of
the mass term, $Z^{\mu}{}_{\nu}:=\left(g^{-1}\bar{f}\right)^{\mu}{}_{\nu}$,
the action of PMG is given by

$\displaystyle S$ $\displaystyle=\frac{M_{\rm Pl}^{2}}{2}\int
d^{4}x\sqrt{-g}\Bigg{[}GR+\frac{3G_{X}^{2}}{2G}\partial_{\mu}X\partial^{\mu}X$
$\displaystyle+m^{2}U\left(X,[Z],\left[Z^{2}\right],\left[Z^{3}\right]\right)\Bigg{]}+S_{\mathrm{m}}\left[g_{\mu\nu},\psi\right]\,,$
(3)

where $G$ is the arbitrary dimensionless function of $X$, $G_{X}:=\partial
G/\partial X$, $m$ is a constant mass parameter, $U$ is the arbitrary
dimensionless function of $X,[Z],[Z^{2}]$, and $[Z^{3}]$, and $S_{\rm m}$ is
the action of the matter field $\psi$ which minimally couples to the physical
metric $g_{\mu\nu}$. The existence of $X$-dependence in the action (3) as well
as the projection tensor (2) manifestly break the translation invariance of
the Stückelberg fields while the global Lorentz invariance is retained. In
other words, thanks to the violation of the internal translational symmetry,
PMG generally includes the non-minimal coupling of $\phi^{a}$ to the Ricci
scalar. Despite this fact, stable cosmological solutions have not been yet
found in the case of $G\neq 1$ Manita and Kimura (2022). For this reason, we
focus on the case with $G=1$ in the following part of this paper.

The modified Einstein equation is given by varying the action with respect to
$g^{\mu\nu}$,

$\displaystyle M_{\rm Pl}^{2}G_{\mu\nu}=T^{(\rm
mass)}_{\mu\nu}+T_{\mu\nu}^{(\rm m)}\,,$ (4)

where $T^{(\rm mass)}_{\mu\nu}$ is the effective energy-momentum tensor of the
graviton mass term that is defined by

$\displaystyle T^{(\rm mass)}_{\mu\nu}$ $\displaystyle=m^{2}M_{\rm
Pl}^{2}\Big{(}{1\over
2}g_{\mu\nu}U-U_{[Z]}{\bar{f}}_{\mu\nu}-2U_{[Z^{2}]}Z^{\rho}_{~{}(\mu}{\bar{f}}_{\nu)\rho}$
$\displaystyle-3U_{[Z^{3}]}Z^{\rho}_{~{}\sigma}Z^{\sigma}_{~{}(\mu}{\bar{f}}_{\nu)\rho}\Big{)}\,,$
(5)

and $T^{(\rm m)}_{\mu\nu}$ is an energy-momentum tensor for the matter field.
We assume a non-relativistic perfect fluid for the matter field,

$\displaystyle T^{(\rm m)}_{\mu\nu}=\rho u_{\mu}u_{\nu}\,,$ (6)

where $\rho$ is the energy density and $u^{\mu}$ is a four-velocity of the
matter field. Requiring the conservation of the energy-momentum tensor,

$\displaystyle\nabla^{\mu}T^{(\rm m)}_{\mu\nu}=0\,,$ (7)

the divergence of the Einstein equation (4) gives

$\displaystyle\nabla^{\mu}T^{(\rm mass)}_{\mu\nu}=0\,.$ (8)

This conservation equation of the effective energy-momentum tensor of the
graviton mass provides the equation of motion of the Stückelberg fields, which
can be obtained by varying the action with respect to $\phi^{a}$.

##  III Cosmological perturbation

In this section, we revisit linear perturbations around a cosmological
background in PMG. We introduce the zero-curvature scaling limit for simple
calculations.

###  III.1 zero-curvature limit

In the following, we consider a homogeneous and isotropic space-time with
positive spatial curvature, i.e., an open ($\kappa>0$) FLRW, given by

$\displaystyle g_{\mu\nu}dx^{\mu}dx^{\nu}$
$\displaystyle=-dt^{2}+a(t)^{2}\Omega_{ij}dx^{i}dx^{j}\,,$ (9)

where

$\displaystyle\Omega_{ij}$
$\displaystyle=\delta_{ij}-\frac{\kappa\,x^{i}x^{j}}{1+\kappa\,x^{k}x^{k}}\,.$
(10)

The unique configuration for Stückelberg fields to preserve homogeneity and
isotropy is provided with Gumrukcuoglu  _et al._ (2011)

$\displaystyle{\phi}^{0}$
$\displaystyle=f(t)\,\sqrt{1+\kappa(x^{2}+y^{2}+z^{2})}\,,$
$\displaystyle{\phi}^{i}$ $\displaystyle=f(t)\,\sqrt{\kappa}\,x^{i}\,,$ (11)

and the ${\bar{f}}$-metric is then given by

$\displaystyle\bar{f}_{\mu\nu}dx^{\mu}dx^{\nu}$ $\displaystyle=\kappa
f^{2}\Omega_{ij}dx^{i}dx^{j}\,.$ (12)

Throughout this present paper, for simplicity, we use the zero-curvature
scaling limit, which is first introduced in de Rham  _et al._ (2014c). We
begin with a redefinition of the function $f(t)$ in (11) as

$\displaystyle f(t)=\frac{\alpha}{\sqrt{\kappa}}+\chi(t)\,,$ (13)

where $\alpha$ is a constant. Taking the $\kappa\to 0$ limit after the
redefinition (13), the line elements for the $g$-metric and the
$\bar{f}$-metric reduce to

$\displaystyle\lim_{\kappa\to 0}{g}_{\mu\nu}dx^{\mu}dx^{\nu}$
$\displaystyle=-dt^{2}+a(t)^{2}\delta_{ij}dx^{i}dx^{j}\,,$ (14)
$\displaystyle\lim_{\kappa\to 0}{\bar{f}}_{\mu\nu}dx^{\mu}dx^{\nu}$
$\displaystyle=\alpha^{2}\delta_{ij}dx^{i}dx^{j}\,.$ (15)

On the other hand, under the redefinition (13), the Lorentz-invariant quantity
$X$ becomes

$\displaystyle
X=-\frac{\alpha^{2}}{\kappa}-\frac{2\alpha}{\sqrt{\kappa}}\chi(t)-\chi(t)^{2}\,,$
(16)

and thus diverges in the limit of $\kappa\to 0$. To avoid this problem, we
redefine $X$ as

$\displaystyle\tilde{X}:=-\frac{\sqrt{\kappa}}{2\alpha}\left(X+\frac{\alpha^{2}}{\kappa}\right)\,,$
(17)

and then $\tilde{X}$ converges to a finite value since $\lim_{\kappa\to
0}\tilde{X}=\chi(t)$. The mass potential $U$ is an arbitrary function of $X$
and the trace of $Z^{n}$, and we can thus redefine it as
$\tilde{U}(\tilde{X},[Z],[Z^{2}],[Z^{3}]):=U(X,[Z],[Z^{2}],[Z^{3}])$.
Therefore, introducing the redefinitions of (13) and (17), the action (3) as
well as the equations of motion (4) still remain finite in the zero-curvature
limit, and the equations of motion for background fields and perturbations in
this limit should agree with ones where the spatial curvature can be ignored.

###  III.2 Background and linear perturbations

In the present paper, we adopt the Newtonian gauge whose metric is given by

$\displaystyle g_{\mu\nu}dx^{\mu}dx^{\nu}=$
$\displaystyle-\left(1+2\Phi(t,\bf{x})\right)dt^{2}$
$\displaystyle+a(t)^{2}\left(1+2\Psi(t,\bf{x})\right)\delta_{ij}d{\bf
x}^{i}d{\bf x}^{j}\,.$ (18)

We define the Stückelberg perturbations through the coordinate transformation,

$\displaystyle x^{\mu}\to x^{\mu}+\delta^{\mu}_{a}\Pi^{a}\,,$ (19)

and then the Stückelberg fields become

$\displaystyle\phi^{0}$
$\displaystyle=\frac{\alpha}{\sqrt{\kappa}}+\chi\left(t+\Pi^{0}\right)+\mathcal{O}(\sqrt{\kappa})\,,$
(20) $\displaystyle\phi^{i}$
$\displaystyle=\alpha\left(x^{i}+\partial^{i}\Pi\right)+\mathcal{O}(\sqrt{\kappa})\,.$
(21)

Consequently, we obtain the perturbed $\bar{f}$-metric in the zero-curvature
limit, which is given by

$\displaystyle\lim_{\kappa\to 0}\bar{f}_{\mu\nu}dx^{\mu}dx^{\nu}$
$\displaystyle~{}=\alpha^{2}\partial^{i}\dot{\Pi}\partial_{i}\dot{\Pi}dt^{2}+2\alpha^{2}(\partial_{i}\dot{\Pi}+\partial^{k}\dot{\Pi}\partial_{k}\partial_{i}\Pi)dtdx^{i}$
$\displaystyle~{}~{}~{}+\alpha^{2}(\delta_{ij}+2\partial_{i}\partial_{j}\Pi+\partial_{i}\partial^{k}\Pi\partial_{k}\partial_{j}\Pi)dx^{i}dx^{j}\,.$
(22)

The matter perturbation is provided by $\rho(t,{\bf
x})=\rho(t)+\delta\rho(t,{\bf x})$, which gives

$\displaystyle T^{({\rm m})0}{}_{~{}0}=-(\rho+\delta\rho)\,,$ (23)
$\displaystyle T^{({\rm m})0}{}_{~{}i}=T^{({\rm m})i}{}_{~{}0}=T^{({\rm
m})i}{}_{~{}j}=0\,.$ (24)

###  III.3 Background and linear equations

By substituting (14) and (15) into (4) and (7), we obtain the background
equations,

$\displaystyle H^{2}-\frac{\rho+\rho_{g}}{3M_{\rm Pl}^{2}}$
$\displaystyle=0\,,$ (25)
$\displaystyle\dot{H}+\frac{\rho+\rho_{g}+p_{g}}{2M_{\rm Pl}^{2}}$
$\displaystyle=0\,,$ (26) $\displaystyle\dot{\rho}+3H\rho$
$\displaystyle=0\,,$ (27)

where we have defined the effective energy density and pressure of the
graviton mass term as

$\displaystyle\rho_{g}$ $\displaystyle:=-{1\over 2}m^{2}M_{\rm Pl}^{2}U\,,$
(28) $\displaystyle p_{g}$ $\displaystyle:={1\over 2}m^{2}M_{\rm
Pl}^{2}\Big{(}U-2\xi^{2}U_{[Z]}-4\xi^{4}U_{[Z^{2}]}-6\xi^{6}U_{[Z^{3}]}\Big{)}\,,$
(29)

and $\xi:=\alpha/a$. Similarly, from (8), we obtain the background Stückelberg
field equation,

$\displaystyle\dot{\rho}_{g}+3H(\rho_{g}+p_{g})$ $\displaystyle=0\,.$ (30)

At the linear level, $(0,0)$ and traceless components of the Einstein equation
(4) give

$\displaystyle\frac{2}{a^{2}}\partial^{2}\Psi+6H\left(H\Phi-\dot{\Psi}\right)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle~{}~{}~{}~{}~{}+\frac{(\rho_{g}+p_{g})\left(\partial^{2}\Pi-3\Psi\right)}{M_{\rm
Pl}^{2}}+\frac{\delta\rho}{M_{\rm Pl}^{2}}$ $\displaystyle=0\,,$ (31)
$\displaystyle\partial^{2}(\Phi+\Psi+a^{2}M_{\rm GW}^{2}\Pi)$
$\displaystyle=0\,,$ (32)

where we define the tensor graviton mass as

$\displaystyle M_{\rm GW}^{2}={2\over M_{\rm
Pl}^{2}}\Big{[}\rho_{g}+p_{g}-2m^{2}M_{\rm
Pl}^{2}\xi^{4}\left(U_{[Z^{2}]}+3\xi^{2}U_{[Z^{3}]}\right)\Big{]}\,.$ (33)

The time and space components of the linearized Stückelberg equations are
given by

$\displaystyle{\cal C}(-3\Psi+\partial^{2}\Pi+3H\Pi^{0})$ $\displaystyle=0\,,$
(34) $\displaystyle\partial_{i}\Big{[}2HM_{\rm Pl}^{2}M_{\rm
GW}^{2}\partial^{2}\Pi-\dot{p}_{g}(\partial^{2}\Pi-3\Psi)~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle+3H(\rho_{g}+p_{g})\Phi-{\cal
C}(\partial^{2}\Pi+3H\Pi^{0}-3\Psi)$
$\displaystyle+a^{2}H(2\dot{\rho}_{g}-3\dot{p}_{g})\dot{\Pi}-3a^{2}H(\rho_{g}+p_{g})\ddot{\Pi}\Big{]}$
$\displaystyle=0\,,$ (35)

where

$\displaystyle{\cal C}$ $\displaystyle:=m^{2}M_{\rm
Pl}^{2}\dot{\chi}^{2}U_{\chi\chi}\,.$ (36)

Integrating out $\Pi^{0}$ by using (34), the equation for $\Pi$, (35), becomes

$\displaystyle 2HM_{\rm Pl}^{2}M_{\rm
GW}^{2}\partial^{2}\Pi-\dot{p}_{g}(\partial^{2}\Pi-3\Psi)+3H(\rho_{g}+p_{g})\Phi$
$\displaystyle~{}~{}+a^{2}H(2\dot{\rho}_{g}-3\dot{p}_{g})\dot{\Pi}-3a^{2}H(\rho_{g}+p_{g})\ddot{\Pi}=0\,.$
(37)

###  III.4 Subhorizon expansion

Hereafter, we consider the spherically symmetric perturbations: ${\cal
E}(t,\bm{x})\to{\cal E}(t,r)$, where ${\cal E}=\\{\Phi,\Psi,\Pi,\Pi^{0}\\}$
and $r=\sqrt{x^{2}+y^{2}+z^{2}}$, and we use the following notation

$\displaystyle\partial_{i}=\frac{x_{i}}{r}\frac{\partial}{\partial r}\,,$ (38)

for simplicity. Now, we would like to derive equations applicable for
subhorizon scales. To this end, we perform subhorizon expansion assuming that
the time scale of the evolution of perturbations is roughly the Hubble time,
i.e., $\dot{{\cal E}}\sim H{\cal E}$, and we also expand the perturbations as
${\cal E}={\cal E}^{(1)}+{\cal E}^{(2)}+\cdots$, where ${\cal E}^{(n)}={\cal
O}(\epsilon^{n})$ and $\epsilon$ is a small expansion parameter. In order to
extract the relevant terms in the subhorizon expansion, we additionally assume
Kenna-Allison  _et al._ (2020b)

$\displaystyle r={\cal O}(\epsilon^{1/2})\,,\quad\partial_{i}={\cal
O}(\epsilon^{-1/2})\,.$ (39)

This leads that the second spatial derivative of the perturbations at the
leading order are the same order of the background, that is,
$\partial^{2}{\cal E}^{(1)}={\cal O}(\epsilon^{0})$, and we therefore need to
impose $\delta\rho={\cal O}(\epsilon^{0})$ in order to reproduce the standard
Poisson equation at the linear level, that is, $\delta\rho\sim a^{-2}M_{\rm
Pl}^{2}\partial^{2}\Psi^{(1)}$.

Based on these assumptions, the linear equations (31), (32) and (37), of order
${\cal O}(\epsilon^{0})$ are respectively given by

$\displaystyle\frac{2}{a^{2}}\partial^{2}\Psi^{(1)}+\frac{\delta\rho}{M_{\rm
Pl}^{2}}+\frac{\rho_{g}+p_{g}}{M_{\rm Pl}^{2}}\partial^{2}\Pi^{(1)}$
$\displaystyle=0\,,$ (40)
$\displaystyle\partial^{2}\Psi^{(1)}+\partial^{2}\Phi^{(1)}+a^{2}M_{\rm
GW}^{2}\partial^{2}\Pi^{(1)}$ $\displaystyle=0\,,$ (41)
$\displaystyle\partial^{2}\Pi^{(1)}$ $\displaystyle=0\,,$ (42)

where we have eliminated the background contributions by using (25)-(27) and
(30).111The background equations can be safely used at this point since the
solutions for the perturbations ${\cal E}^{(1)}$ are sufficiently small
compared with the background as one can see from (43). Imposing that the
perturbations vanish at $r\to\infty$, the solution of these equations is given
by

$\displaystyle\Pi^{(1)}=\frac{Q(t)}{r}\,,\quad\Psi^{(1)}=\frac{G{\cal
M}}{ar}\,,\quad\Phi^{(1)}=-\frac{G{\cal M}}{ar}\,,$ (43)

where ${\cal M}$ is an enclosed mass defined by

$\displaystyle{\cal M}(t,r):=4\pi a^{3}\int_{0}^{r}drr^{2}\delta\rho(r)\,,$
(44)

and $Q$ is the integration constant. The gravitational potentials in (43) are
identical to ones in general relativity, and the fifth force due to the scalar
graviton is absent at the leading order in the subhorizon expansion.

Equations for the next order perturbations $\partial^{2}{\cal E}^{(2)}={\cal
O}(\epsilon^{1})$ are given by

$\displaystyle\frac{2}{a^{2}}\partial^{2}\Psi^{(2)}+6H(H\Phi^{(1)}-\dot{\Psi}^{(1)})$
$\displaystyle\qquad+\frac{\rho_{g}+p_{g}}{M_{\rm
Pl}^{2}}(\partial^{2}\Pi^{(2)}-3\Psi^{(1)})=0\,,$ (45)
$\displaystyle\partial^{2}\left(\Psi^{(2)}+\Phi^{(2)}+a^{2}M_{\rm
GW}^{2}\Pi^{(2)}\right)=0\,,$ (46) $\displaystyle 2HM_{\rm Pl}^{2}M_{\rm
GW}^{2}\partial^{2}\Pi^{(2)}-\dot{p}_{g}(\partial^{2}\Pi^{(2)}-3\Psi^{(1)})$
$\displaystyle\qquad+3H(\rho_{g}+p_{g})\Phi^{(1)}+a^{2}H(2\dot{\rho}_{g}-3\dot{p}_{g})\dot{\Pi}^{(1)}$
$\displaystyle\qquad-3a^{2}H(\rho_{g}+p_{g})\ddot{\Pi}^{(1)}=0\,.$ (47)

The solution of these equations can be obtained recursively. From (43), we get

$\displaystyle\partial^{2}\Pi^{(2)}=$
$\displaystyle\frac{3[H(\rho_{g}+p_{g})-\dot{p}_{g}]}{2HM_{\rm Pl}^{2}M_{\rm
GW}^{2}-\dot{p}_{g}}\frac{G{\cal M}}{ar}$
$\displaystyle+(\text{terms~{}with~{}}\dot{Q},\ddot{Q})\,,$ (48)
$\displaystyle\partial^{2}\Psi^{(2)}=$
$\displaystyle\frac{3a^{2}H}{2}\Bigg{[}\frac{\rho_{g}+p_{g}}{M_{\rm
Pl}^{2}}\frac{2M_{\rm Pl}^{2}M_{\rm GW}^{2}-(\rho_{g}+p_{g})}{2HM_{\rm
Pl}^{2}M_{\rm GW}^{2}-\dot{p}_{g}}+6H\Bigg{]}$
$\displaystyle\times\frac{G{\cal
M}}{ar}+(\text{terms~{}with~{}}\dot{Q},\ddot{Q})\,,$ (49)
$\displaystyle\partial^{2}\Phi^{(2)}=$
$\displaystyle-3a^{2}\Bigg{[}3H^{2}+M_{\rm
GW}^{2}\frac{H(\rho_{g}+p_{g})-\dot{p}_{g}}{2HM_{\rm Pl}^{2}M_{\rm
GW}^{2}-\dot{p}_{g}}$ $\displaystyle+\frac{H(\rho_{g}+p_{g})}{2M_{\rm
Pl}^{2}}\frac{2M_{\rm Pl}^{2}M_{\rm GW}^{2}-(\rho_{g}+p_{g})}{2HM_{\rm
Pl}^{2}M_{\rm GW}^{2}-\dot{p}_{g}}\Bigg{]}\frac{G{\cal M}}{ar}$
$\displaystyle+(\text{terms with }\dot{Q},\ddot{Q})\,.$ (50)

We separated $Q$-dependent terms, which are not explicitly written here. As we
will see later, the matching condition to the nonlinear solution leads to
$Q=0$, and thus $Q$-dependent terms vanish. Note that the solution for
$\Pi^{(2)}$ completely agrees with the result in Manita and Kimura (2022), and
one immediately notice that
$\Pi^{(2)}/\Pi^{(1)}\sim\Psi^{(2)}/\Psi^{(1)}\sim\Phi^{(2)}/\Phi^{(1)}\sim
r/(aH)^{-1}$, which is extremely small inside the horizon.

##  IV Nonlinear perturbations

In this section, we investigate the effect of the nonlinearity of the
additional scalar mode at small scales. The basic idea is similar to the case
of the Vainshtein mechanism since the mass potential in the action (3) yields
higher-order (spatial) derivative interactions such as galileon theories.
However, due to the presence of $[Z]$, $[Z^{2}]$, and $[Z^{3}]$ in an
arbitrary function of the mass potential, the standard perturbative analysis
leads to the infinite number of nonlinear terms responsible for the Vainshtein
screening. To evade this problem, we, hereafter, consider a general polynomial
up to quadratic order for $Z$ as the mass potential:

$\displaystyle U(X,[Z],[Z^{2}],[Z^{3}])$
$\displaystyle=a_{0}(X)+a_{1}(X)[Z]+a_{2}(X)[Z]^{2}+a_{3}(X)[Z^{2}]\,.$ (51)

###  IV.1 Leading order

In order to extract relevant nonlinear terms such as $(\partial^{2}{\cal
E})^{n}$, we again use the similar assumption introduced in Sec. III.4, that
is,

$\displaystyle{\cal E}=\\{\Phi,\Psi,\Pi,\Pi^{0}\\}={\cal O}(\epsilon),$
$\displaystyle\partial_{i}={\cal O}(\epsilon^{-1/2}),\quad\delta\rho={\cal
O}(1)\,,$ (52)

together with $\dot{{\cal E}}\sim H{\cal E}$ and the subhorizon limit,
$r\ll(aH)^{-1}$ without expanding ${\cal E}$ itself. Note that ${\cal E}$ used
hereafter in this subsection corresponds to ${\cal E}^{(1)}$ in the Sec.
III.4. Then, ${\cal O}(1)$ equations contain all the nonlinear terms with
$(\partial^{2}{\cal E})^{n}$ including the background and linear
contributions. Similarly to the linear case, we assume the background
equations (25)-(27) and (30) are satisfied, which can be easily verified after
deriving nonlinear solutions. These equations are shown explicitly in Appendix
A.

The Stückelberg field equation can be easily obtained by combining (67) and
(68) and eliminating the third spatial derivative $\Pi^{\prime\prime\prime}$,

$\displaystyle(1+\Pi^{\prime\prime}){\cal P}_{1}{\cal P}_{2}=0\,,$ (53)

where

$\displaystyle{\cal P}_{1}$ $\displaystyle={\cal C}_{1}+{\cal
C}_{2}\left[\left(\frac{\Pi^{\prime}}{r}\right)^{2}+2\frac{\Pi^{\prime}}{r}\right]+{\cal
C}_{3}\left[(\Pi^{\prime\prime})^{2}+2\Pi^{\prime\prime}\right]\,,$ (54)
$\displaystyle{\cal P}_{2}=~{}$ $\displaystyle{\cal
C}_{4}\left(\Pi^{\prime\prime}+2\frac{\Pi^{\prime}}{r}\right)+4{\cal
C}_{5}\left(\frac{\Pi^{\prime}}{r}\right)^{2}$
$\displaystyle-2\left(\frac{5{\cal C}_{4}}{4}-2{\cal
C}_{5}\right)(\Pi^{\prime\prime})^{2}$ $\displaystyle-\left(\frac{{\cal
C}_{4}}{4}-{\cal
C}_{5}\right)\left[\left(\frac{\Pi^{\prime}}{r}\right)^{4}+4\left(\frac{\Pi^{\prime}}{r}\right)^{3}\right]$
$\displaystyle-\left(\frac{3{\cal C}_{4}}{4}-{\cal
C}_{5}\right)\left[(\Pi^{\prime\prime})^{4}+4(\Pi^{\prime\prime})^{3}\right]$
$\displaystyle+\left(\frac{5{\cal C}_{4}}{2}-2{\cal
C}_{5}\right)\left[(\Pi^{\prime\prime})^{2}+2(\Pi^{\prime\prime})\right]\left[\left(\frac{\Pi^{\prime}}{r}\right)^{2}+2\frac{\Pi^{\prime}}{r}\right]\,.$
(55)

We have utilized the symbol ′ to represent the derivative with respect to $r$.
We have also defined ${\cal C}_{n}$ as

$\displaystyle{\cal C}_{1}$ $\displaystyle=6M_{\rm Pl}^{2}M_{\rm
GW}^{2}+\rho_{g}+9p_{g}\,,$ (56) $\displaystyle{\cal C}_{2}$
$\displaystyle=-3M_{\rm Pl}^{2}M_{\rm GW}^{2}+10\rho_{g}+18p_{g}\,,$ (57)
$\displaystyle{\cal C}_{3}$ $\displaystyle=3(3M_{\rm Pl}^{2}M_{\rm
GW}^{2}-4\rho_{g})\,,$ (58) $\displaystyle{\cal C}_{4}$
$\displaystyle=12\left[H(\rho_{g}+9p_{g})+3\dot{p}_{g}\right]\,,$ (59)
$\displaystyle{\cal C}_{5}$ $\displaystyle=H(12M_{\rm Pl}^{2}M_{\rm
GW}^{2}+5\rho_{g}+93p_{g})$ $\displaystyle~{}~{}+6M_{\rm Pl}^{2}M_{\rm
GW}\dot{M}_{\rm GW}+27\dot{p}_{g}\,.$ (60)

As boundary conditions, we impose that

  * •

the solution is regular at $r=0$

  * •

the solution obtained from (53) connects to the linear solution (43) at large
distances.

Note that the first condition is imposed since the Stückelberg equations (67)
and (68) have no source term. Manifestly, one of the solutions of (53) given
by $1+\Pi^{\prime\prime}=0$ does not satisfy the regularity condition, and we
therefore disregard this solution. To see other solutions, we first expand the
scalar graviton $\Pi^{\prime}$ around $r=0$ as

$\displaystyle\Pi^{\prime}=\sum_{n=0}^{\infty}\alpha_{n}r^{n}\,.$ (61)

In the case of ${\cal P}_{1}=0$, all coefficients except $\alpha_{1}$ are
zero, thus there is an exact solution proportional to $r$, implying that this
solution never matches the linear solution (43). In addition, the Stückelberg
field equation (68) never hold for this solution unless $6M_{\rm Pl}^{2}M_{\rm
GW}^{2}+\rho_{g}+9p_{g}=0$, which cannot be realized in a general cosmological
background. Therefore, this is not a physical solution, and we also disregard
this solution.

As for ${\cal P}_{2}=0$, only $\alpha_{1}$ is the non-vanishing coefficient,
and there are three branches of solutions:

$\displaystyle\Pi^{\prime}=\begin{cases}-r\,,\\\ -2r\,,\\\ 0\,,\end{cases}$
(62)

which are also the exact solutions of ${\cal P}_{2}=0$ as well as (67) and
(68). The first two solutions proportional to $r$ does not satisfy the
matching condition with the linear solution (43) at large scales, and we again
disregard these solutions. Therefore, $\Pi^{\prime}=0$ is the only solution
satisfying all the conditions, and the matching condition leads to the
vanishing of the integration constant in the linear solution of $\Pi^{\prime}$
given in (43), that is,

$\displaystyle Q=0\,.$ (63)

This indicates that the scalar graviton never appears even at small scales in
the leading order of the quasi-static approximation. Therefore, we proceed to
investigate the subleading order equations.

###  IV.2 Subleading order

In the previous subsection, we extracted the relevant nonlinear terms from the
equation of motion under the assumption (52), but this results in
$\Pi^{\prime}=0$ for all scales within the horizon. Thus, as in the case of
the linear perturbations, we need to further expand the perturbations as
${\cal E}={\cal E}^{(1)}+{\cal E}^{(2)}+\cdots$ with ${\cal E}^{(n)}={\cal
O}(\epsilon^{n})$. Note that the analysis done in the previous subsection is
for ${\cal E}^{(1)}$. Then, all the equations of motion for ${\cal E}^{(2)}$
completely coincide with ones in the linear perturbations, i.e., the solutions
are simply given by (48)-(50). In fact, the scalar graviton completely
decouples with the metric perturbations $\Psi$ and $\Phi$ in the decoupling
limit, and thus the gravitational potentials are exactly the same as ones in
Newtonian case, whose results are shown in the Appendix. B.

##  V Summary

In the present paper, we consider both linear and nonlinear perturbations in a
spherically symmetric configuration around the homogeneous and isotropic
background in PMG. In order to simplify the analysis, we utilized the zero-
curvature limit, which is applicable where the spatial curvature can be
ignored. Having in mind that nonlinear derivative interactions appearing in
the galileon-like theories are also present in PMG, we derived the nonlinear
equations for all perturbations and found the analytical solutions. The
condition for the nonlinear solution to be connected to the linear solution at
large scales requires the vanishing of the integration constant for the
leading solution of $\Pi$ defined in (43), that is, $Q=0$. This indicates that
the scalar graviton should be zero at the lowest order of the subhorizon
expansion. As we can easily see in (65) and (66), the modified Einstein
equations remains the same as in general relativity.

We have also confirmed that, even at the subleading order in the subhorizon
expansion, the nonlinear effect due to the scalar graviton is absent at all
scales, and the linear contributions dominate at all scales, which are highly
suppressed as can be seen in the Poisson equation. In fact, one of the PPN
parameters $\gamma:=-\Psi/\Phi$ can be roughly estimated as

$\displaystyle 1-\gamma$ $\displaystyle\sim\frac{a^{2}M_{\rm
GW}^{2}\Pi^{(2)}}{\Phi^{(1)}}\sim r^{2}/(aH)^{-2}\,,$ (64)

where we have assumed $M_{\rm GW}\sim H$ and $\rho_{g}\sim p_{g}\sim M_{\rm
Pl}^{2}H^{2}$. This implies that the correction of $\gamma$ to general
relativity is highly suppressed well inside the horizon. For example, at the
solar system scale, this is estimated as $1-\gamma\sim 10^{-29}$, and is much
smaller than the constraint by the Cassini spacecraft $1-\gamma\lesssim
10^{-5}$ Bertotti  _et al._ (2003).

In summary, the effect of modification due to the scalar graviton is strongly
suppressed in the subhorizon region. As we show in the Appendix B, this result
is consistent with the decoupling limit analysis, and the scalar mode
interacts with neither the tensor mode and external matter fields.

###### Acknowledgements.

We thank Tsutomu Kobayashi and Takahiro Tanaka for insightful comments. Y.M.
acknowledges the xPand package Pitrou  _et al._ (2013) that was used to
explicitly confirm calculations. This work is supported by the establishment
of university fellowships towards the creation of science technology
innovation (Y.M.), Japan Society for the Promotion of Science (JSPS) Overseas
Challenge Program for Young Researchers (Y.M.), JSPS Grants-in-Aid for
Scientific Research No. JP22K03605 (R.K.), and Waseda University Grant for
Special Research Project No. 2022C-632 (S.P.).

##  Appendix A Nonlinear equations at leading order

In this section, we show the leading-order nonlinear equations derived under
the assumption (52) in a spherically symmetric setup. With the order ansatz
(52), we extract all ${\cal O}(\epsilon^{0})$ terms from the Einstein and
Stückelberg equations, which are given by

$\displaystyle{\cal E}^{00}=$
$\displaystyle-\frac{2}{a^{2}}\partial^{2}\Psi-\frac{\delta\rho}{M_{\rm
Pl}^{2}}-\frac{3{\cal C}_{1}-2{\cal C}_{3}}{27M_{\rm Pl}^{2}}[\Pi]-\frac{{\cal
C}_{2}}{18M_{\rm Pl}^{2}}[\Pi]^{2}$ $\displaystyle-\frac{{\cal C}_{1}-{\cal
C}_{2}}{18M_{\rm Pl}^{2}}[\Pi^{2}]+\frac{3{\cal C}_{2}-2{\cal C}_{3}}{54M_{\rm
Pl}^{2}}\left([\Pi^{3}]+\frac{1}{4}[\Pi^{4}]\right)$
$\displaystyle-\frac{{\cal C}_{2}}{18M_{\rm
Pl}^{2}}\left([\Pi][\Pi^{2}]+\frac{1}{4}[\Pi^{2}]^{2}\right)\,,$ (65)
$\displaystyle{\cal E}^{\rm tl}=$
$\displaystyle\frac{1}{a^{2}}\partial^{2}(\Phi+\Psi)+\frac{6{\cal
C}_{1}-3{\cal C}_{2}-2{\cal C}_{3}}{27M_{\rm Pl}^{2}}[\Pi]$
$\displaystyle+\frac{{\cal C}_{1}-4{\cal C}_{2}+{\cal C}_{3}}{9M_{\rm
Pl}^{2}}[\Pi^{2}]+\frac{{\cal C}_{2}}{18M_{\rm
Pl}^{2}}(4[\Pi]^{2}+3[\Pi^{2}])$ $\displaystyle-\frac{2({\cal C}_{2}-2{\cal
C}_{3})}{27M_{\rm Pl}^{2}}\left([\Pi^{3}]+\frac{1}{4}[\Pi^{4}]\right)$
$\displaystyle+\frac{2{\cal C}_{2}}{9M_{\rm
Pl}^{2}}\left([\Pi][\Pi^{2}]+\frac{1}{4}[\Pi^{3}]\right)\,,$ (66)
$\displaystyle{\cal E}^{\rm St,0}=$ $\displaystyle-2{\cal C}_{4}[\Pi]+2(5{\cal
C}_{4}-6{\cal C}_{5})[\Pi^{2}]$ $\displaystyle-(5{\cal C}_{4}-4{\cal
C}_{5})\left([\Pi]^{2}+[\Pi^{2}][\Pi]+\frac{1}{4}[\Pi^{2}]^{2}\right)$
$\displaystyle+(11{\cal C}_{4}-12{\cal
C}_{5})\left([\Pi^{3}]+\frac{1}{4}[\Pi^{4}]\right)$
$\displaystyle+\frac{4}{3}{\cal E}^{{\rm
St},i}\dot{\Pi}_{i}(1+\hat{x}^{k}\hat{x}^{l}\Pi_{kl})^{-1}\,,$ (67)
$\displaystyle{\cal E}^{{\rm St},i}=$ $\displaystyle 6{\cal
C}_{1}\Pi^{ij}{}_{j}+3(2{\cal C}_{1}-{\cal C}_{2}+2{\cal
C}_{3})\Pi^{i}{}_{j}\Pi^{jk}{}_{k}$ $\displaystyle-3({\cal C}_{2}-2{\cal
C}_{3})\Pi^{jk}\Pi^{i}{}_{jk}+8{\cal C}_{3}\Pi^{ij}\Pi^{kl}\Pi_{jkl}$
$\displaystyle-(3{\cal C}_{2}-2{\cal
C}_{3})\Big{(}(\Pi^{2})_{jk}\Pi^{ijk}+\Pi^{ij}(\Pi^{2})^{kl}\Pi_{jkl}$
$\displaystyle+(\Pi^{3})^{i}{}_{j}\Pi^{jk}{}_{k}+(\Pi^{2})^{ij}\Pi^{kl}\Pi_{jkl}\Big{)}$
$\displaystyle-2(3{\cal C}_{2}-4{\cal
C}_{3})(\Pi^{2})^{i}{}_{j}\Pi^{jk}{}_{k}$ $\displaystyle+3{\cal
C}_{2}\Big{(}2[\Pi]\Pi^{ij}{}_{j}+2\Pi^{ij}[\Pi]\Pi_{jk}{}^{k}+[\Pi^{2}]\Pi^{ij}{}_{j}$
$\displaystyle+\Pi^{ij}[\Pi^{2}]\Pi^{k}{}_{kj}+2(\Pi^{2})^{ij}\Pi^{kl}\Pi_{jkl}\Big{)}\,,$
(68)

where $\hat{x}^{i}$ is defined by $\hat{x}^{i}:=x^{i}/r$,
$\Pi_{i}=\partial_{i}\Pi$, $\Pi_{ij}=\partial_{i}\partial_{j}\Pi$,
$(\Pi^{2})_{ij}=\Pi_{ik}\Pi^{k}_{~{}j}$, and so on. Note that the Stückelberg
equations, (67) and (68), does not contain the gravitational potentials $\Psi$
and $\Phi$ and the source term.

##  Appendix B Decoupling limit

In this appendix, we derive the effective action in the so-called decoupling
limit Luty  _et al._ (2003) to see the role of the scalar graviton. Hereafter,
we consider $G=1$.

In order to take the decoupling limit, we, for simplicity, utilize the zero-
curvature scaling limit defined through (13). Then the effective action in the
zero-curvature limit is provided by

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int{\rm
d}^{4}x\sqrt{-g}\frac{M_{\rm Pl}^{2}}{2}\Biggl{[}\,R+m^{2}U\big{(}{\cal
X},[{\cal Z}],[{\cal Z}^{2}],[{\cal Z}]^{2}\big{)}\Biggr{]}$ (69)
$\displaystyle+S_{\rm m}[g,\psi]\,,$

where ${\cal X}={\tilde{\phi}}^{a}{\tilde{\phi}}_{a}$ and ${\cal
Z}=g^{\mu\alpha}\delta_{ij}^{(3)}\partial_{\alpha}{\tilde{\phi}}^{i}\partial_{\nu}{\tilde{\phi}}^{j}$.
Here, the indices $i$ and $j$ are the spacial indices and $\delta_{ij}^{(3)}$
is the Kronecker delta in 3D space. In a flat FLRW metric, one can easily
revive the same background equations (25)-(27) with the Stückelberg
configuration, ${\tilde{\phi}}^{0}=\chi(t)$ and ${\tilde{\phi}}^{i}=x^{i}$,
from the action (69). In the following, we choose the mass potential as

$\displaystyle U=a_{1}({\cal X})[{\cal Z}]+a_{2}({\cal X})[{\cal
Z}]^{2}+a_{3}({\cal X})[{\cal Z}^{2}]\,.$ (70)

Now we decompose the Stückelberg field ${\tilde{\phi}}^{a}$ around unitary
gauge, ${\tilde{\phi}}^{a}=\delta^{a}_{~{}\mu}x^{\mu}$,

$\displaystyle{\tilde{\phi}}^{a}=x^{a}-V^{a}-\partial^{a}\varphi\,,$ (71)

into the vector mode $V^{a}$ and the scalar mode $\varphi^{a}$, and we perturb
the fluctuation of the metric around a Minkowski metric,
$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}/M_{\rm Pl}$. Then we take the decoupling
limit noting that

$\displaystyle m\to 0,\qquad M_{\rm Pl}^{2}\to\infty,$ (72)

keeping

(73)

The effective Lagrangian in the decoupling limit describes physics at
distances in the range $\Lambda_{2}^{-1}<r<m^{-1}$. Since $m\sim H_{0}$ in
order to explain the accelerated expansion of the universe at present Manita
and Kimura (2022), the decoupling limit approach overlaps the analysis using
the subhorizon approximation introduced in Sec. II and IV. Then the decoupling
limit Lagrangian ignoring the vector mode is given by

$\displaystyle{\cal L}_{\rm DL}$ $\displaystyle=$ $\displaystyle-{1\over
4}h^{\mu\nu}{\cal E}^{\alpha\beta}_{\mu\nu}h_{\alpha\beta}+h_{\mu\nu}T_{(\rm
m)}^{\mu\nu}$ (74)
$\displaystyle+a_{1}\,\Lambda_{2}^{4}\left[3+\delta_{ij}^{(3)}\left(\partial^{\mu}\partial_{i}\varphi\partial_{\mu}\partial_{j}\varphi-2\partial_{i}\partial_{j}\varphi\right)\right]$
$\displaystyle+F_{\rm
NL}(\Lambda_{2},\partial_{\mu}\varphi,\partial_{\mu}\partial_{i}\varphi)\,,$

where ${\cal E}^{\mu\nu\alpha\beta}$ is the linearized Einstein-Hilbert
kinetic operator defined as

$\displaystyle{\cal E}^{\mu\nu}_{~{}~{}\alpha\beta}$ $\displaystyle=$
$\displaystyle\left[\eta^{(\mu}_{~{}\alpha}\eta^{\nu)}_{~{}\beta}-\eta^{\mu\nu}\eta_{\alpha\beta}\right]\square-2\partial^{(\mu}\partial_{(\alpha}\eta^{\nu)}_{~{}\beta)}$
(75)
$\displaystyle+\partial^{\mu}\partial^{\nu}\eta_{\alpha\beta}+\partial_{\alpha}\partial_{\beta}\eta^{\mu\nu}\,.$

The decoupling limit Lagrangian does not contain higher-order time
derivatives, i.e., each scalar mode $\varphi$ has up to first derivative with
respect to time, and this clearly shows the absence of the BD ghost.
Differently from the dRGT case, the cutoff energy scale is
$\Lambda_{2}>\Lambda_{3}\equiv(m^{2}M_{\rm Pl})^{1/3}$. Furthermore, the
nonlinear interactions $F_{\rm NL}$ depend on only $\varphi$ and its
derivatives, which means the tensor mode and scalar mode are completely
decoupled in the decoupling limit. As a result, the matter field $T_{(\rm
m)}^{\mu\nu}$ simply couples to graviton only through $h_{\mu\nu}$, and the
gravitational potential is given by the Newtonian one. For this reason, the
vDVZ discontinuity is manifestly absent. This agrees with our results in the
analysis obtained in Sec. IV.

## References

  * de Rham and Gabadadze (2010) C. de Rham and G. Gabadadze, Phys. Rev. D 82, 044020 (2010), arXiv:1007.0443 [hep-th] .
  * de Rham  _et al._ (2011) C. de Rham, G. Gabadadze, and A. J. Tolley, Phys. Rev. Lett. 106, 231101 (2011), arXiv:1011.1232 [hep-th] .
  * Boulware and Deser (1972) D. G. Boulware and S. Deser, Phys. Rev. D 6, 3368 (1972).
  * De Felice  _et al._ (2012) A. De Felice, A. E. Gumrukcuoglu, and S. Mukohyama, Phys. Rev. Lett. 109, 171101 (2012), arXiv:1206.2080 [hep-th] .
  * Gumrukcuoglu  _et al._ (2011) A. E. Gumrukcuoglu, C. Lin, and S. Mukohyama, JCAP 11, 030 (2011), arXiv:1109.3845 [hep-th] .
  * Hinterbichler (2013) K. Hinterbichler, JHEP 10, 102 (2013), arXiv:1305.7227 [hep-th] .
  * Kimura and Yamauchi (2013) R. Kimura and D. Yamauchi, Phys. Rev. D 88, 084025 (2013), arXiv:1308.0523 [gr-qc] .
  * de Rham  _et al._ (2014a) C. de Rham, A. Matas, and A. J. Tolley, Class. Quant. Grav. 31, 165004 (2014a), arXiv:1311.6485 [hep-th] .
  * Gao (2014) X. Gao, Phys. Rev. D 90, 064024 (2014), arXiv:1403.6781 [hep-th] .
  * Folkerts  _et al._ (2011) S. Folkerts, A. Pritzel, and N. Wintergerst, (2011), arXiv:1107.3157 [hep-th] .
  * Huang  _et al._ (2013) Q.-G. Huang, K.-C. Zhang, and S.-Y. Zhou, JCAP 08, 050 (2013), arXiv:1306.4740 [hep-th] .
  * D’Amico  _et al._ (2013) G. D’Amico, G. Gabadadze, L. Hui, and D. Pirtskhalava, Phys. Rev. D 87, 064037 (2013), arXiv:1206.4253 [hep-th] .
  * De Felice and Mukohyama (2016) A. De Felice and S. Mukohyama, Phys. Lett. B 752, 302 (2016), arXiv:1506.01594 [hep-th] .
  * De Felice  _et al._ (2021) A. De Felice, F. Larrouturou, S. Mukohyama, and M. Oliosi, JCAP 04, 015 (2021), arXiv:2012.01073 [gr-qc] .
  * De Felice  _et al._ (2017a) A. De Felice, S. Mukohyama, and M. Oliosi, Phys. Rev. D 96, 024032 (2017a), arXiv:1701.01581 [hep-th] .
  * De Felice  _et al._ (2017b) A. De Felice, S. Mukohyama, and M. Oliosi, Phys. Rev. D 96, 104036 (2017b), arXiv:1709.03108 [hep-th] .
  * Hassan and Rosen (2012) S. F. Hassan and R. A. Rosen, JHEP 02, 126 (2012), arXiv:1109.3515 [hep-th] .
  * de Rham  _et al._ (2014b) C. de Rham, L. Keltner, and A. J. Tolley, Phys. Rev. D 90, 024050 (2014b), arXiv:1403.3690 [hep-th] .
  * de Rham  _et al._ (2014c) C. de Rham, M. Fasiello, and A. J. Tolley, Int. J. Mod. Phys. D 23, 1443006 (2014c), arXiv:1410.0960 [hep-th] .
  * Gumrukcuoglu  _et al._ (2021) A. E. Gumrukcuoglu, R. Kimura, M. Kenna-Allison, and K. Koyama, JCAP 09, 023 (2021), arXiv:2107.01423 [hep-th] .
  * Kenna-Allison  _et al._ (2020a) M. Kenna-Allison, A. E. Gümrükçüoglu, and K. Koyama, Phys. Rev. D 101, 084014 (2020a), arXiv:1912.08560 [hep-th] .
  * Gumrukcuoglu  _et al._ (2020) A. E. Gumrukcuoglu, R. Kimura, and K. Koyama, Phys. Rev. D 101, 124021 (2020), arXiv:2003.11831 [gr-qc] .
  * Kenna-Allison  _et al._ (2020b) M. Kenna-Allison, A. E. Gumrukcuoglu, and K. Koyama, Phys. Rev. D 102, 103524 (2020b), arXiv:2009.05405 [gr-qc] .
  * Will (2014) C. M. Will, Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Vainshtein (1972) A. I. Vainshtein, Phys. Lett. B 39, 393 (1972).
  * van Dam and Veltman (1970) H. van Dam and M. J. G. Veltman, Nucl. Phys. B 22, 397 (1970).
  * Zakharov (1970) V. I. Zakharov, JETP Lett. 12, 312 (1970).
  * Babichev  _et al._ (2009) E. Babichev, C. Deffayet, and R. Ziour, JHEP 05, 098 (2009), arXiv:0901.0393 [hep-th] .
  * Berezhiani  _et al._ (2013) L. Berezhiani, G. Chkareuli, and G. Gabadadze, Phys. Rev. D 88, 124020 (2013), arXiv:1302.0549 [hep-th] .
  * Manita and Kimura (2022) Y. Manita and R. Kimura, Phys. Rev. D 105, 084038 (2022), arXiv:2112.13855 [gr-qc] .
  * Bertotti  _et al._ (2003) B. Bertotti, L. Iess, and P. Tortora, Nature 425, 374 (2003).
  * Pitrou  _et al._ (2013) C. Pitrou, X. Roy, and O. Umeh, Class. Quant. Grav. 30, 165002 (2013), arXiv:1302.6174 [astro-ph.CO] .
  * Luty  _et al._ (2003) M. A. Luty, M. Porrati, and R. Rattazzi, JHEP 09, 029 (2003), arXiv:hep-th/0303116 .