#  Nondivergent deflection of light around a photon sphere of a compact object

Ryuya Kudo  Hideki Asada  Graduate School of Science and Technology, Hirosaki
University, Aomori 036-8561, Japan

###### Abstract

We demonstrate that the location of a stable photon sphere (PS) in a compact
object is not always an edge such as the inner boundary of a black hole
shadow, whereas the location of an unstable PS is known to be the shadow edge
notably in the Schwarzschild black hole. If a static spherically symmetric
(SSS) spacetime has the stable outermost PS, the spacetime cannot be
asymptotically flat. A nondivergent deflection is caused for a photon
traveling around a stable PS, though a logarithmic divergent behavior is known
to appear in most of SSS compact objects with an unstable photon sphere. The
reason for the nondivergence is that the closest approach of a photon is
prohibited in the immediate vicinity of the stable PS when the photon is
emitted from a source (or reaches a receiver) distant from a lens object. The
finite gap size depends on the receiver and source distances from the lens as
well as the lens parameters. The mild deflection angle of light can be
approximated by an arcsine function. A class of SSS solutions in Weyl gravity
exemplify the nondivergent deflection near the stable outer PS.

###### pacs:

04.40.-b, 95.30.Sf, 98.62.Sb

##  I Introduction

Since the first measurement by Eddington and his collaborators Eddington , the
gravitational deflection of light has offered us a powerful tool for tests of
gravitational theories including the theory of general relativity as well as
for astronomical probes of dark matter. The Event Horizon Telescope (EHT) team
has recently succeeded in taking a direct image of the immediate vicinity of
the central black hole candidate of M87 galaxy EHT . In addition, the same
team has just reported measurements of linear polarizations around the same
black hole candidate EHT2021a , which has led to an estimation of the mass
accretion rate EHT2021b . These observations have increased our renewed
interest in the strong deflection of light in the strong gravity region.

The strong deflection of light by Schwarzschild black hole was pointed out by
Darwin Darwin . This phenomena is closely related with a photon sphere (PS)
that, with a horizon, features black holes and other compact objects Perlick ;
CVE ; Bozza2002 ; Hod2013 ; Sanchez ; DFR ; WLG ; SYG ; CGP ; Tsukamoto2017 ;
GW2016 ; Ohgami2015 ; Koga2016 ; Koga2018 ; Koga2019 ; Koga2020 ; Koga2021 ;
Tsukamoto2020 ; Shiromizu ; Yoshino2017 ; Yoshino2020a ; Yoshino2020b ;
Yang2020 ; Tsukamoto2021a ; Izumi2021 ; Tsukamoto2021b ; CG2016 ; Mishra2019 ;
CG2021 ; Guo2021 ; Konoplya2021 ; Gan2021 ; Ghosh2021 ; Siino2021a ;
Siino2021b . A photon surface for a less symmetric case is a generalization of
PSs CVE .

Many years later, Bozza showed that the strong deflection behavior in
Schwarzschild black hole can be well described as the logarithmic divergence
in the deflection angle Bozza2002 . Such a strong deflection as the
logarithmic divergence occurs also in other exotic objects such as wormholes.
For instance, Tsukamoto conducted several extensions of Bozza method for the
strong deflection of light Tsukamoto2017 ; Tsukamoto2020 , in which the
logarithmic behavior is shown to be a quite general feature for a static and
spherically symmetric (SSS) compact object that has a PS, e.g. Ellis
wormholes. The logarithmic behavior in the strong deflection for a finite-
distance receiver and source has been recently confirmed Takizawa2021  by
solving the exact gravitational lens that stands even for an asymptotically
nonflat spacetime Bozza2008 ; Takizawa2020a ; Takizawa2020b .

There exists a single PS outside the horizon of the Schwarzschild black hole,
while Ellis wormhole has a PS without horizons. For both cases, the number of
PSs is one. Tsukamoto obtained the logarithmic behavior for such a SSS
spacetime with a single PS, which was assumed to be unstable. In both of
Schwarzschild black hole and Ellis wormhole, there exists only the unstable
PS. Cunha et al. have recently proven that ultra-compact objects have an even
number of PSs, one of which is stable Cunha . Regarding the interesting
theorem, Hod has found an exception for horizonless spacetimes that possess no
stable PS Hod2018 . What happens, if the PS is stable and a light ray is
deflected around the stable PS?

The main purpose of the present paper is to study the deflection of light
around a PS when the PS is stable in a SSS spacetime. This situation has not
often been considered in detail e.g. Darwin ; Bozza2002 ; Tsukamoto2017 ;
Tsukamoto2021a ; Tsukamoto2021b , except for Hasse and Perlick (2002) HP  that
provided a theorem on a connection among three properties of (1) the presence
of a PS for a saddle point case and a stable case as well as an unstable one,
(2) the centrifugal force reversal and (3) infinitely many images in any SSS
spacetime. Hasse and Perlick (2002), however, did not calculate the deflection
angle of light, because they focused on the multiple imaging HP . We shall
show that, in stead of the logarithmic type of the strong deflection, a mild
deflection is caused near the stable PS.

This paper is organized as follows. In Section II, we reexamine the deflection
angle of light when there exists a stable outer PS in a SSS spacetime. In
Section III, we study how light is deflected in the presence of the stable
outer PS. It is shown that the deflection behavior is not too strong to make
the logarithmic divergence but mild enough to be approximated by an arcsine
function. As an example for such a mild deflection of light, we consider a
class of Weyl gravity model in Section IV. Section V concludes the present
paper. Throughout this paper, we use the unit of $G=c=1$.

##  II Deflection angle integral

###  II.1 SSS spacetime

We consider a SSS spacetime, for which the metric reads

$ds^{2}=-A(r)dt^{2}+B(r)dr^{2}+C(r)(d\theta^{2}+\sin^{2}\theta d\phi^{2}),$
(1)

where $A(r)>0$, $B(r)>0$, and $C(r)>0$ are assumed to be finite. If the SSS
spacetime possesses a horizon, we focus on the outside of the horizon. We do
not assume the asymptotic flatness of the spacetime. Actually, the present
paper considers an asymptotically nonflat case. A lemma on this issue is given
at the end of the present section.

###  II.2 Photon orbits

Without loss of generality, we can consider a photon orbit on the equatorial
plane $(\theta=\pi/2)$ because of the spherical symmetry of the spacetime. On
the equatorial plane in the SSS spacetime, a light ray has two constants of
motion. One is the specific energy $E\equiv A(r)\dot{t}$ and the other is the
specific angular momentum $L\equiv C(r)\dot{\phi}$, where the overdot denotes
the derivative with respect to the affine parameter $\lambda$ along the light
ray.

By using the two constants $E$ and $L$, the impact parameter of light becomes
$b=L/E$. Without loss of generality, we assume $b>0$.

In terms of $b$, the null condition $ds^{2}=0$ is rearranged as the orbit
equation

$\dot{r}^{2}+V(r)=0,$ (2)

where $V(r)$ is defined as

$V(r)\equiv-\frac{L^{2}F(r)}{B(r)C(r)}.$ (3)

Here,

$F(r)\equiv\frac{C(r)}{A(r)b^{2}}-1.$ (4)

The closest approach of a light ray is denoted as $r_{0}$, which satisfies
$V(r_{0})=0$ from the definition of $r_{0}$. Henceforth, evaluation at
$r=r_{0}$ is indicated by the subscript $0$. For instance, $V(r_{0})=0$ is
equivalent to $F_{0}=0$, where $F_{0}\equiv F(r_{0})$. By combining $F_{0}=0$
and Eq. (4), we obtain

$b=\sqrt{\frac{C_{0}}{A_{0}}}.$ (5)

This offers a relation between $b$ and $r_{0}$.

Following Hasse and Perlick HP , we use a particular form of the potential
$\tilde{V}(r)$ that is defined as

$\displaystyle\tilde{V}(r)\equiv\frac{A(r)}{C(r)},$ (6)

which is conformally invariant. Thereby Eq. (2) is rewritten as

$A(r)B(r)\dot{r}^{2}+L^{2}\tilde{V}(r)-E^{2}=0.$ (7)

This simplifies the analysis of the photon sphere and its linear perturbation,
because $V(r)$ depends on $b$ as well as $r$, whereas $\tilde{V}(r)$ does not
include $b$. The potential $\tilde{V}(r)$ plays the central role in the proof
of a theorem that clarifies a connection among three properties of the
presence of a PS, the centrifugal force reversal and infinitely many images in
any SSS spacetime HP .

For the later convenience, we write down the first and second derivatives of
$\tilde{V}(r)$,

$\displaystyle\tilde{V}^{{}^{\prime}}(r)=$
$\displaystyle-\frac{A(r)D(r)}{C(r)},$ (8)
$\displaystyle\tilde{V}^{{}^{\prime\prime}}(r)=$
$\displaystyle-\frac{A(r)}{C(r)}\left[D^{{}^{\prime}}(r)-D(r)^{2}\right],$ (9)

where the prime denotes the differentiation with respect to $r$, and the
functions $D$ is defined as

$\displaystyle D(r)$
$\displaystyle\equiv\frac{C^{\prime}(r)}{C(r)}-\frac{A^{\prime}(r)}{A(r)}.$
(10)

From Eq. (7), we obtain

$\displaystyle[A(r)B(r)]^{{}^{\prime}}\dot{r}^{2}+2A(r)B(r)\ddot{r}+L^{2}\tilde{V}^{{}^{\prime}}(r)=0.$
(11)

Eqs. (7) and (11) tell that a photon orbit is a circle if and only if
$r=r_{m}$ satisfies $\tilde{V}(r_{m})=(b_{m})^{-2}$ and
$\tilde{V}^{{}^{\prime}}(r_{m})=0$, where $r_{m}$ denotes the radius of the PS
and the subscript $m$ denotes the value at $r=r_{m}$. From
$\tilde{V}(r_{m})=(b_{m})^{-2}$, the impact parameter for the PS is obtained
as

$\displaystyle b_{m}=\sqrt{\frac{C_{m}}{A_{m}}},$ (12)

where Eq. (6) is used.

###  II.3 Classification of PS stability

We consider a small displacement $\delta r$ around the PS orbit,
$r=r_{m}+\delta r$. At the linear order in $\delta r$, Eq, (11) gives

$\displaystyle\frac{d^{2}}{d\lambda^{2}}(\delta r)$
$\displaystyle=-\frac{L^{2}\tilde{V}^{{}^{\prime\prime}}(r_{m})}{2A_{m}B_{m}}\delta
r$ $\displaystyle=\frac{L^{2}D^{{}^{\prime}}_{m}}{2B_{m}C_{m}}\delta r,$ (13)

where Eq. (9) and $D(r_{m})=0$ from $\tilde{V}^{{}^{\prime}}(r_{m})=0$ are
used.

The linear stability of the perturbed orbit is determined only by the sign of
$D_{m}^{{}^{\prime}}$, because $A(r)>0$, $B(r)>0$, $C(r)>0$ and $L\neq 0$. A
PS is stable if $D_{m}^{{}^{\prime}}<0$, whereas it is unstable for
$D_{m}^{{}^{\prime}}>0$.

From Eq. (10), we obtain

$\displaystyle
D^{{}^{\prime}}(r)=\frac{C^{\prime\prime}(r)}{C(r)}-\frac{A^{\prime\prime}(r)}{A(r)}-D(r)\left(\frac{C^{\prime}(r)}{C(r)}+\frac{A^{\prime}(r)}{A(r)}\right).$
(14)

On the PS, this is reduced to

$D^{\prime}_{m}=\frac{C_{m}^{\prime\prime}}{C_{m}}-\frac{A_{m}^{\prime\prime}}{A_{m}}.$
(15)

Bozza and Tsukamoto assumed $D^{\prime}_{m}>0$, which means an unstable PS
Bozza2002 ; Tsukamoto2017 . Henceforth, we focus on a stable case
$D^{\prime}_{m}<0$. Such an unusual case is realized in a class of Weyl
gravity model as shown in Section VI.

###  II.4 Total angle integral

From Eq. (2), we obtain

$\displaystyle\frac{dr}{d\phi}=\pm\sqrt{\frac{C(r)F(r)}{B(r)}},$ (16)

for which we choose the plus sign without loss of generality. Integrating this
from a source (denoted as S) to a receiver (denoted as R) leads to the total
change in the longitudinal angle.

