# Casimir-Polder interactions of Rydberg atoms
with graphene-based van der Waals heterostructures

K. Wongcharoenbhorn1, C. Koller2, T.M. Fromhold1, W. Li1 1School of Physics
and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK  2University
of Applied Sciences Wiener Neustadt, Johannes Gutenberg-Straße 3, 2700 Wiener
Neustadt, Austria

###### Abstract

We investigate the thermal Casimir-Polder (CP) potential of 87Rb atoms in
Rydberg $n$S-states near single- and double-layer graphene. The dependence of
the CP potential on parameters such as atom-surface distance, temperature,
principal quantum number $n$ and graphene Fermi energy are explored. Through
large scale numerical simulations, we show that, in the non-retarded regime,
the CP potential is dominated by the non-resonant and evanescent-wave terms
which are monotonic, and that, in the retarded regime, the CP potential
exhibits spatial oscillations. We identify that the most important
contributions to the resonant component of the CP potential come from the
$n$S-$n$P and $n$S-$(n-1)$P transitions. Scaling of the CP potential as a
function of the principal quantum number and temperature is obtained. A
heterostructure comprising hexagonal boron nitride layers sandwiched between
two graphene layers is also studied. When the boron nitride layer is
sufficiently thin, the CP potential can be weakened by changing the Fermi
energy of the top graphene layer. Our study provides insights for
understanding and controlling CP potentials experienced by Rydberg atoms near
single and multi-layer graphene-based van der Waals heterostructures.

##  I Introduction

The development of ultracold-atom physics as platforms for chip-based matter
wave manipulation [1], high-accuracy time keeping systems [2, 3], quantum
computing and simulation [4, 5, 6, 7] is an active research field.
Understanding atom-surface interactions is essential for achieving near-
surface atom trapping, as required for the operation of micro-fabricated atom
chips. There is a substantial body of research on trapping ground-state atoms
in metallic-wire-based atom chips [8, 9, 10, 11, 12, 13, 14] as well as on the
interactions of ground-state atoms with various structures of metallic and
perfect conductors [15, 16, 17]. A conclusive review of atom-surface physics
can be found in [18]. It is known that metallic-wire-based atom chips generate
spatially rough trapping potential due to imperfections in the wires [19, 20,
21, 22, 23], high Johnson-noise currents and strong Casimir-Polder (CP)
attractive interactions between the atoms and the chip, causing, for example,
tunnelling losses [24, 25, 26]. In order to enhance the functionality of such
atom chips, different materials are needed for the current-carrying wires.
Recent studies have shown that two-dimensional (2D) materials could offer
desirable properties for overcoming the limitations of using metallic
conductors as current-carrying wires [20, 27].

There is also enduring interest in the technological applications of two-
dimensional (2D) materials, including graphene, for display devices [28],
flexible sensors [29, 30, 31], photo detectors [32, 33, 34, 35, 36], vertical
field-effect transistors [37, 38], and atom chips [27]. New properties and
methods of cooling or patterning 2D materials have been studied [39, 40, 41,
42, 43, 44, 45, 46, 47]. As current-carrying wires in atom chips, graphene has
desirable electronic properties: it has a very low density of electronic
states, high carrier mobility and a linear band structure with zero band gap
in the vicinity of the Dirac points [48], which leads to Johnson noise and CP
attraction far below those typically found for metallic conductors on bulk
substrates [19, 49, 20, 27]. Importantly, smooth trapping potentials can be
obtained using graphene. For example, graphene made from a helium-ion beam
lithography technique has edge roughness of order
$5\text{\,}\mathrm{n}\mathrm{m}$ [21], while the surface roughness of graphene
encapsulated in hexagonal boron nitride (hBN) is on the order of
$12\text{\,}\mathrm{p}\mathrm{m}$ [50]. The coupling of atoms with graphene’s
surface plasmons might be tuneable via changing the Fermi energy of graphene
as the plasmon frequency is proportional to the fourth root of its electronic
density [51, 52, 53, 54, 55]. Previous studies have shown that this could be
used to tailor the CP interactions of atoms trapped near graphene [56, 57, 58,
59, 60, 61, 62, 63, 64, 65, 66].

Figure 1: Schematic diagram of an $b$-layer system, where each layer is
designated by index $q=1,2,...,b$ and is characterized by thickness, $t_{q}$,
permeability, $\mu_{q}$, and permittivity, $\epsilon_{q}$. An atom (red
circle) is located at $\mathbf{r}_{0}=(x_{0},y_{0},z_{0})$ in layer 1 which is
vacuum. Note that, in the Green’s-function calculations, all layers are
infinitely extended in the $x-y$ plane, and that $t_{1}$ and $t_{b}$ are also
infinite.

In recent years, the study of interactions between Rydberg atoms and surfaces
has attracted much attention. Such studies could help us understand the atom-
surface interaction for highly excited atomic states, and open new quantum
technological applications [67], for example by integrating Rydberg atoms with
condensed matter quantum materials like graphene. Rydberg atoms are highly
excited atoms with large principal quantum number, i.e. $n\gg 1$. The size of
Rydberg atoms is proportional to $n^{2}$, and can be a micron when $n\sim 90$,
resulting in weakly bound valence electrons and high electric polarizability,
which in turn causes strong interactions with nearby surfaces [68]. The
lifetime associated with spontaneous emission, i.e. at $0\text{\,}\mathrm{K}$,
is proportional to $n^{3}$. The long lifetime of Rydberg states allows us to
exploit Rydberg atoms for quantum computing and simulation [69, 70]. Atomic
transition frequencies between adjacent Rydberg states are typically in the
microwave-terahertz regions, which are in the same window as thermal energies
at room temperature; Rydberg atoms interact resonantly with thermal photons,
leading to enhanced CP interactions compared to ground-state atoms [71].
Strong coupling between Rydberg atoms and surface plasmon polaritons or
surface phonon polaritons have also been studied both experimentally and
theoretically [72, 73, 74, 75]. There have been previous studies of the CP
interactions of excited two-level or realistic Rydberg atoms near mirrors,
metallic and dielectric surfaces, and metamaterials [76, 77, 78, 79, 80, 81,
82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96], as well as
investigations of the CP interactions of a laser-driven atom and a surface
[97]. The CP interactions between excited molecules or ions and metallic or
dielectric bodies have also been studied [98, 84, 99]. However, CP
interactions between Rydberg atoms and 2D nano-structures such as graphene
have not been fully explored. To the best of our knowledge, the Rydberg atom-
graphene interaction has only been studied in Ref. [100], which focuses on the
zero-temperature limit and shows that graphene can shield the CP force
emanating from a metallic substrate when graphene-substrate distances are
larger than $4\text{\,}\mu$m.

Integrating Rydberg atoms with graphene could lead to promising and novel
quantum devices [101, 102, 103]. In this work, we focus on investigating the
CP interactions of a 87Rb atom in Rydberg $n$S states with graphene. Two
notable features, which are not found in ground-state atoms in the zero-
temperature limit, emerge: stimulated atomic transitions due to thermal
photons give rise to resonant interactions, which can be further distinguished
into attractive and repulsive potentials, and the large sizes of Rydberg atoms
give rise to quadrupole interactions at short atom-surface distances.

The paper is organized as follows. In Sec. II, we provide the quantum field
theoretical description of the CP potential for atoms in Rydberg states and a
brief description of Rydberg atoms. In Sec. III, we first compare the CP
potential of a 87Rb atom in Rydberg states near a suspended single layer of
graphene whose conductivity is described by the Kubo model, and a
$1$-$\mu$m-thick gold sheet whose optical properties are described by the
Drude model. We also consider another model of graphene, with the results
presented, which includes non-local effects such as lattice defects. In Sec.
IV, we present detailed calculations of the CP potential of single-layer
graphene, which allow us to see various characteristics of Rydberg atoms near
graphene that are different from those of ground-state atoms, especially
effects arising from resonant interactions. In Sec. V, we find simple fitted
empirical functions, which capture the scaling of the CP potential with
principal quantum number and the temperature dependence of the interaction.
Finally, in Sec. VI, we extend our study into graphene-based van der Waals
heterostructures comprising two layers of graphene separated by air or
hexagonal boron-nitride (hBN). This allows us to change the spacing between
the graphene layers and change their chemical potentials in order to tune the
interactions between the trapped atoms and the heterostructures. We conclude
in Sec. VII.

##  II Casimir-Polder Potential near Planar Structures

In this section, we will present a theoretical description of the CP energy
shift of highly excited Rydberg atoms near planar layered structures.
Considering only ultracold atoms in Bose-Einstein condensates allows us to
disregard the velocity-dependent CP potential since the atoms are not moving
at a relativistic speed [104]. Further details of the formalism used can be
found in, for example, Refs. [89, 105, 106, 107]. The CP potential arises from
the coupling between an atom and the surrounding body-modified electromagnetic
radiation, described by the coupling Hamiltonian $\hat{H}_{\textrm{AF}}$. The
interaction Hamiltonian in the case of a Rydberg atom can be split into dipole
interactions (first term) and quadrupole interactions (second term):

$\hat{H}_{\textrm{AF}}=-\hat{\textbf{d}}\cdot\hat{\textbf{E}}(\mathbf{r}_{0})-\hat{\textbf{Q}}\bullet\big{[}\bm{\nabla}\otimes\hat{\textbf{E}}(\mathbf{r}_{0})\big{]},$
(1)

where $\hat{\textbf{d}}$ and $\hat{\textbf{Q}}$ are the atomic dipole moment
and quadrupole moment operators, respectively, with
$\hat{\textbf{E}}(\mathbf{r}_{0})$ being the electromagnetic field at the
position of the atom, $\mathbf{r}_{0}$, $\bullet$ denotes the Frobenius inner
product and $\otimes$ tensor product. The energy-level shift up to the second
order, for an atom in state $\ket{u}$ and the body-modified electromagnetic
field in state $\ket{v}$, is [89]

$\delta
E_{u}=\bra{u,v}\hat{H}_{\textrm{AF}}\ket{u,v}+\sum_{u^{\prime},v^{\prime}\neq
u,v}\frac{\absolutevalue{\bra{u,v}\hat{H}_{\textrm{AF}}\ket{u^{\prime},v^{\prime}}}^{2}}{E_{u,v}-E_{u^{\prime},v^{\prime}}},$
(2)

where $E_{u,v}$ are the unperturbed energy eigenvalues of the atom-field
system: $u$ and $v$ are the quantum numbers of the Rydberg state and the
photon field, respectively. This energy shift $\delta E_{u}$ can be expressed
in terms of the Green’s tensor, which can be decomposed into bulk and
scattering parts. It follows that the energy shift can be split into two
components: the position-independent self-energy (associated with the bulk
part of the Green’s tensor) similar to the Lamb shift and the position-
dependent component (associated with the scattering part of the Green’s
tensor), which is, namely, the CP potential.

For an atom, at position $\mathbf{r}_{0}=(x_{0},y_{0},z_{0})$ in the planar
system shown in Fig. 1, in an incoherent superposition of internal-energy
eigenstates $\ket{u}$ (specified by the principal quantum number $n$, the
orbital angular momentum quantum number $l$, the total angular momentum
quantum number $j$, and the $z$-component of the total angular momentum
quantum number $m$) with probabilities $p_{u}$ as described by a density
matrix

$\hat{\sigma}=\sum_{u}p_{u}\ket{u}\bra{u},$ (3)

the total thermal CP potential at environment temperature $T$ can be written
as [105, 107]

$\begin{split}U_{\mathrm{CP}}(\mathbf{r}_{0})&=\sum_{u}p_{u}U_{u}(\mathbf{r}_{0}),\\\
U_{u}(\mathbf{r}_{0})&=U_{u}^{\mathrm{nres}}(\mathbf{r}_{0})+U_{u}^{\mathrm{res}}(\mathbf{r}_{0}).\end{split}$
(4)

