# Honeycomb anti-dot artificial lattice
as a prototypical correlated Dirac fermions system

A. Biborski<EMAIL_ADDRESS>Academic Centre for Materials and
Nanotechnology, AGH University of Krakow, Al. Mickiewicza 30, 30-059 Krakow,
Poland  M. Zegrodnik<EMAIL_ADDRESS>Academic Centre for
Materials and Nanotechnology, AGH University of Krakow, Al. Mickiewicza 30,
30-059 Krakow, Poland  P. Wójcik<EMAIL_ADDRESS>Faculty of
Physics and Applied Computer Science, AGH University of Krakow, Al.
Mickiewicza 30, 30-059 Krakow, Poland  M. P. Nowak<EMAIL_ADDRESS>Academic Centre for Materials and Nanotechnology, AGH University of Krakow,
Al. Mickiewicza 30, 30-059 Krakow, Poland

###### Abstract

We study theoretically the electronic properties of the artificial quantum dot
honeycomb lattice defined in a two-dimensional electron gas, focusing on the
possibility of achieving a regime in which electronic correlations play a
dominant role. At first we establish a non-interacting model compatible with
recently studied experimentally devices. According to the values of the
obtained electron-electron interaction integrals, we postulate that the
inclusion of inherent electron-gas self-screening is indispensable to
reconstruct the experimental observations. Applying the Thomas-Fermi type of
screening, we show that the radius of the anti-dot is crucial to achieve a
correlated state in which phenomena like antiferromagentic ordering and
interaction-induced insulating state appear. We estimate the conditions for
which the electronically correlated state in an artificial honeycomb lattice
can be realized.

##  I  Introduction

There is a long-standing debate about the role of electron-electron
interactions in two-dimensional Dirac fermion systems Kotov  _et al._ (2012);
Tang  _et al._ (2018); Das Sarma  _et al._ (2007). In particular, for
graphene, massless fermionic quasiparticles are believed to be strongly
renormalized by the presence of electronic interactions Kotov  _et al._
(2012). However, this renormalization influences the Fermi velocity and thus
keeps the semi-metallic character of the system instead of inducing the Mott
insulator phase. The unscreened long-range interactions in two-dimensional
systems may also play a competitive role with respect to Mott localization
Tang  _et al._ (2018); Wu and Tremblay (2014); Assaad and Herbut (2013);
Sorella  _et al._ (2012). However, short-range Hubbard interactions $U$ are
believed to be responsible for the formation of a gap in strongly correlated
systems when the kinetic energy scale given by the hopping amplitude $|t|$ is
noticeably smaller than $U$. In the absence of long-range interactions, the
gap opening in the Dirac fermion systems is supposed to arise for
$U/|t|\gtrsim 3.8\div 4.2$ Meng  _et al._ (2010); Sorella  _et al._ (2012);
Raczkowski  _et al._ (2020). This quantum phase transition from the semi-metal
to the antiferromagnetic insulator at half-filling has been pointed out in a
number of theoretical studiesTang  _et al._ (2018); Wu and Tremblay (2014);
Otsuka  _et al._ (2016); Liebsch and Wu (2013); Assaad and Herbut (2013); Meng
_et al._ (2010); Sorella  _et al._ (2012); Arya  _et al._ (2015); Raczkowski
_et al._ (2020).

The electronic properties of honeycomb systems can be conveniently studied
using the so-called _quantum simulators_ such as artificial lattices of
quantum dots (ALDs) defined in a two-dimensional electron gas (2DEG). This
approach provides the possibility of tuning the properties of the system by
controlling the size of the dots and the spacing between them, as well as the
depth of their confinement potential Saleem  _et al._ (2022). The emergence of
Dirac systems in ALD has been experimentally demonstrated by means of
photoluminescence measurements in a nanodevice constituted in the so-called
anti-dot (AD) architecture Du  _et al._ (2021). The AD approach provides the
opportunity to form quasi-2D bulk systems, contrary to quantum dot assemblies
fabricated by electrostatic gating, where the number of trapping centers
typically does not exceed a dozen. Nevertheless, the photoluminescence spectra
presented by Du et al. Du  _et al._ (2021) indicate the presence of the
spectral dublet related to excitations associated with van Hove singularities,
which maps accurately on the model of non-interacting electrons. Thus, at
least for the AD radius and lattice spacing that was examined, this device
cannot be regarded as a strongly correlated system. Therefore, the natural
question arises Saleem  _et al._ (2022), if tuning of these two parameters may
result in forcing ALD to be in the correlated state, characterized, e.g., by
formation of the Mott gap, emergence of antiferromagnetic (AF) order, etc. The
artificial honeycomb lattices formed in 2DEG have been theoretically studied
in view of their electronic properties by a priori assuming an idealized
confinement similar to graphene Kylänpää  _et al._ (2016); Li  _et al._
(2020); Räsänen  _et al._ (2012) or ADs Saleem  _et al._ (2022); Krix  _et
al._ (2022); Li  _et al._ (2017). These valuable pioneering papers provide
evidence of the influence of electron-electron interactions on the ALD band
structure in the framework of Density Functional Theory Räsänen  _et al._
(2012), possible mechanism of pairing Li  _et al._ (2020), or the stability of
the Dirac cone with respect to the shape of the confining potential Kylänpää
_et al._ (2016) or disorder Tkachenko  _et al._ (2015).

In this work, we investigate conditions for realization of a correlated state
in an artificial Dirac lattice in realistic structures. For concreteness, we
take the device as recently realized by Du et al. Du  _et al._ (2021) as a
starting point and investigate the electronic properties of 2DEG ALD,
establishing the conditions under which those systems transit into the
correlated phase. We present the AD-ALD model elaborated within multiscale
simulations. Namely, first, we construct the mean-field model exploiting the
Schrödinger-Poisson scheme. Subsequently, we utilize the resulting single-
particle picture for construction of the Wannier basis by means of the
projection method and subsequently calculate the electron-electron interaction
amplitudes. We discuss the role of electronic screening, which we include by
applying the Thomas-Fermi model with the screening length characteristic for
GaAs-based 2DEG. Inclusion of screening is indispensable for the
reconstruction of experimental finding, as bare (non-screened) amplitudes
would imply a strongly correlated system for the lattice spacing and AD radius
for which a semimetallic character has been reported Du  _et al._ (2021).
Eventually, we elaborate the interacting Hamiltonian in the second
quantization formulation and solve it for the half-filled case by means of the
Variational Monte Carlo (VMC) method. These simulations provide evidence of
the transition from the semimetallic state to the AF phase and indicators of
the emergence of a Mott gap with an increase in AD diameter.

##  II Model and Mean Field approach

###  II.1 Device and AD periodic potential

We consider an artificial honeycomb lattice defined in a GaAs/AlGaAs
heterostructure Du  _et al._ (2021) (Fig. 1(a)). ADs are assembled in a
triangular lattice (Fig.1(b)) spanned by vectors $\mathbf{R}_{1}=(L,0)$ and
$\mathbf{R}_{2}=(L/2,L\sqrt{3}/2)$ where $L=|\mathbf{R}_{1}|=70$ nm. The
points related to
$\mathbf{R}_{ij}=i\times\mathbf{R}_{1}+j\times\mathbf{R}_{2}$ are associated
with maxima of potential $V_{r_{0}}(\mathbf{r})$ resulting from the patterned
etching of the top layers of GaAs/AlGaAs. Here, we model the periodic
potential assuming a Gaussian contraction, which is given as

$\displaystyle
V_{r_{0}}(\mathbf{r})=\sum_{i,j}V_{0}\exp{\Bigg{[}-\Bigg{(}\frac{||\mathbf{R}_{ij}-\mathbf{r}||}{r_{0}}\Bigg{)}^{2}\Bigg{]}},$
(1)

where $\mathbf{r}=(x,y)$, $V_{0}$, and $r_{0}$ are related to the maximum
height of the potential and the radius of the antidot, respectively. The
Gaussian functional form is believed to properly describe the trapping
potential that comprises quantum dots of diameter less than $\sim 100$ nm in
2DEG formed in GaAs/AlGaAs heterostructures Bednarek  _et al._ (2003); Ciurla
_et al._ (2002). As the estimated radius of the etched antidots in the
experiment Du  _et al._ (2021) is $r_{0}\approx 20\pm 5$ nm, we inspect the
values of $r_{0}$ in a similar range. Note that since $r_{0}\approx L/3$, the
local minima in the landscape of $V_{r_{0}}(\mathbf{r})$ present in Fig. 1(c)
correspond to dots that form a graphene-like lattice.

Figure 1: (a) Cross section of the considered heterostructure. (b) Two
dimensional plot of electrostatic potential induced by etching in 2DEG. Filled
blue circles indicate examplary locations of potential minima.(c) The
electrostatic potential profile in 2DEG along the direction marked as dashed
line in (b). Arrows refer to the potential minima.

The potential amplitude $V_{0}$ in Eq. (1) is determined by means of the
standard Schrödinger-Poisson approach used separately for two cross sections
depicted in Fig. 1(a) by dashed lines. Assuming translation symmetry in the
$x-y$ plane, $\Psi(\mathbf{r})=\phi_{n}(z)\exp(ik_{x}x+ik_{y}y)$, the
Schrödinger equation can be reduced to the 1D form

$\displaystyle\left(-\frac{1}{2}\frac{d}{dx}\frac{1}{m^{*}_{\perp}(z)}\frac{d}{dx}+V_{r_{0}}(z)+\frac{\hbar^{2}k^{2}}{2m^{*}_{\parallel}}\right)\phi_{n}(z)=E_{n}\phi_{n}(z),$
(2)

where $k^{2}=k_{x}^{2}+k_{y}^{2}$, $m^{*}_{\perp(\parallel)}$ is the effective
mass in the direction perpendicular (parallel) to the layers and the potential
$V_{r_{0}}(z)$ is the sum of two components: (i) $V_{band}(z)$ related to the
discontinuity of the conduction band at the GaAs/AlGaAs interfaces and (ii)
the electrostatic potential $V_{e}(z)$ that includes the electron-electron
interaction and the electric field from the gate. The latter can be determined
at the mean-field level from the Poisson equation given by

$\displaystyle\epsilon_{0}\frac{d}{dz}\epsilon_{r}(z)\frac{d}{dz}V_{e}(z)=-\rho(z),$
(3)

with the charge density $\rho(z)=\rho_{d}(z)+\rho_{e}(z)$, where $\rho_{d}(z)$
is the doping profile and $\rho_{e}(z)$ is the electron distribution. The
electron density is obtained by

$\displaystyle\rho_{e}(z)=2\sum_{n}|\phi_{n}(z)|^{2}f_{2D}(E_{n}-\mu),$ (4)

where the factor $2$ accounts for the spin degeneracy, $T$ is the temperature,
$\mu$ is the chemical potential and $f_{2D}$ is the Fermi-Dirac distribution
integrated over the $k_{x}$ and $k_{y}$ components of the wave vector,

$\displaystyle
f_{2D}(E_{n},\mu)=\frac{m^{*}_{\parallel}k_{B}T}{2\pi}\ln\left(1+\exp\left(-\frac{E_{n}-\mu}{k_{B}T}\right)\right).$
(5)

