# The Casimir effect for stack of graphenes

Natalia Emelianova<EMAIL_ADDRESS>CMCC, Universidade Federal do
ABC, Avenida dos Estados 5001, CEP 09210-580, SP, Brazil  Rashid Kashapov
<EMAIL_ADDRESS>Eidos-Robotics LLC , Peterburgskaya 50, 420107,
Kazan, Russia  Nail Khusnutdinov<EMAIL_ADDRESS>CMCC,
Universidade Federal do ABC, Avenida dos Estados 5001, CEP 09210-580, SP,
Brazil

###### Abstract

We consider a stack of parallel sheets composed of conducting planes with
tensorial conductivities. Using the scattering matrix approach, we derive
explicit formulas for the Casimir energy of two, three, and four planes, as
well as a recurrence relation for arbitrary planes. Specifically, for a stack
of graphene, we solve the recurrence relations and obtain formulas for the
Casimir energy and force acting on the planes within the stack. Moreover, we
calculate the binding energy in the graphene stack with graphite interplane
separation, which amounts to $E_{ib}=9.9$ meV/atom. Notably, the Casimir force
on graphene sheets decreases rapidly for planes beyond the first one. In
particular, for the second graphene layer in the stack, the force is $35$
times smaller than that experienced by the first layer.

###### pacs:

03.70.+k, 03.50.De, 68.65.Pq

##  I Introduction

The Casimir effect [1], which was originally considered for two perfect slabs,
now plays an important role in various phenomena in physics, chemistry, and
biology. For instance, books [2, *Milton:2001:CEPMZE,
*Parsegian:2006:VdWFHBCEP, *Milonni:1994:QVItQE] and recent reviews [6] have
highlighted its significance in these fields. The Casimir force, which can
have different values and signs depending on the types of materials used, such
as graphene [7], topological insulators [8, *Lu:2018:cdWtfbatis], chiral
metamaterials [10], Weyl semimetals [11,
*RodriguezLopez:2020:socoriciotiaiws], the shape of material boundaries, and
the external conditions including temperature, chemical potential, and
magnetic field, has demonstrated its importance in chemistry [13], biophysics
[14, *Pawlowski:2013:qcemufobscm], and in layered systems such as graphite.

Since the 1970s, van der Waals/Casimir energy and force have been studied for
multilayered periodic systems, including a periodic stack of dielectric
materials [16, *Ninham:1970:vdWFaTF]. The Lifshitz formula has been
generalized for layered dielectric systems using nonstandard recursion
relations for Fresnel coefficients and the fluctuation-dissipative theorem
[18, *Tomas:2012:cealm]. The Casimir force acting on a plane in a stack
perpendicular to the $z$-axis was calculated as the difference in the $(z,z)$
component of the regularized stress tensor of neighbouring planes. In Ref.
[20], a plane in a stack was considered as a piston for a five-plane system.
The Casimir force for a five-layered magnetodielectric planar system was
discussed in Ref. [21] as an ideal system for detection of the temperature
dependence of the Casimir force, and a path-integral approach was used in Ref.
[22, *Amooghorban:2011:cfmmbgl] to calculate the force for a magnetodielectric
layered planar system. Different geometries of multilayered systems were
considered in Ref. [24] in the framework of normal-mode techniques, and a
modal approach was developed for layered materials in Ref. [25]. Finally, Ref.
[26] explored the thermal and electrostatic manipulation of the Casimir force
in graphene-dielectric multilayers and demonstrated the possibility of
consistent modulation of the Casimir pressure by changing the number of
graphene sheets in the stack.

The study presented herein deals with the calculation of Casimir energy for a
layered system with a different shape, i.e., a set of cylindrical concentric
shells, which has been previously discussed in Ref. [27]. In Ref. [28,
*Khusnutdinov:2015:cefasocp, *Khusnutdinov:2016:cpefasocp,
*Kashapov:2016:tcefpls, *Khusnutdinov:2018:tcaciinp2dm,
*Khusnutdinov:2019:cei2dmss], the layered system of conductive planes was
analysed in details. Two models of conductivity were taken into consideration,
namely the constant conductivity and the Drude–Lorentz 7-oscillator model of
graphene conductivity, based on the conductivity of graphite [34].

It is important to note that, in all of the aforementioned references, the
scalar type of materials was considered. Conversely, the conductivity of
graphene has a tensorial form [7, 35, *Bordag:2016:ECefdg]. For two planes
with matrix Fresnel reflection coefficients, the formula for the Casimir
energy, including some variations and a correct derivation, can be found in
Ref. [37, *Reynaud:2010:sacf, *Ingold:2015:cesa] and Ref. [40], respectively.
In a system consisting of two planes without virtual photon production, this
formula can be represented as an integral over imaginary frequency
$\xi=-\mathrm{i}\omega$,

$E_{2}=\Re\iintop\ilimits@\frac{\mathrm{d}^{2}k}{(2\pi)^{3}}\int_{0}^{\infty}\mathrm{d}\xi\ln\det\left[\mathbf{I}-e^{-2dk_{E}}\bm{r}^{\prime}_{1}\bm{r}_{2}\right],$
(1)

where $\bm{r}^{\prime}_{1}$ and $\bm{r}_{2}$ are the Fresnel reflection
matrices for the first and second planes, $d$ represents the interplane
distance, and $k_{E}=\sqrt{\xi^{2}+\bm{k}^{2}}$. The generalization of this
expression for more plains is not obvious and requires careful derivation,
which is the main focus of this paper.

In Sec. II, the scattering problem for a layered system containing
$\mathsf{n}$ parallel conductive planes with tensorial conductivity is
formulated. A general expression of Casimir’s energy is obtained, and it is
calculated for $\mathsf{n}=2,3$, and $4$ planes with different conductivities
and interplane distances.

In Sec. II, recurrent relations for the energy are derived, and a relation for
energies for $\mathsf{n},\mathsf{n}-1,\mathsf{n}-2$ planes is obtained. In the
case of a stack of graphenes with equal interplane distance, this recurrent
relation can be solved in manifest form. By using eigenvalues of the
reflection Fresnel matrix, the matrix form of solutions is represented as a
sum of scalar forms over eigenvalues. Analogous to Casimir’s energy, the
Casimir force acting on a specific plane in the stack is determined.

Finally, Sec. IV is dedicated to a particular case of a stack of graphenes,
where the Casimir energy and force are numerically calculated at zero
temperature, taking into account the conductivity tensor determined within the
framework of 3D QED [7].

Throughout the paper, the units where $\hbar=c=1$ are used.

##  II Scattering problem

Let us consider the scattering problem for the system of $\mathsf{n}$
conductive planes, which are perpendicular to the axis $z$ in the points
$z_{1}<z_{2}<\cdots<z_{\mathsf{n}}$. This system divides all space on
$\mathsf{n}+1$ domains,
$(D_{1}|z_{1}|D_{2}|z_{2}|\ldots,|z_{\mathsf{n}}|D_{\mathsf{n}+1})$. For a
plane $j$ at position $z=z_{j}$ (see Fig. 1) we have the following scattering
problem

$\begin{pmatrix}\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{j}\\\
\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{j+1}\end{pmatrix}=\bm{\mathcal{S}}_{j}\begin{pmatrix}\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{j}\\\
\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{j+1}\end{pmatrix},$ (2)

where

$\bm{\mathcal{S}}_{j}=\begin{pmatrix}\bm{r}_{j}&\bm{t}^{\prime}_{j}\\\
\bm{t}_{j}&\bm{r}^{\prime}_{j}\end{pmatrix},$ (3)

is the scattering $4\times 4$ matrix on the plain $j$
($j=1,\ldots,\mathsf{n}$). Here, $\bm{r}_{j}$ and $\bm{t}_{j}$ are the
$2\times 2$ reflection and transmission coefficients of plane $j$,
correspondingly. The
$\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{j}(\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{j})$
is 2 component vector $(E_{x},E_{y})^{T}$ in the domain $D_{j}$. The movement
to the right $\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{j}$ (left
$\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{j}$), is defined by factor
$e^{+\mathrm{i}k_{z}z}(e^{-\mathrm{i}k_{z}z})$,

Figure 1: The scattering on the plane $j$ is shown. The waves moving from the
left (right) to the right (left) are defined by factor
$e^{+\mathrm{i}k_{z}z}(e^{-\mathrm{i}k_{z}z})$.

To simplify notation, we do not write out the factors with a position of
planes. They may be restored by replacements

$\bm{r}_{j}\to\bm{r}_{j}e^{2\mathrm{i}k_{z}z_{j}},\
\bm{r}^{\prime}_{j}\to\bm{r}^{\prime}_{j}e^{-2\mathrm{i}k_{z}z_{j}},\
\bm{t}_{j}\to\bm{t}_{j},\ \bm{t}^{\prime}_{j}\to\bm{t}^{\prime}_{j}.$ (4)

The scattering matrix of total system, $\bm{\mathcal{S}}$, is defined by
relation

$\begin{pmatrix}\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{1}\\\
\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{\mathsf{n}+1}\end{pmatrix}=\bm{\mathcal{S}}\begin{pmatrix}\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{1}\\\
\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{\mathsf{n}+1}\end{pmatrix}=\begin{pmatrix}\mathbf{R}&\mathbf{T}^{\prime}\\\
\mathbf{T}&\mathbf{R}^{\prime}\end{pmatrix}\begin{pmatrix}\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{1}\\\
\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{\mathsf{n}+1}\end{pmatrix}.$ (5)

To obtain matrix $\bm{\mathcal{S}}$, we represent the set of connected
scattering equations (2) into form $\mathbf{M}\cdot\mathbf{E}=0$, where

$\mathbf{M}=\begin{pmatrix}\mathbf{A}_{1}&\mathbf{C}_{1}&0&0&\cdots&0&0\\\
0&\mathbf{A}_{2}&\mathbf{C}_{2}&0&\cdots&0&0\\\
0&0&\mathbf{A}_{3}&\mathbf{C}_{3}&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\ddots&\vdots\\\
0&0&0&0&\cdots&\mathbf{A}_{\mathsf{n}+1}&\mathbf{C}_{\mathsf{n}+1}\end{pmatrix},\
\mathbf{E}=\begin{pmatrix}\mathbf{E}_{1}\\\ \mathbf{E}_{2}\\\
\mathbf{E}_{3}\\\ \vdots\\\ \mathbf{E}_{\mathsf{n}+1}\end{pmatrix},$ (6)

and

$\mathbf{A}_{i}=\begin{pmatrix}\bm{r}_{i}&-\mathbf{I}\\\
\bm{t}_{i}&0\end{pmatrix},\mathbf{C}_{i}=\begin{pmatrix}0&\bm{t}^{\prime}_{i}\\\
-\mathbf{I}&\bm{r}^{\prime}_{i}\end{pmatrix},\mathbf{E}_{i}=\begin{pmatrix}\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{i}\\\
\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{i}\end{pmatrix}.$ (7)

Here, $\mathbf{I}$ is the identity matrix $2\times 2$. To obtain a relation
between $\mathbf{E}_{1}$ and $\mathbf{E}_{\mathsf{n}+1}$, we make elementary
transformations of lines in the matrix $\mathbf{M}$ to the form where the last
line has the first and last elements, only. To do this, we sequentially will
make zeros blocks, starting with the third line.

Let us illustrate these transformations by using an example with three planes.
We multiply the second line of $\mathbf{M}$ on the left on
$-\mathbf{A}_{3}\mathbf{C}_{2}^{-1}$ and add it to third line of the matrix:

$\mathbf{M}=\begin{pmatrix}\mathbf{A}_{1}&\mathbf{C}_{1}&0&0\\\
0&\mathbf{A}_{2}&\mathbf{C}_{2}&0\\\
0&-\mathbf{A}_{3}\mathbf{C}_{2}^{-1}\mathbf{A}_{2}&0&\mathbf{C}_{3}\end{pmatrix}.$
(8)

Next, we multiply on the left the first line on
$\mathbf{A}_{3}\mathbf{C}_{2}^{-1}\mathbf{A}_{2}\mathbf{C}_{1}^{-1}$ and add
it to the third line and then multiply the third line on
$\mathbf{C}_{3}^{-1}$:

$\mathbf{M}=\begin{pmatrix}\mathbf{A}_{1}&\mathbf{C}_{1}&0&0\\\
0&\mathbf{A}_{2}&\mathbf{C}_{2}&0\\\
\mathbf{C}_{3}^{-1}\mathbf{A}_{3}\mathbf{C}_{2}^{-1}\mathbf{A}_{2}\mathbf{C}_{1}^{-1}\mathbf{A}_{1}&0&0&\mathbf{I}\end{pmatrix}.$
(9)

The last line corresponds to equations

$\mathbf{K}_{3}\begin{pmatrix}\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{1}\\\
\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{1}\end{pmatrix}=\begin{pmatrix}\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{4}\\\
\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{4}\end{pmatrix},$ (10)

where

$\mathbf{K}_{3}=(\mathbf{C}_{3}^{-1}\mathbf{A}_{3})(\mathbf{C}_{2}^{-1}\mathbf{A}_{2})(\mathbf{C}_{1}^{-1}\mathbf{A}_{1})=\mathbf{B}_{3}\mathbf{B}_{2}\mathbf{B}_{1},$
(11)

and $\mathbf{B}_{i}=\mathbf{C}_{i}^{-1}\mathbf{A}_{i}$ (the general structure
of these matrices may be found in Appendix A). For a system with $\mathsf{n}$
planes we obtain

$(-1)^{\mathsf{n}+1}\mathbf{K}_{\mathsf{n}}\begin{pmatrix}\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{1}\\\
\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{1}\end{pmatrix}=\begin{pmatrix}\stackrel{{{}_{\rightarrow}}}{{\mathbf{E}}}_{\mathsf{n}+1}\\\
\stackrel{{{}_{\leftarrow}}}{{\mathbf{E}}}_{\mathsf{n}+1}\end{pmatrix},$ (12)

where

$\mathbf{K}_{\mathsf{n}}=\mathbf{B}_{\mathsf{n}}\mathbf{B}_{\mathsf{n}-1}\cdots\mathbf{B}_{1},$
(13)

are the matrices $4\times 4$. For $\mathsf{n}=0$, there is no scattering at a
whole, and therefore $\mathbf{K}_{0}=\mathbf{I}_{4}$, the identity matrix
$4\times 4$. By using simple algebra, we can represent this system in the form
(5) and obtain the scattering matrix of the total system in terms of the
matrix $\mathbf{K}$ (we omit index $\mathsf{n}$ in $\mathbf{K}_{\mathsf{n}}$
for simplicity)

$\displaystyle\mathbf{R}$
$\displaystyle=-\mathbf{K}_{22}^{-1}\mathbf{K}_{21},\
\mathbf{R}^{\prime}=\mathbf{K}_{12}\mathbf{K}_{22}^{-1},\
\mathbf{T}^{\prime}=(-1)^{\mathsf{n}+1}\mathbf{K}_{22}^{-1},\allowdisplaybreaks$
$\displaystyle\mathbf{T}$
$\displaystyle=(-1)^{\mathsf{n}+1}\left(\mathbf{K}_{11}-\mathbf{K}_{12}\mathbf{K}_{22}^{-1}\mathbf{K}_{21}\right).$
(14)