$\displaystyle I_{F}(r_{S},r_{R},r_{0})$
$\displaystyle\equiv\sum_{i=S,R}\int_{r_{0}}^{r_{i}}\frac{dr}{\sqrt{\frac{C(r)F(r)}{B(r)}}}.$
(17)

Note that a conventional method discusses the total integral $I$ for the
asymptotic receiver and source ($r_{R}\to\infty$ and $r_{S}\to\infty$) e.g.
Darwin ; Bozza2002 ; Tsukamoto2017 , for which $I-\pi$ gives the deflection
angle of light. On the other hand, the present paper considers finite distance
between the receiver and source. In order to clarify this difference, we use
$I_{F}$ in stead of $I$. Rigorously speaking, $I_{F}-\pi$ is not the
deflection angle but $I_{F}$ is the dominant component of the deflection
angle. This issue is beyond the scope of this paper. See e.g. References
Ishihara2016 ; Ishihara2017 ; Ono2017 ; Ono2019U  on how to define
geometrically the deflection angle for the finite-distance receiver and
source.

Following Bozza and Tsukamoto Bozza2002 ; Tsukamoto2017 , we introduce a
variable as

$z\equiv 1-\frac{r_{0}}{r},$ (18)

to rewrite Eq. (17) as

$I_{F}(z_{S},z_{R},r_{0})=\sum_{i=S,R}\int_{0}^{z_{i}}f(z,r_{0})dz,$ (19)

where we define $z_{i}\equiv 1-r_{0}/r_{i}$ $(i=R,S)$,

$\displaystyle H(z,r_{0})$ $\displaystyle\equiv\frac{CF}{B}(1-z)^{4},$ (20)

and

$\displaystyle f(z,r_{0})$
$\displaystyle\equiv\frac{r_{0}}{\sqrt{H(z,r_{0})}}.$ (21)

In the next subsection, we shall examine the integrand in Eq. (19).

###  II.5 Analysis in the vicinity of PS

By noting $F_{0}=0$, the function $H$ is expanded around $r=r_{0}$ ($z=0$) as

$\displaystyle H(z,r_{0})=\sum_{n=1}^{\infty}c_{n}(r_{0})z^{n},$ (22)

where