Equations (2) and (3) are solved numerically using the finite difference
method with Dirichlet boundary conditions. The existence of 2DEG in the main
quantum well results from the Si delta doping located $30$ nm above. The value
of $n_{d}$ is chosen to correspond to the electron concentration in 2DEG at
the level of $n_{el}\approx 4.5\times 10^{10}\text{ cm}^{-2}$, i.e., about two
electrons per unit cell, as provided in the experiment Du  _et al._ (2021).
Then, for a chosen $n_{d}$ the self-consistent energy profile $V_{2}(z)$ is
determined for the cross section with the etched region (see Fig. 1(a)) and
the potential amplitude $V_{0}$ is estimated based on the adiabatic
approximation according to the formula

$\displaystyle V_{0}\approx E_{2}^{1}-E_{1}^{1},$ (6)

where $E_{2}^{1}$ and $E_{1}^{2}$ are energies related to the lowest lying
states for the etched and non-etched cases, respectively. The procedure is
repeated until self-consistency is reached, which we consider to occur when
the potential variation between two consecutive iterations is less than
$10^{-7}$ eV. Calculations have been carried out for the material parameters
corresponding to GaAs (AlGaAs): $m^{*}_{\perp}($GaAs$)=0.067$,
$m^{*}_{\perp}($AlxGa1-xAs$)=0.067+0.083x$,
$m^{*}_{\parallel}($GaAs$)=m^{*}_{\parallel}($AlxGa1-xAs$)=0.067$,
$\epsilon_{r}($GaAs$)=13.18$, $\epsilon_{r}($AlxGa1-xAs$)=13.18+3.12x$ and the
conduction band minima $E_{c}=E_{g}/2$, where the energy gap
$E_{g}($GaAs$)=1.43$ eV, $E_{g}($AlxGa1-xAs$)=1.43+1.247x$. We set $x=0.3$ and
$T=4$ K.

Figure 2: The potential profiles $V_{1}(z)$ and $V_{2}(z)$. Note, that
$V_{2}(z)$ (etched case) is evaluated for Si doping obtained for the unetched
case. The solid black line indicates the carrier density $\rho_{e}(z)$
computed for potential $V_{1}(z)$, which corresponds to $n_{el}=\int_{z\in
2DEG}\rho_{e}(z)dz\approx 4.5\times 10^{10}$ cm-2. 111Note that we excluded
the bottom layer of $n$-doped AlGaAs as detailed information regarding the
character of doping has not been provided Du  _et al._ (2021)

The resulting potential profiles $V_{1}(z)$ (non-etched case) and $V_{2}(z)$
(depth of etching $35$ nm) are presented in Fig. 2. The local maximum of the
AD potential is estimated to be $V_{0}\approx 600$ meV according to Eq. 6. We
confirm the reliability of our procedure by comparing the calculated Fermi
energy $\mu$ as a function of $n_{el}$. In Fig. 3 we present $\mu(n_{el})$
obtained for the potential profile $V_{1}(z)$, as well as for the available
experimental data Du  _et al._ (2021) and the dependence corresponding to the
ideal 2DEG model. The Fermi energy for the considered concentration range is
$2$ meV $\lesssim\mu\lesssim 2.5$ meV. Remarkably, the energies evaluated
within Schrödinger-Poisson scheme are in almost perfect agreement with those
obtained for the ideal 2DEG with values only about $5\%$ lower than the
experimental ones.

Figure 3: The Fermi energy as a function of the electron concentration for the
unpatterned 2DEG calculated within Schrödinger-Poisson scheme and for the
ideal 2DEG (red circles and black solid line, respectively) compared to the
experimental data (blue diamonds).

The estimated value of $V_{0}$ is finally applied to Eq. 1, which defines the
planar 2D potential for a given parameter $r_{0}$. We focus on AD radii $15$
nm $\leq r_{0}\leq 21$ nm. These values provide $V_{r_{0}}$ that exposes
_local_ maxima at
$\mathbf{R}_{lmax}=p/2\times\mathbf{R}_{1}+q/2\times\mathbf{R}_{2}$ and minima
at $\mathbf{R}_{lmin}=u\times\mathbf{R}_{1}+v\times\mathbf{R}_{2}$ where
$p,q\in\mathbb{Z}$ and $u,v\in\\{p/3,2q/3\\}$ (Fig. 4(a)). In particular, the
points defined by the vectors $\mathbf{R}_{lmin}^{uv}$ correspond to the
honeycomb lattice. The depth of the local minima
$V_{AG}=V_{r_{0}}(\mathbf{R}_{lmax})-V_{r_{0}}(\mathbf{R}_{lmin})$ depends on
$r_{0}$. For the lower range of $r_{0}$ considered, $V_{AG}$ is only of the
order of a few microelectronvolts, while for $r_{0}=20$ nm $V_{AG}\approx 26$
meV, as derived from the data shown in Fig. 4(b). Thus, it may be expected
that a sufficiently large $r_{0}$ enforces a stronger confinement of electrons
and, in turn, their localization in honeycomb trapping centers, which gives
the opportunity for emerging strongly correlated state.

Figure 4: (a) Cross section of the potential $V_{r_{0}}$ for different values
of $r_{0}$. The increase in $r_{0}$ results in deeper local minima
corresponding to the formation of a graphene-like lattice. (b) The depth of
the local minima $V_{AG}$ as a function of $r_{0}$.

###  II.2 Artificial lattice: Mean-field approach

The elaboration of $V_{r_{0}}(z)$ eventually allows us to perform an analysis
of the electronic properties of the formed artificial honeycomb lattice.
First, we provide a mean-field approach in which we decouple the modulated
2DEG from the rest of the system. It is done by factorizing a single-electron
wave function $\Psi_{n}(\mathbf{r},z)=\psi_{n}(\mathbf{r})\phi_{1}(z)$. Within
this assumption, to estimate whether the distribution of electronic density
(at a concentration of two electrons per elementary cell) for the estimated
potential is honeycomb-like, we solve again the Schrödinger-Poisson
problem,this time, in $x-y$ plane for $\psi_{n}(\mathbf{r})$, i.e.,

$\displaystyle\Big{[}-\frac{\hbar^{2}}{2m^{*}_{||}}\nabla^{2}_{r}+V_{r_{0}}(\mathbf{r})+V_{c}(\mathbf{r})\Big{]}\psi_{n}(\mathbf{r})=E_{n}\psi(\mathbf{r})$
(7a)
$\displaystyle\nabla^{1}_{r}V_{c}(\mathbf{r})=-\frac{\rho_{r}(\mathbf{r})}{\epsilon_{0}\epsilon_{r}}$
(7b)
$\displaystyle\rho_{r}(\mathbf{r})=\frac{1}{d}\sum_{n,\sigma}^{N_{el}/2}|\psi_{n\sigma}(\mathbf{r})|^{2},$
(7c)

where $N_{el}=2$ is the number of electrons in the unit cell and
$\sigma=\\{\uparrow,\downarrow\\}$ stands for the $S^{z}$ spin component. The
factor $1/d$ in Eq. 7c refers to _mean_ electron concentration (per spin) with
respect to $z$ direction in 2DEG. Eventually the estimated _total_ mean field
2D potential at given $r_{0}$ is modeled as a sum of $V_{r_{0}}(\mathbf{r})$
and $V_{c}(\mathbf{r})$ resulting from the self-consistent numerical solution
of Eqs. 7(a-c), that is,

$\displaystyle\tilde{V}(\mathbf{r})\equiv
V_{r_{0}}(\mathbf{r})+V_{c}(\mathbf{r}).$ (8)

Eventually, $V_{r_{0}}$ is used to solve the two-dimensional mean field model
defined in Eqs. (7a-c). In our calculations we take a single unit cell and
assume periodic boundary conditions. A computational procedure is performed on
the triangular mesh consisting of $70\times 70$ nodes. We consider the half-
filled case, that is, $\int_{\Omega}\rho_{r}(\mathbf{r})d\mathbf{r}=2$, where
$\Omega$ is an area of the unit cell.

This self-consistent mean-field (MF) approach reveals a meaningful
renormalization of the density distribution compared to that of the free
electron (non-interacting) case. Namely, the inclusion of Coulomb interactions
leads to the more uniform _smearing_ of carriers throughout $\Omega$ as can be
deduced from Fig. 5. Particularly, this tendency is meaningful for the lower
values of $r_{0}$, as $V_{AG}$ increases with radius of AD.

Figure 5: Cross sections of the resulting electronic density along
$\mathbf{R}_{1}+\mathbf{R}_{2}$ vector for the three representative values of
$r_{0}$.

##  III Wannier basis and Single-particle Hamiltonian

The task of constructing the Wannier basis is crucial in view of the reliable
determination of both the single-particle Hamiltonian in Tight Binding
Approximation (TBA) and the electron-electron interaction amplitudes. In the
considered system, the minima of the trapping centers are expected to be
shallow, specifically for the semimetallic case, and the distance between the
neighboring dots in the honeycomb lattice is $\sim r_{0}$. Here, the Wannier
functions are obtained by application of the _Projection Method_ Marzari  _et
al._ (2012) adapted to a 2D case. This method is based on projecting smooth
trial functions $g_{n}(\mathbf{r})$ centered at positions that are supposed to
correspond to the maxima of the electronic density onto Bloch states
represented by $\psi_{n,\mathbf{k}}(\mathbf{r})$ (where $n$ refers to the band
index). Subsequently, the procedure that provides unitary transformation is
carried out, which eventually results in well-localized Wannier states when
$\\{g_{n}(\mathbf{r})\\}$ are correctly guessed. Namely, we obtain
$\psi_{n,\mathbf{k}}(\mathbf{r})$ by numerical diagonalization of the single-
particle mean-field Hamiltonian $\mathcal{\hat{H}}_{0}$ defined as

$\displaystyle\mathcal{\hat{H}}_{0}=-\frac{\hbar^{2}\mathbf{k}^{2}}{2m^{*}_{||}}+\tilde{V}(\mathbf{r}).$
(9)

Next, we take the auxiliary basis $\\{\ket{\phi_{n\mathbf{k}}}\\}$ defined as

$\displaystyle\ket{\phi_{n\mathbf{k}}}=\sum_{m}\ket{\psi_{m\mathbf{k}}}\braket{\psi_{m\mathbf{k}}}{g_{n}},$
(10)

where the summation index $m$ enumerates the considered bands, to construct
the Löwdin-orthonormalized Bloch-like states Marzari  _et al._ (2012)

$\displaystyle\ket{\tilde{\psi}_{n\mathbf{k}}}=\sum_{m}\ket{\phi_{m\mathbf{k}}}\big{(}S_{\mathbf{k}}^{-1/2}\big{)}_{mn},$
(11)

where
$\big{(}S_{\mathbf{k}}\big{)}_{mn}=\braket{\phi_{m\mathbf{k}}}{\phi_{n\mathbf{k}}}$.
Eventually, Wannier states $\\{w_{n\mathbf{R}_{ij}}\\}$ are obtained by
applying the standard Fourier transform, i.e.,

$\displaystyle\ket{w_{n\mathbf{R}_{ij}}}=\frac{1}{2}\frac{L^{2}\sqrt{3}}{(2\pi)^{2}}\int_{BZ}d\mathbf{k}\exp\big{(}{-i\mathbf{k}\cdot\mathbf{R_{ij}}}\big{)}\ket{\tilde{\psi}_{n\mathbf{k}}}.$
(12)