For $\mathsf{n}=0$ we obtain $\mathbf{R}=\mathbf{R}^{\prime}=0$ and
$\mathbf{T}=\mathbf{T}^{\prime}=\mathbf{I}$ as should be the case. The
determinant of this block-matrix reads

$\det\mathbf{K}=\det(\mathbf{K}_{22})\det(\mathbf{K}_{11}-\mathbf{K}_{12}\mathbf{K}_{22}^{-1}\mathbf{K}_{21})=\frac{\det\mathbf{T}}{\det\mathbf{T}^{\prime}}.$
(15)

As noted in Ref. [40], due to symmetry, the scattering in opposite directions
has to be taken into account. All quantities in the opposite direction we will
mark by asterisk over $\accentset{*}{g}=g(-k_{z})$. In Ref. [40], it was shown
that

$\displaystyle\accentset{*}{\mathbf{R}}^{\prime}$
$\displaystyle=(\mathbf{R}^{\prime}-\mathbf{T}\mathbf{R}^{-1}\mathbf{T}^{\prime})^{-1},\
\accentset{*}{\mathbf{R}}=(\mathbf{R}-\mathbf{T}^{\prime}\mathbf{R}^{\prime-1}\mathbf{T})^{-1},\allowdisplaybreaks$
$\displaystyle\accentset{*}{\mathbf{T}}^{\prime}$
$\displaystyle=(\mathbf{T}-\mathbf{R}^{\prime}\mathbf{T}^{\prime-1}\mathbf{R})^{-1},\
\accentset{*}{\mathbf{T}}=(\mathbf{T}^{\prime}-\mathbf{R}\mathbf{T}^{-1}\mathbf{R}^{\prime})^{-1},$
(16)

where the right hand side is calculated for $+k_{z}$. Changing $k_{z}\to-
k_{z}$ we obtain the inverse formulas

$\displaystyle\mathbf{R}^{\prime}$
$\displaystyle=(\accentset{*}{\mathbf{R}}^{\prime}-\accentset{*}{\mathbf{T}}\accentset{*}{\mathbf{R}}^{-1}\accentset{*}{\mathbf{T}}^{\prime})^{-1},\
\mathbf{R}=(\accentset{*}{\mathbf{R}}-\accentset{*}{\mathbf{T}}^{\prime}\accentset{*}{\mathbf{R}}^{\prime-1}\accentset{*}{\mathbf{T}})^{-1},\allowdisplaybreaks$
$\displaystyle\mathbf{T}^{\prime}$
$\displaystyle=(\accentset{*}{\mathbf{T}}-\accentset{*}{\mathbf{R}}^{\prime}\accentset{*}{\mathbf{T}}^{\prime-1}\accentset{*}{\mathbf{R}})^{-1},\
\mathbf{T}=(\accentset{*}{\mathbf{T}}^{\prime}-\accentset{*}{\mathbf{R}}\accentset{*}{\mathbf{T}}^{-1}\accentset{*}{\mathbf{R}}^{\prime})^{-1}.$
(17)

These relations valid for a scattering matrix of each plane, too.

The Casimir energy is defined by the logarithm of the determinant of the
scattering matrix, and it may be represented in the following form [40]

$\det\bm{\mathcal{S}}=\frac{\det\mathbf{R}}{\det\accentset{*}{\mathbf{R}}^{\prime}}=\frac{\det\mathbf{T}^{\prime}}{\det\accentset{*}{\mathbf{T}}^{\prime}}.$
(18)

By using relation (14) we obtain the following formula (see Appendix A)

$\det\bm{\mathcal{S}}=\frac{\det\mathbf{T}^{\prime}}{\det\accentset{*}{\mathbf{T}}^{\prime}}=\frac{\det\accentset{*}{\mathbf{K}}_{22}}{\det\mathbf{K}_{22}}.$
(19)

The contribution to the Casimir energy,

$\ln\det\bm{\mathcal{S}}=\ln\det\accentset{*}{\mathbf{K}}_{22}-\ln\det\mathbf{K}_{22},$
(20)

is the difference of scatterings in opposite directions.

Without changing the Casimir energy, we can multiply the matrix
$\mathbf{K}_{\mathsf{n}}$ on the arbitrary non-degenerate matrix $\mathbf{W}$,
which has no dependence on the positions of planes. We use this freedom and
define the matrix $\bm{\mathcal{D}}_{n}=\mathbf{K}_{22}\mathbf{W}$ in that way
that $\bm{\mathcal{D}}_{\mathsf{n}}=\mathbf{I}+\ldots$. Then,

$\det\bm{\mathcal{S}}=\frac{\det\accentset{*}{\bm{\mathcal{D}}}_{\mathsf{n}}}{\det\bm{\mathcal{D}}_{\mathsf{n}}},$
(21)

and the Casimir energy reads

$E_{\mathsf{n}}=\Re\iintop\ilimits@\frac{\mathrm{d}^{2}k}{(2\pi)^{3}}\int_{0}^{\infty}\mathrm{d}\xi\ln\det\bm{\mathcal{D}}_{\mathsf{n}},$
(22)

where $\bm{\mathcal{D}}_{\mathsf{n}}$ is calculated at the imaginary axis
$\omega=\mathrm{i}\xi$.

The special cases with $\mathsf{n}=1,2,3,4$ planes are considered in Appendix
B. They are

$\displaystyle\bm{\mathcal{D}}_{1}$
$\displaystyle=\mathbf{I},\allowdisplaybreaks$
$\displaystyle\bm{\mathcal{D}}_{2}$
$\displaystyle=\bm{\mathcal{D}}_{2}^{(21)}=\mathbf{I}-\bm{r}_{2}\bm{r}^{\prime}_{1},\allowdisplaybreaks$
$\displaystyle\bm{\mathcal{D}}_{3}$
$\displaystyle=\bm{\mathcal{D}}_{2}^{(32)}\bm{t}^{\prime-1}_{2}\bm{\mathcal{D}}_{2}^{(21)}\bm{t}^{\prime}_{2}-\bm{r}_{3}\bm{t}_{2}\bm{r}^{\prime}_{1}\bm{t}^{\prime}_{2},\allowdisplaybreaks$
$\displaystyle\bm{\mathcal{D}}_{4}$
$\displaystyle=\bm{\mathcal{D}}_{2}^{(43)}\bm{t}^{\prime-1}_{3}\bm{\mathcal{D}}_{2}^{(32)}\bm{t}^{\prime-1}_{2}\bm{\mathcal{D}}_{2}^{(21)}\bm{t}^{\prime}_{2}\bm{t}^{\prime}_{3}-\bm{\mathcal{D}}_{2}^{(43)}\bm{t}^{\prime-1}_{3}\bm{r}_{3}\bm{t}_{2}\bm{r}^{\prime}_{1}\bm{t}^{\prime}_{2}\bm{t}^{\prime}_{3}\allowdisplaybreaks$
$\displaystyle-\bm{r}_{4}\bm{t}_{3}\bm{r}^{\prime}_{2}\bm{t}^{\prime-1}_{2}\bm{\mathcal{D}}_{2}^{(21)}\bm{t}^{\prime}_{2}\bm{t}^{\prime}_{3}-\bm{r}_{4}\bm{t}_{3}\bm{t}_{2}\bm{r}^{\prime}_{1}\bm{t}^{\prime}_{2}\bm{t}^{\prime}_{3}.$
(23)

The positions of planes are restored by replacement (4). Here
$\bm{\mathcal{D}}_{2}^{(il)}=\mathbf{I}-\bm{r}_{i}\bm{r}^{\prime}_{l}$ is the
matrix for two planes $i$ and $l$ (numeration of planes starts from one for
the first plane and up to $\mathsf{n}$ for the last one). We observe in the
above expressions that the matrices $\bm{t}$ appear for internal planes, only.
The first terms are additive contributions to the energy; the rest terms
describe the non-additivity of the Casimir energy. For example, the first term
for $\mathsf{n}=3$ planes gives the following contribution

$\ln\det(\bm{\mathcal{D}}_{2}^{(32)}\bm{t}^{\prime-1}_{2}\bm{\mathcal{D}}_{2}^{(21)}\bm{t}^{\prime}_{2})=\ln\det\bm{\mathcal{D}}_{2}^{(32)}+\ln\det\bm{\mathcal{D}}_{2}^{(21)},$
(24)

which is the sum of the contribution of pairs of neighboring planes – first
with the middle and the middle with the last one. The second term in
$\bm{\mathcal{D}}_{3}$ contains reflection matrices of the first and last
planes and transmission matrix of the middle plane. It looks like a
contribution due to the interaction of the first and last terms through the
middle.

##  III Recurrent relations

To obtain a recurrent relation for $\bm{\mathcal{D}}_{\mathsf{n}}$ we use Eq.
(13). The matrix $\mathbf{K}_{\mathsf{n}}$ has the following form

$\mathbf{K}_{\mathsf{n}+1}=\mathbf{B}_{\mathsf{n}+1}\mathbf{K}_{\mathsf{n}}.$
(25)

By using this relation for index $\mathsf{n}$, $\mathsf{n}+1$ and the manifest
form of matrix $\mathbf{B}_{\mathsf{n}}$ (78) we obtain the recurrent relation
for block $(2,2)$ the matrix $\mathbf{K}_{\mathsf{n}}$

$\mathbf{K}_{\mathsf{n}+2}^{22}=\bm{\phi}_{\mathsf{n}}\mathbf{K}_{\mathsf{n}+1}^{22}+\bm{\psi}_{\mathsf{n}}\mathbf{K}_{\mathsf{n}}^{22},$
(26)

where

$\displaystyle\bm{\phi}_{\mathsf{n}}$
$\displaystyle=-\bm{t}^{\prime-1}_{\mathsf{n}+2}\left[\mathbf{I}-\bm{r}_{\mathsf{n}+2}\left(\bm{r}^{\prime}_{\mathsf{n}+1}-\bm{t}_{\mathsf{n}+1}\bm{r}^{-1}_{\mathsf{n}+1}\bm{t}^{\prime}_{\mathsf{n}+1}\right)\right],\allowdisplaybreaks$
$\displaystyle\bm{\psi}_{\mathsf{n}}$
$\displaystyle=-\bm{t}^{\prime-1}_{\mathsf{n}+2}\bm{r}_{\mathsf{n}+2}\bm{t}_{\mathsf{n}+1}\bm{r}^{-1}_{\mathsf{n}+1}.$
(27)

We define matrices

$\bm{\mathcal{D}}_{\mathsf{n}}=(-1)^{\mathsf{n}}\bm{t}^{\prime}_{\mathsf{n}}\mathbf{K}_{\mathsf{n}}\bm{t}^{\prime}_{1}\bm{t}^{\prime}_{2}\cdots\bm{t}^{\prime}_{\mathsf{n}-1}=\mathbf{I}+\ldots,\
\bm{\mathcal{D}}_{1}=\bm{\mathcal{D}}_{0}=\mathbf{I}.$ (28)

In terms of these matrices, we obtain the relation we need

$\bm{\mathcal{D}}_{\mathsf{n}+2}=\bm{\Phi}_{\mathsf{n}}\bm{\mathcal{D}}_{\mathsf{n}+1}\bm{t}^{\prime}_{\mathsf{n}+1}+\bm{\Psi}_{\mathsf{n}}\bm{\mathcal{D}}_{\mathsf{n}}\bm{t}^{\prime}_{\mathsf{n}}\bm{t}^{\prime}_{\mathsf{n}+1},$
(29)

where

$\displaystyle\bm{\Phi}_{\mathsf{n}}$
$\displaystyle=(\mathbf{I}-\bm{r}_{\mathsf{n}+2}\bm{r}^{\prime}_{\mathsf{n}+1})\bm{t}^{\prime-1}_{\mathsf{n}+1}+\bm{r}_{\mathsf{n}+2}\bm{t}_{\mathsf{n}+1}\bm{r}^{-1}_{\mathsf{n}+1},\allowdisplaybreaks$
$\displaystyle\bm{\Psi}_{\mathsf{n}}$
$\displaystyle=-\bm{r}_{\mathsf{n}+2}\bm{t}_{\mathsf{n}+1}\bm{r}^{-1}_{\mathsf{n}+1}\bm{t}^{-1}_{\mathsf{n}}.$
(30)

The scattering matrix for conductive plane $s$, with conductivity tensor
$\bm{\eta}_{s}=2\pi\bm{\sigma}_{s}$ was found in Ref. [40]:

$\displaystyle\bm{r}_{s}$
$\displaystyle=e^{2\mathrm{i}k_{z}z_{s}}\bm{\alpha}_{s},\
\bm{r}^{\prime}_{s}=e^{-2\mathrm{i}k_{z}z_{s}}\bm{\alpha}_{s},\
\bm{t}_{s}=\mathbf{I}+\bm{\alpha}_{s},\
\bm{t}^{\prime}_{s}=\bm{t}_{s},\allowdisplaybreaks$
$\displaystyle\accentset{*}{\bm{r}}_{s}$
$\displaystyle=-e^{-2\mathrm{i}k_{z}z_{s}}(\mathbf{I}+2\bm{\alpha}_{s})^{-1}\bm{\alpha}_{s},\
\accentset{*}{\bm{r}}^{\prime}_{s}=-e^{2\mathrm{i}k_{z}z_{s}}(\mathbf{I}+2\bm{\alpha}_{s})^{-1}\bm{\alpha}_{s},\allowdisplaybreaks$
$\displaystyle\accentset{*}{\bm{t}}_{s}$
$\displaystyle=(\mathbf{I}+2\bm{\alpha}_{s})^{-1}(\mathbf{I}+\bm{\alpha}_{s}),\accentset{*}{\bm{t}}^{\prime}_{s}=\accentset{*}{\bm{t}}_{s},$
(31)

where

$\displaystyle\bm{\alpha}_{s}$
$\displaystyle=-\frac{\omega^{2}\bm{\eta}_{s}-\bm{k}\otimes\bm{k}\bm{\eta}_{s}+\mathbf{I}k_{z}\omega\det\bm{\eta}_{s}}{\omega^{2}\operatorname{tr}\bm{\eta}_{s}-(\bm{k}\bm{k}\bm{\eta}_{s})+k_{z}\omega(1+\det\bm{\eta}_{s})},\allowdisplaybreaks$
$\displaystyle\ [\bm{\alpha}_{s}]^{i}_{j}$
$\displaystyle=-\frac{\omega^{2}[\bm{\eta}_{s}]^{i}_{j}-k^{i}k^{l}[\bm{\eta}_{s}]_{lj}+\delta^{i}_{j}k_{z}\omega\det\bm{\eta}_{s}}{\omega^{2}\operatorname{tr}\bm{\eta}_{s}-(\bm{k}\bm{k}\bm{\eta}_{s})+k_{z}\omega(1+\det\bm{\eta}_{s})},$
(32)