$\displaystyle c_{1}(r_{0})$ $\displaystyle=\frac{C_{0}D_{0}r_{0}}{B_{0}},$
(23) $\displaystyle c_{2}(r_{0})$
$\displaystyle=\frac{C_{0}}{B_{0}}\left[D_{0}r_{0}\left(\frac{C^{\prime}_{0}}{C_{0}}-\frac{B^{\prime}_{0}}{B_{0}}-3\right)+\frac{r_{0}^{2}}{2}(D_{0}^{2}+D^{\prime}_{0})\right],$
(24) $\displaystyle c_{3}(r_{0})$
$\displaystyle=\frac{C_{0}}{B_{0}}\left[D_{0}r_{0}\left(\frac{3B^{\prime}_{0}r_{0}}{B_{0}}-\frac{3C^{\prime}_{0}r_{0}}{C_{0}}-\frac{B^{\prime}_{0}C^{\prime}_{0}r_{0}^{2}}{B_{0}C_{0}}-3\right)\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{r_{0}^{2}}{2}(D_{0}^{2}+D^{\prime}_{0})\left(\frac{C^{\prime}_{0}r_{0}}{C_{0}}-\frac{B^{\prime}_{0}r_{0}}{B_{0}}-2\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.+\frac{r_{0}^{3}}{6}\left(D_{0}^{3}+3D^{\prime}_{0}D_{0}+D^{\prime\prime}_{0}\right)\right].$
(25)

When the closest approach is located near the PS ($r_{0}=r_{m}$), Eqs. (23) -
(25) become

$\displaystyle c_{1}(r_{m})$ $\displaystyle=0,$ (26) $\displaystyle
c_{2}(r_{m})$ $\displaystyle=\frac{r_{m}^{2}C_{m}D_{m}^{\prime}}{2B_{m}},$
(27) $\displaystyle c_{3}(r_{m})$
$\displaystyle=\frac{r_{m}^{2}C_{m}}{2B_{m}}\left[D_{m}^{\prime}\left(\frac{C_{m}^{\prime}r_{m}}{C_{m}}-\frac{B_{m}^{\prime}r_{m}}{B_{m}}-2\right)+\frac{r_{m}D_{m}^{\prime\prime}}{3}\right].$
(28)

From Eqs. (26) - (28), we find

$H(z,r_{m})=c_{2}(r_{m})z^{2}+O(z^{3}).$ (29)

If the PS is unstable (stable), namely $D^{{}^{\prime}}_{m}>0$
($D^{{}^{\prime}}_{m}<0$), then, $c_{2}(r_{m})>0$ ($<0$). As pointed by
Tsukamoto Tsukamoto2017 , the angle integral $I_{F}$ by Eq. (17) is divergent
logarithmically. On the other hand, the stable PS case ($c_{2}(r_{m})<0$) is
investigated below in detail.

Before closing this section, we mention a relation of the emergence of the
stable outermost PS (SOPS) and the spacetime asymptotic flatness.

Lemma
If a SSS spacetime has the SOPS, the spacetime cannot be asymptotically flat,
for which $\tilde{V}^{{}^{\prime}}(r)=[A(r)/C(r)]^{{}^{\prime}}$ is positive
everywhere outside the SOPS.

Proof
We denote the radius of the SOPS as $r_{SOPS}$. At the location of the SOPS,
$\tilde{V}^{{}^{\prime}}(r_{SOPS})=0$ and
$\tilde{V}^{{}^{\prime\prime}}(r_{SOPS})>0$. There exist no PSs outside of the
SOPS, because the SOPS is the outermost PS. Therefore,
$\tilde{V}^{{}^{\prime\prime}}(r)>0$ for $r>r_{SOPS}$. This means that
$\tilde{V}^{{}^{\prime}}(r)$ is an increasing function of $r$ when
$r>r_{SOPS}$. Hence, $\tilde{V}^{{}^{\prime}}(r)>0$ for $r>r_{SOPS}$.

Here, we add an assumption that the spacetime were asymptotically flat. Then,
we can employ a coordinate system in which Eq. (1) approaches the Minkowski
metric in the polar coordinates asymptotically as $r\to\infty$. Namely,
$A(r)\to 1+O(1/r)$ and $C(r)\to r^{2}+O(r)$, which lead to
$A^{{}^{\prime}}(r)\to O(1/r^{2})$ and $C^{{}^{\prime}}(r)\to 2r+O(r^{0})$.
Thereby,

$\displaystyle\tilde{V}^{{}^{\prime}}(r)$
$\displaystyle=\left(\frac{A(r)}{C(r)}\right)^{{}^{\prime}}$
$\displaystyle=-\frac{2}{r^{3}}+O\left(\frac{1}{r^{4}}\right)$
$\displaystyle\to 0,$ (30)

in the limit as $r\to\infty$. This means that $\tilde{V}^{{}^{\prime}}(r)$ has
an extremum between $r_{SOPS}$ and $r=\infty$, because of its continuity. This
contradicts with that $\tilde{V}^{{}^{\prime}}(r)$ is an increasing function
for $r>r_{SOPS}$. Therefore, the spacetime is not asymptotically flat. A proof
of the lemma is thus completed.

According to Reference Guo2021 , if an axisymmetric, stationary and
asymptotically flat spacetime posses light rings (LRs), the outest LR is
unstable. This means that if a SSS spacetime with PSs is asymptotically flat,
the outest PS is unstable. The contraposition of this is that, if the outest
PS in a SSS spacetime with PSs is stable, the spacetime is not asymptotically
flat. This proves a part of the above lemma but tells nothing about the
positivity of $\tilde{V}^{{}^{\prime}}(r)$ for $r>r_{SOPS}$.

On the other hand, it is clear that the asymptotic flatness is allowed, if the
outermost PS (OPS) is unstable ($\tilde{V}^{{}^{\prime\prime}}(r_{OPS})<0$).
In the rest of this paper, we consider that a receiver and source are located
outside of a stable outer PS. However, it is not specified below whether or
not it is the outermost.

Before closing the section, let us briefly mention the stability of a
spacetime that admits a stable PS. Horizonless ultracompact objects with a
stable PS may suffer from instabilities due to slowly (at most
logarithmically) decaying of perturbations leading to the formation of a
trapped surface Cardoso2014 ; Keir ; Cunha ; AMY . In section VI, therefore,
we consider a black hole model with a stable PS in Weyl gravity. On the other
hand, the above lemma may suggest another possibility compatible with the
instability arguments in References Cardoso2014 ; Keir ; Cunha ; AMY . One
such candidate is a compact object without a black hole horizon but with a
stable PS and a cosmological horizon that is consistent with the asymptotic
nonflatness. Its stability issue is beyond the scope of the present paper.

##  III Mild deflection near the stable outer PS

###  III.1 Stability classification of the outer PS

In the neighborhood of the closest approach $z\sim 0$, higher order terms of
$H(z,r_{m})$ in Eq. (22) are negligible compared with $z$ and $z^{2}$ terms.
The dominant part of $f(z,r_{0})$ for $r_{0}\sim r_{m}$ thus becomes

$f_{D}(z,r_{0})\equiv\frac{r_{0}}{\sqrt{h(z)}},$ (31)

where $h(z)$ is defined as

$\displaystyle h(z)=c_{1}(r_{0})z+c_{2}(r_{0})z^{2}.$ (32)

Associated with $f_{D}(z,r_{0})$, we define the dominant part of the total
angle integral as

$I_{FD}(z_{R},z_{S},r_{0})\equiv\sum_{i=S,R}\int_{0}^{z_{i}}f_{D}(z,r_{0})dz.$
(33)

Before starting calculations of the angle integral, first we investigate a
photon orbit in the stable PS case. If $c_{1}(r_{0})<0$, $h(z)$ is always
negative for $z>0$. This is in contradiction with the nonnegativity of $h(z)$.
Hence, the case of $c_{1}(r_{0})<0$ is discarded. If $c_{1}(r_{0})=0$, then,
$h(z)=c_{2}(r_{0})z^{2}$. The nonnegativity of $h(z)$ together with
$c_{2}(r_{0})<0$ admits only $z=0$. This orbit is a circle. We do not discuss
this special case any more. For the last case $c_{1}(r_{0})>0$, the
nonnegativity of $h(z)$ provides a nontrivial situation; the allowed region
for $z$ is

$\displaystyle 0\leq z\leq-\frac{c_{1}(r_{0})}{c_{2}(r_{0})}.$ (34)

This means that only bound orbits are allowed, whereas the scattering orbits
are prohibited. Note that $h(z)=0$ does not make the integral in Eq. (33)
divergent. Indeed, $h(z)\sim c_{1}(r_{0})z$ near a point $z=0$ and $h(z)\sim-
c_{1}(r_{0})[z+c_{1}(r_{0})/c_{2}(r_{0})]$ near a point $z\sim-
c_{1}(r_{0})/c_{2}(r_{0})$. At the two points, therefore, the integral of
$h(z)$ is not divergent as $\sim\int z^{-1/2}dz\sim 2\sqrt{z}$. The two points
$z=0$ and $z\sim-c_{1}(r_{0})/c_{2}(r_{0})$ are the periastron and apastron,
respectively.

###  III.2 Angle integral for the stable outer PS

Henceforth, we focus on the case of $c_{1}(r_{0})>0$, for which the receiver
and source positions satisfy

$0<z_{i}\leq-\frac{c_{1}(r_{0})}{c_{2}(r_{0})}.$ (35)

Eq. (33) is integrated as

$\displaystyle I_{FD}(z_{R},z_{S},r_{0})=$
$\displaystyle\frac{r_{0}}{\sqrt{-c_{2}(r_{0})}}$
$\displaystyle\times\left[\pi-\sum_{i=S,R}\arcsin\left(1+\frac{2c_{2}(r_{0})}{c_{1}(r_{0})}z_{i}\right)\right].$
(36)

It follows that the arcsine term in Eq. (36) is well defined for the allowed
region by Eq. (35), because

$-1\leq 1+\frac{2c_{2}(r_{0})}{c_{1}(r_{0})}z_{i}<1.$ (37)

###  III.3 $I_{F}$ in terms of $b$

In most of lens studies including References Bozza2002 ; Tsukamoto2017 , it is
convenient to express the deflection angle in terms of the impact parameter in
stead to $r_{0}$, mainly because $b\sim D_{L}\theta$ for the lens distance
$D_{L}$ and the image angle direction $\theta$.

Therefore, we look for an approximate expression of $b$. Near the PS,
$c_{1}(r_{0})$ and $b(r_{0})$ are Taylor-expanded as

$\displaystyle c_{1}(r_{0})$
$\displaystyle=\frac{C_{m}r_{m}D^{\prime}_{m}}{B_{m}}(r_{0}-r_{m})+O((r_{0}-r_{m})^{2}),$
(38) $\displaystyle b(r_{0})$
$\displaystyle=b_{m}+\frac{b_{m}D^{\prime}_{m}}{4}(r_{0}-r_{m})^{2}+O((r_{0}-r_{m})^{3}).$
(39)

From Eq. (38), we find

$\displaystyle r_{0}<r_{m},$ (40)

where $c_{1}(r_{0})>0$ and $D^{{}^{\prime}}_{m}<0$ are used. This inequality
means that the light ray passes by the slight inside of the PS. This is
because the PS is stable. This unusual behavior of the photon orbit implies
also

$\displaystyle b<b_{m},$ (41)

from Eq. (39). See Figure 1 for the photon orbit behavior near the PS.

Figure 1:  Schematic figure of photon orbits near a PS. The PS is denoted by a
solid circle with radius $r_{m}$. The dotted circle has radius $r_{0}$, which
is slightly smaller than $r_{m}$. Namely, $r_{m}-r_{0}$ is taken as a
perturbation around the PS. The initial position of the photon is $r_{0}$, at
which the initial radial velocity is vanishing. The thick dashed blue (in
color) arrow denotes a photon motion when the PS is unstable. The orbital
deviation grows because of being unstable. Therefore, it is difficult for the
inner photon to escape to a far region. This is why the earlier papers focused
on $r_{0}>r_{m}$ when an unstable PS is assumed Bozza2002 ; Tsukamoto2017 . On
the other hand, the thick solid red (in color) arrow denotes a photon motion
when the PS is stable. The orbital deviation gets smaller because of being
stable. It is thus possible that the inner photon ($r_{0}<r_{m}$) crosses the
PS to reach a distant observer.

By combining Eqs. (38) and (39), we obtain, near the PS, an approximate
relation between $c_{1}(r_{0})$ and $b$ as

$\displaystyle
c_{1}(r_{0})\approx\frac{2r_{m}C_{m}\sqrt{-D^{\prime}_{m}}}{B_{m}}\left(1-\frac{b}{b_{m}}\right)^{\frac{1}{2}}.$
(42)

By substituting this into Eq. (36), we obtain

$\displaystyle I_{FD}(z_{R},z_{S},b)$
$\displaystyle\simeq\frac{r_{m}}{\sqrt{-c_{2}(r_{m})}}$
$\displaystyle~{}\times\left[\pi-\sum_{i=S,R}\arcsin\left\\{1-\frac{r_{m}\sqrt{-D^{\prime}_{m}}}{2}\left(1-\frac{b}{b_{m}}\right)^{-\frac{1}{2}}z_{i}\right\\}\right].$
(43)

By direct calculations for a photon traveling near the PS, Eq. (35) is
rearranged as

$0<z_{i}\leq\frac{4}{r_{m}\sqrt{-D^{\prime}_{m}}}\left(1-\frac{b}{b_{m}}\right)^{1/2},$
(44)

for which the arcsine function in Eq. (43) is well defined.

In order to obtain Eq. (43) as an approximate estimation of $I_{F}$ in terms
of $b$, we have used the Taylor series expansion as Eq. (39), which is valid
if $b_{m}>b_{m}|D^{{}^{\prime}}_{m}|(r_{0}-r_{m})^{2}/4$. This requires that
the closest approach of light and the PS are close enough to satisfy

$\displaystyle|r_{0}-r_{m}|<\frac{2}{\sqrt{|D^{{}^{\prime}}_{m}|}}.$ (45)

Yet, a lower bound on $|r_{0}-r_{m}|$ exists in the neighborhood of the PS as
shown below.

###  III.4 Discontinuity between the closest approach and the stable PS

Eq. (44) suggests a proposition on the existence of a gap between the allowed
closest approach and the stable PS.

Proposition
In a SSS spacetime which possess a stable outer photon sphere, the closest
approach of a photon from a source (or to a receiver) located at finite
distance from the lens object is not allowed in the infinitesimal neighborhood
of the stable PS.

Proof
This proposition can be proven by contradiction as follows. We consider a
receiver (or a source) at finite distance from the lens, namely $r_{i}>r_{0}$,
which leads to $z_{i}>0$. We assume the closest approach of a photon orbit
were in the infinitesimal neighborhood of the stable PS.

On the stable PS, $F_{m}=F^{{}^{\prime}}_{m}=0$, while
$F^{{}^{\prime\prime}}_{m}<0$. The last inequality comes from the unstable
condition $D^{{}^{\prime}}_{m}<0$. The function $F$ near the PS is thus
Taylor-expanded as

$\displaystyle
F_{0}=\frac{1}{2}F^{{}^{\prime\prime}}_{m}(r_{0}-r_{m})^{2}+O\left((r_{0}-r_{m})^{3}\right).$
(46)

From Eq. (20), $H$ for $r_{0}\approx r_{m}$ is expanded as

$\displaystyle
H=\frac{C_{m}F^{{}^{\prime\prime}}_{m}}{2B_{m}}(1-z)^{4}(r_{0}-r_{m})^{2}+O\left((r_{0}-r_{m})^{3}\right),$
(47)

where Eq. (46) is used.

By using $B_{m}>0$, $C_{m}>0$ and $F^{{}^{\prime\prime}}_{m}<0$ for Eq. (47),
we find $H<0$ if $r_{0}-r_{m}$ is sufficiently small. $H<0$ contradicts with
the existence of the photon orbit. Our proof is finished.

The above proposition prohibits the closest approach in the infinitesimal
neighborhood of the stable PS. However, it does not tell about a size of the
gap between the allowed closest approach and the stable PS. In order to
discuss the gap size, we use Eq. (44), which is rewritten as

$\displaystyle\frac{b}{b_{m}}\leq
1-\frac{r_{m}^{2}z_{i}^{2}(-D_{m}^{{}^{\prime}})}{16}.$ (48)

Note that Eq. (44) is based on a quadratic approximation up to $z^{2}$.

For finite $z_{i}$ and $D_{m}^{{}^{\prime}}<0$, Eq. (48) demonstrates that $b$
is not allowed in the infinitesimal neighborhood of $b_{m}$. From Eq. (48),
the upper bound on $b$ is

$\displaystyle
b_{max}=b_{m}-\frac{b_{m}r_{m}^{2}z_{L}^{2}(-D_{m}^{{}^{\prime}})}{16}.$ (49)

The gap size $\Delta b\equiv b_{m}-b_{max}$ is thus given by

$\displaystyle\Delta
b=\frac{b_{m}r_{m}^{2}z_{L}^{2}(-D_{m}^{{}^{\prime}})}{16},$ (50)

where $z_{L}$ is the larger one of $z_{R}$ and $z_{S}$, namely the value of
$z$ for the more distant one of the receiver and source.

A separate treatment is needed for the marginal PSs ($D^{\prime}_{m}=0)$ e.g.
CT ; Tsukamoto2020 .

###  III.5 The dominant part and the remainder

Before closing this section, we shall confirm that $I_{FR}\equiv I_{F}-I_{FD}$
is really the remainder. It is written simply as

$I_{FR}(z_{R},z_{S},r_{0})\equiv\sum_{i=R,S}\int_{0}^{z_{i}}f_{R}(z,r_{0})dz,$
(51)

where $f_{R}(z,r_{0})\equiv f(z,r_{0})-f_{D}(z,r_{0})$.

From Eq. (51), we obtain, for $z_{i}<1$,

$\displaystyle
I_{FR}(z_{R},z_{S},r_{0})=\sum_{i=R,S}[f_{R}(z_{i},r_{0})z_{i}+O(z_{i}^{2})].$
(52)

For $r_{0}\sim r_{m}$, $f(z,r_{m})\sim f_{D}(z,r_{m})$ when $z<1$. Hence, by
using the Taylor-expansion method, one may think that
$f_{R}(z,r_{m})=f(z,r_{m})-f_{D}(z,r_{m})=O(z)$. However, let us more
carefully examine an asymptotic expansion of $f_{R}(z,r_{0})$ for $z<1$
footnote-expansion . By straightforward calculations, the asymptotic expansion
of $f(z,r_{0})$ for small $z$ is obtained as

$\displaystyle f(z,r_{0})=$
$\displaystyle\frac{r_{0}}{\sqrt{zc_{1}(r_{0})}}-\frac{r_{0}c_{2}(r_{0})}{2c_{1}(r_{0})^{3/2}}\sqrt{z}+O(z^{3/2}).$
(53)

Similarly, the asymptotic expansion of $f_{D}(z,r_{0})$ is

$\displaystyle f_{D}(z,r_{0})=$
$\displaystyle\frac{r_{0}}{\sqrt{zc_{1}(r_{0})}}-\frac{r_{0}c_{2}(r_{0})}{2c_{1}(r_{0})^{3/2}}\sqrt{z}+O(z^{3/2}).$
(54)

By bringing together Eqs. (53) and (54), we find

$\displaystyle f_{R}(z,r_{0})=O(z^{3/2}).$ (55)

By using this for Eq. (52), we find

$\displaystyle I_{FR}(z_{R},z_{S},r_{0})=O(z_{R}^{5/2},z_{S}^{5/2}).$ (56)

On the other hand, Eq. (36) is expanded, for small $z_{R}$ and $z_{S}$, as

$\displaystyle I_{FD}(z_{R},z_{S},r_{0})$
$\displaystyle=\frac{2r_{0}}{\sqrt{c_{1}(r_{0})}}(\sqrt{z_{R}}+\sqrt{z_{S}})+O(z_{R},z_{S}),$
(57)

where $\arcsin(1-\varepsilon)=\pi/2-\sqrt{2\varepsilon}+O(\varepsilon)$ is
used.

For $z_{R},z_{S}\ll 1$, Eqs. (56) and (57) lead to
$I_{FR}(z_{R},z_{S},r_{0})\ll I_{FD}(z_{R},z_{S},r_{0})$. Therefore, we can
safely say that $I_{FD}$ is the dominant part and $I_{FR}$ is the remainder.

##  IV Example: A class of Weyl gravity model

In the Weyl gravity model, Mannheim and Kazanas found a class of SSS solutions
MK . The metric reads

$ds^{2}=-g(r)dt^{2}+g(r)^{-1}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2}),$ (58)

where

$g(r)=1-3\beta\gamma-\frac{\beta(2-3\beta\gamma)}{r}+\gamma r-\kappa r^{2},$
(59)

and we consider $g(r)>0$, namely the outside of the horizon.

The allowed region for the existence of the stable outer PS is TH

$\displaystyle\beta>0,$ (60) $\displaystyle\gamma<0,$ (61)
$\displaystyle\kappa<0.$ (62)

In this section, we focus on this parameter region.

Figure 2:  Effective potential $\tilde{V}(r)$ in the Weyl gravity model. The
solid blue (in color) curve corresponds to $\gamma=-0.5$ and $\kappa=-0.05$ in
the unit of $\beta=1$, which leads to the SOPS at $r_{m}=7$. From
$(b_{m})^{-2}=\tilde{V}(r_{m})=0.0193878$, $b_{m}=7.18184$ is obtained.
$\tilde{V}(r)$ does not vanish for large $r$, because the spacetime is not
asymptotically flat.

We solve $D(r)=0$ to find two roots as $r=3\beta$ and $r=3\beta-2/\gamma$.
They are corresponding to PSs. One is a stable PS and the other is an unstable
one.

In the present case of $\gamma<0$, the stable outer PS is located at

$\displaystyle r_{m}=3\beta-\frac{2}{\gamma},$ (63)

because

$\displaystyle D_{m}^{\prime}=$ $\displaystyle-\frac{2}{r_{m}^{2}g(r_{m})}$
$\displaystyle=$ $\displaystyle-\frac{2b_{m}^{2}}{r_{m}^{4}}$ $\displaystyle<$
$\displaystyle 0.$ (64)