For each considered $r_{0}$, the numerical representation of the Hamiltonian
given in Eq. 9 is implemented in KWANT package Groth  _et al._ (2014) and
diagonalized in the momentum space for the particular wave vector
$\mathbf{k}$. We probe $32\times 32$ equally spaced points from the first
Brillouin in the unit cell spanned by vectors
$\mathbf{G}_{1}=\frac{2\pi}{L}\Big{(}1,-\frac{\sqrt{3}}{3}\Big{)}$ and
$\mathbf{G}_{2}=\frac{2\pi}{L}\Big{(}0,\frac{2\sqrt{3}}{3}\Big{)}$ that
conform to the standard relation
$\mathbf{G}_{i}\cdot\mathbf{R}_{j}=2\pi\delta_{ij}$. We choose Gaussian trial
functions $g_{n}(\mathbf{r})$,

$\displaystyle
g_{n}(\mathbf{r})\equiv\frac{1}{2\pi\sigma^{2}}\exp{\Bigg{[}-\Bigg{(}\frac{|\mathbf{r}-\mathbf{R}_{lmin}^{n}|}{\sigma}\Bigg{)}^{2}\Bigg{]}}$
(13)

with $n\in\\{\alpha,\beta\\}$, where $\alpha$ and $\beta$ correspond to two
minima present in the confining potential surface at
$\mathbf{R}^{\alpha}_{lmin}=L(1/2,\sqrt{3}/6)$ and
$\mathbf{R}^{\beta}_{lmin}=L(1,\sqrt{3}/3)$ respectively. Thus, we construct
the Wannier basis properly describing the two lowest bands of the system. The
presented results correspond to $\sigma=40$ nm, although our test calculations
have not revealed any significant differences for $\sigma$, which are $20$ and
$60$ nm. The integration manifold in Eq. 10 has been chosen to be the area
corresponding to $5\times 5$ unit cells to ensure the proper decay of
$g_{n}(\mathbf{r})$ for the assumed value of $\sigma$.

Figure 6: 2D plot of Wannier functions $w_{\alpha R_{00}}(x,y)$ for $r_{0}=16$
nm (a), $r_{0}=18$ nm (b) and $r_{0}=20$ nm (c). In (d) their cross sections
along $x=35$ nm (dashed vertical lines in (a-c)) are presented as a function
of $y$. The disappearance of local minima with an increasing value of $r_{0}$
can be identified in terms of isoline inspection (a-c), which is accompanied
by an increase in maximal amplitude at the localization center (a-d).

In Fig. 6(a-c) we show representative plots of the Wannier functions for
$n=\alpha$ and different values of the parameter $r_{0}$. Note that they
preserve the three-fold rotational symmetry $\mathbf{C_{3}}$ as desired. The
same feature holds for $w_{\beta R_{ij}}(\mathbf{r})$ since it is reproducable
from $w_{\alpha R_{ij}}(\mathbf{r})$ by means of the following transformation,

$\displaystyle w_{\beta
R_{lm}}(\mathbf{r})=\mathbf{\hat{T}}_{\mathbf{R}_{l-i,m-j}}\mathbf{\hat{O}}_{\mathbf{R}_{ij}+\mathbf{R}_{2}-\mathbf{R}_{1}}w_{\alpha
R_{ij}}(\mathbf{r}),$ (14)

where $\mathbf{\hat{T}}$ is the translation by the vector
$\mathbf{R}_{l-i,m-j}$ and $\mathbf{\hat{O}}$ is operator of reflection with
respect to the direction given by
${\mathbf{R}_{ij}+\mathbf{R}_{2}-\mathbf{R}_{1}}$. Therefore, we disregard the
analysis for the case $n=\beta$ as the one with $n=\alpha$ is fully
representative; however, we emphasize that Eq. 14 holds (within numerical
precision) also from the perspective of our numerical calculations. The
validity of the procedure has been confirmed by inspecting the orthogonality,
which, by construction, should be fulfilled exactly. We have found that this
requirement is met within the numerical precision.

From the data contained in Fig. 6 one may conclude that the spatial extent of
the resulting Wannier states shrinks with increasing $r_{0}$, which can be
regarded as a natural consequence of mutual enforcement of trapping. However,
since these states are obtained with electron-electron interactions at the
level of the mean-field approach, we do not find it trivial, since the
confining role of the local minima in $V_{r_{0}}$ is compatible with respect
to $V_{c}$ as comes from Fig. 5. Furthermore, in addition to enhanced
localization, the nodal lines that reside in the vicinity of the location of
the nearest neighbor (nn) centers and the three local minima enclosed by them
disappear for $r_{0}$ approaching $20$ nm.

In the next step, evaluated Wannier functions are subsequently utilized for
the construction of TBA Hamiltonian given explicitely in second quantization
fromalism as

$\displaystyle\mathcal{\hat{H}}_{TB}=\sum_{\mu\nu,\sigma}t_{\mu\nu}\hat{c}_{\mu,\sigma}^{\dagger}\hat{c}_{\nu,\sigma},$
(15)

The $\mu,\nu\in\\{(i,j),n\\}$ in Eq. 15 are defined by relation
$\mathbf{R_{\mu}}\equiv\mathbf{R}_{ij}$, where $n$ number kind of Wannier
function, $\hat{c}_{\mu,\sigma}^{\dagger}$($\hat{c}_{\mu,\sigma}$) creates
(annihilates) an electron of spin $\sigma$ associated with appropriate
$w_{\mu}(\mathbf{r})$, and,

$\displaystyle
t_{\mu\nu}\equiv\Big{\langle}w_{\mu}(\mathbf{r})\Big{|}-\frac{\hbar^{2}}{2m^{*}_{||}}\nabla^{2}_{\mathbf{r}}+\tilde{V}(\mathbf{r})\Big{|}w_{\nu}(\mathbf{r})\Big{\rangle}.$
(16)

For the sake of clarity, we introduce the following notation,

$\displaystyle t_{\mu\nu}\equiv
t_{I(\mu,\nu)}:I(\mu,\nu)=Z\in\\{0,1,2,3,..\\},$ (17)

where $I(\mu,\nu)$ maps the pair of position-orbital indices into the proper
natural number $Z$ labeling $Z$-th as the nearest neighbor corresponding to
the pair $(\mu,\nu)$.

The notation is completed by assigning

$\displaystyle t_{\mu\mu}=t_{0}\equiv\epsilon_{0},$ (18)

as is the convention commonly used.

Figure 7: The single particle amplitudes as a function of $r_{0}$ up to the
$5$th-nearest neighbor. The hoppings referring to $Z\geq 5$ play a marginal
role in view of the resulting band structure for each of the examined $r_{0}$.

In Fig. 7 we present values of single particle amplitudes up to $Z=5$ as a
function of $r_{0}$. When the size of the trapping center shrinks, the value
of $\epsilon_{0}$ increases. The opposite holds for $|t_{1}|$, which is,
however, of the order of $10^{-1}$ meV in the entire examined range of
$r_{0}$. The hopping values $t_{2}$ are positive and of an order of magnitude
lower than $|t_{1}|$. Note that we have computed the hopping integrals up to
$Z=17$ and found their rapid decay with increasing $Z$, that is, $|t_{Z\geq
5}|\lesssim 10^{-3}$ meV for each value of $r_{0}$.

Figure 8: Dispersion relations for the two lowest bands obtained for
$r_{0}=15$ nm. MF and TB refer to $\mathcal{\hat{H}}_{0}$ and
$\mathcal{\hat{H}}_{TB}$ respectively. The values of
$\epsilon_{1,2}(\mathbf{k})$ resulting from the diagonalization of
$\mathcal{\hat{H}}_{TB}$ in a momentum space are in good agreement with those
that come from the solution of the eigenproblem for $\mathcal{\hat{H}}_{0}$.
Bands represented by dotted black lines correspond to the free-electron
approach.

Finally, the obtained sets $\\{t_{i}\\}_{r_{0}}$ are utilized to formulate
$\mathcal{\hat{H}}_{TB}$ for each value of $r_{0}$. In Fig. 8 we present a
spectrum for $r_{0}=15$ nm that gives $\epsilon{1,2}(\mathbf{k})$ which is in
good agreement with experimental evidence Du  _et al._ (2021). That is, the
formation of the _Dirac cone_ in the vicinity of $\mathbf{K}$ can be clearly
identified. Notably, the difference
$\Delta\varepsilon_{M}\equiv\varepsilon_{2}(\mathbf{M})-\varepsilon_{1}(\mathbf{M})$
is of particular interest, since as $\mathbf{M}=(\pi,\pi/\sqrt{3})/L$ it can
be confronted against the energies of van Hove singularity peaks in the
photolouminescence spectra. As comes from Fig. 8,
$\Delta\varepsilon_{M}\approx 1$ meV is close to the experimental value $\sim
0.9$ meV. Namely, $\varepsilon_{1,2}$ resulting from both approaches are
nearly identical, which validates the procedure covering Wannier basis
elaboration, as well as the subsequent calculation of hopping parameters
presented in the further part of the paper.

The dispersion relation obtained by the free electron approach provides a
significantly different $\Delta\varepsilon_{M}$, which is $\sim 0.6$ meV.
Also, in this case, both bands are narrower compared to the MF approach. We do
not find perfect agreement of
$\Delta\varepsilon_{M-\Gamma}\equiv\varepsilon_{1}(\mathbf{M})-\varepsilon_{1}(\mathbf{\Gamma})\approx
0.6$ meV with respect to the value $\approx 1$ meV provided by Du et al.
However, from the data contained in their work it is difficult to extract the
precise value of $\varepsilon(\mathbf{\Gamma})$ for carrier concentrations
that result in the presence of ALD. Also, since they indicate that the
emergence of the honeycomb lattice results from the filling of the lower band,
resulting in the occupation of states for which
$\varepsilon({\mathbf{k_{F}}})$ is slightly below $\varepsilon({\mathbf{M}})$,
we consider this disagreement to be a minor issue.

Importantly, with increasing $r_{0}$ we observe a decreasing value of
$\Delta\varepsilon_{M}$, as well as a narrowing of the entire band structure,
reflecting a stronger confinement of electrons in trapping centers. Therefore,
it may be expected that at some $r_{0}$ the enhancement of carrier
localization causes electron-electron interactions to dominate over kinetic
energy and, consequently, drives the system to the strongly correlated regime.
To clarify this presumption, we further study the interaction amplitudes,
construct a full electronic Hamiltonian, and find its approximate ground state
by means of VMC method for the selected sets of interactions.

##  IV Interacting picture

###  IV.1 Electronic interactions and Hamiltonian

The determination of Wannier functions and $\mathcal{\hat{H}}_{TBA}$ allows to
evaluate the complete electronic Hamiltonian $\mathcal{\hat{H}}$,

$\displaystyle\mathcal{\hat{H}}=\mathcal{\hat{H}}_{TB}+\mathcal{\hat{H}}_{e-e},$
(19)

where $\mathcal{\hat{H}}_{e-e}$ contains electron-electron interaction terms.
Although the finite width of 2DEG affects $\mathcal{\hat{H}}_{TB}$ only by the
additive constant (related to confinement in the $z$ direction), it is not the
case for the interacting part Biborski  _et al._ (2021). Thus,
$\mathcal{\hat{H}}_{e-e}$ takes the following form