and $\bm{k}=(k^{1},k^{2})=(k_{1},k_{2})$,
$k_{z}=\sqrt{\omega^{2}-\bm{k}^{2}}$. All of these matrices commutate with
each other.

At imaginary axes, $\omega\to\mathrm{i}\xi\ (k_{z}\to\mathrm{i}k_{E})$ we
obtain

$\bm{\alpha}_{s}=-\frac{\xi^{2}\bm{\eta}_{s}+\bm{k}\otimes\bm{k}\bm{\eta}_{s}+\mathbf{I}k_{E}\xi\det\bm{\eta}_{s}}{\xi^{2}\operatorname{tr}\bm{\eta}_{s}+(\bm{k}\bm{k}\bm{\eta}_{s})+k_{E}\xi(1+\det\bm{\eta}_{s})},$
(33)

where $k_{E}=\sqrt{\xi^{2}+\bm{k}^{2}}$.

For identical planes $\bm{\alpha}_{s}=\bm{\alpha}$ for all $s$, and

$\displaystyle\bm{\mathcal{D}}_{\mathsf{n}+2}$
$\displaystyle=\left[\mathbf{I}+e^{-2k_{E}d_{\mathsf{n}+2,\mathsf{n}+1}}(\mathbf{I}+2\bm{\alpha})\right]\bm{\mathcal{D}}_{\mathsf{n}+1}\allowdisplaybreaks$
$\displaystyle-e^{-2k_{E}d_{\mathsf{n}+2,\mathsf{n}+1}}(\mathbf{I}+\bm{\alpha})^{2}\bm{\mathcal{D}}_{\mathsf{n}},$
(34)

where $d_{i,j}=z_{j}-z_{i}$ is the distance between planes $i$ and $j$.

For identical planes $\bm{\alpha}_{s}=\bm{\alpha}$ and equal interplane
distances $d_{i,j}=d$ we obtain at the imaginary axis

$\bm{\mathcal{D}}_{\mathsf{n}+2}=\left[\mathbf{I}+e^{-2k_{E}d}(\mathbf{I}+2\bm{\alpha})\right]\bm{\mathcal{D}}_{\mathsf{n}+1}-e^{-2k_{E}d}(\mathbf{I}+\bm{\alpha})^{2}\bm{\mathcal{D}}_{\mathsf{n}},$
(35)

and

$\displaystyle\bm{\mathcal{D}}_{2}$
$\displaystyle=\mathbf{I}-e^{-2k_{E}d}\bm{\alpha}^{2},\allowdisplaybreaks$
$\displaystyle\bm{\mathcal{D}}_{3}$
$\displaystyle=(\mathbf{I}-e^{-2k_{E}d}\bm{\alpha}^{2})^{2}-e^{-4k_{E}d}\bm{\alpha}^{2}(\mathbf{I}+\bm{\alpha})^{2},\allowdisplaybreaks$
$\displaystyle\bm{\mathcal{D}}_{4}$
$\displaystyle=(\mathbf{I}-e^{-2k_{E}d}\bm{\alpha}^{2})^{3}-2e^{-4k_{E}d}(\mathbf{I}-e^{-2k_{E}d}\bm{\alpha}^{2})\bm{\alpha}^{2}(\mathbf{I}+\bm{\alpha})^{2}\allowdisplaybreaks$
$\displaystyle-e^{-6k_{E}d}\bm{\alpha}^{2}(\mathbf{I}+\bm{\alpha})^{4}.$ (36)

To solve the recurrent relations (35) we use a generation function method
extended for matrix-valued coefficients. Let us consider the matrix-valued
recurrent relations

$\bm{\mathcal{D}}_{\mathsf{n}+2}=\bm{u}\bm{\mathcal{D}}_{\mathsf{n}+1}+\bm{v}\bm{\mathcal{D}}_{n},$
(37)

with commutative matrices $\bm{u}$ and $\bm{v}$, $[\bm{u},\bm{v}]=0$, and a
generation function

$\mathbf{G}=\sum_{s=0}^{\infty}\bm{\mathcal{D}}_{s}z^{s}.$ (38)

Here,

$\bm{u}=\mathbf{I}+e^{-2k_{E}d}(\mathbf{I}+2\bm{\alpha}),\bm{v}=-e^{-2k_{E}d}(\mathbf{I}+\bm{\alpha})^{2}.$
(39)

Taking into account the recurrent relation (37), we obtain

$\mathbf{G}=\bm{\mathcal{D}}_{0}+(\bm{\mathcal{D}}_{1}-\bm{u}\bm{\mathcal{D}}_{0})z+\bm{u}\mathbf{G}z+\bm{v}\mathbf{G}z^{2}.$
(40)

Therefore,

$\mathbf{G}=\frac{\bm{\mathcal{D}}_{0}+(\bm{\mathcal{D}}_{1}-\bm{u}\bm{\mathcal{D}}_{0})z}{\mathbf{I}-\bm{u}z-\bm{v}z^{2}}.$
(41)

Then, we expand over $z$

$(\mathbf{I}-\bm{u}z-\bm{v}z^{2})^{-1}=\sum_{s=0}^{\infty}\mathbf{M}_{s}z^{s},$
(42)

where

$\mathbf{M}_{s}=\sum_{l=0}^{\left[\frac{s}{2}\right]}\bm{u}^{s-2l}\bm{v}^{l}C^{s-l}_{l}.$
(43)

Here, $[\ldots]$ is integer part, and $C^{s}_{l}$ is the binomial coefficient

$C^{s}_{l}=\frac{s!}{l!(s-l)!}.$ (44)

We set $\mathbf{M}_{-1}=0$, and from Eq. (43) we obtain
$\mathbf{M}_{0}=\mathbf{I},\mathbf{M}_{1}=\bm{u}$.

Taking into account this expansion, we have relation

$\mathbf{G}=\sum_{s=0}^{\infty}z^{s}\left[\bm{\mathcal{D}}_{0}\mathbf{M}_{s}+(\bm{\mathcal{D}}_{1}-\bm{u}\bm{\mathcal{D}}_{0})\mathbf{M}_{s-1}\right]=\sum_{s=0}^{\infty}z^{s}\bm{\mathcal{D}}_{s}.$
(45)

The solution of this relation reads

$\bm{\mathcal{D}}_{s}=\bm{\mathcal{D}}_{0}\mathbf{M}_{s}+(\bm{\mathcal{D}}_{1}-\bm{u}\bm{\mathcal{D}}_{0})\mathbf{M}_{s-1}.$
(46)

For $s=0,1$ these relations become identities.

Therefore,

$\bm{\mathcal{D}}_{\mathsf{n}}=\mathbf{M}_{\mathsf{n}}-\mathbf{M}_{\mathsf{n}-1}e^{-2k_{E}d}(\mathbf{I}+2\bm{\alpha}).$
(47)

Here, $\bm{\mathcal{D}}_{0}=\bm{\mathcal{D}}_{1}=\mathbf{I}$ and
$\bm{\mathcal{D}}_{2},\bm{\mathcal{D}}_{3}$ are given by Eq. (36).

If the matrices $\bm{u}$ and $\bm{v}$ are numbers, $u,v$, then the sum (43)
maybe calculated in close form

$M_{\mathsf{n}}=\frac{\left(u+\sqrt{u^{2}+4v}\right)^{\mathsf{n}+1}-\left(u-\sqrt{u^{2}+4v}\right)^{\mathsf{n}+1}}{2^{\mathsf{n}+1}\sqrt{u^{2}+4v}}.$
(48)

If the $\bm{u}$ and $\bm{v}$ are diagonal matrices

$\bm{u}=\operatorname{diag}(u_{1},u_{2}),\
\bm{v}=\operatorname{diag}(v_{1},v_{2}),$ (49)

then the matrix $\mathbf{M}_{s}$ is diagonal, too, with diagonal elements

$M_{\mathsf{n}}^{(i)}=\frac{\left(u_{i}+\sqrt{u_{i}^{2}+4v_{i}}\right)^{\mathsf{n}+1}-\left(u_{i}-\sqrt{u_{i}^{2}+4v_{i}}\right)^{\mathsf{n}+1}}{2^{\mathsf{n}+1}\sqrt{u_{i}^{2}+4v_{i}}}.$
(50)

The matrix $\bm{\alpha}$ may be diagonalized with eigenvalues
$r_{\mathrm{te}},r_{\mathrm{tm}}$, which correspond to the reflection
coefficients of TE and TM modes [41]. Therefore, the matrices $\bm{u}$ and
$\bm{v}$ (39) are diagonalized, too, and, the matrices $\bm{\mathcal{D}}_{n}$
are diagonal with eigenvalues

$D_{\mathsf{n}}^{\mathsf{x}}=M_{\mathsf{n}}^{\mathsf{x}}-e^{-2k_{E}d}M_{\mathsf{n}-1}^{\mathsf{x}}(1+2r_{\mathsf{x}}),$
(51)

where

$\displaystyle M_{\mathsf{n}}^{\mathsf{x}}$
$\displaystyle=\frac{\left(u_{\mathsf{x}}+\sqrt{u_{\mathsf{x}}^{2}+4v_{\mathsf{x}}}\right)^{\mathsf{n}+1}-\left(u_{\mathsf{x}}-\sqrt{u_{\mathsf{x}}^{2}+4v_{\mathsf{x}}}\right)^{\mathsf{n}+1}}{2^{\mathsf{n}+1}\sqrt{u_{\mathsf{x}}^{2}+4v_{\mathsf{x}}}},\allowdisplaybreaks$
$\displaystyle u_{\mathsf{x}}$
$\displaystyle=1+e^{-2k_{E}d}(1+2r_{\mathsf{x}}),v_{\mathsf{x}}=-e^{-2k_{E}d}(1+r_{\mathsf{x}})^{2},$
(52)

and $\mathsf{x}=\mathrm{te},\mathrm{tm}$.

The matrix $\bm{\alpha}$ has the following eigenvalues

$r_{\mathrm{te}}=-\left(1+\frac{k_{E}}{\eta\xi}\right)^{-1},\
r_{\mathrm{tm}}=-\left(1+\frac{\xi}{\eta k_{E}}\right)^{-1},$ (53)

for the constant conductivity case $\bm{\eta}=2\pi\bm{\sigma}=\eta\mathbf{I}$,
and

$r_{\mathrm{te}}=-\left(1+\frac{k_{E}}{\eta_{\mathrm{gr}}\tilde{k}\Phi}\right)^{-1},\
r_{\mathrm{tm}}=-\left(1+\frac{\tilde{k}}{\eta_{\mathrm{gr}}k_{E}\Phi}\right)^{-1},$
(54)

for a graphene at zero temperature case,

$\bm{\eta}=\eta_{\mathrm{gr}}\frac{\tilde{k}}{\xi}\left(\mathbf{I}-v_{F}^{2}\frac{\bm{k}\otimes\bm{k}}{\tilde{k}^{2}}\right)\Phi\left(\frac{\tilde{k}}{2m}\right),$
(55)

where $\eta_{\mathrm{gr}}=2\pi\sigma_{\mathrm{gr}}=\pi e^{2}/2=0.0114$,
$v_{F}$ is the Fermi velocity, and

$\Phi(y)=\frac{2}{\pi y}\left\\{1+\frac{y^{2}-1}{y}\arctan y\right\\},\
\tilde{k}=\sqrt{\xi^{2}+v_{F}^{2}\bm{k}^{2}}.$

To compare with the constant conductivity case [29] we define a new variable
$t_{\mathsf{x}}$ by relation

$r_{\mathsf{x}}=-\frac{t_{\mathsf{x}}}{1+t_{\mathsf{x}}}.$ (56)

For the constant conductivity case,

$\ t_{\mathrm{te}}=\frac{\eta\xi}{k_{E}},\ t_{\mathrm{tm}}=\frac{\eta
k_{E}}{\xi},$ (57)

and

$\ t_{\mathrm{te}}=\frac{\eta_{\mathrm{gr}}\tilde{k}\Phi}{k_{E}},\
t_{\mathrm{tm}}=\frac{\eta_{\mathrm{gr}}k_{E}\Phi}{\tilde{k}},$ (58)

for the graphene’s sheets. With this definition ($z=k_{E}d$)

$M_{\mathsf{n}}^{\mathsf{x}}=\frac{e^{-\mathsf{n}z}}{(1+t_{\mathsf{x}})^{\mathsf{n}}}\frac{1-f_{\mathsf{x}}^{2(\mathsf{n}+1)}}{f_{\mathsf{x}}^{\mathsf{n}}(1-f_{\mathsf{x}}^{2})},$
(59)

where

$f_{\mathsf{x}}=\sqrt{(\cosh z+t_{\mathsf{x}}\sinh z)^{2}-1}+\cosh
z+t_{\mathsf{x}}\sinh z.$ (60)

Taking these relations into account, we obtain

$D_{\mathsf{n}}^{\mathsf{x}}=-\frac{e^{-\mathsf{n}z}}{f^{\mathsf{n}-1}_{\mathsf{x}}(1+t_{\mathsf{x}})^{\mathsf{n}}}\left(e^{-z}\frac{1-f^{2\mathsf{n}}_{\mathsf{x}}}{1-f^{2}_{\mathsf{x}}}(1-t_{\mathsf{x}})-\frac{1-f^{2(\mathsf{n}+1)}_{\mathsf{x}}}{f_{\mathsf{x}}(1-f^{2}_{\mathsf{x}})}\right).$
(61)

After some algebra, we can transform this expression to the form obtained in
Ref. [29] for the case of scalar constant conductivity case

$D_{\mathsf{n}}^{\mathsf{x}}=-\frac{e^{-(\mathsf{n}-1)z}}{f^{\mathsf{n}-2}_{\mathsf{x}}(1+t_{\mathsf{x}})^{\mathsf{n}}}\left(e^{-z}\frac{1-f^{2(\mathsf{n}-1)}_{\mathsf{x}}}{1-f^{2}_{\mathsf{x}}}-\frac{1+t_{\mathsf{x}}}{f_{\mathsf{x}}}\frac{1-f^{2\mathsf{n}}_{\mathsf{x}}}{1-f^{2}_{\mathsf{x}}}\right).$
(62)

Therefore, we obtain that for a stack of $\mathsf{n}$ identical conductive
planes with conductivity tensor $\bm{\sigma}$ and with identical interplane
distance $d$, the Casimir energy has the following form

$E_{\mathsf{n}}=E_{\mathsf{n}}^{\mathrm{te}}+E_{\mathsf{n}}^{\mathrm{tm}}=\iintop\ilimits@\frac{\mathrm{d}^{2}k}{(2\pi)^{3}}\int_{0}^{\infty}\mathrm{d}\xi\left(\ln
D_{\mathsf{n}}^{\mathrm{te}}+\ln D_{\mathsf{n}}^{\mathrm{tm}}\right),$ (63)

