# Proper ab-initio dielectric function of 2D materials and their polarizable
thickness

Lorenzo Sponza  Laboratoire d’Etude des Microstructures, ONERA-CNRS, UMR104,
Université Paris-Saclay, BP 72, 92322 Châtillon Cedex, France  François
Ducastelle  Laboratoire d’Etude des Microstructures, ONERA-CNRS, UMR104,
Université Paris-Saclay, BP 72, 92322 Châtillon Cedex, France

###### Abstract

In this paper we derive a formalism allowing us to separate inter-layer
contributions to the polarizability of a periodic array of 2D materials from
intra-layer ones. To this aim, effective profile functions are introduced.
They constitute a tight-binding-like layer-localized basis involving two
lengths, the effective thickness $d$ characteristic of the 2D material and the
inter-layer separation $L$. The method permits, within the same formalism,
either to compute the single-layer dielectric function from an ab-initio
periodic calculation (top-down strategy) or to stack several 2D materials to
generate a finite-thickness van der Waals heterostructure (bottom-up
strategy).

Most ab initio codes assume periodic boundary conditions in the 3D space. This
framework is inappropriate to simulate isolated systems, interfaces, defects,
amorphous materials and any non-periodic or partially periodic system.
Isolated 2D sheets fall in this category, with a crystalline structure
periodic in one plane (the $xy$-plane, in the following) and isolated along
the vertical direction $z$. Similarly, the simulation of multilayers or van
der Waals heterostructures is often prohibitive because of the number of atoms
involved, even when lattice mismatch or misalignment between the constituent
layers are neglected. In this paper we address both problems within the same
formalism in a fully ab initio approach.

The first problem is about the calculation of single-layer properties from
3D-periodic ab initio calculations (top-down strategy). In this framework,
fictitious Coulomb interactions arise between the replicas of the system. The
brute force approach consists in creating large simulation cells (supercell:
SC) enclosing the planar unitary cell plus an amount of vacuum large enough to
separate the periodic replicas of the system. A vacuum of $\sim$10 Å is often
required to converge regular ground-state simulations. Though, when it comes
to computing the polarizability $\chi$, for instance in the random phase
approximation (RPA), local dipoles are created. These charge inhomogeneities
worsen the convergence problem because the artifact interactions are
unscreened and hence long-ranged. Furthermore $\chi$ is normalized with
respect to the volume of the simulation cell, and hence depends on the amount
of vacuum which is arbitrary Cudazzo  _et al._ (2011); Sponza  _et al._
(2016); Tian  _et al._ (2020). In our opinion, this renormalization problem
has been too much overlooked in the past; mainly because it is absent in
ground-state calculations. As a result, in the SC scheme, excited state
calculations seldom converge with the amount of vacuum.

However, state-of-the-art calculations are nowadays done within the Coulomb
truncation (CT) scheme Ismail-Beigi (2006); Rozzi  _et al._ (2006); Hüser  _et
al._ (2013); Cudazzo  _et al._ (2016). It consists in letting the Coulomb
interaction vanish at a distance $L/2$ from the sheet, where $L$ is the height
of the simulation cell. This scheme kills by construction the spurious
interactions and permits to reduce $L$ to values comparable to those needed in
SC ground-state simulations. Though, we want to stress that $L$ is still an
arbitrary quantity even in the CT scheme.

The second problem consists in stacking different 2D layers to get finite-
thickness homo- or hetero-structures (bottom-up strategy). A solution to the
two problems must rely on the correct description of the inter-layer
interaction and also on the explicit account for the (effective) thickness of
the layers.

In this Letter we present a general formalism allowing us to separate
analytically the intra-layer from inter-layer contributions to the
polarizability $X$ of a periodic array of layers. This is done by defining
profile functions while computing the polarizability of a single layer. In
this way, the volume-normalization and the fictitious-interaction problems are
treated on the same footing, completing and generalizing recent developments
Nazarov  _et al._ (2014); Nazarov (2015); Latini  _et al._ (2015); Thygesen
(2017); Meckbach  _et al._ (2018); Tian  _et al._ (2020); Gjerding  _et al._
(2020) concerning the top-down strategy. Then, we consider the general case of
a finite-thickness heterostructure and we use the same general framework to
implement a bottom-up strategy similar to the Quantum Electrostatic
Heterostructure (QEH) model recently put forward Andersen  _et al._ (2015);
Thygesen (2017); Gjerding  _et al._ (2020); Cavalcante  _et al._ (2018). With
the intent of removing any arbitrariness from the calculation, we also provide
a recipe to compute the profile functions ab initio and use them to calculate
the macroscopic dielectric function of a single layer. This development allows
us to demonstrate that for the dielectric function to be meaningful, it is
mandatory to take into account the finite thickness of the film.

#### I: Single-layer polarizability

We start from the case of an isolated single layer without any periodicity
along the $z$ axis. We introduce a mixed space representation
$\left|{\mathbf{k}_{\parallel},z}\right)$ for which a generic quantity depends
on the in-plane momentum
($\mathbf{k}_{\parallel}=\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}$, sum
of the crystal momentum $\mathbf{q}_{\parallel}$ and of a reciprocal lattice
vector $\mathbf{G}_{\parallel}$ 111Note that in this notation, the continuous
part never changes. So if
$\mathbf{k}_{\parallel}=\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}$, then
$\mathbf{k}_{\parallel}^{\prime}=\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}^{\prime}$
and analogously $k_{z}^{\prime}=q_{z}+G_{z}^{\prime}$.), and on the real-space
coordinate $z$. We proceed assuming that, we can separate in-plane ($\varphi$)
from out-of-plane ($\xi$) components of the Bloch wavefunction:
$\phi_{i\bm{\kappa}_{\parallel}}(\mathbf{r})=\varphi_{i\bm{\kappa}_{\parallel}}(\bm{\rho})\,\xi(z)$,
where $\bm{\rho}$ is the in-plane component of $\mathbf{r}$. The irreducible
polarizability $\chi^{0}$ is calculated as usual as a sum over all
independent-particle transitions, the oscillator strengths of which involve
matrix elements of the form
$\int\phi^{*}_{i\bm{\kappa}_{\parallel}}(\mathbf{r})e^{i(\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel})\cdot\bm{\rho}}\phi_{j\bm{\kappa}_{\parallel}+\mathbf{q}_{\parallel}}(\mathbf{r})\,d\bm{\rho}$.
This makes appear the normalized effective profile functions
$\theta^{0}=\xi^{2}$ which are positive, real and even functions of $z$ and
localized around $z=0$ (i.e. they vanish for $|z|>d/2$ for a given distance
$d$). This leads to the expression for the irreducible polarizability of the
single layer:

$\chi^{0}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega,z,z^{\prime})=\theta^{0}(z)\chi^{0}_{2D\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)\theta^{0}(z^{\prime})\,.$
(1)

The functions $\theta^{0}$ play here the role of effective polarizable
profiles and they can be associated to the thickness of the isolated layer
through the characteristic length $d$. This thickness is defined more
precisely below. It turns out that the numerical extraction of the profile
function shows a mild dependence on the in-plane momentum, and in the
following we include such a dependence. Then,
$\theta^{0}_{\mathbf{k}_{\parallel}}(z)$ is characterized by an effective
thickness $d_{\mathbf{k}_{\parallel}}$. Since our factorization hypothesis is
done at the wave-function level, we assume no $\omega$-dependence of
$\theta^{0}$.

The random phase approximation (RPA) to the polarizability of the slab $\chi$
satisfies a Dyson equation of the form $\chi=\chi^{0}+\chi^{0}v\chi$. As a
consequence of the $\theta^{0}$ functions appearing in (1), one can show that
$\chi=\theta^{0}\chi_{2D}\theta^{0}$ where the extension in $(z,z^{\prime})$
is carried only by the $\theta^{0}$ functions, while the polarization is
embodied in
$\chi_{2D\,\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)$
which solves the 2D Dyson equation:

$\chi_{2D}=\chi^{0}_{2D}+\chi^{0}_{2D}\mathcal{V}^{\text{in}}\chi_{2D}\,.$ (2)

Here $\mathcal{V}^{\text{in}}$ is the layer-projected Coulomb interaction;

$\mathcal{V}^{\text{in}}(\mathbf{k}_{\parallel})=\frac{2\pi}{|\mathbf{k}_{\parallel}|}\int\theta^{0}_{\mathbf{k}_{\parallel}}(z)e^{-|\mathbf{k}_{\parallel}|\,|z-z^{\prime}|}\theta^{0}_{\mathbf{k}_{\parallel}^{\prime}}(z^{\prime})\,dz\,dz^{\prime}\,.$
(3)

For further details on the entire derivation, see appendix A of the
Supplemental Material.

#### II: Periodic array of layers

In order to study the polarizability $X$ of an array of layers with period $L$
along $z$, we introduce a layer-projected representation
$\left|{\mathbf{k}_{\parallel},n}\right)$ where $n$ is the index of the layer.
We define
$X_{\mathbf{G}_{\parallel},\mathbf{G}_{\parallel}^{\prime},nm}(\mathbf{q}_{\parallel},\omega)=\left({\mathbf{k}_{\parallel},n}\right|X(\omega)\left|{\mathbf{k}_{\parallel}^{\prime},m}\right)$
and the projection from
$\left({\mathbf{k}_{\parallel},z}\left|{\mathbf{k}_{\parallel}^{\prime},n}\right.\right)=\sqrt{d_{\mathbf{k}_{\parallel}}}\theta^{n}_{\mathbf{k}_{\parallel}}(z)\delta_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}$,
where
$\theta^{n}_{\mathbf{k}_{\parallel}}(z):=\theta^{0}_{\mathbf{k}_{\parallel}}(z-nL)$
is the profile function localized on the $n$th layer. Details on the
representations can be found in Appendix B of the Supplemental Material. Since
$\chi^{0}$ does not imply any long-range term, the overlap of successive
layers can be neglected if we further assume
$L>d_{\mathbf{k}_{\parallel}}\,\forall\,\mathbf{k}_{\parallel}$. So $X^{0}$
becomes the sum over the periodically repeated $\chi^{0}$s (equation (1)):

$\displaystyle
X^{0}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega,z,z^{\prime})=$
$\displaystyle\qquad=$
$\displaystyle\sum_{n}\theta^{n}_{\mathbf{k}_{\parallel}}(z)\chi^{0}_{2D\,\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)\theta^{n}_{\mathbf{k}_{\parallel}^{\prime}}(z^{\prime})$
$\displaystyle\qquad=$
$\displaystyle\sqrt{d_{\mathbf{k}_{\parallel}}d_{\mathbf{k}_{\parallel}^{\prime}}}\sum_{n,m}\theta^{n}_{\mathbf{k}_{\parallel}}(z)X^{0}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime},nm}(\mathbf{q}_{\parallel},\omega)\theta^{m}_{\mathbf{k}_{\parallel}^{\prime}}(z^{\prime})\;,$

depending on the representation employed (mixed-space above, or layer-
projected below). By equating these two expressions one gets
$\sqrt{d_{\mathbf{k}_{\parallel}}d_{\mathbf{k}_{\parallel}^{\prime}}}X^{0}_{\mathbf{k}_{\parallel}\mathbf{k}_{\parallel}^{\prime}\,nm}=\delta_{nm}\chi^{0}_{2D\,\mathbf{k}_{\parallel}\mathbf{k}_{\parallel}^{\prime}}$.
Using this equality and by defining
$\hat{\chi}^{(0)}_{\mathbf{k}_{\parallel}\mathbf{k}_{\parallel}^{\prime}}:=\sqrt{d_{\mathbf{k}_{\parallel}}d_{\mathbf{k}_{\parallel}^{\prime}}}X^{(0)}_{\mathbf{k}_{\parallel}\mathbf{k}_{\parallel}^{\prime}}$
and
$\hat{v}(\mathbf{k}_{\parallel})=v(\mathbf{k}_{\parallel})/d_{\mathbf{k}_{\parallel}}$,
one arrives to the layer-projected Dyson equation:

$\hat{\chi}_{nm}=\hat{\chi}^{0}_{nm}+\sum_{p,s}\hat{\chi}^{0}_{np}\hat{v}_{ps}\hat{\chi}_{sm}\,.$
(5)

We notice that
$\hat{v}_{nn}(\mathbf{k}_{\parallel})=\left({\mathbf{k}_{\parallel},n}\right|\hat{v}\left|{\mathbf{k}_{\parallel},n}\right)\equiv\mathcal{V}^{\text{in}}(\mathbf{k}_{\parallel})$
as in (3). Using the definitions introduced above, the $n$th layer term reads
$\hat{\chi}_{nn}=\chi^{0}_{2D}+\chi^{0}_{2D}\mathcal{V}^{\text{in}}\chi^{0}_{2D}+O(2)$.
We conclude that $\hat{\chi}_{nn}\equiv\chi_{2D}$ which is layer-independent
and solves the single-layer equation (2). Note that all matrices are written
in $\mathbf{G}_{\parallel}$ space, so that in-plane local fields are fully
taken into account.

