# Optical response of the tightbinding model on the Fibonacci chain

Hiroki Iijima  Yuta Murakami  Akihisa Koga  Department of Physics, Tokyo
Institute of Technology, Meguro, Tokyo 152-8551, Japan

###### Abstract

We theoretically study the optical conductivity of the tightbinding model
which has two types of the hopping integrals arranged in the Fibonacci
sequence. Due to the lack of the translational symmetry, many peak structures
appear in the optical conductivity as well as the density of states. When the
ratio of two hopping integrals is large, the self-similar structure appears in
the optical conductivity. This implies that the optical response between the
high-energy bands is related to that within the low-energy bands, which should
originate from critical behavior in the wave functions. The effects of
disorders on the optical conductivity are also analyzed in order to show the
absence of the self-similarity in the tightbinding model with the random
sequence.

##  I INTRODUCTION

Quasicrystal has been attracting much interests since the first discovery of
the quasicrystalline phase in the Al-Mn alloy Shechtman  _et al._ (1984).
Experimental and theoretical efforts have been made to explore and understand
physical properties inherent in quasicrystals Elser and Henley (1985); Tsai
_et al._ (1987, 1988, 1990); Goldman and Kelton (1993); Tsai  _et al._ (2000);
Kamiya  _et al._ (2018); Iwasaki  _et al._ (2019); Kimura  _et al._ (1989);
Inaba (1996); Tamura  _et al._ (2021). Among them, the Au-Al-Yb alloy with
Tsai-type clusters Ishimasa  _et al._ (2011) exhibits interesting properties
at low temperatures. In the quasicrystal Au51Al34Yb15, quantum critical
behavior appears, while heavy fermion behavior appears in the approximant
Au51Al35Yb14 Deguchi  _et al._ (2012). This distinct behavior stimulates
theoretical investigations on the quasiperiodic structures in correlated
electron systems Takemori and Koga (2015); Takemura  _et al._ (2015); Watanabe
and Miyake (2013); Otsuki and Kusunose (2016); Andrade  _et al._ (2015); Koga
and Tsunetsugu (2017); Takemori  _et al._ (2020); Sakai and Koga (2021);
Ghadimi  _et al._ (2021). The optical response characteristic of the
quasicrystals have also been observed. In the Al-Cu-Fe quasicrystal, the
linear-$\omega$ dependence in the optical conductivity, which is distinct from
the conventional impurity scattering, has been observed Homes  _et al._
(1991). Furthermore, the direction dependent optical conductivity has been
reported in the Al-Co-Cu quasicrytal Basov  _et al._ (1994). These studies
suggest that the exotic optical properties arise from the quasiperiodic
structure. Although the low-frequency behavior has theoretically been examined
for quasiperiodic systems Aldea and Dulea (1988); Sánchez and Wang (2004);
Varma  _et al._ (2016); Tsunetsugu and Ueda (1991), the optical response at
finite frequencies has not been discussed in detail. An important point is
that the spatial features of the initial and final states play a crucial role
for the optical process. It is known that spatially extended states
characterized by the momenta are realized in the periodic system, while
critical states are, in general, realized in the quasiperiodic one Tsunetsugu
_et al._ (1986); Kohmoto and Banavar (1986); Fujiwara  _et al._ (1989);
Kohmoto  _et al._ (1987). Therefore, it is instructive to clarify the optical
response inherent in quasiperiodic systems by comparing them to systems with
distinct properties of eigenstates.

Motivated by this, we treat the tightbinding model which has two types of the
hopping integrals arranged in the Fibonacci sequence, as a simple model. It is
known that each eigenstate for the system shows critical behavior Kohmoto  _et
al._ (1987); Kohmoto and Banavar (1986); Fujiwara  _et al._ (1989); Kohmoto
_et al._ (1983); Macé  _et al._ (2017), which are characterized by
multifractal properties and the power law decay of the amplitude of the wave
functions in the real space. This eigenstate property is distinct from those
for the periodic systems where the spatially extended states are realized. We
then discuss the optical response in the tightbinding model on the Fibonacci
chain, examining the matrix elements of the current operator, which play a
crucial role for the optical conductivity.

The paper is organized as follows. In Sec. II, we introduce the tightbinding
model on the Fibonacci chain and derive the expression of the optical
conductivity in this system. In Sec. III, we discuss the optical response
inherent in the Fibonacci chain, comparing with that in the approximants. The
effect of the disorders is also addressed. A summary is given in the last
section.

##  II MODEL AND METHOD

We consider the tightbinding model to study the optical conductivity inherent
in the Fibonacci chain. The Hamiltonian is given as

$\displaystyle\hat{H}(t)$
$\displaystyle=-\sum_{n}\left(v_{n}e^{-iqL_{n}A(t)}\hat{c}_{n}^{{\dagger}}\hat{c}_{n+1}+\mathrm{H.c.}\right),$
(1)

where $\hat{c}_{n}(\hat{c}_{n}^{{\dagger}})$ is the annihilation (creation)
operator of a spinless fermion at the $n$th site, $v_{n}$ ($L_{n}$) denotes
the hopping integral (lattice spacing) between the $n$th and $(n+1)$th sites,
and $q$ is the charge of the fermion. The site-independent vector potential
$A(t)$ leads to the uniform electric field ${\cal E}(t)=-\partial_{t}A(t)$.

Here, we introduce the Fibonacci sequence ${\cal S}$ into the Hamiltonian. It
is known that the Fibonacci sequence is generated by means of the substitution
rule for two letters $L$ and $S$: $L\rightarrow LS$ and $S\rightarrow L$.
Applying the substitution rule to the initial sequence ${\cal S}_{1}=S$
iteratively, we obtain sequences $\\{L,LS,LSL,LSLLS,\cdots\\}$. The $i$th
sequence ${\cal S}_{i}$ is composed of $F_{i}$ letters, where $F_{i}$ is the
Fibonacci number. We deal with the tightbinding model with the total number of
sites $N=F_{i}$ and the hopping integral is given as $v_{n}=v_{L}(v_{S})$ if
the $n$th letter of the Fibonacci sequence ${\cal S}_{i}$ is $L$ ($S$).

In the paper, we study the linear optical response for the tightbinding model
on the Fibonacci chain. The optical conductivity is given as
$\sigma(\omega)=J(\omega)/{\cal E}(\omega)$, where $J(\omega)$ and ${\cal
E}(\omega)$ are the Fourier components of the current $J(t)$ and electric
field ${\cal E}(t)$. The current operator is split into two parts as
$\hat{J}(t)=\hat{j}_{1}+\hat{j}_{2}(t)$, with

$\displaystyle\hat{j}_{1}=$ $\displaystyle-
iq\sum_{n}v_{n}L_{n}(\hat{c}_{n}^{{\dagger}}\hat{c}_{n+1}-\hat{c}_{n+1}^{{\dagger}}\hat{c}_{n}),$
(2) $\displaystyle\hat{j}_{2}(t)=$
$\displaystyle-q^{2}\sum_{n}v_{n}L_{n}^{2}(\hat{c}_{n}^{{\dagger}}\hat{c}_{n+1}+\hat{c}_{n+1}^{{\dagger}}\hat{c}_{n})A(t).$
(3)

By means of the Kubo formula, the optical conductivity is expressed as

$\displaystyle\sigma(\omega)$
$\displaystyle=\frac{1}{N}\sum_{a,b}\frac{f(E_{b})-f(E_{a})}{i(\omega+i\delta)}\frac{|\bra{b}\hat{j}\ket{a}|^{2}}{\omega-(E_{b}-E_{a})+i\delta}$
$\displaystyle-\frac{1}{N}\sum_{a}\frac{q^{2}f(E_{a})}{i(\omega+i\delta)}\bra{a}\sum_{n}v_{n}L_{n}^{2}(\hat{c}^{{\dagger}}_{n}\hat{c}_{n+1}+\hat{c}^{{\dagger}}_{n+1}\hat{c}_{n})\ket{a},$
(4)

where $\ket{a}$ is a single electron eigenstate of the model with $A(t)=0$,
$E_{a}$ is the corresponding energy, $f(x)$ is the Fermi distribution
function, and $\delta$ is infinitesimal. The first term represents the optical
transition, and the other is the so-called Drude part. For the sake of
simplicity, we set $L_{n}=1$ in order to focus on the effects of the Fibonacci
structure in the hopping integrals Sánchez and Wang (2004); Aldea and Dulea
(1988).

In the following, we mainly consider the system with $N=F_{19}=4181$ under the
periodic boundary condition. This system is large enough to take the
quasiperiodic structure into account. In fact, we have confirmed that the
obtained results are essentially the same as those with $N=17711$. Here,
setting the Fermi energy at $E=0$ and $v_{L}$ as the unit of energy, we
discuss how the quasiperiodic structure affects the optical linear response.

##  III RESULT

We study the optical response in the tightbinding model with the Fibonacci
structure. First, we briefly discuss the one-particle states in the model
without the external electric field.

Figure 1: Density of states in the tightbinding model on the Fibonacci chain
(a) and the approximant (b) with $v_{S}=2.5$ and $A(t)=0$. Inset in (b) shows
the dispersion relation of the model for the approximant. The vertical axes
are in an arbitrary unit.