where function $D_{\mathsf{n}}$ is given by Eq. (62).

For an infinite number of planes, $\mathsf{n}\to\infty$, the Casimir energy is
divergent, but the energy per unit plane,
$\overline{E}_{n}=E_{\mathsf{n}}/\mathsf{n}$, is finite it is given by the
same expression with replacement $D_{\mathsf{n}}^{\mathsf{x}}\to
D^{\mathsf{x}}$, where

$D^{\mathsf{x}}=\lim_{\mathsf{n}\to\infty}\sqrt[\mathsf{n}]{D_{\mathsf{n}}^{\mathsf{x}}}=\frac{e^{-z}f_{\mathsf{x}}}{1+t_{\mathsf{x}}}.$
(64)

Let us analyze the limits of small and large separations. To make this, we
change integrand variables, $k_{i}\to k_{i}/d$, and $\xi\to\xi/d$, to make
integrand dimensionless. The only place with $d$-dependence is the argument of
function $\Phi=\Phi(\tilde{k}/2md)$.

For small separations, $d\to 0$, we observe that $\Phi\to 1$ and the integrand
does not depend on the interplane distance and therefore $E_{\mathsf{n}}\sim
1/d^{3}$. In the zero order of the Fermi velocity, $v_{F}=0$, the relations
(58) are transformed to the constant conductivity case (57) and the energy
coincides with that obtained in Ref. [29] for the constant conductivity case.
The Drude–Lorentz model considered in Ref. [29] gives $1/d^{5/2}$ dependence.

For large separations, $d\to\infty$, the function $\Phi\approx
8\tilde{k}/(3\pi md)\to 0$. Expanding integrand in Eq. (63) over
$t_{\mathsf{x}}\ll 1$ we obtain

$E_{\mathsf{n}}=-\frac{\eta_{\mathrm{gr}}^{2}(n-1)}{225d^{3}(md)^{2}}\sim\frac{1}{d^{5}}.$
(65)

The dependence $1/(m^{2}d^{5})$ is observed in Ref. [7]. For the case of
graphene/perfect metal the energy $\sim 1/(md^{4})$

To obtain a force acting on the plane $s$ in the stack, we consider a stack of
graphene with equal distances except for plane $s$, which are shifted on the
distance $\varepsilon$ (see Fig. 2) and calculate the derivative

$\bm{\mathcal{D}}^{\prime}_{n}=\left.\frac{\partial\bm{\mathcal{D}}_{n}}{\partial\varepsilon}\right|_{\varepsilon=0}.$
(66)

For this system, the recurrent relation (III) is valid

$\displaystyle\bm{\mathcal{D}}_{l}$
$\displaystyle=\bm{u}\bm{\mathcal{D}}_{l-1}+\bm{v}\bm{\mathcal{D}}_{l-2},$
$\displaystyle 1\leqslant l\leqslant s-1,\allowdisplaybreaks$
$\displaystyle\bm{\mathcal{D}}_{s}$
$\displaystyle=\bm{u}(\varepsilon)\bm{\mathcal{D}}_{s-1}+\bm{v}(\varepsilon)\bm{\mathcal{D}}_{s-2},$
$\displaystyle l=s,\allowdisplaybreaks$ $\displaystyle\bm{\mathcal{D}}_{s+1}$
$\displaystyle=\bm{u}(-\varepsilon)\bm{\mathcal{D}}_{s}+\bm{v}(-\varepsilon)\bm{\mathcal{D}}_{s-1},$
$\displaystyle l=s+1,\allowdisplaybreaks$ $\displaystyle\bm{\mathcal{D}}_{l}$
$\displaystyle=\bm{u}\bm{\mathcal{D}}_{l-1}+\bm{v}\bm{\mathcal{D}}_{l-2},$
$\displaystyle l\geqslant s+2,$ (67)

where,

$\displaystyle\bm{u}$
$\displaystyle=\mathbf{I}+e^{-2k_{E}d}(\mathbf{I}+2\bm{\alpha}),\bm{v}=-e^{-2k_{E}d}(\mathbf{I}+\bm{\alpha})^{2},\allowdisplaybreaks$
$\displaystyle\bm{u}(\varepsilon)$
$\displaystyle=\left[\mathbf{I}+e^{-2k_{E}(d+\varepsilon)}(\mathbf{I}+2\bm{\alpha})\right],\bm{v}(\varepsilon)=-e^{-2k_{E}(d+\varepsilon)}(\mathbf{I}+\bm{\alpha})^{2}.$

The recurrent relations for derivatives have the following form

$\displaystyle\bm{\mathcal{D}}^{\prime}_{l}$ $\displaystyle=0,$ $\displaystyle
1\leqslant l\leqslant s-1,\allowdisplaybreaks$
$\displaystyle\bm{\mathcal{D}}^{\prime}_{s}$
$\displaystyle=\bm{u}^{\prime}\bm{\mathcal{D}}_{s-1}+\bm{v}^{\prime}\bm{\mathcal{D}}_{s-2},$
$\displaystyle l=s,\allowdisplaybreaks$
$\displaystyle\bm{\mathcal{D}}^{\prime}_{s+1}$
$\displaystyle=-\bm{u}^{\prime}\bm{\mathcal{D}}_{s}+\bm{u}\bm{\mathcal{D}}_{s}-\bm{v}^{\prime}\bm{\mathcal{D}}_{s-1},$
$\displaystyle l=s+1,\allowdisplaybreaks$
$\displaystyle\bm{\mathcal{D}}^{\prime}_{l}$
$\displaystyle=\bm{u}\bm{\mathcal{D}}^{\prime}_{l-1}+\bm{v}\bm{\mathcal{D}}^{\prime}_{l-2},$
$\displaystyle l\geqslant s+2,$ (68)

where

$\bm{u}^{\prime}=-2k_{E}e^{-2k_{E}d}(\mathbf{I}+2\bm{\alpha}),\bm{v}^{\prime}=2k_{E}e^{-2k_{E}d}(\mathbf{I}+\bm{\alpha})^{2}.$
(69)

The generation function for this recurrent relation starts from power $z^{s}$

$\mathbf{G}=\sum_{l=s}^{\infty}\bm{\mathcal{D}}^{\prime}_{l}z^{l}=z\bm{u}\mathbf{G}+z^{2}\bm{v}\mathbf{G}+\bm{\mathcal{D}}^{\prime}_{s}z^{s}+(\bm{\mathcal{D}}^{\prime}_{s+1}-\bm{u}\bm{\mathcal{D}}^{\prime}_{s})z^{s+1}.$
(70)

Therefore,

$\mathbf{G}=\frac{\bm{\mathcal{D}}^{\prime}_{s}z^{s}+(\bm{\mathcal{D}}^{\prime}_{s+1}-\bm{u}\bm{\mathcal{D}}^{\prime}_{s})z^{s+1}}{\mathbf{I}-z\bm{u}+z^{2}\bm{v}}.$
(71)

Then we use expansion for denominator (42) and obtain

$\bm{\mathcal{D}}^{\prime}_{\mathsf{n}}=\bm{\mathcal{D}}^{\prime}_{s}\mathbf{M}_{\mathsf{n}-s}+(\bm{\mathcal{D}}^{\prime}_{s+1}-\bm{u}\bm{\mathcal{D}}^{\prime}_{s})\mathbf{M}_{\mathsf{n}-s-1}.$
(72)

For all planes $l\leqslant s$ we have to use the solution (47).

Let us obtain now an expression for force acting on the graphene with number
$s$ in the stack of $\mathsf{n}$ graphenes

$\displaystyle F_{s,\mathsf{n}}$ $\displaystyle=-\left.\frac{\partial
E_{\mathsf{n}}}{\partial\varepsilon}\right|_{\varepsilon=0}=-\lim_{\varepsilon\to
0}\frac{E_{\mathsf{n}}(\varepsilon)-E_{\mathsf{n}}(0)}{\varepsilon}\allowdisplaybreaks$
$\displaystyle=-\Re\lim_{\varepsilon\to
0}\iintop\ilimits@\frac{\mathrm{d}^{2}k}{(2\pi)^{3}}\int_{0}^{\infty}\mathrm{d}\xi\frac{\ln\det\bm{\mathcal{D}}_{\mathsf{n}}(\varepsilon)-\ln\det\bm{\mathcal{D}}_{\mathsf{n}}(0)}{\varepsilon}\allowdisplaybreaks$
$\displaystyle=-\Re\iintop\ilimits@\frac{\mathrm{d}^{2}k}{(2\pi)^{3}}\int_{0}^{\infty}\mathrm{d}\xi\operatorname{tr}\left(\bm{\mathcal{D}}_{\mathsf{n}}^{-1}\bm{\mathcal{D}}^{\prime}_{\mathsf{n}}\right)\allowdisplaybreaks$
$\displaystyle=-\Re\sum_{\mathsf{x}=\mathrm{te},\mathrm{tm}}\iintop\ilimits@\frac{\mathrm{d}^{2}k}{(2\pi)^{3}}\int_{0}^{\infty}\mathrm{d}\xi
D^{\mathsf{x}}_{\mathsf{n}}{}^{-1}{D^{\mathsf{x}}_{\mathsf{n}}}^{\prime}.$
(73)

Direct calculation gives

$D^{\mathsf{x}}_{\mathsf{n}}{}^{-1}{D^{\mathsf{x}}_{\mathsf{n}}}^{\prime}\stackrel{{\scriptstyle\textrm{def}}}{{=}}G_{s,\mathsf{n}}^{\mathsf{x}}=\frac{2zt_{\mathsf{x}}^{2}e^{-z}f(f^{2(s-1)}-f^{2(\mathsf{n}-s)})}{e^{-z}f(1-f^{2(\mathsf{n}-1)}_{\mathsf{x}})-(1+t_{\mathsf{x}})(1-f^{2\mathsf{n}}_{\mathsf{x}})}.$
(74)

This expression is the same as obtained in the scalar case in Ref. [31]. For
infinite number of graphenes in stack, $\mathsf{n}\to\infty$, and finite $s$
we obtain

$G_{s}^{\mathsf{x}}=\frac{2zt_{\mathsf{x}}^{2}f^{2(1-s)}}{1-e^{z}f(1+t_{\mathsf{x}})}.$
(75)

Figure 2: The plane $l=s$ in the stack of $\mathsf{n}$ planes is shifted on
$\varepsilon$ to the right.

Therefore, we obtain that the expression for force has the same form as that
obtained for the scalar constant conductivity case but with corresponding
reflection coefficients for TE and TM modes.

##  IV The stack of graphenes

The numerical evaluation of the Casimir energy for a stack of $\mathsf{n}$
graphenes at zero temperature is plotted in Figs. 3, 4 as a function of
parameter $p=md$ and number of graphene in the stack. We evaluate the function
$\mathcal{E}_{\mathsf{n}}=E_{\mathsf{n}}/\mathsf{n}E_{Cas}$ – the energy per
unit graphene in units of the Casimir energy $E_{Cas}=-\pi^{2}/720d^{3}$ for
two ideal planes.

Figure 3: The Casimir energy per unit graphene $\mathcal{E}_{\mathsf{n}}$ as a
function of interplane distance. Here, $\mathsf{n}=2,4,6,\infty$ from the
bottom upwards. The insert shows the figure’s neighborhood of origin.

For $m=0.1$ eV and interplane distance in graphite $d_{c}=0.3345$ nm, the
parameter $p_{c}=1.69\cdot 10^{-4}$. With this value of mass gap, the domain
of distances in Fig. 3 is $d\in[0,100d_{c}]$. The energy
$\mathcal{E}_{\mathsf{n}}$ is lowered by $16$ times in this interval. For
small separation, we observe the $1/d^{3}$ dependence for the Casimir energy
$E_{\mathsf{n}}$ as was proved above. For large separation, we observe
$1/d^{5}$ dependence of the energy $E_{\mathsf{n}}$.

The $\mathsf{n}$-dependence of $\mathcal{E}_{\mathsf{n}}$ is plotted in Fig. 4
for different interplane distances. The energy decreases in accuracy with
increasing interplane distance, as should be the case (the insert in Fig. 4).

Figure 4: The Casimir energy per unit graphene $\mathcal{E}_{\mathsf{n}}$ as a
function of the number of planes in the stack. Here, $m=0.1$ eV,
$d=d_{c},10d_{c},20d_{c}$ from the top down, where $d_{c}=0.3345$ is
interplane distance in graphite. The insert shows the energy for a stack of an
infinite number of graphene as a function of interplane distance up to
$50d_{c}$.

It is useful to compare these calculations with the case of the constant
conductivity case with universal conductivity of graphene
$\eta_{\mathrm{gr}}$, which was considered in detail in Ref. [29]. The
constant conductivity leads to $1/d^{3}$ dependence and therefore
$\mathcal{E}_{\mathsf{n}}$ does not depend on the interplane distance. The
better agreement appears for small distances, where for the graphene case the
$\mathcal{E}_{\mathsf{n}}$ is approximately constant (see Fig. 3). The case of
constant conductivity gives a greater value $\sim 10\%$.

Let us estimate the binding energy due to the Casimir energy for graphite. For
graphite interplane separation $d_{c}=0.3345$ nm, we obtain the following
energy per unit plane in an infinite stack of graphene
$\overline{E}_{\infty}=59.23\ \textrm{erg}/\textrm{cm}^{2}$. The binding
energy $E_{ib}=\overline{E}_{\infty}/d_{c}\rho_{c}$, where $\rho_{c}=2.23\
\textrm{g}/\textrm{cm}^{3}$ is the graphite mass density. For these values, we
obtain $E_{ib}=9.9$ meV/atom, which is by $10$% smaller than for constant
conductivity case [29]. From the first principle, the cohesion energies are
$24-26$ meV/atom [42] and $24$ meV/atom [43]. The experimental data gives
cohesion $35\pm 10,\ 15$ meV/atom [44] and $61\pm 5$ meV/atom [45]. Most
likely, the Casimir energy gives an essential contribution to the binding
energy.

The numerical evaluations of force acting on the graphene $s$ in the stack of
$\mathsf{n}$ graphene are plotted in Figs. 5, 6. We evaluate the function
$\mathcal{F}_{s,\mathsf{n}}=F_{s,\mathsf{n}}/F_{Cas}$, where
$F_{Cas}=-\pi^{2}/240d^{4}$ is the Casimir force for two ideal planes.

Figure 5: The Casimir force $\mathcal{F}_{s,\mathsf{n}}$ acting on the
graphene $s$ in the stack of $10$ and $20$ planes.

From Fig. 5, we observe that the value of the force falls very quickly,
beginning at the first graphene sheet. The force acting on the second graphene
is in $35$ time is smaller than the force acting on the first one. Fig. 6
shows that the force very quickly becomes the force for infinite number of
graphene – already the stack of $10$ graphene gives the same force as for an
infinite stack.