Let us now focus on the off-diagonal elements of (5). We start splitting the
Coulomb interaction
$\hat{v}_{nm}=\mathcal{V}^{\text{in}}_{nm}+\mathcal{V}^{\text{off}}_{nm}$ with
the definitions
$\mathcal{V}^{\text{in}}_{nm}=\delta_{nm}\mathcal{V}^{\text{in}}$, and
$\mathcal{V}^{\text{off}}_{nm}=\hat{v}_{nm}$ if $n\neq m$ and 0 otherwise.
Equation (5) can now be split into two equations:

$\displaystyle\chi_{2D,nm}$ $\displaystyle=$
$\displaystyle\delta_{nm}\left[\chi_{2D}^{0}\left(1-\mathcal{V}^{\text{in}}\,\chi_{2D}^{0}\right)^{-1}\right]=\delta_{nm}\chi_{2D}$
(6) $\displaystyle\hat{\chi}_{nm}$ $\displaystyle=$
$\displaystyle\delta_{nm}\chi_{2D}+\sum_{q}\chi_{2D}\mathcal{V}^{\text{off}}_{nq}\,\hat{\chi}_{qm}\,.$
(7)

The latter equation is the main result: It expresses the total polarizability
of the array ($\hat{\chi}$) in terms of single-layer polarizabilities
($\chi_{2D}$) interacting through an effective inter-layer Coulomb interaction
($\mathcal{V}^{\text{off}}$). This result permits to develop the bottom-up and
the top-down strategies presented previously. The advantage of our approach is
to combine them within a general and concise formalism based in particular on
the profile functions $\theta^{0}$ and their associate polarizable thickness.
They are crucial quantities for a correct calculation of the dielectric
function. We will show how to compute them from first principles, hence
removing any arbitrariness in ab initio calculations of single layers.

With the intent of connecting our derivation to standard ab initio output,
usually expressed in the $\left|{\mathbf{k}_{\parallel},k_{z}}\right)$
representation, we rewrite (7) in reciprocal space. This is done in detail in
Appendix C of the Supplemental Material; below we report just the main results
for $q_{z}=0$. $\omega$ and $\mathbf{q}=\mathbf{q}_{\parallel}$ variables are
dropped for notational convenience:

$\displaystyle X_{\mathbf{G}\mathbf{G}^{\prime}}$ $\displaystyle=$
$\displaystyle
L^{-1}\vartheta_{\mathbf{k}_{\parallel}}(G_{z})\hat{\chi}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}\vartheta_{\mathbf{k}_{\parallel}^{\prime}}(G_{z}^{\prime})$
(8)
$\displaystyle\chi_{2D,\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}$
$\displaystyle=$
$\displaystyle\hat{\chi}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}-\chi_{2D,\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime\prime}}\mathcal{V}^{\text{off}}(\mathbf{k}_{\parallel}^{\prime\prime})\hat{\chi}_{\mathbf{G}_{\parallel}^{\prime\prime}\mathbf{G}_{\parallel}^{\prime}}$
(9)
$\displaystyle\mathcal{V}^{\text{off}}(\mathbf{k}_{\parallel}^{\prime\prime})$
$\displaystyle=$ $\displaystyle
2\mathcal{V}^{\text{in}}(\mathbf{k}_{\parallel}^{\prime\prime})\left[e^{|\mathbf{k}_{\parallel}^{\prime\prime}|L}-1\right]^{-1}\,.$
(10)

In (8), $\vartheta$ is the Fourier transform along $z$ of the profile function
$\theta^{0}$. A sum over $\mathbf{G}_{\parallel}^{\prime\prime}$ is understood
in (9).

Here we present the practical implementation of the top-down method. We keep
assuming $q_{z}=0$. The first step is to extract the profile function
$\theta^{0}$ from the first-principle periodic calculation of $X$. By sampling
the column $X_{\mathbf{G},0}(\mathbf{q}_{\parallel},\omega=0)$ along $G_{z}$
at fixed $\mathbf{G}_{\parallel}$, one can extract
$\theta^{0}_{\mathbf{k}_{\parallel}}$ from first principles:

$\vartheta_{\mathbf{k}_{\parallel}}(G_{z})=\frac{X_{\mathbf{G}\,0}(\mathbf{q}_{\parallel},0)}{X_{\mathbf{G}_{\parallel}\,0}(\mathbf{q}_{\parallel},0)}\;,$
(11)

from which we obtain $\theta_{\mathbf{k}_{\parallel}}^{0}(z)$, the Fourier
transform of $\vartheta_{\mathbf{k}_{\parallel}}(k_{z})$ which, for large $L$,
can be calculated through a discrete sum
$\theta_{\mathbf{k}_{\parallel}}^{0}(z)\simeq\frac{1}{L}\sum_{G_{z}}e^{iG_{z}z}\vartheta_{\mathbf{k}_{\parallel}}(G_{z})$.
Note that this direct relation between the profile function and the
polarizability confirms the close relation identified already by Tian and
coworkers Tian  _et al._ (2020). However, our ab initio method prevents the
use of an arbitrary thickness parameters to be fixed from fitting procedures
Meckbach  _et al._ (2018) or from physical considerations Latini  _et al._
(2015).

Once the $\vartheta$s extracted, one can invert equation (8) at $\omega=0$
thus getting $\hat{\chi}$ from first principles. Next, we calculate the inter-
layer Coulomb interaction $\mathcal{V}^{\text{off}}$ as defined in (10). This
can be done numerically (NUM) relying on (3), or can be approximated. Within
the _perfect 2D_ approximation (2DA),
$\theta^{0}_{\mathbf{k}_{\parallel}}(z)=\delta(z)$. In the _sharp slab_
approximation (SSA)
$\theta^{0}_{\mathbf{k}_{\parallel}}(z)=d^{-1}\left[H\left(z+d/2\right)-H\left(z-d/2\right)\right]$
where $H(x)$ is the Heaviside step function and $d$ is an effective momentum-
independent thickness. Depending on the the approximation chosen,

$\mathcal{V}^{\text{in}}(\mathbf{k}_{\parallel})=\left\\{\begin{array}[]{ll}2\pi/|\mathbf{k}_{\parallel}|&\text{2DA}\\\
\frac{4\pi}{|\mathbf{k}_{\parallel}|^{3}d^{2}}\left(|\mathbf{k}_{\parallel}|d-1+e^{-|\mathbf{k}_{\parallel}|d}\right)&\text{SSA}\\\
\text{computed as in~{}\eqref{eq:vin}}&\text{NUM}\end{array}\right.$

Other approximations are possible (for instance the 2D Ohno potential as in
Meckbach  _et al._ (2018)). Putting it into (10) and solving (9), one finally
gets $\chi_{2D}$. We stress that, if $X$ has been computed in the CT scheme,
then $\mathcal{V}^{\text{off}}\equiv 0$ and $\hat{\chi}\equiv\chi_{2D}$.

In the bottom-up method we calculate the global polarizability of an
heterostructure composed of $N$ layers of different 2D materials. Equation (7)
can be easily generalized to this case and becomes:

$\hat{\chi}_{nm}=\delta_{nm}\chi_{2D,n}+\sum_{p=1}^{N}\chi_{2D,n}\mathcal{V}^{\text{off}}_{np}\,\hat{\chi}_{pm}\,,$
(12)

with

$\mathcal{V}^{\text{off}}_{nm}(\mathbf{k}_{\parallel})=\int\int\theta^{0}_{\mathbf{k}_{\parallel}}(z-z_{n})v(\mathbf{k}_{\parallel},z,z^{\prime})\theta^{0}_{\mathbf{k}_{\parallel}^{\prime}}(z^{\prime}-z_{m})dzdz^{\prime}\,,$
(13)

where the differences with (7) are that (i) the sum is limited to $N$ layers,
(ii) each layer has its own $\chi_{2D,n}$, and (iii) the generic layer $n$
centered at $z_{n}$ is not necessary in a periodic array. The bottom-up
equation (12) can be applied to particular cases. For instance, in a
multilayer made of the same 2D material $\chi_{2D,n}=\chi_{2D}\,\forall n$.
Also the bulk can be reconstructed by adding in the latter case $N=\infty$ and
$z_{n}=nL$, which actually boils down in solving the pristine equation (7).

#### III: The dielectric function

The inverse microscopic dielectric function of a single layer $\epsilon^{-1}$
relates the total potential $U^{\text{tot}}$ to an external potential
$U^{\text{ext}}$ according to the definition
$U^{\text{tot}}=\epsilon^{-1}U^{\text{ext}}$. The macroscopic average of these
fields is obtained through their projection on the single layer. To do so, we
project them on the $\left|{\mathbf{k}_{\parallel},n=0}\right)$ representation
(see Supp Mat Appendix B), getting
$U^{j}_{M}(\mathbf{q}_{\parallel},\omega)=\sqrt{d_{\mathbf{k}_{\parallel}}}\int\theta^{0}_{\mathbf{k}_{\parallel}}(z)U^{j}_{0}(\mathbf{q}_{\parallel},\omega,z)dz$,
for $j=(\text{ext})$ or $(\text{tot})$. Next, we define the macroscopic
dielectric function as the ratio
$\epsilon_{M}=U^{\text{ext}}_{M}/U^{\text{tot}}_{M}$. If we assume from the
beginning that $U^{\text{ext}}$ is macroscopic (i.e. it vanishes where
$\mathbf{G}\neq 0$), we make use of (1) and of $\epsilon^{-1}=1+v\chi$, we
obtain:

$\epsilon_{M}(\mathbf{q}_{\parallel},\omega)=1/\left[1+\mathcal{V}^{\text{in}}(\mathbf{q}_{\parallel})\chi_{2D\,00}(\mathbf{q}_{\parallel},\omega)\right].$
(14)

As before, $\mathcal{V}^{\text{in}}$ can be computed using the profile
functions or using some approximations (e.g. 2DA or SSA). We stop here to
stress an important point. The Coulomb term $\mathcal{V}^{\text{in}}$
appearing in (14) embodies the interaction between the electrons of the system
and an external charge confined in the slab. This is a proper definition of
the average internal dielectric constant, adopted also in Nazarov  _et al._
(2014); Nazarov (2015); Latini  _et al._ (2015). This differs from that used
by Qiu _et al._ Qiu  _et al._ (2016), who defined it as the ratio between the
screened and the bare interaction between two charges. When describing energy
loss spectra or optical responses, the pertinent macroscopic dielectric
function couples the genuine response of the $2D$ system
$\chi_{2D\,00}(\mathbf{q}_{\parallel},\omega)$ to an external field. Such a
coupling is always $4\pi/|\mathbf{q}_{\parallel}|^{2}$ and not
$\mathcal{V}^{\text{in}}(\mathbf{q}_{\parallel})$ as above. In other terms,
the correct Coulomb potential a priori is not the same in (2) and in (14).
This is also discussed by Nazarov in the case of electron loss spectroscopy
Nazarov (2015). Actually many standard 3D calculations provide values of
$\tilde{\epsilon}_{M}=1/(1+vX)$, where $v$ is either
$4\pi/|\mathbf{q}_{\parallel}|^{2}$ (in the SC scheme) or its truncated
version
$4\pi\left(1-e^{-|\mathbf{q}_{\parallel}|L/2}\right)/|\mathbf{q}_{\parallel}|^{2}$
(CT). So, instead of solving the 2D Dyson equation (9), it is then possible to
extract directly the macroscopic polarisability within the CT scheme using the
fact that by definition $\vartheta(0)=1$, which implies that
$LX_{00}(\mathbf{q}_{\parallel},\omega)=\hat{\chi}_{00}(\mathbf{q}_{\parallel},\omega)$
:

$\chi_{2D,00}(\mathbf{q}_{\parallel},\omega)=\frac{L|\mathbf{q}_{\parallel}|^{2}}{4\pi\left(1-e^{-|\mathbf{q}_{\parallel}|L/2}\right)}\left[\frac{1}{\tilde{\epsilon}_{M}(\mathbf{q}_{\parallel},\omega)}-1\right]\,,$
(15)

which can be finally plugged into (14). No additional approximation has been
done with respect to the full approach presented so far. Instead, some
approximations have to be done when deriving a similar post-processing
treatment for SC data.

Taking the macroscopic limit of equation (9), which becomes hence a scalar
equation relating $\chi_{2D,00}$ to $\hat{\chi}_{00}$, we adopt the 2DA for
$\mathcal{V}^{\text{off}}$, and then:

$\frac{1}{\chi_{2D,00}(\mathbf{q}_{\parallel},\omega)}\approx\frac{4\pi}{|\mathbf{q}_{\parallel}|^{2}}\left[\frac{|\mathbf{q}_{\parallel}|}{e^{|\mathbf{q}_{\parallel}|L}-1}+\frac{\tilde{\epsilon}_{M}(\mathbf{q}_{\parallel},\omega)}{L\left(1-\tilde{\epsilon}_{M}(\mathbf{q}_{\parallel},\omega)\right)}\right].$
(16)

An equation similar to (9), has been derived by Nazarov Nazarov  _et al._
(2014); Nazarov (2015) under the 2DA. He then derived the same post-processing
equation (16) and plugged it in the 2DA version of (14). Our method
generalizes Nazarov’s work, making it applicable to any ab initio scheme (SC
or CT) and extends it beyond the 2DA.