Figure 1(a) shows the density of states (DOS) for the tightbinding model on
the Fibonacci chain with $v_{S}=2.5$. Many delta-function peaks appear in the
DOS, which is known to be purely singular continuous like the Cantor set
Kohmoto  _et al._ (1983); Tang and Kohmoto (1986); Kohmoto  _et al._ (1987);
Sütő (1989). To plot a delta-function peak in DOS, we practically use the
Gaussian with a small width in Figure 1. We also use the same strategy to plot
delta-function peaks in the following. The energy levels are roughly
classified into lower (L), middle (M), and upper (U) bands [see Fig. 1]. When
one focuses on the M band, there exist three smaller bands, suggesting the
nested structure in the DOS Ninomiya (1986). In fact, the DOS scaled by $R$ is
in a good agreement with the original one, where $R=E^{U}_{max}/E^{M}_{max}$
and $E^{\alpha}_{max}$ is the maximum energy in the $\alpha(=U,M)$ band. It is
also known that single-particle states of this model are critical and the
decay for each wave function is slower than exponential one Kohmoto  _et al._
(1987); Fujiwara  _et al._ (1989); Kohmoto and Banavar (1986); Kohmoto  _et
al._ (1983); Macé  _et al._ (2017). In contrast to the Fibonacci case, the
smooth DOS should appear in the periodic system. For comparison, we consider
the tightbinding model on the approximant, where the shorter Fibonacci
sequence ${\cal S}_{4}(=LSL)$ is periodically arranged. This is the minimal
approximant with three bands in common with the Fibonacci chain, as shown in
Fig. 1(b). Since each single-particle state in the approximant is
characterized by the momentum, its wave function is spatially extended, in
contrast to the Fibonacci case. When one considers the approximant with the
longer sequence, the corresponding DOS approaches the Fibonacci one and
critical behavior should appear in the wave function.

Figure 2: (Color online) Real part of the optical conductivity for the
Fibonacci chain (a) and approximant (b) with $v_{S}=2.5$. The dashed line in
(a) shows the optical conductivity scaled by $R$. The areas shaded in yellow,
green, and blue represent different types of excitations, see the text for
details. The vertical axes are in an arbitrary unit.

Now we discuss the optical conductivity in these tightbinding models with
distinct behaviors in the one-particle states. Figure 2 shows the real part of
the optical conductivities in the systems with $v_{S}=2.5$. It is clarified
that several peaks appear in the Fibonacci case, while few structures appear
in the approximant. The peaks in the former can be explained by the energy
spectrum in the DOS, where there exist L, M, and U bands [see Fig. 1(a)]. When
the half-filled model is considered, the optical response is categorized into
three types of excitations; those within the M band, between the M band and
the other bands, and between the L and U bands. In fact, these excitations are
clearly separated, which are shown as the yellow, green, and blue areas in
Fig. 2(a). The areas shaded in yellow, green, and blue represent the
excitations within the M band, between the M band and the other bands, and
between the L band and the U band, respectively. Furthermore, a large
intensity in the optical conductivity appears, which corresponds to the
optical transition between the energy levels with the large DOS. By contrast,
a rather simple structure appears in the approximant with the similar DOS, as
shown in Fig. 2(b). It is clarified that no peak structure appears in the low
frequency region $0<\omega<1$ and in higher frequency region, peaks are not
widely distributed, in contrast to the Fibonacci case. Furthermore, the gap
structure in the optical conductivity is not directly related to the DOS.
Thus, the DOS is not enough to explain the optical response for the
approximants. This originates from the existence of the translational symmetry
in the approximant, which leads to the selection rule for the optical
response. Namely, the optical transition between distinct wave numbers is
forbidden [see Fig. 2(b)]. Only the interband excitations contribute to the
optical transitions with finite frequencies. Namely, the transition within the
M band is absent, although there is a signal from the second term of eq. (II)
at $\omega=0$ , i.e. the Durde peak.

Figure 3: Real part of the optical conductivity as a function of $v_{S}$ for
the Fibonacci chain (a) and the approximant (b). The color scale is in an
arbitrary unit.

The $v_{S}$ dependent conductivity is shown in Fig. 3. When $v_{S}\lesssim 2$,
the energy ranges for intra- and interband excitations overlap and their
intensities are rather small. Therefore, it is hard to see the characteristic
feature of the Fibonacci chain. By contrast, the large $v_{S}$ leads to
different behavior in the optical conductivity. In the case, two interband
excitations are almost proportional to $2v_{S}$ and $v_{S}$, which allows us
to distinguish these excitations in the optical conductivity. With the large
$v_{S}$, the intensities of the optical conductivity are clear both in
Fibonacci and approximant cases. On the other hand, qualitatively distinct
behavior appears in each system: in the Fibonacci case, several peak
structures in the intra- and interband exciations appear even when $v_{S}$ is
large.

Figure 4: (a) Matrix elements of the current operator
$|\braket{b}{\hat{j}}{a}|$ for the Fibonacci chain. (b) Magnified figures of
the rectangular part in (a).

In order to obtain further insight into many spiky behavior in the optical
conductivity for the Fibonacci chain, we focus on the matrix elements of the
current operator $\bra{b}\hat{j}\ket{a}$ in eq. (II), where the indices $a$
and $b$ are listed in ascending order by the eigenenergies. We note that no
degeneracy appears in the energy levels and the matrix elements are uniquely
determined. Figure 4 shows the absolute value of the matrix elements of the
current operator for the Fibonacci chain. We find the dense structure in the
matrix elements although some of them is invisible in the figure. This implies
that the optical transition between an arbitrary pair of states is allowed.
This should be consistent with the fact that there exists no translational
symmetry and each eigenstate in the tightbinding model is critical. By this
reason, such optical properties realized in the quasiperiodic systems do not
depend on the details of band structure. In the tightbinding model on the
Fibonacci sequence ${\cal S}_{n}$, the M band is composed of the $F_{n-3}$
energy levels, and the others are composed of the $F_{n-2}$ levels Zheng
(1987). The matrix elements should be classified into nine groups, and each
has the dense structure. One of the most interesting points is that self-
similar behavior appears. When one focuses on the matrix elements for the M
band, which is marked by the green square in Fig. 4(a), its pattern is similar
to the original one. This self-similar pattern enriches the peak structure in
the optical conductivity. Namely, some intraband excitations in the M band are
related to the interband excitation between L and U bands. In the sense, the
self-similar behavior in the matrix elements are reflected in the optical
conductivity. In fact, the optical conductivity scaled by $R$ is in a good
agreement with the original one, as shown in Fig. 2(a). This is in contrast to
the conventional periodic system, where the optical transition is restricted
between the states with the same momentum (optical selection rule).

Up to now, we have discussed the optical responses for the tightbinding models
on the Fibonacci and its approximant, and have clarified that the structure in
the current operator leads to spiky and self-similar behavior in the optical
conductivity for the Fibonacci system. It is naively expected that disordered
systems should have a dense structure in the matrix elements. Therefore, it
should be instructive to clarify the effect of the disorders for the optical
response in the Fibonacci system. In this study, the bond disorder for the
Fibonacci sequence is introduced. A disorder is created by picking a random
site $i$ and swapping its connecting hoppings $v_{i}$ and $v_{i-1}$;
$v_{i}\leftrightarrow v_{i-1}$. We note that the tightbinding model with large
disorders is not reduced to the random-hopping model since the numbers of
$v_{L}$ and $v_{S}$ are conserved. Here, we clarify how the self-similar
structure in the optical response is affected by the introduction of the
disorders.

Figure 5: (a) DOS of the disordered system with $v_{S}=2.5$ when $d=0.01$ and
$d=0.1$. The real part of the optical conductivity for $v_{S}=2.5$ when
(b)$d=0.1$, (c)$d=0.01$. The areas shaded in grey depict the optical
conductivity of the Fibonacci chain as shown in Fig. 2(a). The vertical axes
are in an arbitrary unit.

Figure 5(a) shows the DOS for the Fibonacci chains with $d=0.01$, and $0.1$,
where $d(=m/N)$ is the density of the disorders and $m$ is the number of swap
operations. The results are evaluated by means of, at least, a thousand
independent random samples and standard deviations are invisible in the
figure. When the disorders are introduced in the one-dimensional chain, it may
be regarded as the system composed of Fibonacci chains with finite length. In
the case, the average of the length is roughly given by $1/d$ and thereby the
characteristic feature of the Fibonacci chain still remains when $d<0.01$.
Increasing the disorders, the singular continuous spectrum inherent in the
Fibonacci chain smears and the spectrum have the intensity in the broader
energy range, as shown in Fig. 5(a). In Fig. 5(b)(c), we show the optical
conductivity in the system with $d=0.01$, and $0.1$. We find that, increasing
$d$, the peak structures which is clear in the Fibonacci case smear eg. at
$\omega\approx 0.6,2.2,2.7,4.5,6.0$. Some peaks instead emerge with increasing
$d$, eg. at $\omega\approx 1.7,4.0$.

Figure 6: The matrix elements of the current operator
$|\braket{b}{\hat{j}}{a}|$ for the disordered Fibonacci chain with $v_{S}=2.5$
and $d=0.1$.

This makes the self-similar structure in the conductivity, which is clarified
by rescaling the frequency, less likely with increasing $d$. This is
consistent with the fact that self-similar pattern smears in the matrix
elements of the current operator, as shown in Fig. 6.

##  IV CONCLUSION

We have investigated the optical response for the one-dimensional tightbinding
model on the Fibonacci chain. We have calculated the optical conductivity in
terms of the Kubo formula and have clarified that it reflects the singular
continuous DOS inherent in the Fibonacci chain. This is contrast to those for
the approximants, where the optical transition is restricted due to the
existence of the translational symmetry. By examining the matrix elements of
the current operator carefully, we have found the self-similar structure in
the optical conductivity, where it is well scaled by the value $R$. We have
also discussed the effects of disorders in the Fibonacci chain and clarified
that the self-similar structures in optical conductivity smear. From these
analyses, we conclude that the self-similar feature observed in the optical
conductivity as well as the DOS originates from the quasiperiodic structure
inherent in the Fibonacci sequence. Since the Fibonacci lattice can be
realized for cold atoms on the optical lattice Singh  _et al._ (2015) and for
photons in a photonic waveguide array Verbin  _et al._ (2015), these systems
are potential playgrounds for our theoretical prediction. In particular, for
the cold atom systems, one can effectively realize the electric field Jotzu
_et al._ (2014). Our results suggest that physical observables such as optical
conductivity can show peculiar structures originating from quasiperiodic
lattice structures. This result shall stimulate further experimental and
theoretical studies to explore characteristic self-similar structures in
optical conductivity in various quasiperiodic systems. To undertand the origin
and the generality of such structures, detailed analyses of wave functions
and/or various types of quasi-periodic systems are required. Furthermore,
nonlinear optical responses have been actively studied in the context of
condensed matter physics. In future, it is also important to clarify how
quasiperiodic structures affect nonlinear optical properties.