The other root $r=3\beta$ is the radius of the unstable inner PS.

Note that $r_{m}$ is larger than $3\beta$ because of $\gamma<0$. See Figure 2
for the effective potential $V(r)$ in the Weyl gravity model with the stable
outer PS.

There is a constraint on $\kappa$ as footnote-Horne

$\displaystyle\kappa<-\frac{\gamma^{2}(1-\beta\gamma)}{(2-3\beta\gamma)^{2}}.$
(65)

Here, $b_{m}$ denotes the critical impact parameter corresponding to the
stable PS. It is obtained as

$\displaystyle
b_{m}=\frac{2-3\beta\gamma}{\sqrt{-\kappa(2-3\beta\gamma)^{2}-\gamma^{2}(1-\beta\gamma)}},$
(66)

where the inside of the square root is always positive when $\kappa$ satisfies
Eq. (65).

Figure 3:  The integrands of Eqs. (19) and (33). The solid blue (in clor)
curve denotes $f$ ($=d\phi/dz$) in the integral $I_{F}$, while the dashed red
(in color) curve denotes $f_{D}$ in $I_{FD}$. The horizontal axis means $z$.
We assume $\kappa=-0.05$ in the unit of $\beta=1$ that is roughly
corresponding to the conventional unit as mass = 1 in the Schwarzschild or
Kottler spacetime. The top and bottom panels assume $\gamma=-0.5$ and
$\gamma=-0.1$, respectively, each of which leads to $r_{m}=7.0$ and
$b_{m}=7.18185$, and $r_{m}=23$ and $b_{m}=4.56815$, respectively. In both
cases, the solid and dashed curves are overlapped better for smaller $z$. For
instance, the difference between them is about 10 percents at $z\sim 0.2$. A
significant deviation at large $z$ is due to a departure from the quadratic
approximation of $H(z)$ in $z$.

From Eqs. (43) and (56), we obtain

$\displaystyle I_{FD}(z_{S},z_{R},b)$
$\displaystyle\simeq\pi-\sum_{i=S,R}\arcsin\left\\{1-\frac{b_{m}}{\sqrt{2}r_{m}}\left(1-\frac{b}{b_{m}}\right)^{-1/2}z_{i}\right\\},$
(67) $\displaystyle I_{FR}(z_{S},z_{R},b)$
$\displaystyle\simeq-\beta\gamma\left(\frac{b_{m}}{2\sqrt{2}r_{m}}\right)^{3/2}\left(1-\frac{b}{b_{m}}\right)^{-3/4}\sum_{i=S,R}z_{i}^{5/2},$
(68)

where straightforward calculations at $O(z^{5/2})$ are done in Eq. (56).

When $z_{T}\equiv z_{R}=z_{S}\ll(1-b/b_{m})^{1/2}$, Eq. (67) provides an
approximate expression of the deflection angle as

$\displaystyle\Delta\phi=$ $\displaystyle I_{F}$ $\displaystyle=$
$\displaystyle
2\sqrt{\frac{\sqrt{2}b_{m}}{r_{m}}\left(1-\frac{b}{b_{m}}\right)^{-1/2}z_{T}}$
$\displaystyle+O\left(\frac{b_{m}}{r_{m}}\left(1-\frac{b}{b_{m}}\right)^{-1/2}z_{T}\right).$
(69)

where we use $\arcsin(1-\varepsilon)=\pi/2-\sqrt{2\varepsilon}+O(\varepsilon)$
for $\varepsilon<1$.

Figure 4:  The total angle integrals $I_{F}$ and $I_{FD}$. The solid blue (in
color) curve denotes the present approximation of $I_{FD}$ by Eq. (67) and the
dashed red (in color) curve denotes numerical calculations of $I_{F}$ in Eq.
(19), where we assume the same values for $\kappa$ and $\gamma$ in Figure 3.
For the simplicity, we choose $z\equiv z_{R}=z_{S}$, which is denoted by the
horizontal axis. For $z\ll 1$, the two curves are very close to each other,
because the quadratic approximation works well especially near $z\sim 0$. The
two curves are close to each other especially for smaller $z$. For instance,
the difference between them is roughly 20 percents at $z\sim 0.2$. A
significant deviation at large $z$ reflects a departure from the quadratic
approximation of $H(z)$ in $z$.

It is natural that $I_{FR}$ is much smaller than $I_{FD}$, as discussed in
Section III. Eq. (67) shows the mild deflection in terms of the arcsine
function. See Figure 4 for a comparison of Eq. (67) and numerical calculations
for $I_{FD}$. As shown in Figure 3, $f$ are dependent on rather strongly on
$\gamma$. The difference between the numerical $f$ and approximate one becomes
significant as $z$ is larger. On the other hand, the closest approach and its
vicinity make a dominant contribution to the total angle integral $I_{F}$. In
Figure 4, therefore, $I_{F}$ shows much weaker dependence on $\gamma$ and $z$.

Finally, we discuss the gap size. By using Eq. (64) for Eq. (50), we obtain

$\displaystyle\Delta b=\frac{b_{m}^{3}z_{L}^{2}}{8r_{m}^{2}}.$ (70)

For $b_{m}\sim r_{m}$, furthermore, it becomes simply

$\displaystyle\Delta b\approx\frac{1}{8}b_{m}z_{L}^{2}.$ (71)

In a SSS spacetime, the photon sphere was thought to be always an edge such as
the inner boundary of a black hole shadow. However, Eqs. (70) and (71) provide
a counterexample, when the PS is stable.

##  V Conclusion

We demonstrated that the location of a stable PS in a compact object is not
always an edge such as the inner boundary of a black hole shadow. We showed
also that a SSS spacetime cannot be asymptotically flat, when the SOPS exists
in the spacetime.

We proved a proposition that the closest approach of a photon is prohibited in
the immediate vicinity of the stable PS when the photon is emitted from a
source (or reaches a receiver) distant from a lens object. We discussed the
gap size. Because of the existence of the gap, the mild deflection is caused
for a photon traveling around the stable PS which exists in a class of SSS
spacetimes.

Finally, we used a class of SSS solutions in Weyl gravity in order to
exemplify the mild deflection near the stable outer PS. It would be
interesting to examine whether a light ray is mildly deflected around a stable
photon surface in a less symmetric or dynamical spacetime. It is left for
future.

###### Acknowledgements.

We would like to thank Mareki Honma for the conversations on the EHT method
and technology. We wish to thank Keita Takizawa, Naoki Tsukamoto, Chul-Moon
Yoo, Kenichi Nakao, Tomohiro Harada and Masaya Amo for the useful discussions.
We thank Yuuiti Sendouda, Ryuichi Takahashi, Masumi Kasai, Kaisei Takahashi,
and Yudai Tazawa for the useful conversations. This work was supported in part
by Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for
Scientific Research, No. 20K03963 (H.A.), in part by Ministry of Education,
Culture, Sports, Science, and Technology, No. 17H06359 (H.A.).

## References

  * (1) F. W. Dyson, A. S. Eddington, and C. Davidson, Phil. Trans. R. Soc. A 220, 291 (1920).
  * (2) K. Akiyama et al. (Event Horizon Telescope Collaboration), Astrophys. J. 875, L1 (2019); Astrophys. J. 875, L2 (2019); Astrophys. J. 875, L3 (2019); Astrophys. J. 875, L4 (2019); Astrophys. J. 875, L5 (2019); Astrophys. J. 875, L6 (2019).
  * (3) K. Akiyama et al. (Event Horizon Telescope Collaboration), Astrophys. J. 910, L12 (2021).
  * (4) K. Akiyama et al. (Event Horizon Telescope Collaboration), Astrophys. J. 910, L13 (2021).
  * (5) C. Darwin, Proc. Roy. Soc. Lond. A 249,180 (1959).
  * (6) V. Bozza, Phys. Rev. D 66, 103001 (2002).
  * (7) N. Tsukamoto, Phys. Rev. D 95, 064035 (2017).
  * (8) N. Tsukamoto, Phys. Rev. D 102, 104029 (2020).
  * (9) N. Tsukamoto, Phys. Rev. D 104, 064022 (2021).
  * (10) N. Tsukamoto, Phys. Rev. D. 104, 124016 (2021).
  * (11) V. Perlick, Living Rev. Relativity 7, 9 (2004).
  * (12) N. G. Sanchez, Phys. Rev. D 18, 1030 (1978).
  * (13) C. M. Claudel, K. S. Virbhadra, and G. F. R. Ellis, J. Math. Phys. 42, 818 (2001).
  * (14) S. Hod, Phys. Lett. B 727, 345 (2013).
  * (15) Y. Decanini, A. Folacci, and B. Raffaelli, Phys. Rev. D 81, 104039 (2010).
  * (16) S. W. Wei, Y. X. Liu, and H. Guo, Phys. Rev. D 84, 041501(R) (2011).
  * (17) I. Z. Stefanov, S. S. Yazadjiev, and G. G. Gyulchev, Phys. Rev. Lett. 104, 251103 (2010).
  * (18) M. Cvetic, G. W. Gibbons, and C. N. Pope, Phys. Rev. D 94, 106005 (2016).
  * (19) G. W. Gibbons, and C. M. Warnick, Phys. Lett. B 763, 169 (2016).
  * (20) T. Ohgami, and N. Sakai, Phys. Rev. D 91, 124020 (2015).
  * (21) Y. Koga and T. Harada, Phys. Rev. D 94, 044053 (2016).
  * (22) Y. Koga and T. Harada, Phys. Rev. D 98, 024018 (2018).
  * (23) Y. Koga, Phys. Rev. D 99, 064034 (2019).
  * (24) Y. Koga, Phys. Rev. D 101, 104022 (2020).
  * (25) Y. Koga, T. Igata, and K. Nakashi, Phys. Rev. D 103, 044003 (2021).
  * (26) T. Shiromizu, Y. Tomikawa, K. Izumi, and H. Yoshino, Prog. Theor. Exp. Phys. 2017, 033E01 (2017).
  * (27) H. Yoshino, K. Izumi, T. Shiromizu, and Y. Tomikawa, Prog. Theor. Exp. Phys. 2017, 063E01 (2017).
  * (28) H. Yoshino, K. Izumi, T. Shiromizu, and Y. Tomikawa, Prog. Theor. Exp. Phys. 2020, 023E02 (2020).
  * (29) H. Yoshino, K. Izumi, T. Shiromizu, and Y. Tomikawa, Prog. Theor. Exp. Phys. 2020, 053E01 (2020).
  * (30) R. Yang, and H. Lu, Eur. Phys. J. C 80, 949 (2020).
  * (31) K. Izumi, Y. Tomikawa, T. Shiromizu, and H. Yoshino, Prog. Theor. Exp. Phys. 2021, 083E02 (2021).
  * (32) C. Cederbaum, and G. J. Galloway, Class. Quantum Grav. 33, 075006 (2016).
  * (33) A. K. Mishra, S. Chakraborty, and S. Sarkar, Phys. Rev. D 99, 104080 (2019).
  * (34) C. Cederbaum, and G. J. Galloway, J. Math. Phys. 62, 032504 (2021).
  * (35) M. Guo, and S. Gao, Phys. Rev. D 103, 104031 (2021).
  * (36) R. A. Konoplya, and A. Zhidenko, Phys. Rev. D 103, 104033 (2021).
  * (37) Q. Gan, P. Wang, H. Wu, and H. Yang, Phys. Rev. D 104, 024003 (2021).
  * (38) R. Ghosh, and S. Sarkar, Phys. Rev. D 104, 044019 (2021).
  * (39) M. Siino, Class. Quantum Grav. 38, 025005 (2021).
  * (40) M. Siino, ArXiv:2107.06551.
  * (41) K. Takizawa, and H. Asada, Phys. Rev. D 103, 104039 (2021).
  * (42) V. Bozza, Phys. Rev. D 78, 103005 (2008).
  * (43) K. Takizawa, T. Ono, and H. Asada, Phys. Rev. D 101, 104032 (2020).
  * (44) K. Takizawa, T. Ono, and H. Asada, Phys. Rev. D 102, 064060 (2020).
  * (45) P. V. P. Cunha, E. Berti and C. A. R. Herdeiro, Phys. Rev. Lett. 119, 251102 (2017).
  * (46) S. Hod, Phys. Lett. B 776, 1 (2018).
  * (47) W. Hasse and V. Perlick, Gen. Rel. Grav. 34, 415, (2002).
  * (48) A. Ishihara, Y. Suzuki, T. Ono, T. Kitamura and H. Asada, Phys. Rev. D 94, 084015 (2016).
  * (49) A. Ishihara, Y. Suzuki, T. Ono and H. Asada, Phys. Rev. D 95, 044017 (2017).
  * (50) T. Ono, A. Ishihara, and H. Asada, Phys. Rev. D 96, 104037 (2017).
  * (51) T. Ono, and H. Asada, Universe, 5(11), 218 (2019).
  * (52) V. Cardoso, L. C. B. Crispino, C. F. B. Macedo, H. Okawa, and P. Pani, Phys. Rev. D 90, 044069 (2014).
  * (53) J. Keir, Class. Quantum Grav. 33, 135009 (2016).
  * (54) A. Addazi, A. Marciano, and N. Yunes, Eur. Phys. J. C 80, 36 (2020).
  * (55) T. Chiba, and M. Kimura, Prog. Theor. Exp. Phys. 4, 043E01 (2017).
  * (56) There is a prohibited region around the stable PS as shown in Section IIID. Because the existence of this gap region, the functions $H$ and $f$ are not differentiable at $r=r_{m}$. As a result, $H$ and $f$ cannot be Taylor-expanded around $r_{m}$. In stead of Taylor expansions, therefore, we look for asymptotic expansions of $I_{F}$, $I_{FD}$ and $I_{FR}$ through $f$, $f_{D}$ and $f_{R}$.
  * (57) P. D. Mannheim, and D. Kazanas, Astrophys. J. 342, 635 (1989).
  * (58) Turner and Horne realized that a bound orbit of a photon is possible in this parameter region for the SSS solution in Weyl gravity TH .
  * (59) G. E. Turner, and K. Horne, Class. Quantum Grav. 37, 095012 (2020).