Figure 6: The Casimir force $\mathcal{F}_{s,\mathsf{n}}$ acting on the first
and second (insert figure) graphene in the stack of $\mathsf{n}$ graphene.

##  V Conclusion

We considered here the Casimir energy for a stack of conductive planes and a
force acting on the plane on this stack. To calculate the Casimir force, the
scattering theory is used. Starting from the scattering on the single plane,
we found the scattering matrix for the set of planes and found recurrent
relation (29) for the Casimir energy (22). In the particular case of a stack
of graphene with equal interplane distance, the recurrent relation (36) may be
solved in the manifest form (47). The reflection Fresnel matrix may be
diagonalized with eigenvalues corresponding to the TE and TM electromagnetic
modes. The resulting expression for the Casimir energy is the sum of
contributions of these modes (63). The function $D_{\mathsf{n}}$ has the same
form (62) which was found in Ref. [29] for the constant conductivity case. We
have to use, only the corresponding reflection coefficients for graphene (54).

For small interplane separations $md\ll 1$, the Casimir energy is divergent as
$1/d^{3}$ (without mass gap parameter $m$). In the opposite limit of large
distances, the energy falls as $1/(m^{2}d^{5})$ in accordance to Ref. [7]. The
binding energy of infinite stack graphene with interplane separation of
graphite, $d_{c}$, is $E_{ib}=9.9$ meV/atom, which is by $10$% smaller than
for constant conductivity case considered in Ref. [29]. The obtained value of
binding energy depends only weakly on the mass gap parameter because the case
of graphite interplane separation corresponds to the case $md_{c}\ll 1$ (see
the insert in Fig. 3).

The Casimir force acting on the graphene $s$ in the stack of $\mathsf{n}$
graphene is a sum of two contributions (73), where function $G_{s,\mathsf{n}}$
has the same form (74) as for constant conductivity case found in Ref. [29]
with corresponding reflections coefficients for TE and TM modes. The force
acting on graphene falls very quickly – for the second plane it is $35$ times
smaller than for the first graphene (see Fig. 5). The extremes graphenes
squeezes the stack. The force acting on graphene in a stack is quick,
$\mathsf{n}\sim 10$ (see Fig. 6), becomes the same as for an infinite stack.

The influence of non-zero temperature and chemical potential can be taken into
account by employing the corresponding conductivity obtained in Ref. [35,
*Bordag:2016:ECefdg]. We will analyze this influence in forthcoming papers.

###### Acknowledgements.

NK was supported in part by the grants 2022/08771-5, 2021/10128-0 of São Paulo
Research Foundation (FAPESP).

##  Appendix A Some relations

I) First of all we prove the formula (19). By using relations (16), (17) and
(14) one has

$\displaystyle\det\bm{\mathcal{S}}$
$\displaystyle=\frac{\det\mathbf{R}}{\det\accentset{*}{\mathbf{R}}^{\prime}}=\frac{\det\mathbf{R}}{\det(\mathbf{R}^{\prime}-\mathbf{R}\mathbf{R}^{-1}\mathbf{T}^{\prime})^{-1}}\allowdisplaybreaks$
$\displaystyle=\frac{\det\left[-\mathbf{K}^{-1}_{22}\mathbf{K}_{21}\right]}{\det\left[\mathbf{K}_{12}\mathbf{K}^{-1}_{22}+\left(\mathbf{K}_{11}-\mathbf{K}_{12}\mathbf{K}_{22}^{-1}\mathbf{K}_{21}\right)\mathbf{K}_{21}^{-1}\mathbf{K}_{22}\mathbf{K}_{22}^{-1}\right]^{-1}}\allowdisplaybreaks$
$\displaystyle=\frac{\det\left[-\mathbf{K}^{-1}_{22}\mathbf{K}_{21}\right]}{\det\left[\mathbf{K}_{11}\mathbf{K}_{21}^{-1}\right]^{-1}}=\frac{\det\left[-\mathbf{K}^{-1}_{22}\mathbf{K}_{21}\right]}{\det\left[\mathbf{K}_{21}\mathbf{K}_{11}^{-1}\right]}=\frac{\det\mathbf{K}_{11}}{\det\mathbf{K}_{22}}.$
(76)

II) The different forms of the matrices $\mathbf{B}_{i}$.

By direct solving the relation $\mathbf{C}^{-1}_{i}\mathbf{C}_{i}=\mathbf{I}$,
we obtain

$\mathbf{C}^{-1}_{i}=\begin{pmatrix}\bm{r}^{\prime}_{i}\bm{t}^{\prime-1}_{i}&-\mathbf{I}\\\
\bm{t}^{\prime-1}_{i}&0\end{pmatrix}.$ (77)

Therefore, we can represent the matrix $\mathbf{B}_{i}$ in different forms

$\mathbf{B}_{i}=\begin{pmatrix}-\accentset{*}{\bm{t}}^{\prime-1}_{i}&-\bm{r}^{\prime}_{i}\bm{t}^{\prime-1}_{i}\\\
\bm{t}^{\prime-1}_{i}\bm{r}_{i}&-\bm{t}^{\prime-1}_{i}\end{pmatrix}=\begin{pmatrix}-\accentset{*}{\bm{t}}^{\prime-1}_{i}&\accentset{*}{\bm{t}}^{\prime-1}_{i}\accentset{*}{\bm{r}}_{i}\\\
\accentset{*}{\bm{r}}^{\prime}_{i}\accentset{*}{\bm{t}}^{\prime-1}_{i}&-\bm{t}^{\prime-1}_{i}\end{pmatrix}=\begin{pmatrix}-\accentset{*}{\bm{t}}^{\prime-1}_{i}&\accentset{*}{\bm{t}}^{\prime-1}_{i}\accentset{*}{\bm{r}}^{\prime}_{i}\\\
\bm{t}^{\prime-1}_{i}\bm{r}_{i}&-\bm{t}^{\prime-1}_{i}\end{pmatrix},$ (78)

taking into account relations (16) and (17). Therefore, the matrices have the
following structure

$\mathbf{B}_{i}=\begin{pmatrix}\bm{b}^{i}_{11}&\bm{b}^{i}_{12}\\\
\accentset{*}{\bm{b}}^{i}_{12}&\accentset{*}{\bm{b}}^{i}_{11}\end{pmatrix}=\begin{pmatrix}\accentset{*}{\bm{b}}^{i}_{22}&\accentset{*}{\bm{b}}^{i}_{21}\\\
\bm{b}^{i}_{21}&\bm{b}^{i}_{22}\end{pmatrix}.$ (79)

With replacement (4) the dependence of the position appears

$\mathbf{B}_{i}=\begin{pmatrix}\bm{b}^{i}_{11}&\bm{b}^{i}_{12}e^{-2\mathrm{i}k_{z}z_{i}}\\\
\bm{b}^{i}_{21}e^{2\mathrm{i}k_{z}z_{i}}&\bm{b}^{i}_{22}\end{pmatrix}.$ (80)

Some properties:

$\mathbf{B}_{i}^{-1}=\begin{pmatrix}-\bm{t}^{-1}_{i}&\bm{t}^{-1}_{i}\bm{r}^{\prime}_{i}\\\
-\bm{r}_{i}\bm{t}^{-1}_{i}&-\accentset{*}{\bm{t}}^{-1}_{i}\end{pmatrix},\
\det\mathbf{B}_{i}=\frac{\det\bm{t}_{i}}{\det\bm{t}^{\prime}_{i}}.$ (81)

Therefore, due to relations (79) the matrix $\mathbf{K}$ has the following
structure

$\mathbf{K}=\begin{pmatrix}\mathbf{K}_{11}&\mathbf{K}_{12}\\\
\accentset{*}{\mathbf{K}}_{12}&\accentset{*}{\mathbf{K}}_{11}\end{pmatrix}=\begin{pmatrix}\accentset{*}{\mathbf{K}}_{22}&\accentset{*}{\mathbf{K}}_{21}\\\
\mathbf{K}_{21}&\mathbf{K}_{22}\end{pmatrix}.$ (82)

The matrices $\mathbf{B}_{i}$ maybe represented in the following forms

$\mathbf{B}_{i}=\bm{\tau}_{i}\bm{\rho}_{i}=\bm{\rho}^{\prime}_{i}\bm{\tau}_{i},$
(83)

where

$\bm{\rho}_{i}=\begin{pmatrix}-\mathbf{I}&\accentset{*}{\bm{r}}_{i}\\\
\bm{r}_{i}&-\mathbf{I}\end{pmatrix},\bm{\rho}^{\prime}_{i}=\begin{pmatrix}-\mathbf{I}&-\bm{r}^{\prime}_{i}\\\
-\accentset{*}{\bm{r}}^{\prime}_{i}&-\mathbf{I}\end{pmatrix},\bm{\tau}_{i}=\begin{pmatrix}\accentset{*}{\bm{t}}^{\prime-1}_{i}&0\\\
0&\bm{t}^{\prime-1}_{i}\end{pmatrix}.$ (84)

It easy to see that

$\bm{\rho}_{i+1}\bm{\rho}^{\prime}_{i}=\begin{pmatrix}\mathbf{I}-\accentset{*}{\bm{r}}_{i+1}\accentset{*}{\bm{r}}^{\prime}_{i}&\bm{r}^{\prime}_{i}-\accentset{*}{\bm{r}}_{i+1}\\\
\accentset{*}{\bm{r}}^{\prime}_{i}-\bm{r}_{i+1}&\mathbf{I}-\bm{r}_{i+1}\bm{r}^{\prime}_{i}\end{pmatrix}.$
(85)

##  Appendix B Different number of planes

###  B.1 The case without planes, $\mathsf{n}=0$

Without planes there is no scattering and $\mathbf{K}_{0}=\mathbf{I}_{4}$ is
diagonal matrix $4\times 4$. We define $\bm{\mathcal{D}}_{0}=\mathbf{I}$.

###  B.2 Single plane, $\mathsf{n}=1$

The matrix $\mathbf{K}_{1}=\mathbf{B}_{1}$ has the following form (see Eq.
(78))

$\mathbf{K}_{1}=\begin{pmatrix}-\accentset{*}{\bm{t}}^{\prime-1}_{1}&-\bm{r}^{\prime}_{1}\bm{t}^{\prime-1}_{1}\\\
\bm{t}^{\prime-1}_{1}\bm{r}_{1}&-\bm{t}^{\prime-1}_{1}\end{pmatrix}.$ (86)

Therefore,

$\bm{\mathcal{D}}_{1}=-\mathbf{K}_{1}^{22}\bm{t}^{\prime}_{1}=\mathbf{I}.$
(87)

This expression does not depend on the position of plane and gives no
contribution to the Casimir energy.

###  B.3 Two planes, $\mathsf{n}=2$

In this case

$\mathbf{K}_{2}=\mathbf{B}_{2}\mathbf{B}_{1}=\bm{\tau}_{2}\begin{pmatrix}\mathbf{I}-\accentset{*}{\bm{r}}_{2}\accentset{*}{\bm{r}}^{\prime}_{1}&\bm{r}^{\prime}_{1}-\accentset{*}{\bm{r}}_{2}\\\
\accentset{*}{\bm{r}}^{\prime}_{1}-\bm{r}_{2}&\mathbf{I}-\bm{r}_{2}\bm{r}^{\prime}_{1}\end{pmatrix}\bm{\tau}_{1}.$
(88)

Therefore,

$\bm{\mathcal{D}}_{2}=\mathbf{I}-\bm{r}_{2}\bm{r}^{\prime}_{1}.$ (89)

###  B.4 Three planes, $\mathsf{n}=3$

For three planes one has

$\mathbf{K}_{3}=\mathbf{B}_{3}\mathbf{B}_{2}\mathbf{B}_{1}=-\bm{\tau}_{3}\begin{pmatrix}\mathbf{I}-\accentset{*}{\bm{r}}_{3}\accentset{*}{\bm{r}}^{\prime}_{2}&\bm{r}^{\prime}_{2}-\accentset{*}{\bm{r}}_{3}\\\
\accentset{*}{\bm{r}}^{\prime}_{2}-\bm{r}_{3}&\mathbf{I}-\bm{r}_{3}\bm{r}^{\prime}_{2}\end{pmatrix}\begin{pmatrix}\accentset{*}{\bm{t}}^{\prime-1}_{2}&\accentset{*}{\bm{t}}^{\prime-1}_{2}\bm{r}^{\prime}_{1}\\\
\bm{t}^{\prime-1}_{2}\accentset{*}{\bm{r}}^{\prime}_{1}&\bm{t}^{\prime-1}_{2}\end{pmatrix}\bm{\tau}_{1}.$
(90)

Then, by using the relations

$\accentset{*}{\bm{r}}^{\prime}\accentset{*}{\bm{t}}^{\prime-1}=-\bm{t}^{\prime-1}\bm{r},\
\accentset{*}{\bm{t}}^{\prime-1}\accentset{*}{\bm{r}}=-\bm{r}^{\prime}\bm{t}^{\prime-1},$
(91)

and Eqs. (16), (17) we obtain

$\bm{\mathcal{D}}_{3}=\bm{\mathcal{D}}_{2}^{(32)}\bm{t}^{\prime-1}_{2}\bm{\mathcal{D}}_{2}^{(21)}\bm{t}^{\prime}_{2}-\bm{r}_{3}\bm{t}_{2}\bm{r}^{\prime}_{1}\bm{t}^{\prime}_{2},$
(92)

where $\bm{\mathcal{D}}_{2}^{(ij)}=\mathbf{I}-\bm{r}_{i}\bm{r}^{\prime}_{j}$
is for two planes $i$ and $j$. In manifest form

$\bm{\mathcal{D}}_{3}=\mathbf{I}-\bm{r}_{3}\bm{r}^{\prime}_{2}-\bm{t}^{\prime-1}_{2}\bm{r}_{2}\bm{r}^{\prime}_{1}\bm{t}^{\prime}_{2}+\bm{r}_{3}\bm{r}^{\prime}_{2}\bm{t}^{\prime-1}_{2}\bm{r}_{2}\bm{r}^{\prime}_{1}\bm{t}^{\prime}_{2}-\bm{r}_{3}\bm{t}_{2}\bm{r}^{\prime}_{1}\bm{t}^{\prime}_{2}.$
(93)

By using replacement (4) we observe that the second term $\sim
e^{2\mathrm{i}k_{z}d_{3,2}}$, the third term $\sim
e^{2\mathrm{i}k_{z}d_{2,1}}$, and the last two terms $\sim
e^{2\mathrm{i}k_{z}d_{3,1}}$, where $d_{i,k}$ is distance between planes $i$
and $k$.

A comment is in order. The first term in expression obtained (92) gives
directly the sum of energies of pars. Indeed,

$\ln\det(\bm{\mathcal{D}}_{2}^{(32)}\bm{t}^{\prime-1}_{2}\bm{\mathcal{D}}_{2}^{(21)}\bm{t}^{\prime}_{2})=\ln\det\bm{\mathcal{D}}_{2}^{(32)}+\ln\det\bm{\mathcal{D}}_{2}^{(21)}.$
(94)