$\displaystyle\mathcal{\hat{H}}_{e-e}=\frac{1}{2}\sum_{\begin{subarray}{c}\mu,\nu,\\\
\gamma,\tau\end{subarray}}\sum_{\sigma,\sigma^{\prime}}V_{\mu\nu\gamma\tau}\hat{c}^{\dagger}_{\mu,\sigma}\hat{c}^{\dagger}_{\nu,\sigma^{\prime}}\hat{c}_{l,\tau,\sigma^{\prime}}\hat{c}_{\gamma,\sigma},$
(20a) $\displaystyle
V_{\mu\nu\gamma\tau}=\Big{\langle}\widetilde{w}_{\mu}(\mathbf{r},z)\widetilde{w}_{\nu}(\mathbf{r^{\prime}},z^{\prime})\Big{|}\hat{V}_{e-e}\Big{|}\widetilde{w}_{\gamma}(\mathbf{r},z)\widetilde{w}_{\tau}(\mathbf{r^{\prime}},z^{\prime})\Big{\rangle},$
(20b) $\displaystyle\widetilde{w}_{\mu}(\mathbf{r},z)\equiv
w_{\mu}(\mathbf{r})\phi(z).$ (20c)

Note, that we take $\phi(z)=\phi_{1}(z)$, i.e., the lowest energy state
obtained within Schrödinger-Poisson scheme as described in Sec. II A. The
electron-electron interaction operator $\hat{V}_{e-e}$ is taken in the Yukawa
potential form, to account for the screening effects among electrons in 2DEG
(see Section Results), that is,

$\displaystyle\hat{V}_{e-e}=\frac{e^{2}\exp{\Big{(}-q_{TF}\sqrt{|\mathbf{r}-\mathbf{r}^{\prime}|^{2}+(z-z^{\prime})^{2}}\Big{)}}}{4\pi\epsilon_{0}\epsilon_{r}\sqrt{|\mathbf{r}-\mathbf{r}^{\prime}|^{2}+(z-z^{\prime})^{2}}},$
(21)

where $e$ is the electron charge and $q_{TF}$ is the estimated length of the
Thomas-Fermi wave vector.

Evaluation of microscopic parameters ($t_{\mu\nu}$,$V_{\mu\nu\gamma\tau}$)
eventually leads to formulation of the Hamiltonian $\mathcal{\hat{H}}$, which
can be diagonalized only by approximate methods. Here, we exploit the
Variational Monte Carlo technique Biborski  _et al._ (2021) as will be
presented at the end of this section.

The leading elements of the electron-electron interaction tensor
$V_{\mu\nu\gamma\tau}$ provided in Eq. 20 can be obtained numerically.
However, to perform this task, the length of the Thomas-Fermi wave vector
$q_{TF}$ in Eq. 21—which tunes the magnitude of screening—needs to be
reasonably estimated.

The screening, directly related to the density of the electron gas, together
with the local confinement of the carriers, may play a key role in the tuning
of the magnitude of electron-electron interactions Knörzer  _et al._ (2022).
In our case, the Wigner-Seitz radius $r_{s}^{nD}$ (where $D$ stands for
dimensionality) at $n_{el}\approx 4.5\times 10^{10}$ cm-2 is
$r_{s}^{2D}\approx 2.66$ when treating the system as strictly two-dimensional,
and $r_{s}^{3D}\approx 2.36$ when one considers it as quasi-two-dimensional,
that is, as being a layer of width $d=25$ nm. Here, both values of $r_{s}$ are
accidentally close to each other. As they are $\sim 2.5$, one concludes that
electron-electron interactions and kinetic energy are of similar magnitude
Kotov  _et al._ (2012); Ortiz  _et al._ (1999). Thus, disregarding the
modulation described by $V_{r_{0}}$ for a while, the system may be safely
considered to be the one in which the influence of Coulomb interactions
between carriers is moderate with respect to its electronic properties, namely
$1\lesssim r_{s}\lesssim 10$. The system falls into the class of the so-called
_intermediate regime_ Ortiz  _et al._ (1999) with respect to the relation
between electron-electron interactions and their kinetic energy.

Also, as for $r_{0}=15$ nm $V_{r_{0}}$ exhibits only shallow local minima and
the resulting density profiles are characterized by the nearly homogeneous
distribution of electrons (excluding areas of ADs) with relatively weak local
maxima located at honeycomb lattice sites, it seems reasonable to assume
metallic-like screening for this case. Therefore, increasing $r_{0}$ should
also affect the screening length. However, assuming its value as that in the
non- or weakly confinig regime, one only overestimates the magnitude of
$q_{TF}$ and when interactions become significantly greater than hopping
amplitudes, one may deduce that the strongly correlated regime is achieved.
Therefore, such a strategy provides only an overestimation of the minimal
values of $r_{0}$ for which ALD should exhibit typical phenomena for the
strongly correlated system.

Figure 9: The amplitudes of Hubbard $U$ on-site interaction as a function of
$r_{0}$ for the selected values of $\lambda_{TF}$.

We inspect the values of $1/q_{TF}\equiv\lambda_{TF}\in[2\text{
nm},10\text{nm}]$—where $\lambda_{TF}$ is the Thomas–Fermi length, which is
supported by the following estimations. First, we may consider the $q_{TF}$
resulting from the ideal 2DEG approach assuming $m^{*}=0.067m_{el}$ and
$\epsilon_{r}\approx 13$, that is, the specific values of the effective mass
and dielectric constant of the GaAs material, respectively. In this case,
since $q_{TF}$ does not depend on the carrier concentration,
$\lambda_{TF}=a_{B}^{*}/2\approx 5$ nm, where $a_{B}^{*}$ is the effective
Bohr radius. On the other hand, taking the 3D approach and
$n_{el}^{3D}=4.5\times 10^{10}/d\text{ cm}^{-2}$, one is left with
$\lambda_{TF}\approx 9$ nm, thus a value nearly twice higher than for the 2D
case. Eventually, $\lambda_{TF}$ can be estimated on the basis of a quasi-2D
model (slab) of 2DEG elaborated by Moreno and Méndez-Moreno Moreno and Méndez-
Moreno (2001). Their approach leads to $\lambda_{TF}\approx 2$ nm for the slab
width of $d=25$ nm, i.e. a value substantially lower than that of the ideal
two-dimensional electron gas. As the formulation of a general description of
screening and its role in a quasi-two-dimensional system is a complex problem
that is beyond the scope of this work, we examine the values of
$\lambda_{TF}$, relying on the estimations presented above.

Figure 10: The ratio between Hubbard on-site repulsion $U$ and the absolute
value of $t_{1}$. The dashed area corresponds to $U_{c}/|t_{1}|\approx 3.8\div
4.2$ for which the AFMI phase is supposed to emerge when considering the
Hubbard model on the honeycomb lattice Sorella  _et al._ (2012); Meng  _et
al._ (2010).

We compute interaction amplitudes of the form $V_{\mu\nu\mu\nu}$, that is, we
take into account the on-site Hubbard interaction $U\equiv V_{0,0}$ (applying
the same rule of indexing as for hoppings—see Eq. 17), as well as intersite
interactions $K_{0,Z}$ with $Z$ up to $5$. The calculation procedure has been
carried out using the Cuba library Hahn (2005); Biborski  _et al._ (2021)
designed for the integration of multi-variable functions. In Fig. 9 we present
values of $U$ as a function of radius $r_{0}$ for the selected screening
lengths $\lambda_{TF}$. According to the increasing localization of Wannier
functions with $r_{0}$, the on-site repulsion amplitudes also increase. For
example, for $\lambda_{TF}=2.5$ nm and $r_{0}=15$ nm $U\approx 0.6$ meV, while
for the same screening length but for $r_{0}=21$ nm, we obtain $U\approx 1.3$
meV. Importantly, for $\lambda_{TF}\gtrsim 5$ nm, the ratio $U/|t_{1}|\gtrsim
3.8$ for $r_{0}=15$ nm, as can be seen in Fig. 10, which may suggest that the
system is in the vicinity of the transition point between the semi-metallic
and antiferromagnetic Mott insulator (AFMI) phases Sorella  _et al._ (2012).
However, since we identify the radius $r_{0}=15$ as the one that corresponds
to the experimental setup for which there is no evidence of AFMI formation Du
_et al._ (2021), we conclude that $\lambda_{TF}$ should be considered as
$\lesssim 5$ nm in the framework of the elaborated model.

Note that since $|t_{1}|$ is one order of magnitude greater than $t_{2}$ for
the entire range of $r_{0}$, the estimate of critical radius $r_{0}^{c}$ for
which AFMI possibly emerges by inspecting $U/|t_{1}|$ is justified, since in
this view the electronic properties of the system should be similar to those
of the Hubbard model on the honeycomb lattice. Thus, taking $\lambda_{TF}=2$
nm, we find that the transition should occur at $r_{0}^{c}\lesssim 20$ nm.

Next, we examine the intersite amplitudes $K_{0,Z}$. We find that for
$\lambda_{TF}\leq 5$ nm, they are at least two orders of magnitude smaller
than the corresponding $U$ for each considered $r_{0}$. Therefore, the
interplay between $U-K$ is unlike to cause the formation of a charge density
wave (CDW), which is believed to compete with the semimetallic and AFMI phases
when the magnitude of intersite interactions is significant Wu and Tremblay
(2014). Eventually, according to the values of the dominating interaction
amplitudes, we expect a transition between the semi-metallic and AFMI phases
when the value of $r_{0}$ increases. We illustrate this trend by the results
obtained from our VMC calculations for $\lambda=2.5$ nm as presented below.

###  IV.2 Ground state of interacting Hamiltonian by means of VMC method

We formulate a minimalistic model for the interacting electrons on honeycomb
lattice described by the Hubbard Hamiltonian with hoppings extending up to the
$5th$ nn. Namely, we take into account only the on-site $U$ interactions, as
the inter-site interactions are—as mentioned—at least of two orders of
magnitude weaker. Calculations have been carried out for the lattice
consisting of $N=12\times 12$ unit cells, that is, for $288$ sites within
imposed periodic boundary conditions. The simulations have been carried out
using the mVMC software Misawa  _et al._ (2019) with the help of self-
elaborated, re-usable and generic Hamiltonian input generator Biborski (2023).
The trial state is in Pfaffian form, i.e. product of antiparallel spin pairing
terms, supplied with the on-site Gutzwiller and inter-site long-range Jastrow
projectors, which are responsible for capturing the electronic correlation
effects Biborski  _et al._ (2021). Variational parameters corresponding to
pairings, Gutzwiller and Jastrow projectors are chosen with periodicity
defined by a supercell consisting of $2\times 2$ unit cells. According to the
spin polarization of the system, we impose the constraint
$\sum_{i}^{2N}\langle\hat{S}^{z}_{i}\rangle=0$, that is, the total $z$
component of the spin is zero.

In the following, we present the calculations obtained for $\lambda_{TF}=2.5$
nm, emphasizing that it can be regarded as overestimated for the strongly
confined regime. However, as will become clear, it allows us to reveal typical
features of the strongly correlated system when $r_{0}>19$ nm.

First, we analyze the average double occupancy defined as

$\displaystyle\langle\hat{d}\rangle\equiv\frac{1}{2N}\sum_{i=1}^{2N}\langle\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}\rangle,$
(22)

where $\hat{n}_{i\sigma}$ is the particle number operator associated with the
site labeled by $i$ and $\sigma=\\{\uparrow,\downarrow\\}$ is the spin
component $z$ of the electron that resides at the site $i$.

Figure 11: The double occupancy $\langle\hat{d}\rangle$ as a function of the
AD radius $r_{0}$ (main figure) and with respect to $U/|t_{1}|$ and $t2/t1$ in
the inset figures. The estimated statistical errors are smaller than the
symbol size. The dashed area in lower-left inset corresponds to $U_{c}\approx
3.8\div 4.2$, i.e., estimated critical value for which the AFMI state emerges.

