# Understanding resistance oscillation in CsV3Sb5 superconductor

Jung Hoon Han  Department of Physics, Sungkyunkwan University, Suwon 16419,
South Korea  Patrick A. Lee  Department of Physics, Massachusetts Institute of
Technology, Cambridge, Massachusetts 02139, USA

###### Abstract

A recent demonstration of the periodic oscillation of resistance in the thin
film of CsV3Sb5 superconductor with a hole in the film suggests that
charge-$4e$ and charge-$6e$ Cooper pairs may have condensed in this compound.
While exciting, such interpretation calls for a precise determination of the
effective area for the passage of Cooper pairs from one end of the lead to the
other. Unlike the traditional Little-Parks effect where the rim around the
hole is thin, the effective hole area is not obviously defined for the “thick-
rim geometry” adopted in Ref. Ge et al., 2022. Here, we note that the
experiment was conducted in a regime where the superconctivity is strongly
fluctuating, which motivates an analysis based on the spacetime formulation of
the time-dependent Ginzburg-Landau theory. We argue that under appropriate
conditions, the optimal semi-classical path is not the geometrically shortest
one, but the one that moves along the edge of the hole and takes advantage of
the reduced fluctuations at the boundary of the hole. The condition for the
cross-over from the geometrically shortest path to the path that sticks to the
wall is clarified. In such a scenario, the geometric area of the hole indeed
emerges as the effective area for the flux, providing a theoretical
justification to the interpretation given in Ref. Ge et al., 2022. The
conclusion of our analysis may have implication for similar experiments in
other superconductors where the geometry of the device is not obviously that
of the Little-Parks experiment employing the thin wall, but of the thick-rim
type such as used in Ref. Ge et al., 2022.

##  I Introduction

An explosion of activity has taken place around the discovery of vanadium-
based kagome metals and the superconducting phase in one of the compounds,
CsV3Sb5 (CVS) Ortiz et al. (2019, 2020). The multiple Fermi pocket structure
in the normal state of CVS compound raises the possibility of pair density
wave (PDW) - superconducting orders at a finite momentum ${\bf Q}$ vector, as
recently detected experimentally Chen et al. (2021); Zhou and Wang (2021).
Furthermore, for symmetry reasons, PDW’s at several ${\bf Q}$ vectors may
appear simultaneously. As analylzed some time ago Agterberg and Tsunetsugu
(2008); Berg et al. (2009); Radzihovsky and Vishwanath (2009); Agterberg et
al. (2011), it is possible that above some temperature the PDW loses long
range order, but a composite order $\Delta$ made up of the product of several
primary PDW at different ${\bf Q}$ vectors according to

$\displaystyle\Delta=\prod_{i=1}^{N}\Delta_{{\bf Q}_{i}}$ (1)

may survive as the dominant order. Each ${\bf Q}_{i}$ here is the finite-
momentum pairing field in the PDW phase and the order parameter $\Delta$
carries momentum $\sum_{i}{\bf Q}_{i}$. The effective charge associated with
$\Delta$ is $N(2e)$. In the case of kagome metal, the likely scenario is the
condensation of three ($N=3$) ${\bf Q}$ vectors related by 120∘ forming a
charge-$6e$ ‘Cooper molecule’ with total momentum $\sum_{i}{\bf Q}_{i}=0$. One
telltale demonstration of such Cooper molecule phase will be the resistance
oscillation near or slightly above the critical temperature in the spirit of
the Little-Parks experiment with flux periods of $h/(2Ne)$ Little and Parks
(1962); Parks and Little (1964). A recent demonstration of resistance
oscillation in the CVS superconductor with flux periods going from $h/2e$ to
$h/4e$ to $h/6e$ Ge et al. (2022) upon increasing the temperature is evidence
that this scenario may have already been realized.

In the original Little-Parks experiment, a thin cylinder with well-defined
radius was used to measure the resistance oscillation near the critical
temperature and the results were interpreted as periodic variations of the
free energy as a function of flux quantum $h/2e$. The penetration depth
$\lambda$ in the Little-Parks setup is considered much larger than the
thickness of the cylinder because of vanishing superconducting order. We may
call this the ’thin-rim geometry’. In Ref. Ge et al., 2022 a hole was created
at the center of the thin-film device similar to the Little-Parks, but the
conducting region is not a thin cylinder of well-defined radius, as can be
seen in the schematic reproduction of the experimental setup in Fig. 1. This
situation may be dubbed the ’thick-rim geometry’.

Figure 1: Schematic geometry of the device used in the experiment of Ref. Ge
et al., 2022. The white region is the hole area. Magnetic field is
perpendicular to the plane.

Furthermore, the resistance oscillations with flux periods of $h/4e$ and
$h/6e$ were observed at temperatures well above the zero-resistance
temperature $T_{c}$ where the magnetic field is expected to fully penetrate
the device. Since the Cooper pair experiences a varying range of enclosed flux
as it goes from one lead to the other in the thick-rim geometry of such
device, how does one decide on the effective area for the flux quantum? To
understand the conundrum, we refer to Ref. Kogan et al., 2004 where the free
energy has been calculated on the annular geometry with different inner ($a$)
and outer ($b$) radii. In the absence of any vortex (which is what we assume
as well), the free energy oscillation is governed by an effective area
(adapted from Eq. (39) of Ref. Kogan et al., 2004)

$\displaystyle A_{\rm eff}=\frac{\pi(b^{2}-a^{2})}{2\ln(b/a)},$ (2)

reducing to the well-known result $\pi a^{2}$ only in the $b\rightarrow a$
limit. In general this effective area is somewhere between $\pi a^{2}$ and
$\pi b^{2}$. For example, if $b=1.3a$ ($b=2a$) the effective area is predicted
to be 1.32 (2.16) times the area of the hole. This is the area that is
expected to determine the period of the original Little-Park experiment. Thus
it is somewhat surprising that in the experiment of Ref. Ge et al., 2022 the
effective area was deduced to be close to the area of the hole, in all three
temperature regimes where the flux quantizations are found to be in multiples
of $h/2e$, $h/4e$, or $h/6e$. In fact, the estimated area was slightly less
than the geometric area of the hole they have in the device. It seems that the
use of such area is not well justified a priori given the large width of the
conducting region. Vortices will not play a significant role in the experiment
either. The thickness ($t$) of the devices used in Ref. Ge et al., 2022 was in
the range of 10-20 nanometers Ge et al. (2022) while the penetration depth is
essentially infinite. This is the situation studied by Pearl Pearl (1964);
Kogan et al. (2001) where the screening of the external magnetic field takes
place over the Pearl length $\Lambda=2\lambda^{2}/t\gg\lambda$, larger than
the device dimension in Ref. Ge et al., 2022. (The $T=0$ estimate of the
penetration depth is already in excess of 300 nm Duan et al. (2021).) For the
problem under consideration, then, it is safely assumed that there is little
effect of vortex fluctuation and the magnetic field is nearly uniform. If so,
the question of a well-defined effective area for the flux quantization looms
large.

An important hint comes from the observation that in the experiment of Ref. Ge
et al., 2022 the periodicity of the oscillations is best defined at the
highest temperatures, where the resistivity is about half of the normal state
and the $h/6e$ flux quantum is observed. This places us firmly in the
fluctuating superconducting regime, i.e., the charge $6e$ order is far from
achieving long range order. The system is best treated as a metal with some
superconducting fluctuations. In contrast, at low temperatures where the
theory of Ref. Kogan et al., 2004 is expected to apply, only a few poorly
formed oscillations has been reported, so that the period cannot be accurately
determined. This motivates us to consider the fluctuating superconducting
regime as opposed to the ordered regime considered in Ref. Kogan et al., 2004.

It is the purpose of this paper to focus on the fluctuation regime and address
these questions in the framework of linearized time-dependent Ginzburg-Landau
(TDGL) theory Skocpol and Tinkham (1975). We do not address the origin or the
stability of $2Ne$-charged objects, but regard the effective charge $6e$ as an
input in the TDGL. Instead, what we address is a problem of general relevance
in all superconducting devices that are thin enough to be in the Pearl regime
and devoid of screening, at the same time being thick-rimmed to defy a naive
definition of the effective area as the geometric area of the hole.

##  II Formulation

The linearized time-dependent Ginzburg-Landau (TDGL) equation near $T_{c}$ is
Skocpol and Tinkham (1975)

$\displaystyle\left[\overline{\alpha}-\frac{\hbar^{2}}{2m}({\bm{\nabla}}-ie^{*}{\bf
A})^{2}\right]\Psi=-\overline{\alpha}\tau\frac{\partial\Psi}{\partial t}.$ (3)

The effective charge $e^{*}=6e$ is assumed in the TDGL. The effective mass $m$
refers to that of charge-$6e$ composite. The symbol $m^{*}$ will be reserved
for another definition of mass - see Eq. (8) - that will play a crucial role
on the analysis. The gap parameter $\overline{\alpha}$ is proportional to
$T-T_{c}$ and related to the coherence length $\xi$ by

$\displaystyle\overline{\alpha}=\frac{\hbar^{2}}{2m\xi^{2}}$ (4)

while the relaxation time $\tau$ is given by Skocpol and Tinkham (1975)

$\displaystyle\tau=\frac{\hbar}{8k_{B}(T-T_{c})}.$ (5)

In this paper we shall take $\tau$ as a phenomenological constant. The
coherence length of the bulk CVS superconductor is estimated to be 20-40 nm
Ortiz et al. (2020); Yu et al. (2021); Liang et al. (2021), comparable to the
thickness of the device used in Ref. Ge et al., 2022.

The TDGL equation of Eq. (3) can be solved in general form by means of path
integral technique. First one takes out $\overline{\alpha}$ through the
substitution

$\displaystyle\Psi=e^{-t/\tau}\psi,$ (6)

as the new $\psi$ satisfies a different TDGL equation

$\displaystyle\frac{\partial\psi}{\partial t}$
$\displaystyle=-\frac{\xi^{2}}{\tau}(-i{\bm{\nabla}}-e^{*}{\bf A})^{2}\psi$
$\displaystyle=-\frac{\hbar}{2m^{*}}(-i{\bm{\nabla}}-e^{*}{\bf A})^{2}\psi.$
(7)

This resembles the Schrödinger equation in imaginary time, and the effective
mass

$\displaystyle m^{*}\equiv\frac{\hbar\tau}{2\xi^{2}},$ (8)

differing from the mass of the charge-$6e$ Cooper pair $m$, shows up. As
mentioned in the introduction, vortex effects will be considered weak and
negligible.

Equation (7) is a diffusion equation with the vector potential ${\bf A}\neq
0$. When ${\bf A}=0$,

$\displaystyle\frac{\partial\psi}{\partial
t}=\frac{\hbar}{2m^{*}}{\bm{\nabla}}^{2}\psi\equiv D{\bm{\nabla}}^{2}\psi$ (9)

is solved by the diffusion kernel (in two dimensions)

$\displaystyle\psi({\bf r}t)$ $\displaystyle=\frac{1}{4\pi
Dt}\exp\left(-\frac{{\bf r}^{2}}{4Dt}\right),$ $\displaystyle\langle
x^{2}\rangle$ $\displaystyle=\langle y^{2}\rangle=2Dt=\frac{\hbar t}{m^{*}}.$
(10)

For ${\bf A}\neq 0$, the partial differential equation (7) is solved in
general as

$\displaystyle\psi({\bf r}t)=\int U({\bf r}t;{\bf r}_{0}t_{0})\psi({\bf
r}_{0}t_{0})d{\bf r}_{0}$ (11)