###### Acknowledgements.

Parts of the numerical calculations are performed in the supercomputing
systems in ISSP, the University of Tokyo. This work is supported by Grant-in-
Aid for Scientific Research from JSPS, KAKENHI Grant Nos. JP20K14412,
JP20H05265, JP21H05017 (Y.M.) and JP19H05821, JP18K04678, JP17K05536 (A.K.),
and JST CREST Grant No. JPMJCR1901 (Y.M.).

## References

  * Shechtman  _et al._ (1984) D. Shechtman, I. Blech, D. Gratias, and J. W. Cahn, Phys. Rev. Lett. 53, 1951 (1984).
  * Elser and Henley (1985) V. Elser and C. L. Henley, Phys. Rev. Lett. 55, 2883 (1985).
  * Tsai  _et al._ (1987) A.-P. Tsai, A. Inoue, and T. Masumoto, Japanese Journal of Applied Physics 26, L1505 (1987).
  * Tsai  _et al._ (1988) A.-P. Tsai, A. Inoue, and T. Masumoto, Japanese Journal of Applied Physics 27, L1587 (1988).
  * Tsai  _et al._ (1990) A. P. Tsai, A. Inoue, Y. Yokoyama, and T. Masumoto, Materials Transactions, JIM 31, 98 (1990).
  * Goldman and Kelton (1993) A. I. Goldman and R. F. Kelton, Rev. Mod. Phys. 65, 213 (1993).
  * Tsai  _et al._ (2000) A. P. Tsai, J. Q. Guo, E. Abe, H. Takakura, and T. J. Sato, Nature 408, 537 (2000).
  * Kamiya  _et al._ (2018) K. Kamiya, T. Takeuchi, N. Kabeya, N. Wada, T. Ishimasa, A. Ochiai, K. Deguchi, K. Imura, and N. K. Sato, Nature Communications 9, 154 (2018).
  * Iwasaki  _et al._ (2019) Y. Iwasaki, K. Kitahara, and K. Kimura, Phys. Rev. Materials 3, 061601 (2019).
  * Kimura  _et al._ (1989) K. Kimura, H. Iwahashi, T. Hashimoto, S. Takeuchi, U. Mizutani, S. Ohashi, and G. Itoh, Journal of the Physical Society of Japan 58, 2472 (1989).
  * Inaba (1996) A. Inaba, Philosophical Magazine Letters 74, 381 (1996).
  * Tamura  _et al._ (2021) R. Tamura, A. Ishikawa, S. Suzuki, T. Kotajima, Y. Tanaka, T. Seki, N. Shibata, T. Yamada, T. Fujii, C.-W. Wang, M. Avdeev, K. Nawa, D. Okuyama, and T. J. Sato, Journal of the American Chemical Society 143, 19938 (2021).
  * Ishimasa  _et al._ (2011) T. Ishimasa, Y. Tanaka, and S. Kashimoto, Philosophical Magazine 91, 4218 (2011).
  * Deguchi  _et al._ (2012) K. Deguchi, S. Matsukawa, N. K. Sato, T. Hattori, K. Ishida, H. Takakura, and T. Ishimasa, Nature Materials 11, 1013 (2012).
  * Takemori and Koga (2015) N. Takemori and A. Koga, J. Phys. Soc. Jpn. 84, 023701 (2015).
  * Takemura  _et al._ (2015) S. Takemura, N. Takemori, and A. Koga, Phys. Rev. B 91, 165114 (2015).
  * Watanabe and Miyake (2013) S. Watanabe and K. Miyake, J. Phys. Soc. Jpn. 82, 083704 (2013).
  * Otsuki and Kusunose (2016) J. Otsuki and H. Kusunose, J. Phys. Soc. Jpn. 85, 073712 (2016).
  * Andrade  _et al._ (2015) E. C. Andrade, A. Jagannathan, E. Miranda, M. Vojta, and V. Dobrosavljević, Phys. Rev. Lett. 115, 036403 (2015).
  * Koga and Tsunetsugu (2017) A. Koga and H. Tsunetsugu, Phys. Rev. B 96, 214402 (2017).
  * Takemori  _et al._ (2020) N. Takemori, R. Arita, and S. Sakai, Phys. Rev. B 102, 115108 (2020).
  * Sakai and Koga (2021) S. Sakai and A. Koga, Mat. Trans. 62, 380 (2021).
  * Ghadimi  _et al._ (2021) R. Ghadimi, T. Sugimoto, K. Tanaka, and T. Tohyama, Phys. Rev. B 104, 144511 (2021).
  * Homes  _et al._ (1991) C. C. Homes, T. Timusk, X. Wu, Z. Altounian, A. Sahnoune, and J. O. Ström-Olsen, Phys. Rev. Lett. 67, 2694 (1991).
  * Basov  _et al._ (1994) D. N. Basov, T. Timusk, F. Barakat, J. Greedan, and B. Grushko, Phys. Rev. Lett. 72, 1937 (1994).
  * Aldea and Dulea (1988) A. Aldea and M. Dulea, Phys. Rev. Lett. 60, 1672 (1988).
  * Sánchez and Wang (2004) V. Sánchez and C. Wang, Phys. Rev. B 70, 144207 (2004).
  * Varma  _et al._ (2016) V. K. Varma, S. Pilati, and V. E. Kravtsov, Phys. Rev. B 94, 214204 (2016).
  * Tsunetsugu and Ueda (1991) H. Tsunetsugu and K. Ueda, Phys. Rev. B 43, 8892 (1991).
  * Tsunetsugu  _et al._ (1986) H. Tsunetsugu, T. Fujiwara, K. Ueda, and T. Tokihiro, Journal of the Physical Society of Japan 55, 1420 (1986).
  * Kohmoto and Banavar (1986) M. Kohmoto and J. R. Banavar, Phys. Rev. B 34, 563 (1986).
  * Fujiwara  _et al._ (1989) T. Fujiwara, M. Kohmoto, and T. Tokihiro, Phys. Rev. B 40, 7413 (1989).
  * Kohmoto  _et al._ (1987) M. Kohmoto, B. Sutherland, and C. Tang, Phys. Rev. B 35, 1020 (1987).
  * Kohmoto  _et al._ (1983) M. Kohmoto, L. P. Kadanoff, and C. Tang, Phys. Rev. Lett. 50, 1870 (1983).
  * Macé  _et al._ (2017) N. Macé, A. Jagannathan, P. Kalugin, R. Mosseri, and F. Piéchon, Phys. Rev. B 96, 045138 (2017).
  * Tang and Kohmoto (1986) C. Tang and M. Kohmoto, Phys. Rev. B 34, 2041 (1986).
  * Sütő (1989) A. Sütő, Journal of Statistical Physics 56, 525 (1989).
  * Ninomiya (1986) T. Ninomiya, J. Phys. Soc. Jpn. 55, 3709 (1986).
  * Zheng (1987) W. M. Zheng, Phys. Rev. A 35, 1467 (1987).
  * Singh  _et al._ (2015) K. Singh, K. Saha, S. A. Parameswaran, and D. M. Weld, Phys. Rev. A 92, 063426 (2015).
  * Verbin  _et al._ (2015) M. Verbin, O. Zilberberg, Y. Lahini, Y. E. Kraus, and Y. Silberberg, Phys. Rev. B 91, 064201 (2015).
  * Jotzu  _et al._ (2014) G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Nature 515, 237 (2014).


# Thin-shell wormholes in (2+1)-dimensional $F(R)$ theories

Cecilia Bejarano, Ernesto F. Eiroa, Griselda Figueroa-Aguirre
Instituto de Astronomía y Física del Espacio (IAFE, CONICET-UBA),
Casilla de Correo 67, Sucursal 28, 1428, Buenos Aires, Argentina e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We construct a broad family of thin-shell wormholes with circular symmetry in
(2+1)-dimensional $F(R)$ theories of gravity, with constant scalar curvature
$R$. We study the stability of the static configurations under perturbations
preserving the symmetry. We present examples of charged thin-shell wormholes
which are asymptotically anti-de Sitter at both sides of the throat. We show
that stable solutions are possible when suitable values of the parameters are
taken.

##  1 Introduction

Traversable Lorentzian wormholes [1] have been widely studied in the context
of General Relativity (GR) and also in alternative models of gravity. They can
represent shortcuts in the same spacetime or connect two different universes.
A quite interesting feature is the fact that they are singularity-free
gravitational objects; moreover, wormhole throats can replace black hole
singularities in solutions within the framework of some modified gravity
theories (see [2] and references therein).

It is well established that wormhole spacetimes within GR always require the
presence of exotic matter –which violates at least one of the energy
conditions– somewhere. This is not necessarily the case outside of GR, so it
is worth examining the possible violation of the energy conditions by
wormholes in different theories of gravity. The so-called thin-shell wormholes
have drawn the attention of the gravitational community, since the exotic
matter can be reduced to a minimal amount and be highly confined. In addition
to this, the junction procedure also allows to analyze the wormhole dynamics
by focusing on the thin shell where the two geometries are joined [3, 4].
Mathematically speaking, these wormholes are built by cutting and pasting two
geometries at the throat (located at the shell), where the corresponding
junction conditions are accomplished; a geodesically complete manifold is
obtained after the matching. This formalism is also applied, for instance, to
study thin layers of matter [5] around vacuum (namely, bubbles) or black
holes.

The most straightforward generalization of the Einstein-Hilbert dynamical
action raises by taking an arbitrary function $F(R)$ of the Ricci scalar $R$,
which is usually named $F(R)$ gravity [6]. This theory, as well as many other
alternatives beyond GR, tries to deal with some of the shortcomings of GR such
that the physical meaning of the spacetime singularities and the phenomena
associated to the dark sector of the matter and the energy content of the
Universe. From the very beginning, $F(R)$ gravity has generated a broad
interest; hence, there are many publications in this gravitational framework.
Here, we just mention some papers related to the study of black holes by
considering constant [7] and non-constant [8] scalar curvature, branes [9],
and traversable wormholes [10].