#### IV: Applications

Figure 1: Profile function at selected parallel momenta for hBN single layer.
A vertical shift is applied to ease the reading.

We present here results obtained by applying the top-down strategy to
hexagonal boron nitride (hBN). All calculations have been done with ABINIT
Gonze  _et al._ (2016) (wavefunctions), DP dp  (the $X^{0}$) and an in-home
code implementing our technique. The computational details are presented in
Appendix E.

The profile functions $\theta_{\mathbf{k}_{\parallel}}(z)$ extracted as in
equations (11) are reported in Figure 1 for some in-plane momenta (on the
right of each profile). Oscillations come from the numerical evaluation of the
Fourier transform, which is the reason for the sharp spike at very small
$\mathbf{q}_{\parallel}$. We verified however that our results (and in
particular $\mathcal{V}^{\text{in}}$) are stable even in presence of these
numerical issues. Notice at low $q$ the dip reflecting the $\pi_{z}$-like
electronic density. Just for comparison, we report also the full with at half
maximum of each profile as an estimate of the effective thickness
$d_{\mathbf{k}_{\parallel}}$. We see that this is a bit lower than the inter-
layer distance in bulk hBN and shrinks for larger $\mathbf{k}_{\parallel}$.

Results of the top-down strategy are plotted in Figure 2, which reports
$\epsilon_{M}(\mathbf{q},\omega=0)$ for $\mathbf{q}=\mathbf{q}_{\parallel}$ on
the $\Gamma-M$ direction of the Brillouin zone.

Let us first discuss the standard methods. In the lower part of the figure we
report $\tilde{\epsilon}_{M}$, i.e. the output of standard ab initio
calculations, in the CT scheme (downward-pointing triangles) and in the SC
scheme (upward-pointing ones). Calculations have been carried out in cells
with different heights: $L=15$ Å (black curves) and $L=30$ Å (red ones). In
the CT calculations, the truncation appears both in $X=X^{0}+X^{0}vX$ and in
$\tilde{\epsilon}_{M}=1+vX$ as it is the most common implementation. The
figure clearly shows that the CT reproduces the correct long-range limit of
$\epsilon_{M}$, while in the short range it coincides with the SC simulation.
As a consequence, the CT results have the expected shape as a function of
$\mathbf{q}$. All these are very well-known properties of the CT scheme.

However, the figure highlights also that, except for the $\mathbf{q}=0$ point,
the value of $\tilde{\epsilon}_{M}$ does actually depend on $L$, which is
still an arbitrary quantity even in this simulation scheme. This arbitrariness
can spoil, for instance, the results of quasiparticle GW corrections or
excitonic calculations because it would affect the screening between the
electron and the hole (on this subject see Refs Tian  _et al._ (2020); Qiu
_et al._ (2016)). Even accepting the argument of an error cancellation between
the GW self-energy and the excitonic kernel (which is not perfect, anyway) one
should conclude that only the position of the excitonic peak would be right,
whereas both the quasiparticle gap and the excitonic binding energy would be
unreliable. Instead, our method (full circles) gives results that are
independent of the vacuum, the small discrepancies in the two cells being due
to differences at the level of the $X^{0}$, and it does reproduce the expected
limits of $\epsilon_{M}$ and has the right functional shape. The comparison
with standard calculations shows that, even when a reasonable amount of vacuum
($L=15$ Å) is used, the CT gives results that systematically underestimate
$\epsilon_{M}$ at finite $\mathbf{q}$.

Now, let us discuss the impact of different approximations in (14). Full
circles correspond to results obtained with the full method, which means that
we solved equation (9) and evaluated numerically $\mathcal{V}^{\text{in}}$
both in (10) and (14). The empty squares correspond to the 2DA version of
(14), which is the approximation adopted also by Nazarov Nazarov  _et al._
(2014); Nazarov (2015). The same $\chi_{2D}$ of the full method has been used.
The inset shows clearly that neglecting completely the thickness of the layer
in (14) leads to unphysical results with a dramatic divergence at finite
$\mathbf{q}$. One may account for the thickness of the slab by associating a
thickness $d$ to the layer in the SSA. The results corresponding to this
approximation are reported in full squares, where we employed a thickness
$d=3.33$ Å, which corresponds to the inter-layer distance in bulk hBN. At
small $\mathbf{q}$ the full method and the SSA one give very similar results,
but they start differing at larger $\mathbf{q}$ where the SSA $\epsilon_{M}$
decreases faster, so it systematically underestimates the dielectric function
at large $\mathbf{q}_{\parallel}$. Besides this, the thickness used in the SSA
is somewhat arbitrary whereas our full method is completely ab initio.

Let us now discuss the performances of the post-processing formulae (15) and
(16). Results are reported as stars and crosses respectively and they refer to
calculations where the SSA has been adopted in (14) with $d$=3.33Å. We
preferred to exemplify the performances of the post-processing formulae with
SSA instead of NUM calculations because in the latter case one should pass
through the extraction of $\theta^{0}$, which some times may require an
involved input-output handling. It is clear that either approach leads to
results that are identical to those obtained solving the full Dyson equation
(9) within the same approximation, which validates the post-processing formula
in both the SL and CT framework.

Finally, it is worth stressing that equation (16) assumes the 2DA in (10).
This indicates that while this approximation is justified at the $\chi^{2D}$
level, it leads instead to dramatically wrong results when employed in (14),
as we demonstrated above. This observation is related to the fundamental
difference between the $\mathcal{V}^{\text{in}}$ appearing in (14) and in
(9)-(10) discussed in the previous section.

Figure 2: Dielectric function of hBN single layer at finite $\mathbf{q}$. Cell
height $L$=15 Å (black) and $L$=30 Å (red). Standard calculations (CT =
downward and SC = upward triangles), compared with (14) evaluated with
different approximations for $\mathcal{V}^{\text{in}}$ (circles and squares),
and using (15) (stars) and (16) (crosses). Inset: Divergence of (14) in the
2DA.

To summarize, we have derived a series of equations permitting to split the
polarizability of a layered material into intra- and inter-layer
contributions. An important step of this derivation is the definition of the
layer-localized profile functions $\theta^{n}$. This introduces two
characteristic lengths in the derivation, the effective thickness $d$ (which
is related to the localization of the basis and can be computed from first
principles) and the vertical periodicity of the crystal (i.e. the height of
the simulation cell) $L$. In the context of the top-down strategy, we have
shown how to obtain the polarizability of an isolated sheet of matter starting
from the polarizability of a periodic array of identical layers as computed
with standard ab initio codes. The case of hBN has been discussed in some
detail. The same formalism can be applied to the bottom-up strategy, which
consists in calculating the global response of multilayers or even layered
bulk materials from the responses of their constituent films.

Funding for this work came from the European Union’s Horizon 2020 research and
innovation program under grand agreement N° 881603 (Graphene Flagship core 3).

##  I Appendix A: Derivation of the single-layer equations

We work in the mixed space representation
$\left|{\mathbf{k}_{\parallel},z}\right)$ where
$\mathbf{k}_{\parallel}=\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}$ is the
planar component of the crystal momentum, expressed as the sum of a continuous
part $\mathbf{q}_{\parallel}$ defined inside the first Brillouin zone and a
reciprocal lattice vector $\mathbf{G}_{\parallel}$. Note that, throughout all
derivations, the short wavelength contribution never changes, i.e.
$\mathbf{k}_{\parallel}^{\prime}=\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}^{\prime}$
and analogously $k_{z}^{\prime}=q_{z}+G_{z}^{\prime}$. Instead, the variable
$z$ spans the perpendicular direction in the real space. Working in the mixed
space representation allows us to let physical quantities such as the atomic
wavefunctions or the polarizability extend perpendicularly to the layer
without assuming any periodicity in that direction. The passage from the real-
space representation $\left|{\bm{\rho},z}\right)$ to the mix-space one is the
in-plane Fourier transformation:

$\left({\bm{\rho},z}\left|{\mathbf{k}_{\parallel},z^{\prime}}\right.\right)=\frac{1}{\sqrt{\mathcal{A}}}\,e^{i\mathbf{k}_{\parallel}\cdot\bm{\rho}}\,\delta(z-z^{\prime})\,,$

where $\mathcal{A}$ is the surface of the $xy$ plane, and $\bm{\rho}=(x,y)$
groups the in-plane real-space coordinates. In the system of units we used,
$1/4\pi\epsilon_{0}=\hbar=e=1$, so the real-space representation of the
Coulomb interaction is
$v(\mathbf{r},\mathbf{r}^{\prime})=1/|\mathbf{r}-\mathbf{r}^{\prime}|$. In the
mixed space representation, it reads:

$v(\mathbf{k}_{\parallel},z,z^{\prime})=\frac{2\pi
e^{-|\mathbf{k}_{\parallel}|\,|z-z^{\prime}|}}{|\mathbf{k}_{\parallel}|}\,.$
(17)

We now make the assumption that the electronic wavefunctions
$\phi_{i\bm{\kappa}_{\parallel}}(\mathbf{r})$ are separable:

$\phi_{i\bm{\kappa}_{\parallel}}(\mathbf{r})=\varphi_{i\bm{\kappa}_{\parallel}}(\bm{\rho})\xi(z)\quad\,.$
(18)

Note that the crystal momentum is only in-plane because we are in the
framework of the isolated layer. The irreducible polarizability $\chi^{0}$ is
calculated as usual as a sum over all independent-particle transitions.

Owing to the factorization hypothesis (18),the numerator of $\chi^{0}$ is the
product of two matrix elements of the form
$\int\phi^{*}_{i\bm{\kappa}_{\parallel}}(\mathbf{r})e^{-i\left(\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}\right)\cdot\bm{\rho}}\phi_{j\bm{\kappa}_{\parallel}+\mathbf{q}_{\parallel}}(\mathbf{r})\,d\bm{\rho}=\xi^{2}(z)M_{ij\bm{\kappa}_{\parallel}}(\mathbf{q}_{\parallel},\mathbf{G}_{\parallel})$.
So two profile functions $\theta^{0}(z)=\xi^{2}(z)$ can be factored out from
the sum over all transitions, which leads to the expression :

$\chi^{0}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega,z,z^{\prime})=\theta^{0}(z)\chi^{0}_{2D,\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)\theta^{0}(z^{\prime})\,.$
(19)

In evaluating numerically the profile functions $\theta^{0}(z)$, we observed a
mild dependence on the in-plane momentum $\mathbf{k}_{\parallel}$, so we
introduce a momentum-dependent profile function
$\theta^{0}_{\mathbf{k}_{\parallel}}(z)$ which will be used from now on. The
profile functions $\theta^{0}_{\mathbf{k}_{\parallel}}(z)$ are defined as
positive, real and even functions of $z$. They are localized around $z=0$,
which means that they are negligible for $|z|>d_{\mathbf{k}_{\parallel}}/2$
where $d_{\mathbf{k}_{\parallel}}$ is a sufficiently large distance from the
slab. This characteristic distance will be actually defined more rigorously in
the layer-projected representation (see Appendix B). Also, we take the profile
functions normalized $\int\theta^{0}_{\mathbf{k}_{\parallel}}(z)dz=1$, which
implies that the dimension of $\theta^{0}_{\mathbf{k}_{\parallel}}(z)$ is the
inverse of a length.

The RPA polarizability of the slab $\chi$ satisfies a Dyson equation of the
form $\chi=\chi^{0}+\chi^{0}v\chi$. To better appreciate its structure, let us
write down the zeroth and the first order of it, while making use of
definition (19). To simplify the notation, we will drop all the frequency and
momentum variables, but it is understood that all polarizabilities are
matrices in the $\mathbf{G}_{\parallel}$ space and hence all the in-plane
local fields are taken into account. Expanding the Dyson equation order by
order we get:

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle\chi(z,z^{\prime})=\theta^{0}(z)\chi_{2D}^{0}\,\theta^{0}(z^{\prime})+$\mbox{}\hfil\\\
&\qquad+\theta^{0}(z)\chi_{2D}^{0}\mathcal{V}^{\text{in}}\chi_{2D}^{0}\,\theta^{0}(z^{\prime})+O(2)\,.\end{split}$
(20)

In the expression above we have introduced a layer-projected Coulomb
interaction

$\mathcal{V}^{\text{in}}(\mathbf{k}_{\parallel})=\int\int\theta^{0}_{\mathbf{k}_{\parallel}}(z)v(\mathbf{k}_{\parallel},z,z^{\prime})\,\theta^{0}_{\mathbf{k}_{\parallel}}(z^{\prime})\,dz\,dz^{\prime}\,,$
(21)

where $v(\mathbf{k}_{\parallel},z,z^{\prime})$ is defined as in (17). It is
easy to convince oneself that every order in $v$ is actually localized around
$z=0$ because each order is sandwiched by two profile functions left out of
the integrals. On the contrary, all the other profile functions are integrated
with the Coulomb interaction $v(\mathbf{k}_{\parallel},z_{1},z_{2})$ making
$\mathcal{V}^{\text{in}}$ appear at every order. Following this argument, one
can define the 2D reducible polarizability $\chi_{2D}$ such that the RPA
polarizability of the slab reads:

$\chi_{\mathbf{G}_{\parallel}\mathbf{G}^{\prime}_{\parallel}}(\mathbf{q}_{\parallel},z,z^{\prime},\omega)=\theta^{0}_{\mathbf{k}_{\parallel}}(z)\chi_{2D,\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)\theta^{0}_{\mathbf{k}_{\parallel}}(z^{\prime}).$
(22)

The planar polarizability $\chi_{2D}$ solves the RPA Dyson equation:

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle\chi_{2D\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)=\chi^{0}_{2D\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)+$\mbox{}\hfil\\\
&+\sum_{\mathbf{G}_{\parallel}^{\prime\prime}}\chi^{0}_{2D\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime\prime}}(\mathbf{q}_{\parallel},\omega)\mathcal{V}^{\text{in}}(\mathbf{k}_{\parallel}^{\prime\prime})\chi_{2D\mathbf{G}_{\parallel}^{\prime\prime}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)\,,\end{split}$
(23)

which depends only on the in-plane components of the crystal momentum.

##  II Appendix B: Different representations

In the paper we work mostly with three different representations: the mixed
space representation (MS) $\left|{\mathbf{k}_{\parallel},z}\right)$, the
layer-projected representation (LP) $\left|{\mathbf{k}_{\parallel},n}\right)$
and the standard reciprocal space representation (RS)
$\left|{\mathbf{k}_{\parallel},k_{z}}\right)$. Many results require to move
from one representation to the other. In defining the LP basis, we project it
on the MS basis introduced before. We can make the choice of taking the LP
basis functions proportional to the profile functions
$\theta^{0}_{\mathbf{k}_{\parallel}}$ introduced above. If we do this, and we
note that $\theta^{0}_{\mathbf{k}_{\parallel}}$ has the dimension of an
inverse length because of its normalization in space, then we can define

$\left({\mathbf{k}_{\parallel},z}\left|{\mathbf{k}_{\parallel}^{\prime},n}\right.\right)=\sqrt{d_{\mathbf{k}_{\parallel}}}\theta^{n}_{\mathbf{k}_{\parallel}}(z)\delta_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}\,,$
(24)

where $n\in\mathbb{Z}$ is the index of the layer and the basis functions
$\theta^{n}_{\mathbf{k}_{\parallel}}(z):=\theta^{0}_{\mathbf{k}_{\parallel}}(z-nL)$
satisfy the orthonormality relations:

$\displaystyle 1$ $\displaystyle=$
$\displaystyle\int\theta^{n}_{\mathbf{k}_{\parallel}}(z)dz\quad\text{ and }$
(25) $\displaystyle\delta_{np}$ $\displaystyle=$ $\displaystyle
d_{\mathbf{k}_{\parallel}}\int\theta^{n}_{\mathbf{k}_{\parallel}}(z)\theta^{p}_{\mathbf{k}_{\parallel}}(z)dz\,.$
(26)

With these definitions we create an orthonormal basis of layer-localized
functions, in analogy to the tight-binding formalism where site-localized
atomic wavefunctions constitute the basis for crystal properties. Besides, the
last equation is also the definition of a characteristic length
$d_{\mathbf{k}_{\parallel}}$, which can be associate to an effective
polarizable thickness of the layer.

Many ab initio codes assume a periodicity of the simulation cell in all
directions, so they employ the RS representation
$\left|{\mathbf{k}_{\parallel},k_{z}}\right)$. Let us now evaluate its
projection on the other two representations introduced above. Its projection
on the MS is a simple Fourier transformation on the $z$ variable:

$\left({\mathbf{k}_{\parallel},z}\left|{\mathbf{k}_{\parallel}^{\prime},k_{z}}\right.\right)=\frac{1}{\sqrt{NL}}e^{ik_{z}z}\delta_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}$
(27)

where $N$ is the number of layers, so $NL$ is the vertical size of the
crystal. Its projection on the LP is:

$\displaystyle\left({\mathbf{k}_{\parallel},n}\left|{\mathbf{k}_{\parallel}^{\prime},k_{z}}\right.\right)$
$\displaystyle=\sum_{\mathbf{G}_{\parallel}^{\prime\prime}}\int\left({\mathbf{k}_{\parallel},n}\left|{\mathbf{k}_{\parallel}^{\prime\prime},z}\right.\right)\left({\mathbf{k}_{\parallel}^{\prime\prime},z}\left|{\mathbf{k}_{\parallel}^{\prime},k_{z}}\right.\right)dz$
$\displaystyle=\sqrt{\frac{d_{\mathbf{k}_{\parallel}}}{NL}}\delta_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}\int\theta^{n}_{\mathbf{k}_{\parallel}}(z)e^{ik_{z}z}dz$
$\displaystyle=\sqrt{\frac{d_{\mathbf{k}_{\parallel}}}{NL}}\delta_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}e^{iq_{z}nL}\vartheta_{\mathbf{k}_{\parallel}}(k_{z})\,.$
(28)

In the last passage we used the definition of
$\theta^{n}_{\mathbf{k}_{\parallel}}(z)$ and used the fact that
$e^{iG_{z}nL}=1$ for all $n\in\mathbb{Z}$. This defines $\vartheta$ as the
Fourier transform of the $n=0$ profile function:

$\theta^{0}_{\mathbf{k}_{\parallel}}(z)=\frac{1}{2\pi}\int\vartheta_{\mathbf{k}_{\parallel}}(k_{z})e^{ik_{z}z}dk_{z}\approx\frac{1}{L}\sum_{G_{z}}\vartheta_{\mathbf{k}_{\parallel}}(G_{z})e^{iG_{z}z}\,,$
(29)

which is how the Fourier transform is actually calculated using standard
implementations of the discrete Fourier transform. This discrete approximation
is justified under the assumption that $L$ is large. To conclude this part, it
is useful to evaluate the Coulomb interaction in the LP representation:

$\displaystyle
v_{nm}(\mathbf{k}_{\parallel})=\left({\mathbf{k}_{\parallel},n}\right|v\left|{\mathbf{k}_{\parallel},m}\right)=$
$\displaystyle=\sum_{\mathbf{G}_{\parallel}^{\prime}}\int\int\left({\mathbf{k}_{\parallel},n}\left|{\mathbf{k}_{\parallel}^{\prime},z}\right.\right)v(\mathbf{k}_{\parallel}^{\prime},z,z^{\prime})\left(\left.{\mathbf{k}_{\parallel}^{\prime},z}\right|{\mathbf{k}_{\parallel},m}\right)dz\,dz^{\prime}$
$\displaystyle=d_{\mathbf{k}_{\parallel}}\int\int\theta^{n}_{\mathbf{k}_{\parallel}}(z)v(\mathbf{k}_{\parallel},z,z^{\prime})\theta^{m}_{\mathbf{k}_{\parallel}}(z^{\prime})dz\,dz^{\prime}.$
(30)

In the case $m=n$ it easy to show, through some simple changes of variables,
that;

$v_{nn}(\mathbf{k}_{\parallel})=d_{\mathbf{k}_{\parallel}}\int\int\theta^{0}(z)v(\mathbf{k}_{\parallel},z,z^{\prime})\theta^{0}(z^{\prime})dz\,dz^{\prime}\,.$
(31)

Comparing this result with definition (21) and remembering that
$\hat{v}_{nm}(\mathbf{k}_{\parallel}):=v_{nm}(\mathbf{k}_{\parallel})/d_{\mathbf{k}_{\parallel}}$,
one demonstrates that
$\hat{v}_{nn}(\mathbf{k}_{\parallel})\equiv\mathcal{V}^{\text{in}}(\mathbf{k}_{\parallel})$.

##  III Appendix C: Rewriting in reciprocal space

Our main result is the Dyson equation describing the array polarizability
$\hat{\chi}$ as generated by the inter-layer interactions between the single
layers, namely :

$\hat{\chi}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}nm}(\mathbf{q}_{\parallel},\omega)=\chi_{2D,\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}nm}(\mathbf{q}_{\parallel},\omega)+\sum_{\mathbf{G}_{\parallel}^{\prime\prime},p,s}\chi_{2D,\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime\prime}np}(\mathbf{q}_{\parallel},\omega)\mathcal{V}^{\text{off}}_{ps}(\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}^{\prime\prime})\hat{\chi}_{\mathbf{G}_{\parallel}^{\prime\prime}\mathbf{G}_{\parallel}^{\prime}sm}(\mathbf{q}_{\parallel},\omega)\,.$
(32)

In the latter equation,
$\chi_{2D,\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}nm}$ is defined
as the solution of the Dyson equation:

$\chi_{2D,nm}=\hat{\chi}^{0}_{nm}+\sum_{pq}\hat{\chi}^{0}_{np}\mathcal{V}^{\text{in}}_{pq}\,\chi_{2D,qm}$
(33)

resulting from the splitting of the Coulomb interaction into an intra-layer
contribution $\mathcal{V}^{\text{in}}_{nm}=\delta_{nm}\mathcal{V}^{\text{in}}$
and an inter-layer one. We also recall that $\hat{\chi}^{0}$ is a sum over
periodically repeated single-layer contributions, so
$\hat{\chi}^{0}_{nm}=\delta_{nm}\chi^{0}_{2D}$. This observation, together
with definition of $\mathcal{V}^{\text{in}}_{nm}$, allow us to simplify
equation (33) as:

$\chi_{2D,nm}=\delta_{nm}\left[\chi^{0}_{2D}+\chi^{0}_{2D}\mathcal{V}^{\text{in}}\chi_{2D}\right]\,,$
(34)

which then becomes diagonal in the layer-projected basis, and actually layer-
independent.

Let us now turn back to equation (32). Its solution is solved by inverting the
Dyson equation, which gives:

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle\hat{\chi}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime},mn}(\mathbf{q}_{\parallel},\omega)=$\mbox{}\hfil\\\
&\qquad=\sum_{\mathbf{G}_{\parallel}^{\prime\prime}}\mathcal{M}^{-1}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime\prime}nm}(\mathbf{q}_{\parallel},\omega)\chi_{2D\mathbf{G}_{\parallel}^{\prime\prime}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)\,,\end{split}$
(35)

with

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle\mathcal{M}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}nm}(\mathbf{q}_{\parallel},\omega)=$\mbox{}\hfil\\\
&=\delta_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}\delta_{nm}-\chi_{2D\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)\mathcal{V}^{\text{off}}_{nm}(\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}^{\prime})\,.\end{split}$
(36)

To make a connection between this result and the output of a standard ab
initio simulation $X_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)$, we
shall switch to the reciprocal-space representation
$\left|{\mathbf{k}_{\parallel},k_{z}}\right)$:

$\displaystyle
X_{\mathbf{G},\mathbf{G}^{\prime}}(\mathbf{q})=\left({\mathbf{k}_{\parallel},k_{z}}\left|{X}\right|{\mathbf{k}_{\parallel}^{\prime},k_{z}^{\prime}}\right)=$
$\displaystyle=\sum_{\mathbf{G}_{\parallel}^{\prime\prime}\mathbf{G}_{\parallel}^{\prime\prime\prime}}\sum_{m,n}\left({\mathbf{k}_{\parallel},k_{z}}\left|{\mathbf{k}_{\parallel}^{\prime\prime},m}\right.\right)X_{\mathbf{k}_{\parallel}^{\prime\prime}\mathbf{k}_{\parallel}^{\prime\prime\prime},mn}\left({\mathbf{k}_{\parallel}^{\prime\prime\prime},n}\left|{\mathbf{k}_{\parallel}^{\prime},k_{z}^{\prime}}\right.\right)\,.$

Using the projection (28) and remembering that :

$X_{\mathbf{k}_{\parallel}\mathbf{k}_{\parallel}^{\prime}mn}=\sqrt{d_{\mathbf{k}_{\parallel}}d_{\mathbf{k}_{\parallel}^{\prime}}}\hat{\chi}_{\mathbf{k}_{\parallel}\mathbf{k}_{\parallel}^{\prime}mn}\,,$

one finally arrives to:

$X_{\mathbf{G},\mathbf{G}^{\prime}}(\mathbf{q})=\frac{1}{L}\vartheta_{\mathbf{k}_{\parallel}}(k_{z})\hat{\chi}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},q_{z})\vartheta_{\mathbf{k}_{\parallel}^{\prime}}(k_{z}^{\prime})\,,$
(37)

where we have introduced the definition:

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle\hat{\chi}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},q_{z})=$\mbox{}\hfil\\\
&\qquad=\frac{1}{N}\sum_{m,n}e^{-iq_{z}mL}\hat{\chi}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime},mn}(\mathbf{q}_{\parallel})e^{iq_{z}nL}\,.\end{split}$
(38)

The latter quantity is basically a discrete Fourier transform of (35). Since
$\chi_{2D}$ is layer-independent, the solutions (35) and (36) projects simply
as:

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle\hat{\chi}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},q_{z},\omega)=$\mbox{}\hfil\\\
&\qquad=\sum_{\mathbf{G}_{\parallel}^{\prime\prime}}\mathcal{M}^{-1}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime\prime}}(\mathbf{q}_{\parallel},q_{z},\omega)\chi_{2D\mathbf{G}_{\parallel}^{\prime\prime}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)\end{split}$
(39)