Here, $U_{u}(\mathbf{r}_{0})$ is the position-dependent part of $\delta
E_{u}$, $U_{u}^{\mathrm{nres}}(\mathbf{r}_{0})$ is the non-resonant potential
due to virtual photons and $U_{u}^{\mathrm{res}}(\mathbf{r}_{0})$ is the
resonant potential due to real thermal photons.

For the potentials due to dipole interactions, the non-resonant term is given
by [106]

$U_{u}^{\textrm{nres}}(\mathbf{r}_{0})=\mu_{0}k_{B}T\sum_{a=0}^{\infty}{}^{{}^{\prime}}\xi_{j}^{2}\big{[}\bm{\alpha}_{u}(i\xi_{a})\bullet\mathbf{G}^{(\mathrm{s})}(\mathbf{r}_{0},\mathbf{r}_{0},i\xi_{j})\big{]},$
(5)

whereas, the resonant term is given by

$U_{u}^{\textrm{res}}(\mathbf{r}_{0})=-\mu_{0}\sum_{k<u}[N(\omega_{uk})+1]\omega_{uk}^{2}\\\
\times\big{(}\textbf{d}_{uk}\otimes\textbf{d}_{ku}\big{)}\bullet\textrm{Re}\mathbf{G}^{(\mathrm{s})}(\mathbf{r}_{0},\mathbf{r}_{0},\omega_{uk})\\\
+\mu_{0}\sum_{k>u}N(\omega_{ku})\omega_{ku}^{2}\big{(}\textbf{d}_{uk}\otimes\textbf{d}_{ku}\big{)}\bullet\textrm{Re}\mathbf{G}^{(\mathrm{s})}(\mathbf{r}_{0},\mathbf{r}_{0},\omega_{ku}),$
(6)

where $\mu_{0}$ is the permeability of free space, $k_{B}$ the Boltzmann
constant, $\mathbf{G}^{(\mathrm{s})}$ is the scattering Green’s tensor (see
Appendix A for details), and $\textbf{d}_{uk}=\bra{u}\hat{\textbf{d}}\ket{k}$
is the dipole matrix element. The atomic dipole polarizability as a function
of radiation frequency, $\omega$, is defined as

$\bm{\alpha}_{u}(\omega)=\lim_{\epsilon\to 0+}\frac{1}{\hbar}\sum_{k\neq
u}\Bigg{(}\frac{\mathbf{d}_{uk}\otimes\mathbf{d}_{ku}}{\omega_{ku}-\omega-i\epsilon}+\frac{\mathbf{d}_{ku}\otimes\mathbf{d}_{uk}}{\omega_{ku}+\omega+i\epsilon}\Bigg{)},$
(7)

with $\omega_{ku}=(E_{k}-E_{u})/\hbar$ denoting the atomic transition angular
frequencies. The purely imaginary frequencies $\xi_{j}=2\pi
k_{\mathrm{B}}Ta/\hbar$, $a=0,1,2,...$ are the Matsubara frequencies and
$N(\omega)=1/[\mathrm{e}^{\hbar\omega/(k_{\mathrm{B}}T)}-1]$ is the average
thermal photon number in accordance with Bose-Einstein statistics. Following
the separation of the scattering Green’s tensor into propagating-wave and
evanescent-wave components, the resonant term of the CP potential can also be
split into two components as well.

In our calculations, the binding energy of the Rydberg series is given by
$E_{n,l,j}=-\mathrm{Ry}/n^{*2}$ [108], where $\mathrm{Ry}$ is the Rydberg
energy and $n^{*}=n-\delta_{n,l,j}$. Here,
$\delta_{n,l,j}=\delta_{0}+{\delta_{2}/(n-\delta_{2})^{2}}$ is the quantum
defect [109] whose values for 87Rb are tabulated in Table 1.

The electron wavefunction at position $\mathbf{r}$ with respect to the ion
core, $\psi(\mathbf{r})=\psi_{n,l,j,m}(\mathbf{r})$, for the valence electron
is described by the Schrödinger equation

$\Big{[}-\frac{\hbar^{2}}{2\mu}\nabla^{2}-\frac{e^{2}}{4\pi\epsilon_{0}r}\Big{]}\psi(\mathbf{r})=[E_{n,l,j}+V(\mathbf{r})]\psi(\mathbf{r}),$
(8)

where $e$ is the electronic charge, $\mu=Mm_{e}/(M+m_{e})$ is the reduced mass
with $M$ being the mass of the atom, $m_{e}$ is the electronic mass,
$r=\absolutevalue{\mathbf{r}}$ and $V(\mathbf{r})$ is a model potential as
given in Ref. [110], which accounts for the finite size of the core at short
range. Eq. (8) can be separated into angular and radial equations which yield,
in standard notation, spherical harmonics $Y_{l,m_{l}}(\theta,\phi)$ and
radial wavefunctions $R_{n,l,j}(r)$ as solutions, respectively. In this work,
the radial wavefunctions are calculated numerically using the tool provided in
Ref. [110] through Numerov’s method [111].

0State0 |  00000$\delta_{0}$00000 |  00000$\delta_{2}$00000
---|---|---
$n\mathrm{S}_{1/2}$ | 3.1311804 | 0.1784
$n\mathrm{P}_{1/2}$ | 2.6548849 | 0.2900
$n\mathrm{P}_{3/2}$ | 2.6416737 | 0.2950
$n\mathrm{D}_{3/2}$ | 1.34809171 | -0.60286
$n\mathrm{D}_{5/2}$ | 1.34646572 | -0.59600
Table 1: Quantum defects for S, P and D states of Rb atoms [112].

##  III Comparisons between material models

In this section, we will investigate how different conductivity models affect
the CP potential. For graphene, we take the Fermi energy and electron
relaxation rate of graphene to be $E_{F}=$0.1\text{\,}\mathrm{e}\mathrm{V}$$
and $\gamma=$4\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$$, respectively,
corresponding to typical values found both theoretically [113, 114, 115, 116]
and in experiments [51], unless otherwise explicitly stated.

###  III.1 Graphene vs Gold

Figure 2: CP potential energy calculated versus principal quantum numbers for
an 87Rb atom near a $1$-$\mu$m-thick gold sheet (dashed curve) or a suspended
single layer of graphene (solid curve). For each $n$, the atom-surface
separation is taken to be $1/10$ of the wavelengths of the
$n\mathrm{S}$-$(n-1)\mathrm{P}_{1/2}$ transitions. This is to ensure that the
calculations are in the non-retarded regime. For $n<12$, the CP potential for
graphene is weaker than that for gold. For higher $n$, the relative strengths
reverse until $n>20$, when the potential of graphene again becomes stronger
than that of gold.

We first consider a free-standing graphene monolayer, which, from Fig. 1, can
be modeled as a two-layer system, in which the monolayer graphene is located
at the interface between layer 1 and layer 2. Its conductivity is modeled by
the Kubo formula and is not a function of wavevector of impinging radiations
(see Appendix B for the description of the model). Note that throughout this
paper, we will refer to the distance between an atomic core and a surface as
the atom-surface distance. In Fig. 2, we consider the CP potential in the non-
retarded regime by choosing the atom-surface distance to be $1/10$ of the
wavelength of the corresponding $n\mathrm{S}$-$(n-1)\mathrm{P}_{1/2}$
transition for each individual principal quantum number between $6$ and $50$.
The atom-surface distances therefore range from
$136\text{\,}\mathrm{n}\mathrm{m}$ for $n=6$ to
$0.881\text{\,}\mathrm{m}\mathrm{m}$ for $n=50$. For comparison, we also show
the CP potential for an atom near a $1$-$\mu$m-thick free-standing gold sheet
in the same figure. Generally, the CP potential of ground-state atoms near
graphene is weaker than for those near a thick metallic sheet [27]. However,
this is not always the case for Rydberg atoms. We can see that when $n<12$,
the CP potential for graphene is weaker than that for gold. The situation
changes for $n\geq 21$: the CP potential for graphene becomes stronger than
that for gold.

###  III.2 Two models for monolayer graphene conductivity

Another conductivity model is the full non-local model, which includes non-
local effects such as lattice defects. In this model, the conductivity is a
function of both frequency and wavevector of impinging radiations; more
details can be found in Appendix B and in, for example, Ref. [115]. The
comparison of the CP potentials for multiple Rydberg states in the linear-log
scale between the Kubo conductivity and the full non-local conductivity is
shown in panels (a) and (b) of Fig. 3, and the CP potential for the 30S state
is shown in (c). We can see that the potentials from the full non-local model
are slightly weaker than those from the local conductivity model. The two
models give consistent results in terms of magnitudes and spatial dependence.
Due to the excellent agreement between the two models, we will utilize the
Kubo model for shorter computational times in the following sections, unless
otherwise explicitly stated.

Figure 3: CP potentials calculated versus the distance of an 87Rb atom from a
graphene monolayer for two conductivity models of graphene: (a) Kubo
conductivity, (b) full non-local conductivity. (c) shows the CP potential
calculated for the 30S state at $T=$10\text{\,}\mathrm{K}$$ from the full non-
local model (dashed curve) and Kubo model (solid curve).

##  IV Characteristics of the CP potential

In this section, we will explore the dependence of the CP potential on atom-
surface distances, temperature and Rydberg states. Since the resonant part of
the CP potential in Rydberg atoms is enhanced by thermal photons and is
significantly larger than in ground-state atoms, we will also investigate in
detail the contributions from individual atomic transitions. This will provide
insights into the characteristics of the CP potential for Rydberg atoms.

###  IV.1 Characteristic quantities

Limit | Condition | $z$ | $T$
---|---|---|---
Retarded | $z_{0}\omega_{-}/c\gg 1$ | $z_{\omega}\ll z_{0}$ | $T_{z}\ll T_{\omega}$
Non-retarded | $z_{0}\omega_{+}/c\ll 1$ | $z_{0}\ll z_{\omega}$ | $T_{\omega}\ll T_{z}$
Spectroscopic low-T | $k_{B}T\ll\hbar\omega_{-}$ | $z_{\omega}\ll z_{T}$ | $T\ll T_{\omega}$
Spectroscopic high-T | $k_{B}T\gg\hbar\omega_{+}$ | $z_{T}\ll z_{\omega}$ | $T_{\omega}\ll T$
Geometric low-T | $k_{B}T\ll\hbar c/z_{0}$ | $z_{0}\ll z_{T}$ | $T\ll T_{z}$
Geometric high-T | $k_{B}T\gg\hbar c/z_{0}$ | $z_{T}\ll z_{0}$ | $T_{z}\ll T$
Table 2: Three pairs of limits and their associated conditions [106].

There are three pairs of characteristic quantities, of which three are related
to distances and the others to temperatures. The first pair is the geometric
distance $z_{0}$ and the geometric temperature, $T_{z}=\hbar c/(z_{0}k_{B})$,
which is the temperature of radiation whose wavelength is of the order
$z_{0}$. The second pair of parameters is the spectroscopic length,
$z_{\omega}=c/\omega_{\pm}$, which is the measure of the wavelength of the
maximum or minimum of the relevant atomic transition frequencies,
$\omega_{\pm}$, and the spectroscopic temperature
$T_{\omega}=\hbar\omega_{\pm}/k_{B}$. The last pair is the thermal length
$z_{T}=\hbar c/(k_{B}T)$ and the environment temperature $T$. Comparing the
geometric quantities with the spectroscopic quantities allows us to determine
the retarded and non-retarded limits, and comparing the spectroscopic
quantities with the thermal quantities allows us to determine the
spectroscopic high- and low-temperature limits. Finally, comparing the
geometric quantities with the thermal quantities allows us to determine the
geometric high- and low-temperature limits as listed in Table 2. It is
impossible to simultaneously realize the retarded, spectroscopic high-
temperature and geometric low-temperature limit; the same is true for the non-
retarded, spectroscopic low-temperature, and geometric high-temperature limit.