# Assessing the role of quantum effects in 2D heterophase MoTe2 field effect
transistors.

Line Jelver<EMAIL_ADDRESS>CAMD, Dept. of Physics, Technical University of
Denmark, Bldg. 309, DK-2800 Kongens Lyngby, Denmark.  Synopsys QuantumATK,
Fruebjergvej 3, PostBox 4, DK-2100 Copenhagen, Denmark.  Center for Nano
Optics, University of Southern Denmark, DK-5230 Odense M, Denmark.  Ole Hansen
DTU Nanolab, Technical University of Denmark, Bldg. 347, DK-2800 Kongens
Lyngby, Denmark.  Karsten Wedel Jacobsen  CAMD, Dept. of Physics, Technical
University of Denmark, Bldg. 309, DK-2800 Kongens Lyngby, Denmark.

###### Abstract

The two-dimensional transition metal dichalcogenides (TMDs) have been proposed
as candidates for the channel material in future field effect transistor
designs. The heterophase design which utilizes the metallic T- or T’ phase of
the TMD as contacts to the semiconducting H phase channel has shown promising
results in terms of bringing down the contact resistance of the device. In
this work, we use ab-initio calculations to demonstrate how atomic-scale and
quantum effects influence the ballistic transport properties in such
heterophase transistors with channel lengths up to 20 nm. We investigate how
the charge transfer depends on the carrier density both in T’-H MoTe2 Schottky
contacts and planar T’-H-T’ MoTe2 transistors. We find that the size of the
Schottky barrier and the charge transfer is dominated by the local atomic
arrangements at the interface and the doping level. Furthermore, two types of
quantum states have a large influence on the charge transport; interface
states and standing waves in the semiconductor due to quantum confinement. We
find that the latter can be associated with rises in the current by more than
an order of magnitude due to resonant tunneling. Our results demonstrate the
quantum mechanical nature of these 2D transistors and highlight several
challenges and possible solutions for achieving a competitive performance of
such devices.

Field effects transistors (FETs) have through several decades been
continuously optimized and down-scaled. Today, quantum effects have become one
of the limiting factors in the device performance. Tunneling effects between
the source and drain as well as through the gate stack make it challenging to
continue the down-scaling. Further development of the energy efficiency and
computational power is therefore to a larger and larger degree achieved by
optimizing the transistor design and materials.[1]

The 2D transition metal dichalcogenides (TMDs) represent a promising group of
materials which could substitute silicon as the semiconducting channel of the
device. These two dimensional semiconductors have the advantage that they, in
contrast to silicon, can be synthesized down to the monolayer limit while
retaining a high mobility.[2, 3] This allows for the ultimate down-scaling of
the transistor channel with a typical monolayer thickness of only 0.7 nm.
However, two large challenges need to be resolved before these materials can
become a viable choice for future transistor designs. Firstly, a scalable
method for area selective doping of 2D materials is not currently
available.[4] This means that the carrier density in a 2D semiconductor cannot
be controlled in a reproducible manner. Secondly, the interface between the 2D
channel and the 3D source and drain contacts is difficult to fabricate. Metal
deposition on a 2D material can be very destructive and impurities and defects
in the interface result in large contact resistances.[5, 6] A promising
solution to this is to exploit that the TMDs exists in different phases which
exhibit either metallic or semiconducting properties. The group-IV TMDs are
generally most stable in the semiconducting 2H phase but in MoTe2, the
metallic 1T’ phase is meta-stable[7] and can be selectively grown[8, 9, 10,
11, 12, 13, 14] or induced from a 2H sample.[15, 16] The MoTe2 heterophase
transistor design relies on this fact and utilizes the 1T’ phase of the TMD as
the source and drain contact. This approach reduces the contact resistance
significantly since charge transport across the interface between 1T’-MoTe2
and a 3D metal contact seems to be less vulnerable to impurities and
defects.[15, 9, 13, 14, 11, 12, 10] The ultimate prospect of this design is to
combine the monolayer versions of these bulk phases, the H- and T’ phase, to
create a completely two-dimensional transistor.

In our previous work[17], we have used ab-initio calculations to demonstrate
how interface resonances at monolayer (ML) T’-H MoTe2 contacts create peaks in
the transmission such that the Schottky barrier heights extracted by the
thermionic emission model yield significantly different results than the
Schottky barrier seen in the density of states. In this work, we go a step
further and investigate the consequences of such quantum effects not only in
the Schottky contacts but also in T’-H-T’ MoTe2 transistors.

The investigations are conducted in two steps. The first step is to examine
the Schottky barriers and charge transport in two different T’-H MoTe2
contacts as the carrier density is gradually increased. FETs operate by
adjusting an electrostatic field using the gate electrode to draw in or push
out carriers in the semiconductor channel. The same effect can therefore be
created by gradually increasing or decreasing the carrier density in the
semiconductor at a single metal-semiconductor interface. Previous ab-initio
studies by Saha _et al._[18] and Urquiza _et al._[19] have investigated two
different doping densities of T’-H MoS2 Schottky contacts and concluded that
an increased doping density decreases both the barrier height and the
depletion width and our previous work[17] has demonstrated the same effect in
T’-H MoTe2 contacts. However, to our knowledge, this is the first time the
properties of these Schottky contacts have been investigated by gradually
changing the carrier density. This gradual change allows for mimicking the
effects of an applied gate potential in a field effect transistor and the
calculations therefore provide new fundamental understanding of the
electrostatic control over the charge transfer in these 2D contacts.

The second step is to directly study the charge transport in MoTe2 heterophase
FETs. The two contacts investigated in the previous part are combined to form
complete 2D FETs with channel lengths of 10, 18, and 20 nm. In case of the 20
nm channel, the self-consistent transport calculations includes 573 atoms
which takes this ab-initio approach to the limit of what is computationally
practical. Devices on this length scale have not previously been investigated
and help to shed light on the behavior of heterophase FETs in a regime where
the channel lengths are comparable to channel lengths used in commercial
devices.[1]

Our analysis of the Schottky contacts demonstrates that the local atomic
configuration at the interface completely dominates the Schottky barrier
height, depletion width, and charge transfer. Furthermore, we find that
quantum effects play an important part in the charge transfer at the
heterophase devices. Both in case of the contacts and the transistors, quantum
states bridge the Schottky barrier and contribute significantly to the total
current. In the contact, interface states increase the transmission
probability and lowers the effective Schottky barrier height. In the
transistors, standing waves in the channel increase this effect and results in
large current fluctuations in the transfer characteristics, even at a channel
length of 20 nm’s. However, we find that our devices are not competitive with
silicon technology neither in terms of efficiency nor output current. The
largest challenges for this heterophase design is that the Schottky barrier at
the source, which partly determines the efficiency, is controlled by the local
atomic configuration and carrier density in the channel and that the ballistic
ON-current is inhibited due to the restriction of momentum conservation along
the interface.

The paper is organized as follows. In Section I, we explain the details
regarding the applied methods and the setup of the Schottky contacts and
heterophase transistors. In Section II, we demonstrate how the Schottky
contacts respond to increasing the carrier density. This part serves as
background to understand the results of the complete transistor simulations
which are the subject of Section III. After demonstrating the results of the
device simulations, we discuss how these calculations compare to previous
experiments and calculations on the MoTe2 heterophase transistor and how the
performance of these devices could be optimized. This is the subject of
Section IV. Finally, we summarize the main conclusions in Section V.

##  I Method

Figure 1: The T’-H Schottky contact and T’-H-T’ transistor of monolayer MoTe2
on a substrate and back-gate. a shows the computational cell of the contacts
consisting of a central region and a T’- and H-phase electrode. The total cell
size is (30.0, 0.718, 15.0) nm. b shows the computational cell of the 2D
transistors with T’ phase source and drain and a H phase channel. The total
cell size is (10.0+$L_{C}$, 0.718, 15.0) nm with $L_{C}$ being 10.0, 18.0 or
20.0 nm. The middle part shows the atomic configurations of the two
investigated interface geometries seen from the top.

Interfaces between the two phases of MoTe2 are most commonly observed along
the zigzag direction or at a 60 degree angle to this direction.[9, 20, 21]
Here, we choose to investigate interfaces along the zigzag direction. In this
direction, the unit cells of the T’ and H phase TMD can be combined in 8
different ways while conserving the stoichiometry of the two phases. A
stability analysis of these eight interfaces is included in the Supplementary
Information.† The c- and d-type are chosen, since these represent some of the
most stable interfaces, the c-type in Te-rich conditions and the d-type in
intermediate and Mo-rich conditions, meanwhile demonstrating the large
variance in interface dipole represented by the different interface
geometries.

The c- and d-type interfaces are shown in the middle part of Figure 1. Both of
these interfaces have been found by combining the two unit cells along the
zigzag direction and thereby creating a strain of 2.12 % in the y-direction of
the T’ phase. The three unit cells of each phase closest to the interface have
been relaxed to a force tolerance of 0.02 eV/Å. The unit cell is doubled in
the y-direction compared to the conventional unit cells to allow for a
stabilizing distortion at the interface. This distortion is clearly visible at
the d-type interface.

We use density functional theory (DFT)[22, 23] with the Perdew-Burke-Ernzerhof
(PBE) exchange-correlation functional[24] and the non-equilibrium Green’s
function (NEGF) method as implemented in QuantumATK.[25] The wave functions
are expanded as linear combinations of atomic orbitals using PseudoDojo
pseudopotentials[26] using a cut-off energy of 100 Ha. The method is described
in details in our previous work[17] on similar interfaces where we also
justify using PBE to describe these systems. Within this description, the H
phase is a direct bandgap semiconductor with a band gap of 1.03 eV and the T’
phase is gapless. The computational cell of the Schottky contacts is shown on
Figure 1a. The total cell size of the central region is (30.0, 0.718, 15.0) nm
and the k-space grid of the NEGF calculations is (401, 6, 1). Figure 1b shows
the computational cell of the 2D planar FETs which combines a c- and d-type
interface as the drain and source contacts, respectively. The symmetry of the
H phase inhibits using the same interface geometry at both source and drain
which is why the device is a combination of the two different geometries. The
cell sizes of the central region is (10.0+$L_{C}$, 0.718, 15.0) nm where
$L_{C}$ is either 10.0, 18.0 or 20.0 nm and the k-space grids are equivalent
to the Schottky contact calculations.