and

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle\mathcal{M}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},q_{z},\omega)=$\mbox{}\hfil\\\
&\quad=\delta_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}-\chi_{2D\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)\mathcal{V}^{\text{off}}(\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel},q_{z})\,.\end{split}$
(40)

The last step to complete the connection with the reciprocal-space
representation is the expression of
$\mathcal{V}^{\text{off}}(\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel},q_{z})$.
Remembering its definition in the layer-projected basis (it vanishes for
$m=n$), one writes

$\mathcal{V}^{\text{off}}(\mathbf{k}_{\parallel},q_{z})=\frac{1}{N}\sum_{m}\sum_{n\neq
m}e^{-iq_{z}mL}\hat{v}_{mn}(\mathbf{k}_{\parallel})e^{iq_{z}nL}\;,$

and
$\hat{v}_{nm}(\mathbf{k}_{\parallel})=d_{\mathbf{k}_{\parallel}}^{-1}v_{nm}(\mathbf{k}_{\parallel})$
where $v_{nm}(\mathbf{k}_{\parallel})$ is defined in (30). Because of the
localization of the $\theta^{n}$ functions, if $n<m$ then $|z-z^{\prime}|>0$
and vice versa. It is therefore convenient to split $\sum_{n\neq
m}=\sum_{n<m}+\sum_{n>m}$. Next, inside the first sum, we change variables
$z-mL=\zeta$ and $z^{\prime}-nL=\zeta^{\prime}$, while inside the second one,
we put $z-mL=-\zeta$ and $z^{\prime}-nL=-\zeta^{\prime}$. Rearranging all
terms, using the parity of the $\theta^{0}(z)$ functions and the identity
$\hat{v}_{00}=\mathcal{V}^{\text{in}}$, we get

$\begin{split}&\mathcal{V}^{\text{off}}(\mathbf{k}_{\parallel},q_{z})=\\\
&\qquad=\mathcal{V}^{\text{in}}(\mathbf{k}_{\parallel})\frac{1}{N}\sum_{m}\left[\sum_{n<m}e^{-(iq_{z}+|\mathbf{k}_{\parallel}|)(m-n)L}+\right.\\\
&\left.\qquad\quad+\sum_{n>m}e^{-(iq_{z}-|\mathbf{k}_{\parallel}|)(m-n)L}\right]\,.\end{split}$

Since both sums depend only on the difference $p=m-n$, we can change the index
and cast them in the form $\sum_{p=1}^{\infty}e^{-(|\mathbf{k}_{\parallel}|\pm
iq_{z})Lp}=\left(e^{(|\mathbf{k}_{\parallel}|\pm iq_{z})L}-1\right)^{-1}$. The
resulting expression for the inter-layer Coulomb interaction reads

$\begin{split}&\mathcal{V}^{\text{off}}(\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel},q_{z})=\mathcal{V}^{\text{in}}(\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel})\times\\\
&\times\left(\frac{1}{e^{(|\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}|+iq_{z})L}-1}+\frac{1}{e^{(|\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}|-iq_{z})L}-1}\right)\,,\end{split}$
(41)

which, evaluated at $q_{z}=0$, gives the result reported in the main text.

##  IV Appendix D: The dielectric function

###  IV.1 In the layer-projected representation

The definition of the microscopic dielectric function in terms of the total
and external fields is:

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle
U^{\text{tot}}_{\mathbf{G}_{\parallel}}(\mathbf{q}_{\parallel},\omega,z)=$\mbox{}\hfil\\\
&=\sum_{\mathbf{G}_{\parallel}^{\prime}}\int\epsilon^{-1}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega,z,z^{\prime})U^{\text{ext}}_{\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega,z^{\prime})dz^{\prime}\,.\end{split}$
(42)

We split now the total potential into an external and an induced part, and
assume the latter to be the classical Hartree potential, then the equation
above leads to the definition:

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle\epsilon^{-1}_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega,z,z^{\prime})=\delta_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}\delta(z-z^{\prime})+$\mbox{}\hfil\\\
&\qquad+\int
v(\mathbf{k}_{\parallel},z,z^{\prime\prime})\chi_{\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega,z^{\prime\prime},z^{\prime})dz^{\prime\prime}\;,\end{split}$
(43)

where $\chi$ is defined as in (22).

We introduce now the macroscopic average of the fields, i.e. the
$\mathbf{G}_{\parallel}=0$ component projected on the single-layer. To this
aim, we refer to the projection (24), where we take $n=0$ and $N=1$ since we
are working on the single-layer framework. Then, for $\alpha=$(tot) or (ext),

$U^{\alpha}_{M}(\mathbf{q}_{\parallel},\omega)=\sqrt{d_{\mathbf{q}_{\parallel}}}\int\theta^{0}_{\mathbf{q}_{\parallel}}(z)U^{\alpha}_{0}(\mathbf{q}_{\parallel},\omega,z)dz\,,$
(44)

and we define
$U^{\text{tot}}_{M}(\mathbf{q}_{\parallel},\omega)=U^{\text{ext}}_{M}(\mathbf{q}_{\parallel},\omega)/\epsilon_{M}(\mathbf{q}_{\parallel},\omega)$.
We can now plug (42) into (44) and assume that $U^{\text{ext}}$ is
macroscopic, which means that the only non-vanishing term is the
$\mathbf{G}_{\parallel}=0$ one. From the definitions (43) and (22), we get the
relation:

$U^{\text{tot}}_{M}(\mathbf{q}_{\parallel},\omega)=\sqrt{d_{\mathbf{q}_{\parallel}}}\int\int\theta^{0}_{\mathbf{q}_{\parallel}}(z)\left[\delta(z-z^{\prime})+\int
v(\mathbf{q}_{\parallel},z,z^{\prime})\theta^{0}_{\mathbf{q}_{\parallel}}(z^{\prime\prime})\chi_{2D,00}(\mathbf{q}_{\parallel},\omega)\theta^{0}_{\mathbf{q}_{\parallel}}(z^{\prime})dz^{\prime\prime}\right]U^{\text{ext}}_{0}(\mathbf{q}_{\parallel},\omega)dz^{\prime}dz\,,$
(45)

which eventually leads to the single layer macroscopic dielectric function:

$\frac{1}{\epsilon_{M}(\mathbf{q}_{\parallel},\omega)}=1+\mathcal{V}^{\text{in}}(\mathbf{q}_{\parallel})\chi_{2D,00}(\mathbf{q}_{\parallel},\omega)\,;$
(46)

###  IV.2 In the reciprocal-space representation

Let us now repeat the derivation in the
$\left|{\mathbf{k}_{\parallel},k_{z}}\right)$ representation. This will allow
us to show that our approach reproduces some results obtained by other authors
in the sharp slab approximation Latini  _et al._ (2015), but actually permits
to go beyond them generalizing their method. Similarly to what done xbefore,
we use (28) with $n=0$ and $N=1$. Furthermore we set from the beginning
$\mathbf{G}_{\parallel}=0$, we assume $U^{\text{ext}}$ to be macroscopic, and
we drop the frequency variable for notational convenience. We get:

$\displaystyle U^{\text{tot}}_{M}(\mathbf{q}_{\parallel})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{d_{\mathbf{q}_{\parallel}}}{L}}\int\sum_{G_{z}}\vartheta_{\mathbf{q}_{\parallel}}(k_{z})U^{\text{tot}}_{(0,G_{z})}(\mathbf{q}_{\parallel},q_{z})dq_{z}$
$\displaystyle U^{\text{ext}}_{M}(\mathbf{q}_{\parallel})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{d_{\mathbf{q}_{\parallel}}}{L}}\int\vartheta_{\mathbf{q}_{\parallel}}(q_{z})U^{\text{ext}}_{(0,0)}(\mathbf{q}_{\parallel},q_{z})dq_{z}\,.$

We can do the additional assumption that the integrands depend weakly on
$q_{z}$ and hence fix $q_{z}=0$ and replace $\int f(q_{z})dq_{z}\approx
Lf(0)/2\pi$ in both expressions. We get the two definitions:

$\displaystyle U^{\text{tot}}_{M}(\mathbf{q}_{\parallel})$
$\displaystyle\approx$
$\displaystyle\frac{\sqrt{d_{\mathbf{q}_{\parallel}}L}}{2\pi}\sum_{G_{z}}\vartheta_{\mathbf{q}_{\parallel}}(G_{z})U^{\text{tot}}_{(0,G_{z})}(\mathbf{q}_{\parallel},0)$
(47) $\displaystyle U^{\text{ext}}_{M}(\mathbf{q}_{\parallel})$
$\displaystyle\approx$
$\displaystyle\frac{\sqrt{d_{\mathbf{q}_{\parallel}}L}}{2\pi}U^{\text{ext}}_{(0,0)}(\mathbf{q}_{\parallel},0)\,,$
(48)

where we have used the fact that
$\vartheta_{\mathbf{q}_{\parallel}}(0)=1\,\forall\mathbf{q}_{\parallel}$ as a
consequence of its normalization in real space.

In the $\left|{\mathbf{k}_{\parallel},k_{z}}\right)$ representation, the
microscopic dielectric function is defined according to the relation
$U^{\text{tot}}_{\mathbf{G}}(\mathbf{q})=\sum_{\mathbf{G}^{\prime}}\epsilon^{-1}_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q})U^{\text{ext}}_{\mathbf{G}^{\prime}}(\mathbf{q})$.
Inserting this definition into (47), using (48), and using the fact that
$U^{\text{ext}}$ is macroscopic, one gets:

$U^{\text{tot}}_{M}(\mathbf{q}_{\parallel})=\sum_{G_{z}}\vartheta_{\mathbf{q}_{\parallel}}(G_{z})\epsilon^{-1}_{(0,G_{z}),(0,0)}(\mathbf{q}_{\parallel})U^{\text{ext}}_{M}(\mathbf{q}_{\parallel})\,,$

from which one defines:

$\frac{1}{\epsilon_{M}(\mathbf{q}_{\parallel})}=\sum_{G_{z}}\vartheta_{\mathbf{q}_{\parallel}}(G_{z})\epsilon^{-1}_{(0,G_{z}),(0,0)}(\mathbf{q}_{\parallel})\,.$
(49)

This result is a generalization of the _quasi 2D_ dielectric function derived
by Latini et al. Latini  _et al._ (2015). Their result can be retrieved within
the _sharp slab_ approximation, which in the reciprocal space gives
$\vartheta_{\mathbf{q}_{\parallel}}(G_{z})=(G_{z}\frac{d}{2})^{-1}\sin(G_{z}\frac{d}{2})$.

Let us now show that (49) leads to the same expression we derived in the main
text. First of all, let us recall that the dielectric function in the
$\left|{\mathbf{k}_{\parallel},k_{z}}\right)$ representation reads:

$\epsilon^{-1}_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)=\delta_{\mathbf{G}\mathbf{G}^{\prime}}+v(\mathbf{q}+\mathbf{G})\chi_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)\,$
(50)

and that the single layer polarizability is:

$\chi_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)=\frac{1}{L}\vartheta_{\mathbf{k}_{\parallel}}(k_{z})\chi_{2D,\mathbf{G}_{\parallel}\mathbf{G}_{\parallel}^{\prime}}(\mathbf{q}_{\parallel},\omega)\vartheta_{\mathbf{k}_{\parallel}}(k_{z}^{\prime})\,,$
(51)

in analogy to (37). Actually, in equation (37), $L$ is the period of the
periodic array, and $NL$ the height of the crystal. In expression (51) $L$ has
to be interpreted as the size of the box of integration, and $N=1$. However,
in practice, it still corresponds to the height of the simulation cell. In any
case, at the end of the calculation this factor will cancel out, in agreement
with the physical requirement that the single-layer dielectric function does
not depend on the size of the simulation cell.

Using (50) and (51), and remembering that
$\vartheta_{\mathbf{k}_{\parallel}}(0)=1$, equation (49) becomes:

$\frac{1}{\epsilon_{M}(\mathbf{q}_{\parallel},\omega)}=1+\frac{\chi_{2D\,00}(\mathbf{q}_{\parallel},\omega)}{L}\sum_{G_{z}}\vartheta_{\mathbf{q}_{\parallel}}^{2}(G_{z})v(\mathbf{q}_{\parallel},G_{z})\,.$

Finally, to complete the demonstration, we shall evaluate the term
$\sum_{G_{z}}$. This is done by projecting successively on the different
representations (remember $N=1$). We get:

$\displaystyle\sum_{G_{z}}\vartheta_{\mathbf{q}_{\parallel}}^{2}(G_{z})v(\mathbf{q}_{\parallel},G_{z})=$
$\displaystyle=$
$\displaystyle\frac{L}{d_{\mathbf{q}_{\parallel}}}\sum_{G_{z}G_{z}^{\prime}}\left({\mathbf{q}_{\parallel},n=0}\left|{\mathbf{q}_{\parallel},G_{z}}\right.\right)\left({\mathbf{q}_{\parallel},G_{z}}\right|v\left|{\mathbf{q}_{\parallel},G_{z}^{\prime}}\right)\left(\left.{\mathbf{q}_{\parallel},G_{z}^{\prime}}\right|{\mathbf{q}_{\parallel},n=0}\right)$
$\displaystyle=$
$\displaystyle\frac{L}{d_{\mathbf{q}_{\parallel}}}\sum_{G_{z}G_{z}^{\prime}}\int\int\left({\mathbf{q}_{\parallel},0}\left|{\mathbf{q}_{\parallel},z}\right.\right)\left({\mathbf{q}_{\parallel},z}\left|{\mathbf{q}_{\parallel},G_{z}}\right.\right)\left({\mathbf{q}_{\parallel},G_{z}}\right|v\left|{\mathbf{q}_{\parallel},G_{z}^{\prime}}\right)\left({\mathbf{q}_{\parallel},G_{z}^{\prime}}\left|{\mathbf{q}_{\parallel},z^{\prime}}\right.\right)\left({\mathbf{q}_{\parallel},z^{\prime}}\left|{\mathbf{q}_{\parallel},0}\right.\right)dzdz^{\prime}\,.$

We use definitions (24), (27) and (28) and arrive to the expression;

$\sum_{G_{z}}\vartheta_{\mathbf{q}_{\parallel}}^{2}(G_{z})v(\mathbf{q}_{\parallel},G_{z})=\sum_{G_{z}G_{z}^{\prime}}\int\int\theta^{0}_{\mathbf{q}_{\parallel}}(z)e^{iG_{z}z}\frac{4\pi}{|\mathbf{q}_{\parallel}+\hat{z}G_{z}|^{2}}\delta_{G_{z}G_{z}^{\prime}}e^{-iG_{z}^{\prime}z^{\prime}}\theta^{0}_{\mathbf{q}_{\parallel}}(z^{\prime})dz\,dz^{\prime}\,.$
(52)

If we notice that
$v(\mathbf{q}_{\parallel},z,z^{\prime})=\frac{1}{L}\sum_{G_{z}}v(\mathbf{q}_{\parallel}+\hat{z}G_{z})e^{iG_{z}(z-z^{\prime})}$
is a Fourier transform at $q_{z}=0$ and $N=1$, then we are left with L times
the definition (21), so the result becomes:

$\sum_{G_{z}}\vartheta_{\mathbf{q}_{\parallel}}^{2}(G_{z})v(\mathbf{q}_{\parallel},G_{z})=L\mathcal{V}^{\text{in}}(\mathbf{q}_{\parallel})\,,$
(53)

which completes the demonstration.

##  V Appendix E: Computational details

In this section we give the details of the simulations of the hBN monolayer.

All ground state calculations have been carried out with the ABINIT simulation
package Gonze  _et al._ (2016). The parameters of the simulation cell are
$a=2.5$ Å for the in-plane lattice parameter and $L=15$ Å or $L=30$ Å for the
height of the simulation cell. In all cases, the local density approximation
(LDA) was used to model the exchange-correlation potential of the Kohn-Sham
Hamiltonian and the cutoff energy defining the basis set was 30 Ha. The
k-point grid used to sample the electronic density was $15\times 15\times 1$,
instead for the wavefunction calculations, necessary to compute $X^{0}$, a
denser grid of $36\times 36$ k-points in the $(x,y)$-plane was employed. Both
grids were centered in $\Gamma$. No Coulomb truncation has been used in the
ground state calculations.

The wavefunctions have been successively passed to the DP dp  simulation
package with which we computed only the irreducible polarizability $X^{0}$ of
the periodic array. The sum over all transitions included 30 bands in the
$L=15$ Å calculations, and 60 in the $L=30$ Å ones. The cutoff energy for the
representation of the matrix elements was 400 eV, and it was 100 eV for the
dimension of the $X_{\mathbf{G}\mathbf{G}^{\prime}}$. No scissor operator has
been applied, so the resulting gap was the LDA one.

All the RPA algorithms connecting $X^{0}$ to $X$ and $\tilde{\epsilon}_{M}$,
with and without Coulomb truncation, as well as all the manipulations proper
to our method have been implemented in a home-made code which takes $X^{0}$ in
input.

## References

  * Cudazzo  _et al._ (2011) P. Cudazzo, I. V. Tokatly, and A. Rubio, Physical Review B - Condensed Matter and Materials Physics 84, 085406 (2011).
  * Sponza  _et al._ (2016) L. Sponza, J. Goniakowski, and C. Noguera, Physical Review B 93, 195435 (2016).
  * Tian  _et al._ (2020) T. Tian, D. Scullion, D. Hughes, L. H. Li, C.-J. Shih, J. Coleman, M. Chhowalla, and E. J. G. Santos, Nano Letters 20, 841 (2020).
  * Ismail-Beigi (2006) S. Ismail-Beigi, Physical Review B - Condensed Matter and Materials Physics 73, 1 (2006).
  * Rozzi  _et al._ (2006) C. A. Rozzi, D. Varsano, A. Marini, and E. K. U. Gross, Physical Review B , 1 (2006).
  * Hüser  _et al._ (2013) F. Hüser, T. Olsen, and K. S. Thygesen, Physical Review B 87, 235132 (2013).
  * Cudazzo  _et al._ (2016) P. Cudazzo, L. Sponza, C. Giorgetti, L. Reining, F. Sottile, and M. Gatti, Physical Review Letters 116, 066803 (2016).
  * Nazarov  _et al._ (2014) V. U. Nazarov, F. Alharbi, T. S. Fisher, and S. Kais, Physical Review B - Condensed Matter and Materials Physics 89, 195423 (2014).
  * Nazarov (2015) V. U. Nazarov, New Journal of Physics 17, 73018 (2015).
  * Latini  _et al._ (2015) S. Latini, T. Olsen, and K. S. Thygesen, Physical Review B 92, 1 (2015).
  * Thygesen (2017) K. S. Thygesen, 2D Materials 4, 022004 (2017).
  * Meckbach  _et al._ (2018) L. Meckbach, T. Stroucken, and S. W. Koch, Physical Review B 97, 035425 (2018).
  * Gjerding  _et al._ (2020) M. N. Gjerding, L. S. R. Cavalcante, A. Chaves, and K. S. Thygesen, The Journal of Physical Chemistry C 124, 11609 (2020).
  * Andersen  _et al._ (2015) K. Andersen, S. Latini, and K. S. Thygesen, Nano Letters 15, 4616 (2015).
  * Cavalcante  _et al._ (2018) L. S. R. Cavalcante, A. Chaves, B. Van Duppen, F. M. Peeters, and D. R. Reichman, Physical Review B 97, 125427 (2018).
  * Note (1) Note that in this notation, the continuous part never changes. So if $\mathbf{k}_{\parallel}=\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}$, then $\mathbf{k}_{\parallel}^{\prime}=\mathbf{q}_{\parallel}+\mathbf{G}_{\parallel}^{\prime}$ and analogously $k_{z}^{\prime}=q_{z}+G_{z}^{\prime}$.
  * Qiu  _et al._ (2016) D. Y. Qiu, F. H. da Jornada, and S. G. Louie, Physical Review B 93, 235435 (2016).
  * Gonze  _et al._ (2016) X. Gonze, F. Jollet, A. A. F., D. Adams, B. Amadon, T. Applencourt, C. Audouzec, J.-M. Beuken, J. Bieder, A. Bokhanchuk, E. Bousquet, B. F., D. Caliste, M. Côté, F. Dahm, F. Da Pieve, M. Delaveau, M. Di Gennaro, B. Dorado, C. Espejo, G. Geneste, L. Genovese, A. Gerossier, M. Giantomassi, Y. Gillet, D. Hamann, L. He, G. Jomard, J. Laflamme Janssen, S. Le Roux, A. Levitt, A. Lherbier, F. Liu, I. Lukačević, A. Martin, C. Martins, M. Oliveira, S. Poncé, Y. Pouillon, T. Rangel, G.-M. Rignanese, A. Romero, B. Rousseau, O. Rubel, A. Shukri, S. M, M. Torrent, M. Van Setten, B. Van Troeye, M. Verstraete, D. Waroquiers, J. Wiktor, B. Xu, A. Zhou, and J. Zwanziger, Comput. Phys. Commun. 205, 106–310 (2016).
  * (19) https://etsf.polytechnique.fr/Software/Ab_Initio.


# Study on Aspect Ratio Variability toward Robustness of Vision Transformer-
based Vehicle Re-identification

Mei Qiu, Lauren Christopher, and Lingxi Li Mei Qiu, Lauren Christopher, and
Lingxi Li are with the Department of Electrical and Computer Engineering,
Purdue University in Indianapolis, 723 West Michigan Street, SL-160,
Indianapolis, Indiana 46202, USA. Emails:
meiqiu<EMAIL_ADDRESS>
###### Abstract

Vision Transformers (ViTs) have excelled in vehicle re-identification (ReID)
tasks. However, non-square aspect ratios of image or video input might
significantly affect the re-identification performance. To address this issue,
we propose a novel ViT-based ReID framework in this paper, which fuses models
trained on a variety of aspect ratios. Our main contributions are threefold:
(i) We analyze aspect ratio performance on VeRi-776 and VehicleID datasets,
guiding input settings based on aspect ratios of original images. (ii) We
introduce patch-wise mixup intra-image during ViT patchification (guided by
spatial attention scores) and implement uneven stride for better object aspect
ratio matching. (iii) We propose a dynamic feature fusing ReID network,
enhancing model robustness. Our ReID method achieves a significantly improved
mean Average Precision (mAP) of 91.0% compared to the the closest state-of-
the-art (CAL) result of 80.9% on VehicleID dataset.

##  I INTRODUCTION

Figure 1:  Aspect ratio distribution of images in training datasets from ReID
benchmark datasets, VeRi-776 and VehicleID, varies significantly. These
datasets show that a substantial portion of the images are non-square.

A fundamental task in intelligent transportation systems is vehicle re-
identification (Re-ID): identifying vehicles across multiple non-overlapping
cameras [1]. Despite its significance, vehicle Re-ID encounters challenges due
to variations in vehicle appearance across different viewpoints, poses,
illumination, and backgrounds. Deep learning models face the challenge of
extracting discriminative features resistant to viewpoint variations [2, 3].
Both global and local features are essential for generating robust
representations for vehicle pairs [4, 5, 6].

Numerous benchmark datasets contain images from real-world surveillance
scenes, including VeRi-776 [2], PKU-VD [7], VehicleID [8], Vehicle-1M [9], and
VERI-Wild [10], crucial for vehicle Re-ID AI development [11]. State-of-the-
art models leverage the self-attention mechanism. Vision transformers (ViTs)
have demonstrated superior ability in capturing discriminative details
compared to previous CNN-based methods [12]. Additionally, images in different
datasets exhibit various aspect ratios, as illustrated in Fig. 1. VeRi-776 and
VehicleID display different size and shape distributions, posing a significant
challenge in model training.

However, unlike CNN-based models, which can handle varying aspect ratios to
some extent due to their translation invariance and local receptive fields,
vision transformers consider the entire image as a sequence of patches. The
fixed patch size and sequential nature of transformers necessitate careful
consideration of how input images are resized and cropped [13, 14, 15, 16].
Early implementations of vision transformers adopted resizing strategies from
CNNs, often distorting the original aspect ratio and potentially compromising
performance, particularly in tasks reliant on object shape and scale.

To address this, subsequent studies explored padding strategies, adaptive
post-patch extraction, trainable resizing networks, aspect-ratio aware
attention mechanisms, and multi-scale/multi-aspect training [17, 18, 19, 20,
13, 21]. However, these approaches entail computational burdens, data
requirements, and optimization challenges.

In summary, gaps persist in applying Vision Transformers (ViTs) to vehicle Re-
ID. Optimization of scaling and resizing strategies, understanding aspect
ratio effects, and exploring data augmentations such as mixup [22] at the
patch level for ViTs in multi-aspect ratio scenarios are needed. Intra-image
mixup can potentially enhance the model’s ability to learn detailed features,
particularly when the whole image is distorted by resizing with unsatisfactory
aspect ratios. TransReID [12] introduces a jigsaw patch module (JPM) to
enhance feature robustness and discrimination ability, representing a step
towards addressing these challenges. However, this intra-mixup occurs at the
feature level; its potential at the pixel level needs to be comprehensively
explored.

Figure 2:  (Left) Existing Method: Image size is typically fixed and set to a
single size with a square shape. (Right) Our Method: Combined Vision
Transformer (ViT)-based ReID model that dynamically fuses features extracted
from multiple models. Each model is trained on a fixed size and aspect ratio.
Figure 3: The structure of each individual model is designed to adapt to the
dataset’s size and aspect ratio distribution. During the patchification
process, the stride size in the horizontal and vertical directions is
dynamically determined based on the input object’s aspect ratio. Subsequently,
a Patch Mixing (PM) module shuffles and mixes patches from the same image
using an attention-guided strategy. Any Vision Transformer (ViT)-based
architecture can be chosen as the backbone. In this study, we select ViT/B-16.
The features extracted from ViTs are used for the vehicle ReID downstream
task.