###  IV.2 Scaling Relations

Since the main focus of this work is on Rydberg atoms, we will evaluate the
scaling of the CP potential with respect to the principal quantum number $n$
(without showing graphical results). As shown in Table 3, the non-resonant
term (as shown in Eq. (5)) scales with $n^{7}$ since the polarizability scales
with $n^{7}$. Note that the Green’s tensor as a function of Matsubara
frequencies does not scale with the principal quantum number. The resonant
term has four quantities that scale with $n$ and their combined contribution
leads to an overall $n^{7}$, similar to the non-resonant term. Surprisingly,
the total CP potential, $U^{\mathrm{nres}}_{u}+U^{\mathrm{res}}_{u}$, follows
an $n^{4}$ power law due to their opposite signs.

Quantity | Power
---|---
Polarizability $\bm{\alpha}_{u}$ | $n^{7}$
Atomic transition frequencies $\omega_{uk}$ | $n^{-3}$
Thermal photon number $N(\omega_{uk})$ | $n^{3}$
Dipole moment $\textbf{d}_{uk}$ | $n^{2}$
Green’s tensor $\mathbf{G}^{(\mathrm{s})}(\omega_{uk})$ | $n^{6}$
CP potential due to dipole interactions $U_{\mathrm{CP}}$ | $n^{4}$
Table 3: Scaling of the CP potential and its associated quantities.

###  IV.3 CP variations with atom-surface distances

Figure 4: (a) Linear scale plots of the total CP potential (solid-black
curve), the non-resonant contribution (dashed-blue curve), the resonant
contribution due to evanescent waves (dashed-red curve) and the resonant
contribution due to propagating waves (dashed-green curve) of 87Rb in the 15S
state versus $z_{0}$; (b) Log-log scale plots of the total CP potential
calculated versus atom-surface separation. Note that
$\lambda=$139\text{\,}\mu$$m is roughly the wavelength of the 15S-14P
transition, which is the dominant transition in this case. The parameters are:
$T=$10\text{\,}\mathrm{K}$$, $z_{\omega}=$21.6\text{\,}\mu$$m and
$z_{T}=$229\text{\,}\mu$$m.

Since the CP potential arises from the interaction between atoms and
electromagnetic waves, we can expect that the behaviour of the CP potential in
space can be split into the retarded and non-retarded regimes. Mathematically,
the spatial variation of the CP potential is determined by the Green’s
function. At small atom-surface distances, the interactions with the
evanescent waves will dominate and, at large atom-surface distances, the
interactions with propagating waves will dominate. The Green’s function scales
with a $1/z_{0}^{3}$ power law in the non-retarded regime, and hence the same
scaling is found in the CP potential.

As an example, Fig. 4 shows the variation of the CP potential of the 15S state
calculated versus $z_{0}$ and plotted on (a) a linear scale and (b) a log-log
scale. Note that the wavelength of the 15S-14P transition, which is the
dominant transition in this case, is $\lambda\approx$139\text{\,}\mu$$m, and
that the two characteristic lengths are $z_{\omega}\approx$22\text{\,}\mu$$m
and $z_{T}\approx$229\text{\,}\mu$$m. The spatial variation of the non-
resonant contribution is monotonic, similar to that of ground-state atoms. The
figure reveals the dominance of the non-resonant contribution together with
the resonant contribution due to evanescent waves at the distances much
smaller than $z_{\omega}$ (below $6\text{\,}\mu$m). For higher atom-surface
separations, up to about $z_{\omega}$, the positive contribution from the
resonant term (evanescent wave) dominates. From $z_{0}\sim\lambda/2$ onward,
the resonant term due to the propagating-wave interactions comes into play and
the CP potential oscillates around zero. The log-log plot in panel (b) also
shows a clear distinction of the non-retarded regime, in which
$U_{\mathrm{CP}}\propto 1/z_{0}^{3}$.

Figure 5: Characteristic wavelengths of the CP potential calculated versus $n$
for $T=$10\text{\,}\mathrm{K}$$ (blue curve) and $T=$300\text{\,}\mathrm{K}$$
(green curve) calculated from the first cycle that appears after the atom-
surface distance exceeds $\lambda_{n\mathrm{S}-(n-1)\mathrm{P}}/2$. Also shown
are the half wavelengths of the nearest downward (red dashed curve) and upward
(black dashed curved) atomic transitions calculated versus principal quantum
numbers. At $10$ K, the wavelengths of the oscillation of the CP potential for
$n$ roughly below $30$ leans towards the half wavelengths of the downward
atomic transitions to the $(n-1)\mathrm{P}$ states, but they shift towards the
half wavelengths of the upward atomic transitions to the $n\mathrm{P}$ states
for most the states considered here due to increases in thermal photons (see
discussion in Sec. IV.4).

In the retarded regime, the CP potential exhibits oscillatory behavior
determined by the dominant atomic transition frequencies: the wavelength of
the spatial oscillation of the CP potential, $\lambda_{\mathrm{CP}}$, is
roughly half the wavelength of the dominant transition frequencies,
$\lambda_{n\mathrm{S}-n\mathrm{P}}$ and
$\lambda_{n\mathrm{S}-(n-1)\mathrm{P}}$ (see Fig. 5). The oscillation starts
when the atom-surface distance is approximately half the corresponding atomic
transition wavelength, which can be seen, for example, at the distance
indicated by the purple vertical line in Fig. 4. For low-$n$ states, the
dominant transition frequencies are high, the CP potential starts to oscillate
at short atom-surface distances around $100\text{\,}\mu$m; for high-$n$
states, the dominant transition frequencies are low with corresponding
wavelengths of a few mm.

In order to quantify how each intermediate state $\ket{k}$ in Eq. (6)
contributes to the total resonant CP potential, we consider a relative
contribution, $R^{\mathrm{res}}$, which is defined as the absolute value of
each term (specified by $k$) divided by the summation of the absolute value of
all terms in Eq. (6). The relative contributions of the four adjacent energy-
states for 87Rb atom in the 15S state ($z_{\omega}\approx$22\text{\,}\mu$$m)
at $z_{0}=$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$,$10\text{\,}\mu$$m, and
$69\text{\,}\mu$m together with $R^{\mathrm{res}}$ for the 13P and 14P states
are shown versus atom-surface distances in Fig. 6. Since we know that the non-
resonant potential is monotonic, the underlying mechanics of the oscillations
of the CP potential must originate from the resonant term. The signs of the
individual contributions alter with distances and at some points the positive
and negative contributions cancel out each other, resulting in a negligible
potential. We can also see that $z_{\omega}\ll z_{T}$, indicating that this is
the spectroscopic low-temperature regime. As a result, the contributions from
all upward atomic transitions are negligible due to the lack of thermal
photons.

Figure 6: Relative contributions of the four adjacent energy-states to the
resonant term for an 87Rb atom in the 15S state
($z_{\omega}\approx$22\text{\,}\mu$$m) at $z_{0}=$ (a) $2\text{\,}\mu$m, (b)
$10\text{\,}\mu$m and (c) $69\text{\,}\mu$m; (d) $R^{\mathrm{res}}$ versus
atom-surface distances. Red and blue bars indicate negative and positive
contributions, respectively. We can see that the signs of $R^{\mathrm{res}}$
for each transition change with distance, and at the points governed by
$\lambda_{\mathrm{CP}}$ the positive and negative contributions cancel out
each other, resulting in an overall zero potential as shown in Fig. 4. Note
that, at this temperature, the contribution from the 15S-15P upward-transition
is very small due to the lack of thermal photons. The parameters are
$T=$10\text{\,}\mathrm{K}$$ and $z_{T}=$229\text{\,}\mu$$m.

###  IV.4 Temperature and Rydberg state dependence

Thermal effects on the CP potential of ground-state atoms and molecules depend
on the dominant atomic transition frequencies, which determine the non-
retarded and retarded types: the near-surface (order of $\mu$m) CP potentials
of molecules are likely to be non-retarded and hence insensitive to
temperature as the dominant transition frequencies are low
($\sim$1\text{\times}{10}^{13}\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{s}$$),
while the CP potentials for ground-state atoms are strongly retarded and thus
greatly affected by temperature since the dominant transition frequencies are
high
($\sim$1\text{\times}{10}^{15}\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{s}$$)
[117, 71, 118].

First, we consider how the conductivity of graphene depends on temperature.
There is an explicit linear temperature dependence in the non-resonant
potential–see Eq. (5). The temperature dependence for the resonant term, Eq.
(6), is embedded in the thermal photon distribution function, where
$N(\omega_{uk})\approx k_{B}T/\hbar\omega_{uk}$ as long as
$\hbar\omega_{uk}/k_{B}T\ll 1$, i.e. the spectroscopic high-temperature limit.
On the other hand, there are five quantities that depend on the principal
quantum number as mentioned in Section IV.2.

Figure 7: Variations of the CP potential with (a) principal quantum number,
calculated at fixed distance $z_{0}=$10\text{\,}\mu$$m and
$T=$10\text{\,}\mathrm{K}$$ with $T_{z}=$229\text{\,}\mathrm{K}$$ and (b)
temperature, calculated for the 40S state at $z_{0}=$5\text{\,}\mu$$m
($T_{z}=$458\text{\,}\mathrm{K}$$ and
$T_{\omega}\approx$3\text{\,}\mathrm{K}$$). The CP potential is proportional
to $n^{4}$ and $T$.

Fig. 7 shows the variations of the CP potential with (a) principal quantum
number, calculated at $z_{0}=$10\text{\,}\mu$$m and (b) temperature,
calculated for the 40S state at $z_{0}=$5\text{\,}\mu$$m
($T_{z}=$458\text{\,}\mathrm{K}$$ and
$T_{\omega}\approx$3\text{\,}\mathrm{K}$$). Even though, the non-resonant and
resonant potential, each scales with $n^{7}$, the total CP potential obeys an
$n^{4}$ power law. The CP potential changes linearly with temperature. The
reason is that the dominant transition energies of Rydberg atoms between two
nearest energy-states, $(n-1)$P and $n$P, are small compared to thermal
energy, $k_{B}T$, which makes the condition $\hbar\omega_{uk}/k_{B}T\ll 1$ so
the average photon number is reduced to $N(\omega_{uk})\approx
k_{B}T/\hbar\omega_{uk}$, confirming a linear relation with $T$ as shown in
panel (b).

Figure 8: The relative contributions of the intermediate states to the
resonant parts of the CP potential at $5\text{\,}\mu$m for an 87Rb atom in
(from (a)-(e) respectively) the 25S, 30S, 35S, 40S and 45S states. Red and
blue bars indicate negative and positive contributions, respectively. In this
case, $T_{\omega}\approx$14\text{\,}\mathrm{K}$$ for the 25S state and
$\approx$2\text{\,}\mathrm{K}$$ for the 45S state, and
$T_{z}=$458\text{\,}\mathrm{K}$$, indicating the non-retarded limit. Comparing
the bar heights among the $(n-1)$P states and among the $n$P states, we can
see that the contributions from states with energy levels below the target
states get smaller as $n$ increases and vice versa for the higher energy
states. The temperature is $T=$10\text{\,}\mathrm{K}$$.