It is well known that the study of thin shells follows in GR the Darmois-
Israel scheme [11]; its counterpart in $F(R)$ gravity is more restrictive [12,
13]. Besides the continuity of the first fundamental form at the matching
hypersurface, for any nonlinear function $F(R)$ the conditions also require
the continuity of the trace of the extrinsic curvature and the continuity of
the scalar curvature there, with the only exception of quadratic $F(R)$, in
which the Ricci scalar can be discontinuous [13]. In this last case, three new
contributions –an external scalar pressure/tension, an external energy flux
vector, and a double layer energy-momentum distribution– emerge in addition to
the standard energy-momentum tensor [13, 14]. The theories with the most
generic Lagrangian containing terms quadratic in the curvature [15, 16] share
the main attributes with quadratic $F(R)$. Several works can be found within
$F(R)$ gravity where the thin-shell formalism is applied to bubbles, layers of
matter enclosing black holes [17, 18], and also in the construction of
traversable wormholes [19, 20]. The peculiar case of pure double layers in the
quadratic $F(R)$ model is studied in Ref. [21]. A related paper on thin shells
within $F(R,T)$ gravity has recently appeared [22].

The interest of $(2+1)$ spacetimes dwells in some of the features which could
throw light on conceptual issues related to high energy situations [23],
present in the context of black hole singularities, quantum gravity, and
string theory. Different three-dimensional geometries were investigated over
the years in relation to black holes [24] and wormholes [25]. Thin shells of
matter [26] as well as thin-shell wormholes [27] were also widely examined.
Black hole solutions with $(2+1)$ dimensions in $F(R)$ gravity have been found
[28, 29, 30, 31]. In the recent article [32], thin shells of matter in $(2+1)$
dimensions within $F(R)$ theory have been analyzed in detail and examples of
stable configurations under radial perturbations have been found.

In this paper, we use the junction conditions in $F(R)$ gravity to construct
circular thin-shell wormholes in $(2+1)$ dimensions. The scalar curvature $R$
is assumed to be constant. We analyze the stability of the static
configurations under perturbations that preserve the symmetry. We consider
examples of anti-de Sitter spacetimes with a conformally invariant Maxwell
field as a source. In Sects. 2 and 3, we introduce the general formalism,
while in Sect. 4 we show the examples. Finally, in Sect. 5 we discuss the
results obtained. We adopt units such that $c=G=1$, with $c$ the speed of
light and $G$ the gravitational constant.

##  2 Junction formalism in (2+1)-dimensional $F(R)$ gravity

The junction formalism allows the construction of a manifold $\mathcal{M}$ as
the result of the union of two parts $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$
through a hypersurface $\Sigma$, which corresponds to a boundary hypersurface
when the matter content vanishes or to a thin shell of matter otherwise. We
begin with a brief review of this formalism within $F(R)$ theories in order to
apply it to a manifold $\mathcal{M}$ with (2+1) dimensions, so that the
hypersurface $\Sigma$ is one-dimensional. We denote the first fundamental form
(also known as induced metric) on $\Sigma$ by $h^{1,2}_{\mu\nu}$ and the
second fundamental form (or extrinsic curvature) by $K^{1,2}_{\mu\nu}$, where
the superscripts label each part of $\mathcal{M}$. The jump of any quantity
$\Upsilon$ across the hypersurface $\Sigma$ is defined by
$[\Upsilon]\equiv(\Upsilon^{2}-\Upsilon^{1})|_{\Sigma}$. The matching between
$\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ at $\Sigma$, should fulfill the so-
called junction conditions. In $F(R)$ gravity [13], the continuity of the
first fundamental form inherited from both $\mathcal{M}_{1,2}$ is required as
in General Relativity

$[h_{\mu\nu}]=0,$ (1)

but when $F^{\prime\prime\prime}(R)\neq 0$ (where the prime means the
derivative with respect to $R$) one has also to demand two additional
conditions [13]: the continuity of the trace of the second fundamental form,

$[K^{\mu}_{\;\;\mu}]=0,$ (2)

and the continuity of the scalar curvature

$\quad[R]=0.$ (3)

In this case, the dynamical equations at the joining hypersurface, with
$S_{\mu\nu}$ the energy-momentum tensor at $\Sigma$, read

$\kappa
S_{\mu\nu}=-F^{\prime}(R_{\Sigma})[K_{\mu\nu}]+F^{\prime\prime}(R_{\Sigma})[\eta^{\gamma}\nabla_{\gamma}R]h_{\mu\nu},\;\;\;\;n^{\mu}S_{\mu\nu}=0,$
(4)

where $\kappa=8\pi$ and $\nabla$ denotes the covariant derivative. When
$F^{\prime\prime\prime}(R)=0$, which corresponds to quadratic
$F(R)=R-2\Lambda+\alpha R^{2}$ gravity, the continuity of the scalar
curvature, given by Eq. (3), is no longer required. Now, the field equations
at $\Sigma$ take the form [13, 14, 15]

$\kappa
S_{\mu\nu}=-[K_{\mu\nu}]+2\alpha([n^{\gamma}\nabla_{\gamma}R]h_{\mu\nu}-[RK_{\mu\nu}]),\qquad
n^{\mu}S_{\mu\nu}=0.$ (5)

However, in the quadratic case, the energy-momentum tensor is not the only
matter contribution at $\Sigma$, it is also necessary to take into account
three other ones (for more details, see [13, 14, 15]): an external energy flux
vector $\mathcal{T}_{\mu}$

$\kappa\mathcal{T}_{\mu}=-2\alpha\overline{\nabla}_{\mu}[R],\qquad
n^{\mu}\mathcal{T}_{\mu}=0,$ (6)

with $\overline{\nabla}$ the intrinsic covariant derivative on
$(\Sigma,h_{\mu\nu})$; an external scalar pressure or tension $\mathcal{T}$

$\kappa\mathcal{T}=2\alpha[R]K^{\gamma}{}_{\gamma};$ (7)

and a two-covariant symmetric tensor distribution $\mathcal{T}_{\mu\nu}$

$\kappa\mathcal{T}_{\mu\nu}=\nabla_{\gamma}\left(2\alpha[R]h_{\mu\nu}n^{\gamma}\delta^{\Sigma}\right),$
(8)

with $\delta^{\Sigma}$ denoting the Dirac delta with support on $\Sigma$,
having a resemblance with the dipole distributions in classical
electrodynamics [13, 14, 15]. These extra contributions are required in order
to obtain a well defined energy-momentum tensor with null divergence, which is
necessary for local conservation. When $R$ has a null jump at $\Sigma$, all
these extra contributions vanish and there is an ordinary thin shell if
$S_{\mu\nu}\neq 0$. The typical scenario in quadratic $F(R)$ corresponds to
the existence of a double layer besides a thin shell at the joining
hypersurface.

##  3 Wormholes with a circular throat: construction and stability

We are interested in metrics with circular symmetry that, adopting the time
$t_{1,2}$, the radial $r>0$, and the angular $0\leq\theta\leq 2\pi$
coordinates, have the form

$ds^{2}=-A_{1,2}(r)dt_{1,2}^{2}+A_{1,2}(r)^{-1}dr^{2}+r^{2}d\theta^{2},$ (9)

where the subscripts label the regions at the two sides of the matching
hypersurface $\Sigma$, defined as a circle with radius $a$. We take the outer
region $r\geq a$ of each geometry in order to define the manifolds
$\mathcal{M}_{1,2}$. The whole manifold
$\mathcal{M}=\mathcal{M}_{1}\cup\mathcal{M}_{2}$ represents a thin-shell
wormhole, with the throat at $\Sigma$, where the flare-out condition is
satisfied. The radius of the throat is taken large enough to avoid the
presence of event horizons –if the original manifolds have them–, in order to
obtain a traversable wormhole. In $\mathcal{M}$ we adopt the coordinates
$X^{\alpha}_{1,2}=(t_{1,2},r,\theta)$, while on $\Sigma$ we use the
coordinates $\xi^{i}=(\tau,\theta)$, with $\tau$ the proper time. We let the
radius $a(\tau)$ be a function of $\tau$, and we denote its derivative with
respect to $\tau$ by an overdot. The proper time should be the same at both
sides of $\Sigma$, then
$dt_{1,2}/d\tau=\sqrt{A_{1,2}(a)+\dot{a}^{2}}/A_{1,2}(a)$, in which the free
signs are fixed by requiring that the times $t_{1,2}$ and $\tau$ all run into
the future. The first fundamental form at each side of the shell is given by

$h^{1,2}_{ij}=\left.g^{1,2}_{\mu\nu}\frac{\partial
X^{\mu}_{1,2}}{\partial\xi^{i}}\frac{\partial
X^{\nu}_{1,2}}{\partial\xi^{j}}\right|_{\Sigma},$ (10)

while the second fundamental form is determined by

$K_{ij}^{1,2}=-n_{\gamma}^{1,2}\left.\left(\frac{\partial^{2}X^{\gamma}_{1,2}}{\partial\xi^{i}\partial\xi^{j}}+\Gamma_{\alpha\beta}^{\gamma}\frac{\partial
X^{\alpha}_{1,2}}{\partial\xi^{i}}\frac{\partial
X^{\beta}_{1,2}}{\partial\xi^{j}}\right)\right|_{\Sigma},$ (11)

where the unit normals ($n^{\gamma}n_{\gamma}=1$), which are chosen to point
from $\mathcal{M}_{1}$ to $\mathcal{M}_{2}$, read111Note that the unit normal
$n^{\mu}$ to $\Sigma$ is well defined without a jump, but for computational
purposes the expressions at both sides of $\Sigma$ are usually given.

$n_{\gamma}^{1,2}=\pm\left\\{\left.\left|g^{\alpha\beta}_{1,2}\frac{\partial
G}{\partial X^{\alpha}_{1,2}}\frac{\partial G}{\partial
X^{\beta}_{1,2}}\right|^{-1/2}\frac{\partial G}{\partial
X^{\gamma}_{1,2}}\right\\}\right|_{\Sigma},$ (12)