In this work, we are the first to propose that aspect ratio is a key factor
affecting vehicle Re-ID performance and the robustness of feature learning in
Vision Transformers (ViTs). We conduct a series of novel, scientific, and
comprehensive experiments to explore the effects of various aspect ratios on
ViT-based ReID.

To enhance the model’s generality to various aspect ratios in the input, we
dynamically fuse features extracted from several models trained on images with
different aspect ratios, as shown in Fig. 2. Additionally, when training a
single model, we propose a novel intra-image patch mixup (PM) data
augmentation method to improve the model’s learning ability on details and
mitigate overfitting during training. Furthermore, to mitigate the distortion
caused by unsatisfactory resizing, we employ an uneven stride strategy in the
patchify step.

The key contributions of this work are:

  * •

Proposed using aspect ratio as a critical factor impacting vehicle Re-ID
performance and ViT’s feature learning robustness.

  * •

Conducted novel and comprehensive experiments to explore the effects of
various aspect ratios on ViT-based ReID.

  * •

Introducted dynamic feature fusion from models trained on different aspect
ratios to enhance model generality.

  * •

Proposed an intra-image patch mixup (PM) data augmentation method to improve
model learning ability and prevent overfitting.

  * •

Implementated an uneven stride strategy to reduce distortion in unsatisfactory
resizing.

##  II Related Works.

Vision Transformer. The Vision Transformer (ViT) is an adaptation of the
transformer architecture from natural language processing (NLP) tasks [23] to
computer vision. Dosovitskiy et al. introduced the ViT [24], the first to
demonstrate that a pure transformer applied directly to sequences of image
patches can excel in large-scale image classification tasks. Within ViTs, the
multi-head self-attention mechanism enables the model to capture diverse
dependencies in images, including shapes, textures, and contextual
relationships between objects in parallel, facilitating efficient learning of
richer data representations. ViTs and their variants, such as DeiT [25], Swin
Transformer [26], PVT [27], and CPVT [28], prove beneficial across tasks
including image classification, object detection, and Re-ID [29].

Vehicle Re-Identification. As discussed, vehicle Re-ID constitutes the primary
focus of this study. Numerous works leveraging deep learning have achieved
notable performance across various public vehicle Re-ID benchmarks. These
works often utilize either CNN backbones [30, 31, 32] or ViT backbones [12,
33, 34, 35]. Common loss functions employed in deep vehicle Re-ID network
training include cross-entropy loss (ID loss) [36], triple loss [37], and
contrastive loss [38].

Figure 4: Displayed are examples from the VehicleID (first two columns) and
VeRi-776 (last two columns) test datasets, illustrating the effects of the
intra-image patch mixup (PM) data augmentation method. This technique blends
various parts of an image based on attention-driven distances, increasing
image complexity to enhance model robustness and reduce overfitting. The top
row presents images without the PM module, while the bottom row features
images processed with the PM module.

##  III Method

Model Structure. We train one model for each major aspect ratio. The
determination of how many models are needed for training with various aspect
ratios is learned from the data, as illustrated in Fig. 3. For each model with
a fixed aspect ratio, the resized input images undergo augmentation using a
Patch Mixup (PM) module. Subsequently, a chosen Vision Transformer (ViT)
backbone and its pre-trained weights are used for initialization. Features
extracted from the last transformer layer are utilized for the ReID task.

Image Input with Adaptive Size and Shape. We employ a statistical method to
estimate the dataset’s size and aspect ratios. Starting with the mean or
median value of the dataset’s size is a suitable approach. For aspect ratios,
clustering methods such as K-means can be used to generate several clusters
representing the aspect ratios present in the original dataset. Subsequently,
a combination of size and aspect ratio guides the resizing of images for the
training of each specific model. Each model is trained on a fixed size and
aspect ratio.

Patchification with uneven stride. To enhance the model’s ability to learn
spatial relationships from image data, we incorporate uneven strides in
different dimensions, based on the aspect ratio. Specifically, in our
approach, the stride sizes $s_{h}$ and $s_{w}$ can vary, with the stride size
being smaller in the shorter dimension compared to the longer dimension. We
maintain a fixed patch size $p$ of 16. The total patch number $n$ is
calculated as
$\left(\frac{H-p}{s_{h}}+1\right)\times\left(\frac{W-p}{s_{w}}+1\right)$.

Patch Mixup intra image module. We propose a novel intra-image data
augmentation method where each patch of an image has a probability to mix
adaptively with another randomly chosen patch from the same image. Mixup
weights are determined based on the spatial distance between these two patches
in pixels, with closer patches assigned higher mixup weights. Let $A$ denote
the original patches with $n$ patches, and their positions are denoted by the
coordinates of their top-left corners $(x_{i},y_{i})$ for the $i$-th patch.
After patch partitioning, all patch indexes are fixed, with each patch having
a unique index ranging from $0$ to $n-1$. The Euclidean distance between the
centers of two patches indexed by $i$ and $j$ can be computed as:

$d(i,j)=\sqrt{(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}}$ (1)

The patch distance matrix $D$ for all patches can then be represented as:

$D=\begin{bmatrix}d(1,1)&d(1,2)&\cdots&d(1,n)\\\
d(2,1)&d(2,2)&\cdots&d(2,n)\\\ \vdots&\vdots&\ddots&\vdots\\\
d(n,1)&d(n,2)&\cdots&d(n,n)\end{bmatrix}$ (2)

Given a distance matrix $D$ where each element $D_{ij}$ represents the
distance between patch $i$ and patch $j$, the attention scores matrix $S$ can
be computed as:

$S=\frac{1}{1+D*p}$ (3)

In this formula, the attention score $S_{ij}$ increases as the distance
$D_{ij}$ decreases, indicating that closer patches have more influence on each
other. After shuffling patches, we generate a new patch-to-patch mixup weights
matrix $S^{\prime}$ based on the patches’ indexes, as shown in Equation 4,
where each value of $S^{\prime}$ is fetched from $S$ using the indexes of the
patches:

$S^{\prime}=S[A_{k},B_{k}^{\prime}]$ (4)

Here, $k$ is the patch index of a patch from the original patches $A$, and
$k^{\prime}$ is the patch index of the selected patch from the shuffled
patches $B$ for mixup processing. Assuming the original patches $A$ are
generated using non-overlapping patchification, and its patch distance matrix
$D$ is calculated based on Equation 1 and 2, with $D=16\times D$ as the patch
is a fixed-size square. The original spatial attention scores and mixup
weights matrix based on the shuffled patches are calculated using Equations 3
and 4. The final patches of this image $F$ are given by: $F=S^{\prime}\times
B+(1-B)\times A$. Several image samples before and after Patch Mixing (PM)
data augmentation are shown in Fig. 4, where it can be observed that after PM,
detailed information from different parts of the same vehicle is mixed. The
ablation study result is shown in Table I. The pseudocode of this method is
shown in Algorithm 1.

Algorithm 1 Patch Mixing guided by Spatial Attention Scores in Vision
Transformer

1:An input image $I$ of dimensions $H\times W\times C$, a patch size $P$, and
stride $S$ for overlapping.

2:$B\leftarrow\text{batch size of }I$

3:$N\leftarrow\text{number of patches along height}$

4:$M\leftarrow\text{number of patches along width}$

5:$D\leftarrow\text{patch distance matrix of size }(N\times M)\times(N\times
M)$

6:$S\leftarrow\text{compute attention scores from }D$

7:for $b=1$ to $B$ do

8: $P_{b}\leftarrow\text{extract patches from image }I_{b}\text{ with
overlap}$

9: $\pi\leftarrow\text{random permutation of }\\{1,\ldots,N\times M\\}$

10: $P^{\prime}_{b}\leftarrow P_{b}[\pi]\text{ shuffled patches using
permutation }\pi$

11: $S^{\prime}_{b}\leftarrow S[\pi]\text{ adjust attention scores with }\pi$

12: for $i=1$ to $N\times M$ do

13: $\lambda\leftarrow S^{\prime}_{b_{i,\pi(i)}}\text{ attention weight for
current patch}$

14: $P^{\prime\prime}_{b_{i}}\leftarrow(1-\lambda)\cdot P_{b_{i}}+\lambda\cdot
P^{\prime}_{b_{\pi(i)}}$

15: end for

16: $O_{b}\leftarrow\text{reconstruct output from mixed patches
}P^{\prime\prime}_{b}$

17:end for

18:return $O$ $\triangleright$ Return the batch of mixed images

ReID Task. We optimize the network by constructing ID loss and triplet loss
for global features. The ID loss $L_{ID}$ is the cross-entropy loss without
label smoothing. For a triplet set $\\{a,p,n\\}$, the triplet loss $L_{T}$
with a soft-margin. The triplet loss $L_{T}$ with soft-margin is defined as:

$L_{T}=\log(1+\exp(||\mathbf{f}_{a}-\mathbf{f}_{p}||^{2}-||\mathbf{f}_{a}-\mathbf{f}_{n}||^{2}))$

where:

  * •

$\mathbf{f}_{a}$, $\mathbf{f}_{p}$, and $\mathbf{f}_{n}$ are the feature
embeddings of the anchor, positive, and negative samples, respectively.

  * •

$||\cdot||^{2}$ denotes the squared Euclidean distance.

Inference Phase. As shown in Fig. 2, we use a Dynamic Feature Fusing strategy
during inference to fuse features from models trained on multiple inputs. The
output feature vector can be represented as:

$\textbf{F\\_out}=w_{1}\cdot\textbf{f}_{1}+w_{2}\cdot\textbf{f}_{2}+w_{3}\cdot\textbf{f}_{3}$

or as a weighted concatenation:

$\textbf{F\\_out}=\left[w_{1}\cdot\textbf{f}_{1};w_{2}\cdot\textbf{f}_{2};w_{3}\cdot\textbf{f}_{3}\right]$

where $w_{1}$, $w_{2}$, and $w_{3}$ are weights. For the adaptive weight
assignment $w$:

$w=\begin{cases}1.3,&\text{if }|\text{model\\_ar}-\text{image\\_ar}|\leq
0.3\\\ 1.0,&\text{if }0.3<|\text{model\\_ar}-\text{image\\_ar}|\leq 0.6\\\
0.9,&\text{otherwise}\end{cases}$

We investigated the ReID performance by summing or concatenating features, and
the results are shown in Table I.

##  IV Experiments

###  IV-A Experiment Settings

Dataset. Two popular vehicle Re-ID benmark datasets: VeRi-776 [2] and
VehicleID [8], are used in this work.

Evaluation Metrics. In line with the Re-ID field, we employ mAP (mean Average
Precision) and CMC (Cumulative Matching Characteristic) as our performance
estimation metrics. Specifically, we consider Rank-1 (R1), Rank-5 (R5), and
Rank-10 (R10) accuracy.

Dataset |  | Image Size/aspect ratio
---
(w/PM or w/o PM)
Stride Size | Fusing Strategy | mAP(%) | R1(%) | R5(%) | R10(%)
| [224,224]/1.0 (w/o) | [16,16] | - | 74.9 | 95.0 | 97.6 | 98.6
| [224,224]/1.0 (w/) | [16,16] | - | 77.8 | 96.1 | 98.8 | 99.1
| [224,212]/0.95 (w/o) | [16,12] | - | 75.5 | 95.2 | 97.9 | 98.5
| [224,212]/0.95 (w/) | [16,12] | - | 73.2 | 93.5 | 97.4 | 98.4
| [224,298]/1.33 (w/o) | [12,16] | - | 78.5 | 96.5 | 98.0 | 99.0
| [224,298]/1.33 (w/) | [12,16] | - | 79.4 (4.5$\uparrow$) | 96.1 | 98.2 | 99.0
|  | Fusing three models
---
(w/o)
- | Weighted Sum | 78.1 | 96.1 | 98.1 | 98.7
|  Fusing three models (w/) | - | Weighted Sum | 81.4(6.5$\uparrow$) | 97.0(2.0$\uparrow$) | 98.5 | 99.1
|  | Fusing three models
---
(w/o)
- |  | Weighted
---
Concatenate
78.0 | 95.9 | 98.1 | 98.6
VeRi776 |  | Fusing three models
---
(w/)
- |  | Weighted
---
Concatenate
81.4 | 97.0 | 98.5 | 99.0
| [384,384]/1.0 (w/o) | [16,16] | - | 90.2 | 85.0 | 97.5 | 99.3
| [384,384]/1.0 (w/) | [16,16] | - | 90.8(0.6$\uparrow$) | 85.7 | 97.8 | 99.0
| [384,308]/0.80 (w/o) | [16,12] | - | 89.9 | 84.1 | 97.6 | 99.2
| [384,308]/0.80 (w/) | [16,12] | - | 89.6 | 83.8 | 97.1 | 98.8
| [384,396]/1.03 (w/o) | [12,16] | - | 87.8 | 81.7 | 96.1 | 98.6
| [384,396]/1.03 (w/) | [12,16] | - | 90.2 | 85.2 | 97.1 | 98.8
|  | Fusing three models
---
(w/o)
- | Weighted Sum | 90.1 | 84.8 | 97.3 | 98.9
|  Fusing three models (w/) | - | Weighted Sum | 91.0(0.8$\uparrow$) | 86.3(1.3$\uparrow$) | 97.4 | 99.0
|  | Fusing three models
---
(w/o)
- |  | Weighted
---
Concatenate
90.5 | 85.4 | 97.4 | 99.0
VehicleID |  | Fusing three models
---
(w/)
- |  | Weighted
---
Concatenate
90.7 | 85.3 | 97.7 | 99.1