We now consider how the resonant contributions change with the principal
quantum number. In Fig. 8, we plot the contributions of the intermediate
states to the resonant parts of the CP potential in the 25S, 30S, 35S, 40S and
45S states for the environment temperature $T=$10\text{\,}\mathrm{K}$$. Note
that the systems are in the non-retarded limit, and the spectroscopic
temperature for the 25S state is approximately $14\text{\,}\mathrm{K}$,
slightly higher than the environment temperature, whereas
$T_{\omega}\approx$2\text{\,}\mathrm{K}$$ for the 45S state, indicating the
spectroscopic high-temperature limit. We can see that the contributions from
the intermediate states $\ket{k}$ to the resonant parts also change when the
principal quantum number of the target state $\ket{u}$ changes. As shown in
Fig. 8, the (downward-transition) contributions from the intermediate states
with energy levels below the target states (below and including $(n-1)$P
state) get smaller as $n$ increases and vice versa for the intermediate states
with higher energies (i.e. upward-transition contributions), which is captured
in panel (f). In other words, the system approaches the spectroscopic high-
temperature limit as $n$ increases, resulting in an increase in the repulsive
resonant potential. This behaviour of the resonant parts will affect the
scaling of the CP potential with $n$ and also enhances the temperature
dependence of the upward-transition contributions as we shall see in more
detail in Sec. V.

We now proceed to look into detail at how the resonant contributions from the
two nearest energy states change with temperature. For the 20S state, there is
a slight increase in the positive contribution from the 20P state as $T$
increases, while the opposite happens to the 19P state, as shown in Fig. 9(a),
resulting in the evanescent-wave potential being more positive as $T$
increases. Fig. 9(b) shows the non-resonant and resonant contributions due to
evanescent and propagating waves. The attractive contribution from the non-
resonant part is stronger than the repulsive resonant potential.

Figure 9: (a) Relative contributions from the 19P and 20P states; (b) the non-
resonant, evanescent-wave, and propagating-wave components of the CP potential
for an atom in the 20S state at $5\text{\,}\mu$m
($T_{\omega}=$35\text{\,}\mathrm{K}$$). Panel (a) shows that the positive
contribution from the 20P state increases, while the negative contribution
from the 19P state decreases as $T$ increases, resulting in a more positive
evanescent-wave potential as $T$ increases, as shown in panel (b).

###  IV.5 Effects of changing Fermi energy

In this section, we will investigate how the CP potential of single-layer
graphene depends on its Fermi energy, which influences the conductivity of the
graphene sheet and, hence, its reflection coefficients and underlying Green’s
tensor. We expect stronger non-resonant contributions to the CP potential as
we increase the Fermi energy since we are adding more conduction electrons.
Regarding the resonant part, the dominant atomic transition frequencies of the
atomic states usually considered are in the microwave-terahertz spectral
regions, which are small compared to the typical Fermi energy. Consequently,
as discussed in Appendix B, we expect that graphene’s optical conductivity
will be dominated by the intraband-process term, which is almost linear in
$E_{F}$ and is a few orders of magnitude higher than the universal AC
conductivity of graphene $\sigma_{0}$, as shown in Fig. 10. However, the CP
potential does not vary linearly with the Fermi energy as we can see in Fig.
11, where we show the CP potential calculated at (a) $z_{0}=$2\text{\,}\mu$$m
and (b) $z_{0}=$10\text{\,}\mu$$m in the 30S state. In (a), at
$T=$50\text{\,}\mathrm{K}$$, the potential peaks at
$E_{F}=$0\text{\,}\mathrm{e}\mathrm{V}$$ then decreases and rises again as the
Fermi energy deviates from $0$ eV. This is not the case for
$T=$300\text{\,}\mathrm{K}$$ and $T=$400\text{\,}\mathrm{K}$$: the potential
then gets more attractive as the doping level increases. In (b), all three
curves exhibit similar behaviors: the potential is the most attractive at
$E_{F}=$0\text{\,}\mathrm{e}\mathrm{V}$$.

Figure 10: The real and imaginary parts of the optical conductivity of
graphene, calculated at
$\omega=$9.88\text{\times}{10}^{11}\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{s}$$,
which is the atomic transition frequency between the 30S and 30P state. The
conductivity increases almost linearly as the Fermi energy increases. The
temperature is $T=$300\text{\,}\mathrm{K}$$. Figure 11: CP potential
calculated versus Fermi energy, $E_{F}$, at atom-surface distances (a)
$z_{0}=$2\text{\,}\mu$$m and (b) $z_{0}=$10\text{\,}\mu$$m for three different
temperatures for a 87Rb atom in the 30S state.

##  V Fitted empirical functions

In this section, we will try to find fitted empirical functions that follow
the power laws discussed in the preceding sections. Let us begin by writing
the CP potential in the form:

$U_{\mathrm{CP}}=-\frac{C_{\alpha}}{z_{0}^{\alpha}},$ (9)

where $C_{\alpha}$ is a positive dispersion coefficient and $\alpha$ is a
scaling power to be evaluated.

###  V.1 Low-temperature limit

We first consider the non-retarded, spectroscopic low-temperature limit in
order to determine the scaling relation of the CP potential that largely
depends on the atom-surface separation and on the principal quantum number.
Figure 12 shows the CP potential at $T=$10\text{\,}\mathrm{K}$$ for (a) a
$1\text{\,}\mu$m-thick gold sheet, (b) graphene modeled by the Kubo
conductivity and (c) graphene modeled by the full non-local conductivity.
Panel (d) shows the dispersion coefficients, $C_{3}$, versus the principal
quantum numbers. From the slopes of the CP potential in panels (a)-(c),
$\alpha$ is determined to be $3$ and it follows from the slopes of $C_{3}$ in
panel (d) that $C_{3}\propto n^{4}$ approximately. To be more precise,
assuming $C_{3}=q_{1}n^{4}+q_{2}n^{3}$ ($q_{1}$ and $q_{2}$ are coefficients),
we find, for the gold sheet,

$\frac{C_{3}}{(\mathrm{Hz\cdot
m^{3}})}=($1.936\text{\times}{10}^{-16}\text{\,}$)n^{4}-($1.893\text{\times}{10}^{-15}\text{\,}$)n^{3},$
(10)

for graphene modeled by the Kubo conductivity,

$\frac{C_{3}}{(\mathrm{Hz\cdot
m^{3}})}=($1.923\text{\times}{10}^{-16}\text{\,}$)n^{4}-($1.840\text{\times}{10}^{-15}\text{\,}$)n^{3},$
(11)

and for graphene modeled by the full non-local model,

$\frac{C_{3}}{(\mathrm{Hz\cdot
m^{3}})}=($1.924\text{\times}{10}^{-16}\text{\,}$)n^{4}-($1.848\text{\times}{10}^{-15}\text{\,}$)n^{3}.$
(12)

Alternatively, if we assume that $C_{3}$ is proportional to a single power of
$n$, we find, for the gold sheet,

$\frac{C_{3}}{(\mathrm{Hz\cdot
m^{3}})}=($3.394\text{\times}{10}^{-17}\text{\,}$)n^{4.397},$ (13)

for graphene modeled by the Kubo conductivity,

$\frac{C_{3}}{(\mathrm{Hz\cdot
m^{3}})}=($3.543\text{\times}{10}^{-17}\text{\,}$)n^{4.385},$ (14)

and for graphene modeled by the full non-local model,

$\frac{C_{3}}{(\mathrm{Hz\cdot
m^{3}})}=($3.517\text{\times}{10}^{-17}\text{\,}$)n^{4.387}.$ (15)

We can see from the above equations that both the Kubo and full non-local
models give quantitatively similar results.

Figure 12: Panels (a)-(c) show the CP potential calculated, versus atom-
surface separation, and plotted on a log-log scale for a $1$-$\mu$m-thick gold
sheet, graphene modeled by the Kubo conductivity and graphene modeled by the
full non-local conductivity, respectively. Panel (d) shows the relations
between the dispersion coefficients, $C_{3}$ (obtained from Eq. (9)), and the
principal quantum numbers of the atomic states considered, calculated at
$z_{0}=$10\text{\,}\mu$$m, with two insets zoomed near $n=30$ and $n=40$.

###  V.2 Including temperature-dependent effects

Figure 13: (a, b) The CP potential and the dispersion coefficient calculated
versus temperature for a graphene-monolayer modeled by the Kubo conductivity;
(c, d) the polynomial coefficients of the fitted empirical formula. The
spacing of $C_{3}$ follows an $n^{4}$ power law, while its slope follows an
$n^{7}$ power law. The parameter is $z_{0}=$5\text{\,}\mu$$m.

In this subsection, we try to include the scaling law of the CP potential with
temperature in an attempt to find a more complete fitted empirical formula
which can describe the CP potential of a 87Rb atom near a single graphene
sheet in the non-retarded regime and for various $n$. This is done by
calculating the CP potential versus temperature to obtain the dispersion
coefficient, $C_{3}$, as a function of temperature, $T$, and principal quantum
number, $n$, in the form of a linear equation $C_{3}(n,T)=p_{1}(n)T+p_{2}(n)$.
Figs. 13(a) and (b) show the CP potential near graphene vs $T$ for the 20S,
25S, 30S, 35S and 40S states and its associated dispersion coefficient,
$C_{3}$. The $U_{\mathrm{CP}}(T)$ curves show a linear relationship between
the potential and temperature. However, the slopes of the $U_{\mathrm{CP}}(T)$
and $C_{3}(T)$ curves become steeper as $n$ increases since the positive
resonant contributions increase as $n$ increases, as shown in Fig. 8. The non-
linear relations between the coefficients, $p_{1}$ and $p_{2}$, and $n$ are
shown in Figs. 13(c) and (d). We can fit them with polynomial equations of
degree $7$ and $4$ as follows (see Sec. IV.2):

$\frac{p_{1}(n)}{(\mathrm{Hz\cdot
m^{3}}/K)}=\Big{[}-$4\text{\times}{10}^{-25}\text{\,}$n^{7}-$9.38\text{\times}{10}^{-15}\text{\,}$\Big{]},$
(16)

$\frac{p_{2}(n)}{(\mathrm{Hz\cdot
m^{3}})}=\Big{[}$1.866\text{\times}{10}^{-16}\text{\,}$n^{4}-$1.614\text{\times}{10}^{-15}\text{\,}$n^{3}\Big{]}.$
(17)

We may write the fitted empirical function for the CP potential for graphene-
monolayer in the spectroscopic high-temperature and non-retarded regime as

$U_{\mathrm{CP}}(n,T,z_{0})=-\frac{p_{1}(n)T+p_{2}(n)}{z_{0}^{3}}.$ (18)

###  V.3 Accuracy of the fitted empirical functions

We now check the accuracy of our fitted empirical formula, Eq.(18), compared
with the general formula, Eq. (4), by plotting the CP potential for atom-
surface distances in the range $1\text{\,}\mu$m to $10\text{\,}\mu$m and
calculating the absolute values of the energy differences. Fig. 14 shows the
CP potential calculated by the general formula (blue solid curves) and the
empirical formula (orange dashed curves), together with their corresponding
energy differences, $\absolutevalue{\Delta U}$, at (a, c)
$T=$10\text{\,}\mathrm{K}$$ and (b, d) $T=$300\text{\,}\mathrm{K}$$. We can
see that the empirical formula gives results that are off by
$\sim$1\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ at $1\text{\,}\mu$m. Beyond
this point, the differences are of order smaller than Megahertz. The empirical
formula provides the least accurate results for the 20S state with a relative
error as high as $10\%$. As for the other two states, the relative error is
far below $1\%$.

Figure 14: The CP potential calculated by the general formula (blue solid
curves) and by the fitted empirical function Eq.(18) (orange dashed curves) at
(a) $T=$10\text{\,}\mathrm{K}$$ and (b) $T=$300\text{\,}\mathrm{K}$$ and their
corresponding differences (c, d). The fitted function provides the least
accurate results for the 20S state with a relative error as high as $10\%$. As
for the other two states, the relative error is far below $1\%$. The parameter
is $z_{0}=$5\text{\,}\mu$$m.

##  VI Heterostructures containing two graphene layers