where $t_{0}$ is some earlier time to $t$, ${\bf r}_{0}$ covers the area of
the device (excluding the hole). The propagator $U$ is given in the path
integral form

$\displaystyle U({\bf r}t;{\bf r}_{0}t_{0})$ $\displaystyle=\int[{\cal D}{\bf
r}]\exp\left[i\frac{2\pi}{\Phi^{*}_{0}}\int_{{\bf r}_{0}}^{\bf r}{\bf A}\cdot
d{\bf r}-\frac{m^{*}}{2\hbar}\int_{t_{0}}^{t}\dot{\bf r}^{2}dt\right],$ (12)

where $\Phi^{*}_{0}=h/e^{*}$ is the flux quantum pertinent to our problem. The
sum $[{\cal D}{\bf r}]$ spans all the spacetime paths that begin at $({\bf
r}_{0}t_{0})$ and end at $({\bf r}t)$. The Aharonov-Bohm (AB) phase factor
appears despite the fact that we are solving the diffusion (imaginary-time
Schrödinger) equation. It is a common feature of geometric terms such as the
AB phase to appear as a complex phase even in the imaginary-time evolution of
the path integral.

The product of two wave functions at unequal times and positions is
$\Psi^{*}({\bf r}_{1}t_{1})\Psi({\bf r}_{2}t_{2})$ where

$\displaystyle\Psi({\bf r}_{2}t_{2})$
$\displaystyle=e^{-(t_{2}\\!-\\!t_{0})/\tau}\int U({\bf r}_{2}t_{2};{\bf
r}_{0}t_{0})\psi({\bf r}_{0}t_{0})d{\bf r}_{0}$ $\displaystyle\Psi^{*}({\bf
r}_{1}t_{1})$ $\displaystyle=e^{-(t_{1}\\!-\\!t_{0})/\tau}\int U^{*}({\bf
r}_{1}t_{1};{\bf r}^{\prime}_{0}t_{0})\psi^{*}({\bf r}^{\prime}_{0}t_{0})d{\bf
r}^{\prime}_{0}.$ (13)

Taking the initial time $t_{0}=t_{1}$ in both path integrals without loss of
generality gives

$\displaystyle\Psi({\bf r}_{2}t_{2})\Psi^{*}({\bf
r}_{1}t_{1})=e^{-(t_{2}-t_{1})/\tau}\times$ $\displaystyle~{}~{}~{}~{}~{}\int
U({\bf r}_{2}t_{2};{\bf r}_{0}t_{1})\psi({\bf r}_{0}t_{1})\psi^{*}({\bf
r}_{1}t_{1})d{\bf r}_{0}.$ (14)

The non-local conductivity $\sigma({\bf r}_{1},{\bf r}_{2})$ where ${\bf
r}_{1}$ and ${\bf r}_{2}$ are the two leads coupled to the superconducting
device is obtained from the current-current correlation function $\langle{\bf
J}({\bf r}_{1}t_{1})\cdot{\bf J}({\bf r}_{2}t_{2})\rangle$ Glazman et al.
(1992). In turn, the current-current correlation function follows from the
two-particle density matrix

$\displaystyle\Gamma({\bf r}^{\prime}_{1}t_{1},{\bf r}^{\prime}_{2}t_{2};{\bf
r}_{1}t_{1},{\bf r}_{2}t_{2})$ $\displaystyle~{}~{}=\langle\Psi^{*}({\bf
r}^{\prime}_{1}t_{1})\Psi({\bf r}_{1}t_{1})\Psi^{*}({\bf
r}^{\prime}_{2}t_{2})\Psi({\bf r}_{2}t_{2})\rangle$ (15)

by taking spatial derivatives

$\displaystyle\langle{\bf J}({\bf r}_{1}t_{1})\cdot{\bf J}({\bf
r}_{2}t_{2})\rangle\sim$ $\displaystyle~{}~{}({\bf D}_{1}\\!-\\!{\bf
D}^{*}_{1^{\prime}})\cdot({\bf D}_{2}\\!-\\!{\bf
D}^{*}_{2^{\prime}})\Gamma({\bf r}^{\prime}_{1}t_{1},{\bf
r}^{\prime}_{2}t_{2};{\bf r}_{1}t_{1},{\bf r}_{2}t_{2})|_{{\bf r}={\bf
r}^{\prime}}.$ (16)

The $\langle\cdots\rangle$ implies thermal and other averages. The covariant
derivatives are

$\displaystyle{\bf D}_{1,2}$ $\displaystyle={\bm{\nabla}}_{1,2}-ie^{*}{\bf
A}({\bf r}_{1,2})$ $\displaystyle{\bf D}^{*}_{1,2}$
$\displaystyle={\bm{\nabla}}_{1,2}+ie^{*}{\bf A}({\bf r}_{1,2}).$ (17)

We can express the two-point correlation function, using Eq. (14),

$\displaystyle\Psi^{*}({\bf r}^{\prime}_{1}t_{1})\Psi({\bf
r}_{1}t_{1})\Psi^{*}({\bf r}^{\prime}_{2}t_{2})\Psi({\bf r}_{2}t_{2})$
$\displaystyle~{}~{}~{}=e^{-2(t_{2}-t_{1})/\tau}\int U({\bf r}_{2}t_{2};{\bf
r}_{0}t_{1})U^{*}({\bf r}^{\prime}_{2}t_{2};{\bf
r}^{\prime}_{0}t_{1})\psi({\bf r}_{0}t_{1})\psi^{*}({\bf
r}^{\prime}_{1}t_{1})\psi^{*}({\bf r}^{\prime}_{0}t_{1})\psi({\bf
r}_{1}t_{1})d{\bf r}_{0}d{\bf r}^{\prime}_{0}.$ (18)

This is an exact representation of the two-point correlation function based on
the path-integral solution of TDGL. Eventually, one takes ${\bf
r}_{1,2}\rightarrow{\bf r}^{\prime}_{1,2}$ and the integrals over ${\bf
r}_{0}$ and ${\bf r}^{\prime}_{0}$ pick up the largest contributions from
${\bf r}_{0}\approx{\bf r}^{\prime}_{1}$ and ${\bf r}^{\prime}_{0}\approx{\bf
r}_{1}$. Under the approximation the two-particle correlation function
simplifies to

$\displaystyle\Gamma({\bf r}^{\prime}_{1}t_{1},{\bf r}^{\prime}_{2}t_{2};{\bf
r}_{1}t_{1},{\bf r}_{2}t_{2})$ $\displaystyle\rightarrow
e^{-2(t_{2}-t_{1})/\tau}\langle U({\bf r}_{2}t_{2};{\bf
r}^{\prime}_{1}t_{1})U^{*}({\bf r}^{\prime}_{2}t_{2};{\bf
r}_{1}t_{1})|\psi({\bf r}^{\prime}_{1}t_{1})|^{2}|\psi({\bf
r}_{1}t_{1})|^{2}\rangle$ $\displaystyle\propto
e^{-2(t_{2}-t_{1})/\tau}\langle U({\bf r}_{2}t_{2};{\bf
r}^{\prime}_{1}t_{1})U^{*}({\bf r}^{\prime}_{2}t_{2};{\bf
r}_{1}t_{1})\rangle.$ (19)

The first propagator $U({\bf r}_{2}t_{2};{\bf r}^{\prime}_{1}t_{1})$ is the
sum over all paths from $({\bf r}^{\prime}_{1}t_{1})$ to $({\bf r}_{2}t_{2})$
while the second propagator $U^{*}({\bf r}^{\prime}_{2}t_{2};{\bf
r}_{1}t_{1})$ is (complex conjugate of) the sum over the paths from $({\bf
r}_{1}t_{1})$ to $({\bf r}^{\prime}_{2}t_{2})$. In the next section, we
analyze the two-particle correlator in the semi-classical approximation with
the geometry of Fig. 1, similar to the one adopted in Ref. Ge et al., 2022, in
mind.

##  III Optimal path in the punctured geometry

In the punctured geometry resembling the device of the experiment as shown in
Fig. 2(a), paths can be divided as $+$ and $-$ depending on whether it passes
over or under the inner square. There are four possibilities of paths for the
product of propagators, $(++),(+-),(-+),(--)$. Keeping in mind that eventually
we will take ${\bf r}^{\prime}_{1,2}={\bf r}_{1,2}$, the complex phase factors
cancel out in the $(++)$ and $(--)$ cases but gives rise to the Aharonov-Bohm
(AB) phase in the $(+-)$ and $(-+)$ cases. For the interpretation of the
experiment in Ref. Ge et al., 2022 it suffices to consider the two terms

$\displaystyle\langle U({\bf r}_{2}t_{2};{\bf r}^{\prime}_{1}t_{1})U^{*}({\bf
r}^{\prime}_{2}t_{2};{\bf r}_{1}t_{1})\rangle$
$\displaystyle\rightarrow\langle U_{+}({\bf r}_{2}t_{2};{\bf
r}^{\prime}_{1}t_{1})U^{*}_{-}({\bf r}^{\prime}_{2}t_{2};{\bf
r}_{1}t_{1})\rangle$ $\displaystyle~{}~{}~{}+\langle U_{-}({\bf
r}_{2}t_{2};{\bf r}^{\prime}_{1}t_{1})U^{*}_{+}({\bf r}^{\prime}_{2}t_{2};{\bf
r}_{1}t_{1})\rangle.$ (20)

Figure 2: (a) Schematic geometry of the device used in Ref. Ge et al., 2022.
The leads shown in Fig. 1 are reduced to two point ${\bf r}_{1}$ and ${\bf
r}_{2}$ now. The superconducting region is confined between the inner and the
outer rectangles. The hole region is colored. Two geometrically shortest paths
connecting the leads at ${\bf r}_{1}$ and ${\bf r}_{2}$ are shown as $({\bf
r}_{1};{\bf r}_{3\pm};{\bf r}_{4\pm};{\bf r}_{2})$. Lengths of various line
segments are designated as $a,b,c,d$. (b) Geometrically shortest paths (in
blue) and the alternative, edge Aharonov-Bohm path (in red). (c) The
difference between blue and red paths is highlighted as the travel through the
hypotenuse or the base + perpendicular sides of the right triangle.

The propagators can be evaluated semi-classically. The exponential weight of
the kinetic energy implies that shorter paths carry the larger weight. Taking
${\bf r}_{1}$ and ${\bf r}_{2}$ to be the opposite centers at the edges of the
outer square, the shortest paths from ${\bf r}_{1}$ to ${\bf r}_{2}$ are the
ones going just above and just below the inner square as shown in Fig. 2. The
shortest paths are divided into three linear segments $({\bf r}_{1}t_{1};{\bf
r}_{3\pm}t_{3})$, $({\bf r}_{3\pm}t_{3};{\bf r}_{4\pm}t_{4})$, $({\bf
r}_{4\pm}t_{4};{\bf r}_{2}t_{2})$, and the corresponding classical action
$S_{c}=(m^{*}/2)\int\dot{\bf r}^{2}~{}dt$ evaluated for each:

$\displaystyle
S_{13}=\frac{m^{*}}{2}\frac{c^{2}}{t_{3}\\!-\\!t_{1}},S_{34}=\frac{m^{*}}{2}\frac{d^{2}}{t_{4}\\!-\\!t_{3}},S_{42}=\frac{m^{*}}{2}\frac{c^{2}}{t_{2}\\!-\\!t_{4}}.$
(21)