with $G(r)\equiv r-a$ (null at $\Sigma$) and the upper and the lower signs
corresponding to $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, respectively. On the
hypersurface $\Sigma$, we prefer to work in the orthonormal basis
$\\{e_{\hat{\tau}}=e_{\tau},e_{\hat{\theta}}=a^{-1}e_{\theta}\\}$, which
allows an straightforward interpretation of the results. Then, for the general
geometries (9), the first fundamental form is
$h^{1,2}_{\hat{\imath}\hat{\jmath}}=\mathrm{diag}(-1,1)$, the unit normals are

$n_{\gamma}^{1,2}=\pm\left(-\dot{a},\frac{\sqrt{A_{1,2}(a)+\dot{a}^{2}}}{A_{1,2}(a)},0\right),$
(13)

and the only non-null components of the extrinsic curvature at each side of
$\Sigma$ result

$\begin{split}K_{\hat{\tau}\hat{\tau}}^{1,2}&=\pm\frac{A^{\prime}_{1,2}(a)+2\ddot{a}}{2\sqrt{A_{1,2}(a)+\dot{a}^{2}}},\\\
K_{\hat{\theta}\hat{\theta}}^{1,2}&=\mp\frac{1}{a}\sqrt{A_{1,2}(a)+\dot{a}^{2}},\end{split}$
(14)

from which we calculate their jumps

$\begin{split}[K_{\hat{\tau}\hat{\tau}}]&=-\frac{A^{\prime}_{2}(a)+2\ddot{a}}{2\sqrt{A_{2}(a)+\dot{a}^{2}}}-\frac{A^{\prime}_{1}(a)+2\ddot{a}}{2\sqrt{A_{1}(a)+\dot{a}^{2}}},\\\
[K_{\hat{\theta}\hat{\theta}}]&=\frac{1}{a}\sqrt{A_{2}(a)+\dot{a}^{2}}+\frac{1}{a}\sqrt{A_{1}(a)+\dot{a}^{2}}.\end{split}$
(15)

The continuity of the trace of the extrinsic curvature at the shell, i.e. Eq.
(3), gives the following equation

$\frac{2\ddot{a}+A_{2}^{\prime}(a)}{2\sqrt{A_{2}(a)+\dot{a}^{2}}}+\frac{2\ddot{a}+A_{1}^{\prime}(a)}{2\sqrt{A_{1}(a)+\dot{a}^{2}}}+\frac{1}{a}\left(\sqrt{A_{1}(a)+\dot{a}^{2}}+\sqrt{A_{2}(a)+\dot{a}^{2}}\right)=0,$
(16)

which, for the static configurations with a shell radius $a_{0}$, takes the
form

$\frac{A_{2}^{\prime}(a_{0})}{2\sqrt{A_{2}(a_{0})}}+\frac{A_{1}^{\prime}(a_{0})}{2\sqrt{A_{1}(a_{0})}}+\frac{1}{a_{0}}\left(\sqrt{A_{1}(a_{0})}+\sqrt{A_{2}(a_{0})}\right)=0.$
(17)

We now have all the ingredients to find the matter content at the matching
hypersurface and then to perform the stability analysis of the thin shell. The
energy-momentum tensor in the orthonormal basis has the form
$S_{{}_{\hat{\imath}\hat{\jmath}}}={\rm diag}(\sigma,p)$, where $\sigma$ is
the energy density and $p=p_{\hat{\theta}}$ is the transverse pressure. In
what follows, the scalar curvature is constant in each of the regions at the
sides of the throat, then $[\eta^{\gamma}\nabla_{\gamma}R]=0$, which
simplifies Eqs. (4) and (5), while the so-called [13] brane tension
$\lambda=F^{\prime\prime}(R_{\Sigma})[\eta^{\gamma}\nabla_{\gamma}R]$ has a
null value at $\Sigma$.

###  3.1 General $F(R)$ gravity with $[R]=0$

For an arbitrary non-quadratic $F(R)$ theory of gravity, there are three
junction conditions which must be satisfied, that is Eqs. (1), (2), and (3).
Then, for a constant value $R_{0}$ at both sides of the throat, from Eq. (4)
we obtain the energy density and the pressure on the shell

$\displaystyle\sigma$
$\displaystyle=\frac{F^{\prime}(R_{0})}{\kappa}\left(\frac{2\ddot{a}+A_{2}^{\prime}(a)}{2\sqrt{A_{2}(a)+\dot{a}^{2}}}+\frac{2\ddot{a}+A_{1}^{\prime}(a)}{2\sqrt{A_{1}(a)+\dot{a}^{2}}}\right),$
(18) $\displaystyle p$
$\displaystyle=-\frac{F^{\prime}(R_{0})}{a\kappa}\left(\sqrt{A_{2}(a)+\dot{a}^{2}}+\sqrt{A_{1}(a)+\dot{a}^{2}}\right).$
(19)

By using the expression (16), we rewrite the energy density as

$\sigma=-\frac{F^{\prime}(R_{0})}{a\kappa}\left(\sqrt{A_{2}(a)+\dot{a}^{2}}+\sqrt{A_{1}(a)+\dot{a}^{2}}\right),$
(20)

so that we can relate it with the pressure by $\sigma-p=0$. Since the
stability analysis is done with respect to the static configurations, we show
the corresponding expressions of the energy density and the pressure

$\displaystyle\sigma_{0}$
$\displaystyle=\frac{F^{\prime}(R_{0})}{\kappa}\left(\frac{A_{2}^{\prime}(a_{0})}{2\sqrt{A_{2}(a_{0})}}+\frac{A_{1}^{\prime}(a_{0})}{2\sqrt{A_{1}(a_{0})}}\right),$
(21) $\displaystyle p_{0}$
$\displaystyle=-\frac{F^{\prime}(R_{0})}{a_{0}\kappa}\left(\sqrt{A_{2}(a_{0})}+\sqrt{A_{1}(a_{0})}\right).$
(22)

Now we have that $\sigma_{0}-p_{0}=0$. In quadratic $F(R)$, when $[R]=0$ all
these equations are also valid with $F^{\prime}(R_{0})=1+2\alpha R_{0}$, as it
can be easily seen from Eq. (5); the extra contributions proportional to $[R]$
vanish in this case.

###  3.2 Quadratic $F(R)$ gravity with $[R]\neq 0$

As it is was stated in Sect. 2, the particular case of quadratic $F(R)$
gravity is less restricted: the condition about the continuity of the scalar
curvature (3) is no longer necessary, so we can take constant values
$R_{1}\neq R_{2}$ at the sides of the throat. Then, we get from the dynamical
equations (5) for this case

$\displaystyle\sigma$ $\displaystyle=\frac{1+2\alpha
R_{2}}{\kappa}\left(\frac{2\ddot{a}+A_{2}^{\prime}(a)}{2\sqrt{A_{2}(a)+\dot{a}^{2}}}\right)+\frac{1+2\alpha
R_{1}}{\kappa}\left(\frac{2\ddot{a}+A_{1}^{\prime}(a)}{2\sqrt{A_{1}(a)+\dot{a}^{2}}}\right),$
(23) $\displaystyle p$ $\displaystyle=-\frac{1+2\alpha
R_{2}}{\kappa}\left(\frac{\sqrt{A_{2}(a)+\dot{a}^{2}}}{a}\right)-\frac{1+2\alpha
R_{1}}{\kappa}\left(\frac{\sqrt{A_{1}(a)+\dot{a}^{2}}}{a}\right).$ (24)

The energy-momentum tensor should be completed by the extra contributions.
Namely, the external scalar pressure/tension

$\mathcal{T}=\frac{2\alpha
R_{2}}{\kappa}\left(\frac{2\ddot{a}+A_{2}^{\prime}(a)}{2\sqrt{A_{2}(a)+\dot{a}^{2}}}+\frac{\sqrt{A_{2}(a)+\dot{a}^{2}}}{a}\right)+\frac{2\alpha
R_{1}}{\kappa}\left(\frac{2\ddot{a}+A_{1}^{\prime}(a)}{2\sqrt{A_{1}(a)+\dot{a}^{2}}}+\frac{\sqrt{A_{1}(a)+\dot{a}^{2}}}{a}\right),$
(25)

from which we obtain the relation $\sigma-p=\mathcal{T}$; while the external
energy flux vector results $\mathcal{T}_{\mu}=0$ and the two-covariant
symmetric tensor distribution $\mathcal{T}_{\mu\nu}$ is proportional to
$(2\alpha[R]/\kappa)[R]h_{\mu\nu}$. The static values of the energy density
and the pressure read

$\displaystyle\sigma_{0}$ $\displaystyle=\frac{1+2\alpha
R_{2}}{\kappa}\left(\frac{A_{2}^{\prime}(a_{0})}{2\sqrt{A_{2}(a_{0})}}\right)+\frac{1+2\alpha
R_{1}}{\kappa}\left(\frac{A_{1}^{\prime}(a_{0})}{2\sqrt{A_{1}(a_{0})}}\right),$
(26) $\displaystyle p_{0}$ $\displaystyle=-\frac{1+2\alpha
R_{2}}{\kappa}\left(\frac{\sqrt{A_{2}(a_{0})}}{a_{0}}\right)-\frac{1+2\alpha
R_{1}}{\kappa}\left(\frac{\sqrt{A_{1}(a_{0})}}{a_{0}}\right);$ (27)

the extra contributions adopt the form

$\mathcal{T}_{0}=\frac{2\alpha
R_{2}}{\kappa}\left(\frac{A_{2}^{\prime}(a_{0})}{2\sqrt{A_{2}(a_{0})}}+\frac{\sqrt{A_{2}(a_{0})}}{a_{0}}\right)+\frac{2\alpha
R_{1}}{\kappa}\left(\frac{A_{1}^{\prime}(a_{0})}{2\sqrt{A_{1}(a_{0})}}+\frac{\sqrt{A_{1}(a_{0})}}{a_{0}}\right),$
(28)

which satisfies $\sigma_{0}-p_{0}=\mathcal{T}_{0}$; while
$\mathcal{T}_{\mu}^{(0)}=0$ and $\mathcal{T}_{\mu\nu}^{(0)}$ is proportional
to $(2\alpha[R]/\kappa)h_{\mu\nu}$.