In Fig. 11 we show $\langle\hat{d}\rangle$ as a function of $r_{0}$. Its value
decreases with an increasing AD radius. This behavior is driven by the mutual
enhancement of $U$, as well as the narrowing of the occupied band. We
performed auxiliary calculations for $r_{0}=17.15,17.25$ and $17.5$ nm to
identify an abrupt drop in $\langle\hat{d}\rangle$ in vicinity of
$r_{0}\approx 17$ nm, which is the characteristic indicator of transition to
the insulating state driven by electronic correlations.

The emergence of AFMI state is also supported by the analysis of the spin-spin
correlation functions. Namely, we investigate the AF spin order parameter
provided by Sorella et al. Sorella  _et al._ (2012), but defined for the
$z$-component of $\hat{S}_{i}$, i.e.,

$\displaystyle
m_{s^{z}}^{2}\equiv\frac{1}{4N^{2}}\Bigg{\langle}\Big{[}\sum_{ij}\big{(}\hat{S}^{z}_{i,\alpha}-\hat{S}^{z}_{j,\beta}\big{)}\Big{]}^{2}\Bigg{\rangle}.$
(23)

Figure 12: The spin order parameter $m_{s^{z}}^{2}$ defined in Eq. 23 as a
function of the AD radius (main figure). The insets contain $m_{s^{z}}^{2}$
dependencies on $U/|t_{1}|$ and $t_{2}/t_{1}$. The dashed area in upper-left
inset corresponds to $U_{c}\approx 3.8\div 4.2$, i.e., estimated critical
value for which the AFMI state appears with increasing value of $U$. Estimated
statistical errors are smaller than the size of the symbols.

In Fig. 12 we present $m_{s^{z}}^{2}$ analogously to the double occupancies
shown in Fig. 11. It is clearly visible that the AF order for $r_{0}\leq 17$
nm is marginal since $m_{s^{z}}^{2}\propto 10^{-3}$. However, for $r_{0}>17$
nm, we observe its radical increase, namely, for $r_{0}=18$ nm, which is the
next case considered, it is $\sim 0.14$.

Recapitulating, when $r_{0}$ is translated into the ratio $U/|t_{1}|$ (see the
insets in Figs. 11, 12), simultaneous, anomal behaviour of both
$\langle\hat{d}\rangle$ and $m_{s^{z}}^{2}$ can be observed for
$r_{0}=17.15\div 17.25$ nm. These values of $r_{0}$ correspond to
$U/|t_{1}|=4.13\div 4.38$, the range which nearly overlap with $3.8\div 4.2$,
i.e., critical values for which transition to AFMI takes place, as it has been
estimated in the framework of Quantum Monte-Carlo approaches Sorella  _et al._
(2012); Meng  _et al._ (2010); Raczkowski  _et al._ (2020).

##  V Summary and conclusions

In this paper, we studied the anti-dot honeycomb lattice model in view of the
importance of electron-electron interactions. Within the assumed form of the
quantum dot potential and the multilevel computational scheme, we provided
evidence of the possibility of the formation of an antiferromagnetic Mott
insulator. Our modeling scheme reveals that the radius of the antidot is
crucial in this context. In addition, the inclusion of metallic screening
among electrons seems to be essential for the strength of the electron-
electron interactions, specifically when the radius of the antidot is small.
Note that we did not have to take into account the screening of the metallic
gates that induce the periodic potential Knörzer  _et al._ (2022); Biborski
_et al._ (2021); Byrnes  _et al._ (2008) as the honeycomb structure here
results from the triangular etching pattern of the upper layers. Also, we do
not take into account any kind of scattering centers, such as ionic impurities
or any other kind of disorder in the system. Therefore, the whole screening
considered here has its origin in the _metallic_ nature of 2DEGs Kollar and
Sachdev (2002). In this spirit our supporting VMC analysis reveals possible
mapping of interacting electrons in AD honeycomb lattice to the Hubbard model.
Our study provides a straightforward way to exploit the devices that are
currently fabricated for the realization of strongly correlated artificial
Dirac systems. The experimental confirmation of our predictions would open up
the opportunity for a better understanding of the electronic properties of
graphene-like systems in a controllable manner. Moreover, we indicate that in
addition to the correlation effects studied here, the AD honeycomb lattice
formed in 2DEG may also serve as a simulation platform for strain-induced
opening of the gap in honeycomb systems Cocco  _et al._ (2010); Rut  _et al._
(2023), since strain can be emulated by the distorted QD assembly. Thus, the
proposed modeling scheme may also be considered useful from this point of
view.

##  VI Acknowledgement

This work was supported by National Science Centre (NCN) agreement number
UMO-2020/38/E/ST3/00418.

## References

  * Kotov  _et al._ (2012) V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and A. H. Castro Neto, Rev. Mod. Phys. 84, 1067 (2012).
  * Tang  _et al._ (2018) H.-K. Tang, J. N. Leaw, J. N. B. Rodrigues, I. F. Herbut, P. Sengupta, F. F. Assaad, and S. Adam, Science 361, 570 (2018), https://www.science.org/doi/pdf/10.1126/science.aao2934 .
  * Das Sarma  _et al._ (2007) S. Das Sarma, E. H. Hwang, and W.-K. Tse, Phys. Rev. B 75, 121406 (2007).
  * Wu and Tremblay (2014) W. Wu and A.-M. S. Tremblay, Phys. Rev. B 89, 205128 (2014).
  * Assaad and Herbut (2013) F. F. Assaad and I. F. Herbut, Phys. Rev. X 3, 031010 (2013).
  * Sorella  _et al._ (2012) S. Sorella, Y. Otsuka, and S. Yunoki, Scientific Reports 2, 992 (2012).
  * Meng  _et al._ (2010) Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, and A. Muramatsu, Nature 464, 847 (2010).
  * Raczkowski  _et al._ (2020) M. Raczkowski, R. Peters, T. T. Phùng, N. Takemori, F. F. Assaad, A. Honecker, and J. Vahedi, Phys. Rev. B 101, 125103 (2020).
  * Otsuka  _et al._ (2016) Y. Otsuka, S. Yunoki, and S. Sorella, Phys. Rev. X 6, 011029 (2016).
  * Liebsch and Wu (2013) A. Liebsch and W. Wu, Phys. Rev. B 87, 205127 (2013).
  * Arya  _et al._ (2015) S. Arya, P. V. Sriluckshmy, S. R. Hassan, and A.-M. S. Tremblay, Phys. Rev. B 92, 045111 (2015).
  * Saleem  _et al._ (2022) Y. Saleem, A. Dusko, M. Cygorek, M. Korkusinski, and P. Hawrylak, Phys. Rev. B 105, 205105 (2022).
  * Du  _et al._ (2021) L. Du, Z. Liu, S. J. Wind, V. Pellegrini, K. W. West, S. Fallahi, L. N. Pfeiffer, M. J. Manfra, and A. Pinczuk, Phys. Rev. Lett. 126, 106402 (2021).
  * Kylänpää  _et al._ (2016) I. Kylänpää, F. Berardi, E. Räsänen, P. García-González, C. A. Rozzi, and A. Rubio, New J. Phys. 18 (2016), 10.1088/1367-2630/18/8/083014.
  * Li  _et al._ (2020) T. Li, J. Ingham, and H. D. Scammell, Phys. Rev. Res. 2, 043155 (2020).
  * Räsänen  _et al._ (2012) E. Räsänen, C. A. Rozzi, S. Pittalis, and G. Vignale, Phys. Rev. Lett. 108, 246803 (2012).
  * Krix  _et al._ (2022) Z. E. Krix, H. D. Scammell, and O. P. Sushkov, Phys. Rev. B 105, 075120 (2022).
  * Li  _et al._ (2017) Y.-M. Li, X. Zhou, Y.-Y. Zhang, D. Zhang, and K. Chang, Phys. Rev. B 96, 035406 (2017).
  * Tkachenko  _et al._ (2015) O. A. Tkachenko, V. A. Tkachenko, I. S. Terekhov, and O. P. Sushkov, 2D Materials 2, 014010 (2015).
  * Bednarek  _et al._ (2003) S. Bednarek, B. Szafran, K. Lis, and J. Adamowski, Phys. Rev. B 68, 155333 (2003).
  * Ciurla  _et al._ (2002) M. Ciurla, J. Adamowski, B. Szafran, and S. Bednarek, Physica E Low Dimens. Syst. Nanostruct. 15, 261 (2002).
  * Note (1) Note that we excluded the bottom layer of $n$-doped AlGaAs as detailed information regarding the character of doping has not been provided Du  _et al._ (2021).
  * Marzari  _et al._ (2012) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
  * Groth  _et al._ (2014) C. W. Groth, M. Wimmer, A. R. Akhmerov, and X. Waintal, New J. Phys. 16, 063065 (2014).
  * Biborski  _et al._ (2021) A. Biborski, M. P. Nowak, and M. Zegrodnik, Phys. Rev. B 104, 245430 (2021).
  * Knörzer  _et al._ (2022) J. Knörzer, C. J. van Diepen, T.-K. Hsiao, G. Giedke, U. Mukhopadhyay, C. Reichl, W. Wegscheider, J. I. Cirac, and L. M. K. Vandersypen, Phys. Rev. Res. 4, 033043 (2022).
  * Ortiz  _et al._ (1999) G. Ortiz, M. Harris, and P. Ballone, Phys. Rev. Lett. 82, 5317 (1999).
  * Moreno and Méndez-Moreno (2001) M. Moreno and R. M. Méndez-Moreno, International Journal of Quantum Chemistry 82, 269 (2001), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.1049 .
  * Hahn (2005) T. Hahn, Comput. Phys. Commun. 168, 78 (2005).
  * Misawa  _et al._ (2019) T. Misawa, S. Morita, K. Yoshimi, M. Kawamura, Y. Motoyama, K. Ido, T. Ohgoe, M. Imada, and T. Kato, Comput. Phys. Commun. 235, 447 (2019).
  * Biborski (2023) A. Biborski,  (2023), 10.5281/zenodo.8172389.
  * Byrnes  _et al._ (2008) T. Byrnes, N. Y. Kim, K. Kusudo, and Y. Yamamoto, Phys. Rev. B 78, 075320 (2008).
  * Kollar and Sachdev (2002) M. Kollar and S. Sachdev, Phys. Rev. B 65, 121304 (2002).
  * Cocco  _et al._ (2010) G. Cocco, E. Cadelano, and L. Colombo, Phys. Rev. B 81, 241412 (2010).
  * Rut  _et al._ (2023) G. Rut, M. Fidrysiak, D. Goc-Jagło, and A. Rycerz, International Journal of Molecular Sciences 24 (2023), 10.3390/ijms24021509.


# The Einstein Toolkit: A Student’s Guide

Nicholas Choustikov<EMAIL_ADDRESS>Institute of Astronomy, University of
Cambridge, Cambridge, CB3 0HA, United Kingdom

(October 21, 2020)

###### Abstract

The Einstein Toolkit represents a unique opportunity for students to explore
the world of numerical relativity, without the need for high-level computing
power or knowledge of the mathematics behind the simulations themselves. This
document aims to act as a resource for students aiming to learn about and get
involved with Numerical Relativity, focusing specifically on the modelling of
Binary Black Hole Mergers in order to produce sample Gravitational Waveforms
as would be measured by Gravitational Wave detectors such as LIGO.