The total action $S=S_{13}+S_{34}+S_{42}$ is minimized when we choose the time
interval proportionately with distance,
$t_{3}-t_{1}:t_{4}-t_{3}:t_{2}-t_{4}=c:d:c$ and equals

$\displaystyle S({\bf r}_{1},{\bf
r}_{2})=\frac{m^{*}}{2}\frac{(2c+d)^{2}}{t_{2}-t_{1}}.$ (22)

The semi-classical propagators are

$\displaystyle\langle U_{\pm}({\bf r}_{2}t_{2};{\bf
r}^{\prime}_{1}t_{1})U^{*}_{\mp}({\bf r}^{\prime}_{2}t_{2};{\bf
r}_{1}t_{1})\rangle$
$\displaystyle\sim\langle\exp\left(i\frac{2\pi}{\Phi^{*}_{0}}\bigl{[}\int_{{\bf
r}^{\prime}_{1}}^{{\bf r}_{2}}{\bf A}\cdot d{\bf r}_{\pm}-\int_{{\bf
r}_{1}}^{{\bf r}^{\prime}_{2}}{\bf A}\cdot d{\bf r}_{\mp}\bigr{]}\right)$
$\displaystyle\times\exp\left(-S({\bf r}^{\prime}_{1},{\bf
r}_{2})/\hbar-S({\bf r}_{1},{\bf r}^{\prime}_{2})/\hbar\right)\rangle.$ (23)

Calculating the current-current correlation function using Eq. (16) is
straightforward and yields the oscillatory factor
$\cos(2\pi\Phi_{g}/\Phi^{*}_{0})$ where $\Phi_{g}$ is the flux enclosed inside
the geometric area spanned by the shortest path. This area is obviously larger
than that of the inner square and therefore disagrees with an interpretation
that associates the AB phase with the flux through the hole. Worse yet, the
flux $\Phi_{g}$ depends explicitly on the geometric dimension of $a$ so that
in a different device with the larger width $a$, a totally different estimate
of $\Phi_{g}$ is inevitable. If such were the case, an interpretation of the
resistance oscillation as due to flux periodicity is highly unreliable.

The semi-classical path consistent with the conventional AB interpretation
would be the one where both ${\bf r}_{3\pm}$ and ${\bf r}_{4\pm}$ hit the
center of the inner square, and the remainder of the path goes around its
perimeter as shown by red in Fig. 2(b). We refer to such alternative path as
the edge Aharonov-Bohm (edge-AB) path since it goes around the edges of the
hole and picks out its area. The difference between the edge-AB path and the
geometrically preferred path is whether one goes from ${\bf r}_{1}$ to ${\bf
r}_{3\pm}$ through the hypotenuse of the right triangle or through the base +
perpendicular edges of it, as shown in Fig. 2(c). Designating the lengths of
three edges as $a,b,c$, a simple back-of-the-envelop calculation shows that
the action for the two paths (called $c$-path and $ab$-path) are respectively

$\displaystyle S_{c}$ $\displaystyle=\frac{m^{*}}{2}\frac{c^{2}}{t}$
$\displaystyle S_{a}+S_{b}$
$\displaystyle=\frac{m^{*}}{2}\frac{a^{2}}{t^{\prime}}+\frac{m^{*}}{2}\frac{b^{2}}{t-t^{\prime}}$
$\displaystyle\rightarrow\frac{m^{*}}{2}\frac{(a+b)^{2}}{t},$ (24)

the last line following upon minimizing $S_{a}+S_{b}$ with respect to the
intermediate time $t^{\prime}$. Since $(a+b)^{2}>c^{2}$, we obtain
$S_{a}+S_{b}>S_{c}$ as expected. This simple result is as expected, but poses
a challenge in interpreting the experiment Ge et al. (2022). To get a
heightened sense of the paradox, imagine the limit $a\rightarrow\infty$ such
that the area enclosed by the $c$-path becomes very large. The flux enclosed
inside the geometric area will vary rapidly with the magnetic field, and the
oscillation will become unobservable. A more attractive scenario is that the
particle somehow chooses the edge-AB path and maintain the same area for
enclosing the flux quantum regardless of the size of $a$. But how?

Below we propose a quick fix to the conundrum and discuss its physical origin
in the following section. Instead of the bare action used in Eq. (28), adopt
the following effective action

$\displaystyle S_{a}$
$\displaystyle=\frac{m^{*}}{2}\frac{a^{2}}{t^{\prime}}+\alpha
a^{2}t^{\prime},$ $\displaystyle S_{b}$
$\displaystyle=\frac{m^{*}}{2}\frac{b^{2}}{t-t^{\prime}}+\beta
b^{2}(t-t^{\prime}),$ $\displaystyle S_{c}$
$\displaystyle=\frac{m^{*}}{2}\frac{c^{2}}{t}+\alpha c^{2}t.$ (25)

The extra $\propto t$ term in the effective action can be derived in a
standard manner as a result of fluctuation around the classical path (see
following section). The coefficients $\alpha$ and $\beta$ are in principle
different because in the case of $b$-path, the presence of the edge prohibits
the fluctuation to roughly half the amount compared to $a$\- or $c$-path where
the fluctuation is uninhibited. For now we treat them as parameters
$\alpha,\beta>0$ and minimize the new action $S_{a}+S_{b}$ with respect to
$t^{\prime}$, yielding

$\displaystyle\frac{m^{*}}{2}\frac{a^{2}}{t^{\prime
2}}-\frac{m^{*}}{2}\frac{b^{2}}{(t-t^{\prime})^{2}}=\alpha a^{2}-\beta b^{2}.$
(26)

An approximate solution is found by parameterizing
$t^{\prime}=t[a/(a+b)-\epsilon]$ and $t-t^{\prime}=t[b/(a+b)+\epsilon]$ and
finding, to first order in $\epsilon$:

$\displaystyle\epsilon\approx\frac{(\alpha a^{2}-\beta
b^{2})ab}{(a+b)^{4}}\frac{t^{2}}{m^{*}}.$ (27)

We get $\epsilon>0$ when $\alpha a^{2}>\beta b^{2}$. This makes sense, as the
particle wants to spend more time ($\epsilon>0$) in the $b$-region where the
linear-in-$t$ part of the action is less than that in the $a$-region. To the
same order in $\epsilon$ we find

$\displaystyle S_{a}+S_{b}-S_{c}$
$\displaystyle\approx\frac{m^{*}ab}{t}-\left(\alpha
ab+(\alpha-\beta)\frac{b^{3}}{a+b}\right)t$ $\displaystyle-\frac{(\alpha
a^{2}-\beta b^{2})^{2}ab}{(a+b)^{4}}\frac{t^{3}}{m^{*}}.$ (28)

For $\alpha>\beta$ the second term which is linear in $t$ is negative and will
dominate the first term for large enough $t$, making the $ab$-path more
favorable than the $c$-path despite its longer trajectory. Similarly the third
term is negative if $\beta$ is small enough compared with $\alpha$, only
speeding up the crossover from $c$\- to $ab$-path. We shall later show that in
fact $\beta$ appraches zero. This will be our route to stabilizing the long
path that hugs the edge. The expression derived in Eq. (28) invites the
definition of a time scale $t_{c}$:

$\displaystyle t_{c}=\sqrt{m^{*}/\alpha}.$ (29)

For $t\gtrsim t_{c}$, the linear-in-$t$ correction terms start to dominate
over the $\sim 1/t$ term and $S_{a}+S_{b}$ becomes less than $S_{c}$, favoring
the longer $ab$-path over the $c$-path on the hypotenuse. The existence of the
time scale $t_{c}$ beyond which the edge-AB path is favored provides
justification to the AB interpretation using the hole area and is the key
assertion of our analysis.

Figure 3: Fluctuating paths give rise to area fluctuations. $\Delta A>0$
($\Delta A<0$) fluctuations are shown as red (blue) areas.

##  IV Effective Action

Now we give physical justification to the $\propto t$ term in the effective
action, estimate the coefficients $\alpha,\beta$, and carve out the conditions
under which the crossover of paths from short to long can take place. As
schematically shown in Fig. 3, we consider fluctuations away from the
classical paths and their effects on the action at the Gaussian level. A well-
known correction to the classical action $S\sim(m^{*}/2)(l^{2}/t)$ is the
prefactor $t^{-d/2}$ where $d$ is the spatial dimension. Such correction,
though, will modify the analysis (which concerns quantities in the exponent)
only by logarithmic amount and not play a vital role. A more pertinent
correction comes from changes in the AB phase $\oint{\bf A}\cdot d{\bf r}$. As
shown in Fig. 3, fluctuation of the classical path implies the change in the
area, giving rise to the change in the phase factor

$\displaystyle\exp\left(\frac{2\pi iB}{\Phi^{*}_{0}}\Delta A\right)$ (30)

where $\Delta A$ is the area enclosed between the classical and the
fluctuating (non-classical) path - red and blue region in Fig. 3. For the
‘free’ region the area $\Delta A$ takes both positive and negative values.
This is not the case, however, if the classical path is against the wall, i.e.
it travels along the edges of the hole. In this case the area fluctuation is
always of one sign: $\Delta A>0$. For the two paths shown in Fig. 2(c), the
difference arises in going from ${\bf r}_{2}$ to ${\bf r}_{3\pm}$ and from
${\bf r}_{4\pm}$ to ${\bf r}_{2}$. In both cases, one only needs to compare
the effective actions for the hypotenuse and the base+perpendicular paths of
Fig. 2(c).

If the time $t$ is not too long, the typical fluctuating path will be ‘smooth’
and unlikely to cross the classical path. In that case the area $\Delta A$ can
be estimated reasonably well by the length $l$ of the classical path
multiplied by the width in the transverse direction. Given that the governing
equation (9) is a diffusion equation, the condition for this is $t<t_{D}$
where the diffusion time is given by

$\displaystyle t_{D}\sim l^{2}/D.$ (31)

The size of the transverse fluctuation can be estimated by using the one-
dimensional diffusion equation (as we are interested in the fluctuation in the
transverse direction to the classical path),

$\displaystyle\frac{\partial\psi}{\partial
t}=\frac{\hbar}{2m^{*}}\frac{\partial^{2}\psi}{\partial y_{\perp}^{2}}\equiv
D\frac{\partial^{2}\psi}{\partial y_{\perp}^{2}},$ (32)

for the transverse coordinate $y_{\perp}$. Using its well-known solution,
$\psi(y_{\perp},t)=(4\pi Dt)^{-1/2}\exp(-y_{\perp}^{2}/4Dt)$, $\langle
y_{\perp}^{2}\rangle=2Dt$, we estimate

$\displaystyle\Delta y_{\perp}$ $\displaystyle\sim\sqrt{2Dt}=\sqrt{\hbar
t/m^{*}}=\xi\sqrt{2t/\tau},$ $\displaystyle\Delta A$ $\displaystyle\sim
l\xi\sqrt{2t/\tau}$ (33)

where $l$ is the length of a particular classical trajectory and $t$ is the
time taken to traverse it.

Evaluation of the two-particle correlator, Eq. (23), ultimately depends on
summing over various paths, and this implies summation over various areas
$\Delta A$. We therefore introduce the Gaussian distribution function for the
area

$\displaystyle P(\Delta A)=\frac{e^{-(\Delta
A)^{2}/2\sigma_{A}^{2}}}{\sqrt{2\pi\sigma_{A}^{2}}},\left(\sigma_{A}=l\xi\sqrt{\frac{2t}{\tau}}\right),$
(34)