###  3.3 Stability

In order to analyze the stability of the static circular throats, we employ
the standard potential analogy method. By using
$\ddot{a}=(1/2)d(\dot{a}^{2})/da$, we can rewrite Eq. (16) to obtain the
equivalent equation $az^{\prime}(a)+z(a)=0$ where
$z=\sqrt{A_{2}(a)+\dot{a}^{2}}+\sqrt{A_{1}(a)+\dot{a}^{2}}$, which represent a
differential equation for $\dot{a}^{2}$ in terms of an effective potential
$\dot{a}^{2}=-V(a)$ given by

$V(a)=-\frac{a^{2}\left(A_{1}(a)-A_{2}(a)\right)^{2}}{4a_{0}^{2}\left(\sqrt{A_{1}\left(a_{0}\right)}+\sqrt{A_{2}\left(a_{0}\right)}\right)^{2}}-\frac{a_{0}^{2}\left(\sqrt{A_{1}\left(a_{0}\right)}+\sqrt{A_{2}\left(a_{0}\right)}\right)^{2}}{4a^{2}}+\frac{A_{1}(a)+A_{2}(a)}{2}.$
(29)

Since $V(a_{0})=0$ and $V^{\prime}(a_{0})=0$, the stability of the
configurations under radial perturbations is determined by the positive sign
of $V^{\prime\prime}(a_{0})$, which reads

$\displaystyle V^{\prime\prime}(a_{0})$ $\displaystyle=$
$\displaystyle-\frac{3\left(\sqrt{A_{1}\left(a_{0}\right)}+\sqrt{A_{2}\left(a_{0}\right)}\right)^{2}}{2a_{0}^{2}}-\frac{\left(\sqrt{A_{1}\left(a_{0}\right)}-\sqrt{A_{2}\left(a_{0}\right)}\right)^{2}}{2a_{0}^{2}}$
(30)
$\displaystyle-\frac{\left(A_{1}^{\prime}(a_{0})-A_{2}^{\prime}(a_{0})\right)^{2}}{2\left(\sqrt{A_{1}\left(a_{0}\right)}+\sqrt{A_{2}\left(a_{0}\right)}\right)^{2}}-\frac{2\left(A_{1}(a_{0})-A_{2}(a_{0})\right)\left(A_{1}^{\prime}(a_{0})-A_{2}^{\prime}(a_{0})\right)}{a_{0}\left(\sqrt{A_{1}\left(a_{0}\right)}+\sqrt{A_{2}\left(a_{0}\right)}\right)^{2}}$
$\displaystyle-\frac{\left(A_{1}(a_{0})-A_{2}(a_{0})\right)\left(A_{1}^{\prime\prime}(a_{0})-A_{2}^{\prime\prime}(a_{0})\right)}{2\left(\sqrt{A_{1}\left(a_{0}\right)}+\sqrt{A_{2}\left(a_{0}\right)}\right)^{2}}+\frac{A_{1}^{\prime\prime}(a_{0})+A_{2}^{\prime\prime}(a_{0})}{2}.$

Now we are in position to explore some concrete examples.

##  4 Anti-de Sitter thin-shell wormholes with charge

In $F(R)=R+f(R)$ gravity222Here use the nomenclature in which the first term
corresponds to General Relativity. coupled to nonlinear electrodynamics, the
$(2+1)$-dimensional action reads

$I=\frac{1}{16\pi}\int
d^{3}x\sqrt{-g}\left(R+f(R)+\left(-F_{\alpha\beta}F^{\alpha\beta}\right)^{s}\right),$
(31)

with $s$ an arbitrary positive nonlinearity parameter ($s\neq 1/2$) and
$F_{\alpha\beta}F^{\alpha\beta}$ the Maxwell invariant, where
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the
electromagnetic tensor field, being $A_{\mu}$ the gauge potential. The
corresponding energy-momentum tensor has the form
$T_{\mu\nu}=1/(4\pi)\left[-sF_{\mu\gamma}F_{\nu}^{\;\;\gamma}(-F_{\alpha\beta}F^{\alpha\beta})^{s-1}-(1/4)g_{\mu\nu}(-F_{\alpha\beta}F^{\alpha\beta})^{s}\right]$.
In the case with $s=3/4$, the conformally invariant Maxwell field as a source
is obtained [28], having a traceless $T_{\mu\nu}$. The field equations in the
metric formalism are

$R_{\mu\nu}(1+f^{\prime}(R))-\frac{1}{2}g_{\mu\nu}(R+f(R))+(g_{\mu\nu}\nabla_{\gamma}\nabla^{\gamma}-\nabla_{\mu}\nabla_{\nu})f^{\prime}(R)=8\pi
T_{\mu\nu},$ (32)
$\partial_{\mu}\left(\sqrt{-g}F^{\mu\nu}(-F_{\alpha\beta}F^{\alpha\beta})^{-1/4}\right)=0.$
(33)

By considering a constant scalar curvature $R_{0}$, these field equations have
a solution of the form (9), with the metric function [28]

$A(r)=-M-\frac{\left(2\mathcal{Q}^{2}\right)^{3/4}}{2\left(1+f^{\prime}(R_{0})\right)r}-\frac{r^{2}R_{0}}{6},$
(34)

being $M$ the mass and $\mathcal{Q}$ the charge. The only non-null independent
component of the electromagnetic field is $F_{tr}=\mathcal{Q}/r^{2}$. When
$\mathcal{Q}=0$, the well-known vacuum static BTZ geometry [24] of General
Relativity is obtained, which is also a solution in $F(R)$ gravity [29] with
constant scalar curvature. The spacetime is asymptotically anti-de Sitter
since it requires that $R_{0}<0$. The geometry has a curvature singularity at
$r=0$ because the Kretschmann scalar diverges there [28]. The trace of the
field equations

$R_{0}\left(1+f^{\prime}(R_{0})\right)-\frac{3}{2}\left(R_{0}+f(R_{0})\right)=0,$
(35)

allows to define the effective cosmological constant $\Lambda_{e}$

$R_{0}=\frac{3f(R_{0})}{2f^{\prime}(R_{0})-1}\equiv 6\Lambda_{e}.$ (36)

We can also define an effective charge by

$Z=\frac{\left(2\mathcal{Q}^{2}\right)^{3/4}}{2\left(1+f^{\prime}(R_{0})\right)},$
(37)

which can be positive or negative, depending on the sign of
$F^{\prime}(R_{0})=1+f^{\prime}(R_{0})$. The effective Newton constant
$G_{\mathrm{eff}}=G/F^{\prime}(R)=1/F^{\prime}(R)$ is positive when
$F^{\prime}(R)>0$, preventing in this case the graviton to be a ghost [6] (see
also Ref. [33]). Note that the sign of $Z$ is fixed by the choice of the
$F(R)$ theory and the value $R_{0}$, while the squared charge
$\mathcal{Q}^{2}$ can only modify the absolute value of $Z$. As usual, the
radii of the horizons are determined by solving $A(r)=0$ and taking the real
and positive solutions333For $Z\neq 0$ the analytic expressions for the
solutions of the cubic equations are cumbersome, hence it is meaningless to
show them here.. Then, there are three possible cases:

  * $-$

If $Z\geq 0$, there is only a solution with $r_{h}$ the radius of the event
horizon; in particular, for $Z=0$ (which corresponds to a quadratic equation)
one gets $r_{h}=\sqrt{-6M/R_{0}}$.

  * $-$

If $Z_{c}<Z<0$ ($Z_{c}=(-2M/3)\sqrt{-2M/R_{0}}$) there are two solutions, one
corresponds to the radius of the event horizon $r_{h}$, and the other to the
radius of the inner horizon $r_{i}<r_{h}$; when $Z=Z_{c}$ both horizons merge
into one.

  * $-$

If $Z<Z_{c}$, there are no horizons and the singularity at the origin is
naked.

In the following, we study some examples by taking into account the metric
function (34). We construct static thin-shell wormholes with radius $a_{0}$ by
applying the junction formalism, and we analyze their stability under radial
perturbations which preserve the symmetry. We also determine the kind of
matter (ordinary or exotic) located at the thin-shell by checking the weak
energy condition (WEC)444In the orthonormal basis, it takes the form
$\sigma_{0}\geq 0$ and $\sigma_{0}+p_{0}\geq 0$..

###  4.1 Wormhole symmetric across the throat

Figure 1: Wormhole with the same value of the scalar curvature $R_{0}$, mass
$M$, and charge $\mathcal{Q}$ at the sides of the throat with radius $a_{0}$,
in a general $F(R)$ theory. The parameter
$Z=\left(2\mathcal{Q}^{2}\right)^{3/4}/\left(2\left(1+f^{\prime}(R_{0})\right)\right)$
has the same sign as $F^{\prime}(R_{0})=1+f^{\prime}(R_{0})$ (see text). The
meshed zones represent normal matter (exotic otherwise) and the gray ones have
no physical meaning (see text). Left: $M=0.5$; right: $M=1$. The solid lines
correspond to the stable configurations while the dotted lines to the unstable
ones.

We proceed with the construction of a thin-shell wormhole symmetric across the
throat, that is, with equal values of the scalar curvature
$R_{1}=R_{2}=R_{0}$, the mass $M_{1}=M_{2}=M$, and the charge
$\mathcal{Q}_{1}=\mathcal{Q}_{2}=\mathcal{Q}$, so that the metric functions,
given by Eq. (34), are the same $A_{1}(r)=A_{2}(r)=A(r)$ in both regions. When
$Z\geq Z_{c}=(-2M/3)\sqrt{-2M/R_{0}}$, the value of the throat radius $a_{0}$
is taken larger than $r_{h}$ in order to remove the region inside the horizon
corresponding to the original geometry, while for $Z<Z_{c}$ any value
$a_{0}>0$ removes the naked singularity. The proper matching at $\Sigma$
requires that $a_{0}$ can only take values satisfying Eq. (17). The energy
density $\sigma_{0}$ and the pressure $p_{0}$ are obtained by replacing the
metric functions in Eqs. (21) and (22), respectively. A configuration is
stable under radial perturbations when the second derivative of the potential
evaluated at $a_{0}$, obtained by replacing the metric functions and their
derivatives in Eq. (30), is positive.