The feasibility of using the Einstein Toolkit on Amazon’s Web Services (AWS)
has also been explored for the first time, confirming a potential - yet
expensive - method for students or citizen scientists to gain access to high-
level computing power sufficient to play with and run their own numerical
relativity simulations.

###### Contents

  1. 1 Introduction
  2. 2 Requirements and Feasibility
  3. 3 Installing and Building the Toolkit
  4. 4 Interacting with the Toolkit
  5. 5 Example Parameter Files
  6. 6 Diving into a Black Hole Binary
    1. 6.1 The Merger
    2. 6.2 The Model
  7. 7 Using Toolkit Data
  8. 8 Conclusion

##  1 Introduction

The Einstein Toolkit (ET) is a community-driven platform of core computational
tools [1] which provide a powerful tool to perform high-level numerical-
relativity and relativistic astrophysical simulations, even while operating
within the constraints of desktop-level computational power [1, 2]. Such a
toolkit appeals therefore as a particularly useful tool for students beginning
to learn and use Numerical Relativity (NR) within the context of university-
level physics, as well as to researches carrying out advanced simulations.

Therefore, this tutorial has been written with a keen undergraduate student in
mind. The intention being to set up said student with the ability to utilize
the toolkit, as well as to understand what goes into parameter files used for
simulations - whether for a research project or out of personal curiosity. As
a result, this tutorial will not require a strong understanding of General
Relativity (GR), though a basic understanding of the astrophysical phenomena
being studied is incredibly useful.

This tutorial will first go through the entire process of installing, building
and then using the Einstein Toolkit - assuming little knowledge of shell
scripting - or of the toolkit itself. In section 6 it will then go through a
sample parameter file (.rpar), allowing the reader an opportunity to become
familiar with parameter files as well as the Cactus Thorn structure [2]. The
process of accessing and parsing data produced following a successful run is
then discussed in section 7. For more experienced readers, it is possible to
follow the Einstein Toolkit’s own tutorial, available with an executable
jupyter notebook available at https://einsteintoolkit.org/documentation/new-
user-tutorial.html.

##  2 Requirements and Feasibility

Naturally, the equations being solved in NR are computationally expensive [3].
As a result, despite the fact that the ET has been well-optimized for even
desktop-level computational power, there are certain computational
requirements.

In principle, the Toolkit can be built and run on systems with one processor
and 4GB of RAM. As will be discussed, often high-level, accurate simulations
can easily become computationally expensive to the point of effective
impossibility on such systems, though alternatives are suggested.

Furthermore, in order to ease installation and use of the toolkit, it is
recommended to use a Linux distribution. All examples in this tutorial were
performed using a clean system running the minimal Ubuntu 20.04 distribution
[4]. For such a clean system, the following packages are required (with their
most-recent versions) in order to build and run the ET, with a sample of the
required bash code shown below:

Package | Installation Code
---|---
Git | sudo apt install git
Subversion | sudo apt-get install subversion
Curl | sudo apt install curl
Python2/3 | sudo apt install python2/3
GCC | sudo apt-get install build-essential
GFortran | sudo apt-get install gfortran
MPI | sudo apt-get install mpi

Next is a question of feasibility. As mentioned before, many high-level
simulations will require literally thousands of core-hours, and as a result
can quickly grow to become too computationally-expensive to run feasibly on a
personal-computer. For instance, using the models that will be discussed in
sections 5 and 6 as examples, the TOV Star Simulation requires approximately
10 core-minutes, while the Binary Black Hole (BBH) Merger requires 8700 core-
hours. This clearly renders it rather expensive to run on a laptop. Therefore,
in cases where the user does not have access to a computing cluster, the best
option - despite being fiscally expensive - might potentially be to make use
of the Amazon Web Services (AWS) cloud-computing service.

Available to members, the AWS Elastic Computing Cloud (EC2) is a service which
provides highly-customizable computing capacity in the cloud. This is done, by
first choosing the ’shape’ of a desired virtual machine (VM). This decides the
CPU and RAM size available to the machine, from a list of options pre-
determined by AWS. Following this, an image (custom or default) is chosen
which will run on the VM. For the purpose of this project, a public image
(freely available) running a clean installation of Ubuntu 20.04 [4] was
chosen, to be run on an m5.16xlarge instance. This setup offered performance
equivalent to 64 CPUs with 256 GiB memory [5], at the cost of $3.072 per hour.
This though can quickly become equally prohibitively expensive. In order to
combat this, a system was setup with a large detachable volume which could be
attached to the primary instance while calculations were being run, and then
attached to a much smaller (and cheaper) instance in order to then carry out
data analysis.

This option can equally become particularly expensive, but it certainly offers
an opportunity for students with external funding to play with computing power
previously only available in dedicated computing clusters.

During the process of writing this document, several models of BBH mergers,
with varying mass ratios, initial separations and intrinsic spins were started
and run through the first 20,000 iterations. These yielded the expected
behaviour and waveforms up to this point in the model. None were however run
to completion due to a lack of sufficient funding.

##  3 Installing and Building the Toolkit

Once all necessary modules have been acquired, the following steps should be
taken in order to install the ET. It is important to note, that in cases where
only one user will be utilising the toolkit, that user itself should run all
setup code. If on the other hand multiple users will be interacting with the
toolkit, then the root user should not be used. This is specifically important
due to the fact that by default, all simulation data will be stored in the
user’s home directory, as described in section 7.

To start the process, the following three commands should be called, in the
directory within which the toolkit will be installed:

⬇

1 curl -kl0
https://raw.githubusercontent.com/gridaphobe/CRL/ET_2020_05/GetComponents

2 chmod a+x GetComponents

3 ./GetComponents \--parallel
https://bitbucket.org/einsteintoolkit/manifest/raw/ET_2020_05/einsteintoolkit.th

This code downloads the GetComponents package, which is then executed in order
to install the ET itself. Following this, chmod allows the package to be
executed, and the third line executes it, providing the einsteintoolkit.th
file as an argument, to be used as the thornlist. Once this is done, in
principle all component thorns of the ET have been acquired, they just need to
be configured for a given machine, such that they will be able to interact
with machine-specific queuing systems.

This is principally done using the Simulation Factory (”SimFactory”), which
itself is configured as shown:

⬇

1 cd Cactus

2 ./simfactory/bin/sim setup-silent

Once this step has been completed, SimFactory has been configured to run using
the machine’s default settings, which can be accessed and edited in
./simfactory/mdb/machines/<hostname>.ini. Any changes here should be carried
out with care, as they may lead to warnings or errors when simulations are
started-up. Furthermore, it is important to note that by default the toolkit
is configured to run with only $1$ node.

Now, the ET must be built. Assuming that there are no issues with the
SimFactory setup, and that all required packages are functioning properly, the
next command should completely build the toolkit.

⬇

1 ./simfactory/bin/sim build -j<N> --thornlist ../einsteintoolkit.th

where $N$ is the number of processors used when building, and should be
changed to suit the given machine. It should be noted that even with multiple
processors, this step can take a very long time. For example, with $N=2$, it
took approximately three and a half hours to complete, while with $N=64$, it
took fewer than 15 minutes to completely build.

Assuming this step finished with no issues, the toolkit has now been built,
and is ready to run several test examples, as described in section 5.

##  4 Interacting with the Toolkit

The toolkit is made up of a series of Cactus components, or thorns, which for
all intents and purposes can be thought of as modules which can be called upon
to accomplish various tasks. The complete thornlist along with all
documentation is available at
https://einsteintoolkit.org/documentation/ThornGuide.php.

Tasks achieved by these thorns include multiple vacuum spacetime solvers,
relativistic hydrodynamic solvers, initial data components, along with others
for data output and analysis.

Therefore, in order to run a simulation, it is necessary to write a parameter
file (.par or .rpar) which tells the toolkit which thorns to use, along with
the desired parameters - all working together to produce the required model.
As a result, it quickly becomes clear that dependant on the system being
investigated, different thorns will prove far more efficient and useful -
especially compared to thorns which on the surface appear to accomplish the
same task.

Once a parameter file has been written, all interaction with the ET can be
done using SimFactory. In principle, the user will instruct SimFactory to
create and run a simulation, following a given parameter file. This can be
achieved using the following command:

⬇

1 lamboot

2 ./simfactory/bin/sim create-run <SIM_NAME> <PARFILE> <N> <WALLTIME>

It is important to note that when using a computing cluster, it is necessary
to call lamboot every time an instance running the ET is started, as otherwise
errors will be encountered when the simulation is started.

This code will create a simulation under the name <SIM_NAME>, and then
proceeds to run it, using <N> processors with a run time of <WALLTIME>.

More useful, is an equivalent command to create and submit a given simulation:

⬇

1 ./simfactory/bin/sim create-submit <SIM_NAME> <PARFILE> <N> <WALLTIME>

this will add the simulation to the queue, telling SimFactory to begin running
it as soon as there are sufficient resources available. It is also
particularly useful, as in doing so the simulation will run in the background
- therefore allowing the user to do something else in the meantime. Such as:

⬇

1 ./simfactory/bin/sim list-simulations

2 ./simfactory/bin/sim show-output <SIM_NAME> \--follow

3 ./simfactory/bin/sim stop <SIM_NAME>

as expected, the first command will list all simulations in the toolkit’s
memory, along with their statuses (RUNNING/FINISHED). The next command will
show the full terminal output of the simulation, with live output. This can be
stopped with <Ctrl + C> without stopping the simulation itself. Alternatively,
the \--follow parameter can be removed, in which case the entire output
history of the simulation up to that point will be printed. This allows the
user to see how the simulation is proceeding. Finally, the third command will
terminate the simulation inelegantly, if required.

There are a plethora of other commands available for working with SimFactory,
and these can be found and explored on the SimFactory documentation [6].

##  5 Example Parameter Files

As is historically known, often the most simple test of a running system is to
produce a Hello World output. For this purpose, there is HelloWorld.par, which
is found within the toolkit at
$\sim$/Cactus/arrangemets/CactusExamples/HelloWorld/par/HelloWorld.par. On
running this ”simulation”, if the output ”Hello World” is seen printed
anywhere, it can be considered a successful test - and as a result the toolkit
is considered to be running correctly.

Following this, a further - more relevant - example simulation is that of a
simple TOV Neutron Star (NS).

Here, the simulation evolves a single, stable, non-rotating NS using a
solution to the Tolman-Oppenheimer-Volkoff equation [7]

$\frac{dP}{dr}=-\frac{Gm}{r^{2}}\rho\bigg{(}1+\frac{P}{\rho
c^{2}}\bigg{)}\bigg{(}1+\frac{4\pi
r^{3}P}{mc^{2}}\bigg{)}\bigg{(}1-\frac{2Gm}{rc^{2}}\bigg{)}^{-1}$ (1)

as the initial condition. Furthermore, the model assumes a polytropic equation
of state $P=K\rho^{\Gamma}$, with $\Gamma=2$ and $K=100$.

Despite the fact that this is a fairly low-resolution model, it serves as a
reasonable representation of a non-rotating NS with a mass of $M=1.4M_{\odot}$
and as a result a reasonable test for the ET, without being too
computationally intensive. The simulation’s maximum density ($\rho_{c}$) has
been plotted in figure 1. Here, the observed oscillations are produced by
interactions between the star and the artificial atmosphere [1].

The parameter file used, along with a further description and documentation is
available at https://einsteintoolkit.org/gallery/ns/index.html.