and sum over $\Delta A$ with this probability. For the area fluctuation
against the edge where $\Delta A>0$, the normalization factor is only half as
large: $\sqrt{2\pi\sigma_{A}^{2}}\rightarrow\sqrt{2\pi\sigma_{A}^{2}}/2$. The
following exact integral identities are useful in estimating the effect of
area fluctuation both without and with the wall:

$\displaystyle\int_{-\infty}^{\infty}\exp\left(2\pi i\frac{B\Delta
A}{\Phi^{*}_{0}}\right)P(\Delta A)d(\Delta A)$ $\displaystyle=e^{-\alpha
l^{2}t}$ $\displaystyle\int_{0}^{\infty}\exp\left(2\pi i\frac{B\Delta
A}{\Phi^{*}_{0}}\right)P(\Delta A)d(\Delta A)$ $\displaystyle=e^{-\alpha
l^{2}t}\left[1+i{\rm erfi}(\sqrt{\alpha t}l)\right].$ (35)

In the case of the first line, one can think of the $\alpha l^{2}t$ in the
exponent as an addition to the classical action. Comparing this to the
effective action in Eq. (25), we realize that this $\alpha$ is precisely the
phenomenological $\alpha$ there, but now it is explicitly given by

$\displaystyle\alpha=(2\pi B\xi/\Phi^{*}_{0})^{2}/\tau.$ (36)

The imaginary error function erfi($x$) in the second line of the identity has
the property

$\displaystyle{\rm
erfi}(x)\equiv\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}e^{y^{2}}dy=e^{x^{2}}D(x),$
(37)

and the Dawson function $D(x)$ has the asymptotic properties

$\displaystyle D(x)$ $\displaystyle\sim x~{}~{}~{}(x\ll 1)$ $\displaystyle
D(x)$ $\displaystyle\sim 1/2x~{}~{}~{}(x\gg 1),$ (38)

with the maximum value less than $1$ at $x\approx O(1)$. Using the Dawson
function one can re-write the half-infinite integral in the second line of Eq.
(35) as

$\displaystyle\int_{0}^{\infty}\exp\left(2\pi i\frac{B\Delta
A}{\Phi^{*}_{0}}\right)P(\Delta A)d(\Delta A)$ $\displaystyle~{}~{}=e^{-\alpha
l^{2}t}+iD(\sqrt{\alpha t}l).$ (39)

If we could ignore the Dawson function and keep the first term only, it means
the parameter $\beta$ equals $\alpha$. If instead the Dawson function
dominates at long times, the decay with time is only algebraic, and $\beta$ is
effectively zero. The crossover from $\beta\approx\alpha$ to $\beta\approx 0$
takes place around the time $t_{\beta}$ which is estimated by $\alpha
l^{2}t_{\beta}\approx 1$, or in physical terms,

$\displaystyle\frac{t_{\beta}}{\tau}=\left(\frac{\Phi^{*}_{0}}{2\pi B\xi
l}\right)^{2}.$ (40)

We return to Eq. (28) for $S_{a}+S_{b}-S_{c}$ and insert $\alpha$ from Eq.
(36). We consider time $t>t_{\beta}$ (to be justified later) so that we can
set $\beta=0$, to find

$\displaystyle S_{a}+S_{b}-S_{c}$
$\displaystyle\approx\frac{m^{*}ab}{t}\Bigl{(}1-\left[1+\frac{b^{2}}{a(a+b)}\right]\left(\frac{t}{t_{c}}\right)^{2}$
$\displaystyle~{}~{}-\frac{a^{4}}{(a+b)^{4}}\left(\frac{t}{t_{c}}\right)^{4}\Bigr{)}.$
(41)

Accordingly, the crossover from the $c$-path to the $ab$-path takes place when

$\displaystyle\frac{t}{t_{c}}\gtrsim\sqrt{\frac{a^{2}+ab}{a^{2}+ab+b^{2}}}.$
(42)

The factor on the r.h.s. is of order one for typical aspect ratios $a:b$, and
we obtain that for $t\gtrsim t_{c}$ the path crossover indeed takes place. We
also need the estimate of the crossover time $t_{c}$ in terms of the
relaxation time $\tau$ and other physical parameters

$\displaystyle\frac{t_{c}}{\tau}=\frac{\sqrt{m^{*}/\alpha}}{\tau}=\frac{1}{2\sqrt{2}\pi}\frac{\Phi^{*}_{0}}{B\xi^{2}}.$
(43)

Comparing with Eq. (40), we find the ratio

$\displaystyle\frac{t_{c}}{t_{\beta}}=\sqrt{2}\pi\frac{Bl^{2}}{\Phi^{*}_{0}}.$
(44)

Since $l$ is roughly half the linear dimension of the hole, we get
$Bl^{2}\approx BA/4$, $A$ being the area of the hole. Thus for a field $B_{1}$
defined by $B_{1}A=\Phi^{*}_{0}$ which corresponds to the first $AB$
oscillation, $t_{c}$ is already comparable or larger than $t_{\beta}$ and this
ratio increases for larger $B$. This justifies our use of $\beta=0$.

The analysis so far only concerned a part of the overall trajectory going from
${\bf r}_{1}$ to ${\bf r}_{2}$. For a complete analysis let us go back the
full edge-AB trajectory ${\bf r}_{1}$ to ${\bf r}_{2}$ and estimate its
effective action:

$\displaystyle S_{13}$
$\displaystyle=\frac{m^{*}}{2}\frac{(a+b)^{2}}{t_{3}\\!-\\!t_{1}}+\alpha\frac{a^{3}}{a+b}(t_{3}-t_{1}),$
$\displaystyle S_{34}$
$\displaystyle=\frac{m^{*}}{2}\frac{d^{2}}{t_{4}\\!-\\!t_{3}},$ $\displaystyle
S_{42}$
$\displaystyle=\frac{m^{*}}{2}\frac{(a+b)^{2}}{t_{2}\\!-\\!t_{4}}+\alpha\frac{a^{3}}{a+b}(t_{2}-t_{4}).$
(45)

The subscript in $S_{ij}$ means the trajectory from ${\bf r}_{i}$ to ${\bf
r}_{j}$. Note that the $ab$-path is already assumed in writing down $S_{13}$
and $S_{42}$ above. The linear-$t$ term is absent in $S_{34}$ since the
fluctuation there is against the wall and we already showed $\beta=0$ in that
case. An estimate of $t_{3},t_{4}$ minimizing the total action
$S=S_{13}+S_{34}+S_{42}$ is given by

$\displaystyle t_{3}-t_{1}$ $\displaystyle=t_{2}-t_{4}=\frac{a+b}{l_{\rm
eff}}(t_{2}-t_{1}),$ $\displaystyle t_{4}-t_{3}$
$\displaystyle=\frac{d}{l_{\rm eff}}(t_{2}-t_{1}),$ (46)

where $l_{\rm eff}=2(a+b)+d$. The total action is

$\displaystyle S$ $\displaystyle=S_{13}+S_{34}+S_{42}$
$\displaystyle=\frac{m^{*}}{2}\frac{l_{\rm
eff}^{2}}{t_{2}\\!-\\!t_{1}}+\left(\alpha\frac{2a^{3}}{l_{\rm
eff}}+\frac{1}{\tau}\right)(t_{2}-t_{1}),$ (47)

where we now include the exponential time factor $e^{-(t_{2}-t_{1})/\tau}$ as
part of the effective action.

The conductivity $\sigma({\bf r}_{1},{\bf r}_{2})$ follows from taking the
Fourier transform of the current-current correlation $\langle{\bf J}({\bf
r}_{1}t_{1})\cdot{\bf J}({\bf r}_{2}t_{2})\rangle$ with respect to
$e^{-i\omega(t_{2}-t_{1})}$ and taking $\omega\rightarrow 0$ limit. In turn,
$\langle{\bf J}({\bf r}_{1}t_{1})\cdot{\bf J}({\bf r}_{2}t_{2})\rangle$ is
dominated by the path that maximizes the propagator, Eq. (23), or minimizes
the total action $S$ in Eq. (47). The minimization of $S$ with respect to
$t_{2}-t_{1}$ gives

$\displaystyle t_{2}-t_{1}=t^{*}=\frac{t_{c}}{\sqrt{2\left[2a^{3}/l_{\rm
eff}^{3}+1/\alpha\tau l_{\rm eff}^{2}\right]}}$ (48)

where $t_{c}$ is the crossover time in Eq. (43). Among the two terms in the
denominator,

$\displaystyle\frac{1}{\alpha\tau l_{\rm
eff}^{2}}=\left(\frac{\Phi^{*}_{0}}{2\pi B\xi l_{\rm
eff}}\right)^{2}\approx\left(\frac{B_{1}}{B}\right)^{2}\left(\frac{A}{2\pi\xi
l_{\rm eff}}\right)^{2}.$ (49)

We estimate $\xi=2\times 10^{-8}$m and $l_{\rm eff}=5\times 10^{-7}$m. By way
of comparison, the area of the hole $A$ in the experiment is estimated to be
$2.76\times 10^{-14}$m2 , which corresponds to a linear dimension of about
$1.7\times 10^{-7}$m. Thus $l_{\rm eff}$ is about 3 times the dimension of the
hole. With the $2\pi$ factor in our favor, we find that this expression is
comparable or less than one already at $B=B_{1}$, the field of the first $AB$
oscillation and will decrease rapidly at larger magnetic field. On the other
hand, the first term in the denominator of Eq. (48) is independent of $B$. For
this reason we keep only the first term in the denominator and write

$\displaystyle t^{*}\approx t_{c}\sqrt{\frac{l_{\rm eff}^{3}}{4a^{3}}}.$ (50)

The value of $t^{*}$ is then used to estimate the time to go from ${\bf
r}_{1}$ to ${\bf r}_{3\pm}$ as

$\displaystyle t_{3}-t_{1}=\frac{a+b}{l_{\rm
eff}}t^{*}=t_{c}\sqrt{\frac{(2a+2b+d)(a+b)^{2}}{4a^{3}}}.$ (51)

Recall, as we argued earlier, that the time $t_{3}-t_{1}$ to complete the
$ab$-path needs to be larger than $t_{c}$ in order to justify the crossover
from $c$-path to the $ab$-path. Now let’s estimate $t_{3}-t_{1}$ explicitly.
For $a:b:d=1:1:1$, for instance, we get $t_{3}-t_{1}=\sqrt{5}t_{c}>t_{c}$. For
$a:b:d=1:1:2$ we get an even larger $t_{3}-t_{1}=\sqrt{6}t_{c}$. The current
device setup Ge et al. (2022) is therefore consistent with the statement that
the time to travel between ${\bf r}_{1}$ and ${\bf r}_{3\pm}$ is sufficiently
long to justify the choice of $ab$-path over the $c$-path. Note further than
this condition is geometric, i.e. does not depend on the strength of the
magnetic field $B$.

For the $b$\- and $d$-trajectories, we already concluded that there are no
corrections to the classical action because $\beta$ is effective zero.
However, the area covered by a given trajectory always satisfies $\Delta A>0$
and this might give rise to the increase in the effective area from the
geometric area of the hole, $2b\times d$. Contrary to the naive expectation,
this is not the case. This can be seen easily from the expression in Eq. (39),
which contains no oscillatory terms. To be more explicit, for the segments of
the path hugging the edge, we may use the half-integral identities in Eqs.
(39) and (38) to obtain the compound action