TABLE I: The basic results of our method. The blue color marks the baseline
result we used for comparison and the red color marks the best performance.

Implementation Details. After analyzing the data size and aspect ratio
distributions from the VeRi-776 and VehicleID training datasets, we
standardized the resized image heights to 224 pixels for VeRi-776 and 384
pixels for VehicleID. We used a uniform patch size of $16\times 16$ for both
training and testing phases, with stride sizes set to 16 for the longer
dimension and 12 for the shorter. The aspect ratios were $[1.0,0.95,1.33]$ for
VeRi-776 and $[1.0,0.80,1.03]$ for VehicleID. In the testing stage, we used
the entire VeRi-776 test dataset and the largest VehicleID subset, containing
800 vehicles. Image preprocessing included 50% random horizontal flipping,
padding, cropping, and erasing.

All experiments were conducted on four NVIDIA RTX A6000 GPUs using PyTorch
with FP16 training. We used an SGD optimizer with a momentum of 0.1 and a
weight decay of 1e-4, maintaining equal weights of 1.0 for ID and triplet
losses. The batch size was set at 128, with four images per ID, across 120
training epochs. Initial learning rates were set at 0.035 for VeRi-776 and
0.045 for VehicleID, both decreasing linearly.

###  IV-B Results

Major Results of Our Method. Our experiments (Table I) compare baseline
performance with different input shapes and the absence of patch mixup (PM)
data augmentation. On VeRi-776, non-square input of $224\times 298$ without PM
outperforms square input of $224\times 224$ by 4.5% in mAP. With sum-based
weighted feature fusion from three models, the mAP improves by over 6.5%.
Similarly, on VehicleID, non-square input of $384\times 396$ without PM
achieves 0.6% higher mAP than square input of $384\times 384$, while feature
fusion raises mAP by over 0.8%.

Comparsion with the State-of-the-Art Methods. In Tables II and III, we compare
our best model’s ReID performance with three state-of-the-art methods (RPTM
[39], CAL [40], and TransReID [12]) on VeRi-776 and VehicleID large test
datasets. Our model, using ViT-B/16 backbone, outperforms pure ViT-B/16 by 23%
mAP and 0.5% R1 on VeRi-776. On VehicleID, our approach achieves 91.0% mAP and
97.4% R5, surpassing RPTM by 10.5% mAP and 1.1% R5. These results demonstrate
the significant potential for improving ReID performance by solely modifying
ViTs’ input without altering their architectures.

Figure 5: ReID Performance of Square Input on VeRi-776.

Figure 6: ReID Performance of Square Input on VehicleID.

Backbone |  | Input
---
Size
Method | mAP | R1 | R5
ResNet-101 | 240x240 | RPTM[39] | 88.0 | 97.3 | 98.4
ResNet-50 | 384x192 | CAL[40] | 74.3 | 95.4 | 97.9
ViT/B-16 | 384x128 | TransReID[12] | 82.0 | 97.1 | -
ViT/B-16 | 256x128 |  | ViT/B-16
---
Baseline[12]
78.2 | 96.5 | -
ViT/B-16 | Fused | Ours | 81.5 | 97.0 | 98.5

TABLE II: Comparison with state-of-the-art methods of Re-ID on VeRi-776

Backbone |  | Input
---
Size
Method |  | mAP
---
(%)
R1(%) | R5(%)
ResNet-101 | 240x240 | RPTM[39] | 80.5 | 92.9 | 96.3
ResNet-50 | 384x192 | [40] | 80.9 | 75.1 | 88.5
ViT/B-16 | 256x128 | TransReID[12] | - | 85.2 | 97.5
ViT/B-16 | 256x128 |  | ViT/B-16
---
Baseline[12]
- | 83.5 | 96.7
ViT/B-16 | Fused | Ours | 91.0 | 86.3 | 97.4

TABLE III: Comparison with state-of-the-art methods of Re-ID on VehicleID
large dataset

Ablation Study: Square input with different sizes. We assess the impact of
input size on ReID performance using square shapes of $224\times 224$,
$256\times 256$, $384\times 384$, and $416\times 416$ for both VeRI-776 and
VehicleID datasets. All the trainings are conducted without using the PM
module. The results are shown in Fig. 6 and 6. On both datasets, the ReID
performance does not linearly increase with the input size. The highest ReID
accuracy is observed when the input size is 256 pixels for VeRi-776 and 384
pixels for VehicleID.

Figure 7: Average inference Time of single Image from Test Dataset.  Figure 8:
Attention map on VeRi-776 without and with PM module. The first row shows the
results without using PM module, and the second row shows the results using PM
module.  Figure 9: Attention map on VehicleID without and with PM module. The
first row shows the results without using PM module, and the second row shows
the results using PM module.

###  IV-C Results Analysis

Patch Number affects Average Inference Time. Fig. 7 illustrates the
correlation between patch count and inference time per object. Increased
patches escalate time costs, indicating optimization opportunities. Despite
longer inference times from model fusion, mitigation is possible through mask
or attention strategies to trim patch count during training, thus curtailing
parameter size.

How does patch mixup (PM) work? Grad-CAM visualizes attention (Fig. 8 and 9).
Intra-image Patch Mixup improves object focus, especially in VeRi-776,
enhancing robustness to view and aspect ratios, notably in VehicleID.
Utilizing self-attention, our method enhances spatial relations, improving
global feature learning in ViTs.

##  V Conclusion

We explore aspect ratio’s impact on ViT-based vehicle ReID. Fusion of varied
aspect ratio models boosts robustness and Re-ID performance. Our intra-image
Patch Mixup augmentation enhances generalization. Outperforming baselines on
VeRi-776 and VehicleID, though fusion may raise inference time, network
pruning can be used to mitigate this. Future work aims at efficient ReID
models for diverse aspect ratios.

## References

  * [1] Hongbo Wang  “A survey of vehicle re-identification based on deep learning”  In _IEEE Access_ 7 IEEE, 2019, pp. 172443–172469
  * [2] Zhedong Zheng  “VehicleNet: Learning robust visual representation for vehicle re-identification”  In _IEEE Transactions on Multimedia_ 23 IEEE, 2020, pp. 2683–2693
  * [3] Xiaobo Chen  “Global-Local Discriminative Representation Learning Network for Viewpoint-Aware Vehicle Re-identification in Intelligent Transportation”  In _IEEE Transactions on Instrumentation and Measurement_ IEEE, 2023
  * [4] Jinjia Peng  “Learning multi-region features for vehicle re-identification with context-based ranking method”  In _Neurocomputing_ 359 Elsevier, 2019, pp. 427–437
  * [5] Xinyu Zhang  “Part-guided attention learning for vehicle instance retrieval”  In _TITS_ 23.4 IEEE, 2020, pp. 3048–3060
  * [6] Jianyang Gu  “An efficient global representation constrained by Angular Triplet loss for vehicle re-identification”  In _Pattern Analysis and Applications_ 24 Springer, 2021, pp. 367–379
  * [7] Ke Yan  “Exploiting multi-grain ranking constraints for precisely searching visually-similar vehicles”  In _ICCV_ , 2017, pp. 562–570
  * [8] Xinchen Liu  “A deep learning-based approach to progressive vehicle re-identification for urban surveillance”  In _ECCV_ , 2016, pp. 869–884  Springer
  * [9] Haiyun Guo  “Learning coarse-to-fine structured feature embedding for vehicle re-identification”  In _AAAI_ 32.1, 2018
  * [10] Yihang Lou  “Veri-wild: A large dataset and a new method for vehicle re-identification in the wild”  In _ICCV_ , 2019, pp. 3235–3243
  * [11] Zakria  “Trends in vehicle re-identification past, present, and future: A comprehensive review”  In _Mathematics_ 9.24 MDPI, 2021, pp. 3162
  * [12] Shuting He  “Transreid: Transformer-based object re-identification”  In _ICCV_ , 2021, pp. 15013–15022
  * [13] Junjie Ke  “Musiq: Multi-scale image quality transformer”  In _ICCV_ , 2021, pp. 5148–5157
  * [14] Wanli Liu  “Is the aspect ratio of cells important in deep learning? A robust comparison of deep learning methods for multi-scale cytopathology cell image classification: From convolutional neural networks to visual transformers”  In _Computers in biology and medicine_ 141 Elsevier, 2022, pp. 105026
  * [15] Xiaofeng Mao  “Towards robust vision transformer”  In _CVPR_ , pp. 12042–12051
  * [16] Mostafa Dehghani  “Patch n’Pack: NaViT, a Vision Transformer for any Aspect Ratio and Resolution”  In _arXiv preprint arXiv:2307.06304_ , 2023
  * [17] Zhuofan Xia  “Vision transformer with deformable attention”  In _CVPR_ , 2022, pp. 4794–4803
  * [18] Pengyuan Lv  “SCViT: A spatial-channel feature preserving vision transformer for remote sensing image scene classification”  In _IEEE Transactions on Geoscience and Remote Sensing_ 60 IEEE, 2022, pp. 1–12
  * [19] Hyeonbin Hwang  “Vision transformer equipped with neural resizer on facial expression recognition task”  In _ICASSP_ , 2022, pp. 2614–2618  IEEE
  * [20] Hancheng Zhu  “ARET-IQA: An Aspect-Ratio-Embedded Transformer for Image Quality Assessment”  In _Electronics_ 11.14 MDPI, 2022, pp. 2132
  * [21] Siyuan Li  “Multi-aspect convolutional-transformer network for SAR automatic target recognition”  In _Remote Sensing_ 14.16 MDPI, 2022, pp. 3924
  * [22] H Zhang  “mixup: Beyond empirical risk management”  In _ICLR_ , 2018, pp. 1–13
  * [23] Ashish Vaswani  “Attention is all you need”  In _Advances in neural information processing systems_ 30, 2017
  * [24] Alexey Dosovitskiy  “An image is worth 16x16 words: Transformers for image recognition at scale. arXiv 2020”  In _arXiv preprint arXiv:2010.11929_ , 2010
  * [25] Hugo Touvron  “Training data-efficient image transformers & distillation through attention”  In _International conference on machine learning_ , 2021, pp. 10347–10357  PMLR
  * [26] Ze Liu  “Swin transformer: Hierarchical vision transformer using shifted windows”  In _ICCV_ , 2021, pp. 10012–10022
  * [27] Wenhai Wang  “Pyramid vision transformer: A versatile backbone for dense prediction without convolutions”  In _ICCV_ , 2021, pp. 568–578
  * [28] Xiangxiang Chu  “Conditional positional encodings for vision transformers”  In _arXiv preprint arXiv:2102.10882_ , 2021
  * [29] Salman Khan  “Transformers in vision: A survey”  In _ACM computing surveys (CSUR)_ 54.10s ACM New York, NY, 2022, pp. 1–41
  * [30] Kaiming He  “Deep residual learning for image recognition”  In _CVPR_ , 2016, pp. 770–778
  * [31] Raja Muhammad Saad Bashir  “Vr-proud: Vehicle re-identification using progressive unsupervised deep architecture”  In _Pattern Recognition_ 90 Elsevier, 2019, pp. 52–65
  * [32] Geoffrey Roman-Jimenez  “Improving vehicle re-identification using CNN latent spaces: Metrics comparison and track-to-track extension”  In _IET Computer Vision_ 15.2 Wiley Online Library, 2021, pp. 85–98
  * [33] Jiawei Lian  “Transformer-based attention network for vehicle re-identification”  In _Electronics_ 11.7 MDPI, 2022, pp. 1016
  * [34] Hao Luo  “An empirical study of vehicle re-identification on the AI City Challenge”  In _CVPR_ , 2021, pp. 4095–4102
  * [35] Ran Wei  “Transformer-Based Domain-Specific Representation for Unsupervised Domain Adaptive Vehicle Re-Identification”  In _TITS_ 24.3 IEEE, 2022, pp. 2935–2946
  * [36] Zhedong Zheng  “A discriminatively learned cnn embedding for person reidentification”  In _TOMM_ 14.1 ACM New York, NY, USA, 2017, pp. 1–20
  * [37] Hao Liu  “End-to-end comparative attention networks for person re-identification”  In _TIP_ 26.7 IEEE, 2017, pp. 3492–3506
  * [38] Raia Hadsell  “Dimensionality reduction by learning an invariant mapping”  In _CVPR_ 2, 2006, pp. 1735–1742  IEEE
  * [39] Adhiraj Ghosh  “Relation preserving triplet mining for stabilising the triplet loss in re-identification systems”  In _WACV_ , 2023, pp. 4840–4849
  * [40] Yongming Rao  “Counterfactual attention learning for fine-grained visual categorization and re-identification”  In _ICCV_ , 2021, pp. 1025–1034