The second term breaks the additivity of the Casimir energy. It contains
refraction coefficients the first and last plane without middle plane.

###  B.5 Four planes

We have

$\displaystyle\mathbf{K}_{4}$
$\displaystyle=\bm{\tau}_{4}\begin{pmatrix}\mathbf{I}-\accentset{*}{\bm{r}}_{4}\accentset{*}{\bm{r}}^{\prime}_{3}&\bm{r}^{\prime}_{3}-\accentset{*}{\bm{r}}_{4}\\\
\accentset{*}{\bm{r}}^{\prime}_{3}-\bm{r}_{4}&\mathbf{I}-\bm{r}_{4}\bm{r}^{\prime}_{3}\end{pmatrix}\bm{\tau}_{3}\bm{\tau}_{2}\begin{pmatrix}\mathbf{I}-\accentset{*}{\bm{r}}_{2}\accentset{*}{\bm{r}}^{\prime}_{1}&\bm{r}^{\prime}_{1}-\accentset{*}{\bm{r}}_{2}\\\
\accentset{*}{\bm{r}}^{\prime}_{1}-\bm{r}_{2}&\mathbf{I}-\bm{r}_{2}\bm{r}^{\prime}_{1}\end{pmatrix}\bm{\tau}_{1}.$
(95)

Let us express now this expression to that without asterisk. We use relations
(91). Then,

$\displaystyle\bm{\mathcal{D}}_{4}$
$\displaystyle=\bm{\mathcal{D}}_{2}^{(43)}\bm{t}^{\prime-1}_{3}\bm{\mathcal{D}}_{2}^{(32)}\bm{t}^{\prime-1}_{2}\bm{\mathcal{D}}_{2}^{(21)}\bm{t}^{\prime}_{2}\bm{t}^{\prime}_{3}-\bm{\mathcal{D}}_{2}^{(43)}\bm{t}^{\prime-1}_{3}\bm{r}_{3}\bm{t}_{2}\bm{r}^{\prime}_{1}\bm{t}^{\prime}_{2}\bm{t}^{\prime}_{3}\allowdisplaybreaks$
$\displaystyle-\bm{r}_{4}\bm{t}_{3}\bm{r}^{\prime}_{2}\bm{t}^{\prime-1}_{2}\bm{\mathcal{D}}_{2}^{(21)}\bm{t}^{\prime}_{2}\bm{t}^{\prime}_{3}-\bm{r}_{4}\bm{t}_{3}\bm{t}_{2}\bm{r}^{\prime}_{1}\bm{t}^{\prime}_{2}\bm{t}^{\prime}_{3}.$
(96)

The first term describes additive part of the Casimir energy and rest terms
give non-additivity of the energy.

## References

  * Casimir [1948] H. B. G. Casimir, Kon. Ned. Akad. Wetensch. Proc. 51, 793–795 (1948).
  * Bordag  _et al._ [2009a] M. Bordag, G. L. Klimchitskaya, U. Mohideen, and V. M. Mostepanenko,  _Advances in the Casimir Effect_ (2009) p. 768.
  * Milton [2001] K. A. Milton,  _Casimir Effect: Physical Manifestation of Zero-Point Energy._ (World Scientific, Singapore, 2001) p. 301.
  * Parsegian [2006] A. V. Parsegian,  _Van der Waals Forces. A Handbook for Biologists, Chemists, Engineers, and Physicists_ (Cambridge University Press, 2006) p. 380.
  * Milonni [1994] P. W. Milonni,  _The Quantum Vacuum. An Introduction to Quantum Electrodynamics._ (Academic Press, New York, 1994) p. 522.
  * Woods  _et al._ [2016] L. M. Woods, D. A. R. Dalvit, A. Tkatchenko, P. Rodriguez-Lopez, A. W. Rodriguez, and R. Podgornik, Rev. Mod. Phys. 88, 45003 (2016).
  * Bordag  _et al._ [2009b] M. Bordag, I. V. Fialkovsky, D. M. Gitman, and D. V. Vassilevich, Phys. Rev. B 80, 245406 (2009b).
  * Grushin and Cortijo [2011] A. G. Grushin and A. Cortijo, Phys. Rev. Lett. 106, 20403 (2011).
  * Lu [2018] B.-S. Lu, Phys. Rev. B 97, 045427 (2018).
  * Zhao  _et al._ [2009] R. Zhao, J. Zhou, T. Koschny, E. N. Economou, and C. M. Soukoulis, Phys. Rev. Lett. 103, 103602 (2009).
  * Wilson  _et al._ [2015] J. H. Wilson, A. A. Allocca, and V. Galitski, Phys. Rev. B 91, 235115 (2015).
  * Rodriguez-Lopez  _et al._ [2020] P. Rodriguez-Lopez, A. Popescu, I. Fialkovsky, N. Khusnutdinov, and L. M. Woods, Commun. Mater. 1, 14 (2020).
  * Sheehan [2009] D. P. Sheehan, J. Chem. Phys 131, 104706 (2009).
  * Machta  _et al._ [2012] B. B. Machta, S. L. Veatch, and J. P. Sethna, Phys. Rev. Let. 109, 138101 (2012).
  * Pawlowski and Zielenkiewicz [2013] P. H. Pawlowski and P. Zielenkiewicz, J. Membrane Biol. 246, 383–389 (2013).
  * Ninham and Parsegian [1970a] B. W. Ninham and V. A. Parsegian, J. Chem. Phys. 53, 3398–3402 (1970a).
  * Ninham and Parsegian [1970b] B. W. Ninham and V. A. Parsegian, J. Chem. Phys. 52, 4578 (1970b).
  * Tomaš [2002] M. S. Tomaš, Phys. Rev. A 66, 052103 (2002).
  * Tomaš [2012] M.-S. Tomaš, Int. J. Mod. Phys.: Conf. Ser. 14, 561–565 (2012).
  * Teo [2010] L. P. Teo, Phys. Rev. A 81, 032502 (2010).
  * Ellingsen and Brevik [2007] S. A. Ellingsen and I. Brevik, J. Phys. A: Math. Theor. 40, 3643–3664 (2007).
  * Kheirandish  _et al._ [2011] F. Kheirandish, M. Soltani, and J. Sarabadani, Ann. Phys. 326, 657–667 (2011).
  * Amooghorban  _et al._ [2011] E. Amooghorban, M. Wubs, N. A. Mortensen, and F. Kheirandish, Phys. Rev. A 84, 013806 (2011).
  * Sernelius [2014] B. E. Sernelius, Phys. Rev. B 90, 155457 (2014).
  * Davids  _et al._ [2010] P. S. Davids, F. Intravaia, F. S. S. Rosa, and D. A. R. Dalvit, Phys. Rev. A 82, 062111 (2010).
  * Abbas  _et al._ [2017] C. Abbas, B. Guizal, and M. Antezza, Phys. Rev. Lett. 118, 126101 (2017).
  * Tatur and Woods [2008] K. Tatur and L. M. Woods, Phys. Lett. Sect. A Gen. At. Solid State Phys. 372, 6705–6710 (2008).
  * Khusnutdinov and Kashapov [2015] N. Khusnutdinov and R. Kashapov, Theor. Math. Phys. 183, 491 (2015).
  * Khusnutdinov  _et al._ [2015] N. Khusnutdinov, R. Kashapov, and L. M. Woods, Phys. Rev. D 92, 045002 (2015).
  * Khusnutdinov  _et al._ [2016] N. Khusnutdinov, R. Kashapov, and L. M. Woods, Phys. Rev. A 94, 012513 (2016).
  * Kashapov  _et al._ [2016] R. Kashapov, N. Khusnutdinov, and L. M. Woods, Int. J. Mod. Phys. A 31, 1641028 (2016).
  * Khusnutdinov  _et al._ [2018] N. Khusnutdinov, R. Kashapov, and L. M. Woods, 2D Mater. 5, 35032 (2018).
  * Khusnutdinov and Woods [2019] N. Khusnutdinov and L. M. Woods, JETP Lett. 110, 183 (2019).
  * Djurišić and Li [1999] A. B. Djurišić and E. H. Li, J. Appl. Phys. 85, 7404–7410 (1999).
  * Fialkovsky  _et al._ [2011] I. V. Fialkovsky, V. N. Marachevsky, and D. V. Vassilevich, Phys. Rev. B 84, 035446 (2011).
  * Bordag  _et al._ [2016] M. Bordag, I. Fialkovskiy, and D. Vassilevich, Phys. Rev. B 93, 075414 (2016).
  * Bordag [1995] M. Bordag, J. Phys. A. Math. Gen. 28, 755–765 (1995).
  * Reynaud  _et al._ [2010] S. Reynaud, A. Canaguier-Durand, R. Messina, A. Lambrecht, and P. A. Maia Neto, Int. J. Mod. Phys. A 25, 2201–2211 (2010).
  * Ingold and Lambrecht [2015] G.-L. Ingold and A. Lambrecht, Am. J. Phys. 83, 156–162 (2015).
  * Fialkovsky  _et al._ [2018] I. Fialkovsky, N. Khusnutdinov, and D. Vassilevich, Phys. Rev. B 97, 165432 (2018).
  * Khusnutdinov and Emelianova [2019] N. Khusnutdinov and N. Emelianova, Int. J. Mod. Phys. A 34, 1950008 (2019).
  * Schabel and Martins [1992] M. C. Schabel and J. L. Martins, Phys. Rev. B 46, 7185–7188 (1992).
  * Rydberg  _et al._ [2003] H. Rydberg, M. Dion, N. Jacobson, E. Schroder, P. Hyldgaard, S. I. Simak, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 91, 126402 (2003).
  * Benedict  _et al._ [1998] L. X. Benedict, N. G. Chopra, M. L. Cohen, A. Zettl, S. G. Louie, and V. H. Crespi, Chem. Phys. Lett. 286, 490–496 (1998).
  * Zacharia  _et al._ [2004] R. Zacharia, H. Ulbricht, and T. Hertel, Phys. Rev. B 69, 155406 (2004).


# Fit-free determination of scale invariant equations of state:
application to the 2D Bose gas across the Berezinksii–Kosterlitz–Thouless
transition

Rémi Desbuquois  Institute for Quantum Electronics, ETH Zurich, 8093 Zurich,
Switzerland  Tarik Yefsah  MIT-Harvard Center for Ultracold Atoms, Research
Laboratory of Electronics, Department of Physics, Massachusetts Institute of
Technology, Cambridge, Massachusetts 02139, USA  Lauriane Chomaz  Laboratoire
Kastler Brossel, CNRS, UPMC, ENS, 24 rue Lhomond, F-75005 Paris, France
Collège de France, 11 Place Marcelin Berthelot, 75231 Paris Cedex 05, France
Christof Weitenberg  Institut für Laserphysik, Universität Hamburg, Luruper
Chaussee 149, D-22761 Hamburg, Germany  Laura Corman  Laboratoire Kastler
Brossel, CNRS, UPMC, ENS, 24 rue Lhomond, F-75005 Paris, France  Sylvain
Nascimbène  Laboratoire Kastler Brossel, CNRS, UPMC, ENS, 24 rue Lhomond,
F-75005 Paris, France  Jean Dalibard  Laboratoire Kastler Brossel, CNRS, UPMC,
ENS, 24 rue Lhomond, F-75005 Paris, France  Collège de France, 11 Place
Marcelin Berthelot, 75231 Paris Cedex 05, France

###### Abstract

We present a general “fit-free” method for measuring the equation of state
(EoS) of a scale-invariant gas. This method, which is inspired from the
procedure introduced by Ku et al. [Science 335, 563 (2012)] for the unitary
three-dimensional Fermi gas, provides a general formalism which can be readily
applied to any quantum gas in a known trapping potential, in the frame of the
local density approximation. We implement this method on a weakly-interacting
two-dimensional Bose gas in the vicinity of the Berezinskii-Kosterlitz-
Thouless transition, and determine its EoS with unprecedented accuracy in the
critical region. Our measurements provide an important experimental benchmark
for classical field approaches which are believed to accurately describe
quantum systems in the weakly interacting but non-perturbative regime.

###### pacs:

03.75.-Hh, 05.70.Ce, 67.85.-d

Homogeneous matter at thermal equilibrium is described by an equation of state
(EoS), _i.e.,_ a functional relation between thermodynamic variables of the
system. While the EoS is analytically known for ideal gases, one must resort
to approximations or numerical calculations to determine the EoS of
interacting fluids, which must then be confronted to experiments. Thanks to a
precise control of temperature, confining potential and interaction strength,
cold atomic gases constitute a system of choice for the experimental
determination of quantum matter EoS Ho and Zhou (2009). While performed on
atomic systems, such measurements often provide crucial insight on generic
physical problems, well beyond the atomic physics perspective. Prominent
examples are the recent measurements of the EoS of atomic Fermi gases
Nascimbène et al. (2010); Navon et al. (2010); Ku et al. (2012); Van Houcke et
al. (2012), which provided a precious quantitative support for our
understanding of strongly interacting fermions at low temperature. Another
important paradigm accessible to atomic gases is found in two-dimensional
quantum systems, where the low temperature state is established via a defect-
driven transition. This generic phenomenon of two-dimensional systems is
described by the celebrated Berezinskii-Kosterlitz-Thouless (BKT) theory, with
a scope that ranges from superconductivity to quantum Hall bilayer physics to
high energy physics.

In this context, the weakly interacting two-dimensional Bose gas is of
particular interest as it supports the fundamental principles of the BKT
theory, while allowing for a simplified theoretical description. Indeed, for
small enough inter-particle interactions, the dynamics of the two-dimensional
Bose gas is well captured by a classical field model Prokof’ev et al. (2001);
Prokof’ev and Svistunov (2002), which is itself described by a dimensionless
coupling constant and exhibits scale-invariance 111Beyond the classical field
level the scale invariance is broken Vogt et al. (2012), but the deviation can
be neglected for small enough interaction strengths Olshanii et al. (2010);
Langmack et al. (2012). In general, scale invariance occurs in any fluid where
no explicit energy/length scale is associated to the (binary) interaction
potential. For the weakly interacting two-dimensional case, the length scale
introduced by the 3D scattering length is normalized by the extension of the
system in the third dimension, and this dimensionless ratio characterizes the
effective 2D interaction strength. Scale invariance is also found in the
unitary Fermi gas, where the scattering length describing $s-$wave
interactions diverges (for a review, see Giorgini et al. (2008)). This
property considerably simplifies the EoS structure, as general dimensionless
quantities such as the phase space density $\mathcal{D}$ which usually depend
separately on the chemical potential $\mu$ and the temperature $T$ can only be
expressed as the ratio $\mu/k_{\rm B}\,T$ owing to the absence of other energy
scales Hung et al. (2011); Yefsah et al. (2011).

