# Towards high-throughput many-body perturbation theory: efficient algorithms
and automated workflows

Miki Bonacci<EMAIL_ADDRESS>FIM Department, University of Modena
and Reggio Emilia, Via Campi 213/a, Modena, Italy  S3 Center, Istituto
Nanoscienze, CNR, Via Campi 213/a, Modena, Italy  Junfeng Qiao  Theory and
Simulation of Materials (THEOS) and National Centre for Computational Design
and Discovery of Novel Materials (MARVEL), École Polytechnique Fédérale de
Lausanne, CH-1015 Lausanne, Switzerland  Nicola Spallanzani  S3 Center,
Istituto Nanoscienze, CNR, Via Campi 213/a, Modena, Italy  Antimo Marrazzo
Dipartimento di Fisica, Università di Trieste, I-34151 Trieste, Italy
Giovanni Pizzi  Theory and Simulation of Materials (THEOS) and National Centre
for Computational Design and Discovery of Novel Materials (MARVEL), École
Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland  Laboratory
for Materials Simulations (LMS), Paul Scherrer Institut (PSI), CH-5232
Villigen PSI, Switzerland  Elisa Molinari  FIM Department, University of
Modena and Reggio Emilia, Via Campi 213/a, Modena, Italy  S3 Center, Istituto
Nanoscienze, CNR, Via Campi 213/a, Modena, Italy  Daniele Varsano  S3 Center,
Istituto Nanoscienze, CNR, Via Campi 213/a, Modena, Italy  Andrea Ferretti  S3
Center, Istituto Nanoscienze, CNR, Via Campi 213/a, Modena, Italy  Deborah
Prezzi  S3 Center, Istituto Nanoscienze, CNR, Via Campi 213/a, Modena, Italy

###### Abstract

The automation of ab initio simulations is essential in view of performing
high-throughput (HT) computational screenings oriented to the discovery of
novel materials with desired physical properties. In this work, we propose
algorithms and implementations that are relevant to extend this approach
beyond density functional theory (DFT), in order to automate many-body
perturbation theory (MBPT) calculations. Notably, a novel algorithm pursuing
the goal of an efficient and robust convergence procedure for GW and BSE
simulations is provided, together with its implementation in a fully automated
framework. This is accompanied by an automatic GW band interpolation scheme
based on maximally-localized Wannier functions, aiming at a reduction of the
computational burden of quasiparticle band structures while preserving high
accuracy. The proposed developments are validated on a set of representative
semiconductor and metallic systems.

##  I Introduction

Computational HT screening is nowadays a consolidated approach to materials
discovery curtarolo_high-throughput_2013 ; vecchio_high-throughput_2021 ;
luo_highthroughput_2021 , as a complementary and accelerated tool with respect
to experimental efforts. In the last decade, seminal works in this field have
addressed, among many other topics, the discovery of novel 2D materials
ashton_topology-scaling_2017 ; cheon_data_2017 ; Mounet_2018 ; choud_2d ;
MC2DB ; haastrup_computational_2018 ; marrazzo_relative_2019 , the
identification of optimal new lithium-ion battery anodes kirklin_high-
throughput_2013 ; zhang_computational_2019 , thermoelectric
chen_understanding_2016 ; bhattacharya_high-throughput_2015 , photocatalysts
castelli_computational_2012  and photovoltaic light harvesting
yu_identification_2012 ; yan_solar_2017 ; kuhar_high-throughput_2018
materials. The success of these studies relies on the development of different
software and implementations lejaeghere_reproducibility_2016  that were able
to encode complex domain-specific knowledge into automated and robust
workflows – enforcing rigorous computational protocols
lejaeghere_reproducibility_2016 ; SSSP  and managing all the steps concerning
a simulation – which thus require the least possible human intervention
maffioletti_gc3pie:_2012 ; curtarolo_aflow:_2012 ; jain_fireworks:_2015 ;
PIZZI2016218 ; hjorth_larsen_atomic_2017 ; mathew_atomate:_2017 ;
mortensen_myqueue:_2020 ; huber_aiida_2020 ; uhrin_workflows_2021 .

Concerning the electronic structure field, most of these works and
implementations are based on DFT, which allows one to compute total energies,
optimized geometries, and other ground state properties of materials with
predictive accuracy. However, different approaches are required for the
accurate prediction of excited-state properties of materials, such as
quasiparticle (QP) band structures and absorption spectra, which are typically
crucial for the description of active processes in modern optoelectronic
technologies, like photovoltaics, photocatalysis, light-emitting diodes
(LEDs), photodetectors and solar cells graph_opto ; MXenes_opto_review ;
pathways_photovoltaics ; zhu_fotoc ; photodetectors_1D_CdSe ;
xia_ultrafast_2009 ; lee_highdetectivity_2019 . In this context, MBPT and
Green’s function methods represent the state-of-the-art tools, where charged
(electronic quasi-particle levels) and neutral excitations (optical
properties, electron energy loss spectra) can be obtained by means of the GW
approximation and the Bethe-Salpeter equation (BSE), respectively
martin_reining_ceperley_2016 .

To date, a limited number of attempts have been made toward automation
van_setten_automation_2017  and HT screening haastrup_computational_2018 ;
rasmussen_towards_2021  based on these MBPT approaches, mainly because of
their conceptual and computational complexity. Indeed, depending on the
specific physical problem, different levels of theory might be adopted
huser_quasiparticle_2013 , further branching off depending on the chosen
approximations and implementations rangel_repr_2020 ; stankovski_g_2011 .
Furthermore, even for the simplest approximations, these calculations require
the control over a much larger parameter space with respect to DFT, with
parameters that are often interdependent and might change depending on the
specific implementations adopted, but the choice of which is always crucial to
obtain reliable results. From a computational point of view, this type of
simulations are constituted by a chain of distinct steps, often relying on the
usage of different software tools, each of them with their own specificity,
e.g., in terms of memory and parallelization requirements. Memory requirements
are much heavier than in standard DFT simulations even for moderate system
size, and often require massive usage of parallel computing resources.
Calculations often fail due to memory overflow and have to be restarted with
careful choice of parameters. All of these problems make the application of
MBPT-based approaches a complex and difficult task per se, and its automation
still an open challenge.

Building on pioneering works in the field van_setten_automation_2017 ;
rasmussen_towards_2021 , we here focus on the development of an improved
algorithm aiming at an efficient and computationally cost-effective management
of the choice of converged parameters for accurate GW-BSE calculations. The
algorithm is implemented in the AiiDA framework huber_aiida_2020 ;
uhrin_workflows_2021 , a platform that is routinely used for HT studies
Mounet_2018 ; mercado_silico_2018 ; prandini_precision_2018 ; Vitale2020  and
that incorporates the ADES model for Automation, Data, Environment and Sharing
PIZZI2016218 . As detailed below, this implementation allows us to encode an
efficient error handling, memory and parallelization management, and logic
computational flows within automated python workflows. Moreover, it guarantees
a seamless interoperability of different software codes that tackle the
different steps usually involved in MBPT simulations, i.e., the preliminary
DFT part (here using Quantum ESPRESSO 2009qe ; 2017qe ), the GW-BSE
calculations (Yambo yambo1 ; yambo2 ), and any required post-processing. In
particular, for the latter, we here introduce a scheme based on maximally
localized Wannier functions marzari_maximally_2012  for the automatic
interpolation of GW band structures, which interfaces the Wannier90Pizzi2020
and Yambo projects. All of these developments point to a drastic reduction of
both human and computational efforts, key issues for enabling HT studies. In
addition, by incorporating the different domain-specific scientific and
computational competences into robust and reliable workflows, we aim at making
accurate GW-BSE calculations available for the materials science community at
large, including non-experts in the field (e.g., via graphical user interfaces
aiidalab ), similarly to what has already happened with DFT.

The manuscript is organised as follows. In the Results and Discussion Section,
we first introduce a model for the convergence surface in the N-dimensional
space of the GW (BSE) variables, and present our improved algorithm for
efficiently retrieving converged values for the main (interdependent)
parameters. The implementation of this algorithm within the aiida-yambo plugin
is then described, followed by the presentation of the novel aiida-yambo-
wannier90 plugin that encodes the GW band interpolation based on
Wannierization. Both these implementations are validated for selected
prototypical systems. Additional details on simulations are provided in the
Methods section. In the remainder of this section, we instead introduce the
main concepts and quantities related to the GW and BSE schemes, which will be
useful to properly understand the subsequent Sections.

### GW approximation

Accurate electronic band structures of materials can be computed within the
MBPT framework by correcting the Kohn-Sham (KS) DFT eigenvalues with a self-
energy term $\Sigma$ by means of the $GW$ approximation, i.e. $\Sigma$ is
approximated with the first term of the perturbation series expansion in terms
of the screened Coulomb interaction $W$ Hedin_1965 . We hereafter consider the
simplest and most widespread implementation of GW, i.e. the so-called
$G_{0}W_{0}$ approach, where $G_{0}$ is the KS independent-particle one-body
Green’s function and $W$ is computed within the random-phase approximation
(RPA) onida_electronic_2002 . Nevertheless, the convergence algorithm
presented here can be used also for more sophisticated flavours of the theory,
e.g. the self-consistent GW. Under these assumptions and by considering a
plane-wave expansion, the self-energy term $\Sigma_{n\textbf{k}}$ for a band
$n$ at a given k-point – written as the sum of the Fock exchange
($\Sigma^{x}$) and the frequency-dependent correlation ($\Sigma^{c}$) terms –
is given by:

$\Sigma^{x}_{n\textbf{k}}=-\sum_{m}^{\text{occ}}\int\frac{d\textbf{q}}{(2\pi)^{3}}\sum_{\textbf{G}}^{G^{x}_{\text{cut}}}v_{\bf{G}}(\textbf{q})|\rho_{nm}(\textbf{k},\textbf{q},\textbf{G})|^{2}f_{m,\textbf{k}-\textbf{q}}$
(1)

and

$\begin{split}\Sigma^{c}_{n\textbf{k}}(\omega)=&-i\sum_{m}^{N_{b}}\int\frac{d\textbf{q}}{(2\pi)^{3}}\sum_{\textbf{G}\textbf{G}^{\prime}}^{G_{\text{cut}}}\rho_{nm}(\textbf{k},\textbf{q},\textbf{G})\rho^{*}_{nm}(\textbf{k},\textbf{q},\textbf{G}^{\prime})\times\\\
&\int
d\omega^{\prime}W_{\bf{GG^{\prime}}}(\textbf{q},\omega^{\prime})\times\\\
&\Bigg{[}\frac{f_{m,\textbf{k}-\textbf{q}}}{\omega-\omega^{\prime}-\epsilon_{m,\textbf{k}-\textbf{q}}-i\eta}+\frac{1-f_{m,\textbf{k}-\textbf{q}}}{\omega-\omega^{\prime}-\epsilon_{m,\textbf{k}-\textbf{q}}+i\eta}\Bigg{]},\end{split}$
(2)

where $v_{\bf{G}}(\textbf{q})$ and
$W_{\bf{GG^{\prime}}}(\textbf{q},\omega^{\prime})$ are the bare and screened
Coulomb interaction, respectively. The generalized dipole matrix elements
$\rho_{nm}(\textbf{k},\textbf{q},\textbf{G})$ are defined as:

$\rho_{nm}(\textbf{k},\textbf{q},\textbf{G})=\bra{n\textbf{k}}e^{i(\textbf{q}+\textbf{G})\cdot\textbf{r}}\ket{m\textbf{k}-\textbf{q}},$
(3)

where $\epsilon^{\text{KS}}_{n\textbf{k}}$ and $|n\mathbf{k}\rangle$ are the
corresponding KS eigenvalues and eigenvectors. Here, $W$ can be written in
terms of the reducible polarizability $\chi$ ($W=v+v\chi v$), which is in turn
computed by solving a Dyson equation for the RPA irreducible polarizability
$\chi^{0}$:

$\chi_{\bf{G}\bf{G}^{\prime}}(\textbf{q},\omega)=\sum^{G_{\text{cut}}}_{\bf{G}^{\prime\prime}}\left[I-v(\bf{q})\chi^{0}(\bf{q},\omega)\right]_{\bf{G}\bf{G}^{\prime\prime}}^{-1}\chi^{0}_{\bf{G}^{\prime\prime}\bf{G}^{\prime}}(\bf{q},\omega),$
(4)

where

$\chi^{0}_{\textbf{G}\textbf{G}^{\prime}}(\textbf{q},\omega)=2\sum_{nm}^{N_{b}}\int_{BZ}\frac{d\textbf{k}}{(2\pi)^{3}}\text{
}\rho^{*}_{mn}(\textbf{k},\textbf{q},\textbf{G})\rho_{mn}(\textbf{k},\textbf{q},\textbf{G}^{\prime})\left[\frac{f_{n\textbf{k}-\textbf{q}}(1-f_{m\textbf{k}})}{\omega+\epsilon_{n\textbf{k}-\textbf{q}}-\epsilon_{m\textbf{k}}}-\frac{f_{n\textbf{k}-\textbf{q}}(1-f_{m\textbf{k}})}{\omega+\epsilon_{m\textbf{k}}-\epsilon_{n\textbf{k}-\textbf{q}}}\right].$
(5)

In practice, the above quantities are computed by introducing specific
approximations that crucially impact both the accuracy and the computational
and memory costs of the calculations. In particular, Eq. (4) as well as Eq.
(2) involve a discrete, finite number of G vectors, set by a cutoff parameter
$G_{\text{cut}}$, which defines the size of the response matrix. Moreover,
with the parameter $N_{b}$, we introduce a truncation over the KS states
summation for both response $\chi^{0}_{\textbf{G}\textbf{G}^{\prime}}$ (Eq.
(5)) and self-energy $\Sigma^{c}$ (Eq. (2)), which should in principle include
all occupied and an infinite number of empty states. Finally, all the
integrals in reciprocal space are computed on a discrete k-points (q-points)
grid whose size, $N_{\textbf{k}}$, defines the accuracy of the sampling of the
Brillouin zone (BZ). All of these three parameters, $G_{\text{cut}}$, $N_{b}$,
and $N_{\textbf{k}}$, need to be increased till the desired convergence is
reached.

One of the major obstacles to automate the convergence procedure lies in the
interdependence of the first two parameters, $G_{\text{cut}}$ and $N_{b}$,
such that their convergence has to be performed jointly, as thoroughly
discussed elsewhere stankovski_g_2011 ; gao_speeding_2016 ;
van_setten_automation_2017 . Indeed, Eqs. 4 and 5 contains a summation over
both empty states and reciprocal lattice vectors (G), and the expression of
the generalized dipole terms in Eq. (3) is such that matrix elements with
large G are governed by high-energy KS states van_setten_automation_2017 .
Furthermore, the interdependence is non-trivial, given the presence of an
inversion in Eq. (4) that further enters in the evaluation of the correlation
self-energy, Eq. (2). Given the lack of an efficient “recipe” to carry out
this non-trivial, coupled convergence, one has to resort to an iterative
convergence of each of the two parameters by fixing the other, in an
alternating way. This procedure, in addition of being tedious, is
computationally very expensive, sometimes representing the most cumbersome
part of a GW calculation.

### Bethe-Salpeter Equation

Starting from the GW quasi-particles, the solution of the BSE
strinati1988application ; onida_electronic_2002  can give access to optical
properties of materials via the macroscopic dielectric function:

$\begin{split}&\epsilon_{M}(\omega,\mathbf{q})=1-\frac{2}{VN_{q}}v(\mathbf{q})\sum_{\lambda}\Phi_{\lambda}(\mathbf{q})\\\
&\times\left[\frac{1}{E^{\lambda}(\mathbf{q})-(\omega+i\eta)}+\frac{1}{E^{\lambda}(\mathbf{q})+(\omega+i\eta)}\right],\end{split}$
(6)

where $V$ is the volume of the unit cell; $N_{q}$ is the number of q-points
sampling the Brillouin zone (BZ); $E^{\lambda}(\mathbf{q})$ is the eigenvalue
of the exciton $\lambda$ at momentum
$\textbf{q}=\textbf{k}-\textbf{k}^{\prime}$, and the corresponding exciton
oscillator strength $\Phi_{\lambda}(\textbf{q})$ is defined as:

$\Phi_{\lambda}(\textbf{q})=\left|\sum_{(vc,\textbf{k})}\bra{v\textbf{k}-\textbf{q}}e^{-i\textbf{q}\cdot\textbf{r}}\ket{c\textbf{k}}A^{\lambda}_{vc\textbf{k}}(\textbf{q})\right|^{2}.$
(7)

Equation (7) contains a summation over valence and conduction bands ($v,c$),
and k-point mesh, which are the main parameters to converge for a BSE
calculation. The terms $A^{\lambda}_{vc\textbf{k}}(\textbf{q})$ represent the
weight of each electron-hole transition contributing the exciton $\lambda$, as
resulting from the solution of the BSE. The summation over bands mainly
defines the range of energies under investigation. The k-point mesh is
connected to the accuracy with which we describe the exciton composition in
terms of single-particle $v\rightarrow c$ transitions over the entire BZ, and
it is usually significantly larger than the one needed to converge the
corresponding GW band structure. Additional convergence parameters are the
number of G vectors used for expanding the KS wave-functions in the transition
matrix elements, as defined in Eq. (3), and for Fast-Fourier-Transform (FFT)
operations, as well as the plane-wave expansion of the BSE kernel (both direct
and exchange terms). These parameters are usually inherited from GW
convergences.

##  II Results and discussion

### Description of the convergence surface

The above-described coupled convergence of the parameters, combined with a
much worse scaling than DFT and more computationally and memory demanding
calculations, call for efficient procedures to describe and explore the
convergence space in GW-BSE simulations. A possible strategy is to describe
the convergence space in terms of an analytic function of the parameters
schindlmayr_analytic_2013 ; klimes_predictive_2014 ; maggio+gwsmall ;
van_setten_automation_2017 ; rangel_repr_2020 ; rasmussen_towards_2021 . For a
general (N+1) dimensional space, a model convergence surface $f(\textbf{x})$
that represents the value of a given observable $E($x$)$ (e.g., quasiparticle
energies or excitonic eigenvalues) as a function of the $N$ parameters
$\textbf{x}=[x_{1},...,x_{N}]$ can be defined as:

$f(\textbf{x})=\prod_{i}^{N}\left(\frac{A_{i}}{x_{i}^{\alpha_{i}}}+b_{i}\right),$
(8)

where $A_{i}$, $b_{i}$ and $\alpha_{i}$ are free fitting parameters, and
$B=\prod_{i}^{N}b_{i}$ is the extrapolated converged value. Of course, the
accuracy of the latter depends on the actual region of the parameter space
explored for the evaluation of the convergence behaviour, i.e., how far is the
extrapolated value from the exact one. As such, $B$ might not always be a good
choice for guiding the search of the convergence parameters.

In the following, we introduce conditions on the mixed partial derivatives of
Eq. (8) as a way to address the parameter interdependence. Indeed, the
adoption of an analytical form for the description of the convergence space
has the clear advantage of enabling the calculation of all-order derivatives,
once the fitting parameters are known. On the other hand, not taking directly
into account this interdependence can result in a very tedious convergence
procedure, as it would require the cyclic repetition of multiple univariate
convergences, as mentioned in the Introduction (further details are provided
in Supplemental Material supp-info and Ref. yambo2, ).

For the expression in Eq. (8), the gradient components are:

$f_{x_{i}}^{\prime}(\textbf{x})=-\alpha_{i}\frac{A_{i}}{x_{i}^{\alpha_{i}+1}}\prod_{j\neq
i}^{N}\left(\frac{A_{j}}{x_{j}^{\alpha_{j}}}+b_{j}\right),$ (9)

while second derivatives are:

$f_{x_{i}}^{\prime\prime}(\textbf{x})=\alpha_{i}(\alpha_{i}+1)\frac{A_{i}}{x_{i}^{\alpha_{i}+2}}\prod_{j\neq
i}^{N}\left(\frac{A_{j}}{x_{j}^{\alpha_{j}}}+b_{j}\right),$ (10)
$f_{x_{i},x_{j}}^{\prime\prime}(\textbf{x})=\alpha_{i}\alpha_{j}\frac{A_{i}A_{j}}{x_{i}^{\alpha_{i}+1}x_{j}^{\alpha_{j}+1}}\prod_{k\neq
i,j}^{N}\left(\frac{A_{k}}{x_{k}^{\alpha_{k}}}+b_{k}\right).$ (11)

The asymptotic region of the convergence surface can be determined by
imposing, for each parameter $x_{i},x_{j}$ with $i,j=1,...,N$, two conditions:

$\left\\{\begin{array}[]{ll}|f_{x_{i}}^{\prime}(\textbf{x})|<\Delta_{i}\\\\[7.0pt]
|f_{x_{i},x_{j}}^{\prime\prime}(\textbf{x})|<\Delta_{ij}\end{array}\right.$
(12)

The first condition determines the region in which the convergence surface
becomes flat (thus approaching convergence), whereas the condition on second
partial derivatives ensures that the $N$ parameters are no longer
interdependent. The threshold values $\Delta_{i}$ and $\Delta_{ij}$ can be
tuned according to the desired accuracy. 111According to our experience, the
asymptotic region can be safely identified by choosing $\Delta_{i}=5\cdot
10^{-5}$ and $\Delta_{ij}=1\cdot 10^{-8}$. Once this asymptotic region has
been determined, a guess for the converged value of $f(\textbf{x})$,
$E_{guess}$, is made, and its accuracy is then checked and validated according
to the automated algorithm described below.

### Convergence algorithm

Figure 1:  Flowchart of the convergence algorithm. After generating the grid
for the N-dimensional parameter space, a subset of M simulations is performed.
The results are then fitted to predict the converged parameters. Finally, the
accuracy (Eq. (13)a) and the convergence (Eq. (13)b) of the prediction are
verified, and the procedure iterated, if needed.

Figure 1 schematically depicts our convergence algorithm, which is purposed to
obtain the desired accuracy on GW-BSE results with the least possible number
of calculations. In the following, we consider $N_{b}$ and $G_{\text{cut}}$ as
the two interdependent parameters to be converged simultaneously. We remark
that, while the algorithm is specifically designed to handle coupled
convergences, it can be successfully used to accelerate convergence tests with
respect to any other parameter, such as the BZ k-point mesh or the FFT grids.

The first step (i) consists in the construction of the N-dimensional space of
parameters as a grid of equally spaced points, with spacing and ranges
provided from input. It is worth noting that, to fit the functional form of
Eq. (8), one needs to generate a grid with minimum 3$N$ points, since
$f(\textbf{x})$ contains three free fitting parameters for each of the $N$
dimensions of the parameter space. However, a further reduction on the minimal
grid size (that is, on the minimum number of calculations to perform) can be
obtained by fixing the power-law dependence, $\alpha_{i}$, to a given value,
as suggested in Ref. 38, thus resulting in a minimal grid size of 2N. Usually,
much denser grids are generated, where (ii) M${}_{0}\geq$ 2$N$ calculations
are performed on a subset of the grid points, chosen such as to efficiently
sample the parameter space.

Next (iii), the results of the calculations are fitted by using the expression
Eq. (8). As mentioned above, the power laws are fixed to given values:
$\alpha_{i}\in$ $\\{$1,2$\\}$ $\forall$ i=1, …N, and the one resulting in the
lowest mean squared error is chosen. The asymptotic region is then identified
by computing the first and second order derivatives (Eqs. (9)-(11)), and
imposing the conditions in Eq. (12). Given the asymptotic region, a guess
converged value is selected $E_{guess}$= $f(N^{0}_{b},G^{0}_{cut})$, where
($N^{0}_{b},G^{0}_{cut}$) are the lowest values that can be chosen for the
parameters such that $E_{guess}$ is within a desired convergence threshold
$\Delta$ with respect to the asymptotic region.

To establish the accuracy and convergence of the fitting procedure, (iv)
$E(N^{0}_{b},G^{0}_{cut})$ is computed and compared to the outcome of the fit
$E_{guess}$, by considering the chosen threshold $\Delta$ (see Eq. (13)a
below). If the accuracy condition is satisfied, we need to check the
convergence of the fitting procedure, that is, a new pair of parameters
$(N^{1}_{b},G^{1}_{cut})$ is obtained from the fit by adding the
$(N^{0}_{b},G^{0}_{cut})$ point to the initial M0 grid, and compared to the
previous one (see Eq. (13)b below). The last step is repeated until
convergence is reached. If the accuracy condition is not satisfied, the grid
is instead shifted toward higher values of the parameters, and the steps
(ii)-(iv) are repeated until the two conditions:

$\displaystyle|E(N^{j}_{b},G^{j}_{cut})-E^{fit}(N^{j}_{b},G^{j}_{cut})|<\Delta$
(13a) $\displaystyle(N^{j+1}_{b},G^{j+1}_{cut})=(N^{j}_{b},G^{j}_{cut}),$
(13b)

are simultaneously satisfied for the $j$-th iteration.

### The aiida-yambo plugin and automated workflows

The above convergence algorithm has been implemented in the new version of the
aiida-yambo plugin aiida-yambo-github , which is meant to fully automate GW-
BSE calculations by interfacing the Yambo project yambo1 ; yambo2  and the
AiiDA informatics infrastructure and workflow management system
huber_aiida_2020 ; uhrin_workflows_2021 . The automation concerns input
generation, scheduler submission, and output parsing phases. 222The output
parsing of the aiida-yambo plugin is partially done by using yambopy functions
yambo2 . Thanks to the AiiDA infrastructure, links between single calculations
are managed on the fly by ad-hoc, dynamic workflows (the so-called workchains
in the AiiDA jargon), i.e. their execution path is not fixed, but can depend
on the results of completed calculations. This allows for the implementation
of complex logics, such as those characterizing the convergence algorithm and
GW band interpolation that we propose in this work. Moreover, each
calculation, together with inputs and outputs, is stored in the AiiDA
relational database, thus ensuring data provenance and full reproducibility of
results.

Currently, the aiida-yambo plugin supports quasiparticle ($G_{0}W_{0}$ and
COSHEX Hedin_1965  level) and optical properties (IP-RPA and BSE) simulations,
as well as interfaces with different codes (e.g., Quantum ESPRESSO and
Wannier90). These options are implemented in the YamboCalculation and
YppCalculation classes, which manage individual simulations (including data
interfacing) that can be performed by using the Yambo code. On top of them,
task-specific workflows are implemented, and organized in a modular way, in
order to automate tasks of increasing complexity. In particular, the aiida-
yambo plugin contains three main workflows, each of them targeting a precise
task:

  * •

YamboRestart: automation of error handling and restart for each
YamboCalculation;

  * •

YamboWorkflow: automation of the single GW or BSE flow (composed of several
interlinked steps, explained in the following);

  * •

YamboConvergence: automation of the convergence (composed of multiple
YamboWorkflow runs).

Their nested organization is shown in Fig. 2.

Figure 2: Hierarchical structure of the aiida-yambo workchains. The highest
level workflow is YamboConvergence, which calls multiple YamboWorkflow
workchains. YamboWorkflow comprises all the steps needed to perform individual
GW-BSE calculations from scratch. In case of failures, it calls the
YamboRestart workchain for automatic error handling. The outputs are stored in
the AiiDA database in a human readable fashion, and are easily accessible and
shareable by the user.

The highest level workflow is represented by the YamboConvergence workchain,
which implements the full automation of the convergence algorithm of Sec. II,
thus allowing for all Yambo simulations to be organized on the fly, without
any external user intervention. The user is only requested to provide a python
list containing the information on the parameter space to be explored. An
example of such input reads:

⬇

[

{

'var': [

'BndsRnXp',

'GbndRnge',

'NGsBlkXp'

],

'start': [50, 50, 2],

'stop': [400, 400, 10],

'delta': [50, 50, 2],

'max': [1000, 1000, 36],

'what': ['gap_GG'],

'conv_thr': 0.1,

'conv_thr_units': 'eV',

},

]

where the Yambo variables “BndsRnXp” and “GbndRnge” govern the convergence
over empty states, $N_{b}$, to be carried out jointly with that on the size of
the response matrix, $G_{\text{cut}}$ (“NGsBlkXp” variable). The edges of the
grid (“start” and “stop”) and its spacing (“delta”), together with an upper
bound of the parameter space (“max”), limiting the search to computationally
accessible calculations, are also set by the user. The key “what” indicates
the quantity to be converged – in our example, the direct band gap of the
material at $\Gamma$ point – up to a given convergence threshold $\Delta$
(“conv$\\_$thr” key). The output summarizes the convergence history and allows
the user to easily parse the converged simulation. YamboConvergence allows one
to converge several many-body quantities, like quasiparticle levels, band-
gaps, as well as optical excitation energies. Notably, the convergence block
in Fig. 2 can be skipped if converged parameters are already known.

Each single GW (BSE) calculation is instead automated within YamboWorkflow,
which is the core workchain of the plugin that takes care of performing all
the steps needed in a typical Yambo simulation – from preliminary self-
consistent (SCF) and non-self-consistent (NSCF) DFT calculations to the actual
GW (BSE) calculations, and the related post-processing. The workflow ensures a
robust interoperability between DFT and MBPT codes (Quantum ESPRESSO and
Yambo, respectively), and links subsequent calculations, interfacing the data
automatically. In practice, YamboWorkflow encodes the specific flowchart
underlying each requested calculation, and allows for its dynamic execution
according to the instructions provided in input. This implies performing all
the intermediate steps needed for a specific calculations without the need of
instructing them explicitly, or, on the contrary, to skip some of the
intermediate steps for which parent calculations are available, fully
exploiting the YamboWorkflow provenance information.

To support a restart mechanism in case of code failures, YamboWorkflow takes
advantage of the YamboRestart workchain, a sub-level workflow that encodes an
automatic error handler (inherited from the AiiDA BaseRestartWorkchain class)
which, depending on the encountered failure, automatically instructs a restart
run. For out-of-memory errors or failures connected with insufficient wall-
time requests, YamboRestart automatically resubmits the calculation by
appropriately changing the requested resources (e.g., the maximum wall-time
and the MPI/OpenMP balance); parallelization errors are managed by overwriting
the parallelism variables set in input by the user with the default
parallelism decided on the fly by yambo. In all these cases, an efficient,
CPU-time saving restart mechanism is implemented, which avoids to restart
unfinished runs from scratch by automatically retrieving and enabling the
reuse of stored data files.

As a final issue, we would like to discuss the possibility to develop
protocols for MBPT calculations. Indeed, most of the DFT-based AiiDA plugins
enable the use of protocols huber_common_2021 , that is, the possibility of
creating inputs with pre-populated default values for several parameters. Such
protocols are usually code-agnostic and robust, given the high level of
reproducibility of DFT with different quantum engines
lejaeghere_reproducibility_2016 . Moreover, their reliability is guaranteed by
means of large scale studies spanning systems with a wide variety of
characteristics (i.e., metals, semiconductors, dimensionality and so on)
prandini_precision_2018 .

Concerning MBPT calculations, the possibility to define protocols is still an
open issue van_setten_automation_2017 . First of all, code-agnostic parameters
are not at all easy to be determined as it is for DFT-based codes, because the
MBPT implementations and the subsequent definition of parameters can differ in
very many aspects, as highlighted in the Introduction section. Secondly, the
high computational cost of these simulations has limited so far the number of
systems to be studied extensively, which is crucial to define a reliable
statistics on convergence parameters. Last but not least, DFT-based protocols
usually result in safe but overconverged parameters, an approach that might
lead to unfeasible calculations when moving to the GW-BSE framework.

For all these reasons, we believe that an efficient, fully automated
convergence tool, as the one presented here, is currently the most valuable
solution. Nonetheless, also in view of possible future developments, the
aiida-yambo plugin provides an implementation for a protocol framework for
both GW and BSE simulations, which is currently pre-populated on the basis of
previous experience on a limited subset of systems. Such protocols concern
several parameters connected to the main yambo input variables, such as the
FFT grids (“FFTGvecs”), the summation over empty states (“BndsRnXp” and
“GbndRnge”), the plane-wave expansion for the polarizability (“NGsBlkXp”) and
the BZ k-point sampling. A detailed documentation on these and other aspects
concerning the aiida-yambo plugin is provided elsewhere aiida-yambo-docs .

### Automatic GW bands interpolation: the aiida-yambo-wannier90 plugin

Figure 3: Flowchart of the YamboWannier90WorkChain for the Wannier
interpolation of GW band structures. After performing the Yambo convergence,
the workflow searches for a commensurate k-point meshes for Yambo and
Wannier90, and carries out the corresponding Yambo QP calculation. Given the
QP corrections, the workflow proceeds with the Wannierization and the band
interpolation at DFT and GW level.

GW band interpolation from Wannierization is a crucial task in order to obtain
the most accurate quasiparticle band structure with the lowest computational
cost. This task is encoded in the aiida-yambo-wannier90 plugin aiida-yambo-
wannier90-github ; aiida-yambo-wannier90-docs . Essentially, the plugin
provides a meta-workflow, called YamboWannier90WorkChain, which utilizes the
automation and error handling of the underlying aiida-yambo and aiida-
wannier90-workflows plugins for GW convergence and Wannierization,
respectively. The flowchart of the workflow is summarized in Fig. 3. Starting
from a given crystal structure, the workflow first launches a YamboConvergence
workflow for automatic convergence. Then, it finds the minimal commensurate
mesh with the Wannier90 ones that satisfies the GW convergence conditions (see
below). Thirdly, it runs YamboWorkflow to compute all the quasiparticle
corrections required for the Wannierization on the commensurate mesh, and a
subsequent ypp calculation (by means of YppRestart) to extract the GW
corrections in a Wannier90 eig file format. Fourthly, the workflow Wannierizes
the KS wavefunctions, saving the unitary transformation matrices of maximal
localization, and interpolates the band structure. Finally, the workflow
performs the Wannierization procedure at $G_{0}W_{0}$ level, which consists in
incorporating the GW corrections into the DFT eigenvalues, and interpolating
the band structure by using the DFT Wannierization outcomes.

A crucial step of the workflow is finding a commensurate mesh for both GW QP
calculations and Wannierization. Indeed, the GW mesh resulting from automated
convergence might not always be compatible with the mesh required by Wannier90
to ensure interpolation accuracy. Notably, considering a Monkhorst-Pack (MP)
grid for the Wannierization, the corresponding GW mesh must be an integer
multiple of the MP grid. We here propose a recipe to find the minimal
commensurate meshes for GW and Wannier90 calculations, as depicted in Fig. 4.
Considering $n_{d}$ as the number of k-points chosen by the YamboConvergence
workflow, and $n_{c}$ the number of k-points chosen by the Wannierization
protocol (typically based on a k-point spacing, 0.2 Å-1) Vitale2020  , the
target is to find a new ($n_{d}^{\prime}$, $n_{c}^{\prime}$) such that the
dense mesh $n_{d}^{\prime}=k\cdot n_{c}^{\prime}$, where $k\in\mathbb{N}$,
i.e., natural number. The given input $(n_{d},n_{c})$ restricts the search
space to a sector bounded by $k_{\text{low}}$ and $k_{\text{high}}$ (see Fig.
4), where $k_{\text{low}}=1$
($n_{d}^{\prime}=n_{c}^{\prime}=\max(n_{d},n_{c})$ is always a good solution),
and $k_{\text{high}}=\lceil\frac{n_{d}}{n_{c}}\rceil$, where
$\lceil\cdot\rceil$ indicates the ceiling integer. The search always succeeds
since $s_{\text{low}}=(\max(n_{d},n_{c}),\max(n_{d},n_{c}))$ and
$s_{\text{high}}=(k_{\text{high}}\cdot n_{c},n_{c})$ are already two good
solutions. In fact, the optimal solution is often inside the triangular region
determined by the input $(n_{d},n_{c})$, $s_{\text{low}}$, and
$s_{\text{high}}$. The final solution is chosen according to the $\ell^{1}$
distance to the input, such to ensure the minimal increase in computational
cost. It is also possible to change the metric, e.g., pushing the solution
towards increasing the Wannier mesh or GW mesh, depending on which calculation
is less cumbersome. The aforementioned recipe is repeated for each of the
three dimensions of the MP grid.

Figure 4: Recipe to find commensurate meshes for GW and Wannier90
calculations. Using input meshes (11, 5) as an example, final commensurate
meshes (12,6) are found. The $k_{high}$ and $k_{low}$ lines, that intersect in
the origin, are respectively the imposed upper and lower bound for searching
the commensurate meshes; the orange dots are possible solutions; the red dot
is the chosen solution, which is the closest to the input in the metric of
$\ell^{1}$ norm.

### Validation of the workflows

System | E${}^{G_{0}W_{0}}_{gap}$ (eV) | E${}^{ref}_{gap}$ (eV) | $N_{b}$ |  $G_{\text{cut}}$ (Ry) |  $\rho_{\textbf{k}}$(Å-1) |  $\Delta^{\Gamma}$(meV) | $\Delta^{\Gamma}_{\%}$
---|---|---|---|---|---|---|---
Si | 1.18 | 1.16 rangel_repr_2020  | 400 | 16 | 0.33 | 10 | 0.3
Diamond | 5.82 | 5.63 gao_2015  | 300 | 20 | 1 | 70 | 0.1
ZnO | 2.36 | 2.35 stankovski_g_2011  | 800 | 28 | 0.25 | 10 | 0.4
TiO2 | 3.20 | 3.20 rangel_repr_2020  | 600 | 12 | 0.2 | 10 | 0.3
MoS2 | 2.54 | 2.54 Rasmussen_Q0  | 400 | 8 | 0.25 | 10 | 0.3
hBN bulk | 6.28 | 6.30 huser_quasiparticle_2013  | 800 | 20 | 0.25 | 82 | 1.0
hBN 2D | 6.84 | 7.06 Rasmussen_Q0  | 1200 | 28 | 0.2 | 42 | 0.5
Table 1: $G_{0}W_{0}$ convergence tests on prototypical semiconductors. For
each system, the minimum band gap $E^{G_{0}W_{0}}_{gap}$ (second column) is
computed by using the parameters obtained by the automated convergence
algorithm implemented in the YamboConvergence workchain, and compared with
previous $GW$ results (third column). The considered parameters are the number
of empty states, $N_{b}$, the energy cutoff on the G vectors,
$G_{\text{cut}}$, and the irreducible Brillouin Zone (iBZ) k-points density
$\rho_{\textbf{k}}$, expressed as the maximum distance between adjacent points
along a reciprocal axis. The first two parameters are converged jointly. The
last two columns include the convergence thresholds imposed on the G0W0 band
gap at the $\Gamma$ point, both in absolute ($\Delta^{\Gamma}$) and relative
($\Delta^{\Gamma}_{\%}$) terms.

The proposed convergence algorithm, as implemented in the YamboConvergence
workchain, has been validated by performing convergence studies for the
quasiparticle $G_{0}W_{0}$ gap of a set of prototypical semiconductors:
silicon, diamond, ZnO, rutile TiO2, monolayer MoS2, bulk and monolayer hBN.
The convergence addresses the direct band gap at the $\Gamma$ point with
respect to the two coupled parameters, $N_{b}$ and $G_{\text{cut}}$, and the
$\mathbf{k}$-point grid (in terms of k-points density $\rho_{\textbf{k}}$).
Once convergence is achieved, the minimum band gap $E^{G_{0}W_{0}}_{gap}$ is
also computed. The results are summarized in Table 1, where, together with
$E^{G_{0}W_{0}}_{gap}$, we report the final parameters resulting from the
convergence procedure ($N_{b}$, $G_{\text{cut}}$, and $\rho_{\textbf{k}}$) as
well as the convergence threshold (absolute, $\Delta^{\Gamma}$, and relative,
$\Delta^{\Gamma}_{\%}$) adopted in each case. Our results are found in good
agreement with previous findings, which are also reported in the Table 1.
Deviations with respect to reference results can be ascribed to different GW
implementations and/or different DFT starting points used in the related works
rangel_repr_2020 . Further details are contained in the Supplemental Material
supp-info .

Figure 5: Convergence algorithm applied to monolayer hBN. (a) Convergence of
the direct quasiparticle band gap at the $\Gamma$ point with respect to the
coupled parameters $N_{b}$ and $G_{\text{cut}}$. The blue shaded area
represents the starting grid, while the red shaded area the shifted grid
obtained after the first iteration, according to the flowchart in Fig. 1.
Black squares represent the actual calculations performed by the workflow; the
blue square is the first guess for the converged parameters; the red square
indicates the final, converged point. The colormap specifies the relative
error with respect to the converged point. A maximum absolute error of
$\Delta$ = 42 meV is achieved for ($N_{b}$, $G_{\text{cut}}$) = (1200,28 Ry),
corresponding to a maximum relative error of $\Delta_{\%}$ = 0.5 $\%$. (b)
Convergence of the direct quasiparticle band gap at the $\Gamma$ point with
respect to the k-mesh. The black points represent the actual calculations
performed by the workflow, whereas the blue points are the ones obtained
within the fitting procedure and used to predict the convergence. The final,
converged mesh (red square) is achieved with five simulations.

Figure 5 shows the convergence procedure for monolayer hBN, considering both
the joint convergence with respect to Nb and $G_{\text{cut}}$ (panel a), and
the single-parameter convergence with respect to the k-mesh (panel b).
Starting from an input grid with N${}_{b}\in$ [200,800] and G${}_{cut}\in$
[4,16] Ry (Figure 5a, blue shaded area), a subset of 6 calculations is
performed (black squares). Since the converged guess is above the upper bound
of the parameter space (see “max” variable in the plugin description), a new
shifted grid (orange shaded area) is considered, and a first guess for the
converged parameters is found from the fitting procedure (blue square), now
satisfying the accuracy condition (Eq. (13)a). A new fit is performed
including this additional point, which results in a new converged guess (red
square). The procedure is repeated till the converged result is verified to be
consistent with the prediction within the given threshold (Eq. (13)a) and to
be the true converged point (Eq. (13)b). A similar path is followed for the
k-mesh convergence (Figure 5b): a limited number of calculations is initially
performed (black squares), from which the fitting is evaluated (blue curve),
and the smallest grid compatible with the given threshold is finally selected
(here 8$\times$8$\times$1). We note that, despite the results of the
simulations seem to have an oscillating behaviour with respect to the fitted
curve, the error bar considered here (from the 14$\times$14$\times$1 mesh) is
$\sim 0.13\%$ of the band gap at $\Gamma$, i.e. $\sim 10$ meV, in line with
the accuracy of state-of-the-art GW calculations.

The dependence of the outcome of the algorithm on the input settings (e.g.,
the initial value of the parameters and the boundaries of the convergence
space) is a key issue for evaluating the robustness of the algorithm itself.
Indeed, the convergence procedure has been tested on the starting grid for
monolayer MoS2. By using two different grids, N${}_{b}\in$ [200,800],
G${}_{cut}\in$ [4,20] Ry and N${}_{b}\in$ [200,1200], G${}_{cut}\in$ [8,24]
Ry, the convergence point obtained from the workflow remains the same, i.e.,
(Nb,$G_{\text{cut}}$) = (400,8). Another important issue to evaluate is the
efficiency of the algorithm. We notice that, in the case of the 2D-hBN
convergence shown in the left panel of Figure 5, only 14 calculations where
required to reach convergence. Older implementations yambo2  would have
required $\geq$ 25 simulations to achieve a final result, with over 40%
reduction.

Next, the YamboWannier90WorkChain as been tested on bulk Si and Cu, in order
to validate the automatic Wannier interpolation of GW band structures for both
semiconductors and metals. Results are plotted in Fig. 6, where we compare the
DFT bands, the Wannier interpolated DFT, and the Wannier interpolated GW
bands. For Si (Fig. 6a), the comparison between computed and interpolated DFT
bands shows that the results are almost identical, indicating the accuracy of
the Wannierization of the KS wavefunctions. Moreover, the typical band gap
opening upon inclusion of GW corrections is found when comparing KS and QP
band structures.

For Cu (Fig. 6b), we obtain a discrepancy of $\sim$ 10 meV around the Fermi
energy (here set to zero) between computed and interpolated bands at the DFT
level. Better accuracy can be achieved imposing more stringent values of the
involved parameters. At QP level, the GW correction is very small around the
Fermi level ($\sim$ 37 meV), but still not negligible. Here, the GW
convergence is more stringent than for Si, especially concerning the k-point
mesh. Indeed, denser grids are needed to account for the contribution of
intra-band transitions in the q$\rightarrow$0 limit, which is crucial for
metallic systems but not explicitly included within the plasmon pole
approximation hybertsen_electron_1986 ; godby_metal-insulator_1989 .
Considering the converged parameters, ($N_{b}$, $G_{\text{cut}}$,
$\rho_{\textbf{k}}$) = (400, 18 Ry, 0.2 Å-1), the quasiparticle evaluations
required to interpolate the bands for the minimum converged Wannier90 k-point
mesh (16$\times$16$\times$16) become 2900. This quite large number of QP can
be easily computed using the YamboWorkflow workchain thanks to the possibility
to split the QP calculation in several runs, each of them computing only a
fraction of the GW corrections, and then collecting all the data in a final
database well-suited for the YamboWannier90WorkChain. Since the number of QP
corrections to compute can be quite high, in the Supporting Information
(section S.IIIA) we suggest an effective way to reduce the number of the
required calculations for cases when accurate QP corrections are needed only
in a limited energy region, e.g. around the Fermi energy, and energies outside
of the chosen region can be approximated e.g. through a scissor and stretching
correction.

Figure 6: Wannier interpolation of GW band structures. Band structure of Si
(a) and Copper (b). Interpolated bands are plotted for both DFT (red dashed
line) and GW (green solid line) eigenvalues, as compared to the DFT computed
bands (black dots).

##  III Conclusions

In this work, we have presented the successful design and implementation of
advanced algorithms in state-of-the-art GW (BSE) calculations, that is,
convergence between interdependent parameters, error handling and automatic
band interpolation by means of Wannierization. We validated the tools on
selected cases among semiconductors and metallic systems. The results
contained in this work clearly show the power of these newly developed
workflows for the automated study of excited states properties of materials,
paving the way for achieving high-throughput MBPT studies. Thanks to these
developments and within the next-generation of pre-exascale and exascale
supercomputers, these simulations may become extensively and routinely
performed by the materials-science community in the near future.

##  IV Methods

For all the systems studied here, we used symmetrized geometries in such a way
to reduce the computational cost of simulations. We do not expect relevant
differences in the results obtained with fully-relaxed structures. DFT
simulations were carried out by using the Quantum ESPRESSO simulation package,
which implements plane-wave basis set and pseudopotential approach. The KS-DFT
exchange-correlation functional was approximated using GGA-PBE Perdew-Burke-
Ernzerhof1996PRL , through the optimized norm-conserving Vanderbilt (ONCV)
SG15 oncv1 ; oncv2  pseudopotentials. In the case of ZnO, we adopted Local
Density Approximation (LDA), to compare the results with the existing
literature stankovski_g_2011 ; rangel_repr_2020 , and PseudoDojo
pseudopotentials van_setten_pseudodojo_2018 . GW and BSE results were obtained
by means of the Yambo code. The frequency dependence of the screened
interaction potential was approximated by using the Godby-Needs plasmon pole
approximation godby_metal-insulator_1989  (GNPPA), and the quasiparticle
energies were calculated according to the G0W0 approximation Hedin_1965 ;
onida_electronic_2002 , as implemented in Yambo. The Bruneval-Gonze technique
BG  was used to reduce the number of empty states $N_{b}$ needed for the
construction of the correlation Self-Energy $\Sigma_{c}$ (Eq. (2)). For low-
dimensional systems, spurious interactions between supercell replica were
avoided using a slab truncation of the Coulomb potential rozzi_exact_2006
along the non-periodic direction; its divergences are cured by means of the
Random Integration Method pulci_ab_1998  (RIM), which also accelerates
convergence with respect to the BZ sampling. For 2D systems, specifically, we
adopted a recently developed accelerating technique based on stochastic
integration of the screened potential guand_RIMW , which allows to have GW-
converged results using reduced Monkhorst-Pack k-points grids, just slightly
denser than the DFT one. Finally, Wannierization and band interpolations are
performed by means of the Wannier90 code. All the simulations are performed
using the automated workflows implemented in the aiida-yambo and aiida-yambo-
wannier90 plugins, developed for the AiiDA platform and presented here as part
of the results achieved in this work. Input parameters are generated using the
protocols procedure, as implemented in the corresponding plugins.

## Data availability

The data supporting the findings of this paper are available on the Materials
Cloud talirz_materials_2020  at https://doi.org/10.24435/materialscloud:6w-qh.
Results obtained in this work can be reproduced by means of the example
scripts delivered within the aiida-yambo and aiida-yambo-wannier90 plugins.

## Code availability

All the codes used in this work are fully available to the community by means
of their repositories, and supported by appropriate documentations. The Yambo
code is accessible at https://www.yambo-code.eu/download/. The Quantum
ESPRESSO and Wannier90 codes can be found, respectively, at
https://www.quantum-espresso.org/download and http://www.wannier.org/download.

The AiiDA infrastructure is available at http://www. aiida.net/download. AiiDA
plugins can be downloaded from the corresponding GitHub repositories, already
referenced throughout the manuscript, see Refs. 61; 64.

## Acknowledgements

We acknowledge stimulating discussions with Nicola Marzari, Michael O. Atambo,
Gianluca Prandini, Dario A. L. Valido, Alberto Guandalini and Fulvio Paleari.
This work was supported by: MaX – MAterials design at the eXascale – a
European Centre of Excellence funded by the European Union’s program
H2020-INFRAEDI-2018-1 (Grant No. 824143); the European Union’s Horizon 2020
research and innovation programme (Grant No. 957189, BIG-MAP, also part of the
BATTERY 2030+ initiative, Grant No. 957213); SUPER (Supercomputing Unified
Platform – Emilia-Romagna) from Emilia-Romagna PORFESR 2014-2020 regional
funds; the Swiss National Science Foundation (SNSF) Project Funding (Grant No.
200021E$\\_$206190 “FISH4DIET”); NCCR MARVEL, a National Centre of Competence
in Research, funded by the Swiss National Science Foundation (Grant No.
205602). Computational time on the Marconi100 and Galileo100 machines at
CINECA was provided by the Italian ISCRA program.

## Author contributions

M.B., J.Q., A.M. and N.S. contributed to the development of the aiida-yambo
plugin; M.B. designed, implemented and tested the automatic convergence
algorithm, and the other workflows belonging to the aiida-yambo plugin. J.Q.
and M.B. implemented the aiida-yambo-wannier90 plugin. E.M., A.F., D.V., D.P,
G.P. were responsible for the project supervision and coordination. M.B.,
J.Q., and D.P. wrote the manuscript with contributions from all authors.

The Authors declare no Competing Financial or Non-Financial Interests

## References

  * (1) S. Curtarolo, G. L. W. Hart, M. B. Nardelli, N. Mingo, S. Sanvito, and O. Levy, The high-throughput highway to computational materials design,  Nature Materials 12, 191 (2013).
  * (2) K. S. Vecchio, O. F. Dippo, K. R. Kaufmann, and X. Liu, High-throughput rapid experimental alloy development (HT-READ),  Acta Materialia 221, 117352 (2021).
  * (3) S. Luo, T. Li, X. Wang, M. Faizan, and L. Zhang, High‐throughput computational materials screening and discovery of optoelectronic semiconductors,  WIREs Computational Molecular Science 11 (2021).
  * (4) M. Ashton, J. Paul, S. B. Sinnott, and R. G. Hennig, Topology-Scaling Identification of Layered Solids and Stable Exfoliated 2D Materials,  Physical Review Letters 118, 106101 (2017).
  * (5) G. Cheon, K.-A. N. Duerloo, A. D. Sendek, C. Porter, Y. Chen, and E. J. Reed, Data Mining for New Two- and One-Dimensional Weakly Bonded Solids and Lattice-Commensurate Heterostructures,  Nano Letters 17, 1915 (2017).
  * (6) N. Mounet, M. Gibertini, P. Schwaller, D. Campi, A. Merkys, A. Marrazzo, T. Sohier, I. E. Castelli, A. Cepellotti, G. Pizzi, and et al., Two-dimensional materials from high-throughput computational exfoliation of experimentally known compounds,  Nature Nanotechnology 13, 246–252 (2018).
  * (7) K. Choudhary, I. Kalish, R. Beams, and F. Tavazza, High-throughput Identification and Characterization of Two-dimensional Materials using Density functional theory,  Scientific Reports 7, 5179 (2017).
  * (8) N. Mounet, M. Gibertini, P. Schwaller, D. Campi, A. Merkys, A. Marrazzo, T. Sohier, I. E. Castelli, A. Cepellotti, and G. Pizzi, Two-dimensional materials from high-throughput computational exfoliation of experimentally known compounds,  Materials Cloud Archive 2020.158 (2020).
  * (9) S. Haastrup, M. Strange, M. Pandey, T. Deilmann, P. S. Schmidt, N. F. Hinsche, M. N. Gjerding, D. Torelli, P. M. Larsen, A. C. Riis-Jensen, J. Gath, K. W. Jacobsen, J. Jørgen Mortensen, T. Olsen, and K. S. Thygesen, The Computational 2D Materials Database: high-throughput modeling and discovery of atomically thin crystals,  2D Materials 5, 042002 (2018).
  * (10) A. Marrazzo, M. Gibertini, D. Campi, N. Mounet, and N. Marzari, Relative Abundance of Z2 Topological Order in Exfoliable Two-Dimensional Insulators,  Nano Letters 19, 8431 (2019).
  * (11) S. Kirklin, B. Meredig, and C. Wolverton, High-Throughput Computational Screening of New Li-Ion Battery Anode Materials,  Advanced Energy Materials 3, 252 (2013).
  * (12) Z. Zhang, X. Zhang, X. Zhao, S. Yao, A. Chen, and Z. Zhou, Computational Screening of Layered Materials for Multivalent Ion Batteries,  ACS Omega 4, 7822 (2019).
  * (13) W. Chen, J.-H. Pöhls, G. Hautier, D. Broberg, S. Bajaj, U. Aydemir, Z. M. Gibbs, H. Zhu, M. Asta, G. J. Snyder, B. Meredig, M. A. White, K. Persson, and A. Jain, Understanding thermoelectric properties from high-throughput calculations: trends, insights, and comparisons with experiment,  Journal of Materials Chemistry C 4, 4414 (2016).
  * (14) S. Bhattacharya and G. K. H. Madsen, High-throughput exploration of alloying as design strategy for thermoelectrics,  Physical Review B 92, 085205 (2015).
  * (15) I. E. Castelli, T. Olsen, S. Datta, D. D. Landis, S. Dahl, K. S. Thygesen, and K. W. Jacobsen, Computational screening of perovskite metal oxides for optimal solar light capture,  Energy Environ. Sci. 5, 5814 (2012).
  * (16) L. Yu and A. Zunger, Identification of Potential Photovoltaic Absorbers Based on First-Principles Spectroscopic Screening of Materials,  Physical Review Letters 108, 068701 (2012).
  * (17) Q. Yan, J. Yu, S. K. Suram, L. Zhou, A. Shinde, P. F. Newhouse, W. Chen, G. Li, K. A. Persson, J. M. Gregoire, and J. B. Neaton, Solar fuels photoanode materials discovery by integrating high-throughput theory and experiment,  Proceedings of the National Academy of Sciences 114, 3040 (2017).
  * (18) K. Kuhar, M. Pandey, K. S. Thygesen, and K. W. Jacobsen, High-Throughput Computational Assessment of Previously Synthesized Semiconductors for Photovoltaic and Photoelectrochemical Devices,  ACS Energy Letters 3, 436 (2018).
  * (19) K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I. E. Castelli, S. J. Clark, A. Dal Corso, S. de Gironcoli, T. Deutsch, J. K. Dewhurst, I. Di Marco, C. Draxl, M. Dułak, O. Eriksson, J. A. Flores-Livas, K. F. Garrity, L. Genovese, P. Giannozzi, M. Giantomassi, S. Goedecker, X. Gonze, O. Grånäs, E. K. U. Gross, A. Gulans, F. Gygi, D. R. Hamann, P. J. Hasnip, N. A. W. Holzwarth, D. Iuşan, D. B. Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik, E. Küçükbenli, Y. O. Kvashnin, I. L. M. Locht, S. Lubeck, M. Marsman, N. Marzari, U. Nitzsche, L. Nordström, T. Ozaki, L. Paulatto, C. J. Pickard, W. Poelmans, M. I. J. Probert, K. Refson, M. Richter, G.-M. Rignanese, S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, S. Sharma, F. Tavazza, P. Thunström, A. Tkatchenko, M. Torrent, D. Vanderbilt, M. J. van Setten, V. Van Speybroeck, J. M. Wills, J. R. Yates, G.-X. Zhang, and S. Cottenier, Reproducibility in density functional theory calculations of solids,  Science 351, aad3000 (2016).
  * (20) G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, Precision and efficiency in solid-state pseudopotential calculations,  npj Computational Materials 4, 72 (2018).
  * (21) S. Maffioletti and R. Murri,  GC3Pie: A Python framework for high-throughput computing,  in Proceedings of EGI Community Forum 2012 / EMI Second Technical Conference — PoS(EGICF12-EMITC2), page 143, Munich, Germany, 2012, Sissa Medialab.
  * (22) S. Curtarolo, W. Setyawan, G. L. Hart, M. Jahnatek, R. V. Chepulskii, R. H. Taylor, S. Wang, J. Xue, K. Yang, O. Levy, M. J. Mehl, H. T. Stokes, D. O. Demchenko, and D. Morgan, AFLOW: An automatic framework for high-throughput materials discovery,  Computational Materials Science 58, 218 (2012).
  * (23) A. Jain, S. P. Ong, W. Chen, B. Medasani, X. Qu, M. Kocher, M. Brafman, G. Petretto, G. Rignanese, G. Hautier, D. Gunter, and K. A. Persson, FireWorks: a dynamic workflow system designed for high‐throughput applications,  Concurrency and Computation: Practice and Experience 27, 5037 (2015).
  * (24) G. Pizzi, A. Cepellotti, R. Sabatini, N. Marzari, and B. Kozinsky, AiiDA: automated interactive infrastructure and database for computational science,  Comput. Mater. Sci. 111, 218 (2016).
  * (25) A. Hjorth Larsen, J. Jørgen Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. Bjerre Jensen, J. Kermode, J. R. Kitchin, E. Leonhard Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. Bergmann Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacobsen, The atomic simulation environment—a Python library for working with atoms,  Journal of Physics: Condensed Matter 29, 273002 (2017).
  * (26) K. Mathew, J. H. Montoya, A. Faghaninia, S. Dwarakanath, M. Aykol, H. Tang, I.-h. Chu, T. Smidt, B. Bocklund, M. Horton, J. Dagdelen, B. Wood, Z.-K. Liu, J. Neaton, S. P. Ong, K. Persson, and A. Jain, Atomate: A high-level interface to generate, execute, and analyze computational materials science workflows,  Computational Materials Science 139, 140 (2017).
  * (27) J. Mortensen, M. Gjerding, and K. Thygesen, MyQueue: Task and workflow scheduling system,  Journal of Open Source Software 5, 1844 (2020).
  * (28) S. P. Huber, S. Zoupanos, M. Uhrin, L. Talirz, L. Kahle, R. Häuselmann, D. Gresch, T. Müller, A. V. Yakutovich, C. W. Andersen, F. F. Ramirez, C. S. Adorf, F. Gargiulo, S. Kumbhar, E. Passaro, C. Johnston, A. Merkys, A. Cepellotti, N. Mounet, N. Marzari, B. Kozinsky, and G. Pizzi, AiiDA 1.0, a scalable computational infrastructure for automated reproducible workflows and data provenance,  Sci. Data 7, 300 (2020).
  * (29) M. Uhrin, S. P. Huber, J. Yu, N. Marzari, and G. Pizzi, Workflows in AiiDA: Engineering a high-throughput, event-based engine for robust and modular computational workflows,  Comp. Mat. Sci. 187, 110086 (2021).
  * (30) A. Bablich, S. Kataria, and M. C. Lemme, Graphene and Two-Dimensional Materials for Optoelectronic Applications,  Electronics 5 (2016).
  * (31) X. Zhang, J. Shao, C. Yan, R. Qin, Z. Lu, H. Geng, T. Xu, and L. Ju, A review on optoelectronic device applications of 2D transition metal carbides and nitrides,  Materials and Design 200, 109452 (2021).
  * (32) J. Jean, P. R. Brown, R. L. Jaffe, T. Buonassisi, and V. Bulović, Pathways for solar photovoltaics,  Energy Environ. Sci. 8, 1200 (2015).
  * (33) S. Zhu and D. Wang, Photocatalysis: Basic Principles, Diverse Forms of Implementations and Emerging Scientific Opportunities,  Adv. En. Mater. 7, 1700841 (2017).
  * (34) W. Jin and L. Hu, Review on Quasi One-Dimensional CdSe Nanomaterials: Synthesis and Application in Photodetectors,  Nanomaterials 9 (2019).
  * (35) F. Xia, T. Mueller, Y.-m. Lin, A. Valdes-Garcia, and P. Avouris, Ultrafast graphene photodetector,  Nature Nanotechnology 4, 839 (2009).
  * (36) W. Lee, S. Ha, H. Lee, J. Bae, B. Jang, H. Kwon, Y. Yun, S. Lee, and J. Jang, High‐Detectivity Flexible Near‐Infrared Photodetector Based on Chalcogenide Ag ${}_{\textrm{2}}$ Se Nanoparticles,  Advanced Optical Materials 7, 1900812 (2019).
  * (37) R. M. Martin, L. Reining, and D. M. Ceperley,  Interacting Electrons: Theory and Computational Approaches,  Cambridge University Press, 2016.
  * (38) M. J. van Setten, M. Giantomassi, X. Gonze, G.-M. Rignanese, and G. Hautier, Automation methodologies and large-scale validation for G W : Towards high-throughput G W calculations,  Physical Review B 96, 155207 (2017).
  * (39) A. Rasmussen, T. Deilmann, and K. S. Thygesen, Towards fully automated GW band structure calculations: What we can learn from 60.000 self-energy evaluations,  npj Computational Materials 7, 22 (2021).
  * (40) F. Hüser, T. Olsen, and K. S. Thygesen, Quasiparticle GW calculations for solids, molecules, and two-dimensional materials,  Physical Review B 87, 235132 (2013).
  * (41) T. Rangel, M. Del Ben, D. Varsano, G. Antonius, F. Bruneval, F. H. da Jornada, M. J. van Setten, O. K. Orhan, D. D. O’Regan, A. Canning, A. Ferretti, A. Marini, G.-M. Rignanese, J. Deslippe, S. G. Louie, and J. B. Neaton, Reproducibility in G 0 W 0 calculations for solids,  Computer Physics Communications 255, 107242 (2020).
  * (42) M. Stankovski, G. Antonius, D. Waroquiers, A. Miglio, H. Dixit, K. Sankaran, M. Giantomassi, X. Gonze, M. Côté, and G.-M. Rignanese, G 0 W 0 band gap of ZnO: Effects of plasmon-pole models,  Physical Review B 84, 241201 (2011).
  * (43) R. Mercado, R.-S. Fu, A. V. Yakutovich, L. Talirz, M. Haranczyk, and B. Smit, In Silico Design of 2D and 3D Covalent Organic Frameworks for Methane Storage Applications,  Chemistry of Materials 30, 5069 (2018).
  * (44) G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, Precision and efficiency in solid-state pseudopotential calculations,  npj Computational Materials 4, 72 (2018).
  * (45) V. Vitale, G. Pizzi, A. Marrazzo, J. R. Yates, N. Marzari, and A. A. Mostofi, Automated high-throughput Wannierisation,  npj Comput. Mater. 6, 66 (2020).
  * (46) P. Giannozzi, S. Baroni, N. Bonini, M. Caland ra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials,  J. Phys.: Condens. Matter 21, 395502 (2009).
  * (47) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, J. DiStasio, R. A., A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H. Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H. V. Nguyen, A. Otero-de-la-Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with Quantum ESPRESSO,  J. Phys.: Condens. Matter 29, 465901 (2017).
  * (48) A. Marini, C. Hogan, M. Grüning, and D. Varsano, yambo: An ab initio tool for excited state calculations,  Comput. Phys. Commun. 180, 1392 (2009).
  * (49) D. Sangalli, A. Ferretti, H. Miranda, C. Attaccalite, I. Marri, E. Cannuccia, P. Melo, M. Marsili, F. Paleari, A. Marrazzo, G. Prandini, P. Bonfà, M. O. Atambo, F. Affinito, M. Palummo, A. Molina-Sánchez, C. Hogan, M. Grüning, D. Varsano, and A. Marini, Many-body perturbation theory calculations using the yambo code,  J. Phys.: Condens. Matter 31, 325902 (2019).
  * (50) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Maximally localized Wannier functions: Theory and applications,  Reviews of Modern Physics 84, 1419 (2012).
  * (51) G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, Wannier90 as a community code: new features and applications,  Journal of Physics: Condensed Matter 32, 165902 (2020).
  * (52) A. V. Yakutovich, K. Eimre, O. Schütt, L. Talirz, C. S. Adorf, C. W. Andersen, E. Ditler, D. Du, D. Passerone, B. Smit, N. Marzari, G. Pizzi, and C. A. Pignedoli, AiiDAlab – an ecosystem for developing, executing, and sharing scientific workflows,  Computational Materials Science 188, 110165 (2021).
  * (53) L. Hedin, New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem,  Phys. Rev. 139, 796 (1965).
  * (54) G. Onida, L. Reining, and A. Rubio, Electronic excitations: density-functional versus many-body Green’s-function approaches,  Rev. Mod. Phys. 74, 601 (2002).
  * (55) W. Gao, W. Xia, X. Gao, and P. Zhang, Speeding up GW Calculations to Meet the Challenge of Large Scale Quasiparticle Predictions,  Scientific Reports 6, 36849 (2016).
  * (56) G. Strinati, Application of the Green’s functions method to the study of the optical properties of semiconductors,  La Rivista del Nuovo Cimento (1978-1999) 11, 1 (1988).
  * (57) A. Schindlmayr, Analytic evaluation of the electronic self-energy in the G W approximation for two electrons on a sphere,  Physical Review B 87, 075104 (2013).
  * (58) J. Klimeš, M. Kaltak, and G. Kresse, Predictive G W calculations using plane waves and pseudopotentials,  Physical Review B 90, 075125 (2014).
  * (59) E. Maggio, P. Liu, M. J. van Setten, and G. Kresse, GW100: A Plane Wave Perspective for Small Molecules,  Journal of Chemical Theory and Computation 13, 635 (2017),  PMID: 28094981.
  * (60) See Supplemental Material for details on the older convergence algorithm, convergence plots obtained within this work and additional information in the Wannierization steps.
  * (61) The aiida-yambo code is available at https://github.com/yambo-code/aiida-yambo.
  * (62) S. P. Huber, E. Bosoni, M. Bercx, J. Bröder, A. Degomme, V. Dikan, K. Eimre, E. Flage-Larsen, A. Garcia, L. Genovese, D. Gresch, C. Johnston, G. Petretto, S. Poncé, G.-M. Rignanese, C. J. Sewell, B. Smit, V. Tseplyaev, M. Uhrin, D. Wortmann, A. V. Yakutovich, A. Zadoks, P. Zarabadi-Poor, B. Zhu, N. Marzari, and G. Pizzi, Common workflows for computing material properties using different quantum engines,  npj Computational Materials 7, 136 (2021).
  * (63) The aiida-yambo documentation is available at https://aiida-yambo.readthedocs.io/en/master/.
  * (64) The aiida-yambo-wannier90 code is available at https://github.com/aiidaplugins/aiida-yambo-wannier90.
  * (65) The aiida-yambo-wannier90 documentation is available at https://aiida-yambo-wannier90.readthedocs.io/en/latest/.
  * (66) S.-P. Gao, Band gaps and dielectric functions of cubic and hexagonal diamond polytypes calculated by many-body perturbation theory,  physica status solidi (b) 252, 235 (2015).
  * (67) F. A. Rasmussen, P. S. Schmidt, K. T. Winther, and K. S. Thygesen, Efficient many-body calculations for two-dimensional materials using exact limits for the screened potential: Band gaps of MoS2$,h$-BN, and phosphorene,  Phys. Rev. B 94, 155406 (2016).
  * (68) M. S. Hybertsen and S. G. Louie, Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies,  Phys. Rev. B 34, 5390 (1986).
  * (69) R. W. Godby and R. J. Needs, Metal-insulator transition in Kohn-Sham theory and quasiparticle theory,  Phys. Rev. Lett. 62, 1169 (1989).
  * (70) J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple,  Phys. Rev. Lett. 77, 3865 (1996).
  * (71) D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials,  Phys. Rev. B 88, 085117 (2013).
  * (72) M. Schlipf and F. Gygi, Optimization algorithm for the generation of ONCV pseudopotentials,  Computer Physics Communications 196, 36 (2015).
  * (73) M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table,  Computer Physics Communications 226, 39 (2018).
  * (74) F. Bruneval and X. Gonze, Accurate G W self-energies in a plane-wave basis using only a few empty states: Towards large systems,  Physical Review B 78, 085125 (2008).
  * (75) C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, Exact Coulomb cutoff technique for supercell calculations,  Phys. Rev. B 73, 205119 (2006).
  * (76) O. Pulci, G. Onida, R. Del Sole, and L. Reining, Ab Initio Calculation of Self-Energy Effects on Optical Properties of GaAs(110),  Physical Review Letters 81, 5374 (1998).
  * (77) A. Guandalini, P. D’Amico, A. Ferretti, and D. Varsano,  Efficient GW calculations in two dimensional materials through a stochastic integration of the screened potential, 2022.
  * (78) L. Talirz, S. Kumbhar, E. Passaro, A. V. Yakutovich, V. Granata, F. Gargiulo, M. Borelli, M. Uhrin, S. P. Huber, S. Zoupanos, C. S. Adorf, C. W. Andersen, O. Schütt, C. A. Pignedoli, D. Passerone, J. VandeVondele, T. C. Schulthess, B. Smit, G. Pizzi, and N. Marzari, Materials Cloud, a platform for open computational science,  Scientific Data 7, 299 (2020).


Rajeev Singh

# Formalism of hydrodynamics with spin degrees of freedom

###### Abstract

In this article we review the perfect-fluid hydrodynamics with spin framework
proposed recently. This framework generalises the standard relativistic
hydrodynamics framework to include spin degrees of freedom and provides a
natural method to describe the spin polarization evolution of massive spin 1/2
particles. This formalism is based on the GLW (de Groot - van Leeuwen - van
Weert) energy-momentum tensor and spin tensor. We show here using Bjorken
model that how this spin hydrodynamics framework may be used for the
determination of the observables which describes the particle polarization
measured in the experiment.

##  1 Introduction

Measurements of $\Lambda$ hyperons spin polarization made by the STAR
Collaboration [1, 2, 4, 3] created immense interest in the theoretical physics
community focusing on the interplay between the orbital angular momentum of
the matter and particle polarization in the relativistic heavy-ion collisions
[5, 6, 7, 8, 9, 10, 27, 12, 20, 21, 22, 23, 26, 24, 19, 25, 11, 13, 14, 15,
18, 17, 35, 16, 28, 29, 30, 31, 32, 33, 34]. Thermal-based models describing
the global spin polarization [38, 39, 40, 41, 36, 37], somehow fail to explain
the differential results [4]. These models deal with the spin polarization at
the freeze-out where they assume that thermal vorticity determines spin
polarization of particles [42, 6]. Following Refs. [43, 44, 45, 46], we
investigate the spacetime evolution of particle spin polarization by including
the spin degrees of freedom in the standard hydrodynamic formalism in the
Bjorken hydrodynamical setup [47].

##  2 Hydro-dynamical equations and spin polarization tensor

For particles with spin-${\nicefrac{{1}}{{2}}}$, perfect fluid hydrodynamics
formalism is formed on the conservation laws for charge, energy and linear
momentum and, angular momentum, which is based on GLW [48] energy-momentum
tensor, $T^{\alpha\beta}_{\rm GLW}$, and spin tensor,
$S^{\alpha\beta\gamma}_{\rm GLW}$111Herein we assume small spin polarization
of particles ($|\omega_{\mu\nu}|<1$)., namely [43, 44, 45]

$\displaystyle\quad\partial_{\mu}N^{\mu}=0,\qquad\partial_{\mu}T^{\mu\nu}_{\rm
GLW}=0,\qquad\partial_{\lambda}S_{\rm GLW}^{\lambda,\alpha\beta}=T_{\rm
GLW}^{\beta\alpha}-T_{\rm GLW}^{\alpha\beta},$ (1)

with

$\displaystyle N^{\alpha}=nU^{\alpha},\qquad T^{\alpha\beta}_{\rm
GLW}=(\varepsilon+P)U^{\alpha}U^{\beta}-Pg^{\alpha\beta},$ (2)

with $N^{\alpha}$ being the net baryon charge current and $\varepsilon$, $P$,
$n$ is the energy density, pressure and baryon number density, respectively,
and $U^{\beta}$ is the time-like hydrodynamic flow four-vector. Due to the
symmetricity of energy-momentum tensor in Eq. (1), conservation of the angular
momentum in-turn implies conservation of the spin part separately [46]. The
spin tensor [44] is $S^{\alpha,\beta\gamma}_{\rm GLW}={\cal
C}\left[n_{(0)}(T)U^{\alpha}\omega^{\beta\gamma}+{\cal
A}_{(0)}\,U^{\alpha}U^{\delta}U^{[\beta}\omega^{\gamma]}_{\hphantom{\alpha}\delta}+\,{\cal
B}_{(0)}\,\Big{(}U^{[\beta}\Delta^{\alpha\delta}\omega^{\gamma]}_{\hphantom{\alpha}\delta}+U^{\alpha}\Delta^{\delta[\beta}\omega^{\gamma]}_{\hphantom{\alpha}\delta}+U^{\delta}\Delta^{\alpha[\beta}\omega^{\gamma]}_{\hphantom{\alpha}\delta}\Big{)}\right]$
where ${\cal C}=\cosh(\xi)$, ${\cal B}_{(0)}=-\frac{2}{\hat{m}^{2}}s_{(0)}(T)$
and ${\cal A}_{(0)}=-3{\cal B}_{(0)}+2n_{(0)}(T)$, and $n_{(0)}(T)$,
$s_{(0)}(T)$ is the number density and entropy density of spin-less and
neutral massive Boltzmann particles [49], respectively, with
$\Delta^{\alpha\beta}$ being the spatial projection operator orthogonal to $U$
with $\xi=\mu/T$ (i.e. baryon chemical potential over temperature) and
$\hat{m}=m/T$, $m$ is the particle mass.
Being an asymmetric second rank tensor, the polarization tensor
$\omega_{\alpha\beta}$ can be written as

$\displaystyle\omega_{\alpha\beta}$ $\displaystyle=$
$\displaystyle\kappa_{\alpha}U_{\beta}-\kappa_{\beta}U_{\alpha}+\epsilon_{\alpha\beta\gamma\delta}U^{\gamma}\omega^{\delta},$
(3)

where $\kappa$ and $\omega$ are four-vectors orthogonal to fluid flow four-
vector $U$. For system with boost-invariant flow and transversely homogeneous
set-up the following basis are useful for further calculations

$\displaystyle U^{\alpha}=\left(\cosh(\eta),0,0,\sinh(\eta)\right),\quad
X^{\alpha}=\left(0,1,0,0\right),$ $\displaystyle
Y^{\alpha}=\left(0,0,1,0\right),\quad
Z^{\alpha}=\left(\sinh(\eta),0,0,\cosh(\eta)\right)$ (4)

with $\eta=\frac{1}{2}\ln((t+z)/(t-z))$ being the space-time rapidity. Using
these four-vector basis, $\kappa^{\alpha}$ and $\omega^{\alpha}$ can be
decomposed in the following way where all the scalar coefficients are function
of proper time $(\tau)$ only.

$\displaystyle\kappa^{\alpha}={C}_{\kappa X}(\tau)X^{\alpha}+{C}_{\kappa
Y}(\tau)Y^{\alpha}+{C}_{\kappa
Z}(\tau)Z^{\alpha},\quad\omega^{\alpha}={C}_{\omega
X}(\tau)X^{\alpha}+{C}_{\omega Y}(\tau)Y^{\alpha}+{C}_{\omega
Z}(\tau)Z^{\alpha}$ (5)

Using Eq. (5) in spin conservation law and projecting the resulting tensorial
equation on $U_{\beta}X_{\gamma}$, $U_{\beta}Y_{\gamma}$,
$U_{\beta}Z_{\gamma}$, $Y_{\beta}Z_{\gamma}$, $X_{\beta}Z_{\gamma}$, and,
$X_{\beta}Y_{\gamma}$ we get equations of motion for scalar coefficients which
turn out to evolve independently. Due to the rotational symmetry in the
transverse plane, coefficients ${C}_{\kappa X}$, ${C}_{\kappa Y}$ and
${C}_{\omega X}$, ${C}_{\omega Y}$ follow the same structure of differential
equations, respectively.

##  3 Particle spin polarization at freeze-out

Spin polarization tensor evolution permit us to calculate the mean spin
polarization per particle defined as
$\langle\pi_{\mu}\rangle=E_{p}\frac{d\Pi_{\mu}(p)}{d^{3}p}/E_{p}\frac{d{\cal{N}}(p)}{d^{3}p}$
[46] with

$\displaystyle E_{p}\frac{d\Pi_{\mu}(p)}{d^{3}p}$ $\displaystyle=$
$\displaystyle-\frac{\cosh(\xi)}{(2\pi)^{3}m}\int\Delta\Sigma_{\lambda}p^{\lambda}\,e^{-\beta\cdot
p}\,\tilde{\omega}_{\mu\beta}p^{\beta},$ $\displaystyle
E_{p}\frac{d{\cal{N}}(p)}{d^{3}p}$ $\displaystyle=$
$\displaystyle\frac{4\cosh(\xi)}{(2\pi)^{3}}\int\Delta\Sigma_{\lambda}p^{\lambda}\,e^{-\beta\cdot
p},$ (6)

where $E_{p}\frac{d\Pi_{\mu}(p)}{d^{3}p}$ being the total value of the Pauli-
Lubański vector for particles and $E_{p}\frac{d{\cal{N}}(p)}{d^{3}p}$ being
the momentum density of all particles having momentum $p$. Using canonical
boost, the polarization vector $\langle\pi^{\star}_{\mu}\rangle$ in the
particle rest frame can be obtained, where its longitudinal component is [47]

$\displaystyle\langle\pi^{\star}_{z}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{8m}\Big{[}\left(\frac{m\cosh(y_{p})+m_{T}}{m_{T}\cosh(y_{p})+m}\right)\left[\chi\left({C}_{\kappa
X}p_{y}-{C}_{\kappa Y}p_{x}\right)+2{C}_{\omega Z}m_{T}\right]$ (7)
$\displaystyle+\frac{\chi\,m\,\sinh(y_{p})\left({C}_{\omega
X}p_{x}+{C}_{\omega Y}p_{y}\right)}{m_{T}\cosh(y_{p})+m}\Big{]}$

where
$\chi\left(\hat{m}_{T}\right)=\left(K_{0}\left(\hat{m}_{T}\right)+K_{2}\left(\hat{m}_{T}\right)\right)/K_{1}\left(\hat{m}_{T}\right)$,
$\hat{m}_{T}=m_{T}/T$ and $y_{p}$ is the rapidity.

##  4 Results and Conclusion

For the Bjorken flow, charge current conservation and energy and linear
momentum conservation is written as $\frac{dn}{d\tau}+\frac{n}{\tau}=0$ and
$\frac{d\varepsilon}{d\tau}+\frac{(\varepsilon+P)}{\tau}=0$, respectively. In
Fig.1 (above panel left), proper-time dependence of temperature and baryon
chemical potential is shown, where temperature decreases with proper-time,
whereas ratio of baryon chemical potential over temperature increases with
proper time. In Fig.1 (above panel right), dependence of the $C$ coefficients
on proper time is shown. Using thermodynamical quantities and $C$ coefficients
at freeze-out, different components of the PRF average polarization vector
$\langle\pi^{\star}_{\mu}\rangle$ are obtained, see Fig.1 (below panel). We
notice that $\langle\pi^{\star}_{y}\rangle$ is negative which reflects the
initial system’s spin angular momentum. At midrapidity, the longitudinal
component of polarization vector ($\langle\pi^{\star}_{z}\rangle$) is zero and
$\langle\pi^{\star}_{x}\rangle$ shows quadrupole structure. Our results
presented here does not show quadrupole structure of the longitudinal
polarization observed in experiments because of the symmetries we assumed. We
presented results using the perfect-fluid hydrodynamics framework with spin
for the Bjorken scenario [50, 51]. Our approach show that scalar coefficients
$C$ which describes spin polarization evolve independently with weak proper-
time dependence. Our results may be used to determine the particle spin
polarization at freeze-out and we have shown that particle spin polarization
formed at freeze-out hypersurface depicts the initial polarization direction.
I thank Wojciech Florkowski and Radoslaw Ryblewski for clarifying and
inspiring discussions. This research was supported in part by the Polish
National Science Center Grants No. 2016/23/B/ST2/00717 and No.
2018/30/E/ST2/00432.

Figure 1: Above panel left: Proper-time dependence of temperature $T$ rescaled
by its initial value $T_{0}$ (solid line) and the ratio of $\mu$ (baryon
chemical potential) and $T$ (temperature) rescaled by the initial ratio
$\mu_{0}/T_{0}$ (dotted). Above panel right: Proper-time dependence of the
coefficients ${C}_{\kappa X}$, ${C}_{\kappa Z}$, ${C}_{\omega X}$ and
${C}_{\omega Z}$. Below panel: Components of the PRF average polarization
three-vector of $\Lambda$’s with the initial conditions $\mu_{0}=800$ MeV,
$T_{0}=155$ MeV, $C_{\kappa,0}=(0,0,0)$ and $C_{\omega,0}=(0,0.1,0)$ for mid-
rapidity.

## References

  * [1] L. Adamczyk et al. [STAR Collaboration], Nature 548, 62 (2017).
  * [2] J. Adam et al. [STAR Collaboration], Phys. Rev. C 98, 014910 (2018).
  * [3] T. Niida [STAR Collaboration], Nucl. Phys. A 982, 511 (2019).
  * [4] J. Adam et al. [STAR Collaboration], Phys. Rev. Lett.  123, no. 13, 132301 (2019).
  * [5] F. Becattini and L. Tinti, Annals Phys.  325, 1566 (2010).
  * [6] F. Becattini et al., Annals Phys.  338, 32 (2013).
  * [7] D. Montenegro et al., Phys. Rev. D 96, no. 5, 056012 (2017) Addendum: [Phys. Rev. D 96, no. 7, 079901 (2017)].
  * [8] D. Montenegro, L. Tinti and G. Torrieri, Phys. Rev. D 96, no. 7, 076016 (2017)
  * [9] F. Becattini, W. Florkowski and E. Speranza, Phys. Lett. B 789, 419 (2019).
  * [10] B. Boldizsár, M. I. Nagy and M. Csanád, Universe 5, no. 5, 101 (2019).
  * [11] S. Y. F. Liu, Y. Sun and C. M. Ko, Phys. Rev. Lett. 125, no.6, 062301 (2020)
  * [12] W. Florkowski et al., Phys. Rev. C 100, no. 5, 054907 (2019).
  * [13] H. Z. Wu et al., Phys. Rev. Research.  1, 033058 (2019).
  * [14] F. Becattini, G. Cao and E. Speranza, Eur. Phys. J. C 79, no. 9, 741 (2019).
  * [15] J. j. Zhang et al., Phys. Rev. C 100, no. 6, 064904 (2019).
  * [16] K. Fukushima and S. Pu, arXiv:2001.00359 [hep-ph].
  * [17] W. Florkowski, A. Kumar and R. Ryblewski, Acta Phys. Polon. B 51, 945-959 (2020).
  * [18] S. Li and H. U. Yee, Phys. Rev. D 100, no. 5, 056022 (2019).
  * [19] K. Hattori, Y. Hidaka and D. L. Yang, Phys. Rev. D 100, no. 9, 096011 (2019).
  * [20] N. Weickgenannt et al., Phys. Rev. D 100, no. 5, 056018 (2019).
  * [21] K. Hattori et al., Phys. Lett. B 795, 100 (2019).
  * [22] V. E. Ambrus, JHEP 08, 016 (2020).
  * [23] X. L. Sheng, L. Oliva and Q. Wang, Phys. Rev. D 101, no.9, 096005 (2020).
  * [24] Y. B. Ivanov, V. D. Toneev and A. A. Soldatov, Phys. Atom. Nucl. 83, no.2, 179-187 (2020).
  * [25] Y. Xie, D. Wang and L. P. Csernai, Eur. Phys. J. C 80, no.1, 39 (2020).
  * [26] G. Y. Prokhorov, O. V. Teryaev and V. I. Zakharov, Phys. Rev. D 99, no. 7, 071901 (2019).
  * [27] G. Y. Prokhorov, O. V. Teryaev and V. I. Zakharov, JHEP 1902, 146 (2019).
  * [28] G. Y. Prokhorov, O. V. Teryaev and V. I. Zakharov, JHEP 03, 137 (2020).
  * [29] D. L. Yang, Phys. Rev. D 98, no. 7, 076019 (2018).
  * [30] Y. Liu and X. Huang, Nucl. Sci. Tech. 31, no.6, 56 (2020).
  * [31] S. Tabatabaee and N. Sadooghi, Phys. Rev. D 101, no.7, 076022 (2020).
  * [32] S. Bhadury et al., [arXiv:2002.03937 [hep-ph]].
  * [33] Y. Liu, K. Mameda and X. Huang, Chin. Phys. C 44, 094101 (2020).
  * [34] D. Yang, K. Hattori and Y. Hidaka, JHEP 20, 070 (2020).
  * [35] X. Deng et al., Phys. Rev. C 101, no.6, 064908 (2020).
  * [36] H. Taya et al. [ExHIC-P], Phys. Rev. C 102, no.2, 021901 (2020).
  * [37] F. Becattini and M. A. Lisa, [arXiv:2003.03640 [nucl-ex]].
  * [38] F. Becattini et al., Phys. Rev. C 95, no. 5, 054902 (2017).
  * [39] I. Karpenko and F. Becattini, Eur. Phys. J. C 77, no. 4, 213 (2017).
  * [40] H. Li et al., Phys. Rev. C 96, no. 5, 054908 (2017).
  * [41] Y. Xie, D. Wang and L. P. Csernai, Phys. Rev. C 95, no. 3, 031901 (2017).
  * [42] F. Becattini, F. Piccinini and J. Rizzo, Phys. Rev. C 77, 024906 (2008).
  * [43] W. Florkowski et al., Phys. Rev. C 97, no. 4, 041901 (2018).
  * [44] W. Florkowski et al., Phys. Rev. D 97, no. 11, 116017 (2018).
  * [45] W. Florkowski, R. Ryblewski and A. Kumar, Prog. Part. Nucl. Phys.  108, 103709 (2019).
  * [46] W. Florkowski, A. Kumar and R. Ryblewski, Phys. Rev. C 98, no. 4, 044906 (2018).
  * [47] W. Florkowski et al., Phys. Rev. C 99, no. 4, 044910 (2019).
  * [48] S. R. De Groot, W. A. Van Leeuwen, C. G. Van Weert, Relativistic Kinetic Theory, Principles and Applications, Amsterdam, North-Holland, 1980.
  * [49] W. Florkowski, Phenomenology of Ultra-Relativistic Heavy-Ion Collisions, Singapore: World Scientific, 2010.
  * [50] J. D. Bjorken, Phys. Rev. D 27, 140 (1983).
  * [51] L. Rezzolla and O. Zanotti, Relativistic Hydrodynamics, Oxford University Press, 2013.