Figure 1: Maximum internal density within a non-rotating neutron star of mass
$M=1.4M_{\odot}$ plotted as a function of time. Model was produced using the
Einstein Toolkit, by numerically evolving a system from initial conditions
produced as a spherically symmetric solution to the TOV equation for a compact
neutron star.

##  6 Diving into a Black Hole Binary

###  6.1 The Merger

Since the advent of effective numerical relativity, one particularly
successful topic of research has been that of BBH mergers. Such systems are
typically produced by either dynamical means, whereby dynamical interactions
between black holes in dense clusters can lead to binaries forming [8], or by
non-dynamical means, whereby binaries made up of two massive stars survive as
a close binary due to their tidally locked high spin, until they both collapse
into black holes [9].

These systems evolve as follows. First, the two black holes exist as a
classical binary, with the two compact bodies orbiting each other. As they
orbit, Gravitational Waves (GW) are emitted, acting as an energy-loss function
[10]. This leads to the two bodies orbiting closer to and faster around each
other, in the inspiral phase. Eventually, the two collide and merge -
producing the GW ’chirp’ [11] \- lasting a fraction of a second - followed by
the ringdown phase [12]. Measurements of the GWs produced throughout the
system’s evolution can yield valuable information about the two bodies
involved, as well as the exact nature of their orbit. This data can therefore
provide an extraordinary opportunity to test GR under such high-energy
conditions [13].

In principle, the inspiral and ringdown can be modelled using Post-Newtonian
approximations [14] and black hole perturbation theory [15] respectively.
However in order to fully model the entire evolution of a given system, the
full equations of GR must be solved, and as a result NR truly comes into its
own.

Obtaining accurate models of BBH mergers has become incredibly important, with
the recent construction of gravitational wave observatories such as LIGO [16]
and VIRGO [17] paving the way for several observational runs for which sample
waveforms are naturally needed. One of the most important results came in 2015
with the detection of GW150914 [11], the first successful detection of
gravitational waves due to a BBH merger, coming a century after Einstein’s
prediction. This system involved two black holes with initial inferred masses
of $36M_{\odot}$ and $29M_{\odot}$ at a distance of
$410\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$ from Earth [18]. It is a model of
this very event which is to be discussed here.

###  6.2 The Model

The simulation discussed in this section can be found at
https://einsteintoolkit.org/gallery/bbh/index.html, with the simulation’s
parameter file available at
https://bitbucket.org/einsteintoolkit/einsteinexamples/raw/master/par/GW150914/GW150914.rpar
[19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29].

Investigating the parameter file, it is clear to see that the file’s preamble
defines the system’s initial data, with the following values:

Variable | Value | Meaning
---|---|---
$D$ | $10$ | Initial puncture separation
$q$ | $36.0/29.0$ | Mass ratio $(q\geq 1)$
M | $1.0$ | Total mass
$\vec{\chi}_{p}$ | $[0,0,0.31]$ | Primary BH spin
$\vec{\chi}_{m}$ | $[0,0,-0.46]$ | Secondary BH spin
$p_{r}$ | $-0.00084541527$ | Radial linear momentum
$p_{\phi}$ | $0.09530152297$ | Azimuthal linear momentum

It is important to note, that all calculations within this parameter file -
and by extension - within most of the toolkit, are carried out in Geometrized
units [30], whereby the convention is taken that $G=c=1$, therefore implying
that all distances and times are given in units of mass. Furthermore, the
convention is typically taken that the total mass ($M=m_{p}+m_{m}$) is equal
to one. The mass ratio ($q=\frac{m_{p}}{m_{m}}$) is one of great importance,
influencing the behaviour of the inspiral and subsequent merger in a profound
manner. Here it is defined to be greater than one. $\chi_{p}\text{ and
}\chi_{m}$ are the three-dimensional dimensionless intrinsic spins of the two
black holes respectively, defined by $\vec{\chi}=\vec{J}/M^{2}$ [31], where
$\vec{J}$ is the black hole’s angular momentum and it is known that
$0\leq|\vec{\chi}|\leq 1$.

Following this, the remainder of the preamble sets up various necessary values
such as the estimated duration of the merger, time steps of each iteration as
well as setting up the grid structure used during the simulation itself.

The parameter file then lists all thorns required for the simulation,
activating and subsequently specifying operational parameters. These are
described within the thorn guide, provided at
https://einsteintoolkit.org/documentation/ThornGuide.php. Here, a few
important thorns will be discussed in order to paint a picture of how the
simulation will proceed.

Several of the thorns, including Coordinates, Interpolate2 and
GlobalDerivative are parts of the llama multi-block infrastructure [21]. The
purpose of these thorns is to work in unison to allow the use of multiple
curvilinear coordinate grid patches for Cactus, using grids provided by
Carpet. Here, Thornburg04 is a system consisting of 7 patches. These include a
central cubic Cartesian grid capable of adaptive mesh refinement, surrounded
by 6 ”inflated cube” spherical grids, together covering a spherical shell.
WaveExtractCPM is also particularly important, providing a method to extract
gauge-invariant perturbations of the Schwarzschild metric away from the binary
itself. These data then allow for the study of gravitational waves being
emitted throughout the simulation’s evolution.

As mentioned above, Carpet is a mesh refinement driver [26] which is able to
run with multiple grid patches. It runs by actively altering the refinement of
areas of the grid, providing greater sensitivity and accuracy to sensitive or
turbulent regions - as compared to other, less interesting ones.

The AHFinderDirect thorn quickly and accurately locates apparent horizons of
black holes based on provided initial guesses [25]. The local apparent horizon
is the surface boundary which separates outward directed light rays travelling
outwards from outward directed light rays travelling inwards. It therefore
acts as the local notion of a black hole’s boundary, giving a good description
of the puncture’s position.

TwoPunctures is a thorn which provides initial data for a system with two
black holes [27] (two punctures) located along the x-axis, with opposite
linear momenta along the y-axis. It does so by utilising a spectral method in
solving the necessary Hamiltonian constraint equations as described by [27].

The QuasiLocalMeasures thorn allows for the real-time calculation of quasi
local measures such as masses, linear and angular momenta on a closed two-
dimensional surface [29]. These data are then outputted to allow the
simulation’s progress to be tracked. Following this, the ML_BSSN thorn is used
to evolve the simulation, by iteratively solving a BSSN-type formulation of
the Einstein Field Equations [23], using a collection of fourth-order accurate
integration methods.

Unsurprisingly, the Output section gives details of which data is produced, at
which intervals. Furthermore, it is specified which format is used to store
the data. Specifically, IOASCII and IOHDF5 are used, with the rationale
described in the next section.

##  7 Using Toolkit Data

Simulations such as the BBH Merger above produce a huge amount of data, with
the full simulation producing several terabytes of data [19]. However, even
with lightweight data (comprising a small number of iterations) plenty of
information can be gleamed, as shown in the figures below.

In general, all data produced in simulations using the Einstein Toolkit are
saved in one of two formats: ASCII (.asc) and HDF5 (.h5). Both formats have
their advantages and disadvantages and are therefore used, though for
different data types.

Data stored as .asc files are immediately accessible on all computing
platforms, appearing as a typical .txt file. These can immediately be parsed
by python, using a tool such as the pandas package [32]. Although this is
doable, it is recommended to use a short script to convert the data into a
.csv (comma separated values) file, for ease of use:

⬇

1 #!/bin/sh

2 cat puncturetracker-pt_loc..asc | grep "^[^#;]" > foo1.txt

3 sed -e ’s/ /\t/g’ foo1.txt > foo2.txt

4 sed -e ’s/\t/,/g’ foo2.txt > puncturetracker-pt_loc.csv

5

6 echo done

Here, data produced by the PunctureTracker thorn is converted into a .csv
format, while simultaneously formating the data to be comma-delimited and
removing all commented lines. This modified file can then be read in using the
following lines in python:

⬇

1 import pandas as pd

2 data = pd.read_csv("puncturetracker-pt_loc.csv", names=headers)

3 time = pd.DataFrame.to_numpy(data[’time’])

where headers is a 1x52 array containing the data column headers as described
in the original puncturetracker-pt_loc..asc file. The final line then extracts
the ’time’ column, converting it into a numpy array - which can then easily be
used for any further analysis.

On the other hand, data saved as a .h5 file are stored in a hierarchical
format, which allows for the storage of huge amounts of data with a high
degree of compression. Furthermore, this format allows for very fast
Input/Ouput (I/O) as compared to ASCII. However, in order to visually inspect
and edit data stored in a .h5 format, it is necessary to use an interface tool
such as HDFView [33].

In order to parse .h5 data using python, the h5py package [34] is required.
The following code will then allow you to read in data:

⬇

1 import h5py

2

3 path = ’GW150914_28/output-0000/GW150914_28/mp_psi4.h5’

4 filename = l2_m2_r500.00

5 hf = h5py.File(path, ’r’)

6 data = hf[filename]

This code will open the example simulation’s Weyl Scalar ($\Psi_{4}$) [35]
information and extract data for the (2,2) spherical harmonic mode. The
resultant array will have three columns - giving the Real and Imaginary
components as a function of time.

Following this, it is clearly possible to plot and visualize the simulation
data. Here, several examples are shown using the above BBH simulation data
[19]. Data was parsed using code similar to that which was shown above, and
all post-processing and plotting were done using the freely accessible
packages numpy, scipy and matplotlib.pyplot [36, 37, 38].

Figure 2: Several plots of data produced in a simulation of the binary black
hole merger GW150914 using the Einstein Toolkit. Top left: Movement of the
black holes’ centroid throughout the simulation is shown. The inspiral stage
is clearly seen, with the rate of convergence increasing as the two black
holes approach one another. Top right: Here the precise radial separation
between the two black holes is shown as a function of time. This data was then
fitted, with the result having the expected form. Bottom: Raw Weyl Scalar data
for the dominant (2,2) spherical harmonic mode, as measured at a distance of
$r=500M$. This shows the expected increase in both frequency and amplitude
culminating in the merger itself.

Figure 2 shows several plots of the raw simulation data. These include the
motion of the system’s centroid, the change in radial separation as a function
of time as well as raw data from the Weyl Scalar for the dominant (2,2)
spherical harmonic mode. The first of these clearly demonstrates the expected
quasicircular nature of the inspiral, with an innermost stable circular orbit
of $r_{\text{ISCO}}\approx 6.6M$. The second shows a rapid increase in the
rate of radial approach as a function of separation. This data was fitted to
an equation of the form: $r(t)=r_{0}(1-t/t_{\text{coalesce}})^{b}$ with the
linear regression algorithm yielding $b=0.2361\pm 0.0003$. For binary systems
it is possible to show that the averaged power emitted by gravitational waves
is given by [39]

$\frac{dE}{dt}=-\frac{32}{5}\frac{\mu^{2}M^{3}}{r^{5}}F(e)=\frac{32}{5}\frac{M_{1}^{2}M_{2}^{2}(M_{1}+M_{2})}{r^{5}}F(e),$
(2)

where the factor

$F(e)=(1-e^{2})^{-7/2}\bigg{(}1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\bigg{)}.$
(3)

Here, taking the Newtonian approximation that $E=\frac{-M_{1}M_{2}}{2r}$ and
that $e\approx 0$ gives

$\frac{dr}{dt}=-\frac{64}{5}\frac{M_{1}M_{2}(M_{1}+M_{2})}{r^{3}},$ (4)