Some representative results are shown in Fig. 1. In all plots, the meshed
zones represent normal matter satisfying the weak energy condition (exotic
matter otherwise) and the gray areas have no physical meaning, corresponding
to the removed part of the original manifold. The solid lines correspond to
stable static solutions, while the dotted lines to the unstable ones. In the
left plot we have adopted $M=0.5$ and in the right one $M=1$. We see that
there exist two solutions made of normal matter for a short range of
$\sqrt{|R_{0}|}Z$, corresponding to negative values of $Z$ thus requiring the
presence of ghost fields. The solution with the largest radius is stable while
the other one is unstable under radial perturbations. A modification in the
value of the mass $M$ only results in a change of scale, without affecting the
qualitative behavior of the solutions.

###  4.2 Wormhole asymmetric in the mass and the charge

Figure 2: Wormhole with the same value $R_{0}$ of the scalar curvature, but
different values of mass $M_{1}\neq M_{2}$ at the sides of the throat with
radius $a_{0}$, in any $F(R)$ theory. The parameters
$Z_{1,2}=\left(2\mathcal{Q}_{1,2}^{2}\right)^{3/4}/F^{\prime}(R_{0})$, with
$\mathcal{Q}_{1,2}$ the values of charge, have the same sign as
$F^{\prime}(R_{0})$ (see text). The solid and dotted lines, and the style of
the different zones have the same meaning as in Fig. 1. In all plots
$M_{1}=0.5$ and $M_{2}=1$. Upper row, left: $Z_{1}=0.4Z_{2}$; upper row,
right: $Z_{1}=0.7Z_{2}$; lower row, left: $Z_{1}=Z_{2}$; lower row, right:
$Z_{1}=1.3Z_{2}$.

Now we construct a wormhole with the same value $R_{1}=R_{2}=R_{0}$ of the
scalar curvature, but different masses $M_{1}$ and $M_{2}$, and charges
$\mathcal{Q}_{1}$ and $\mathcal{Q}_{2}$ that can be equal or different at the
sides of the shell. The metric functions $A_{1}(r)$ and $A_{2}(r)$ are given
by Eq. (34), in which we use the definitions
$Z_{1,2}=\left(2\mathcal{Q}_{1,2}^{2}\right)^{3/4}/(2\left(1+f^{\prime}(R_{0})\right))$.
In order to prevent the presence of horizons and singularities, we take the
radius of the throat $a_{0}$ larger than the horizon radii $r_{h}^{(1)}$ and
$r_{h}^{(2)}$ –when present– of the original manifolds. The solutions with
radius $a_{0}$ have to satisfy the Eq. (17), and the energy density
$\sigma_{0}$ and the pressure $p_{0}$ result by replacing the metric functions
in Eqs. (21) and (22), respectively. A configuration is stable under radial
perturbations when the second derivative of the potential
$V^{\prime\prime}(a_{0})$, shown in Eq. (30), is positive.

In Fig. 2 we show the most illustrative results. Again, in all plots the
meshed zones represent normal matter (exotic otherwise) and the gray areas
have no physical meaning, corresponding to the removed parts of the original
manifolds. The solid lines correspond to stable static solutions, while the
dotted lines to the unstable ones. In all plots the values of the masses are
$M_{1}=0.5$ and $M_{2}=1$, while the relationship between the charges is
$Z_{1}=\eta Z_{2}$, with $\eta$ taking the values $0.4$, $0.7$, $1$, and
$1.3$, respectively. We can see that in all cases there exist two solutions
made of normal matter for a short range of negative values of the charge
$Z_{2}$. This means that solutions are only possible under the presence of
ghosts. Regarding their stability, in all cases the solution with the smallest
radius is unstable, while the one with the largest one is stable. When the
value of $\eta$ increases, the range of $Z_{2}$ for which the solutions exist
becomes smaller. A modification in the values of the non-null masses $M_{1}$
and $M_{2}$ causes a change of scale without affecting the qualitative
behavior of the solutions.

###  4.3 Wormhole asymmetric in the scalar curvature

Figure 3: Wormhole with different values $R_{1}\neq R_{2}$ of the scalar
curvature, and the same values of mass $M$ and charge $\mathcal{Q}$, at the
sides of the throat with radius $a_{0}$, only in quadratic $F(R)$ theory. The
parameter $Z=\left(2\mathcal{Q}^{2}\right)^{3/4}/\left(2\left(1+2\alpha
R_{A}\right)\right)$, with $R_{A}=(R_{1}+R_{2})/2$, has the same sign as
$F^{\prime}(R_{A})=1+2\alpha R_{A}$ (see text). The solid and dotted lines,
and the style of the different zones have the same meaning as in Fig. 1. In
all plots $M=1$. Upper row: $R_{1}=0.9R_{A}$ and $R_{2}=1.1R_{A}$; lower row:
$R_{1}=0.8R_{A}$ and $R_{2}=1.2R_{A}$. Left column: $\alpha R_{A}=-0.1$ for
$Z>0$ and $\alpha R_{A}=-0.9$ for $Z<0$; right column: $\alpha R_{A}=-0.2$ for
$Z>0$ and $\alpha R_{A}=-0.8$ for $Z<0$

.

As a final example, we analyze a thin-shell wormhole with equal values of the
mass $M_{1}=M_{2}=M$ and the charge
$\mathcal{Q}_{1}=\mathcal{Q}_{2}=\mathcal{Q}$ for both regions, but with
different values of the scalar curvature $R_{1}\neq R_{2}$, so that we have to
work in quadratic $F(R)$. The metric functions $A_{1}(r)$ and $A_{2}(r)$ are
then given by Eq. (34), with $R_{1}$ or $R_{2}$ as appropriate. For
convenience, we define the average value $R_{A}=(R_{1}+R_{2})/2$, which will
be useful in the presentation of the results. We introduce the parameter
$Z=\left(2\mathcal{Q}^{2}\right)^{3/4}/\left(2\left(1+2\alpha
R_{A}\right)\right)$, which has the same sign as $F^{\prime}(R_{A})=1+2\alpha
R_{A}$. With this definition, we have that $Z_{1,2}=Z\left(1+2\alpha
R_{A}\right)/\left(1+2\alpha R_{1,2}\right)$. As it was previously explained
in detail, the throat radius $a_{0}$ should be larger than the radius of any
of the horizons $r_{h}^{(1)}$ and $r_{h}^{(2)}$ of the original manifolds in
order to remove them and the regions inside, and also satisfy Eq. (17). The
energy density $\sigma_{0}$, the pressure $p_{0}$, and the external scalar
pressure or tension $\mathcal{T}_{0}$ are calculated by replacing the metric
functions in Eqs. (26), (27), and (28), respectively; while –as shown above–
there is a null external energy flux vector $\mathcal{T}_{\mu}$ and the double
layer energy-momentum distribution $\mathcal{T}_{\mu\nu}^{(0)}$ is
proportional to $(2\alpha[R]/\kappa)h_{\mu\nu}$. A configuration is stable
under radial perturbations when the second derivative of the potential, shown
in Eq. (30), fulfill the inequality $V^{\prime\prime}(a_{0})>0$.

Some representative results are displayed in Fig. 3. The styles of the
different regions and the lines have the same meaning explained above. We have
taken $M=1$ in all plots, while the values of the scalar curvature are
$R_{1}=0.9R_{A}$ and $R_{2}=1.1R_{A}$ in the upper row, and $R_{1}=0.8R_{A}$
$R_{2}=1.2R_{A}$ in the lower one. In the left column we have adopted $\alpha
R_{A}=-0.1$ for $Z>0$ and $\alpha R_{A}=-0.9$ for $Z<0$; while in the right
one the values are $\alpha R_{A}=-0.2$ for $Z>0$ and $\alpha R_{A}=-0.8$ for
$Z<0$. In the plots, we find two solutions for a short range of negative
values of $\sqrt{|R_{A}|}Z$, the one with the largest radius is stable under
radial perturbations while the other is unstable; both of them are made of
normal matter with the presence of ghost fields. We can also see that the
larger the absolute value of the difference between the scalar curvature
$R_{1}$ and $R_{2}$ across the throat, the smaller the range of values of
$\sqrt{|R_{A}|}Z$ where these solutions exist. A similar behavior can be seen
when increasing the difference between the values of the theory parameter
$\alpha$ at the sides of the shell, which reduces the range of values of
$\sqrt{|R_{A}|}Z$ where the solutions can be found. A variation in the mass
$M$, not illustrated in the figure for brevity, leads again only to a scale
change.

##  5 Discussion

In this work, we have presented a wide class of (2+1)-dimensional thin-shell
wormholes with circular symmetry, within the framework of $F(R)$ theories of
gravity with constant scalar curvature $R$. We have analyzed the matter
content at the throat where the thin shell is located and we have studied the
stability of the static configurations under perturbations preserving the
symmetry.