In this section, we study the CP potential of heterostructures containing two
separated layers of graphene, which provide a wider range of tuneable optical
properties than single-layer graphene (SLG) [119, 120, 121]. Hexagonal boron
nitride (hBN) is commonly used to integrate with graphene to form complex
structures; we can use it, for example, as a substrate, a tunnel barrier
[122], or an encapsulating layer [50, 123]. We can grow hBN on graphene
vertically [124] or laterally [125].

Let us consider a simple van der Waals (vdW) heterostructure comprising an hBN
layer with permittivity $\epsilon_{\mathrm{hBN}}=3.58$ [126] sandwiched
between two graphene monolayers with spacing $d$ between them. Bare graphene
layers (i.e. without an hBN spacer) in vacuum will also be considered as a
comparison. Figure 15 shows the CP potential of the said structures at
$T=$300\text{\,}\mathrm{K}$$ and $z_{0}=$2\text{\,}\mu$$m for a 87Rb atom in
the 30S state. In (a), the spacing $d$ is varied between $1$ nm and
$1\text{\,}\mu$m; the potential becomes more negative when we increase the
spacing between two graphene layers until
$d\approx$11\text{\,}\mathrm{n}\mathrm{m}$$ for the graphene-vacuum-graphene
structure and $d\approx$73\text{\,}\mathrm{n}\mathrm{m}$$ for the graphene-
hBN-graphene structure, when the potential gets less negative again. In (b),
we show that the potential of the structure containing two layers of graphene
approaches that of the single-layered one as $d$ is increased to
$1\text{\,}\mu$m.

Figure 15: CP potential of a 87Rb atom in the 30S state (a) versus the spacing
$d$ of heterostructures containing two graphene layers, separated by vacuum
and by hBN; (b) versus $z_{0}$ for a single-layer graphene (SLG, black dashed
curve) and a graphene-vacuum-graphene structure with spacing
$d=$11\text{\,}\mathrm{n}\mathrm{m}$$ and $d=$1\text{\,}\mu\mathrm{m}$$. The
parameters are: T = $300\text{\,}\mathrm{K}$,
$E_{F}=$0.1\text{\,}\mathrm{e}\mathrm{V}$$, $z_{0}=$2\text{\,}\mu\mathrm{m}$$
and $\epsilon_{\mathrm{hBN}}=3.58$.

In order to show the relationship between $z_{0}$ and $d$, in Fig. 16 we plot
the color map of the differences (in kHz) between the CP potential of SLG and
that of graphene-vacuum-graphene structure on a linear-log scale. We can see
once again that the behavior of the two-graphene-layered structure approaches
that of SLG as the spacing $d$ increases and that the atom effectively
experiences the two-graphene-layered structure as SLG when the atom-surface
distance is large enough.

Figure 16: Color map of the energy difference (in kHz) between the CP
potential of single-layer graphene and graphene-vacuum-graphene structure with
spacing $d$ calculated versus $z_{0}$ and $d$ (on a linear-log scale) for a
87Rb atom in the 30S state. The plot shows that the behavior of the two layers
of graphene approaches that of SLG as the spacing $d$ increases. The
parameters are: T = $300\text{\,}\mathrm{K}$,
$E_{F}=$0.1\text{\,}\mathrm{e}\mathrm{V}$$.

Now let us consider how the CP potential changes when we simultaneously vary
the Fermi levels of the top and bottom graphene layers. As a background, for a
ground-state atom near two layers of graphene, the CP potential generally
becomes more attractive as graphene sheets are electrically doped [100].
Figure 17 shows the CP potential (in MHz) versus the Fermi energies $E_{F}$ of
the top and bottom layer of the graphene-hBN-graphene structure for a 87Rb
atom in the 30S state at T = $300\text{\,}\mathrm{K}$, taking
$z_{0}=$2\text{\,}\mu$$m and $d=$10\text{\,}\mathrm{n}\mathrm{m}$$. We can see
that in the middle of the colour map where the Fermi energies of both layers
are zero, the potential is still stronger than in most regions of the
parameter plane. Moreover, it is very surprising that the potential is the
most attractive when $E_{F}$ of the top layer is zero. This cannot simply be
explained by considering the conductivity alone; in addition the non-resonant
and resonant terms are enhanced differently when the Fermi energies are
varied.

Lastly, let us investigate how the CP potential depends on the spacing $d$ and
on the Fermi energy of the top layer when the Fermi energy of the bottom layer
is fixed. Figure 18 shows the CP potential (in MHz) of the graphene-vacuum-
graphene structure calculated versus the spacing $d$ between the graphene
layers (in a logarithmic scale) and the Fermi energy $E_{F}$ of the top layer
when the Fermi energy of the bottom layer is set to
$0\text{\,}\mathrm{e}\mathrm{V}$ at $T$ = $300\text{\,}\mathrm{K}$ and
$z_{0}=$2\text{\,}\mu$$m. We may divide the behavior of the CP potential
landscape into three regimes: the first is the SLG-like regime, in which the
separation between the two graphene layers is large. In this regime, the
potential becomes more negative when the Fermi energy is increased. The second
regime is when $d\sim 400$-$700\text{\,}\mathrm{n}\mathrm{m}$; the potential
then barely changes with Fermi energy. The third regime is when
$d<$400\text{\,}\mathrm{n}\mathrm{m}$$. In this regime, the potential is the
most attractive when $E_{F}=$0\text{\,}\mathrm{e}\mathrm{V}$$.

Figure 17: The CP potential (in MHz) calculated versus the Fermi energies
$E_{F}$ of the top and bottom layer for a 87Rb atom in the 30S state. Here, we
assume that an hBN slab is sandwiched between graphene sheets. Surprisingly,
the potential is not the weakest when both layers are undoped. The parameters
are: T = $300\text{\,}\mathrm{K}$, $z_{0}=$2\text{\,}\mu$$m and
$d=$10\text{\,}\mathrm{n}\mathrm{m}$$. Figure 18: The CP potential (in MHz)
calculated (and plotted on a log-linear scale ) versus the spacing $d$ between
graphene layers and the Fermi energy $E_{F}$ of the top layer for a 87Rb atom
in the 30S state when the Fermi energy of the bottom layer
$=$0\text{\,}\mathrm{e}\mathrm{V}$$. For $d$ below approximately
$400\text{\,}\mathrm{n}\mathrm{m}$, the potential is the most attractive when
$E_{F}=0$ then gets weaker as $E_{F}$ moves away from $0$ symmetrically.
Thereafter, the behavior starts to become similar to that of a single layer of
graphene. The parameters are: T = $300\text{\,}\mathrm{K}$,
$z_{0}=$2\text{\,}\mu$$m.

##  VII Conclusions

In summary, we have calculated and analyzed the CP potential of a rubidium
Rydberg atom positioned near single-layer and double-layer graphene. These
calculations used a Green’s-function method in the framework of macroscopic
quantum electrodynamics. The optical conductivity of the layer(s) was modeled
by the local Kubo equation, including non-local effects. Together, the atomic
electric polarizability, the dipole matrix elements of the atom, and the
electromagnetic reflection coefficients of the surface layer(s), determine the
CP potential. Since Rydberg atoms have higher electric polarizabilities than,
and their dipole matrix elements between adjacent atomic states exceed, those
of ground-state atoms, their CP potential interaction with graphene-based
multilayers is complex.

We have shown that, at $T=$300\text{\,}\mathrm{K}$$ in the non-retarded
regime, the CP potential of monolayer graphene is weaker than that for a
${1}$-$\mu$m-thick gold sheet for low-$n$ states, but the CP potential of
graphene becomes increasingly attractive when $n>20$. Regarding the different
models of graphene’s conductivity, the near-surface potential calculated using
the local conductivity is slightly more attractive than that calculated from
the full non-local conductivity models for typical doping levels
($E_{F}=$0.1\text{\,}\mathrm{e}\mathrm{V}$$) and $T=$10\text{\,}\mathrm{K}$$.
In general, in the non-retarded limit, the CP potential is determined by both
the monotonic attractive non-resonant potential and the evanescent-wave
resonant potential. In the retarded limit, the CP potential is dominated by
the resonant potential alone and spatially oscillates with a periodicity that
approximately equals the half wavelength of the nearest downward transition at
low temperature and the nearest upward transition at high temperature. The
spatial oscillations in the CP potential start to occur when the atom-surface
separation begins to exceed these wavelengths. Thermal effects come into play
when the thermal energy is resonant with the atomic transition energies. In
the non-retarded, spectroscopic low-temperature limit, the atomic transitions
are dominated by downward transitions, which give rise to an attractive
resonant potential. In contrast, in the spectroscopic high-temperature limit,
there are enough thermal photons to stimulate upward transitions, resulting in
a repulsive potential. The spectroscopic high-temperature limit can be easily
realized by increasing the principal quantum number even below room
temperature. Doping graphene generally results in a more attractive CP
potential at high temperatures.

Finally, we have investigated heterostructures containing two graphene sheets
with varying inter-layer separation and Fermi energies. We found that the
effects of changing the spacing and Fermi energies are inter-related. When the
spacing is small ($d\ll$1\text{\,}\mu$$m), doping the top layer either
positively or negatively weakens the CP potential. By contrast, when the
spacing is large the behavior approaches that of single-layer graphene.

Possible future work could be done on multiple-layer structures comprising
graphene and other 2D materials. The complex interference of electromagnetic
waves within and near such structures could greatly affect the resonant CP
potential. Studying the many-body effects of Rydberg atoms near graphene-based
heterostructures might also be interesting since parameters such as the layer
chemical potentials, inter-layer spacing, and/or the number and type of 2D
layers can be altered to manipulate the interaction with, and behavior of,
nearby trapped atoms.

Acknowledgements: This work is supported by the EPSRC through Grant No.
EP/R04340X/1 via the QuantERA project “ERyQSenS”.

##  Appendix A Green’s tensor

An electric field created by a radiating electric dipole can be described by a
classical Green’s tensor whose form for planar multilayered systems is well-
known. When analysing the resonant part of the CP potential, it is useful to
split the equal-position scattering Green’s tensor into evanescent-wave and
propagating-wave components as follows [106, 116]:

$\mathbf{G}^{(\mathrm{s})}(\mathbf{r}_{0},\mathbf{r}_{0},\omega)=\mathbf{G}^{(\mathrm{s})}_{\mathrm{evan}}(\mathbf{r}_{0},\mathbf{r}_{0},\omega)+\mathbf{G}^{(\mathrm{s})}_{\mathrm{prop}}(\mathbf{r}_{0},\mathbf{r}_{0},\omega).$
(19)

The evanescent-wave component takes the form

$\mathbf{G}^{(\mathrm{s})}_{\mathrm{evan}}(\mathbf{r}_{0},\mathbf{r}_{0},\omega)=\frac{1}{8\pi}\int_{0}^{\infty}\mathrm{d}\kappa^{\perp}\mathrm{e}^{-2\kappa^{\perp}z_{0}}\\\
\times\Big{[}\mathbf{M}_{\alpha}r_{s}^{(\mathrm{s})}(k^{\parallel},\omega)+\frac{c^{2}}{\omega^{2}}\mathbf{M}_{\beta}r_{p}^{(\mathrm{s})}(k^{\parallel},\omega)\Big{]},$
(20)

and the propagating-wave component takes the form

$\mathbf{G}^{(\mathrm{s})}_{\mathrm{prop}}(\mathbf{r}_{0},\mathbf{r}_{0},\omega)=\frac{i}{8\pi}\int_{0}^{\omega/c}\mathrm{d}k^{\perp}\mathrm{e}^{2ik^{\perp}z_{0}}\\\
\times\Big{[}\mathbf{M}_{\alpha}r_{s}^{(\mathrm{s})}(k^{\parallel},\omega)+\frac{c^{2}}{\omega^{2}}\mathbf{M}_{\beta}r_{p}^{(\mathrm{s})}(k^{\parallel},\omega)\Big{]},$
(21)