To simulate the electrostatic environment in a FET, both the 2D Schottky
contacts and transistors have been placed on top of a thin metal back gate and
a dielectric region of 5 nm with a dielectric constant of 3.9 corresponding to
SiO2. Both are described as a classical continuum. The lowest atom of the
system is placed 2.4 Å above the dielectric region. This value is chosen since
it lies close to the experimental value of the distance between 2H MoS2 and
SiO2[27] and since we have found that a variation between 2.4 and 4 Å does not
change the barrier height or depletion width. The electrostatic effect of the
substrate is to screen the interface dipole which increases the charge
transfer since the potential difference between the metal and semiconductor is
conserved despite the presence of the substrate. This increases the depletion
width but conserves the barrier height.

We use a continuous doping model to add charge carriers to the semiconductor
and we apply a p-type doping since MoTe2 most commonly exhibit a p-type
intrinsic doping in ambient conditions.[9, 28] The electrons per atom are
explicitly modified by scaling the atomic densities.[29] We choose to dope all
the atoms in the system to avoid any artificial effects of defining a sharp
interface between the two phases. The extra holes in the metallic phase have a
negligible effect.

##  II Schottky contacts

Figure 2: The Schottky barrier evolution with varying levels of p-doping. a
and c show the valence band edge of the H phase for the c- and d-type
contacts, respectively. The grey line indicates the first Mo-atom belonging to
the H phase. b and d show the barrier heights as a function of the doping
level (circles) along with the maximum of the valence band (black line). The
density of states at the position of the grey line on a and c at each doping
level are included as a grey-scale contour to illustrate the density of
interface states between the Fermi level and valence band maximum.

Our starting point is to investigate the two Schottky contacts. The results
from these investigations will help to analyze the results of the simulations
on the complete transistor setup. The p-doping level is varied between
$10^{12}$ and $1.5\times 10^{13}$ cm-2 and the effects on the Schottky barrier
height, depletion width, and charge transport are investigated. The Schottky
barrier height is found from the projected density of states (DOS) as the
distance from the Fermi level to the minima of the valence band. The valence
band edge is found by defining the gap region as energies at which the
projected DOS is negligible. The evolution of the valence band positions of
the two contacts, when the doping is increased, are shown on Figure 2a and 2c.
The Fermi level is placed at the energy zero point and the x-coordinate of the
leftmost Mo-atom, which belonged to the H-phase before relaxation, is marked
by a vertical grey line. The c-type contact only shows band bending for the
three lowest doping levels. From then on, the valence bands are flat and the
barrier height coincides with the maximum of the valence band. The d-type
contact show band bending for all the investigated doping levels.

The barrier height is shown together with the valence band maxima as a
function of the doping level on Figure 2b and 2d. The size of the barrier
differs significantly between the two interface geometries and so does the
behavior with the doping density. The barrier height does not vary smoothly
with the doping level but rather jumps between energy plateaus. Our previous
work[17] and the work by Liu _et al._[30] have demonstrated that these
heterophase MoTe2 contacts host interface states which increases the tunneling
probability. At the interfaces investigated here, some of these states bridge
the Schottky barrier and result in both increased tunneling probability and an
effectively lower barrier than what the strength of the interface dipole
dictates. This effect is illustrated on Figure 3. The interface dipole is
created due to the depletion of holes in the valence band and the
corresponding charge accumulation in the metal which screens the depletion
charge. The band bending therefore reflects the strength of the interface
dipole and is an electrostatic effect. The density of states, on the other
hand, is a representation of the available quantum states which electrons and
holes can occupy. The barrier height extracted from the density of states
therefore represents an effective barrier for quantum transport which does not
necessarily coincide with the electrostatic barrier. This effective barrier
will in the following be referred to as the DOS barrier, $\Phi^{DOS}$, in
agreement with the definition in our previous work.[17]

Figure 3: The effect of interface states on the DOS barrier, $\Phi^{DOS}$,
between H- and T’ phase MoTe2. The electrostatic field from the interface
dipole causes the valence band to bend downwards creating a barrier height
represented by the dashed black line. The DOS associated with interface states
bridges this electrostatically defined barrier reducing the effective barrier
for quantum transport illustrated by the green line.

The density of interface states have been included on Figure 2b and 2d as a
contour plot. The contour represent the DOS at the position of the vertical
grey line on Figure 2a and 2c. The states of interest lie at energies between
the bottom of the barrier and the maximum of the valence band. Energies below
this range host the valence band states and energies above cannot contribute
to transport since these lie within the band gap of the H phase electrode. The
effect of the interface DOS is most easily seen in the case of the d-type
contact. The barrier bottom is seen to jump between windows of low interface
DOS around -0.17 and -0.26 eV. Note, that the slope of the interface states on
this plot are much lower than the slope of the valence band maxima. This
demonstrates that the interface states originate from the T’-phase rather than
the H-phase which was the case for the interface resonances in our previous
study as well. The interface DOS illustrated on the contour plot can be
associated both with regular metal induced gap states, i.e. the exponential
tail of bulk states in the metal, and such interface resonances.

Figure 4: The electron density at the two Schottky contacts (a and b) and the
depletion width (c and d) of the d-type interface with varying p-doping
following the color-scale on Figure 2. The depletion width is compared to the
predictions by Ilatikhameneh _et al._[31] (dashed line) and Nipane _et
al._[32] (dotted line). d shows the depletion width as a function of the work-
function difference divided by the doping density and a linear fit to the
calculated depletion widths (solid line).

These results illustrate that the barrier height depends non-trivially on the
doping level due to the presence of interface states. We will now investigate
how the interface dipole and depletion width vary with the doping level. In
these monolayer contacts, the electronic density and the potential depend a
lot on the atomic rearrangements at the phase boundary. To visualize how the
interface dipole changes with the doping level, the electronic density of a
contact with a low doping level, $N_{2D}^{ref}=1\times 10^{11}$ cm-2, is used
as a reference. Almost no charge transfer is possible in the reference system
since less than 0.02 holes are available in the entire semiconductor region.
After the subtraction of the reference density, the doping difference between
the reference device and this device is subtracted in order to remove the
background doping density,

$\displaystyle\Delta n_{e}(x)$
$\displaystyle=\frac{1}{W}\int[n_{e}(x,y,z)-n_{e}^{ref}(x,y,z)]\,\mathrm{d}y\,\mathrm{d}z$
$\displaystyle-(N_{2D}^{dop}-N_{2D}^{ref}).$ (1)

This difference in electronic density can be seen on Figure 2a and 2b for the
c-type and d-type contacts, respectively. The difference in the size of the
dipole is quite clear. The c-type contact shows much less variation of the
electronic density through the device which is a footprint of a very small
charge accumulation in the depletion region and corespondingly a small barrier
height. The d-type contact shows a much larger variation and the interface
dipole is seen to become more narrow in space as the doping is increased. The
minimum of the density can be used to illustrate the smooth variation of the
electronic density with the doping level and is shown on the insets.

Since the depletion width vanish already at a doping level of $4\times
10^{12}$ cm-2 in the c-type contact, we only consider how the depletion width
varies with the doping level of the d-type contact. The depletion width can be
found from the bending of the valence bands on Figure 2c. The values are
compared on Figure 4c and 4d to the expressions derived using classical
electrostatics by Ilatikhameneh _et al._[31] and Nipane _et al._[32],

$\displaystyle x_{D}^{\mathrm{Ilatikhameneh}}$
$\displaystyle=\frac{\pi\varepsilon_{\mathrm{eff}}\Delta V}{\ln(4)qN_{2D}},$
$\displaystyle x_{D}^{\mathrm{Nipane}}$
$\displaystyle=\frac{\pi^{2}\varepsilon_{\mathrm{eff}}\Delta V}{8GqN_{2D}}.$

$\Delta V$ is the work function difference between the two materials in the
junction, $\varepsilon_{\mathrm{eff}}$ is the effective dielectric constant
which is determined by the dielectrics surrounding the 2D material which in
this case is
$\varepsilon_{\mathrm{eff}}=(\varepsilon_{above}+\varepsilon_{below})/2=2.45\varepsilon_{0}$,
$G\approx 0.915$ is Catalan’s constant, and $q$ is the elementary charge.
These depletion widths have been calculated using the work function of the H-
and T’ phase electrode at each doping level with computational parameters
matching those used for the electrodes of the devices.

Comparing to the ab-initio results, it is clear that the depletion width can
not be predicted from the difference in workfunctions and doping level alone.
First of all, the values of the two different phase boundaries, as already
mentioned, differ a lot and secondly, the depletion width of the d-type
contact is considerably shorter than predicted by the classical models.
However, as demonstrated by Figure 4d, the classically predicted linear
scaling with $\Delta V/N_{2D}$ seems to be fulfilled.

Figure 5: Total and spectral current though the Schottky contacts at varying
levels of p-doping. a shows the total current through the two contacts as a
function of the doping level and at a bias of $V_{bias}=10$ mV. b shows the
spectral current of the d-type contact at varying doping levels. The color-
scale follows Figure 2.

Finally, we investigate how the charge transport though the contact changes
with the carrier density. We apply a small positive bias of $V_{bias}=10$ mV
across the contacts. This drives a hole current from the metal to the
semiconductor, the current can be calculated using the Landauer-Büttiker
formula,

$\displaystyle I$ $\displaystyle=\frac{2q}{h}\int T(E,\mu_{L},\mu_{R})\times$
$\displaystyle\left[f\left(\frac{E-\mu_{L}}{k_{B}T}\right)-f\left(\frac{E-\mu_{R}}{k_{B}T}\right)\right]\mathrm{d}E.$
(2)

$h$ is Plack’s constant, $k_{B}$ is the Boltzmann constant,
$T(E,\mu_{L},\mu_{R})$ is the transmission, $f$ is the Fermi-Dirac
distribution, $\mu_{L}$ and $\mu_{R}$ are the chemical potentials of the two
electrodes, defined as $\mu_{L}-\mu_{R}=qV_{bias}$, and $T$ is the
temperature. We use $T=300$ K to find the current at room temperature. The
resulting currents through the two contacts are shown on Figure 5a. The
current has a relatively smooth dependence on the carrier density and the
magnitude differs about two orders of magnitude between the two geometries at
low doping levels but reaches a comparable value at high doping levels. Note,
that the current is shown on a logarithmic scale. The variation of the
spectral current through the d-type contact is shown on Figure 5b. The
spectral current is the term inside the integral in eq (2). Two peaks are
present at the energies where a high and relatively constant interface DOS can
be identified on Figure 2d at all the investigated carrier densities. Those
states therefore contribute significantly to the current regardless of the
value of the Schottky barrier height and we believe that they can be
associated with interface resonances.

To summarize, the Schottky barrier in heterophase MoTe2 Schottky contacts can
be electrostatically tuned to increase the current such that these can be
utilized in a transistor. As expected, the interface dipole and charge
transfer respond in a smooth manner when increasing the carrier density. This
smooth increase is also reflected in the depletion width which scales
reciprocally with the carrier density. However, the size of the barrier
height, the length of the depletion width, and the total magnitude of the
current are determined by atomic-scale and quantum effects. It has been
established several times[33, 18, 19], which we also confirm with these
calculations, that the barrier height in the TMD heterojunctions depends on
the local atomic positions at the phase boundary. From our investigation of
the depletion width, we therefore naturally find that the magnitude of this
quantity, as well as the size of the interface dipole, is highly dependent on
the interface geometry as well. The barrier height behaves non-trivially when
the carrier density is increased since interface states can bridge the
barrier. Some of these interface states increase the transmission probability
significantly and therefore also increase the current trough the contact.

##  III Heterophase FETs

After concluding these initial investigations of the isolated Schottky
contacts, we are ready to simulate a complete transistor setup with source,
drain, channel and back-gate as shown on Figure 1b. We combine the two
Schottky contacts to create transistors with channel lengths of 10, 18 and 20
nm’s and choose a doping level of $N_{A}=5\times 10^{12}$ cm-2.