which can easily be integrated to give
$r(t)=r_{0}(1-t/t_{\text{coalesce}})^{1/4}$. Therefore, the simulation’s
orbital separation decreases slightly faster than this approximation, as would
be expected given the relativistic nature of the system.

Following this, the bottom plot of figure 2 shows the variation in the Weyl
Scalar of the dominant (2,2) spherical harmonic mode, as measured at a
distance of $r=500M$. This demonstrates the expected sharp rises in both
frequency and amplitude, as well as the complete ’chirp’ [11] itself, as the
two black holes merge.

In order to calculate the measurable strain ($h$) - as would be measured by a
ground-based detector such as LIGO - we use the following formalism [40]

$h=h_{+}-ih_{\text{$\times$}}=\int_{-\infty}^{t}dt^{\prime}\int_{-\infty}^{t^{\prime}}dt^{\prime\prime}\Psi_{4},$
(5)

where $h_{+}$ and $h_{\times}$ are the plus and cross polarizations
respectively, as they appear in the transverse-traceless gauge [41]. In order
to carry out this integration, and to account for a linear drift we use the
procedure outlined in [40]. Doing so yields figure 3 below.

Furthermore, for figure 3 we have accounted for the propagation time taken for
the signal to reach the detector at $r=500M$. This was done by approximating
the propagation time by $t\approx r^{*}$, where $r^{*}$ is the Regge-Wheeler
”tortoise coordinate” [42]

$r^{*}=r+2M\ln\bigg{|}\frac{r}{2M}-1\bigg{|}.$ (6)

Figure 3: Plot of the dominant $h_{+}^{2,2}$ strain, as would be measured by a
gravitational wave observatory such as LIGO. Key features include the
constant-amplitude and frequency waves produced by the initial inspiral, the
rapidly increasing frequency and amplitude of the merger itself, as well as
the equally rapid ringdown, showing the aftermath of the binary black hole
merger. The entire signal lasted 0.2s, with a peak gravitational wave
frequency of 250Hz.

Figure 3 shows the build-up and chirp in strain of the $h_{+}^{2,2}$ spherical
harmonic mode, as produced by the GW150914 inspiral, merger and ringdown.
These stages are characterized by the initial approximately constant amplitude
and frequency, followed by a rapid increase in both frequency and amplitude,
followed by an equally rapid drop in both during the ringdown phase.
Accounting for units, the entire signal would have lasted 0.2 seconds, with a
peak amplitude at a frequency of 150Hz - as detected by LIGO [11].

Following the merger, the resultant black hole was inferred to have a total
mass of $62^{+4}_{-4}M_{\odot}$, implying that
$3.0^{+0.5}_{-0.5}M_{\odot}c^{2}$ was radiated as gravitational waves [18],
such as those shown in figure 3. The successful detection of GW150914 was a
landmark event in gravitational wave astronomy. It provided the first evidence
of the existence of binary stellar-mass black hole systems, as well as being
the first direct detection of gravitational waves.

##  8 Conclusion

To conclude, this document has described the process of installing, building
and subsequently using the Einstein Toolkit, from the perspective of a
student. The aim being to dispel some of the ’mystery’, and to show students
that the world of numerical relativity is not only fascinating, but entirely
feasible to start exploring, without the need for a computing cluster, or even
heavy technical skills.

Furthermore, the document has explored several examples of numerical
simulations done using the toolkit, with the focus being on a model of a
binary black hole merger. This was done using system properties inferred by
LIGO of the GW150914 event. Several plots were then produced - with the
process described and explained.

The feasibility of utilising AWS computing infrastructure in order to carry
out computational physics has also been considered. Though an expensive tool,
AWS is able to provide an extraordinary amount of computing power - which a
student would have otherwise found difficult to obtain. This being the case,
more work must to be done to benchmark AWS capabilities, in order to ascertain
the optimal setup for uses of this kind.

As a result, it has been reiterated that the Einstein Toolkit is a powerful
tool, allowing for both high-level simulations as well as experimentation with
numerical relativity ’at home’. It has also been shown that it is not only
entirely feasible for students to use the Toolkit, but equally that it is an
incredible opportunity for those students to learn and gain experience of
using numerical relativity.

## References

  * [1] F. Löffler, J. Faber, E. Bentivegna, T. Bode, P. Diener, R. Haas, I. Hinder, B. C. Mundim, C. D. Ott, E. Schnetter, and et al., “The einstein toolkit: a community computational infrastructure for relativistic astrophysics,” Classical and Quantum Gravity, vol. 29, p. 115001, May 2012.
  * [2] M. ZILHÃO, M. LÖFFLER, FRANK Zilhão, and F. Löffler, “An introduction to the einstein toolkit,” International Journal of Modern Physics A, vol. 28, p. 1340014, Sep 2013.
  * [3] J. Blackman, S. E. Field, M. A. Scheel, C. R. Galley, C. D. Ott, M. Boyle, L. E. Kidder, H. P. Pfeiffer, and B. Szilágyi, “Numerical relativity waveform surrogate model for generically precessing binary black hole mergers,” Physical Review D, vol. 96, Jul 2017.
  * [4] “https://releases.ubuntu.com/20.04/,” Sept. 2020.
  * [5] “https://aws.amazon.com/ec2/instance-types/,” Sept. 2020.
  * [6] E. Schnetter. and M. Thomas, “http://simfactory.org/info/documentation/userguide/use.html,” Sept. 2020.
  * [7] J. R. Oppenheimer and G. M. Volkoff, “On massive neutron cores,” Physical Review, vol. 55, 1939.
  * [8] S. Banerjee, “Stellar-mass black holes in young massive and open stellar clusters and their role in gravitational-wave generation,” Monthly Notices of the Royal Astronomical Society, p. stw3392, Jan 2017.
  * [9] P. Marchant, et al., “A new route towards merging massive black holes,” Astronomy & Astrophysics, vol. 50, 2016.
  * [10] É. É. Flannagan and S. A. Hughes, “The basics of gravational wave theory,” New Journal of Physics, vol. 7, 2005.
  * [11] B. P. Abbott, et al., “Observation of gravitational waves from a binary black hole merger,” Physical Review Letters, 2016.
  * [12] S. Bhagwat, et al, “On choosing the start time of binary black hole ringdowns,” Physical Review D, vol. 97, May 2018.
  * [13] P. Amaro-Seoane, et al., “Research update on extreme-mass-ratio inspirals,” Journal of Physics: Conference Series, vol. 610, 2015.
  * [14] C. M. Will, “On the unreasonable effectiveness of the post-newtonian approximation in gravitational physics,” PNAS, vol. 108, no. 15, pp. 5938–5945, 2011.
  * [15] A. Le Tiec, “The overlap of numerical relativity, perturbation theory and post-newtonian theory in the binary black hole problem,” International Journal of Modern Physics D, vol. 23, p. 1430022, Sep 2014.
  * [16] B. P. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen, R. S. Amin, S. B. Anderson, W. G. Anderson, M. A. Arain, and et al., “Ligo: the laser interferometer gravitational-wave observatory,” Reports on Progress in Physics, vol. 72, p. 076901, Jun 2009.
  * [17] T. Accadia, et al., “Virgo: a laser interferometer to detect gravitational waves,” Journal of Instrumentation, vol. 7, Mar. 2012.
  * [18] B. P. Abbott, et al., “Properties of the binary black hole merger gw150914,” Physical Review Letters, 2016.
  * [19] B. Wardell, I. Hinder, and E. Bentivegna, “Simulation of GW150914 binary black hole merger using the Einstein Toolkit,” Sept. 2016.
  * [20] F. Loffler et al., “The Einstein Toolkit: A Community Computational Infrastructure for Relativistic Astrophysics,” Class. Quant. Grav., vol. 29, p. 115001, 2012.
  * [21] D. Pollney, C. Reisswig, E. Schnetter, N. Dorband, and P. Diener, “High accuracy binary black hole simulations with an extended wave zone,” Phys. Rev., vol. D83, p. 044045, 2011.
  * [22] M. Thomas and E. Schnetter, “Simulation factory: Taming application configuration and workflow on high-end resources,” in Grid Computing (GRID), 2010 11th IEEE/ACM International Conference on, pp. 369 –378, oct. 2010\.
  * [23] J. D. Brown, P. Diener, O. Sarbach, E. Schnetter, and M. Tiglio, “Turduckening black holes: An Analytical and computational study,” Phys. Rev., vol. D79, p. 044023, 2009.
  * [24] S. Husa, I. Hinder, and C. Lechner, “Kranc: a Mathematica application to generate numerical codes for tensorial evolution equations,” Comput. Phys. Commun., vol. 174, pp. 983–1004, 2006.
  * [25] J. Thornburg, “A Fast apparent horizon finder for three-dimensional Cartesian grids in numerical relativity,” Class. Quant. Grav., vol. 21, pp. 743–766, 2004.
  * [26] E. Schnetter, S. H. Hawley, and I. Hawke, “Evolutions in 3-D numerical relativity using fixed mesh refinement,” Class. Quant. Grav., vol. 21, pp. 1465–1488, 2004.
  * [27] M. Ansorg, B. Bruegmann, and W. Tichy, “A Single-domain spectral method for black hole puncture data,” Phys. Rev., vol. D70, p. 064011, 2004.
  * [28] T. Goodale, G. Allen, G. Lanfermann, J. Massó, T. Radke, E. Seidel, and J. Shalf, “The Cactus framework and toolkit: Design and applications,” in Vector and Parallel Processing – VECPAR’2002, 5th International Conference, Lecture Notes in Computer Science, (Berlin), Springer, 2003.
  * [29] O. Dreyer, B. Krishnan, D. Shoemaker, and E. Schnetter, “Introduction to isolated horizons in numerical relativity,” Phys. Rev., vol. D67, p. 024018, 2003.
  * [30] A. L. Myers, “Natural units of general relativity,” 2016.
  * [31] C. S. Reynolds, “Observing black holes spin,” Nature Astronomy, vol. 3, 2019\.
  * [32] W. McKinney, “Data structures for statistical computing in python,” in Proceedings of the 9th Python in Science Conference (S. van der Walt and J. Millman, eds.), pp. 51 – 56, 2010.
  * [33] “https://www.hdfgroup.org/downloads/hdfview/,” Sept. 2020.
  * [34] A. Collette, “https://www.h5py.org/,” Sept. 2020.
  * [35] L. A. Gómez López and G. D. Quiroga, “Asymptotic structure of spacetime and the newman-penrose formalism: a brief review,” Revista Mexicana de Física, vol. 63, 2017.
  * [36] J. D. Hunter, “Matplotlib: A 2d graphics environment,” Computing in Science & Engineering, vol. 9, no. 3, pp. 90–95, 2007.
  * [37] C. R. Harris, et al., “Array programing with numpy,” Nature, vol. 585, 2020.
  * [38] P. Virtanen, et al., “Scipy 1.0: fundamental algorithms for scientific computing in python,” Nature Methods, vol. 17, 2020.
  * [39] M. Celoria, R. Oliveri, A. Sesana, and M. Mapelli, “Lecture notes on black hole binary astrophysics,” 2018.
  * [40] C. Reisswig, D. Pollney, “Notes on the integration of numerical relativity waveforms,” Classical and Quantum Science, 2011.
  * [41] P. Wagle, et al., “Polarization modes of gravitational waves in quadratic gravity,” Physical Review D, 2019.
  * [42] S. J. Crothers, “On the regge-wheeler tortoise and the kruskal-szekerer coordinates,” Progress in Physics, vol. 3, 2006.