$\displaystyle\left(\frac{i}{2\sqrt{\alpha(t_{3}\\!-\\!t_{1})}b}\right)^{2}\left(\frac{i}{2\sqrt{\alpha(t_{4}\\!-\\!t_{3})}d}\right)^{2}\left(\frac{i}{2\sqrt{\alpha(t_{2}\\!-\\!t_{4})}b}\right)^{2}$
$\displaystyle\sim-\frac{1}{(t_{2}-t_{1})^{3}}.$ (52)

This factor, being purely real, will have no contribution to the AB phase
$\exp[(2\pi i/\Phi^{*}_{0})\oint{\bf A}\cdot d{\bf r}]$. Thus, contrary to
naive expectation, the effective area stays the same as the geometric area of
the inner square: $A_{\rm eff}=2b\times d$. The situation is analogous to what
happens in three dimensional quantum oscillations, where the oscillation
period is given by the extremal area, despite a broad distributions in the
cross-section area that deviates from the extremal area with the same sign.
Mathematically this is because the integral in the second line of Eq. (35 has
a sharp cutoff and the effective area is determined by that cut-off.

##  V Discussion

Several time scales were introduced and discussed throughout the paper. Among
them, the diffusion time $t_{D}$ was crucial in allowing us to use the
Gaussian probability function $P(\Delta A)$ for the area fluctuation that
increases diffusively with time $\Delta A\propto\sqrt{t}$, Eq. (33). Since
$D=\hbar/2m^{*}$ and $m^{*}=\hbar\tau/2\xi^{2}$,

$\displaystyle t_{D}$ $\displaystyle=l^{2}/D=(l/\xi)^{2}\tau$ (53)

The other two time scales, $t_{\beta}$ and $t_{c}$, were defined in Eq. (40)
and (43), respectively. Putting all three time scales together gives

$\displaystyle\frac{t_{\beta}}{\tau}:\frac{t_{c}}{\tau}:\frac{t_{D}}{\tau}$
$\displaystyle\approx\left(\frac{\Phi^{*}_{0}}{2\pi B\xi
l}\right)^{2}:\frac{\Phi^{*}_{0}}{2\sqrt{2}\pi
B\xi^{2}}:\left(\frac{l}{\xi}\right)^{2}$
$\displaystyle\approx\left(\frac{\Phi^{*}_{0}}{2\pi
Bl^{2}}\right)^{2}:\frac{\Phi^{*}_{0}}{2\sqrt{2}\pi Bl^{2}}:1.$ (54)

For instance, the condition $t_{D}\gtrsim t_{c}$ amounts to having
$1\gtrsim\Phi^{*}_{0}/(2\sqrt{2}\pi Bl^{2})$, which is nicely fulfilled if we
take $Bl^{2}$ to be the AB flux through the hole, so that $B\gtrsim B_{1}$ ,
i.e. at least one $AB$ oscillations have been seen. The inequality
$t_{c}\gtrsim t_{\beta}$ holds under the similar condition.

As mentioned in the introduction, this paper does not address the microscopic
justification of the charge-$6e$ formation in the CVS superconductor. We take
it for granted and write down the Ginzburg-Landau theory accordingly. The
starting point of our investigation is the observation that the experimental
setup of Ref. Ge et al., 2022 calls for a theoretical justification of the
effective area one should use in estimating the flux quantum. The effective
area used by the authors of the experiment suggests that the area of hole (the
inner square) is playing a dominant role, but on the basis of some earlier
literature, Kogan et al. (2001, 2004) it is not clear whether this is
justified. Our investigation shows that the use of the extremal area - the
area of the inner square - does get justified under the experimentally viable
circumstances $\Phi/\Phi^{*}_{0}\gtrsim 1$ ($\Phi$=flux through the hole) if
the pairing order is fluctuating and has only short range order.

Heuristically, the preferred path that dominate the contribution to the
correlation function is not always the geometrically preferrred path, but the
one that likes to ‘stick to the edge’ because there the fluctuation (and hence
the effective action) is less. This simple picture arising from our analysis
is obviously applicable to the more mundane, charge-$2e$ superconductors under
the similar device setup, as long as we are in the fluctuation regime. We
surmise that this result may be more general, as long as the system is
described by a sum over diffusing Feynman paths. The basic principle rests on
the analog of Eq. (35) which states that the average over a distribution of
$AB$ areas is dominated by the extremal area in the presence of a sharp cut-
off, and not by the average area. Thus the result that the $AB$ oscillation is
given by the area of the hole may also apply to the charge-$e$ transport in a
normal metal in the thick-rim geometry.

###### Acknowledgements.

J.H.H. was supported by NRF-2019R1A6A1A10073079. He also acknowledges
financial support from EPIQS Moore theory centers at MIT and Harvard, and
thanks Manhyung Han for preparing the figures. Hospitality extended by IBS-
CCES at Seoul National University during the completion of this work is
gratefully acknowledged. P.L. acknowledges the support by DOE office of Basic
Sciences Grant No. DE-FG02-03ER46076. We thank Ziqiang Wang for enlightening
discussion and comments on the manuscript.

## References

  * Ge et al. (2022) J. Ge, P. Wang, Y. Xing, Q. Yin, H. Lei, Z. Wang, and J. Wang, _Discovery of charge-4e and charge-6e superconductivity in kagome superconductor csv3sb5_ (2022), URL https://arxiv.org/abs/2201.10352.
  * Ortiz et al. (2019) B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski, M. Bordelon, J. S. Mangum, I. W. H. Oswald, J. A. Rodriguez-Rivera, J. R. Neilson, S. D. Wilson, et al., Phys. Rev. Materials 3, 094407 (2019), URL https://link.aps.org/doi/10.1103/PhysRevMaterials.3.094407.
  * Ortiz et al. (2020) B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, et al., Phys. Rev. Lett. 125, 247002 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.125.247002.
  * Chen et al. (2021) H. Chen, H. Yang, B. Hu, Z. Zhao, J. Yuan, Y. Xing, G. Qian, Z. Huang, G. Li, Y. Ye, et al., Nature 599, 222 (2021), ISSN 1476-4687, URL https://doi.org/10.1038/s41586-021-03983-5.
  * Zhou and Wang (2021) S. Zhou and Z. Wang, _Doped orbital chern insulator, chern fermi pockets, and chiral topological pair density wave in kagome superconductors_ (2021), URL https://arxiv.org/abs/2110.06266.
  * Agterberg and Tsunetsugu (2008) D. F. Agterberg and H. Tsunetsugu, Nature Physics 4, 639 (2008), ISSN 1745-2481, URL https://doi.org/10.1038/nphys999.
  * Berg et al. (2009) E. Berg, E. Fradkin, and S. A. Kivelson, Nature Physics 5, 830 (2009), ISSN 1745-2481, URL https://doi.org/10.1038/nphys1389.
  * Radzihovsky and Vishwanath (2009) L. Radzihovsky and A. Vishwanath, Phys. Rev. Lett. 103, 010404 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.103.010404.
  * Agterberg et al. (2011) D. F. Agterberg, M. Geracie, and H. Tsunetsugu, Phys. Rev. B 84, 014513 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.84.014513.
  * Little and Parks (1962) W. A. Little and R. D. Parks, Phys. Rev. Lett. 9, 9 (1962), URL https://link.aps.org/doi/10.1103/PhysRevLett.9.9.
  * Parks and Little (1964) R. D. Parks and W. A. Little, Phys. Rev. 133, A97 (1964), URL https://link.aps.org/doi/10.1103/PhysRev.133.A97.
  * Kogan et al. (2004) V. G. Kogan, J. R. Clem, and R. G. Mints, Phys. Rev. B 69, 064516 (2004), URL https://link.aps.org/doi/10.1103/PhysRevB.69.064516.
  * Pearl (1964) J. Pearl, Applied Physics Letters 5, 65 (1964), eprint https://doi.org/10.1063/1.1754056, URL https://doi.org/10.1063/1.1754056.
  * Kogan et al. (2001) V. G. Kogan, V. V. Dobrovitski, J. R. Clem, Y. Mawatari, and R. G. Mints, Phys. Rev. B 63, 144501 (2001), URL https://link.aps.org/doi/10.1103/PhysRevB.63.144501.
  * Duan et al. (2021) W. Duan, Z. Nie, S. Luo, F. Yu, B. R. Ortiz, L. Yin, H. Su, F. Du, A. Wang, Y. Chen, et al., Science China Physics, Mechanics & Astronomy 64, 107462 (2021), ISSN 1869-1927, URL https://doi.org/10.1007/s11433-021-1747-7.
  * Skocpol and Tinkham (1975) W. J. Skocpol and M. Tinkham, Reports on Progress in Physics 38, 1049 (1975), URL https://doi.org/10.1088/0034-4885/38/9/001.
  * Yu et al. (2021) F. H. Yu, T. Wu, Z. Y. Wang, B. Lei, W. Z. Zhuo, J. J. Ying, and X. H. Chen, Phys. Rev. B 104, L041103 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.104.L041103.
  * Liang et al. (2021) Z. Liang, X. Hou, F. Zhang, W. Ma, P. Wu, Z. Zhang, F. Yu, J.-J. Ying, K. Jiang, L. Shan, et al., Phys. Rev. X 11, 031026 (2021), URL https://link.aps.org/doi/10.1103/PhysRevX.11.031026.
  * Glazman et al. (1992) L. I. Glazman, F. W. J. Hekking, and A. Zyuzin, Phys. Rev. B 46, 9074 (1992), URL https://link.aps.org/doi/10.1103/PhysRevB.46.9074.


# The transition from non-convex region to hyperbolic region and vice versa in
the Beryllium-Graphene-Beryllium three-layer structure: Plasmonics

Badie Ghavami<EMAIL_ADDRESS>School of Nano Science, Institute for
Research in Fundamental Sciences (IPM), Tehran 19395-5531, Iran  Elnaz
Rostampour  Department of Physics, Urmia University, Urmia, Iran

###### Abstract

Meta-materials are one of the important discussions of condensed matters which
have unusual optical properties. They can be used to regulate light and wave
propagation. Natural hyperbolic and non-convex plasmons are the
characteristics of meta-materials that are observed in these materials. In
this paper, we have studied Beryllium-Graphene-Beryllium (BeCBe) three-layer
optical properties and have shown this material has hyperbolic plasmon in the
THz frequency region and non-convex plasmons. We have shown most of the
plasmon energy propagation along the $q_{y}$ axis for the hyperbolic region.
Also, it can be seen that there is an elliptic case only in 1029.71 THz
frequency. The starting point of the hyperbolic region occurs at the frequency
of 289.31 THz, which is the first and the narrowest hyperbolic window. Due to
the existence of far more free electrons in the BeCBe plasmonic material, the
presence of the visible region in the third hyperbolic window causes a high
plasma frequency.

non-convex plasmons, hyperbolic plasmon, dielectric function, electron energy
loss spectrum

##  I Introduction

The plasmon-polariton is an attractive optical phenomenon found in 2D
materials. By note in plasmonic properties, we have four classes by attention
imaginary conductivity part: isotropic, anisotropic, hyperbolic, and chiral
plasmon-polaritons [1, 2, 3, 4, 5, 6]. In these classes, hyperbolic materials
interest us, because they enable a wide range of applications that include
far-field subwavelength imaging, nanolithography, emission engineering [7],
negative index waveguides [8], subdiffraction photonic funnels [9], and
nanoscale resonators [10]. Hyperbolic materials show dispersion and combine
properties of transparent reflective metals and dielectrics that ellipse
materials don’t show these properties [7]. Therefore, the behavior of
hyperbolic materials causes them to be considered. Hyperbolic materials are
highly anisotropic materials that are principal components of real parts of
their dielectric (imaginary part of conductivity) tensor having opposite signs
[11, 12]. They can naturally support very limited hyperbolic polarization [13,
14, 15]. For example, it could be shown for graphene when the imaginary part
of the surface conductivity is negative (hyperbolic properties), it becomes
transverse electric surface polaritons in the frequency region [16].
Hyperbolic metasurface are different conditions possible, an instance of, it
can be realized for certain LC-type contours in the microwave frequency region
[17, 18] and formed by anisotropic plasmonic particles 2D lattice in the
infrared region (graphene nanoribbons) [19]. In the terahertz (THz) range,
hyperbolic polaritons have observed MgB2, graphite, cuprate, and ruthenate
[20]. Also, in the range of mid-infrared (mid-IR), for instance, hyperbolic
polaritons were observed in the Van der Waals crystal [21]. Provide natural
hyperbolic polaritons [22] and tunable hyperbolic plasmons [23] in the THz and
mid-IR in the 2D materials observed. Black phosphorus (BP) is one of the 2D
materials to study new plasmonic designs which have tunable optical properties
and are flexible [24, 25, 26], because of having high tunable direct bandgap
and carrier mobility in-plane anisotropy. Pure phosphorous didn’t have
hyperbolic plasmons in the near-infrared but they made hyperbolic plasmons
with external modulation such as carbon phosphide [23].
Despite the hyperbolic behavior of the metamaterial in the hyperbolic
metamaterial, the vertical component of the effective permittivity depended on
the number of layers of the graphene sheets [27, 28, 29, 30]. In a multilayer
periodic structure composed of graphene layers in the THz frequency range, the
elliptic to hyperbolic regime existed using an external gate voltage [31]. In
the material examined here, the presence of Be atom makes us not have the
elliptic regime, and there is an elliptic case only in frequency 1029.71 THz.
Unlike the studied case in three layers of graphene, where there was an
elliptical region.
In a multi-layer graphene-dielectric composite material, backward waves and
forward waves were propagated under the hyperbolic and elliptical dispersion
regimes, respectively [32]. The formation of hybrid polaritons in the
multilayer structures composed of graphene and hexagonal boron nitride film
was more in hyperbolic regions which despite the intrinsic natural hyperbolic
behavior of hexagonal boron nitride, there was a more effective hyperbolic
behavior [33]. In the multilayer structure of graphene and Al2O3, elliptic to
hyperbolic dispersion was observed at a wavelength of 4.5 $\mu m$ by an
infrared ellipsometer[34].
This paper is written as follows: In section II, we have briefly explained the
theory of dielectric function and optical properties of 2D materials, and in
section III is explained the computational method by details that made the ab-
initio method as DFT calculation. Section IV becomes a discussion result of
the BeCBe three-layer and investigates plasmonic properties for variable
frequency. Finally, the summary and conclusion are provided in section V.

##  II Computational details

We have calculated the dielectric function in reciprocal space representation
for ground-state electronic structure. Here, We will briefly explain the
theory of calculations. The dielectric function is dependent on wave vector
and frequency which is written as $\epsilon_{G,G^{\prime}}(q,\omega)$, where
$G$, $q$, and $\omega$ are reciprocal lattice vector, momentum, and frequency
of energy transfer in excitation energy respectively. Therefore, the
macroscopic dielectric function is given by
$\epsilon_{M}=[\epsilon_{00}(q,\omega)]^{-1}$ taking into account the local
field correction [35], the absorption spectrum is given by
$\alpha=\Im\epsilon_{M}(q\to 0,\omega)$, and electron energy loss spectrum
(EELS) is gotten by $eels=-\Im[\epsilon^{-1}(q,\omega)]$.
The dielectric tensor of materials is:

$\epsilon(q,\omega)=\begin{pmatrix}\epsilon_{xx}(q,\omega)&0\\\
0&\epsilon_{yy}(q,\omega)\\\ \end{pmatrix}$ (1)

where $\epsilon(q,\omega)=\Re\epsilon(q,\omega)+i\Im\epsilon(q,\omega)$ or the
conductivity tensor equal to:

$\sigma(q,\omega)=\begin{pmatrix}\sigma_{xx}(q,\omega)&0\\\
0&\sigma_{yy}(q,\omega)\\\ \end{pmatrix}.$ (2)

To study surface mode properties, it is considered k-surface,
$\omega(q_{x},q_{y})$=constant, to be shown in Fig.5. The mode wave-front
propagation direction in anisotropic materials is defined by wave-vector
$\bm{q}$ which doesn’t generally match with the energy flow direction and it
is defined by the group velocity $\bm{v_{gr}}=\nabla_{q}\omega(q)$ [36, 1].
The velocity group represents the surface mode in anisotropic materials that
transports energy in the orthogonal direction in the k-surfaces. Thus, the
following equation can be used to explain this phenomenon:

$k_{0}\sqrt{q_{x}^{2}+q_{y}^{2}-k_{0}^{2}}[\frac{c}{2\pi}-\frac{2\pi}{c}\Im\sigma_{xx}\Im\sigma_{yy}]=\\\
\Im\sigma_{xx}(q_{x}^{2}-k_{0}^{2})+\Im\sigma_{yy}(q_{y}^{2}-k_{0}^{2})$ (3)

where $k_{0}=\omega/c$ is the free-space wave number ($\omega$ and $c$ are the
angular frequency and free-space speed of light respectively). In the
following sections, we have investigated anisotropic plasmons and hyperbolic
plasmons which are of the special anisotropic type for BeCBe’s three layers.

##  III Computational Method

The ab-initio method has done with density functional theory(DFT) calculations
which were performed by ’GPAW’ code [37] for 16 bands on the Brillouin zone is
sampled by 14 $\times$ 14 $\times$ 1 k-point mesh of the Monkhorst-Pack [38].
The cutoff energy is obtained at 500 eV. The exchange-correlation(XC) energy
functional ’optPBE-vdW’ [39, 40, 41] is implemented self-consistently in GPAW
considering the Van der Waals potentials. Also, it has been utilized XC
functional with Broyden mixing [42] as the exchange-correlation energy along
with the gride-based projection-augmented wave(PAW) to determine the
electronic and optical properties of the BeCBe three-layer structure. The
structure is optimized by Hellmann-Feynman forces, which are smaller than 0.01
eV/$A^{\circ}$ on each atom. The optical properties are calculated by
calculating the polarizability response function using the framework of the
random phase approximation(RPA). The response function has been evaluated
using the Gritsenko-Leeuwen-Lenthe-Baerends-Solis-Correlation potential(GLLB-
SC) [43] and mesh grid for plasmonic and absorption is utilized 80 $\times$ 70
$\times$ 1 $q$-grid. Finally, we have investigated non-convex and hyperbolic
regions for the electrostatically BeCBe three-layer based on Maxwell’s
boundary condition and the DFT method.

##  IV Result and Discussion

Electronic structure of BeCBe: In Fig.1, the system under consideration
consists of three layers of Beryllium, Graphene, and Beryllium. There is 3.68
$A^{\circ}$ between each layer after optimizing energy. We have studied the
configuration band structure using the DFT method. Fig.2 has shown the band
structure of the system.

Figure 1: (Color online) It is shown (a) the geometric of the BeCBe three
layer, (b) top view structure, and (c) distance between layers is
$3.68A^{\circ}$, respectively. Here, Berelium and Carbon atoms are shown with
green and gray colors, respectively. The unit cell parameters are
$a=2.398A^{\circ}$, $b=4.151A^{\circ}$, $c=25.0A^{\circ}$, and
$\alpha=\beta=\gamma=90^{\circ}$ after optimization.  Figure 2: (Color online)
Band structure of the BeCBe along high symmetry paths of the Brillouin zone.
Fermi energy is $E_{F}=-1.35$ eV.

According to the band structure and density of states (DOS)- Fig.3\- diagrams,
it is observed that the gap energy is zero and this structure has metal
properties.

Figure 3: (Color online) Density of states (DOS) diagram BeCBe. It shows to be
metal bulk.

Near the valence band maximum (VBM) and the conduction band minimum (CBM) at
the $\Gamma$ point and X point, we have found the different effective mass
values ($m^{\star}_{x}$, $m^{\star}_{y}$) for the x-direction and y-direction.
We get effective mass through the following equation [44] for two directions:

$m^{\star}_{x(y)}=\hbar^{2}(\frac{\partial^{2}E}{\partial^{2}k_{x(y)}})^{-1}$
(4)

where fermi energy is $E_{F}=-1.35$ eV, and the anisotropy ratio is defined
by[45]:

$\frac{\omega_{q,x}}{\omega_{q,y}}=\sqrt{\frac{m^{*}_{y}}{m^{*}_{x}}}$ (5)

where $m^{*}_{x(y)}$ and $\omega_{qx(y)}$ are effective mass and plasmon
angular frequency, respectively. This equation is a good scale to determine
plasmonic and polariton-type regimes. Thus, we would have expected to have
plasmon hyperbolic in the system.
In Fig.2 exist two linear bands forming Dirac point which are crossed between
$\Gamma$ to $X$ point above the fermi energy level, further $Y$ to $S$ point
below the fermi energy level. The crossing linear bands are both due to the
$p_{z}$ orbitals. The band gap in paths $\Gamma$ to $X$ point, $Y$ to $\Gamma$
point, and $\Gamma$ to $S$ point are zero and show material properties in
these paths. To paths $X$ to $S$ and $S$ to $Y,$ the band structure has gape
energy as in these paths’ system is a semiconductor. Also, this view can
expect to view plasmonic hyperbolic phenomena that the surface of the system
is shown semiconductor and the bulk structure is metal. In Fig.3, DOS has
shown to be a metal bulk structure.

Figure 4: (Color online) Imaginary and real Dielectric function of the BeCBe
three-layer for two directions x, y. The plot shows plasmonic different
regions in the THz frequency.

Optical properties: We have studied optical properties in our structure. For
this aim, we calculated and illustrated imaginary and real dielectric
functions in two directions x and y. Fig.4 have shown the dielectric function
of the BeCBe structure. The polariton behavior in a metamaterial and
metasurface depends on the real part of the dielectric function sign for each
direction as for the isotropic plasmonic, the anisotropic plasmonic, the
hyperbolic plasmons, and the chiral plasmons are given by:
$\epsilon_{xy}=\epsilon_{yx}=0$, $\epsilon_{xx}=\epsilon_{yy}$[46, 47, 2, 48],
$\epsilon_{xy}=\epsilon_{yx}=0$,
$\Re\epsilon_{xx}\times\Re\epsilon_{yy}>0$[36, 49],
$\epsilon_{xy}=\epsilon_{yx}=0$,
$\Re\epsilon_{xx}\times\Re\epsilon_{yy}<0$[50],
and $\epsilon_{xy}=-\epsilon_{yx}\neq 0$, $\epsilon_{xx}=\epsilon_{yy}$[51],
respectively.

This property can be mentioned based on the tensor of conductive elements:

$\epsilon(q_{i},\omega)=1+j\frac{2\pi q_{i}\sigma_{ii}(\omega)}{\omega}$ (6)

where $i=x,y$ and $j=\sqrt{-1}$.
We show dielectric function vs frequency in Fig. 4. For the free-moving
electrons, the polarization response is in opposite direction to the electric
field therefore the real part of the dielectric tensor of materials is
negative below the plasma frequency. The hyperbolic regime is usually the
component of the negative dielectric tensor in only one or two spatial
directions.

Figure 5: (color online) Diagram of (a) hyperbolic and (b) anisotropic in a
few frequencies. When change sign of the dielectric function real part (the
conductivity imaginary part), there is an interband transition. Figure 6:
(color online) Selected hyperbolic regime. Graph elongation factor vs
frequency (red). Hyperbolas are shown by the angle between the group velocity
and the vertical movement axis.

Optical properties of BeCBe: The imaginary and real parts of the dielectric
function in the $x$ and $y$ directions for BeCBe are drawn in Fig.4, which
have different signs in the specified frequency range. The obtained results
show that BeCBe has three hyperbolic frequency windows at the frequency of
$0-2000$ THz, one of which is wide and the other two are narrow, and
$\Re\epsilon_{xx}<0$, $\Re\epsilon_{yy}>0$ in them. Hyperbolic frequency
windows are shown with blue shadows. The first hyperbolic window is at the
frequency of $289.31-294.64$ THz, the second hyperbolic window is at the
frequency of $316.06-326.83$ THz, and the third hyperbolic window is at the
frequency of $343.05-632.47$ THz. The behavior of the real parts of the
dielectric function is as follows: The value of $\Re\epsilon_{xx}$ at
frequencies less than $632.47$ THz is always negative, while the value of
$\Re\epsilon_{yy}$ at frequencies less than $289.31$ THz is negative and is
placed in the non-convex region.
With the increase of frequency, the value of $\Re\epsilon_{yy}$ is positive
and leads to the first hyperbolic region, the value of which reaches zero at
the frequency of $294.64$ THz, and with the increase of frequency, it becomes
negative and then reaches zero. In this case, the non-convex region is formed.
At the frequencies of $316.06$ THz to $326.83$ THz, the value of
$\Re\epsilon_{yy}$ is positive and the second hyperbolic region is formed,
with the increase of the frequency, the value of $\Re\epsilon_{yy}$ becomes
negative and then reaches zero at the frequency of $343.05$ THz, which forms
the non-convex region. Again, with the increase of frequency, the value of
$\Re\epsilon_{yy}$ is positive up to the frequency of $1023$ THz, and the
value of $\Re\epsilon_{xx}$ becomes positive from the frequency of $632.47$
THz onwards and forms the hyperbolic region up to the frequency of $632.47$
THz. In fact, when $\Re\epsilon_{xx}$ and $\Re\epsilon_{yy}$ have the same
sign, an ellipse region is formed, and when the sign of $\Re\epsilon_{yy}$
changes, a hyperbolic region is formed. The non-convex regions have an almost
eight-like shape. Due to the fact that inductive and capacitive behaviors
occur at the frequencies between the transition points - the hyperbolic
frequencies-, in this material, in the region where $\Re\epsilon_{xx}$ is
negative, we have the inductive response and in the region where
$\Re\epsilon_{yy}$ is positive, we have the capacitive response. The change in
the sign of the dielectric function in this material indicates the change in
the behavior of polaritons. For BeCBe, the first hyperbolic window and the
second hyperbolic window are in the infrared region. The third hyperbolic
window shows the transition from the infrared region to the visible region. In
this material, $\Re\epsilon_{xx}$ reaches zero later than $\Re\epsilon_{yy}$.
The monolayer of metal-free two-dimensional carbon phosphide has one
hyperbolic frequency window and this window is wider compared to the third
BeCBe hyperbolic frequency window [23]. On the other hand, the width of the
third hyperbolic frequency window of BeCBe is more than $MgB_{2}$ [52], thin
films of $WTe_{2}$ [53], and monolayer black phosphorus [54]. So, the order of
the hyperbolic frequency window of these five materials can be written as
Monolayer carbon phosphide $>$ BeCBe $>$ $MgB_{2}$ $>$ monolayer black
phosphorus $>$ thin films of $WTe_{2}$.
Fig.5 shows isofrequency counters of surface plasmon polaritons. BeCBe has
frequencies as $\omega_{h1}=289.31$ THz, $\omega_{h2}=283.999$ THz,
$\omega_{h3}=305.33$ THz, $\omega_{h4}=316.06$ THz, $\omega_{h5}=322.22$ THz,
and $\omega_{h6}=343.05$ THz in the non-convex region. These frequencies are
vertical non-convex contours. In the hyperbolic region, frequency
$\omega_{h\nu p}=436.60$ THz is in hyperbolic form with large openings in
Fig.5. In fig.5, these frequencies are in hyperbolic form with large openings.
In this metamaterial, with the change of the sign of the real part of the
dielectric function, there are interband transitions. The transition
frequencies are as Table 1.

Table 1:  Transmission frequencies (THz) are shown for the hyperbolic plasmonic shape with a large opening. $\omega_{tt1}$ | $\omega_{tt2}$ | $\omega_{tt3}$ | $\omega_{tt4}$ | $\omega_{tt5}$ | $\omega_{tt6}$
---|---|---|---|---|---
289.31 | 294.64 | 316.06 | 326.83 | 343.05 | 632.46

In fig.5, these frequencies are in hyperbolic form with small openings. In
case of transition between interband and intraband, hyperbolic dispersion
occurs.

Figure 7: (color online) The plasmon dispersion for the $q_{x}$ and $q_{y}$
directions. Plasmon dispersion is more in the $y$ direction than in the $x$
direction. Plasmons dispersion in the $x$ direction to $q_{x}=1A^{\circ-1}$
and in the $y$ direction to $q_{y}=1.5A^{\circ-1}$ and then enter the single-
particle excitations region. Here, it is seen $\omega\propto\sqrt{q}$.

In fig.5, with the increase in frequency, the curves change from non-convex to
hyperbolic, from hyperbolic to non-convex, and so on. The degree of elongation
of the isofrequency contours is expressed by the elongation factor. The
elongation factor determines the non-convex and hyperbolic regions. At
frequencies where the elongation factor is positive, we have the non-convex
regime, and at frequencies where the elongation factor is negative, we have
the hyperbolic regime. In Fig. 6, the elongation factor is plotted in terms of
frequency. According to the obtained results, the elongation factor of this
metamaterial has the maximum value in three vertical frequencies as
$\omega_{\nu 1}=282.9$ THz, $\omega_{\nu 2}=310.5$ THz, and $\omega_{\nu
3}=334.5$ THz. These frequencies have a dumbbell shape as shown in Fig.5.
Before this frequency and after this frequency, frequencies are placed in the
non-convex regime. At frequencies around $\omega_{\nu}$
($\Re\epsilon_{yy}\gg\Re\epsilon_{xx}$), the phenomenon of low-loss
canalization occurs along the $q_{y}$ direction.
In the transition frequencies $\omega_{tti}$ for $i=1..6$
($\Re\epsilon_{xx}<0$, $\Re\epsilon_{yy}\approx 0$), we have the phenomenon of
canalization and most of the limited plasmon energy is propagated along the
$q_{y}$ axis. In this case, the direction of group velocity and plasmon energy
propagation coincide. This situation is more in hyperbolic isofrequency than
in non-convex isofrequency. The propagation of surface plasmon polaritons is
between the asymptotes of the hyperbola $q_{y}=\pm q_{x}tan(\theta)$.
Hyperbolas are expressed by the angle between the velocity of the group
$v_{g}$ and the axis of the vertical movement $q_{y}$
($\theta=tan^{-1}\sqrt{\frac{\Re\epsilon_{xx}}{\Re\epsilon_{yy}}}$) as shown
in Fig.6. Non-convex and hyperbolic surface plasmon polaritons are expressed
by equation 3. In the range of $0-632.47$ THz, $\Re\epsilon_{xx}$ and
$\Re\epsilon_{yy}$ do not cut each other
($\Re\epsilon_{xx}\neq\Re\epsilon_{yy}$) and there is no circular crossing
frequency. So the elongation factor near one is in the form of vertical lines.
According to the obtained results, the behavior of polaritons in BeCBe is
anisotropic plasmonic in non-convex regions and hyperbolic plasmons in
hyperbolic regions. The hyperbolic regime indicates directional plasmons. In
Fig.7, plasmons propagate in the $x$ and $y$ directions. Plasmon dispersion is
more in the $y$ direction than in the $x$ direction. $\Re\epsilon_{yy}$ is
equal to zero at different frequencies. Plasmons propagate in the $x$
direction to $q_{x}=1A^{\circ-1}$ and in the $y$ direction to
$q_{y}=1.5A^{\circ-1}$ and then enter the single-particle excitations region.
In BeCBe, small acoustic plasmons with very low intensity exist in both the
$q_{x}$ and $q_{y}$ directions at very low energies. BeCBe exhibits strong
anisotropic plasmon.
As we expect, at low $q_{x}$ and $q_{y}$, with the proximity of the surface
plasmon polariton to the Sommerfeld-Zenneck wave [55], it has a photon-like
behavior, and with the increase of $q_{x}$ and $q_{y}$ and the bending of the
dispersion relation curve, it approaches the asymptotic limit of the surface
plasma frequency.

##  V Conclusions

In summary, hyperbolic and non-convex plasmonic properties are revealed in
this material so that the results show three hyperbolic frequency windows. The
widest hyperbolic frequency window is located in the frequency range of
$343.05-632.47$ THz. We investigate the behavior of polaritons in non-convex
and hyperbolic regions with a two-dimensional description of the dielectric
function. The opposite sign of the real parts of the dielectric function in
both $x$ and $y$ directions indicates hyperbolic plasmon behavior and their
same sign indicates anisotopic plasmonic behavior. Plasmons are dispersed in
both $x$ and $y$ directions for this hyperbolic metamaterial. In the
hyperbolic region, most of the plasmon energy propagates along the $q_{y}$
axis. In the non-convex regions around the frequencies $\omega_{\nu 1}$,
$\omega_{\nu 2}$, and $\omega_{\nu 3}$ the phenomenon of low-loss canalization
in the $q_{y}$ direction is created. The presence of non-Euclidean geometry in
BeCBe makes it more practical.

## References

  * Nemilentsau  _et al._ [2019] A. Nemilentsau, T. Low, and G. Hanson, Chiral and hyperbolic plasmons in novel 2-d materials, in  _Carbon-Based Nanoelectromagnetics_ (Elsevier, 2019) pp. 119–138.
  * Garcia de Abajo [2014] F. J. Garcia de Abajo, Graphene plasmonics: challenges and opportunities, Acs Photonics 1, 135 (2014).
  * Grigorenko  _et al._ [2012] A. N. Grigorenko, M. Polini, and K. Novoselov, Graphene plasmonics, Nature photonics 6, 749 (2012).
  * Gangaraj  _et al._ [2016] S. A. H. Gangaraj, T. Low, A. Nemilentsau, and G. W. Hanson, Directive surface plasmons on tunable two-dimensional hyperbolic metasurfaces and black phosphorus: Green’s function and complex plane analysis, IEEE Transactions on Antennas and Propagation 65, 1174 (2016).
  * Mojarro  _et al._ [2022] M. Mojarro, R. Carrillo-Bastos, and J. A. Maytorena, Hyperbolic plasmons in massive tilted two-dimensional dirac materials, Physical Review B 105, L201408 (2022).
  * Rostampour [2021] E. Rostampour, The effect of nitrogen and boron on magneto-optical and dielectric properties of the azafullerenes c60-nxn (x= b, n; n= 1, 2), Solid State Communications 324, 114145 (2021).
  * Drachev  _et al._ [2013] V. P. Drachev, V. A. Podolskiy, and A. V. Kildishev, Hyperbolic metamaterials: new physics behind a classical problem, Optics express 21, 15048 (2013).
  * Podolskiy and Narimanov [2005] V. A. Podolskiy and E. E. Narimanov, Strongly anisotropic waveguide as a nonmagnetic left-handed system, Physical Review B 71, 201101 (2005).
  * Govyadinov and Podolskiy [2006] A. A. Govyadinov and V. A. Podolskiy, Metamaterial photonic funnels for subdiffraction light compression and propagation, Physical Review B 73, 155108 (2006).
  * Yao  _et al._ [2011] J. Yao, X. Yang, X. Yin, G. Bartal, and X. Zhang, Three-dimensional nanometer-scale optical cavities of indefinite medium, Proceedings of the National Academy of Sciences 108, 11327 (2011).
  * Poddubny  _et al._ [2013] A. Poddubny, I. Iorsh, P. Belov, and Y. Kivshar, Hyperbolic metamaterials, Nature photonics 7, 948 (2013).
  * Ferrari  _et al._ [2015] L. Ferrari, C. Wu, D. Lepage, X. Zhang, and Z. Liu, Hyperbolic metamaterials and their applications, Progress in Quantum Electronics 40, 1 (2015).
  * Gjerding  _et al._ [2017] M. N. Gjerding, R. Petersen, T. G. Pedersen, N. A. Mortensen, and K. S. Thygesen, Layered van der waals crystals with hyperbolic light dispersion, Nature communications 8, 1 (2017).
  * Narimanov and Kildishev [2015] E. E. Narimanov and A. V. Kildishev, Naturally hyperbolic, Nature Photonics 9, 214 (2015).
  * Dai  _et al._ [2014] S. Dai, Z. Fei, Q. Ma, A. Rodin, M. Wagner, A. McLeod, M. Liu, W. Gannett, W. Regan, K. Watanabe, _et al._ , Tunable phonon polaritons in atomically thin van der waals crystals of boron nitride, Science 343, 1125 (2014).
  * Mikhailov and Ziegler [2007] S. A. Mikhailov and K. Ziegler, New electromagnetic mode in graphene, Physical review letters 99, 016803 (2007).
  * Chshelokova  _et al._ [2012] A. V. Chshelokova, P. V. Kapitanova, A. N. Poddubny, D. S. Filonov, A. P. Slobozhanyuk, Y. S. Kivshar, and P. A. Belov, Hyperbolic transmission-line metamaterials, Journal of Applied Physics 112, 073116 (2012).
  * Shchelokova  _et al._ [2014] A. V. Shchelokova, D. S. Filonov, P. V. Kapitanova, and P. A. Belov, Magnetic topological transition in transmission line metamaterials, Physical Review B 90, 115155 (2014).
  * Trushkov and Iorsh [2015] I. Trushkov and I. Iorsh, Two-dimensional hyperbolic medium for electrons and photons based on the array of tunnel-coupled graphene nanoribbons, Physical Review B 92, 045305 (2015).
  * Sun  _et al._ [2014] J. Sun, N. M. Litchinitser, and J. Zhou, Indefinite by nature: from ultraviolet to terahertz, Acs Photonics 1, 293 (2014).
  * Zheng  _et al._ [2019] Z. Zheng, N. Xu, S. L. Oscurato, M. Tamagnone, F. Sun, Y. Jiang, Y. Ke, J. Chen, W. Huang, W. L. Wilson, _et al._ , A mid-infrared biaxial hyperbolic van der waals crystal, Science advances 5, eaav8690 (2019).
  * Tao  _et al._ [2019] J. Tao, L. Wu, G. Zheng, and S. Yu, Cherenkov polaritonic radiation in a natural hyperbolic material, Carbon 150, 136 (2019).
  * Dehdast  _et al._ [2021] M. Dehdast, Z. Valiollahi, M. Neek-Amal, B. Van Duppen, F. Peeters, and M. Pourfath, Tunable natural terahertz and mid-infrared hyperbolic plasmons in carbon phosphide, Carbon 178, 625 (2021).
  * Van Veen  _et al._ [2019] E. Van Veen, A. Nemilentsau, A. Kumar, R. Roldán, M. I. Katsnelson, T. Low, and S. Yuan, Tuning two-dimensional hyperbolic plasmons in black phosphorus, Physical Review Applied 12, 014011 (2019).
  * Huang  _et al._ [2018] X. Huang, Y. Cai, X. Feng, W. C. Tan, D. M. N. Hasan, L. Chen, N. Chen, L. Wang, L. Huang, T. J. Duffin, _et al._ , Black phosphorus carbide as a tunable anisotropic plasmonic metasurface, ACS photonics 5, 3116 (2018).
  * Gomez-Diaz and Alu [2016] J. Gomez-Diaz and A. Alu, Flatland optics with hyperbolic metasurfaces, ACS Photonics 3, 2211 (2016).
  * Gric and Hess [2017] T. Gric and O. Hess, Tunable surface waves at the interface separating different graphene-dielectric composite hyperbolic metamaterials, Optics express 25, 11466 (2017).
  * Sreekanth  _et al._ [2013] K. Sreekanth, A. De Luca, and G. Strangi, Negative refraction in graphene-based hyperbolic metamaterials, Applied Physics Letters 103, 023107 (2013).
  * Zhu  _et al._ [2013] B. Zhu, G. Ren, S. Zheng, Z. Lin, and S. Jian, Nanoscale dielectric-graphene-dielectric tunable infrared waveguide with ultrahigh refractive indices, Optics express 21, 17089 (2013).
  * Xiang  _et al._ [2014] Y. Xiang, X. Dai, J. Guo, H. Zhang, S. Wen, and D. Tang, Critical coupling with graphene-based hyperbolic metamaterials, Scientific reports 4, 5483 (2014).
  * Iorsh  _et al._ [2013] I. V. Iorsh, I. S. Mukhin, I. V. Shadrivov, P. A. Belov, and Y. S. Kivshar, Publisher’s note: Hyperbolic metamaterials based on multilayer graphene structures [phys. rev. b 87, 075416 (2013)], Physical Review B 88, 039904 (2013).
  * Othman  _et al._ [2013] M. A. Othman, C. Guclu, and F. Capolino, Graphene–dielectric composite metamaterials: evolution from elliptic to hyperbolic wavevector dispersion and the transverse epsilon-near-zero condition, Journal of Nanophotonics 7, 073089 (2013).
  * Zhao  _et al._ [2017] B. Zhao, B. Guizal, Z. M. Zhang, S. Fan, and M. Antezza, Near-field heat transfer between graphene/hbn multilayers, Physical Review B 95, 245437 (2017).
  * Chang  _et al._ [2016] Y.-C. Chang, C.-H. Liu, C.-H. Liu, S. Zhang, S. R. Marder, E. E. Narimanov, Z. Zhong, and T. B. Norris, Realization of mid-infrared graphene hyperbolic metamaterials, Nature communications 7, 1 (2016).
  * Adler [1962] S. L. Adler, Quantum theory of the dielectric constant in real solids, Physical Review 126, 413 (1962).
  * Nemilentsau  _et al._ [2016] A. Nemilentsau, T. Low, and G. Hanson, Anisotropic 2d materials for tunable hyperbolic plasmonics, Physical review letters 116, 066804 (2016).
  * Enkovaara  _et al._ [2010] J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen, M. Dułak, L. Ferrighi, J. Gavnholt, C. Glinsvad, V. Haikola, H. Hansen, _et al._ , Electronic structure calculations with gpaw: a real-space implementation of the projector augmented-wave method, Journal of physics: Condensed matter 22, 253202 (2010).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for brillouin-zone integrations, Physical review B 13, 5188 (1976).
  * Dion  _et al._ [2004] M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, Van der waals density functional for general geometries, Physical review letters 92, 246401 (2004).
  * Wellendorff  _et al._ [2012] J. Wellendorff, K. T. Lundgaard, A. Møgelhøj, V. Petzold, D. D. Landis, J. K. Nørskov, T. Bligaard, and K. W. Jacobsen, Density functionals for surface science: Exchange-correlation model development with bayesian error estimation, Physical Review B 85, 235149 (2012).
  * Klimeš  _et al._ [2009] J. Klimeš, D. R. Bowler, and A. Michaelides, Chemical accuracy for the van der waals density functional, Journal of Physics: Condensed Matter 22, 022201 (2009).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Physical review letters 77, 3865 (1996).
  * Gritsenko  _et al._ [1995] O. Gritsenko, R. van Leeuwen, E. van Lenthe, and E. J. Baerends, Self-consistent approximation to the kohn-sham exchange potential, Physical Review A 51, 1944 (1995).
  * Wasserman [2005] A. Wasserman, Effective masses, in  _Encyclopedia of Condensed Matter Physics_, edited by F. Bassani, G. L. Liedl, and P. Wyder (Elsevier, Oxford, 2005) pp. 1–5.
  * Frenzel  _et al._ [2017] A. J. Frenzel, C. C. Homes, Q. D. Gibson, Y. Shao, K. W. Post, A. Charnukha, R. J. Cava, and D. N. Basov, Anisotropic electrodynamics of type-ii weyl semimetal candidate wte 2, Physical Review B 95, 245140 (2017).
  * Scholz  _et al._ [2013] A. Scholz, T. Stauber, and J. Schliemann, Plasmons and screening in a monolayer of mos 2, Physical Review B 88, 035135 (2013).
  * Wang  _et al._ [2015] Y. Wang, J. Z. Ou, A. F. Chrimes, B. J. Carey, T. Daeneke, M. M. Alsaif, M. Mortazavi, S. Zhuiykov, N. Medhekar, M. Bhaskaran, _et al._ , Plasmon resonances of highly doped two-dimensional mos2, Nano letters 15, 883 (2015).
  * Hanson [2008] G. W. Hanson, Dyadic green’s functions and guided surface waves for a surface conductivity model of graphene, Journal of Applied Physics 103, 064302 (2008).
  * Low  _et al._ [2014] T. Low, R. Roldán, H. Wang, F. Xia, P. Avouris, L. M. Moreno, and F. Guinea, Plasmons and screening in monolayer and multilayer black phosphorus, Physical review letters 113, 106802 (2014).
  * Gomez-Diaz  _et al._ [2015] J. S. Gomez-Diaz, M. Tymchenko, and A. Alu, Hyperbolic plasmons and topological transitions over uniaxial metasurfaces, Physical review letters 114, 233901 (2015).
  * Kumar  _et al._ [2016] A. Kumar, A. Nemilentsau, K. H. Fung, G. Hanson, N. X. Fang, and T. Low, Chiral plasmon in gapped dirac systems, Physical Review B 93, 041413 (2016).
  * Gao  _et al._ [2019] H. Gao, X. Zhang, W. Li, and M. Zhao, Tunable broadband hyperbolic light dispersion in metal diborides, Optics Express 27, 36911 (2019).
  * Wang  _et al._ [2020] C. Wang, S. Huang, Q. Xing, Y. Xie, C. Song, F. Wang, and H. Yan, Van der waals thin films of wte2 for natural hyperbolic plasmonic surfaces, Nature communications 11, 1 (2020).
  * Wang  _et al._ [2021] F. Wang, C. Wang, A. Chaves, C. Song, G. Zhang, S. Huang, Y. Lei, Q. Xing, L. Mu, Y. Xie, _et al._ , Prediction of hyperbolic exciton-polaritons in monolayer black phosphorus, Nature communications 12, 1 (2021).
  * Schumacher  _et al._ [1998] D. Schumacher, C. Rea, D. Heitmann, and K. Scharnberg, Surface plasmons and sommerfeld–zenneck waves on corrugated surfaces:: Application to high-tc superconductors, Surface science 408, 203 (1998).