Figure 6: Transfer characteristics of the heterophase FETs at $V_{bias}=-10$ mV. The ON-currents and sub-threshold slopes (SS) are indicated. The SS values are found from the fit indicated by the red lines and the sub-threshold regimes are indicated by the shaded red regions. The flat-band condition, $V_{FB}$, is marked by the green lines and the threshold potential at $V_{TH}=-2.7$ V by the blue lines. Purple stars indicate rises in the current due to resonant tunneling. Ref. | $t_{EOT}$ | $t_{ch}$ | $L_{C}$ | $I^{ON}$ | $V_{sd}(I^{ON})$ | $V_{DD}$ | SS |
---|---|---|---|---|---|---|---|---
Zhang  _et al._ [12] | 1.87 nm | 8 nm |  4 $\mu$m |  0.25 $\mu$A/$\mu$m | 0.1 V | 0.8 V | 69 mV/dec |
Sung  _et al._ [9] | 300 nm | 5-6 nm |  4-15 $\mu$m |  15 nA/$\mu$m | 50 mV | 40 V |  |
Ma  _et al._ [10] | 300 nm | 8 nm |  20 $\mu$m |  39.6 nA/$\mu$m | -0.1 V | 100 V |  |
This Work | 5 nm | ML | 10 nm |  1.6 nA/$\mu$m | -10 mV | 5.7 V | 0.23 V/dec |
This Work | 5 nm | ML | 18 nm |  0.79 nA/$\mu$m | -10 mV | 5.1 V | 0.14 V/dec |
This Work | 5 nm | ML | 20 nm |  1.2 nA/$\mu$m | -10 mV | 5.1 V | 0.12 V/dec |
This Work* | 5 nm | ML | 20 nm |  21 nA/$\mu$m | -0.1 V | - | - |
This Work* | 5 nm | ML | 20 nm |  0.15 mA/$\mu$m | -0.7 V | - | - |
IDRS 2022[1] | 1 nm |  6 nm111Fin width of the finFET design. | 16 nm |  912 $\mu$A/$\mu$m | 0.7 V | 0.7 V | 82 mV/dec |
Table 1: Equivalent oxide thicknesses, $t_{EOT}$, channel thicknesses,
$t_{ch}$, channel lengths, $L_{C}$, ON-currents, $I^{ON}$, power supply
voltages, $V_{DD}$, and sub-threshold slopes (SS) of MoTe2 based heterophase
devices compared to the IRDS 2022 goals.[1] *Extrapolated non self-
consistently.

The transfer characteristics of the 2D FETs are shown on Figure 6 at a bias of
$V_{bias}=-10$ mV. This bias results in a device with a d-type barrier at the
source and flat bands at the drain. The transfer characteristics are
calculated by gradually varying the potential between the back-gate and source
electrode and calculate the current using eq. (2). At a gate potential of 3.0
V, in case of the 10 nm channel, and 2.4 V, in case of the other two devices,
the current is switched completely off. Above those potentials, the devices
shift polarity and begin to be dominated by electron transport in the
conduction band. The sub-threshold regimes where the current rises
exponentially are indicated as the shaded red regions on Figure 6. The sub-
threshold slopes in this region are 0.24, 0.14 and 0.12 V/dec for the 10, 18
and 20 nm FETs respectively. The values of the 18 and 20 nm devices agree
reasonably well with the SS value of a heterophase MoS2 device fabricated by
Nourbakhsh _et al._[34] of 120 mV/dec whereas it is considerably larger than
the SS value measured in a MoTe2 heterophase device fabricated by Zhang _et
al._[12] of merely 69 mV/dec. The impressively low SS value in the latter
device is most likely due to the optimized setup used in this case. The
heterophase transistor was placed on the high-$\kappa$ dielectric HfO2 with an
equivalent oxide thickness (EOT) of 1.87 nm’s and the back gate positioned
only below the 2H phase. A low EOT generally both decreases the SS value and
increases the ON-current in a device.[35]

In the Supplementary Information,† we use a classical electrostatic model to
estimate how an increment in the applied gate potential will result in a
potential increment at the position of the semiconductor in a device setup
such as ours. The relation between the two is given by the ideality factor,
$\eta$, which determines the SS value, $\mathrm{SS}=\eta\times 60$ mV/dec. The
model assumes a constant charge density in the channel and estimates the
effect of the geometry of the device which creates a capacitive potential
divider formed by the gate capacitance in series with the capacitance of the
substrate and vacuum region above the transistor. This model predicts SS
values of 0.20, 0.12, and 0.12 V/dec for the 10, 18 and 20 nm channels
respectively. Most of the deviation from the thermodynamic limit of
$\mathrm{SS}=60$ mV/dec can therefore be explained by the setup of the device.
The remaining discrepancy is most likely due to the variation in the charge
density along the channel and the tunneling between the source and drain. The
source-drain tunneling is most important in the 10 nm device which explains
why the SS value of this device is very large.

The flat-band condition, $V_{FB}$, is met when the exponential dependence of
the current breaks off. This point is marked by the vertical green line on
Figure 6. At lower gate potentials, the devices enter the Schottky barrier
regime where the current is conducted as a mixture of thermionic emission and
tunneling of the holes. The threshold potential at $V_{G}=V_{TH}=-2.7$ V is
defined as the point where the valence band edge rises above the chemical
potential of the source and has been identified by calculating the projected
DOS. In this condition, we find ON-currents of 1.6, 0.79 and 1.2 nA/$\mu$m,
respectively. Previous experimental measurements on MoTe2 2H-1T’ heterophase
devices with channel thickness between 5 and 8 nm and channel lengths on the
$\mu$m scale have presented ON-currents around 15 nA/$\mu$m at a bias of 50 mV
and 39.6 nA/$\mu$m at a bias of 0.1 V.[9, 10] By non self-consistently
increasing the bias, we can compare our ON-currents with these values and get
a reasonable agreement. The extrapolated ON-currents of the 20 nm device are
summarized in Table 1 and the extrapolated ON-currents of the two other
devices can be found in the Supplementary Information.† In the table, we also
include the measurements on the optimized device fabricated by Zhang _et
al._[12] which results in an ON-current which is two orders of magnitude
higher than the other measurements and our results.

The energy efficiency of the devices is reflected both in the value of the
sub-threshold slope and in the power supply voltage, $V_{DD}$. We define
$V_{DD}$ as the gate potential difference between the threshold potential and
the potential of minimal current. This value is related to the Schottky
barrier height since a lower barrier at the source will result in a lower
flat-band potential and decrease this value. The 10 nm device has a $V_{DD}$
value of 5.7 V and the 18 and 20 nm devices have $V_{DD}$ values of 5.1 V.
These values are around an order of magnitude higher than what typical silicon
based devices require.[1] The power supply voltage required in the 10 nm
device is larger than in the two other devices due to the proximity of the
contacts. The interface dipole fields at the source and drain interacts which
results in a less efficient gate response and a larger flat-band potential.
The effect is illustrated in the Supplementary Information.† In the 10 nm
device, also the OFF-current is larger than in the other devices since this
device isn’t fully turned off when the switch in polarity occurs. It is
important to note that the switch in polarity is directly related to the
bandgap of the channel. PBE systematically underestimates the bandgap of the
2D TMDs which means that the potential of the polarity switch is
underestimated as well.

Figure 7: Spectral current in units of q/h in the transistor with a channel
length of 20 nm. The two large current peaks marked by the purple stars on
Figure 6 are denoted QR1 and QR2. At the quantum resonances (QR), significant
contributions in the tunneling current occur through these peaks.

The Schottky barrier regime between $V_{TH}$ and $V_{FB}$ is dominated by
large spikes in the current which have been marked by the purple stars on
Figure 6. The spikes in the current occur due to resonant tunneling and the
effect is especially pronounced in the 20 nm device. On Figure 7, the spectral
current of this device is shown at different gate potentials between the first
purple star at $V_{G}=V_{QR1}=1.0$ V and the threshold potential. As the
device is turned on, the peaks associated with the quantized states in the
channel appear and disappear whenever a quantized state contribute to the
tunneling current. At the two potentials corresponding to the purples stars,
large peaks due to resonant tunneling can be identified at low energies.

Figure 8: Quantum effects in MoTe2 heterophase transistors. a shows the
projected DOS of the 20 nm transistor at a gate potential of $V_{G}=0.3$ V.
The chemical potential of the drain, $\mu_{D}$, and the source $\mu_{S}$ are
indicated by the dotted black lines. The interface states are marked by green
contours and some of the standing waves in blue contours. b shows the
transmission spectrum where the peaks due to the standing waves and interface
states have been marked in blue and green. c shows the k-resolved transmission
perpendicular to the transport direction. The most predominant peak in the
transmission spectrum is marked by the blue arrows.

The effect is most pronounced at a gate potential of 0.3 V. On Figure 8a and
8b, the projected DOS at this gate potential is shown along with the
transmission spectrum. Two kinds of quantum phenomena cause the peaks to occur
in the transmission and spectral current. The first type is the interface
states discussed in the previous chapter which have been highlighted in green
on the figure, the second type is standing waves in the channel which form due
to quantum confinement. These have been marked in blue. In these devices,
extraordinary large peaks occur in the transmission when a standing wave
bridge the Schottky barrier. An example of this is highlighted by the blue
arrows on the figure.

Note, that the transmission is quite low even below the bottom of the barrier.
This is due to the restriction of momentum conservation in ballistic
transport. In p-type MoTe2 heterophase contacts, the ballistic transport is
significantly reduced since only very few k-points along the
$Y\rightarrow\Gamma$ path are occupied both in the metallic and semiconducting
phase. Conserving the momentum perpendicular to the transport direction is
therefore only possible for a narrow range of k-points. This is illustrated in
Figure 8c, where the $k_{y}$-dependence of the transmission is illustrated. We
therefore expect that including the effects of electron phonon scattering
might actually increase the device performance.

To summarize, we find that resonant tunneling and momentum conservation play a
very large role for the device performance even in these relatively long
channel devices. Comparing with the behavior of the Schottky contacts, we find
that the standing waves in the channel affect the current to a larger degree
than the effect of the interface resonances. Since the standing waves are a
property of the semiconductor and not the metal, the electrostatic field from
the gate moves these significantly up and down in energy. When a standing wave
bridges the Schottky barrier, the current increases dramatically. This results
in large current fluctuations in the transfer characteristics with peaks in
the current rising more than an order of magnitude. As it was demonstrates on
Figure 8a and 8b, the interface states also contribute to the current in the
transistors but the effect is smaller and relatively constant with the gate
potential since they have origin in the metallic phase.

##  IV Discussion

Previous calculations on similar devices have been performed by Fan _et
al._[36] These devices differ from ours since the T phase rather than T’ phase
was used as contacts, the channel lengths were much smaller, between 5 and 8
nm, the interface was along the armchair direction, and the devices had both a
top and back gate with EOT’s of only 0.5 nm. Furthermore, that study used the
Local Density Approximation (LDA) functional rather than the PBE functional.
The results from those calculations showed ON-currents which are around an
order of magnitude larger than those reported here and a power supply voltage
around 1 V. For WTe2 and WSe2 based devices, the results from Fan _et al._
suggest that devices based on the T phase generally conduct more current than
T’-based devices. The reason for this difference might lie in the difference
in dispersion between the T and T’ phase. Since the T phase is metallic
whereas the T’ phase is semi-metallic, the states are generally more densely
distributed in k-space which allows for a greater probability of ballistic
transport of carriers. Whether the rest of the difference in the ON-current
results from the difference in the device setup or in the method between this
study and ours is difficult to say. However, the behavior of the interfaces
along the zigzag direction seems most relevant since the experimental evidence
suggests that the interfaces occur along this direction.[9, 20, 21]

In Table 1, we have summarized our results and compared to the three
experimental studies mentioned in the previous section and to the
International Roadmap for Devices and Systems (IRDS) requirements[1] of a 2022
high-performance logic device. In summary, our 18 and 20 nm devices show a
performance which is comparable to previous measurements on devices on top of
SiO2 substrates. Comparing to the IRDS goal, the energy efficiency and total
ON-current of the devices are not competitive with the performance of silicon
based devices. The experimental results by Zhang _et al._[12], on the other
hand, suggests that the IRDS requirements might be fulfilled by using a
high-$\kappa$ dielectric as the gate oxide which is also the common practice
in silicon based devices. However, the device dimensions in the study by Zhang
_et al._ is much larger than both our calculations and the IRDS requirements.
This therefore opens the question of whether the heterophase design will be
competitive when scaling down the channel length of the devices. Judging from
these considerations, the idealised monolayer limit with perfect interfaces
might not be preferable since ballistic transport through such geometries are
limited by momentum conservation. However, electron-phonon scattering might
remove some of this effect and since the device performance is dominated by
the local atomic configuration at the interface, the presence of defects in
the interface might also dramatically alter the behavior. Furthermore, as
demonstrated for the Schottky contacts, a high carrier density will smear out
the difference between interface geometries and result in a relatively low
Schottky barrier. A heavy doping of the H phase might therefore be a viable
approach to obtain a competitive design.