where $z_{0}$ is the shortest distance between the surface and the center of
the atom, $r_{s}^{(\mathrm{s})}$ and $r_{p}^{(\mathrm{s})}$ are the Fresnel
reflection coefficients for the $s$\- and $p$-polarized waves, respectively,
and

$\kappa^{\perp}=\Big{(}k^{\parallel
2}-\mu_{1}\epsilon_{1}\frac{\omega^{2}}{c^{2}}\Big{)}^{1/2},\\\
k^{\perp}=\Big{(}\mu_{1}\epsilon_{1}\frac{\omega^{2}}{c^{2}}-k^{\parallel
2}\Big{)}^{1/2}$ (22)

with $k^{\parallel 2}=k_{x}^{2}+k_{y}^{2}$.
The tensors $\mathbf{M}_{\alpha}$ and $\mathbf{M}_{\beta}$ in Eqs. (20) and
(21) are given by

$\displaystyle\mathbf{M}_{\alpha}=$ $\displaystyle\begin{pmatrix}1&0&0\\\
0&1&0\\\ 0&0&0\end{pmatrix},$ (23) $\displaystyle\mathbf{M}_{\beta}=$
$\displaystyle\begin{cases}\begin{pmatrix}\kappa^{\perp 2}&0&0\\\
0&\kappa^{\perp 2}&0\\\ 0&0&2k^{\parallel 2}\end{pmatrix}\text{for
$\mathbf{G}^{(\mathrm{s})}_{\mathrm{evan}}$},\\\ \\\ \begin{pmatrix}-k^{\perp
2}&0&0\\\ 0&-k^{\perp 2}&0\\\ 0&0&2k^{\parallel 2}\end{pmatrix}\text{for
$\mathbf{G}^{(\mathrm{s})}_{\mathrm{prop}}$}.\end{cases}$ (24)

##  Appendix B Graphene’s Optical Properties

###  B.1 Conductivity models

In this section, we present two models for the optical conductivity of
monolayer graphene: (i) based on local Kubo conductivity,
$\sigma^{\mathrm{Kubo}}(\omega)$, which depends only on the angular
frequencies, $\omega$, of the incident electromagnetic radiations and ignores
spatial dispersion in the graphene surface [115, 127, 128] and (ii) using the
full non-local conductivity, $\sigma^{\mathrm{fnl}}(q,\omega)$, which is
derived from the Lindhard polarisation function in random-phase (RPA) and
relaxation-time (RT) approximations [129, 130, 131, 115]. The latter takes
into account spatial dispersion by non-local causes when considering the
interactions of incident photons, surface plasmons with wavenumbers, $q$, and
graphene’s electrons.

The Kubo conductivity can be expressed as the sum of two contributions:
$\sigma_{\mathrm{intra}}(\omega)$, which arises from intraband transition
processes of the electrons and $\sigma_{\mathrm{inter}}(\omega)$, which
describes transitions between the conduction and valence bands; both terms may
be written as follows:

$\sigma_{\mathrm{intra}}(\omega)=\frac{\sigma_{0}}{\pi}\frac{4}{\hbar\gamma-i\hbar\omega}\bigg{[}E_{F}+2k_{B}T\ln\bigg(1+e^{-E_{F}/k_{B}T}\bigg{missing})\bigg{]},$
(25)

$\sigma_{\mathrm{inter}}(\omega)=\sigma_{0}\bigg{[}G(\hbar\omega/2)+i\frac{4\hbar\omega}{\pi}\int_{0}^{\infty}dE\frac{G(E)-G(\hbar\omega/2)}{(\hbar\omega)^{2}-4E^{2}}\bigg{]},$
(26)

in which

$G(X)=\frac{\sinh\bigg{(}\frac{\displaystyle X}{\displaystyle
k_{B}T}\bigg{)}}{\cosh\bigg(\frac{\displaystyle E_{F}}{\displaystyle
k_{B}T}\bigg{missing})+\cosh\bigg(\frac{\displaystyle X}{\displaystyle
k_{B}T}\bigg{missing})},$ (27)

where $\sigma_{0}=e^{2}/(4\hbar)$ is the universal alternating-current
conductivity of graphene, $\gamma$ is the electron relaxation rate in
graphene, $E_{F}$ is the Fermi energy and $T$ is the temperature of the
graphene layer.

As for the other model–the full non-local conductivity, we start by providing
the 2D polarizability in the RPA, $P(q,\omega)=P_{re}(x,y)+iP_{im}(x,y)$, in
the $T=$0\text{\,}\mathrm{K}$$ limit, where $P_{re}(x,y)$ and $P_{im}(x,y)$
are the real and imaginary parts, respectively, in six regions of the
$(q,\hbar\omega)$-space as depicted in Fig. 19 as follows:

Figure 19: Schematic diagram showing six regions in the
$(q,\hbar\omega)$-plane of graphene surface-plasmons (GSP), which are
distinguished by the coupling between the plasmons and excitations of
electron-hole pairs in graphene.

$\bullet$ region 1A

$\displaystyle P_{re}(x,y)=$ $\displaystyle-2t_{1},$ (28) $\displaystyle
P_{im}(x,y)=$
$\displaystyle\phantom{i}\frac{1}{4}t_{1}t_{3}[C_{h}(t_{5})-C_{h}(t_{4})],$
(29)

$\bullet$ region 2A

$\displaystyle P_{re}(x,y)=$
$\displaystyle-2t_{1}+\frac{1}{4}t_{1}t_{3}C(t_{5}),$ (30) $\displaystyle
P_{im}(x,y)=$ $\displaystyle-\frac{1}{4}t_{1}t_{3}C_{h}(t_{4}),$ (31)

$\bullet$ region 3A

$\displaystyle P_{re}(x,y)=$
$\displaystyle-2t_{1}+\frac{1}{4}t_{1}t_{3}[C(t_{4})-C(t_{6})],$ (32)
$\displaystyle P_{im}(x,y)=$ $\displaystyle\phantom{i}0,$ (33)

$\bullet$ region 1B

$\displaystyle P_{re}(x,y)=$
$\displaystyle-2t_{1}+\frac{1}{4}t_{1}t_{2}[C_{h}(t_{4})-C_{h}(t_{5})],$ (34)
$\displaystyle P_{im}(x,y)=$ $\displaystyle\phantom{i}0,$ (35)

$\bullet$ region 2B

$\displaystyle P_{re}(x,y)=$
$\displaystyle-2t_{1}+\frac{1}{4}t_{1}t_{2}C_{h}(t_{4}),$ (36) $\displaystyle
P_{im}(x,y)=$ $\displaystyle\phantom{i}\frac{1}{4}t_{1}t_{2}C(t_{5}),$ (37)

$\bullet$ region 3B

$\displaystyle P_{re}(x,y)=$
$\displaystyle-2t_{1}+\frac{1}{4}t_{1}t_{2}[C_{h}(t_{4})-C_{h}(t_{6})],$ (38)
$\displaystyle P_{im}(x,y)=$ $\displaystyle-\frac{\pi}{4}t_{1}t_{2},$ (39)

where $x=q/k_{F}$, $y=\hbar\omega/E_{F}$, $t_{1}={k_{F}}/{\pi\hbar v_{F}}$,
$t_{2}=x^{2}/\sqrt{y^{2}-x^{2}}$, $t_{3}=x^{2}/\sqrt{x^{2}-y^{2}}$,
$t_{4}=(2+y)/x$, $t_{5}=(2-y)/x$, $t_{6}=(y-2)/x$. Here $v_{F}$ and $k_{F}$
are Fermi velocity and wavenumbers, respectively. The auxiliary functions are
defined as follows

$\displaystyle C_{h}(a)=$ $\displaystyle a\sqrt{a^{2}-1}-\mathrm{arccosh}(a),$
(40) $\displaystyle C(a)=$ $\displaystyle a\sqrt{1-a^{2}}-\mathrm{arccos}(a).$
(41)

The 2D polarizability described above only takes into account intrinsic
mechanisms for the decay of GSP into electron-hole pairs. To include extrinsic
processes such as collisions with lattice defects or impurity scattering, we
extend the previous model by also including scattering events within the
relaxation-time approximation, which allows us to express the 2D
polarizability in the RPA-RT approximation as follows:

$P_{\gamma}(q,\omega)=\frac{(1+i\gamma/\omega)P(q,\omega+\mathrm{i\gamma})}{(1+i\gamma/\omega)\cdot
P(q,\omega+\mathrm{i\gamma})/P(q,0)}.$ (42)

The RPA-RT dielectric function can be written in terms of the 2D polarisation
function as

$\epsilon^{\mathrm{RPA-RT}}(q,\omega)=\epsilon_{r}-v_{q}P_{\gamma}(q,\omega),$
(43)

where $\epsilon_{r}$ is the relative permittivity of the medium in which the
graphene layer is embedded and $v_{q}=e^{2}/2\epsilon_{0}q$ is the Fourier
transform of the Coulomb interaction. Additionally, the longitudinal
conductivity can also be written in terms of the Lindhard polarizability:

$\sigma^{\mathrm{fnl}}(q,\omega)=ie^{2}\frac{\omega}{q}P_{\gamma}(q,\omega).$
(44)

Note that it is this quantity–the conductivity–that will be inserted into our
transfer-matrix calculations of the reflection coefficients of graphene, when
we calculate the CP potential.

Figure 20: Real (a) and imaginary (b) parts of the 2D polarizability of
graphene, $P_{\gamma}(q,\omega)\times h$, plotted in four regions–1A, 2A, 1B
and 2B–of the $(q,\hbar\omega)$ plane. The parameters are
$\gamma=$4\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$$,
$E_{F}=$0.1\text{\,}\mathrm{e}\mathrm{V}$$.

Fig. 20 shows the real and imaginary parts of the 2D polarizability of
graphene, calculated using Eq. (42). Note that the imaginary part in these
regions is always negative.

##  Appendix C Permittivity of metals

The permittivity of gold at radiation frequency $\omega$ is described by the
Drude model [132]

$\epsilon_{\mathrm{metal}}(\omega)=1-\frac{\omega_{\mathrm{p}}^{2}}{\omega^{2}+i\Gamma_{D}\omega},$
(45)

where, for gold,
$\omega_{\mathrm{p}}=$1.35\text{\times}{10}^{16}\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{s}$$
is the plasma frequency and
$\Gamma_{D}=$17.13\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$$ is the electron
scattering rate [133].