The usual method for determining the EoS of a cold atomic gas starts with the
measurement of the density distribution $n(\bm{r})$ in a smoothly varying,
confining potential $V(\bm{r})$. The gas is assumed to be in thermal
equilibrium, and its chemical potential at the trap center is $\mu$. Using the
Local Density Approximation (LDA), the measured ${\cal D}(\bm{r})$ is linked
to that of a uniform fluid with the same interaction strength, the same
temperature, and with chemical potential $\mu(\bm{r})=\mu-V(\bm{r})$ Ho and
Zhou (2009). The values of $T$ and $\mu$ for a given realization of the gas
are obtained by comparing the low-density wings of $n(\bm{r})$ with the known
theoretical result for a dilute fluid. A major drawback of this method is the
fact that any systematic error in the determination of the density, _e.g._ due
to imperfect calibration of the probing system, will lead to inaccurate values
of $\mu$ and $T$, and thus affect the measurement of the EoS ${\cal
D}(\mu,T)$. Recently, an alternative fit-free method that does not suffer from
this limitation has been put forward in Ku et al. (2012) for the measurement
of the EoS of the scale-invariant 3D Fermi gas. It is based on the use of two
specific thermodynamic variables, pressure and compressibility; in addition,
absolute energy scales $T$ and $\mu$ were replaced by a single relative scale
${\rm d}\mu$, which was itself determined by the LDA through ${\rm d}\mu=-{\rm
d}V$.

The purpose of the present Letter is two-fold. First, we describe a method
that generalizes the procedure introduced in Ku et al. (2012), which does not
rely on specific thermodynamic variables but rather provides a generic
formalism that can readily be applied to other quantum systems. Second, we
implement this method on a two-dimensional (Bose) fluid, for which all
information about the spatial density $n(\bm{r})$ is directly accessible from
an image of the cloud. This is to be contrasted with the 3D case studied in Ku
et al. (2012), where one has to resort to an Abel transform to reconstruct the
spatial density profile, before implementing the fit-free determination of the
EoS. The EoS reconstructed here reaches a precision which makes it suitable
for a quantitative comparison to the classical-field Monte Carlo calculation
Prokof’ev and Svistunov (2002) in the critical region, which could not be
conclusive from previous measurements Hung et al. (2011); Yefsah et al.
(2011). Our measurement, with a relative statistical error smaller than 1% on
the detectivity, is found in excellent agreement (better than 5%) with the
prediction obtained from Prokof’ev and Svistunov (2002) at the critical point
and deeper in the superfluid regime. In the normal regime close to the
transition point, we observe a deviation on the order of 15%, which might be
due to beyond classical-field effects.

We start our analysis by considering an atomic gas in thermodynamic
equilibrium confined in a known potential $V(\bm{r})$. The only hypothesis for
the method is the LDA, which entails that $n(\bm{r})$ depends on position only
through the local value of the trapping potential: $n(\bm{r})=n[V(\bm{r})]$.
Although this method is applicable to any dimension, we focus here on the
particular case of the two-dimensional gas for the sake of clarity. Let us
introduce the energy $E[V(\bm{r}))]$ with $\bm{r}=(x,y)$, defined by 222In
dimension $d$, this equation takes the general form $E=\hbar^{2}\,n^{2/d}/m$

$E=\frac{\hbar^{2}}{m}n,$ (1)

which we want to combine with other relevant energies in order to form useful
dimensionless variables. Though no absolute energy scales are readily
available, a relative energy scale is provided by the variation of the
trapping potential ${\rm d}V$. Furthermore, quantities formed in this manner
are directly connected by the LDA to the properties of the uniform gas using
the relation ${\rm d}\mu=-{\rm d}V$. Thus, we define the dimensionless
quantities

$X_{\nu}\equiv
E^{\nu-1}\frac{\partial^{\nu}E}{\partial\mu^{\nu}}=(-1)^{\nu}\,E^{\nu-1}\frac{{\rm
d}^{\nu}E}{{\rm d}V^{\nu}}\,,$ (2)

where $\nu$ is an integer. By convention, $X_{0}=1$ and a negative $\nu$ will
instead correspond to $|\nu|$ successive integrations of $E$ with respect to
$V$, with for example

$X_{-1}=\frac{1}{E^{2}}\int^{\infty}_{V}E(V^{\prime})\;{\rm d}V^{\prime}\,.$
(3)

From a given image of the gas $n(\bm{r})$, one can thus construct all
functions $X_{\nu}(V)$. In the case of a scale invariant system, the knowledge
of a single thermodynamic variable $X_{\nu}$ is sufficient to determine the
state of the fluid, hence the values of all other variables
$X_{\nu^{\prime}}$. In other words, all individual measurements must collapse
on a single line in each plane $\\{X_{\nu},\,X_{\nu^{\prime}}\\}$,
irrespective of their temperature and chemical potential. Such a line is a
valid EoS of the fluid under consideration.

Once the $X_{\nu}$ are known, all other thermodynamic quantities can be
determined, up to an integration constant. In particular, one can derive the
phase-space density $\mathcal{D}$ and the ratio $\alpha=\mu/k_{\rm B}T$. Let
us suppose that a point $(X_{\nu}^{(0)},\,X_{\nu^{\prime}}^{(0)})$ can be
identified in a known portion the EoS, and that it corresponds to the values
$\alpha_{0}$ and $\mathcal{D}_{0}$. The link between the set $\\{X_{\nu}\\}$
and $(\alpha,\mathcal{D})$ is provided by

$\displaystyle\mathcal{D}(X_{\nu}^{(1)})=\mathcal{D}_{0}\exp\left(\int_{X_{\nu}^{(0)}}^{X_{\nu}^{(1)}}\frac{X_{1}}{(\nu-1)X_{1}X_{\nu}+X_{\nu+1}}{\rm
d}X_{\nu}\right),$ (4)
$\displaystyle\alpha(X_{\nu}^{(1)})=\alpha_{0}+\frac{1}{2\pi}\int_{X_{\nu}^{(0)}}^{X_{\nu}^{(1)}}\frac{\mathcal{D}(X_{\nu})}{(\nu-1)X_{1}X_{\nu}+X_{\nu+1}}{\rm
d}X_{\nu}\,.$ (5)

The determination of $(\alpha,\mathcal{D})$ thus requires the knowledge of a
triplet $\\{X_{1},X_{\nu},X_{\nu+1}\\}$. This requirement can be weakened by
choosing $\nu=1$ or $\nu=-1$, in which case only the pairs
$\\{X_{1},X_{-1}\\}$ or $\\{X_{1},X_{2}\\}$ are needed, and
$(\alpha,\mathcal{D})$ are given by

$\displaystyle\mathcal{D}(X_{-1}^{(1)})=\mathcal{D}_{0}\exp\left(\int_{X_{-1}^{(0)}}^{X_{-1}^{(1)}}\frac{X_{1}}{1-2X_{1}X_{-1}}{\rm
d}X_{-1}\right)\,,$ (6)
$\displaystyle\alpha(X_{-1}^{(1)})=\alpha_{0}+\frac{1}{2\pi}\int_{X_{-1}^{(0)}}^{X_{-1}^{(1)}}\frac{\mathcal{D}(X_{-1})}{1-2X_{1}X_{-1}}{\rm
d}X_{-1}\,,$ (7)

or

$\displaystyle\mathcal{D}(X_{1}^{(1)})=\mathcal{D}_{0}\exp\left(\int_{X_{1}^{(0)}}^{X_{1}^{(1)}}\frac{X_{1}}{X_{2}}{\rm
d}X_{1}\right)\,,$ (8)
$\displaystyle\alpha(X_{1}^{(1)})=\alpha_{0}+\frac{1}{2\pi}\int_{X_{1}^{(0)}}^{X_{1}^{(1)}}\frac{\mathcal{D}(X_{1})}{X_{2}}{\rm
d}X_{1}.$ (9)

We illustrate this general procedure with a few examples. For the simple case
of a Maxwell-Boltzmann gas, the EoS in terms of the $X_{\nu}$’s can be
obtained analytically and one gets for example $X_{1}X_{-1}=1$ and
$X_{2}=X_{1}^{2}$. For an interacting 2D gas, the EoS is not known
analytically; however for the bosonic case, it can be approximated in two
limiting cases. For $\mu<0$, the gas is only weakly degenerate and the mean-
field energy of an atom in the gas can be written $2\tilde{g}(\hbar^{2}n/m)$,
where the dimension-less coefficient $\tilde{g}$ (assumed here $\ll 1$)
characterizes the strength of the interaction. The thermodynamics is then well
described by the prediction of the Hartree–Fock theory Hadzibabic and Dalibard
(2011)

$\mathcal{D}=-\ln\left(1-e^{\alpha-\tilde{g}\mathcal{D}/\pi}\right),$ (10)

from which we extracted numerically the values of $X_{-1}$ and $X_{1}$ and
plotted the corresponding EoS in Fig. 1. In the opposite case of a strongly
degenerate gas (chemical potential positive and larger than $k_{\rm B}T$), the
gas is described by the Thomas–Fermi equation
$\mathcal{D}=2\pi\,\alpha/\tilde{g}$. All $X_{\nu}$ are then constant, with
$X_{-1}=\tilde{g}/2$, $X_{1}=1/\tilde{g}$ and $X_{2}=0$.

Figure 1: (Color online) Determination of the EoS with variables $X_{-1}$ and
$X_{1}$, along with known limits for the EoS. The simple cases of the ideal
Bose gas (Boltzmann gas) is shown in blue dashed (dotted) line. The known
limits of the EoS of the weakly interacting 2D Bose gas are indicated by a
black point for the Thomas–Fermi limit and by a black full line for the
Hartree–Fock mean field theory. The red line results from the averaging over
all the separate intensity profiles, with the error bars corresponding to the
standard error introduced by the averaging procedure. The data shown here
contains $\sim$ 100 different values of $X_{-1}$. Inset: distribution of
measured values of $X_{-1}$ and $X_{1}$. The gray level indicates the number
of individual data points falling in each pixel.

We now turn to the practical implementation of this method for processing data
obtained with a quasi-2D rubidium gas. Our experimental preparation follows
the lines detailed in Rath et al. (2010); Yefsah et al. (2011). We start with
a 3D gas of 87Rb atoms, confined in their $F=m_{F}=2$ state in a magnetic
trap. To create a 2D system, we shine an off-resonant blue-detuned laser beam
on the atoms, with an intensity node in the plane $z=0$. The resulting
potential provides a strong confinement perpendicular to this plane, with
oscillation frequency $\omega_{z}/2\pi=1.9\,(2)$ kHz, which decreases at most
by 5% over typical distribution radii. This corresponds to the interaction
strength $\tilde{g}=\sqrt{8\pi}\,a/l_{z}\approx 0.1$, where $a$ is the 3D
scattering length and $l_{z}=\sqrt{\hbar/m\omega_{z}}$ Petrov et al. (2000).
The energy $\hbar\omega_{z}$ is comparable to the thermal energy $k_{\rm B}T$,
which ensures that most of the atoms occupy the ground state of the potential
along $z$ (see Yefsah et al. (2011) and below). The magnetic trap provides a
harmonic confinement in the $xy$ plane, with mean oscillation frequency
$\omega_{r}/2\pi=20.6\,(1)$ Hz.

After letting the cloud equilibrate for 3 s in this combined magnetic and
dipole trap, we measure the density distribution $n(\bm{r})$ by performing
_in-situ_ absorption imaging with a probe beam perpendicular to the plane of
the atoms. In order to overcome the multiple scattering effects stemming from
the high spatial density of our sample Rath et al. (2010); Chomaz et al.
(2012), the intensity of the probe beam $I$ is chosen to be much larger than
the saturation intensity $I_{\rm sat}$ of the $D_{2}$ resonance line of Rb
(typically, $I/I_{\rm sat}\approx 40$). In this regime, the response of
individual atoms is independent of their environment Reinaudi et al. (2007).
Owing to the photonic shot-noise, this method inherently suffers from a low
signal-to-noise ratio. For this reason, we complement it by the conventional
low-intensity absorption imaging ($I/I_{\rm sat}\approx 0.7$). In this regime,
the noise level is much lower, though high atomic densities cannot be
measured. Since both imaging processes destroy the atomic distribution, each
sample is prepared twice and imaged successively in the low and high intensity
regime. For the data analyzed below, we used $\approx 80$ images, with
temperatures ranging from 30 nK to 150 nK, and atoms numbers from 25 000 to
120 000.

In Fig. 2 we show typical density distributions of 2D atomic clouds, together
with the corresponding function $n[V(\bm{r})]$. The cloud (a) exhibits a
significant thermal fraction, contrarily to cloud (b), which is essentially in
the Thomas-Fermi regime. The latter illustrates the power of this fit-free
method, since it can be incorporated as such in our determination of the EoS.
On the opposite, it would be discarded in a conventional approach, owing to
the impossibility to assign it a temperature.

Figure 2: (Color online) (a) and (b) Density distributions of 2D atomic
samples of 87Rb corresponding to a partially degenerate (a) and a strongly
degenerate cloud (b). (c) Corresponding function $n[V(\bm{r})]$resulting from
azimuthal averaging. The distributions are obtained with high intensity
imaging.

Though both choices of variables $(X_{-1},X_{1})$ and $(X_{1},X_{2})$ are in
principle possible, the latter requires the experimental evaluation of a
second-order derivative, which often suffers from a poor signal-to-noise
ratio. By contrast the choice $(X_{-1},X_{1})$, also adopted in Ku et al.
(2012) when writing the EoS in terms of pressure and compressibility, appears
particularly robust 333In our case, the compressibility is given by
$\kappa=X_{1}/(n^{2}\,h^{2}/m)$ and the pressure by
$P=X_{-1}\,n^{2}\,h^{2}/m$. For each image, we perform an azimuthal average
and compute a set of $\approx 70$ data points $(X_{-1},\,X_{1}$), where the
low (high) values of $X_{-1}$ correspond to the high (low) density regions of
the image.

In a first step, we combine all sets obtained from images acquired at various
temperatures and various atom numbers to test the scale invariance. As
explained above, each individual measurement should sit on the same universal
curve in the $(X_{-1},X_{1})$ plane, provided the interaction strength
$\tilde{g}$ is constant. We show in the inset of Fig. 1 the repartition of
data points in the $(X_{-1},X_{1})$ plane, which fall as expected around a
single curve.

After having checked the (approximate) scale invariance of our quasi-2D Bose
gases, we average all data sets $(X_{-1},X_{1})$ obtained from the ensemble of
images and arrive at the determination of the EoS shown in Fig. 1 444We
extract data points with $X_{-1}<2.5$ ($X_{-1}>2.5$) from high (low) intensity
imaging.. In order to re-express this EoS in terms of the more traditional
variables $\alpha$ and $\mathcal{D}$, we now need to apply the transformations
of Eqs. (6) and (7). However, this transformation must be adapted to account
for possible imperfections in the calibration of the detectivity of our
imaging setup. Indeed, as in most cold atoms experiments, we only measure the
density up to a global multiplicative factor $\beta$ 555This factor is taken
into account in Eq. 1 by $E=\beta(\hbar^{2}/m)n$. The definition of $X_{\nu}$
is unchanged., which is defined as the ratio between the unknown actual
absorption cross-section and the ideal one expected for monochromatic probe
light in the absence of stray magnetic fields. Taking this calibration factor
into account amounts to replacing Eqs. (4) and (5) by