##  V Conclusion

Our simulations of the heterophase FET devices demonstrated that quantum
effects play a large role in the performance of these devices even at a
channel length of 20 nm’s. We find that quantum resonances result in resonant
tunneling and large fluctuations in the current in the Schottky barrier regime
of the transfer characteristics. The quantum resonances stem from both
interface states at the contacts and standing waves in the channel. The
interface states are a material specific property of the devices whereas the
standing waves are an effect of the scale and dimensionality of the device and
can be expected to be a general property of 2D transistors. Furthermore, we
find that the restriction of momentum conservation inhibit the ballistic
transport significantly in these idealised devices.

All-in-all, our results demonstrate many challenges related to the heterophase
FET design. Effects from electron-phonon scattering and interface defects
might actually improve the performance and calculations estimating such
effects could shed light on this possibility. Likewise, moving from the
monolayer to the few layer limit or ensuring a high carrier density in the
channel might resolve some of the difficulties since this will somewhat smear
out some of the dependency on the interface geometry. The results by
Nourbakhsh _et al._[34] where MoS2 based heterophase devices with both one and
three layers were investigated showed that the 3-layer device performs better,
which agrees with this conclusion.

###### Acknowledgements.

This work is partly funded by the Innovation Fund Denmark (IFD) under File No.
5189-00082B. We would like to thank Kurt Stokbro and Daniele Stradi for their
contributions.

## References

  * IRD [2020] IEEE International Roadmap for Devices and Systems (2020).
  * Akinwande  _et al._ [2019a] D. Akinwande, C. Huyghebaert, C.-H. Wang, M. I. Serna, S. Goossens, L.-J. Li, H.-S. P. Wong, and F. H. L. Koppens, Graphene and two-dimensional materials for silicon technology, Nature 573, 507 (2019a).
  * Akinwande  _et al._ [2019b] D. Akinwande, C. Huyghebaert, C.-H. Wang, M. I. Serna, S. Goossens, L.-J. Li, H.-S. P. Wong, and F. H. L. Koppens, Graphene and two-dimensional materials for silicon technology, Nature 573, 507 (2019b).
  * Balasubramaniam  _et al._ [2019] B. Balasubramaniam, N. Singh, P. Kar, A. Tyagi, J. Prakash, and R. K. Gupta, Engineering of transition metal dichalcogenide based 2D nanomaterials through doping for environmental applications, Molecular Systems Design and Engineering 4, 804 (2019).
  * Allain  _et al._ [2015] A. Allain, J. Kang, K. Banerjee, and A. Kis, Electrical contacts to two-dimensional semiconductors, Nature Materials 14, 1195 (2015).
  * Kim  _et al._ [2017] C. Kim, I. Moon, D. Lee, M. S. Choi, F. Ahmed, S. Nam, Y. Cho, H. J. Shin, S. Park, and W. J. Yoo, Fermi Level Pinning at Electrical Metal Contacts of Monolayer Molybdenum Dichalcogenides, ACS Nano 11, 1588 (2017).
  * Duerloo  _et al._ [2014] K.-A. N. Duerloo, Y. Li, and E. J. Reed, Structural phase transitions in two-dimensional Mo- and W-dichalcogenide monolayers, Nature Communications 5, 4214 (2014).
  * Yoo  _et al._ [2017] Y. Yoo, Z. P. DeGregorio, Y. Su, S. J. Koester, and J. E. Johns, In-Plane 2H-1T’ MoTe2 Homojunctions Synthesized by Flux-Controlled Phase Engineering, Advanced Materials 29, 1 (2017).
  * Sung  _et al._ [2017] J. H. Sung, H. Heo, S. Si, Y. H. Kim, H. R. Noh, K. Song, J. Kim, C.-S. Lee, S.-Y. Seo, D.-H. Kim, H. K. Kim, H. W. Yeom, T.-H. Kim, S.-Y. Choi, J. S. Kim, and M.-H. Jo, Coplanar semiconductor–metal circuitry defined on few-layer MoTe2 via polymorphic heteroepitaxy, Nature Nanotechnology 12, 1064 (2017).
  * Ma  _et al._ [2019] R. Ma, H. Zhang, Y. Yoo, Z. P. Degregorio, L. Jin, P. Golani, J. Ghasemi Azadani, T. Low, J. E. Johns, L. A. Bendersky, A. V. Davydov, and S. J. Koester, MoTe2 Lateral Homojunction Field-Effect Transistors Fabricated using Flux-Controlled Phase Engineering, ACS Nano 13, 8035 (2019).
  * Zhang  _et al._ [2019a] X. Zhang, Z. Jin, L. Wang, J. A. Hachtel, E. Villarreal, Z. Wang, T. Ha, Y. Nakanishi, C. S. Tiwary, J. Lai, L. Dong, J. Yang, R. Vajtai, E. Ringe, J. C. Idrobo, B. I. Yakobson, J. Lou, V. Gambin, R. Koltun, and P. M. Ajayan, Low Contact Barrier in 2H/1T’ MoTe2 In-Plane Heterostructure Synthesized by Chemical Vapor Deposition, ACS Applied Materials and Interfaces 11, 12777 (2019a).
  * Zhang  _et al._ [2019b] Q. Zhang, X.-F. Wang, S.-H. Shen, Q. Lu, X. Liu, H. Li, J. Zheng, C.-P. Yu, X. Zhong, L. Gu, T.-L. Ren, and L. Jiao, Simultaneous synthesis and integration of two-dimensional electronic components, Nature Electronics 2, 164 (2019b).
  * Xu  _et al._ [2019a] X. Xu, S. Chen, S. Liu, X. Cheng, W. Xu, P. Li, Y. Wan, S. Yang, W. Gong, K. Yuan, P. Gao, Y. Ye, and L. Dai, Millimeter-Scale Single-Crystalline Semiconducting MoTe2 via Solid-to-Solid Phase Transformation, Journal of the American Chemical Society 141, 2128 (2019a).
  * Xu  _et al._ [2019b] X. Xu, S. Liu, B. Han, Y. Han, K. Yuan, W. Xu, X. Yao, P. Li, S. Yang, W. Gong, D. A. Muller, P. Gao, Y. Ye, and L. Dai, Scaling-up Atomically Thin Coplanar Semiconductor–Metal Circuitry via Phase Engineered Chemical Assembly, Nano Letters 19, 6845 (2019b).
  * Cho  _et al._ [2015] S. Cho, S. Kim, J. H. Kim, J. Zhao, J. Seok, D. H. Keum, J. Baik, D.-H. Choe, K. J. Chang, K. Suenaga, S. W. Kim, Y. H. Lee, and H. Yang, Phase patterning for ohmic homojunction contact in MoTe2, Science 349, 625 (2015).
  * Song  _et al._ [2016] S. Song, D. H. Keum, S. Cho, D. Perello, Y. Kim, and Y. H. Lee, Room Temperature Semiconductor-Metal Transition of MoTe2 Thin Films Engineered by Strain, Nano Letters 16, 188 (2016).
  * Jelver  _et al._ [2021] L. Jelver, D. Stradi, K. Stokbro, and K. W. Jacobsen, Schottky barrier lowering due to interface states in 2D heterophase devices, Nanoscale Advances 3, 567 (2021), arXiv:1911.09521 .
  * Saha and Mahapatra [2017] D. Saha and S. Mahapatra, Asymmetric Junctions in Metallic-Semiconducting-Metallic Heterophase MoS2, IEEE Transactions on Electron Devices 64, 2457 (2017).
  * Urquiza and Cartoixà [2020] M. L. Urquiza and X. Cartoixà, Schottky barriers, emission regimes and contact resistances in 2H-1T’ MoS2 lateral metal-semiconductor junctions from first-principles, 2D Materials 7, 045030 (2020).
  * Xu  _et al._ [2020] X. Xu, B. Han, S. Liu, S. Yang, X. Jia, W. Xu, P. Gao, Y. Ye, and L. Dai, Atomic‐Precision Repair of a Few‐Layer 2H‐MoTe2 Thin Film by Phase Transition and Recrystallization Induced by a Heterophase Interface, Advanced Materials 32, 2000236 (2020).
  * Han  _et al._ [2020] B. Han, C. Yang, X. Xu, Y. Li, R. Shi, K. Liu, H. Wang, Y. Ye, J. Lu, D. Yu, and P. Gao, Correlating the electronic structures of metallic/semiconducting MoTe2 interface to its atomic structures, National Science Review , 1 (2020), arXiv:1911.09885 .
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous Electron Gas, Phy. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects, Phys. Rev. 140, A1133 (1965).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Smidstrup  _et al._ [2019] S. Smidstrup, T. Markussen, P. Vancraeyveld, J. Wellendorff, J. Schneider, T. Gunst, B. Verstichel, D. Stradi, P. A. Khomyakov, U. G. Vej-Hansen, M.-E. Lee, S. T. Chill, F. Rasmussen, G. Penazzi, F. Corsetti, A. Ojanperä, K. Jensen, M. L. N. Palsgaard, U. Martinez, A. Blom, M. Brandbyge, and K. Stokbro, QuantumATK: an integrated platform of electronic and atomic-scale modelling tools, Journal of Physics: Condensed Matter 32, 015901 (2019).
  * van Setten  _et al._ [2018] M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table, Computer Physics Communications 226, 39 (2018).
  * Woollam and Somoano [1976] J. A. Woollam and R. B. Somoano, Superconducting critical fields of alkali and alkaline-earth intercalates of mo${\mathrm{s}}_{2}$, Phys. Rev. B 13, 3843 (1976).
  * Liu  _et al._ [2018a] J. Liu, Y. Wang, X. Xiao, K. Zhang, N. Guo, Y. Jia, S. Zhou, Y. Wu, Q. Li, and L. Xiao, Conversion of Multi-layered MoTe2 Transistor Between P-Type and N-Type and Their Use in Inverter, Nanoscale Research Letters 13, 291 (2018a).
  * dop [2021] The details on the doping methods implemented in the atk code can be found in the quantumatk technical notes (2021).
  * Liu  _et al._ [2018b] S. Liu, J. Li, B. Shi, X. Zhang, Y. Pan, M. Ye, R. Quhe, Y. Wang, H. Zhang, J. Yan, L. Xu, Y. Guo, F. Pan, and J. Lu, Gate-tunable interfacial properties of in-plane ML MX2 1T’-2H heterojunctions, Journal of Materials Chemistry C 6, 5651 (2018b).
  * Ilatikhameneh  _et al._ [2017] H. Ilatikhameneh, T. Ameen, F. Chen, H. Sahasrabudhe, G. Klimeck, and R. Rahman, Dramatic Impact of Dimensionality on the Electrostatics of PN Junctions (2017).
  * Nipane  _et al._ [2017] A. Nipane, S. Jayanti, A. Borah, and J. T. Teherani, Electrostatics of lateral p-n junctions in atomically thin materials, Journal of Applied Physics 122, 194501 (2017).
  * Li  _et al._ [2017] P. Li, Y. Wen, X. He, Q. Zhang, C. Xia, Z.-M. Yu, S. A. Yang, Z. Zhu, H. N. Alshareef, and X.-X. Zhang, Evidence for topological type-II Weyl semimetal eWTe2, Nature Communications 8, 2150 (2017).
  * Nourbakhsh  _et al._ [2016] A. Nourbakhsh, A. Zubair, R. N. Sajjad, T. K. Amir, W. Chen, S. Fang, X. Ling, J. Kong, M. S. Dresselhaus, E. Kaxiras, K. K. Berggren, D. Antoniadis, and T. Palacios, MoS2 Field-Effect Transistor with Sub-10 nm Channel Length, Nano Letters 16, 7798 (2016).
  * Houssa  _et al._ [2016] M. Houssa, D. Athanasios, and A. Molle,  _2D materials for nanoelectronics_ (CRC Press, 2016).
  * Fan  _et al._ [2018] Z. Q. Fan, X. W. Jiang, J. Chen, and J. W. Luo, Improving Performances of In-Plane Transition-Metal Dichalcogenide Schottky Barrier Field-Effect Transistors, ACS Applied Materials and Interfaces 10, 19271 (2018).