## References

  * Chen and Fan [2020] X. Chen and B. Fan, Reports on Progress in Physics 83, 076401 (2020).
  * Liu  _et al._ [2018] L. Liu, D.-S. Lü, W.-B. Chen, T. Li, Q.-Z. Qu, B. Wang, L. Li, W. Ren, Z.-R. Dong, and J.-B. e. a. Zhao, Nature Communications 9, 2760 (2018).
  * Ren  _et al._ [2020] W. Ren, T. Li, Q. Qu, B. Wang, L. Li, D. Lü, W. Chen, and L. Liu, National Science Review 7, 1828 (2020).
  * Knight and Walmsley [2019] P. Knight and I. Walmsley, Quantum Sci. Technol. 4, 040502 (2019).
  * Acín  _et al._ [2018] A. Acín, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, D. Esteve, N. Gisin, S. Glaser, F. Jelezko, S. Kuhr, M. Lewenstein, M. Riedel, P. Schmidt, R. Thew, A. Wallraff, I. Walmsley, and F. Wilhelm, New J. Phys. 20, 080201 (2018).
  * Meinert  _et al._ [2020] F. Meinert, C. Hölzl, M. A. Nebioglu, A. D’Arnese, P. Karl, M. Dressel, and M. Scheffler, Phys. Rev. Research 2, 023192 (2020).
  * Geier  _et al._ [2021] S. Geier, N. Thaicharoen, C. Hainaut, T. Franz, A. Salzinger, A. Tebben, D. Grimshandl, G. Zürn, and M. Weidemüller, Science 374, 1149 (2021).
  * Haug  _et al._ [2019] T. Haug, S. Y. Buhmann, and R. Bennett, Phys. Rev. A 99, 012508 (2019).
  * Scheel  _et al._ [2005] S. Scheel, P. K. Rekdal, P. L. Knight, and E. A. Hinds, Phys. Rev. A 72, 042901 (2005).
  * Rekdal  _et al._ [2004] P. K. Rekdal, S. Scheel, P. L. Knight, and E. A. Hinds, Phys. Rev. A 70, 013811 (2004).
  * Folman  _et al._ [2002] R. Folman, P. Krüger, J. Schmiedmayer, J. Denschlag, and C. Henkel, Advances In Atomic, Molecular, and Optical Physics , 263 (2002).
  * Henkel  _et al._ [1999] C. Henkel, S. Pötting, and M. Wilkens, Applied Physics B 69, 379 (1999).
  * Hinds and Hughes [1999] E. A. Hinds and I. G. Hughes, Journal of Physics D: Applied Physics 32, R119 (1999).
  * Jones  _et al._ [2003] M. P. A. Jones, C. J. Vale, D. Sahagun, B. V. Hall, and E. A. Hinds, Phys. Rev. Lett. 91, 080401 (2003).
  * Messina and Passante [2007] R. Messina and R. Passante, Phys. Rev. A 76, 032107 (2007).
  * Vasile and Passante [2008] R. Vasile and R. Passante, Phys. Rev. A 78, 032108 (2008).
  * Ferreri  _et al._ [2019] A. Ferreri, M. Domina, L. Rizzuto, and R. Passante, Symmetry 11 (2019).
  * Laliotis  _et al._ [2021a] A. Laliotis, B.-S. Lu, M. Ducloy, and D. Wilkowski, AVS Quantum Science 3, 043501 (2021a).
  * Sinuco-León  _et al._ [2011] G. Sinuco-León, B. Kaczmarek, P. Krüger, and T. M. Fromhold, Phys. Rev. A 83, 021401 (2011).
  * Sinuco-León  _et al._ [2018] G. Sinuco-León, P. Krüger, and T. Fromhold, Journal of Modern Optics 65, 677 (2018).
  * Aigner  _et al._ [2008] S. Aigner, L. D. Pietra, Y. Japha, O. Entin-Wohlman, T. David, R. Salem, R. Folman, and J. Schmiedmayer, Science 319, 1226 (2008).
  * Japha  _et al._ [2008] Y. Japha, O. Entin-Wohlman, T. David, R. Salem, S. Aigner, J. Schmiedmayer, and R. Folman, Phys. Rev. B 77, 201407 (2008).
  * Schumm  _et al._ [2005] T. Schumm, J. Estève, C. Figl, J.-B. Trebbia, C. Aussibal, H. Nguyen, D. Mailly, I. Bouchoule, C. I. Westbrook, and A. Aspect, The European Physical Journal D 32, 171 (2005).
  * Henkel  _et al._ [2003] C. Henkel, P. Krüger, R. Folman, and J. Schmiedmayer, Applied Physics B: Lasers and Optics 76, 173 (2003).
  * Lin  _et al._ [2004] Y.-j. Lin, I. Teper, C. Chin, and V. Vuletić, Phys. Rev. Lett. 92, 050404 (2004).
  * Laliotis  _et al._ [2021b] A. Laliotis, B.-S. Lu, M. Ducloy, and D. Wilkowski, AVS Quantum Science 3, 043501 (2021b).
  * Wongcharoenbhorn  _et al._ [2021] K. Wongcharoenbhorn, R. Crawford, N. Welch, F. Wang, G. Sinuco-León, P. Krüger, F. Intravaia, C. Koller, and T. M. Fromhold, Phys. Rev. A 104, 053108 (2021).
  * Ergoktas  _et al._ [2021] M. Ergoktas, G. Bakan, E. Kovalska, L. Le Fevre, R. Fields, P. Steiner, X. Yu, O. Salihoglu, S. Balci, V. Fal’ko, K. Novoselov, R. Dryfe, and C. Kocabas, Nature Photonics 15, 493 (2021).
  * Liu  _et al._ [2020a] M. Liu, Z. Li, X. Zhao, R. Young, and I. Kinloch, Nano Letters 21, 833 (2020a).
  * Kim  _et al._ [2020] Y. Kim, T. Kim, J. Lee, Y. S. Choi, J. Moon, S. Y. Park, T. H. Lee, H. K. Park, S. A. Lee, and M. S. e. a. Kwon, Advanced Materials 33, 2004827 (2020).
  * Carvalho  _et al._ [2021] A. F. Carvalho, B. Kulyk, A. J. S. Fernandes, E. Fortunato, and F. M. Costa, Advanced Materials , 2101326 (2021).
  * Mišeikis  _et al._ [2020] V. Mišeikis, S. Marconi, M. A. Giambra, A. Montanaro, L. Martini, F. Fabbri, S. Pezzini, G. Piccinini, S. Forti, B. Terrés, I. Goykhman, L. Hamidouche, P. Legagneux, V. Sorianello, A. C. Ferrari, F. H. L. Koppens, M. Romagnoli, and C. Coletti, ACS Nano 14, 11190 (2020).
  * Gosciniak  _et al._ [2020] J. Gosciniak, M. Rasras, and J. B. Khurgin, ACS Photonics 7, 488 (2020).
  * Gosciniak and Khurgin [2020] J. Gosciniak and J. B. Khurgin, ACS Omega 5, 14711 (2020).
  * Zhu  _et al._ [2021] Y. Zhu, H. Xu, P. Yu, and Z. Wang, Applied Physics Reviews 8, 021305 (2021).
  * AlAloul and Rasras [2021] M. AlAloul and M. Rasras, J. Opt. Soc. Am. B 38, 602 (2021).
  * Britnell  _et al._ [2013] L. Britnell, R. V. Gorbachev, A. K. Geim, L. A. Ponomarenko, A. Mishchenko, M. T. Greenaway, T. M. Fromhold, K. S. Novoselov, and L. Eaves, Nature Communications 4, 1794 (2013).
  * Liu  _et al._ [2020b] L. Liu, Y. Liu, and X. Duan, Science China Information Sciences 63, 201401 (2020b).
  * Ribeiro and Terças [2016] S. Ribeiro and H. Terças, Phys. Rev. A 94, 043420 (2016).
  * Candussio  _et al._ [2021] S. Candussio, M. V. Durnev, S. Slizovskiy, T. Jötten, J. Keil, V. V. Bel’kov, J. Yin, Y. Yang, S.-K. Son, A. Mishchenko, V. Fal’ko, and S. D. Ganichev, Phys. Rev. B 103, 125408 (2021).
  * Tripathi  _et al._ [2021] M. Tripathi, F. Lee, A. Michail, D. Anestopoulos, J. G. McHugh, S. P. Ogilvie, M. J. Large, A. A. Graf, P. J. Lynch, J. Parthenios, K. Papagelis, S. Roy, M. A. S. R. Saadi, M. M. Rahman, N. M. Pugno, A. A. K. King, P. M. Ajayan, and A. B. Dalton, ACS Nano 15, 2520 (2021).
  * Ta  _et al._ [2020] H. Q. Ta, A. Bachmatiuk, R. G. Mendes, D. J. Perello, L. Zhao, B. Trzebicka, T. Gemming, S. V. Rotkin, and M. H. Rümmeli, Advanced Materials 32, 2002755 (2020).
  * Wei  _et al._ [2021] T. Wei, F. Hauke, and A. Hirsch, Advanced Materials 33, 2104060 (2021).
  * Cui  _et al._ [2021a] L. Cui, J. Wang, and M. Sun, Reviews in Physics 6, 100054 (2021a).
  * Liu  _et al._ [2021a] M. Liu, Y. Zhang, G. L. Klimchitskaya, V. M. Mostepanenko, and U. Mohideen, Phys. Rev. B 104, 085436 (2021a).
  * Slizovskiy  _et al._ [2021] S. Slizovskiy, A. Garcia-Ruiz, A. I. Berdyugin, N. Xin, T. Taniguchi, K. Watanabe, A. K. Geim, N. D. Drummond, and V. I. Fal’ko, Nano Letters 21, 6678 (2021).
  * Greenaway  _et al._ [2021] M. T. Greenaway, P. Kumaravadivel, J. Wengraf, L. A. Ponomarenko, A. I. Berdyugin, J. Li, J. H. Edgar, R. K. Kumar, A. K. Geim, and L. Eaves, Nature Communications 12, 6392 (2021).
  * McCann [2012] E. McCann,  _Electronic Properties of Monolayer and Bilayer Graphene_ (Springer Berlin Heidelberg, 2012) pp. 237–275.
  * Fermani  _et al._ [2007] R. Fermani, S. Scheel, and P. L. Knight, Phys. Rev. A 75, 062905 (2007).
  * Thomsen  _et al._ [2017] J. D. Thomsen, T. Gunst, S. S. Gregersen, L. Gammelgaard, B. S. Jessen, D. M. A. Mackenzie, K. Watanabe, T. Taniguchi, P. Bøggild, and T. J. Booth, Phys. Rev. B 96, 014101 (2017).
  * Ju  _et al._ [2011] L. Ju, B. Geng, J. Horng, C. Girit, M. Martin, Z. Hao, H. A. Bechtel, X. Liang, A. Zettl, Y. R. Shen, and F. Wang, Nature Nanotechnology 6, 630 (2011).
  * Yan  _et al._ [2013] H. Yan, T. Low, W. Zhu, Y. Wu, M. Freitag, X. Li, F. Guinea, P. Avouris, and F. Xia, Nature Photonics 7, 394 (2013).
  * Cui  _et al._ [2021b] L. Cui, J. Wang, and M. Sun, Reviews in Physics 6, 100054 (2021b).
  * Haakh  _et al._ [2014] H. R. Haakh, C. Henkel, S. Spagnolo, L. Rizzuto, and R. Passante, Phys. Rev. A 89, 022509 (2014).
  * Zeng and Zubairy [2021] X. Zeng and M. S. Zubairy, Phys. Rev. Lett. 126, 117401 (2021).
  * Cysne  _et al._ [2014] T. Cysne, W. J. M. Kort-Kamp, D. Oliver, F. A. Pinheiro, F. S. S. Rosa, and C. Farina, Phys. Rev. A 90, 052511 (2014).
  * Liu  _et al._ [2021b] M. Liu, Y. Zhang, G. L. Klimchitskaya, V. M. Mostepanenko, and U. Mohideen, Phys. Rev. B 104, 085436 (2021b).
  * Henkel  _et al._ [2018] C. Henkel, G. L. Klimchitskaya, and V. M. Mostepanenko, Phys. Rev. A 97, 032504 (2018).
  * Chaichian  _et al._ [2012] M. Chaichian, G. L. Klimchitskaya, V. M. Mostepanenko, and A. Tureanu, Phys. Rev. A 86, 012515 (2012).
  * Churkin  _et al._ [2010] Y. V. Churkin, A. B. Fedortsov, G. L. Klimchitskaya, and V. A. Yurova, Phys. Rev. B 82, 165433 (2010).
  * Bordag  _et al._ [2006] M. Bordag, B. Geyer, G. L. Klimchitskaya, and V. M. Mostepanenko, Phys. Rev. B 74, 205431 (2006).
  * Cysne  _et al._ [2016] T. P. Cysne, T. G. Rappoport, A. Ferreira, J. M. V. P. Lopes, and N. M. R. Peres, Phys. Rev. B 94, 235405 (2016).
  * Nichols  _et al._ [2016] N. S. Nichols, A. Del Maestro, C. Wexler, and V. N. Kotov, Phys. Rev. B 93, 205412 (2016).
  * Klimchitskaya and Mostepanenko [2020] G. L. Klimchitskaya and V. M. Mostepanenko, Universe 6, 150 (2020).
  * Kaur  _et al._ [2014] K. Kaur, J. Kaur, B. Arora, and B. K. Sahoo, Phys. Rev. B 90, 245405 (2014).
  * Khusnutdinov  _et al._ [2018] N. Khusnutdinov, R. Kashapov, and L. M. Woods, 2D Materials 5, 035032 (2018).
  * Šibalić and Adams [2018] N. Šibalić and C. S. Adams,  _Rydberg Physics_, 2399-2891 (IOP Publishing, 2018).
  * Kohlhoff [2016] M. W. Kohlhoff, The European Physical Journal Special Topics 225, 3061 (2016).
  * Hermann-Avigliano  _et al._ [2014] C. Hermann-Avigliano, R. C. Teixeira, T. L. Nguyen, T. Cantat-Moltrecht, G. Nogues, I. Dotsenko, S. Gleyzes, J. M. Raimond, S. Haroche, and M. Brune, Phys. Rev. A 90, 040502 (2014).
  * Wu  _et al._ [2021] X. Wu, X. Liang, Y. Tian, F. Yang, C. Chen, Y.-C. Liu, M. K. Tey, and L. You, Chinese Physics B 30, 020305 (2021).
  * Ellingsen  _et al._ [2010] S. A. Ellingsen, S. Y. Buhmann, and S. Scheel, Phys. Rev. Lett. 104, 223003 (2010).
  * Törmä and Barnes [2014] P. Törmä and W. L. Barnes, Reports on Progress in Physics 78, 013901 (2014).
  * Kübler  _et al._ [2013] H. Kübler, D. Booth, J. Sedlacek, P. Zabawa, and J. P. Shaffer, Phys. Rev. A 88, 043810 (2013).
  * Ribeiro  _et al._ [2013] S. Ribeiro, S. Y. Buhmann, and S. Scheel, Phys. Rev. A 87, 042508 (2013).
  * Sheng  _et al._ [2016] J. Sheng, Y. Chao, and J. P. Shaffer, Phys. Rev. Lett. 117, 103201 (2016).
  * Barton [1974] G. Barton, Journal of Physics B: Atomic and Molecular Physics 7, 2134 (1974).
  * Power and Thirunamachandran [1982] E. A. Power and T. Thirunamachandran, Phys. Rev. A 25, 2473 (1982).
  * Barton [1988] G. Barton, Physica Scripta T21, 11 (1988).
  * Meschede  _et al._ [1990] D. Meschede, W. Jhe, and E. A. Hinds, Phys. Rev. A 41, 1587 (1990).
  * Hinds and Sandoghdar [1991] E. A. Hinds and V. Sandoghdar, Phys. Rev. A 43, 398 (1991).
  * Bambini and Robinson [1992] A. Bambini and E. J. Robinson, Phys. Rev. A 45, 4661 (1992).
  * Fichet  _et al._ [1995] M. Fichet, F. Schuller, D. Bloch, and M. Ducloy, Phys. Rev. A 51, 1553 (1995).
  * Gorza  _et al._ [2001] M.-P. Gorza, S. Saltiel, H. Failache, and M. Ducloy, The European Physical Journal D 15, 113 (2001).
  * Bushev  _et al._ [2004] P. Bushev, A. Wilson, J. Eschner, C. Raab, F. Schmidt-Kaler, C. Becher, and R. Blatt, Phys. Rev. Lett. 92, 223602 (2004).
  * Buhmann  _et al._ [2004] S. Y. Buhmann, L. Knöll, D.-G. Welsch, and H. T. Dung, Phys. Rev. A 70, 052117 (2004).
  * Gorza and Ducloy [2006] M.-P. Gorza and M. Ducloy, The European Physical Journal D 40, 343 (2006).
  * Buhmann and Welsch [2008] S. Y. Buhmann and D.-G. Welsch, Phys. Rev. A 77, 012110 (2008).
  * Schiefele and Henkel [2010] J. Schiefele and C. Henkel, Phys. Rev. A 82, 023605 (2010).
  * Crosse  _et al._ [2010] J. A. Crosse, S. Å. Ellingsen, K. Clements, S. Y. Buhmann, and S. Scheel, Phys. Rev. A 82, 010901 (2010).
  * Antezza  _et al._ [2014] M. Antezza, C. Braggio, G. Carugno, A. Noto, R. Passante, L. Rizzuto, G. Ruoso, and S. Spagnolo, Phys. Rev. Lett. 113, 023601 (2014).
  * Ribeiro  _et al._ [2015] S. Ribeiro, S. Y. Buhmann, T. Stielow, and S. Scheel, EPL (Europhysics Letters) 110, 51003 (2015).
  * Armata  _et al._ [2016] F. Armata, R. Vasile, P. Barcellona, S. Y. Buhmann, L. Rizzuto, and R. Passante, Phys. Rev. A 94, 042511 (2016).
  * Carvalho  _et al._ [2018] J. C. d. A. Carvalho, P. Pedri, M. Ducloy, and A. Laliotis, Phys. Rev. A 97, 023806 (2018).
  * Yang  _et al._ [2019] B. Yang, B. Zhang, Z. Liu, and H. Yao, Journal of Physics B: Atomic, Molecular and Optical Physics 52, 095501 (2019).
  * Fang  _et al._ [2019] W. Fang, G.-X. Li, J. Xu, and Y. Yang, Opt. Express 27, 37753 (2019).
  * Buhmann  _et al._ [2021] S. Y. Buhmann, S. M. Giesen, M. Diekmann, R. Berger, S. Aull, P. Zahariev, M. Debatin, and K. Singer, New Journal of Physics 23, 083040 (2021).
  * Fuchs  _et al._ [2018] S. Fuchs, R. Bennett, and S. Y. Buhmann, Phys. Rev. A 98, 022514 (2018).
  * Wilson  _et al._ [2003] M. A. Wilson, P. Bushev, J. Eschner, F. Schmidt-Kaler, C. Becher, R. Blatt, and U. Dorner, Phys. Rev. Lett. 91, 213602 (2003).
  * Barcellona  _et al._ [2016] P. Barcellona, R. Passante, L. Rizzuto, and S. Y. Buhmann, Phys. Rev. A 93, 032508 (2016).
  * Ribeiro and Scheel [2013] S. Ribeiro and S. Scheel, Phys. Rev. A 88, 042519 (2013).
  * M.  _et al._ [2011] M. M. M., H. R. Haakh, T. Calarco, C. P. Koch, and C. Henkel, Quantum Information Processing 10, 771 (2011).
  * Wade  _et al._ [2016] C. G. Wade, Šibalić N., N. R. de Melo, J. M. Kondo, C. S. Adams, and K. J. Weatherill, Nature Photonics 11, 40 (2016).
  * Saffman [2016] M. Saffman, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 202001 (2016).
  * Scheel and Buhmann [2009] S. Scheel and S. Y. Buhmann, Phys. Rev. A 80, 042902 (2009).
  * Scheel and Buhmann [2008] S. Scheel and S. Buhmann, Acta Physica Slovaca. 58, 675 (2008).
  * Buhmann [2012] S. Y. Buhmann,  _Dispersion Forces II: Many-Body Effects, Excited Atoms, Finite Temperature and Quantum Friction_ (Springer, Berlin, 2012).
  * Buhmann and Scheel [2008] S. Y. Buhmann and S. Scheel, Phys. Rev. Lett. 100, 253201 (2008).
  * Gallagher [1988] T. F. Gallagher, Reports on Progress in Physics 51, 143 (1988).
  * Seaton [1983] M. J. Seaton, Reports on Progress in Physics 46, 167 (1983).
  * Robertson  _et al._ [2021] E. Robertson, N. Šibalić, R. Potvliege, and M. Jones, Computer Physics Communications 261, 107814 (2021).
  * Numerov [1927] B. Numerov, Astronomische Nachrichten 230, 359 (1927).
  * Li  _et al._ [2003] W. Li, I. Mourachko, M. W. Noel, and T. F. Gallagher, Phys. Rev. A 67, 052502 (2003).
  * Gric [2019] T. Gric, Opt Quant Electron 51, 202 (2019).
  * Andryieuski and Lavrinenko [2013] A. Andryieuski and A. V. Lavrinenko, Optics Express 21, 9144 (2013).
  * Gonçalves and Peres [2016] P. A. D. Gonçalves and N. M. R. Peres,  _An Introduction to Graphene Plasmonics_ (World Scientific, Singapore, 2016).
  * Amorim  _et al._ [2017] B. Amorim, P. A. D. Gonçalves, M. I. Vasilevskiy, and N. M. R. Peres, Applied Sciences 7, 1158 (2017).
  * Ellingsen  _et al._ [2009] S. A. Ellingsen, S. Y. Buhmann, and S. Scheel, Phys. Rev. A 79, 052903 (2009).
  * Klimchitskaya and Mostepanenko [2018] G. L. Klimchitskaya and V. M. Mostepanenko, Phys. Rev. A 98, 032506 (2018).
  * Rodrigo  _et al._ [2017] D. Rodrigo, A. Tittl, O. Limaj, F. J. G. d. Abajo, V. Pruneri, and H. Altug, Light: Science & Applications 6, 10.1038/lsa.2016.277 (2017).
  * Woessner  _et al._ [2014] A. Woessner, M. B. Lundeberg, Y. Gao, A. Principi, P. Alonso-González, M. Carrega, K. Watanabe, T. Taniguchi, G. Vignale, M. Polini, and et al., Nature Materials 14, 421 (2014).
  * Brar  _et al._ [2014] V. W. Brar, M. S. Jang, M. Sherrott, S. Kim, J. J. Lopez, L. B. Kim, M. Choi, and H. Atwater, Nano Letters 14, 3876 (2014), publisher: American Chemical Society.
  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * Han  _et al._ [2019] X. Han, J. Lin, J. Liu, N. Wang, and D. Pan, The Journal of Physical Chemistry C 123, 14797 (2019).
  * Yankowitz  _et al._ [2014] M. Yankowitz, J. Xue, and B. J. LeRoy, Journal of Physics: Condensed Matter 26, 303201 (2014).
  * Wrigley  _et al._ [2021] J. Wrigley, J. Bradford, T. James, T. S. Cheng, J. Thomas, C. J. Mellor, A. N. Khlobystov, L. Eaves, C. T. Foxon, S. V. Novikov, and P. H. Beton, 2D Materials 8, 034001 (2021).
  * Laturia  _et al._ [2018] A. Laturia, M. L. Van de Put, and W. G. Vandenberghe, npj 2D Materials and Applications 2, 10.1038/s41699-018-0050-x (2018).
  * Stauber  _et al._ [2008] T. Stauber, N. M. R. Peres, and A. K. Geim, Phys. Rev. B 78, 085432 (2008).
  * Hanson [2013] G. W. Hanson, Journal of Applied Physics 113, 029902 (2013).
  * Lindhard [1954] J. Lindhard, Dan. Mat. Fys. Medd. 28, 1 (1954).
  * Wunsch  _et al._ [2006] B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New Journal of Physics 8, 318 (2006).
  * Hwang and Das Sarma [2007] E. H. Hwang and S. Das Sarma, Phys. Rev. B 75, 205418 (2007).
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Saunders College Publishing, 1976).
  * Zeman and Schatz [1987] E. J. Zeman and G. C. Schatz, The Journal of Physical Chemistry 91, 634 (1987).