$\displaystyle\mathcal{D}(X_{\nu}^{(1)})=\mathcal{D}_{0}\exp\left(\int_{X_{\nu}^{(0)}/\beta^{\nu}}^{X_{\nu}^{(1)}/\beta^{\nu}}\frac{X_{1}}{(\nu-1)X_{1}X_{\nu}+X_{\nu+1}}{\rm
d}X_{\nu}\right)$ (11)
$\displaystyle\alpha(X_{\nu}^{(1)})=\alpha_{0}+\frac{\beta}{2\pi}\int_{X_{\nu}^{(0)}/\beta^{\nu}}^{X_{\nu}^{(1)}/\beta^{\nu}}\frac{\mathcal{D}(X_{\nu})}{(\nu-1)X_{1}X_{\nu}+X_{\nu+1}}{\rm
d}X_{\nu}$ (12)

where the bounds of the integrals now depend on $\beta$ and where
$X_{\nu}^{(0)}/\beta^{\nu}$ corresponds to the reference values $\alpha_{0}$
and $\mathcal{D}_{0}$. The value of $\beta$ is _a priori_ unknown; however, it
can be determined by performing a least-square fit of the prediction of the
Hartree–Fock mean-field theory to the measured EoS over a proper region of
phase-space densities. In other words, the fact that the value of $\beta$ is
unknown can be handled in this method, provided one has a good theoretical
knowledge of the EoS over a segment in parameter space, rather than at a
single point $(\alpha_{0},\mathcal{D}_{0})$. Here, we consider the EoS to be
well approximated by the Hartree–Fock mean-field theory for $\alpha<0$.

Figure 3: (Color online) (a) Equation of state of the 2D Bose gas, determined
with Eqs. 6 and 7 (red points), with a detailed view of the critical region
around the BKT transition (blue dash-dotted line) in (b). Statistical error
bars are too small to be shown on these plots. We show for comparison the
classical field Monte Carlo prediction $\mathcal{D}_{\rm th}$ Prokof’ev and
Svistunov (2002) in black squares, the Thomas–Fermi limit in black dash-dotted
line and the Hartree–Fock mean-field theory in black full line. We provide a
quantitative estimate of the difference between measurement and prediction in
(c). There, we plot $\mathcal{D}/\mathcal{D}_{\rm th}-1$, where zero indicates
perfect agreement. The error bars result from a boostrap analysis of the
experimental data.

We choose the central point $(\alpha_{0},\mathcal{D}_{0})$ of the fitting
segment at $X_{-1}^{0}=3$, which corresponds to a phase-space density
$\mathcal{D}_{0}=1.45$ and $\alpha_{0}=-0.22$, well within the Hartree–Fock
mean-field regime. This choice is a good compromise between higher values of
$X_{-1}$, which have a significantly higher noise level and lower values of
$X_{-1}$, for which the EoS is not known. The fit of the Hartree–Fock mean-
field theory yields a detectivity $\beta=0.456\,(1)$ 666The error bar on this
value is calculated through a bootstrap analysis. The result is consistent
with our (less precise) previous estimate for the same setup, $\beta=0.40(2)$,
which was based on a different method Yefsah et al. (2011).. The EoS in terms
of the variables $(\alpha,\mathcal{D})$ – obtained after a small correction
due to excited states of the $z$ motion (see below) – is shown in Fig. 3a,
along with the numerical prediction $\mathcal{D}_{\rm th}$ Prokof’ev and
Svistunov (2002). Note that the reconstructed EoS is remarkably smooth and
doesn’t display any particular feature at the transition point. This
observation is also made on the EoS for pressure, entropy and heat capacity
shown in the Supplementary Material. This fact illustrates the “infinite
order” nature of the BKT transition which is not associated to any singularity
of thermodynamic quantities777A singularity however appears in the superfluid
density, which has to be distinguished from the total density, as the former
enters into account only in the response of the system to an externally driven
drag or rotation., as opposed to phase transitions driven by the breaking of a
continuous symmetry, such as the second order lambda transition observed at
MIT Ku et al. (2012). To compare quantitatively the reconstructed EoS with the
numerical prediction, we plot the quantity $\mathcal{D}/\mathcal{D}_{\rm
th}-1$ in Fig. 3, and find it to lie consistently below 15 %, and even below 5
% around the phase transition, which occurs at $\mu_{C}/k_{\rm B}T\approx
0.17$ Prokof’ev and Svistunov (2002). The deviation observed in the
fluctuation region below the critical chemical potential might signal
deviations to the classical field picture which is expected to be accurate for
$\tilde{g}\ll\,1$ Prokof’ev et al. (2001); Prokof’ev and Svistunov (2002).
Theoretically this deviation could be addressed using Quantum Monte Carlo
methods Holzmann and Krauth (2008); Rançon and Dupuis (2012).

We have emphasized that the above method does not require any fit of
individual images, nor the knowledge of their temperature $T$ or chemical
potential $\mu$. Furthermore, once the EoS has been determined, it is possible
to fit the full EoS to the individual images and to infer their $(T,\mu)$. We
find from this _a posteriori_ fit that the temperatures of our gases range
from 30 nK to 150 nK, which is of the same order as the harmonic oscillator
level spacing $\hbar\omega_{z}/k_{\rm B}\sim 100$ nK. Consequently, our
warmest clouds do not strictly lie in the 2D regime; some atoms occupy excited
states of the $z$ motion, which breaks the scale invariance of the EoS.
Fortunately, as a consequence of Bose statistics, this fraction of atoms
remains small even when $k_{\rm B}T\sim\hbar\omega_{z}$. Their contribution to
the total density can then be self-consistently evaluated and removed, in a
similar manner to Tung et al. (2010); Yefsah et al. (2011). The remaining
ground-state density can then be used to refine the determination of the EoS.
In the data presented in Fig. 3, this subtraction of the contribution of the
excited states of the $z$ motion has already been performed. The exact
procedure used for the excited state subtraction is presented in detail in the
Supplemental Material, together with the initial EoS obtained before this
subtraction (the deviations with respect to the final one shown in Fig. 3 are
minute).

In conclusion, we have presented a method to determine the EoS of a scale-
invariant fluid. It is based on the identification of a set of dimensionless
variables $X_{\nu}$ directly related to the measured density of the fluid.
This method does not rely on thermometry of individual images, nor on the
precise calibration of the detectivity, and leads to a strong reduction of the
noise level in the measurement. We have applied it to the case of a weakly
interacting Bose gas and obtained its EoS with a precision of few percents, in
excellent agreement with a theoretical prediction obtained from a classical
Monte Carlo simulation. Using the response of the gas to a gauge field,
originating for example from a rotation, this method could be extended to
access the superfluid faction of the gas along the lines proposed in Ho and
Zhou (2009). In principle, this method is not limited to scale invariant
systems, and could be extended to any situation described by two independent
dimensionless parameters, such as the zero temperature limit of the Fermi gas,
either for a spin-balanced gas with varying interactions Navon et al. (2010),
or for a unitary spin-imbalanced Fermi gas Shin (2008).

###### Acknowledgements.

We thank Jérôme Beugnon, Jason Ho, Boris Svistunov and Martin Zwierlein for
useful discussions and comments. This work is supported by IFRAF, ANR
(ANR-12-BLAN-AGAFON), ERC (Synergy UQUAM). L. Chomaz and L. Corman acknowledge
the support from DGA, and C. Weitenberg acknowledges the support from the EU
(PIEF-GA-2011-299731).

## References

  * Ho and Zhou (2009) T.-L. Ho and Q. Zhou, Nature Physics 6, 131 (2009).
  * Nascimbène et al. (2010) S. Nascimbène, N. Navon, K. J. Jiang, F. Chevy, and C. Salomon, Nature 463, 1057 (2010).
  * Navon et al. (2010) N. Navon, S. Nascimbène, F. Chevy, and C. Salomon, Science 328, 729 (2010).
  * Ku et al. (2012) M. J. H. Ku, A. T. Sommer, L. W. Cheuk, and M. W. Zwierlein, Science 335, 563 (2012).
  * Van Houcke et al. (2012) K. Van Houcke, F. Werner, E. Kozik, N. Prokof’ev, B. Svistunov, M. J. H. Ku, A. T. Sommer, L. W. Cheuk, A. Schirotzek, and M. W. Zwierlein, Nature Physics 8, 366 (2012).
  * Prokof’ev et al. (2001) N. Prokof’ev, O. Ruebenacker, and B. Svistunov, Physical Review Letters 87, 270402 (2001).
  * Prokof’ev and Svistunov (2002) N. V. Prokof’ev and B. V. Svistunov, Phys. Rev. A 66, 043608 (2002).
  * Giorgini et al. (2008) S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).
  * Hung et al. (2011) C.-L. Hung, X. Zhang, N. Gemelke, and C. Chin, Nature 470, 236 (2011).
  * Yefsah et al. (2011) T. Yefsah, R. Desbuquois, L. Chomaz, K. J. Günter, and J. Dalibard, Phys. Rev. Lett. 107, 130401 (2011).
  * Hadzibabic and Dalibard (2011) Z. Hadzibabic and J. Dalibard, Rivista del Nuovo Cimento 34, 389 (2011).
  * Rath et al. (2010) S. P. Rath, T. Yefsah, K. J. Günter, M. Cheneau, R. Desbuquois, M. Holzmann, W. Krauth, and J. Dalibard, Phys. Rev. A 82, 013609 (2010).
  * Petrov et al. (2000) D. S. Petrov, M. Holzmann, and G. V. Shlyapnikov, Phys. Rev. Lett. 84, 2551 (2000).
  * Chomaz et al. (2012) L. Chomaz, L. Corman, T. Yefsah, R. Desbuquois, and J. Dalibard, New Journal of Physics 14, 055001 (2012).
  * Reinaudi et al. (2007) G. Reinaudi, T. Lahaye, Z. Wang, and D. Guéry-Odelin, Optics Letters 32, 3143 (2007).
  * Holzmann and Krauth (2008) M. Holzmann and W. Krauth, Physical Review Letters 100, 190402 (2008).
  * Rançon and Dupuis (2012) A. Rançon and N. Dupuis, Physical Review A 85, 063607 (2012).
  * Tung et al. (2010) S. Tung, G. Lamporesi, D. Lobser, L. Xia, and E. A. Cornell, Phys. Rev. Lett. 105, 230408 (2010).
  * Shin (2008) Y.-I. Shin, Physical Review A 77, 041603 (2008).
  * Vogt et al. (2012) E. Vogt, M. Feld, B. Fröhlich, D. Pertot, M. Koschorreck, and M. Köhl, Phys. Rev. Lett. 108, 070404 (2012).
  * Olshanii et al. (2010) M. Olshanii, H. Perrin, and V. Lorent, Phys. Rev. Lett. 105, 095302 (2010).
  * Langmack et al. (2012) C. Langmack, M. Barth, W. Zwerger, and E. Braaten, Phys. Rev. Lett. 108, 060402 (2012).

Supplemental material

## Removing the contribution for thermally populated transverse states

The method presented in the main text assumes the EoS to be strictly scale
invariant. However, when the temperature $T$ of an individual realization is
on the order of the harmonic oscillator level spacing $\hbar\omega_{z}/k_{\rm
B}\sim 100$ nK, a fraction of the atoms occupy transverse excited states of
the confinement potential. Their contribution can be evaluated and self-
consistently removed in the following manner.

  1. 1.

We calculate the EoS in variables $(X_{-1},X_{1})$ with the method described
in the main text. This process yields a preliminary determination of the EoS,
as well as a value for the detectivity $\beta_{(0)}=0.481(4)$. Both these
measurements contain a systematic error introduced by the thermally populated
transverse levels.

  2. 2.

We measure the temperature and the chemical potential of the atomic
distribution of each individual image by fixing the detectivity $\beta$ to the
value determined above, and by fitting with the EoS measured previously. Since
the EoS contains excited levels in its determination, we expect that the
measured temperature will also be affected. Note that the population of the
excited states is only 10% of the total population at most: we therefore
expect a similarly small shift of the temperature and chemical potential.

  3. 3.

We self-consistently evaluate the contribution of the excited levels to the
total density, assuming the atoms in the excited states of the $z$ motion to
be in the HFMF regime Tung et al. (2010); Yefsah et al. (2011) and subtract
them to obtain an estimate of the population of the ground state. We find that
the phase-space density associated with the excited states is at most 0.5,
which justifies to describe these states with the Hartree–Fock mean-field
approximation.

  4. 4.

With the new estimate of the population of the ground state for each image, we
generate a new EoS with the method outlined in the main text. We also
determine

Figure S1: (Color online) EoS of the 2D Bose gas, determined with Eqs. 6 and
7. The measurement without (with) subtraction of the population of the excited
levels is shown in blue (red) line.

a more accurate value of the detectivity, and find $\beta=$$0.456(1)$.

  5. 5.

In principle, this process can be iterated to obtain an even more accurate
determination of the temperature for each image. However, this does not lead
to an improved measurement.

As shown in Fig. S1, the subtraction of the excited states does not
significantly affect the EoS, though the detectivity is significantly
modified.

Figure S2: (Color online) Equations of state of the 2D Bose gas, derived from
the Hartree–Fock mean-field theory (black line) and reconstructed following
our method (red points), expressed with (a) the pressure $P$, (b) the entropy
per particle $S/N$ and (c) the heat capacity per particle $C/N$.

## EoS in other variables

Owing to the scale invariance of the two dimensional Bose gas, the EoS with
any pair of variables can in principle be deduced from the one in the main
text. We show in Fig. S2 three measurements of such EoS. For example, the
pressure is deduced from the density by $n=(\partial P/\partial\mu)_{T}$ and
is expressed in dimensionless units:

$\mathcal{P}=\frac{P\,\lambda_{\rm T}^{2}}{k_{\rm
B}}=\frac{1}{2\pi}\mathcal{D}^{2}\,X_{-1}.$ (13)

Once the pressure is known, extensive variables such as the entropy and the
heat capacity can be deduced. The former is derived from the entropy per unit
area $s=(\partial P/\partial T)_{\mu}$, and the latter from the internal
energy $C=(\partial U/\partial T)_{\mu}$. Expressed in terms of previously
determined thermodynamic quantities, we have

$\frac{S}{N\,k_{\rm
B}}=2\frac{\mathcal{P}}{\mathcal{D}}-\alpha\quad\mathrm{and}\quad\frac{C}{N\,k_{\rm
B}}=\frac{\mathcal{P}}{\mathcal{D}}-\alpha+\alpha\,X_{-1}\,X_{1}$ (14)