We have considered three examples of wormholes in $F(R)$ gravity coupled to
conformally invariant nonlinear electrodynamics, which are asymptotically
anti-de Sitter at both sides of the throat. In the first example, the
spacetime is symmetric across the throat, in the second one is asymmetric in
the mass and the charge, and in the third one –only in quadratic $F(R)$– is
asymmetric in the scalar curvature. We have obtained that, in all cases,
solutions are present only for negative values of $F^{\prime}(R)$, which means
that the presence of ghost fields is always required. In the first example,
for symmetric wormholes with fixed (negative) scalar curvature $R_{0}$, we
have found two solutions made of normal matter –satisfying WEC– for a short
range of the squared charge $\mathcal{Q}^{2}$. The solution with the largest
radius is stable while the other is unstable under radial perturbations. Any
change in the mass $M$ is reflected as a change of scale conserving the
general behavior of the solutions. In the second example, corresponding to the
wormhole with the same scalar curvature $R_{0}$ at both sides of the shell but
asymmetric in the mass and the charge, we have also found two solutions with
normal matter having the same stability characteristics displayed in the
symmetric case. Adopting $M_{1}<M_{2}$, if we increase the ratio
$\mathcal{Q}_{1}^{2}/\mathcal{Q}_{2}^{2}$, the solutions exist for a smaller
range of values of charge. Instead, when we modify the masses (maintaining
that $M_{1}<M_{2}$) we only see a change of scale without altering the
qualitative behavior of the solutions. Finally, in the third example where the
wormhole is only asymmetric in the scalar curvature, there are two solutions
made of normal matter, with the same stability behavior as in previous
examples. When the absolute value of the difference between the constant
(negative) scalar curvatures $R_{1}$ and $R_{2}$ of the two regions grows, the
range of $\mathcal{Q}^{2}$ where these solutions exist becomes smaller.
Analogous behavior can be observed when we decrease the value of the quadratic
parameter $\alpha$ of the theory; as it happens with the symmetric case, any
change in $M$ results in a change of scale.

Summarizing, we can say that in our examples of symmetric and asymmetric
wormholes across the throat –whether we consider general theories or only
quadratic ones–, a pair of solutions, made of normal matter but requiring the
presence of ghost fields, exists for a certain set of the parameters of the
theory. We have found that under radial perturbations one of these solutions
is stable while the other is unstable.

## Acknowledgments

This work has been supported by CONICET and Universidad de Buenos Aires. C. B.
thanks the partial support of the John Templeton Foundation.

## References

  * [1] M. S. Morris and K. S. Thorne, Am. J. Phys. 56, 395 (1988); M. Visser, Lorentzian Wormholes (AIP Press, NewYork, 1996).
  * [2] C. Bejarano, G. J. Olmo, and D. Rubiera-Garcia, Phys. Rev. D 95, 064043 (2017); C. Bejarano, F. S. N. Lobo, G. J. Olmo, and D. Rubiera-Garcia, Eur. Phys. J. C 77, 776 (2017).
  * [3] E. Poisson and M. Visser, Phys. Rev. D 52, 7318 (1995); E.F. Eiroa and G.E. Romero, Gen. Relativ. Gravit. 36, 651 (2004); E.F. Eiroa, Phys. Rev. D 78, 024018 (2008); N. Montelongo Garcia, F.S.N. Lobo, and M. Visser, Phys. Rev. D 86, 044026 (2012); S.D. Forghani, S. Habib Mazharimousavi, and M. Halilsoy, Eur. Phys. J. C 78, 469 (2018); T. Berry, F.S.N. Lobo, A. Simpson, and M. Visser, Phys. Rev. D 102, 064054 (2020).
  * [4] E.F. Eiroa and C. Simeone, Phys. Rev. D 81, 084022 (2010); 90, 089906(E) (2014); S. Habib Mazharimousavi, M. Halilsoy, and Z. Amirabi, Phys. Rev. D 89, 084003 (2014); E.F. Eiroa, E. Rubín de Celis, and C. Simeone, Eur. Phys. J. C 79, 272 (2019); S.D. Forghani, S. Habib Mazharimousavi, and M. Halilsoy, J. Cosmol. Astropart. Phys. 10, 067 (2019).
  * [5] P.R. Brady, J. Louko, and E. Poisson, Phys. Rev. D 44, 1891 (1991); M. Ishak and K. Lake, Phys. Rev. D 65, 044011 (2002); F.S.N. Lobo and P. Crawford, Class. Quantum Gravity 22, 4869 (2005); E.F. Eiroa and C. Simeone, Phys. Rev. D 83, 104009 (2011).
  * [6] T.P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82, 451 (2010); A. De Felice and S. Tsujikawa, Living Rev. Relativ. 13, 3 (2010); S. Nojiri, S.D. Odintsov, and V.K. Oikonomou, Phys. Rep. 692, 1 (2017).
  * [7] T. Multamäki and I. Vilja, Phys. Rev. D 74, 064022 (2006); S. Capozziello, A. Stabile, and A. Troisi, Class. Quantum Gravity 25, 085004 (2008); A. de la Cruz-Dombriz, A. Dobado, and A.L. Maroto, Phys. Rev. D 80, 124011 (2009); 83, 029903(E) (2011); T. Moon, Y.S. Myung, and E.J. Son, Gen. Relativ. Gravit. 43, 3079 (2011).
  * [8] L. Sebastiani and S. Zerbini, Eur. Phys. J. C 71, 1591 (2011); S. Habib Mazharimousavi, M. Halilsoy, and T. Tahamtan, Eur. Phys. J. C 72, 1851 (2012); E. Elizalde, G.G.L. Nashed, S. Nojiri, and S.D. Odintsov, Eur. Phys. J. C 80, 109 (2020).
  * [9] S. Chakraborty and S. SenGupta, Class. Quantum Gravity 33, 225001 (2016); N. Dimakis, A. Giacomini, and A. Paliathanasis, Phys. Lett. B 777, 361 (2018).
  * [10] A. DeBenedictis and D. Horvat, Gen. Relativ. Gravit. 44, 2711 (2012); T. Harko, F.S.N. Lobo, M.K. Mak, and S.V. Sushkov, Phys. Rev. D 87, 067504 (2013); J.L. Rosa, J.P.S. Lemos, and F.S.N. Lobo, Phys. Rev. D 98, 064054 (2018); H. Golchin and M.R. Mehdizadeh, Eur. Phys. J. C 79, 777 (2019); F.S.N. Lobo, G.J. Olmo, E. Orazi, D. Rubiera-Garcia, and A. Rustam, Phys. Rev. D 102, 104012 (2020).
  * [11] G. Darmois, Mémorial des Sciences Mathématiques, Fascicule XXV, Chap. V (Gauthier-Villars, Paris, 1927); W. Israel, Nuovo Cimento B 44, 1 (1966); 48, 463(E) (1967).
  * [12] N. Deruelle, M. Sasaki, and Y. Sendouda, Prog. Theor. Phys. 119, 237 (2008).
  * [13] J.M.M. Senovilla, Phys. Rev. D 88, 064015 (2013).
  * [14] J.M.M. Senovilla, Class. Quantum Gravity 31, 072002 (2014).
  * [15] B. Reina, J.M.M. Senovilla, and R. Vera, Class. Quantum Gravity 33, 105008 (2016).
  * [16] V.A. Berezin, V.I. Dokuchaev, Y.N. Eroshenko, and A.L. Smirnov, Class. Quantum Gravity, 38, 045014 (2021).
  * [17] E.F. Eiroa and G. Figueroa-Aguirre, Eur. Phys. J. C 78, 54 (2018); E.F. Eiroa and G. Figueroa-Aguirre, Eur. Phys. J. C 79, 171 (2019).
  * [18] E.F. Eiroa and G. Figueroa-Aguirre, Eur. Phys. J. Plus 135, 774 (2020).
  * [19] E.F. Eiroa and G. Figueroa-Aguirre, Eur. Phys. J. C 76, 132 (2016); E.F. Eiroa and G. Figueroa-Aguirre, Phys. Rev. D 94, 044016 (2016).
  * [20] M. Zaeem-ul-Haq Bhatti, A. Anwar, and S. Ashraf, Mod. Phys. Lett. A 32, 1750111 (2017); S. Habib Mazharimousavi, Eur. Phys. J. C 78, 612 (2018); S. Habib Mazharimousavi, M. Halilsoy, and K. Kianfar, Eur. Phys. J. Plus 135, 440 (2020).
  * [21] E.F. Eiroa, G. Figueroa-Aguirre, and J.M.M. Senovilla, Phys. Rev. D 95, 124021 (2017).
  * [22] J.L. Rosa, Phys. Rev. D 103, 104069 (2021).
  * [23] S. Carlip, Living Rev. Relativ. 8, 1 (2005).
  * [24] M. Bañados, C. Teitelboim, and J. Zanelli, Phys. Rev. Lett. 69, 1849 (1992); M. Bañados, M. Henneaux, C. Teitelboim, and J. Zanelli, Phys. Rev. D 48, 1506 (1993).
  * [25] G.P. Perry and R.B. Mann, Gen. Relativ. Gravit. 24, 305 (1992); S.W. Kim, H.J. Lee, S.K. Kim, and J. Yang, Phys. Lett. A 183, 359 (1993); M.S.R. Delgaty and R.B. Mann, Int. J. Mod. Phys. D 4, 231 (1995); W.T. Kim, J.J. Oh, and M.S. Yoon, Phys. Rev. D 70, 044006 (2004); M. Jamil, M.U. Farooq, Int. J. Theor. Phys. 49, 835 (2010).
  * [26] E.F. Eiroa and C. Simeone, Phys. Rev. D 87, 064041 (2013).
  * [27] F. Rahaman, A. Banerjee, and I. Radinschi, Int. J. Theor. Phys. 51, 1680 (2012); A. Banerjee, Int. J. Theor. Phys. 52, 2943 (2013); C. Bejarano, E.F. Eiroa, and C. Simeone, Eur. Phys. J. C 74, 3015 (2014); S. Habib Mazharimousavi, M. Halilsoy, Eur. Phys. J. C, 75, 81 (2015); P. Bhar and A. Banerjee, Int. J. Mod. Phys. D, 24, 1550034 (2015).
  * [28] S.H. Hendi, B. Eslam Panah, and R. Saffari, Int. J. Mod. Phys. D 23, 1450088 (2014).
  * [29] S.H. Hendi, Int. J. Theor. Phys. 53, 4170 (2014).
  * [30] S.H. Hendi, R. Ramezani-Arani, and E. Rahimi, Phys. Lett. B 805 135436 (2020).
  * [31] T. Karakasis, E. Papantonopoulos, Z.-Y. Tang, and B. Wang, Phys. Rev. D 103, 064063 (2021).
  * [32] E. F. Eiroa and G. Figueroa-Aguirre, Phys. Rev. D 103, 044011 (2021).
  * [33] K.A. Bronnikov, M.V. Skvortsova, and A.A. Starobinsky, Grav. Cosmol. 16, 216 (2010).