Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid,
E-28049 Madrid, Spain

# Shift current with Gaussian basis sets & general prescription for maximally-
symmetric summations in the irreducible Brillouin zone

M. A. García-Blázquez<EMAIL_ADDRESS>Departamento de
Física de la Materia Condensada, Universidad Autónoma de Madrid, E-28049
Madrid, Spain  J. J. Esteve-Paredes  Departamento de Física de la Materia
Condensada, Universidad Autónoma de Madrid, E-28049 Madrid, Spain  A. J. Uría
Departamento de Física de la Materia Condensada, Universidad Autónoma de
Madrid, E-28049 Madrid, Spain  J. J. Palacios  Departamento de Física de la
Materia Condensada, Universidad Autónoma de Madrid, E-28049 Madrid, Spain

###### Abstract

The bulk photovoltaic effect is an experimentally verified phenomenon by which
a direct charge current is induced within a non-centrosymmetric material by
light illumination. Calculations of its intrinsic contribution, the shift
current, are nowadays amenable from first-principles employing plane-waves
bases. In this work we present a general method for evaluating the shift
conductivity in the framework of localized Gaussian basis sets that can be
employed in both the length and velocity gauges, carrying the idiosyncrasies
of the quantum-chemistry approach. The (possibly magnetic) symmetry of the
system is exploited in order to fold the reciprocal space summations to the
representation domain, allowing to reduce computation time and unveiling the
complete symmetry properties of the conductivity tensor under general light
polarization.

##  1 Introduction

The generation of a non-oscillating response in a material medium under an
incident electric oscillating field is a general feature that occurs at all
even orders in the perturbative expansion, that is, it is a non-linear optical
phenomenon. For responses that transform as vectors, such as an electric
current, an elemental symmetry analysis shows that all the even order response
tensors must vanish in the presence of inversion symmetry, hence such
frequency-independent quantity can only arise in non-centrosymmetric
materials. In this regard, the emergence of a direct charge current in an
homogeneous material induced by light is known as the bulk photovoltaic effect
(BPVE) 1. It was established by the mid seventies with earlier experimental
reports in ferroelectric materials 2, 3, 4, 5, and continued gathering
attention during the next decades 6, 7, 8, 9, 10, 11, 12, 13, 14. However, the
potential applications in solar cells 15, 16, and the advances in both
experimental facilities and first-principles capabilities have driven a
considerable surge of studies in recent years 17, 18, 19, 20, 21, 22, 23, 24,
25, 26.

The BPVE is part of the total second-order optical response, which in addition
includes second-harmonics and, for polychromatic electric fields,
contributions of mixed frequencies. In turn, the BPVE can be separated into 3
essentially different contributions 17: the shift current, a static and
coherent (stemming from the off-diagonal part of the density matrix) response
that under time-reversal ($\hat{\pazocal{T}}$) symmetry appears only with
linearly polarized light; and two transient contributions that eventually
reach a steady state, namely the injection current, which under
$\hat{\pazocal{T}}$ symmetry appears only with circularly polarized light, and
the ballistic current, which under $\hat{\pazocal{T}}$ symmetry emerges purely
from coherent scattering processes such as electron-phonon or electron-hole
interactions that introduce an imbalance between the carrier generation rates
across the Brillouin zone (BZ). Disregarding excitonic effects, the BPVE in
non-metallic systems occurs at frequencies above the band gap. These
quantities, or equivalently the corresponding third-rank tensors
$\sigma^{a;bc}(\omega)$ as a function of a single frequency, admit expressions
in terms of the quasi-particle properties that are are amenable to numerical
evaluation via quantum mechanical methods. Specifically, these microscopic
expressions can be obtained by diagrammatic approaches for the ballistic
current 27, 28, 29, and solving the density matrix perturbatively 1, employing
Wilson loops 30 or again by diagrammatic techniques 31 for the injection and
shift currents. However, only the latter one is truly intrinsic to the single-
particle system, in the sense that it can be computed exclusively from the
band structure and electronic eigenfunctions without further modelling.

The calculation of the shift current presents some difficulties or subtleties
starting from the choice of gauge for the interaction of electrons with the
field 32, 33, 34, 31. The most generally applicable method, the length gauge,
requires evaluating numerical derivatives with respect to the crystalline
momentum $\bm{k}$ of quantities that are not gauge invariant. In contrast, the
velocity gauge constitutes a more straightforward alternative, although it
carries an (a priori) infinite sum over the electronic states external to the
direct optical transition. Both gauges require the evaluation of the matrix
representation of the velocity operator in the set of crystalline
eigenfunctions, and both are expected to yield equal results in the limit of a
complete basis for describing the latter. There currently exist methods for
evaluating the shift current in the single-particle approximation from
density-functional theory (DFT) 35, tight-binding (including Wannierizations
36, 37) and $\bm{k}\cdot\bm{p}$ band structures 18. Yet, as it is often the
case in physics-leaning studies, the DFT calculations are almost invariably
assumed to employ a plane-wave basis, at least when the velocity operator is
not approximated by the momentum.

In this work, we present a formalism for computing the shift conductivity
tensor in non-metallic crystals, in both length and velocity gauges, from
first-principles employing Gaussian basis sets. It is based on an exact
calculation of the velocity and Berry connection matrix elements through the
analytical evaluation of the real-space integrals involved. The use of a
localized basis presents some advantages and disadvantages with respect to the
plane-waves alternative inherited from the DFT methods: the whole chain of
calculations should generally be faster, the evaluation of position (and by
extension, velocity) matrix elements is straightforward, hybrid functionals
can be used at little cost (which may allow to obtain an accurate gap avoiding
scissor corrections or GW calculations), all-electron calculations can be
performed, and no artificial replication of layers is required in 2D
materials. On the other hand, a customized basis optimization has to be
performed for each system while limited by the superposition error and
diffusive exponents, and errors from the lack of completeness of the basis are
more likely (making the more delocalized unoccupied states particularly hard
to reproduce). It is expected that the reliability of this method is highly
correlated with the ability to properly reproduce the relevant occupied and
unoccopied states (dictated by the frequency range) with a Gaussian basis.

A further benefit of the use of localized bases lies in the guarantee that the
complete symmetry of the system is preserved, in slight contrast with the
maximally-localized Wannier representation. The properties of the
crystallographic point group can then be exploited to reduce the summations
over the BZ that are required for the shift conductivity to properly weighted
sums, which encode the whole (possibly magnetic) symmetry of the system, only
over its irreducible part or representation domain. We present a complete list
of the explicit formulae for each space group including the magnetic
configurations, where the structure type is used to parametrize the
irreducible domain and the (magnetic) point group determines the precise
folding of the $\bm{k}-$resolved conductivity.

##  2 Shift current: definition, considerations & numerical evaluation

A general expression for the total second-order optical response under
homogeneous illumination can be obtained by solving the density matrix in
perturbation theory for the field. In particular, for a uniform polychromatic
electric field
$\tilde{\bm{E}}(t)=\sum_{j}\bm{E}(\omega_{j})e^{-i\omega_{j}t}+\text{c.c.}$,
the shift current is defined as the intrinsic second-order DC component

$J^{a}_{\text{shift}}=2\sum_{j\>|\>\omega_{j}>0}\sum_{b,c}\real\left[\sigma^{a;bc}_{\text{shift}}(\omega_{j})E^{b}(\omega_{j})E^{c}(-\omega_{j})\right]$
(1)

where $a,b,c$ label the spatial components in the chosen coordinate system.

###  2.1 Length gauge

In the length gauge, the electric potential is chosen as
$\hat{V}=\absolutevalue{e}\hat{\bm{r}}\cdot\bm{E}$ and the perturbative
expression for the shift conductivity third-rank tensor in a non-metallic
material ultimately reads 38, 14, 21

$\sigma^{a;bc}_{\text{shift}}(\omega)=-\frac{ig_{s}\pi\absolutevalue{e}^{3}}{2\hbar^{2}V}\sum_{\bm{k}\in\text{BZ}}\sum_{m,n}f_{m,n}(\bm{k})\left[A^{b}_{m,n}(\bm{k})A^{c}_{n,m;a}(\bm{k})-A^{c}_{n,m}(\bm{k})A^{b}_{m,n;a}(\bm{k})\right]\delta(\omega_{m,n}(\bm{k})-\omega)$
(2)

which is valid irrespective of whether time-reversal $\hat{\pazocal{T}}$ is a
symmetry. In this formula

  * •

$m,n$ label the eigenstates of the periodic single-particle Hamiltonian
$\hat{H}\ket{\psi_{n,\bm{k}}}=E_{n,\bm{k}}\ket{\psi_{n,\bm{k}}}$, which
satisfy Bloch’s theorem:
$\bra{\bm{r}}\ket{\psi_{n,\bm{k}}}=\psi_{n,\bm{k}}(\bm{r})=e^{i\bm{k}\bm{r}}u_{n,\bm{k}}(\bm{r})$
with $u_{n,\bm{k}}$ having the periodicity of the direct lattice. $\bm{k}$ is
the crystalline momentum or label of the irreducible representations of the
translation group. $\hbar\omega_{m,n}(\bm{k})\equiv E_{m,\bm{k}}-E_{n,\bm{k}}$
and $f_{m,n}(\bm{k})\equiv f_{m}(\bm{k})-f_{n}(\bm{k})$ is the difference of
Fermi distributions, hereafter taken at zero temperature.

  * •

$A^{b}_{m,n}(\bm{k})=\left.\bra{u_{m,\bm{k}^{\prime}}}i\hat{\partial}_{k^{\prime
b}}\ket{u_{n,\bm{k}^{\prime}}}\right|_{\bm{k}}$ is the $b-$th spatial
component of the Berry connection matrix elements, which satisfy
$A^{b}_{m,n}(\bm{k})=A^{b}_{n,m}(\bm{k})^{*}$.

  * •

$A^{b}_{n,m;a}(\bm{k})=\left.\partial_{k^{\prime
a}}A^{b}_{n,m}(\bm{k}^{\prime})\right|_{\bm{k}}-i[A^{a}_{n,n}(\bm{k})-A^{a}_{m,m}(\bm{k})]A^{b}_{n,m}(\bm{k})$
is the generalized derivative (GD) of the Berry connection.

  * •

$V=N_{\bm{k}}V_{\text{PUC}}$ is the volume (area in 2D, or longitude in 1D) of
the crystal, with $N_{\bm{k}}\to\infty$ the number of terms in the Brillouin
zone (BZ) summation (or discretized integration) and $V_{\text{PUC}}$ the
volume of the primitive unit cell. $\delta(\omega)$ is a nascent (or
broadened) Dirac delta function. $g_{s}=1$ (2) in the presence (absence,
respectively) of spin-dependent terms in the Hamiltonian (excluding the
doubled states in the $m,n$ summations).

  * •

We have introduced a global minus sign in agreement with Ref. 21 and equation
58 of Ref. 14. We note that some authors do not include the $1/2$ factor in
the conductivity, instead cancelling it with the prefactor in (1).

It follows that in general
$\sigma^{a;bc}_{\text{shift}}(\omega)=\sigma^{a;cb}_{\text{shift}}(\omega)^{*}=\sigma^{a;bc}_{\text{shift}}(-\omega)^{*}$
and $\sigma^{a;bb}_{\text{shift}}$ is real. For linearly polarized light,
$E^{b}(\omega_{j})E^{c}(-\omega_{j})$ is real for all components and only
$\real\sigma^{a;bc}_{\text{shift}}$ contributes to $J^{a}_{\text{shift}}$. In
contrast, for circular polarization $E^{b}(\omega_{j})E^{c}(-\omega_{j})$ is
complex for some $b,c$, hence both the real and imaginary parts of the
conductivity may contribute to the current in general. Under
$\hat{\pazocal{T}}$ symmetry, in particular excluding any permanent magnetic
alignment, $\hat{\pazocal{T}}u_{n,\bm{k}}=e^{i\theta_{n,\bm{k}}}u_{n,-\bm{k}}$
for an arbitrary phase $\theta_{n,\bm{k}}$; hence applying the anti-unitary
transformation in the inner products, and noting that
$\hat{\pazocal{T}}i\hat{\bm{\nabla}}_{\bm{k}}\hat{\pazocal{T}}^{-1}=-i\hat{\bm{\nabla}}_{\bm{k}}$

$\displaystyle
A^{b}_{m,n}(\bm{k})=\left.\bra{u_{m,-\bm{k}^{\prime}}}e^{-i\theta_{m,\bm{k}}}-i\hat{\partial}_{k^{\prime
b}}e^{i\theta_{n,\bm{k}}}\ket{u_{n,-\bm{k}^{\prime}}}^{*}\right|_{\bm{k}}=e^{i\theta_{n,m}^{\bm{k}}}\left[A^{b}_{m,n}(-\bm{k})^{*}+\frac{\partial\theta_{n,\bm{k}}}{\partial
k^{b}}\delta_{m,n}\right]\>\Rightarrow$ $\displaystyle
A^{c}_{n,m;a}(\bm{k})=-e^{i\theta_{m,n}^{\bm{k}}}\left[\left.\frac{\partial
A^{c}_{n,m}(\bm{k}^{\prime})^{*}}{\partial k^{\prime
a}}\right|_{-\bm{k}}+i\left(A^{a}_{n,n}(-\bm{k})^{*}-A^{a}_{m,m}(-\bm{k})^{*}\right)A^{c}_{n,m}(-\bm{k})^{*}\right]=-e^{i\theta_{m,n}^{\bm{k}}}A^{c}_{n,m;a}(-\bm{k})^{*}$

where $\theta_{n,m}^{\bm{k}}\equiv\theta_{n,\bm{k}}-\theta_{m,\bm{k}}$. Thus
$f_{m,n}(\bm{k})A^{b}_{m,n}(\bm{k})A^{c}_{n,m;a}(\bm{k})=-\left[f_{m,n}(-\bm{k})A^{b}_{m,n}(-\bm{k})A^{c}_{n,m;a}(-\bm{k})\right]^{*}$
and $\sigma^{a;bc}_{\text{shift}}$ is real with $\hat{\pazocal{T}}$ symmetry.
In this case also
$-A^{c}_{n,m}(\bm{k})A^{b}_{m,n;a}(\bm{k})=A^{c}_{m,n}(\bm{k})A^{b}_{n,m;a}(\bm{k})$
under the BZ summation. The difference between the currents for right and left
circular polarization, which is proportional to
$\imaginary\sigma^{a;bc}_{\text{shift}}$ only, is therefore vanishing; and the
shift current is associated with linearly polarized light. In the event of
$\hat{\pazocal{T}}$ breaking, a circular shift current generally emerges as
well as a non-stationary injection current 21, 23. We remark that the
injection current can be equally computed with the method described in this
work, but we do not show explicit results since the calculation is
straightforward compared to that for the shift current 38, 14, up to an
extrinsic scattering rate that is often set phenomenologically.

###  2.2 Berry connection and velocity in a local basis

The single-particle crystalline eigenstates are generally expanded in a set of
states 111Which in practice is not an actual (complete) basis in this function
space due to its finiteness. Nevertheless, we use this term referring to non-
complete sets according to standard convention. satisfying Bloch’s theorem as

$\ket{\psi_{n,\bm{k}}}=\sum_{\mu}c_{\mu,n}(\bm{k})\ket{\phi_{\mu,\bm{k}}}$

where $\mu$ is in principle a generic label and the coefficients are obtained
from the generalized eigenvalue problem

$\sum_{\mu^{\prime}}\left[H_{\mu,\mu^{\prime}}(\bm{k})-E_{n}(\bm{k})S_{\mu,\mu^{\prime}}(\bm{k})\right]c_{\mu^{\prime},n}(\bm{k})=0$

where the Hamiltonian $H_{\mu,\mu^{\prime}}(\bm{k})$ and overlap
$S_{\mu,\mu^{\prime}}(\bm{k})$ matrix elements are the representations of the
Hamiltonian $\hat{H}$ and identity $\hat{I}$ operator, respectively, in the
set $\set{\ket{\phi_{\mu,\bm{k}}}}_{\mu}$ for each $\bm{k}$. In a local basis,
which is repeated in each unit cell and labelled by the lattice vectors
$\bm{R}$, the Bloch states can in turn be expanded in agreement with Bloch’s
theorem

$\ket{\phi_{\mu,\bm{k}}}=\frac{1}{\sqrt{N}}\sum_{\bm{R}}e^{i\bm{k}\bm{R}}\ket{\varphi_{\mu,\bm{R}}}$

where $N=V/V_{\text{PUC}}\to\infty$ is the number of unit cells in the
crystal. Therefore

$\displaystyle
H_{\mu,\mu^{\prime}}(\bm{k})=\sum_{\bm{R}}e^{i\bm{k}\bm{R}}H_{\mu,\mu^{\prime}}(\bm{R}),\text{
with
}H_{\mu,\mu^{\prime}}(\bm{R})=\bra{\varphi_{\mu,\bm{0}}}\hat{H}\ket{\varphi_{\mu^{\prime},\bm{R}}},$
(3) $\displaystyle
S_{\mu,\mu^{\prime}}(\bm{k})=\sum_{\bm{R}}e^{i\bm{k}\bm{R}}S_{\mu,\mu^{\prime}}(\bm{R}),\text{
with
}S_{\mu,\mu^{\prime}}(\bm{R})=\bra{\varphi_{\mu,\bm{0}}}\ket{\varphi_{\mu^{\prime},\bm{R}}}$

due to the periodicity of $\hat{H}(\bm{r})$ and $\hat{I}$.

The Berry connection matrix elements can then be expressed in the local basis
by inserting the transformation
$u_{n,\bm{k}}(\bm{r})=e^{-i\bm{k}\bm{r}}\sum_{\mu}c_{\mu,n}(\bm{k})\frac{1}{\sqrt{N}}\sum_{\bm{R}}e^{i\bm{k}\bm{R}}\varphi_{\mu,\bm{R}}(\bm{r})$,
yielding after some algebra

$A^{b}_{m,n}(\bm{k})=-i\left.\frac{\partial
c^{*}_{\mu,m}(\bm{k}^{\prime})}{\partial k^{\prime
b}}\right|_{\bm{k}}S_{\mu,\mu^{\prime}}(\bm{k})c_{\mu^{\prime},n}(\bm{k})+c^{*}_{\mu,m}(\bm{k})r^{b}_{\mu,\mu^{\prime}}(\bm{k})c_{\mu^{\prime},n}(\bm{k})$

which can be compactly expressed in matrix form as

$A^{b}(\bm{k})=-i\left.\frac{\partial c^{\dagger}(\bm{k}^{\prime})}{\partial
k^{\prime
b}}\right|_{\bm{k}}S(\bm{k})c(\bm{k})+c^{\dagger}(\bm{k})r^{b}(\bm{k})c(\bm{k})$
(4)

We have introduced the position matrix elements in the Bloch basis

$\bm{r}_{\mu,\mu^{\prime}}(\bm{k})=\sum_{\bm{R}}e^{i\bm{k}\bm{R}}\bm{r}_{\mu,\mu^{\prime}}(\bm{R}),\text{
with
}\bm{r}_{\mu,\mu^{\prime}}(\bm{R})=\bra{\varphi_{\mu,\bm{0}}}\hat{\bm{r}}\ket{\varphi_{\mu^{\prime},\bm{R}}}=\int_{\mathbb{R}^{3}}\varphi_{\mu,\bm{0}}(\bm{r})^{*}\bm{r}\varphi_{\mu^{\prime},\bm{R}}(\bm{r})d^{3}\bm{r}$
(5)

which is well defined in a localized basis, albeit the diagonal components
depend on the origin choice. Indeed, a rigid shift of the form
$\bm{r}\rightarrow\bm{r}+\bm{r}_{0}$ (which does not alter the position
operator $\hat{\bm{r}}$ itself) with $\bm{r}_{0}$ restricted to the unit cell,
results in $A^{b}_{m,n}(\bm{k})\rightarrow
A^{b}_{m,n}(\bm{k})-\bm{r}_{0}\delta_{m,n}$. Nevertheless, it is easy to see
that this arbitrary factor is cancelled in the GD $A^{b}_{m,n;a}$ rendering
the shift conductivity invariant under this choice.

As it can be observed in (4), the calculation of both diagonal and off-
diagonal Berry connections requires evaluating numerical derivatives with
respect to $\bm{k}$. In order to avoid further complications in the term
$\partial_{k^{a}}A^{b}_{n,m}$ of the GD $A^{b}_{n,m;a}$, such as the
introduction of a second grid for derivatives or the increase in the required
significant digits, we employ the identity

$\bm{A}_{n,m}(\bm{k})=\frac{\bm{v}_{n,m}(\bm{k})}{i\omega_{n,m}(\bm{k})}\>,\;\;\;m\neq
n$ (6)

which can be readily obtained by expanding
$\bm{\nabla}_{\bm{k}}\left[\bra{u_{m,\bm{k}}}e^{-i\bm{k}\bm{r}}\hat{H}e^{i\bm{k}\bm{r}}\ket{u_{n,\bm{k}}}\right]=0$
for $m\neq n$. In (6),
$\bm{v}_{n,m}(\bm{k})=\bra{\psi_{n,\bm{k}}}\bm{\hat{v}}\ket{\psi_{m,\bm{k}}}=\bm{v}_{m,n}(\bm{k})^{*}$
are the matrix elements (here in the set of eigenstates) of the velocity
operator $\hat{\bm{v}}\equiv(i/\hbar)[\hat{H},\hat{\bm{r}}]$. It can be shown
that the velocity matrix elements have the following form in the local basis
39

$\bm{v}_{n,m}(\bm{k})=\frac{1}{\hbar}\sum_{\mu,\mu^{\prime}}c^{*}_{\mu,n}(\bm{k})\left[\bm{\nabla}_{\bm{k}}H_{\mu,\mu^{\prime}}(\bm{k})-E_{n}(\bm{k})\bm{\nabla}_{\bm{k}}S_{\mu,\mu^{\prime}}(\bm{k})+i\hbar\omega_{n,m}(\bm{k})\bm{r}_{\mu,\mu^{\prime}}(\bm{k})\right]c_{\mu^{\prime},m}(\bm{k})$
(7)

which is independent on the origin choice. The $\bm{k}-$derivatives in this
expression are all analytical, in particular
$\bm{\nabla}_{\bm{k}}H_{\mu,\mu^{\prime}}(\bm{k})=i\sum_{\bm{R}}e^{i\bm{k}\bm{R}}\bm{R}H_{\mu,\mu^{\prime}}(\bm{R})$
and likewise for $S_{\mu,\mu^{\prime}}(\bm{k})$. Therefore,
$\partial_{k^{a}}A^{b}_{n,m}$ and $A^{b}_{n,m}$ can be computed employing (6)
with numerical derivatives only of the first and zeroth order, respectively.
Inserting these into (2), noting that the $f_{m,n}$ factor forces $m\neq n$
and that the derivatives of $\omega_{n,m}$ cancel out, one obtains

$\displaystyle\sigma^{a;bc}_{\text{shift}}(\omega)=-\frac{ig_{s}\pi\absolutevalue{e}^{3}}{2\hbar^{2}V}\sum_{\bm{k}\in\text{BZ}}\sum_{m,n}\frac{f_{m,n}}{\omega^{2}_{m,n}}\left[v^{b}_{m,n}v^{c}_{n,m;a}-v^{c}_{n,m}v^{b}_{m,n;a}\right]\delta(\omega_{m,n}-\omega)=$
(8)
$\displaystyle-\frac{ig_{s}\pi\absolutevalue{e}^{3}}{2\hbar^{2}V}\sum_{\bm{k}\in\text{BZ}}\sum_{m,n}\frac{f_{m,n}}{\omega^{2}_{m,n}}\left\\{v^{b}_{m,n}\left[\frac{\partial
v^{c}_{n,m}}{\partial
k^{a}}-iv^{c}_{n,m}\left(A^{a}_{n,n}-A^{a}_{m,m}\right)\right]-(b\leftrightarrow
c)^{*}\right\\}\delta(\omega_{m,n}-\omega)$

where $(b\leftrightarrow c)^{*}$ represents the conjugated of the previous
term inside the curly brackets with the components $b$ and $c$ permuted (even
if $b=c$), and the $\bm{k}$ dependence has been omitted for brevity. Note that
with $\hat{\pazocal{T}}$ symmetry the $-(b\leftrightarrow c)^{*}$ term is
equivalent to $+(b\leftrightarrow c)$ under the BZ summation.

A subtle issue in (8) and other equivalent length gauge formulae is the
evaluation of numerical derivatives of quantities that are not gauge-
invariant, in particular the coefficients $c_{\mu,m}$ in the diagonal Berry
connections, see (4), and the velocities $\bm{v}_{n,m}$ for $n\neq m$; which
are respectively determined up to arbitrary phase factors $\theta_{m,\bm{k}}$
and $\theta_{m,n}^{\bm{k}}=\theta_{m,\bm{k}}-\theta_{n,\bm{k}}$. The
continuity of (LABEL:HkSk) in $\bm{k}$ makes all $c(\bm{k})$ (and
$\bm{v}(\bm{k})$ by extension) also continuous, except in general for the
phases. In order to fix a continuous gauge, we impose that

$e^{i\theta_{m,\bm{k}}}\sum_{\mu}c_{\mu,m}(\bm{k})=\absolutevalue{\sum_{\mu}c_{\mu,m}(\bm{k})}\in\mathbb{R}\>,\;\;\forall
m,\bm{k}$

The necessary $\bm{k}-$ derivatives are well defined this way, and the
$\theta_{m,\bm{k}}$ factors are all cancelled in the gauge-invariant (8) by a
similar argument than in the time-reversal case above. The only remaining
caveat is to keep track of the correct band indexing when a degeneracy occurs
between the infinitesimally close $\bm{k}\pm\bm{h}$ points defining the
numerical derivatives. However, this issue may be neglected by mapping the
$\bm{k}-$summation to the interior of the irreducible Brillouin zone (IBZ),
where only accidental degeneracies may occur, see Section 3.

In 2D materials the out of plane tensor components, i.e., involving at least
one index along the non-periodic direction $\bm{z}$, can be computed on equal
footing (and likewise for 1D systems). Regarding the system as a periodic
stacking of layers, $H_{\mu,\mu^{\prime}}(\bm{R}_{z})$ 222In charge-neutral 2D
systems, the Coulomb potential may be replaced by the Parry potential instead
of the usual Ewald potential in 3D 40., $S_{\mu,\mu^{\prime}}(\bm{R}_{z})$,
$\bm{r}_{\mu,\mu^{\prime}}(\bm{R}_{z})$ are exponentially vanishing for inter-
layer vectors $\bm{R}_{z}$ in the limit of large layer separation, making all
$k^{z}-$derivatives null. In this case (4) and (7) reduce to
$A^{z}_{m,n}=-iv_{m,n}^{z}/\omega_{m,n}=\sum_{\mu,\mu^{\prime}}c^{\dagger}_{m,\mu}r^{z}_{\mu,\mu^{\prime}}c_{\mu^{\prime},n}$.
While the $a=z$ component in the shift conductivity tensor may not be of
interest, in some point groups the IBZ summation requires the calculation of
some of these components for $a=x,y$; in which case the numerical derivatives
in (8) are cancelled and the expression is simplified significantly.

An alternative treatment of (2) that is often found in the literature 1, 14,
35 consists on the introduction of the shift vector, which involves the term
$\bm{\partial_{k^{a}}}\Phi^{b}_{n,m}(\bm{k})$ where
$v^{b}_{n,m}=\absolutevalue{v^{b}_{n,m}}e^{-i\Phi^{b}_{n,m}}$. This is
obtained by noting that in the GD
$\partial_{k^{a}}A^{b}_{n,m}=A^{b}_{n,m}\partial_{k^{a}}\log A^{b}_{n,m}$ when
$A^{b}_{n,m}\neq 0$. The term inside the curly brackets in (8) is then
equivalent to

$-iv^{b}_{m,n}v^{c}_{n,m}\left[\frac{\partial(\Phi^{b}_{n,m}+\Phi^{c}_{n,m})}{\partial
k^{a}}+2(A^{a}_{n,n}-A^{a}_{m,m})+i\frac{\partial\log\absolutevalue{v^{c}_{n,m}/v^{b}_{n,m}}}{\partial
k^{a}}\right]$

For linearly polarized light one can always rotate the coordinate system,
initially based on the crystallographic structure, such that $\bm{E}(\omega)$
points along, say, the $b$ direction. Then only the $\sigma^{a;bb}$ component
contributes to the current along $a$ in (1), and the computation of the tensor
is slightly simplified; in particular the modulus term vanishes in the last
expression. This is, however, not advisable for practical calculations because
it requires evaluating the tensor for each field direction with a different
coordinate system, which may also hinder the obtention of the Hamiltonian
matrix elements from the electronic structure code. In this work (8) is
employed instead, since the computational cost is similar.

###  2.3 Evaluation in Gaussian basis sets

The evaluation of (8) from first-principles requires thus the knowledge of the
matrix elements of $\hat{H}$, $\hat{I}$ and $\hat{\bm{r}}$ in the local basis
for a sufficiently large number of lattice vectors The first one,
$H_{\mu,\mu^{\prime}}(\bm{R})$, must be evaluated self-consistently, typically
in a DFT or hybrid DFT-HF (Hartree-Fock) scheme; and is generally expected to
be provided by the corresponding electronic structure code for the chosen
functional. The latter two, $S_{\mu,\mu^{\prime}}(\bm{R})$ and
$\bm{r}_{\mu,\mu^{\prime}}(\bm{R})$ can be manually pre-computed from the
(possibly optimized) atomic structure. If the local functions are harmonic
Gaussian-type orbitals (GTOs), this can be done analytically. In that case
$\mu$ is a multi-index labelling the atoms $a$ (located at $\bm{d}_{a}$) in
the unit cell 333Ghost atoms would be treated analogously as long as they
preserve the space group symmetry., the pair of orbital quantum numbers $l,m$
($m=-l,\dots,l$), the shells $\lambda$ discerning the harmonics with identical
$l$ and, in the presence of spin dependent terms in $\hat{H}$ such as spin-
orbit coupling (SOC) or magnetic ordering, the $m_{s}=\pm 1/2$ spin quantum
number. The contracted real GTOs
$\varphi_{\mu,\bm{R}}:\mathbb{R}^{3}\to\mathbb{R}$ are then defined as 41, 42

$\varphi_{\mu,\bm{R}}(\bm{r})=N_{\lambda,l}\left[\sum_{j}c_{l,m,j}d_{\lambda,j}G(\alpha_{\lambda,j},\bm{r}-\bm{d}_{a}-\bm{R})\right]X_{l,m}(\bm{r}-\bm{d}_{a}-\bm{R})$

where $N_{\lambda,l}$ and $c_{l,m,j}$ are normalization coefficients (see
Appendix E of Ref. 41), $d_{\lambda,j}$ and $\alpha_{\lambda,j}$ are the
selected contraction coefficients and exponents, $G(\alpha,\bm{r})=e^{-\alpha
r^{2}}$ are the Gaussian-type radial functions and $X_{l,m}(\bm{r})$ are the
real solid harmonics. The latter are obtained from the (not normalized)
spherical harmonics $Y_{l,m}(\bm{r})$ as

$X_{l,m}(\bm{r})=\frac{r^{l}}{2}\cdot\left\\{\begin{aligned}
&Y_{l,\absolutevalue{m}}(\bm{r})+Y_{l,-\absolutevalue{m}}(\bm{r})\>,&\text{ if
}m\geq 0\\\
&-i\left[Y_{l,\absolutevalue{m}}(\bm{r})-Y_{l,-\absolutevalue{m}}(\bm{r})\right]\>,&\text{
if }m<0\end{aligned}\right.$

The central integrals that need be evaluated to obtain $S(\bm{k})$ and
$\bm{r}(\bm{k})$ are then

$\displaystyle
I_{3}(\bm{n},\bm{n}^{\prime},\bm{r}_{0},\bm{r}^{\prime}_{0},\alpha,\alpha^{\prime})=I(n_{x},n^{\prime}_{x},x_{0},x^{\prime}_{0},\alpha,\alpha^{\prime})\cdot
I(n_{y},n^{\prime}_{y},y_{0},y^{\prime}_{0},\alpha,\alpha^{\prime})\cdot
I(n_{z},n^{\prime}_{z},z_{0},z^{\prime}_{0},\alpha,\alpha^{\prime})\>,$ (9)
$\displaystyle
I(n_{x},n^{\prime}_{x},x_{0},x^{\prime}_{0},\alpha,\alpha^{\prime})\equiv\int_{\mathbb{R}}(x-x_{0})^{n_{x}}(x-x^{\prime}_{0})^{n^{\prime}_{x}}e^{-\alpha(x-x_{0})^{2}-\alpha^{\prime}(x-x^{\prime}_{0})^{2}}dx=I(n_{x},n^{\prime}_{x},x_{0}-x^{\prime}_{0},0,\alpha,\alpha^{\prime})$

The one-dimensional integrals appearing in $\bm{r}(\bm{k})$ are then computed
as

$\int_{\mathbb{R}}(x-x_{0})^{n}(x-x^{\prime}_{0})^{n^{\prime}}xe^{-\alpha(x-x_{0})^{2}-\alpha^{\prime}(x-x^{\prime}_{0})^{2}}dx=x^{\prime}_{0}I(n_{x},n^{\prime}_{x},x_{0}-x^{\prime}_{0},0,\alpha,\alpha^{\prime})+I(n_{x},n^{\prime}_{x}+1,x_{0}-x^{\prime}_{0},0,\alpha,\alpha^{\prime})$

There are several methods to tabulate the one-dimensional integrals in (9),
e.g., by recursion over $n$ and $n^{\prime}$. In this work we have instead
employed the following master expression, which can be deduced from Ref. 43

$\displaystyle I(n,n^{\prime},x_{0},0,\alpha,\alpha^{\prime})=$ $\displaystyle
e^{-\tilde{\alpha}\>\alpha^{\prime}x_{0}^{2}}\>x_{0}^{n+n^{\prime}}\sqrt{\frac{\pi\tilde{\alpha}}{\alpha}}\>n!\sum_{k=0}^{n}(-1)^{n+k}\frac{(n^{\prime}+k)!}{(n-k)!\>k!}\>\tilde{\alpha}^{n^{\prime}+k}\sum_{h=0}^{\lfloor\frac{n^{\prime}+k}{2}\rfloor}\frac{1}{(n^{\prime}+k-2h)!\>h!}\left(\frac{1}{4\tilde{\alpha}\>\alpha\>x_{0}^{2}}\right)^{h}\>,\;\forall
x_{0}\neq 0$

where $\tilde{\alpha}\equiv\alpha/(\alpha+\alpha^{\prime})$ and

$I(n,n^{\prime},0,0,\alpha,\alpha^{\prime})=\left\\{\begin{aligned}
&\sqrt{\frac{\pi}{2^{n+n^{\prime}}(\alpha+\alpha^{\prime})^{n+n^{\prime}+1}}}(n+n^{\prime}-1)!!\>,&\text{
for }n+n^{\prime}\text{ even}\\\ &0\>,&\text{ for }n+n^{\prime}\text{
odd}\end{aligned}\right.$

We note that only one of $S_{\mu,\mu^{\prime}}(\pm\bm{R})$ (and likewise for
$\bm{r}_{\mu,\mu^{\prime}}(\pm\bm{R})$) needs be computed for each
$\bm{R}\neq\bm{0}$ since

$S_{\mu,\mu^{\prime}}(-\bm{R})=S_{\mu^{\prime},\mu}(\bm{R})\>,\;\;\bm{r}_{\mu,\mu^{\prime}}(-\bm{R})=\bm{r}_{\mu^{\prime},\mu}(\bm{R})-\bm{R}S_{\mu^{\prime},\mu}(\bm{R})$

and only the upper or lower triangle for $\bm{R}=\bm{0}$. If needed, the
number of matrix elements may be further restricted such that only atoms in
the asymmetric unit are considered in, say, the bra. The remaining entries can
then be reconstructed by employing the (spinless) transformation properties of
the real solid harmonics
$\hat{g}^{-1}X_{l,m}(\bm{r})=X_{l,m}(g\bm{r})=\sum_{m^{\prime}}\pazocal{D}^{l}_{m,m^{\prime}}(g)X_{l,m^{\prime}}(\bm{r})$,
and the position operator
$\hat{g}^{-1}\hat{r}^{b}\hat{g}=\sum_{b^{\prime}}\pazocal{D}^{1}_{b,b^{\prime}}(g)\hat{r}^{b^{\prime}}$,
where $g\in O(3)$ and $\pazocal{D}^{l}$ is the representation of $O(3)$ of
angular momentum $l$ 44 ($l=1$ for $\bm{r}$). The result is

$\displaystyle
S_{(\mathbbm{g}a,\lambda,l,m),(\mathbbm{g}^{\prime}a^{\prime},\lambda^{\prime},l^{\prime},m^{\prime})}(\bm{R})=\sum_{m_{1},m_{1}^{\prime}}\pazocal{D}^{l}_{m,m_{1}}(g)^{*}\pazocal{D}^{l^{\prime}}_{m^{\prime},m^{\prime}_{1}}(g)S_{(a,\lambda,l,m_{1}),(\mathbbm{g}^{-1}\mathbbm{g}^{\prime}a^{\prime},\lambda^{\prime},l^{\prime},m_{1}^{\prime})}(g^{-1}\bm{R})\>,$
(10) $\displaystyle
r^{b}_{(\mathbbm{g}a,\lambda,l,m),(\mathbbm{g}^{\prime}a^{\prime},\lambda^{\prime},l^{\prime},m^{\prime})}(\bm{R})=\sum_{m_{1},m_{1}^{\prime}}\pazocal{D}^{l}_{m,m_{1}}(g)^{*}\pazocal{D}^{l^{\prime}}_{m^{\prime},m^{\prime}_{1}}(g)\sum_{b^{\prime}}\cdot$
$\displaystyle\left[\pazocal{D}^{1}_{b,b^{\prime}}(g)r^{b^{\prime}}_{(a,\lambda,l,m_{1}),(\mathbbm{g}^{-1}\mathbbm{g}^{\prime}a^{\prime},\lambda^{\prime},l^{\prime},m_{1}^{\prime})}(g^{-1}\bm{R})+\pazocal{D}^{1}_{b,b^{\prime}}(g^{-1})t^{b^{\prime}}S_{(a,\lambda,l,m_{1}),(\mathbbm{g}^{-1}\mathbbm{g}^{\prime}a^{\prime},\lambda^{\prime},l^{\prime},m_{1}^{\prime})}(g^{-1}\bm{R})\right]$

where $\mathbbm{g}=(g|\bm{t})$ is a general non-symmorphic operation in the
crystallographic point group $F$ (or space group excluding lattice
translations, $G/T$), and $\mathbbm{g}a$ represents the atom located at
$\mathbbm{g}\bm{d}_{a}$. Note that the atom
$\mathbbm{g}^{-1}\mathbbm{g}^{\prime}a^{\prime}$ may require a non-trivial
lattice vector in order to be mapped to the unit cell, thus altering
$g^{-1}\bm{R}$. If SOC is considered, then
$\hat{g}^{-1}X_{l,m,m_{s}}(\bm{r})=\sum_{m^{\prime},m^{\prime}_{s}}\pazocal{D}^{l}_{m,m^{\prime}}(g)\pazocal{D}^{1/2}_{m_{s},m^{\prime}_{s}}(g^{-1})X_{l,m^{\prime},m^{\prime}_{s}}$,
where
$\pazocal{D}^{1/2}(c_{\theta})=\pazocal{D}^{1/2}(ic_{\theta})=e^{-i(\bm{\sigma}\bm{e})\theta/2}$
(in the $\set{\uparrow,\downarrow}$ basis order) is the projective
representation of $O(3)$ of angular momentum $1/2$ which is even under
inversion $i$; $\bm{\sigma}$ being the Pauli vector and $\bm{e}$ the
counterclockwise rotation axis. The previous relations would be modified in
consequence.

In plane wave schemes, $\hat{\bm{v}}$ is sometimes replaced by the momentum
$\hat{\bm{p}}=-i\hbar\hat{\bm{\nabla}}$. This substitution is not exact in HF
or hybrid DFT-HF schemes because $\hat{\bm{r}}$ and the Fock operator do not
commute, or likewise when employing non-local pseudopotentials or including
relativistic effect such as SOC. While the deviations in the final quantities
are often not large 37, the increase in computational cost in the Gaussian
scheme is marginal enough to advise against the use of this approximation in
general, except perhaps for extremely large unit cells. Regardless, the
relevant integrals would be computed as

$\displaystyle\int_{\mathbb{R}}(x-x_{0})^{n_{x}}e^{-\alpha(x-x_{0})^{2}}\frac{\partial}{\partial
x}(x-x^{\prime}_{0})^{n^{\prime}_{x}}e^{-\alpha^{\prime}(x-x^{\prime}_{0})^{2}}dx=$
$\displaystyle
n^{\prime}_{x}I(n_{x},n^{\prime}_{x}-1,x_{0}-x^{\prime}_{0},0,\alpha,\alpha^{\prime})-2\alpha^{\prime}I(n_{x},n^{\prime}_{x}+1,x_{0}-x^{\prime}_{0},0,\alpha,\alpha^{\prime})$

###  2.4 Velocity gauge

Alternatively, the velocity gauge can be employed by imposing the minimal
coupling
$\hat{H}\left(-i\hbar\bm{\nabla}\right)\rightarrow\hat{H}\left(-i\hbar\bm{\nabla}+\absolutevalue{e}\bm{\pazocal{A}}\right)$,
where $\bm{\pazocal{A}}$ is the vector potential. An analogous derivation in
perturbation theory then yields 8, 23, 33, 31

$\sigma^{a;bc}_{\text{total}}(\omega)=-\frac{g_{s}\absolutevalue{e}^{3}}{2\hbar^{2}\omega^{2}V}\sum_{\bm{k}\in\text{BZ}}\sum_{m,n,l}\left[\frac{f_{m,n}v^{b}_{n,m}}{\omega_{m,n}-\omega+i\varepsilon}\left(\frac{v^{a}_{m,l}v^{c}_{l,n}}{\omega_{m,l}+i\varepsilon}-\frac{v^{c}_{m,l}v^{a}_{l,n}}{\omega_{l,n}+i\varepsilon}\right)-(b\leftrightarrow
c)^{*}\right]$

Both the shift and the injection conductivities are encoded in this formula.
The former can be obtained by taking only the imaginary part of the product of
complex denominators in the $\varepsilon\to 0^{+}$ limit, which is equivalent
to taking $\real\sigma^{a;bc}_{\text{total}}$ under $\hat{\pazocal{T}}$
symmetry (i.e., considering linear polarization)

$\sigma^{a;bc}_{\text{shift}}(\omega)=\frac{g_{s}\pi\absolutevalue{e}^{3}}{2\hbar^{2}V}\sum_{\bm{k}\in\text{BZ}}\sum_{m,n}\frac{f_{m,n}}{\omega^{2}_{m,n}}\imaginary\left[v^{b}_{m,n}\sum_{l\neq
m,n}\left(\frac{v^{a}_{n,l}v^{c}_{l,m}}{\omega_{n,l}}-\frac{v^{c}_{n,l}v^{a}_{l,m}}{\omega_{l,m}}\right)+(b\leftrightarrow
c)\right]\delta(\omega_{m,n}-\omega)$ (11)

In contrast with (2) and (8), the evaluation of (11) avoids the numerical
derivatives at the cost of a sum over all states $l$ that are external to the
direct optical transition. Resulting from the completeness relation
$\sum_{n}\ket{\psi_{n,\bm{k}}}\bra{\psi_{n,\bm{k}}}=\hat{I}$, in principle it
must span the whole set of bands (which with localized bases is seldom
demanding, computationally) even if they are not well represented above a
certain window from the Fermi level, but the sum should nevertheless converge
to the correct result when a sufficiently large basis is employed. While the
length gauge explicitly involves only the pair of bands corresponding to the
direct optical transition at the field frequency, a large basis should still
be needed to properly reproduce the conduction bands involved. For grids of
equal size, the evaluation of (11) is more straightforward and less
computationally demanding than (8), however, the assumption of completeness
(which is not strictly true in finite bases) makes the length gauge approach
the most reliable in general.

A similar expression to (11) that is frequently employed in the literature is
obtained within the length gauge by employing the following sum rule for the
GD

$A^{b}_{n,m;a}\rightarrow-\frac{1}{\omega_{n,m}}\left[\left(A^{b}_{n,m}\Delta^{a}_{n,m}+A^{a}_{n,m}\Delta^{b}_{n,m}\right)+\sum_{l\neq
m,n}\left(v^{b}_{n,l}A^{a}_{l,m}-A^{a}_{n,l}v^{b}_{l,m}\right)+iw^{a,b}_{n,m}\right]$
(12)

where
$w^{a,b}_{n,m}\equiv\bra{\psi_{n,\bm{k}}}[\hat{r}^{a},\hat{v}^{b}]\ket{\psi_{m,\bm{k}}}=w^{b,a}_{n,m}$
and
$\bm{\Delta}_{n,m}(\bm{k})\equiv\bm{v}_{n,n}(\bm{k})-\bm{v}_{m,m}(\bm{k})=\bm{\nabla}_{k}\omega_{n,m}(\bm{k})$
444The last identity results from the expansion of
$\bm{\nabla}_{\bm{k}}[\bra{u_{n,\bm{k}}}e^{i\bm{k}\bm{r}}\hat{H}e^{i\bm{k}\bm{r}}\ket{u_{n,\bm{k}}}]=\bm{\nabla}_{\bm{k}}E_{n,\bm{k}}$.
This sum rule can be obtained by expanding
$\partial_{k^{a}}\partial_{k^{b}}[\bra{u_{m,\bm{k}}}e^{-i\bm{k}\bm{r}}\hat{H}e^{i\bm{k}\bm{r}}\ket{u_{n,\bm{k}}}]=0$
and inserting the completeness relation above. Employing (12) in (2) yields

$\displaystyle\sigma^{a;bc}_{\text{shift}}(\omega)\rightarrow-\frac{ig_{s}\pi\absolutevalue{e}^{3}}{2\hbar^{2}V}\sum_{\bm{k}\in\text{BZ}}\sum_{m,n}\frac{f_{m,n}}{\omega^{2}_{m,n}}\cdot$
(13) $\displaystyle\left\\{v^{b}_{m,n}\left[\sum_{l\neq
m,n}\left(\frac{v^{a}_{n,l}v^{c}_{l,m}}{\omega_{n,l}}-\frac{v^{c}_{n,l}v^{a}_{l,m}}{\omega_{l,m}}\right)-\frac{v^{a}_{n,m}\Delta^{c}_{n,m}}{\omega_{n,m}}-iw^{a,c}_{n,m}\right]-(b\leftrightarrow
c)^{*}\right\\}\delta(\omega_{m,n}-\omega)$

where we note that the $\Delta^{a}$ terms cancel out. Clearly, this expression
coincides with (11) in the presence of $\hat{\pazocal{T}}$ symmetry, except
for the last two terms inside the square brackets. The term $w^{a,c}$, which
would clearly vanish in the absence of non-local terms in the Hamiltonian
($\hat{\bm{v}}=\hat{\bm{p}}$), is often computed in tight-binding or Wannier
schemes by differentiating the Hamiltonian matrix 18, 37 but its calculation
in pure DFT is non-trivial and often ignored, leaving the two-band
$\bm{\Delta}$ terms as the only difference in practice between the velocity
gauge (11) and length gauge with sume rule (LABEL:SCts) expressions. While the
contribution of these terms is usually small, in Section 2.5 we show that the
proper agreement of the length gauge formula (2) or (8) is with the velocity
gauge expression (11), at least if one neglects $w^{a,c}$.

###  2.5 First-principles results and discussion

In Figure 1 we show the shift conductivity computed with large-sized Gaussian
basis sets for some representative non-magnetic materials, in both the length
(8) and velocity (11) gauges. The self-consistent electronic structure problem
has been solved with the CRYSTAL code 45, 46, from which the
$H_{\mu,\mu^{\prime}}(\bm{R})$ are readily obtained. The input files for the
self-consistent calculations can be found in the Supporting Information, in
addition to the resulting band structures.

The starting points for the basis sets were the following: in MoS2, def2-QZVP
47 for S and pob-TZVP-rev2 48 for Mo; in GeS, def2-QZVP; in GaAs, m-pVDZ-PP-
Heyd 49; in BaTiO3, def2-QZVP with pseudo-potential (PP) from pob-TZVP-rev2 50
for Ba. In all cases the bases were modified in order to enable (or preserve)
the convergence and obtain a sufficiently accurate band structure for the
conduction bands in the energy ranges displayed in Figure 1, except for GaAs
which already presented a good dispersion with the unmodified Heyd basis. The
standard GGA PBE functional 51 was used in MoS2, GeS and BaTiO3 in order to
facilitate contrasting with the literature, while the short-range corrected
hybrid HSE06 functional 52 was employed in GaAs for the same reason. In the
latter case, the use of a hybrid functional allows to obtain the experimental
band gap at $\Gamma$ of $\sim 1.5\>\text{eV}$ avoiding the use of a scissor
correction or a $GW$ calculation.

The initial grids for the conductivity contained $2000\times 2000$ points in
the BZ for MoS2, $1500\times 1500$ for GeS, $400\times 400\times 400$ for GaAs
and $200\times 200\times 200$ for BaTiO3 555The grid choices were here
influented by benchmarking purposes, and substantially coarser ones will often
yield good results.; and were subsequently restricted to the IBZ as explained
in Section 3, in particular employing (19) due to the absence of magnetism and
the expressions from the list for the corresponding space groups: 187
($D_{3h}$) reduced to 2D for MoS2, 31 ($C_{2v}$) reduced to 2D with $c_{2,x}$
contained in the lattice plane for GeS, 216 ($T_{d}$) for GaAs and 99
($C_{4v}$) for BaTiO3. The delta function in the
$\sigma_{\text{shift}}^{a;bc}$ expressions has been approximated by a narrow
normal distribution
$\delta(x)\sim\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{x^{2}}{2\sigma^{2}}\right]$
with standard deviation $\sigma=20\>\text{meV}$ in all cases. The absence of
(unphysical) rapid fluctuations in the curves indicates that this value is not
small for the chosen $\bm{k}-$grids. The numerical derivatives in the length
gauge expression (8) have been symmetrized as
$\partial_{k^{a}}f(\bm{k})\sim\frac{1}{2h}(f(\bm{k}+h\bm{e}^{a})-f(\bm{k}-h\bm{e}^{a}))$,
with $h=\frac{10^{-7}}{a_{0}}$ in all cases. The direct lattice summations in
(LABEL:HkSk) and (5) have been truncated to the first (by length) 179 vectors
in MoS2, 120 in GeS, 179 in GaAs and 260 in BaTiO3.

The agreement with the results in the literature is generally good 23, 19, 36,
37, 35, 53, specially taking into account that moderate discrepancies can be
found commonly due to the high sensibility of $\sigma_{\text{shift}}$ to the
lattice parameters, atomic coordinates and electronic eigenfunctions 18, 54,
26, in conjunction with the different calculation methods (as outlined in the
previous sections) and convergence parameters such as the BZ grid or the
broadening of the delta functions. The proper description of the eigenstates
in the $\hbar\omega$ energy range is critical to obtain satisfactory results,
hence in principle the largest possible basis set allowed by convergence
should be used, typically a reduced QZVP or augmented TZVP. Nevertheless, in
materials with not particularly delocalized empty conduction states a smaller,
but properly calibrated basis could suffice. It should be kept in mind that
the truncation of the $\bm{R}-$sums may need to be loosed if smaller Gaussian
exponents are introduced.

It can be observed in Figure 1 that the results of both length and velocity
gauges are virtually identical in all cases, even if the external sum in (11)
is truncated in practice by the finiteness of the basis. This is in agreement
with the results of Reference 34 for third-order calculations in graphene.
However, the inclusion of the two-band $\bm{\Delta}$ terms in (LABEL:SCts) by
the sum rule induces a small, but noticeable discrepancy in some of the
$\sigma^{a;bc}_{\text{shift}}$ components, specifically $\sigma^{x;yy}$,
$\sigma^{y;yx}$ in GeS and $\sigma^{z;xx}$ in BaTiO3, while all other
components (including MoS2 and GaAs) were not visibly affected. Most of the
unaltered components have a symmetry reason to remain so: the $\bm{\Delta}$
terms clearly do not contribute to $a=b=c$ components in (LABEL:SCts), nor to
$b=c=z$ components in 2D materials since $\Delta^{z}=0$. Furthermore, it can
also be seen from (LABEL:SCts) that a sum over even permutations of $(a,b,c)$
cancels the $\bm{\Delta}$ terms, which is precisely underlying in GaAs as can
be seen from the folded summation formula for $T_{d}$ in Section 3. Only the
$\sigma^{x;xz}$ component in BaTiO3 cannot be explained by symmetry reasons,
albeit in this case only a term of the form
$v^{z}_{n,m}v^{x}_{m,n}\Delta^{x}_{n,m}$ contributes under $\hat{\pazocal{T}}$
symmetry, which is most likely small for numerical reasons in the dispersion
along $k^{x}$.

(a)

(b)

(c)

(d)

Figure 1:  Shift conductivity tensor in (a) monolayer MoS2, (b) monolayer GeS,
(c) GaAs and (d) BaTiO3 computed in the both the length (8) (solid lines) and
velocity gauge (11) (dashed lines). All linearly independent components are
shown in each case, excluding $\sigma^{z;zx}$ in GeS for an out of plane
current. IBZ summations with time-reversal symmetry (19) have been employed.

(a)

(b)

Figure 2:  Shift conductivity tensor in (a) monolayer GeS and (b) BaTiO3
computed in the length gauge (8) (solid lines) and employing the sum rule
(LABEL:SCts) (dashed lines) excluding the $w^{a,c}$ term. IBZ summations with
time-reversal symmetry (19) have been employed.

##  3 Irreducible Brillouin Zone summation

The proper use of GTOs ensures that the crystalline eigenstates indeed
transform according to the irreducible representations of the space group $G$,
and
$\hat{\mathbbm{g}}\psi_{n,\bm{k}}(\bm{r})=\psi_{n,g\bm{k}}(\bm{r}-\bm{t}-\bm{R})=\psi_{n,\bm{k}}(\mathbbm{g}^{-1}\bm{r})$,
$\forall\mathbbm{g}=(g|\bm{t}+\bm{R})\in G$. Then, since the velocity operator
transforms under the corresponding change of coordinates as
$\hat{\mathbbm{g}}^{-1}\hat{v}^{b}\hat{\mathbbm{g}}=\sum_{b^{\prime}}\pazocal{D}^{1}_{b,b^{\prime}}(g)\hat{v}^{b^{\prime}}$
666Note that the translations commute with $\hat{H}$ as fixed parameters (see
discussion around (10) for context) and likewise for $\hat{\partial}_{k^{b}}$,
the following identities must satisfied, up to an arbitrary phase factor in
the eigenstates that has no impact in the conductivity

$v^{b}_{m,n}(g\bm{k})=\sum_{b,b^{\prime}}\pazocal{D}^{1}_{b,b^{\prime}}(g)v^{b^{\prime}}_{m,n}(\bm{k})\;,\;\partial_{(gk^{b})}=\sum_{b,b^{\prime}}\pazocal{D}^{1}_{b,b^{\prime}}(g)\partial_{k^{b^{\prime}}}\;,\;A^{b}_{m,n}(g\bm{k})=\sum_{b,b^{\prime}}\pazocal{D}^{1}_{b,b^{\prime}}(g)A^{b^{\prime}}_{m,n}(\bm{k})\;,\;\forall
g\in F$ (14)

where $F$ is the (isogonal) point group of the material, formed by
disregarding all translations of the space group. On the other hand, under the
time-reversal operation
$\hat{\pazocal{T}}\hat{v}^{b}\hat{\pazocal{T}}^{-1}=-\hat{v}^{b}$ (odd power
of time units), so that recalling from Section 2.1 the action of
$\hat{\pazocal{T}}$ on the Berry connection, if time-reversal is a symmetry
then

$v^{b}_{m,n}(-\bm{k})=-v^{b}_{n,m}(\bm{k})\;,\;A^{b}_{m,n}(-\bm{k})=A^{b}_{n,m}(\bm{k})$
(15)

where again the arbitrary phase factors have been already cancelled. We note
that (14), (15) hold irrespective of whether SOC is included, since the space
group is the same 777The space group can be regarded as unaltered if
projective representations of the little groups are employed 44, which do not
appear in this derivation regardless. Alternatively, vector (standard)
representations of enlarged little groups can be considered 55. On the other
hand, we disregard the extra spin-only operations of the spin-space group that
may appear without SOC. These may possibly introduce further restrictions, but
would not overrule the present ones in any case. and the representations of
the eigenstates are not being used (but rather of the operators), and the sign
of $\hat{\pazocal{T}}^{2}=\pm\hat{I}$ has no impact. Consequently, the
following analysis is also valid irrespective of SOC.

These transformation properties can now be exploited to reduce the summations
over the BZ in (2), (8) or (LABEL:SCts) to properly weighted sums over the IBZ
or representation domain, from which the whole BZ is reconstructed by applying
point group symmetries (possibly including $\pazocal{T}$). This ensures that
each sampled $\bm{k}$ point provides unique information, reducing the
computation time of the original BZ grid by approximately the order the of the
point group, and avoids the touching of symmetry-enforced degeneracies which
are troublesome for computing the numerical derivatives in the length gauge.
Furthermore, the identification of the finite tensor entries and the linear
dependencies between them is straightforward from the IBZ expressions.

In order to avoid exceptions with duplicated points, we hereafter assume that
the BZ sampling does not contain any high symmetry point (or, to be precise,
any point whose little group is not trivial), which is without loss of
generality in sufficiently fine grids. Several cases are distinguished
depending on time-reversal symmetry, which determines the magnetic point group
$M$.

  * •

Magnetic point group of type I: $M$ contains only unitary operations, i.e.,
$M=F$ and $\pazocal{T}$ is excluded. Then by (14) the first term in (8)
satisfies

$\displaystyle\sum_{\bm{k}\in\text{BZ}}v^{b}_{m,n}(\bm{k})\left.\frac{\partial
v^{c}_{n,m}(\bm{k}^{\prime})}{\partial k^{\prime
a}}\right|_{\bm{k}}=\sum_{g\in
F}\sum_{\bm{k}\in\text{IBZ}}v^{b}_{m,n}(g\bm{k})\left.\frac{\partial
v^{c}_{n,m}(\bm{k}^{\prime})}{\partial k^{\prime a}}\right|_{g\bm{k}}=$
$\displaystyle\sum_{\bm{k}\in\text{IBZ}}\sum_{a^{\prime},b^{\prime},c^{\prime}}\mathbb{D}^{F}_{a,a^{\prime};b,b^{\prime};c,c^{\prime}}v^{b^{\prime}}_{m,n}(\bm{k})\left.\frac{\partial
v^{c^{\prime}}_{n,m}(\bm{k}^{\prime})}{\partial k^{\prime
a^{\prime}}}\right|_{\bm{k}}$

and likewise for the other terms, also in the form of (2) or (LABEL:SCts).
Here we have defined, for any subset $L\in O(3)$,

$\mathbb{D}^{L}_{a,a^{\prime};b,b^{\prime};c,c^{\prime}}\equiv\sum_{g\in
L}\pazocal{D}^{1}_{a,a^{\prime}}(g)\pazocal{D}^{1}_{b,b^{\prime}}(g)\pazocal{D}^{1}_{c,c^{\prime}}(g)$
(16)

Then, writing the $\bm{k}-$integrand explicitly for any of the forms (2), (8)
or (LABEL:SCts),

$\sigma^{a;bc}_{\text{shift}}\equiv-\frac{i}{N_{\bm{k}}}\sum_{\bm{k}\in\text{BZ}}\left[\tilde{\sigma}^{a;bc}_{\text{shift}}(\bm{k})-\tilde{\sigma}^{a;cb}_{\text{shift}}(\bm{k})^{*}\right]$
(17)

it follows that

$\sigma^{a;bc}_{\text{shift}}=-\frac{i}{\absolutevalue{F}N_{\bm{k}}^{\text{IBZ}}}\sum_{\bm{k}\in\text{IBZ}}\sum_{a^{\prime},b^{\prime},c^{\prime}}\mathbb{D}^{F}_{a,a^{\prime};b,b^{\prime};c,c^{\prime}}\left[\tilde{\sigma}^{a^{\prime};b^{\prime}c^{\prime}}_{\text{shift}}(\bm{k})-\tilde{\sigma}^{a^{\prime};c^{\prime}b^{\prime}}_{\text{shift}}(\bm{k})^{*}\right]$
(18)

where we have used that $N_{\bm{k}}=\absolutevalue{F}N_{\bm{k}}^{\text{IBZ}}$,
with $\absolutevalue{F}$ the order of $F$ and $N_{\bm{k}}^{\text{IBZ}}$ the
number of grid points in the IBZ. It is easy to see that
$\mathbb{D}^{F}_{a,a^{\prime};b,b^{\prime};c,c^{\prime}}=0$ for any group with
inversion symmetry $i\in F$, hence $\sigma^{a;bc}_{\text{shift}}=0$ in
agreement with the obvious requirement from perturbation theory.

  * •

Magnetic point group of type II: time-reversal is a symmetry by itself, i.e.,
$M=F+\pazocal{T}F$. This precludes any permanent magnetic ordering. Then by
(14) and (15)

$\displaystyle\sum_{\bm{k}\in\text{BZ}}v^{b}_{m,n}(\bm{k})\left.\frac{\partial
v^{c}_{n,m}(\bm{k}^{\prime})}{\partial k^{\prime
a}}\right|_{\bm{k}}=\frac{1}{g_{i}}\sum_{g\in
F}\sum_{\bm{k}\in\text{IBZ}}\left[v^{b}_{m,n}(g\bm{k})\left.\frac{\partial
v^{c}_{n,m}(\bm{k}^{\prime})}{\partial k^{\prime
a}}\right|_{g\bm{k}}+v^{b}_{m,n}(-g\bm{k})\left.\frac{\partial
v^{c}_{n,m}(\bm{k}^{\prime})}{\partial k^{\prime
a}}\right|_{-g\bm{k}}\right]=$
$\displaystyle\frac{1}{g_{i}}\sum_{\bm{k}\in\text{IBZ}}\sum_{a^{\prime},b^{\prime},c^{\prime}}\mathbb{D}^{F}_{a,a^{\prime};b,b^{\prime};c,c^{\prime}}2i\imaginary\left[v^{b^{\prime}}_{m,n}(\bm{k})\left.\frac{\partial
v^{c^{\prime}}_{n,m}(\bm{k}^{\prime})}{\partial k^{\prime
a^{\prime}}}\right|_{\bm{k}}\>\right]$

where $g_{i}=2$ if $i\in F$ or in the 2D cases with out of plane $c_{2,z}$
rotational symmetry, and $g_{i}=1$ otherwise. The group $F$ is again forced to
be non-centrosymmetric, but time-reversal symmetry effectively halves the IBZ
with respect to the type I case by introducing a relation between the
$\pm\bm{k}$ pairs (except in the cases where $g_{i}=2$). The general shift
current tensor (17) thus satisfies

$\sigma^{a;bc}_{\text{shift}}=\frac{1}{\absolutevalue{F}N^{\text{IBZ}}_{\bm{k}}}\sum_{\bm{k}\in\text{IBZ}}\sum_{a^{\prime},b^{\prime},c^{\prime}}\mathbb{D}^{F}_{a,a^{\prime};b,b^{\prime};c,c^{\prime}}\imaginary\left[\tilde{\sigma}^{a^{\prime};b^{\prime}c^{\prime}}_{\text{shift}}(\bm{k})+\tilde{\sigma}^{a^{\prime};c^{\prime}b^{\prime}}_{\text{shift}}(\bm{k})\right]$
(19)

where $N_{\bm{k}}=2\absolutevalue{F}N_{\bm{k}}^{\text{IBZ}}/g_{i}$. In
agreement with Section 2.1, $\sigma^{a;bc}_{\text{shift}}$ is real in this
case.

  * •

Magnetic point group of type III: exactly half of the unitary operations are
paired with time-reversal, i.e., $M=H+\pazocal{T}(F-H)$ where $H\subset F$ is
an ordinary point group and $F-H$ is thus not a group. In this case (14) is
valid for $g\in H$ whereas (15) must be used in combination with (14) for
$g\in F-H$, and the IBZ is defined by $F$ as in type I. Hence, symmetrizing in
$\pm\bm{k}$,

$\displaystyle\sum_{\bm{k}\in\text{BZ}}v^{b}_{m,n}(\bm{k})\left.\frac{\partial
v^{c}_{n,m}(\bm{k}^{\prime})}{\partial k^{\prime
a}}\right|_{\bm{k}}=\frac{1}{2}\sum_{\bm{k}\in\text{IBZ}}\sum_{\bm{k}^{\prime}=\pm\bm{k}}\left[\sum_{g\in
H}v^{b}_{m,n}(g\bm{k}^{\prime})\left.\frac{\partial
v^{c}_{n,m}(\bm{k}^{\prime\prime})}{\partial k^{\prime\prime
a}}\right|_{g\bm{k}^{\prime}}+\right.$ $\displaystyle\left.\sum_{\tilde{g}\in
F-H}v^{b}_{m,n}(-\tilde{g}\bm{k}^{\prime})\left.\frac{\partial
v^{c}_{n,m}(\bm{k}^{\prime\prime})}{\partial k^{\prime\prime
a}}\right|_{-\tilde{g}\bm{k}^{\prime}}\right]=$
$\displaystyle\sum_{\bm{k}\in\text{IBZ}}\sum_{a^{\prime},b^{\prime},c^{\prime}}\left[\mathbb{D}^{H}_{a,a^{\prime};b,b^{\prime};c,c^{\prime}}v^{b^{\prime}}_{m,n}(\bm{k})\left.\frac{\partial
v^{c^{\prime}}_{n,m}(\bm{k}^{\prime})}{\partial k^{\prime
a^{\prime}}}\right|_{\bm{k}}-\mathbb{D}^{F-H}_{a,a^{\prime};b,b^{\prime};c,c^{\prime}}v^{b^{\prime}}_{n,m}(\bm{k})\left.\frac{\partial
v^{c^{\prime}}_{m,n}(\bm{k}^{\prime})}{\partial k^{\prime
a^{\prime}}}\right|_{\bm{k}}\>\right]$

and, noting that $\mathbb{D}^{F-H}=\mathbb{D}^{F}-\mathbb{D}^{H}$, (17) can be
expressed as

$\displaystyle\sigma^{a;bc}_{\text{shift}}=-\frac{i}{\absolutevalue{F}N^{\text{IBZ}}_{\bm{k}}}\sum_{\bm{k}\in\text{IBZ}}\sum_{a^{\prime},b^{\prime},c^{\prime}}$
(20)
$\displaystyle\left[2\mathbb{D}^{H}_{a,a^{\prime};b,b^{\prime};c,c^{\prime}}\real\left(\tilde{\sigma}^{a^{\prime};b^{\prime}c^{\prime}}_{\text{shift}}(\bm{k})-\tilde{\sigma}^{a^{\prime};c^{\prime}b^{\prime}}_{\text{shift}}(\bm{k})\right)-\mathbb{D}^{F}_{a,a^{\prime};b,b^{\prime};c,c^{\prime}}\left(\tilde{\sigma}^{a^{\prime};b^{\prime}c^{\prime}}_{\text{shift}}(\bm{k})^{*}-\tilde{\sigma}^{a^{\prime};c^{\prime}b^{\prime}}_{\text{shift}}(\bm{k})\right)\right]$

If the crystalline structure is centrosymmetric disregarding magnetism, there
are three cases upon magnetization. The first one is that $i$ is completely
removed by the magnetization, $i\notin F$ and $M$ is of type I or III, thus it
imposes no restrictions on $\sigma^{a;bc}_{\text{shift}}$. The second one is
$i\in H\subset F$, which implies that $\mathbb{D}^{H}=\mathbb{D}^{F}=0$ and
$\sigma^{a;bc}_{\text{shift}}=0$. The third one is $i\in F-H$, which implies
that $\mathbb{D}^{F}=0$ and $\real\sigma^{a;bc}_{\text{shift}}=0$. Therefore,
a centrosymmetric material may host a finite shift current upon magnetization
as long as inversion is not a symmetry by itself but in combination with time-
reversal (often termed $\pazocal{P}\pazocal{T}$ symmetry), a situation that is
most common in antiferromagnets and which by (1) and (20) yields a shift
current under circular (or elliptical) polarization.

Note that (18), (19), (20) also include the case of 2D materials with the
global symmetry $\sigma_{h}\in F$ (mirror plane parallel to the lattice),
since the required $1/2$ factor in the
$\sum_{\bm{k}\in\text{BZ}}\to(1/2)\sum_{g\in F}\sum_{\bm{k}\in\text{IBZ}}$
substitution is absorbed in
$\absolutevalue{F/C_{1h}}=N_{\bm{k}}/N_{\bm{k}}^{\text{IBZ}}$.

Therefore, for a given material one must evaluate (18), (19) or (20) according
to its magnetic point group, which requires computing (16) for the appropriate
subgroup of $O(3)$, and parametrizing the IBZ according to its space group
since the lattice type determines the BZ. In the list at the end of this
section we show, for each (unitary) non-centrosymmetric point group, the
values of (16) in the simplified notation
$aa^{\prime};bb^{\prime};cc^{\prime}$, excluding the vanishing components as
well as those that are redundant due to the permutation properties
$aa^{\prime};bb^{\prime};cc^{\prime}=aa^{\prime};cc^{\prime};bb^{\prime}=bb^{\prime};aa^{\prime};cc^{\prime}=bb^{\prime};cc^{\prime};aa^{\prime}=cc^{\prime};aa^{\prime};bb^{\prime}=cc^{\prime};bb^{\prime};aa^{\prime}$
and the reciprocity relation
$a^{\prime}a;b^{\prime}b;c^{\prime}c=aa^{\prime};bb^{\prime};cc^{\prime}$.
Afterwards the finite components of the $\sigma^{a;bc}_{\text{shift}}$ tensor
and the relations between them are displayed, omitting the persistent
$\sigma^{a;bc}_{\text{shift}}=(\sigma^{a;cb}_{\text{shift}})^{*}$ and
specifying which of (18) or (19) are finite if only one of them is 888For
example, $\imaginary\sigma^{x;yz}_{I}$ indicates that (18) is purely imaginary
and (19) is null for the $a=x,b=y,c=z$ component. Note that
$\sigma^{a;bc}_{I}=\sigma^{a;bc}_{II}$ implies that (18) is real. We omit the
mandatory equalities between (18) and (19) for $b=c$ components.. When
multiple orientations along the Cartesian axes are possible, we indicate the
orientation of the generating operations with the notation $c_{n,a}$ for a
counterclockwise $2\pi/n-$fold rotation with axis along $a$, $\sigma_{a}$ for
a reflection through a plane perpendicular to $a$ and $s_{n,a}$ for the
improper rotation $c_{n,a}\sigma_{a}$. If a different orientation is employed,
then there exists a transformation $r\in O(3)$ of our coordinate system that
matches the chosen orientation, and the new conductivity tensor
$\overline{\sigma}^{a;bc}$ can be computed from the one given here as
$\overline{\sigma}^{a;bc}=\sum_{a^{\prime},b^{\prime},c^{\prime}}\sigma^{a^{\prime};b^{\prime},c^{\prime}}$,
which will often result in a simple permutation of coordinates.

On the other hand, the grid of the IBZ can be obtained (non-uniquely) as a
restriction of a larger grid by imposing a set of linear constraints on the
coefficients $\alpha_{i}$ of the
$\bm{k}=\sum_{i=1}^{\text{dim}}\alpha_{i}\bm{G}_{i}$ points (“dim” being the
dimension of the lattice) in the basis of reciprocal lattice vectors
$\bm{G}_{i}$. Specifically,

$\text{IBZ}=\hat{R}\set{\sum_{i=1}^{\text{dim}}\alpha_{i}\bm{G}_{i}\>}{\>UA^{t}\bm{\alpha}<u}$
(21)

where $U$ and $u$ depend on the space group and are given in the list below
for a particular choice of reciprocal lattice vectors
$\bm{G}^{\text{ref}}_{i}$, which are those of Reference 55 (see tables 3.1 and
3.3 therein) and which we indicate in the list in $2\pi$ units for each
lattice type under the first point group for which they appear. The
coefficients $\alpha_{i}$ should initially be spanning the $(-1,1)$ interval
because, while a range of length 1 for each would suffice, it is not
guaranteed that our specific parametrization yields the same range for each
coefficient 999Note that in this case a $N_{1}\times N_{2}\times N_{3}$ grid
for the $\alpha_{i}$ coefficients corresponds to a
$(N_{1}/2)\times(N_{2}/2)\times(N_{3}/2)$ grid of the BZ. In (21), $A$ and $R$
are introduced to allow for different sets of lattice vectors: $A$ is the
transformation matrix relating both sets of reciprocal vectors as

$\begin{pmatrix}\bm{G}_{1}\\\ \vdots\\\
\bm{G}_{\text{dim}}\end{pmatrix}=A\begin{pmatrix}\bm{G}^{\text{ref}}_{1}\\\
\vdots\\\ \bm{G}^{\text{ref}}_{\text{dim}}\end{pmatrix}$

and $\hat{R}$ is the operator represented by the $R\in SO(\text{dim})$
rotation matrix that relates both BZs, acting on the set of $\bm{k}$ points to
its right. The parametrized IBZs lie within the BZ except for the triclinic
and monoclinic systems. The space groups $G$ are labelled by their
international number, and the different variations of IBZs that cannot be
obtained through rotations are considered.

The list can also be employed with 2D materials in the following way: first
identify the 3D space group that is compatible with the two-dimensional
structure, taking into account that in this case only groups with
$\bm{G}^{\text{ref}}_{3}\parallel z$ and lacking non-symmorphic translations
along $z$ can be compatible. Then, if $F=C_{2}$ or $C_{2v}$, permute two
columns of $U$ to match the $c_{2}$ axis orientation 101010For example, in the
GeS $C_{2v}$ ($c_{2x}$) geometry of the present work, the first and third
columns are permuted., if $F=C_{s}$ ($\sigma_{y}$) do the columns permutation
(123)$\to$(132), if $F=C_{s}$ ($\sigma_{x}$) do the columns permutation
(123)$\to$(231); and if any of these 3 permutations was performed, do the same
permutation of the $(x,y,z)$ indices in $aa^{\prime};bb^{\prime};cc^{\prime}$
and $\sigma^{a;bc}$. Finally, eliminate all rows in $U$ and $u$ whose first
two entries in $U$ are zero. The reduction of IBZs with $\hat{\pazocal{T}}$
has been chosen such that the last step yields the correct IBZ in 2D
(invariant with $c_{2,z}$ or halved without it) with all compatible space
groups.

Either the red or the blue rows in $U$ and $u$ are included for a given
system, in particular the red rows (excluding the blue) are included when the
IBZ is not affected by $\pazocal{T}$, and the blue rows (excluding the red)
when the IBZ is halved by $\pazocal{T}$. This can be determined for the 4
types of magnetic space groups 55, which we denote as $MG$, as an application
of the previous discussion of magnetic point groups, in addition to the
specific formula for $\sigma^{a;bc}_{\text{shift}}$.

  * •

Magnetic space group of type I: $MG=G$ $\rightarrow$ red rows, equation (18).

  * •

Magnetic space group of type II: $MG=G+\pazocal{T}G$ $\rightarrow$ blue rows,
equation (19).

  * •

Magnetic space group of type III: $MG=\tilde{G}+\pazocal{T}(G-\tilde{G})$,
where $\tilde{G}$ is an space group whose point group has half the order of
$G$’s point group, the latter of which determines the IBZ $\rightarrow$ red
rows, equation (20).

  * •

Magnetic space group of type IV: $MG=G+\pazocal{T}(e|\bm{t}_{0})G$, where
$(e|\bm{t}_{0})\notin G$ is a pure translation. The IBZ is defined by $G$, for
our purpose with $\pazocal{T}$ since the extra translation, as in non-
symmorphic groups, is inconsequential $\rightarrow$ blue rows, equation (19).

  *   * •

$F=C_{1}$

$aa^{\prime};bb^{\prime};cc^{\prime}=\delta_{a,a^{\prime}}\delta_{b,b^{\prime}}\delta_{c,c^{\prime}}$
$\text{No restrictions on }\sigma^{a;bc}$

    * $\circ$

$G=1$ $[\Gamma_{t}]$

$\bm{G}_{i}\text{ arbitrary, }(A=R=I),\;U=\begin{pmatrix}-1&0&0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
0&-1&0\\\ 0&1&0\\\ 0&0&-1\\\ 0&0&1\end{pmatrix},\;u=\begin{pmatrix}0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1/2}\\\ 0\\\
1\\\ 0\\\ 1\end{pmatrix}$

  * •

$F=C_{2}$ ($c_{2,z}$)

$zz;zz;zz=2,\;zz;xx;xx=2,\;zz;yy;yy=2,\;xx;yy;zz=2$
$\sigma^{z;zz},\;\sigma^{z;xx},\;\sigma^{z;yy},\;\sigma^{x;xz},\;\sigma^{y;yz},\;\sigma^{x;yz},\;\sigma^{y;zx},\;\sigma^{z;xy}$

    * $\circ$

$G=3-4$ [$\Gamma_{m}$]

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}-1/(b\cdot\tan\gamma)&-1/b&0\\\
1/(a\cdot\sin\gamma)&0&0\\\ 0&0&1/c\end{pmatrix},\;U=\begin{pmatrix}-1&0&0\\\
1&0&0\\\ 0&-1&0\\\ 0&1&0\\\ 0&0&-1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\end{pmatrix},\;u=\begin{pmatrix}0\\\
1/2\\\ 0\\\ 1\\\ 0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1/2}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\end{pmatrix}$

    * $\circ$

$G=5$ [$\Gamma_{m}^{b}$]

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}-1/(b\cdot\tan\gamma)&-1/b&0\\\
1/(a\cdot\sin\gamma)&0&-1/c\\\
1/(a\cdot\sin\gamma)&0&1/c\end{pmatrix},\;\;U,u\text{ as in }G=3$

  * •

$F=C_{s}$ ($\sigma_{z}$)

$xx;xx;xx=2,\;yy;yy;yy=2,\;xx;yy;yy=2,\;yy;xx;xx=2,\;xx;zz;zz=2,\;yy;zz;zz=2$
$\sigma^{x;xx},\;\sigma^{y;yy},\;\sigma^{x;yy},\;\sigma^{y;xx},\;\sigma^{x;xy},\;\sigma^{y;yx},\;\sigma^{x;zz},\;\sigma^{y;zz},\;\sigma^{z;zx},\;\sigma^{z;zy}$

    * $\circ$

$G=6-7$ [$\Gamma_{m}$]

$U=\begin{pmatrix}-1&0&0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
0&-1&0\\\ 0&1&0\\\ 0&0&-1\\\ 0&0&1\end{pmatrix},\;u=\begin{pmatrix}0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1/2}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\ 0\\\
1\\\ 0\\\ 1/2\end{pmatrix}$

    * $\circ$

$G=8-9$ [$\Gamma_{m}^{b}$]

$U,u\text{ as in }G=6$

  * •

$F=D_{2}$

$xx;yy;zz=4$ $\sigma^{x;yz},\;\sigma^{y;zx},\;\sigma^{z;xy}$

    * $\circ$

$G=16-19$ [$\Gamma_{o}$]

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}0&-1/b&0\\\ 1/a&0&0\\\
0&0&1/c\end{pmatrix},\;U=\begin{pmatrix}1&0&0\\\ -1&0&0\\\ 0&-1&0\\\ 0&1&0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}\\\
0&0&1\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\end{pmatrix},\;u=\begin{pmatrix}0\\\
1/2\\\ 0\\\ 1/2\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\ 1/2\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\end{pmatrix}$

    * $\circ$

$G=20-21$ [$\Gamma_{o}^{b}$]. $\eta\equiv\frac{b^{2}-a^{2}}{b^{2}+a^{2}}$

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}1/a&-1/b&0\\\ 1/a&1/b&0\\\
0&0&1/c\end{pmatrix}\text{with }a>b,\;U=\begin{pmatrix}1&1&0\\\ \eta&1&0\\\
0&0&1\\\ 1&-1&0\\\ -1&-1&0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\end{pmatrix},\;u=\begin{pmatrix}1\\\
1/2\\\ 1/2\\\ 0\\\ 0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\end{pmatrix}$

    * $\circ$

$G=22$ [$\Gamma_{o}^{f}$]

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}1/a&1/b&1/c\\\
-1/a&-1/b&1/c\\\ 1/a&-1/b&-1/c\end{pmatrix}$

      * *

Variation 1: $a^{-2}\geq b^{-2}+c^{-2}$.
$\mu_{ab}\equiv\frac{-a^{2}b^{2}+a^{2}c^{2}+b^{2}c^{2}}{a^{2}b^{2}+a^{2}c^{2}+b^{2}c^{2}}$,
$\mu_{ac}\equiv\frac{a^{2}b^{2}-a^{2}c^{2}+b^{2}c^{2}}{a^{2}b^{2}+a^{2}c^{2}+b^{2}c^{2}}$

$U=\begin{pmatrix}1&-\mu_{ab}&\mu_{ab}+\mu_{ac}-1\\\ 1&-1&-1\\\ 1&1&-1\\\
-1&1&1\\\ -1&1&-1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mu_{ab}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mu_{ac}}\end{pmatrix},\;u=\begin{pmatrix}1/2\\\
1\\\ 1\\\ 0\\\ 0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\end{pmatrix}$

      * *

Variation 2: $a^{-2}<b^{-2}+c^{-2}$ and $b^{-2}<c^{-2}+a^{-2}$.
$\mu_{ab}\equiv\frac{-a^{2}b^{2}+a^{2}c^{2}+b^{2}c^{2}}{a^{2}b^{2}+a^{2}c^{2}+b^{2}c^{2}}$,
$\mu_{ac}\equiv\frac{a^{2}b^{2}-a^{2}c^{2}+b^{2}c^{2}}{a^{2}b^{2}+a^{2}c^{2}+b^{2}c^{2}}$

$U=\begin{pmatrix}1&-\mu_{ab}&\mu_{ab}+\mu_{ac}-1\\\ 1&-1&-1\\\ 1&1&-1\\\
1&-1&1\\\ -1&1&1\\\ -1&1&-1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mu_{ab}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mu_{ac}}\end{pmatrix},\;u=\begin{pmatrix}1/2\\\
1\\\ 1\\\ 1\\\ 0\\\ 0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\end{pmatrix}$

    * $\circ$

$G=23-24$ [$\Gamma_{o}^{v}$]. $\eta_{ab}\equiv\frac{a^{2}}{a^{2}+b^{2}}$,
$\eta_{ac}\equiv\frac{a^{2}}{a^{2}+c^{2}}$,
$\eta_{bc}\equiv\frac{b^{2}}{b^{2}+c^{2}}$

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}1/a&0&1/c\\\ 0&-1/b&1/c\\\
1/a&-1/b&0\end{pmatrix}\text{with }c=\max(a,b,c)\>,$
$U=\begin{pmatrix}1&\eta_{ac}&1-\eta_{ac}\\\
\eta_{bc}&2\eta_{bc}-1&\eta_{bc}-1\\\ 1-\eta_{ab}&-\eta_{ab}&1-2\eta_{ab}\\\
1&1&0\\\ 0&1&1\\\ -1&0&-1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-2\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-\eta_{ac}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{bc}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\eta_{bc}-1}\end{pmatrix},\;u=\begin{pmatrix}1/2\\\
1/2\\\ 1/2\\\ 1\\\ 0\\\ 0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\end{pmatrix}$

  * •

$F=C_{2v}$ ($c_{2,z}$, $\sigma_{x}$)

$zz;zz;zz=4,\;zz;xx;xx=4,\;zz;yy;yy=4$
$\sigma^{z;zz},\;\sigma^{z;xx},\;\sigma^{z;yy},\;\sigma^{x;xz},\;\sigma^{y;yz}$

    * $\circ$

$G=25-34$ [$\Gamma_{o}$] $\rightarrow$ Same as $G=16$

    * $\circ$

$G=35-41$ [$\Gamma_{o}^{b}$] $\rightarrow$ Same as $G=20$

    * $\circ$

$G=42-43$ [$\Gamma_{o}^{f}$] $\rightarrow$ Same as $G=22$

    * $\circ$

$G=44-46$ [$\Gamma_{o}^{v}$] $\rightarrow$ Same as $G=23$

  * •

$F=C_{4}$ ($c_{4,z}$)

$zz;zz;zz=4,\;zz;xx;xx=zz;yy;yy=zz;xy;xy=2,\;xx;yy;zz=-xy;yx;zz=2$
$\sigma^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy},\;\sigma^{x;xz}=\sigma^{y;yz},\;\sigma^{x;yz}=-\sigma^{y;zx},\;\imaginary\sigma^{z;xy}_{I}$

    * $\circ$

$G=75-78$ [$\Gamma_{q}$]

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}1/a&0&0\\\ 0&1/a&0\\\
0&0&1/c\end{pmatrix},\;U=\begin{pmatrix}-1&0&0\\\ 1&0&0\\\ 0&-1&0\\\ 0&1&0\\\
0&0&1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\end{pmatrix},\;u=\begin{pmatrix}0\\\
1/2\\\ 0\\\ 1/2\\\ 1/2\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\end{pmatrix}$

    * $\circ$

$G=79-80$ [$\Gamma_{q}^{v}$]

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}0&1/a&1/c\\\ 1/a&0&1/c\\\
1/a&1/a&0\end{pmatrix}$

      * *

Variation 1: $a>c$. $\eta_{ac}\equiv\frac{a^{2}}{a^{2}+c^{2}}$

$U=\begin{pmatrix}\eta_{ac}&1&1-\eta_{ac}\\\ 1&1&2\\\ -1&0&-1\\\ 0&-1&-1\\\
1&\eta_{ac}&1-\eta_{ac}\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-2\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-\eta_{ac}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-2\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-\eta_{ac}}\end{pmatrix},\;u=\begin{pmatrix}1/2\\\
1\\\ 0\\\ 0\\\ 1/2\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\end{pmatrix}$

      * *

Variation 2: $c>a$. $\eta_{ac}\equiv\frac{a^{2}}{a^{2}+c^{2}}$

$U=\begin{pmatrix}\eta_{ac}&1&1-\eta_{ac}\\\ 1&1&2\\\ 1&1&0\\\ -1&0&-1\\\
0&-1&-1\\\ 1&\eta_{ac}&1-\eta_{ac}\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-2\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-\eta_{ac}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-2\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-\eta_{ac}}\end{pmatrix},\;u=\begin{pmatrix}1/2\\\
1\\\ 1\\\ 0\\\ 0\\\ 1/2\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\end{pmatrix}$

  * •

$F=S_{4}$ ($s_{4,z}$)

$zz;xx;xx=zz;yy;yy=-zz;xy;xy=2,\;xx;yy;zz=xy;yx;zz=2$
$\sigma^{z;xx}=-\sigma^{z;yy},\;\sigma^{x;xz}=-\sigma^{y;yz},\;\sigma^{x;yz}=\sigma^{y;zx},\;\sigma^{z;xy}_{I}=\sigma^{z;xy}_{II}$

    * $\circ$

$G=81$ [$\Gamma_{q}$] $\rightarrow$ Same as $G=75$

    * $\circ$

$G=82$ [$\Gamma_{q}^{v}$] $\rightarrow$ Same as $G=79$

  * •

$F=D_{4}$ ($c_{4,z}$, $c_{2,x}$)

$xx;yy;zz=-xy;yx;zz=4,\;$
$\sigma^{x;yz}=-\sigma^{y;zx},\;\imaginary\sigma^{z;xy}_{I}$

    * $\circ$

$G=89-96$ [$\Gamma_{q}$]

$U=\begin{pmatrix}1&-1&0\\\ -1&0&0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}\\\
0&0&1\\\ 0&1&0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\end{pmatrix},\;u=\begin{pmatrix}0\\\
0\\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
1/2\\\ 1/2\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\end{pmatrix}$

    * $\circ$

$G=97-98$ [$\Gamma_{q}^{v}$]

      * *

Variation 1: $a>c$. $\eta_{ac}\equiv\frac{a^{2}}{a^{2}+c^{2}}$

$U=\begin{pmatrix}\eta_{ac}&1&1-\eta_{ac}\\\ 1&1&2\\\ 1&-1&0\\\ -1&0&-1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-2\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-\eta_{ac}}\end{pmatrix},\;u=\begin{pmatrix}1/2\\\
1\\\ 0\\\ 0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\end{pmatrix}$

      * *

Variation 2: $c>a$. $\eta_{ac}\equiv\frac{a^{2}}{a^{2}+c^{2}}$

$U=\begin{pmatrix}\eta_{ac}&1&1-\eta_{ac}\\\ 1&1&2\\\ 1&1&0\\\ 1&-1&0\\\
-1&0&-1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-2\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1-\eta_{ac}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix},\;u=\begin{pmatrix}1/2\\\
1\\\ 1\\\ 0\\\ 0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\end{pmatrix}$

  * •

$F=C_{4v}$ ($c_{4,z}$, $\sigma_{x}$)

$zz;zz;zz=8,\;zz;xx;xx=zz;yy;yy=zz;xy;xy=4$
$\sigma^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy},\;\sigma^{x;xz}=\sigma^{y;yz}$

    * $\circ$

$G=99-106$ [$\Gamma_{q}$] $\rightarrow$ Same as $G=89$

    * $\circ$

$G=107-110$ [$\Gamma_{q}^{v}$] $\rightarrow$ Same as $G=97$

  * •

$F=D_{2d}$ ($s_{4,z}$, $c_{2,x}$)

$xx;yy;zz=xy;yx;zz=4,\;$
$\sigma^{x;yz}=\sigma^{y;zx},\;\sigma^{z;xy}_{I}=\sigma^{z;xy}_{II}$

    * $\circ$

$G=111-118$ [$\Gamma_{q}$] $\rightarrow$ Same as $G=89$

    * $\circ$

$G=119-122$ [$\Gamma_{q}^{v}$] $\rightarrow$ Same as $G=97$

  * •

$F=C_{3}$ ($c_{3,z}$)

$\displaystyle
xx;xx;xx=-xx;xy;xy=3/4,\;yy;yy;yy=-yy;yx;yx=3/4,\;zz;zz;zz=3,\;$
$\displaystyle zz;xx;xx=zz;yy;yy=zz;xy;xy=3/2,\;xx;yy;zz=-xy;yx;zz=3/2$
$\displaystyle\sigma^{x;xx}=-\sigma^{x;yy}=-\sigma^{y;yx}_{I}=-\sigma^{y;yx}_{II},\;\sigma^{y;yy}=-\sigma^{y;xx}=-\sigma^{x;xy}_{I}=-\sigma^{x;xy}_{II},\;\sigma^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy},\;$
$\displaystyle\sigma^{x;xz}=\sigma^{y;yz},\;\sigma^{x;yz}=-\sigma^{y;zx},\;\imaginary\sigma^{z;xy}_{I}$

    * $\circ$

$G=143-145$ [$\Gamma_{h}$]

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}1/\sqrt{3}a&-1/a&0\\\
2/\sqrt{3}a&0&0\\\ 0&0&1/c\end{pmatrix},\;U=\begin{pmatrix}1&0&0\\\ 0&0&-1\\\
0&0&1\\\ 1&2&0\\\ -1&1&0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-2}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix},\;u=\begin{pmatrix}0\\\
1/2\\\ 1/2\\\ 1\\\ 1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\end{pmatrix}$

    * $\circ$

$G=146$ [$\Gamma_{rh}$]

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}0&-2/3a&1/3c\\\
1/\sqrt{3}a&1/3a&1/3c\\\ -1/\sqrt{3}a&1/3a&1/3c\end{pmatrix}$

      * *

Variation 1: $a>\sqrt{2}c$. $\eta_{ac}\equiv\frac{a^{2}-2c^{2}}{a^{2}+c^{2}}$,
$\eta_{a4c}\equiv\frac{2a^{2}-4c^{2}}{a^{2}+4c^{2}}$

$U=\begin{pmatrix}\eta_{a4c}&2&\eta_{a4c}\\\ 0&1&-1\\\
-\eta_{a4c}&-\eta_{a4c}&-2\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\\\
1&-1&0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-2}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{a4c}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{a4c}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\end{pmatrix},\;u=\begin{pmatrix}1\\\
1\\\ 1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\ 0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix}$

      * *

Variation 2: $a<\sqrt{2}c$. $\eta_{ac}\equiv\frac{a^{2}-2c^{2}}{a^{2}+c^{2}}$,
$\eta_{a4c}\equiv\frac{2a^{2}-4c^{2}}{a^{2}+4c^{2}}$

$U=\begin{pmatrix}\eta_{a4c}&2&\eta_{a4c}\\\ 1&1&\eta_{ac}\\\
-\eta_{a4c}&-\eta_{a4c}&-2\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\\\
1&1&1\\\ 1&-1&0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-2}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{a4c}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{a4c}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
-1&-1&-1\\\ -1&-\eta_{ac}&-1\end{pmatrix},\;u=\begin{pmatrix}1\\\ 1\\\ 1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\ 3/2\\\
0\\\ {\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\ 3/2\\\
1\end{pmatrix}$

  * •

$F=D_{3}$ ($c_{3,z}$, $c_{2,x}$)

$xx;xx;xx=-xx;xy;xy=3/2,\;xx;yy;zz=-xy;yx;zz=3$
$\sigma^{x;xx}=-\sigma^{x;yy}=-\sigma^{y;yx}_{I}=-\sigma^{y;yx}_{II},\;\sigma^{x;yz}=-\sigma^{y;zx},\;\imaginary\sigma^{z;xy}_{I}$

    * $\circ$

$G=149-154$ [$\Gamma_{h}$]

$U=\begin{pmatrix}1&0&0\\\ 0&0&-1\\\ 0&0&1\\\ 1&2&0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-2}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix},\;u=\begin{pmatrix}0\\\
1/2\\\ 1/2\\\ 1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix}$

    * $\circ$

$G=155$ [$\Gamma_{rh}$]

      * *

Variation 1: $a>\sqrt{2}c$. $\eta_{ac}\equiv\frac{a^{2}-2c^{2}}{a^{2}+c^{2}}$,
$\eta_{a4c}\equiv\frac{2a^{2}-4c^{2}}{a^{2}+4c^{2}}$

$U=\begin{pmatrix}\eta_{a4c}&2&\eta_{a4c}\\\ 0&1&-1\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{a4c}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{a4c}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-2}\\\
-1&0&1\\\ 1&-1&0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}\end{pmatrix},\;u=\begin{pmatrix}1\\\
1\\\ {\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
0\\\ 0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\end{pmatrix}$

      * *

Variation 2: $a<\sqrt{2}c$. $\eta_{ac}\equiv\frac{a^{2}-2c^{2}}{a^{2}+c^{2}}$,
$\eta_{a4c}\equiv\frac{2a^{2}-4c^{2}}{a^{2}+4c^{2}}$

$U=\begin{pmatrix}\eta_{a4c}&2&\eta_{a4c}\\\ 1&1&\eta_{ac}\\\
-\eta_{a4c}&-\eta_{a4c}&-2\\\ -1&0&1\\\ 1&1&1\\\ 1&-1&0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\eta_{ac}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}\end{pmatrix},\;u=\begin{pmatrix}1\\\
1\\\ 1\\\ 0\\\ 3/2\\\ 0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}3/2}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\end{pmatrix}$

  * •

$F=C_{3v}$ ($c_{3,z}$, $\sigma_{y}$)

$xx;xx;xx=-xx;xy;xy=3/2,\;zz;zz;zz=6,\;zz;xx;xx=zz;yy;yy=zz;xy;xy=3$
$\sigma^{x;xx}=-\sigma^{x;yy}=-\sigma^{y;yx}_{I}=-\sigma^{y;yx}_{II},\;\sigma^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy},\;\sigma^{x;xz}=\sigma^{y;yz}$

    * $\circ$

$G=156-159$ [$\Gamma_{h}$] $\rightarrow$ Same as $G=149$

    * $\circ$

$G=160-161$ [$\Gamma_{rh}$] $\rightarrow$ Same as $G=155$

  * •

$F=C_{6}$ ($c_{6,z}$)

$zz;zz;zz=6,\;zz;xx;xx=zz;yy;yy=zz;xy;xy=3,\;xx;yy;zz=-xy;yx;zz=3$
$\sigma^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy},\;\sigma^{x;xz}=\sigma^{y;yz},\;\sigma^{x;yz}=-\sigma^{y;zx},\;\imaginary\sigma^{z;xy}$

    * $\circ$

$G=168-173$ [$\Gamma_{h}$]

$U=\begin{pmatrix}1&0&0\\\ 0&0&1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\\\
1&2&0\\\ -1&1&0\\\ -1&-1&0\end{pmatrix},\;u=\begin{pmatrix}0\\\ 1/2\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\\\ 1\\\
1\\\ 0\end{pmatrix}$

  * •

$F=C_{3h}$ ($c_{3,z}$, $\sigma_{z}$)

$xx;xx;xx=-xx;xy;xy=3/2,\;yy;yy;yy=-yy;yx;yx=3/2$
$\sigma^{x;xx}=-\sigma^{x;yy}=-\sigma^{y;yx}_{I}=-\sigma^{y;yx}_{II},\;\sigma^{y;yy}=-\sigma^{y;xx}=-\sigma^{x;xy}_{I}=-\sigma^{x;xy}_{II}$

    * $\circ$

$G=174$ [$\Gamma_{h}$]

$U=\begin{pmatrix}1&0&0\\\ 0&0&-1\\\ 0&0&1\\\ 1&2&0\\\ -1&1&0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-2}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix},\;u=\begin{pmatrix}0\\\
0\\\ 1/2\\\ 1\\\ 1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\end{pmatrix}$

  * •

$F=D_{6}$ ($c_{6,z}$, $c_{2,x}$)

$xx;yy;zz=-xy;yx;zz=6$
$\sigma^{x;yz}=-\sigma^{y;zx},\;\imaginary\sigma^{z;xy}_{I}$

    * $\circ$

$G=177-182$ [$\Gamma_{h}$]

$U=\begin{pmatrix}1&0&0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}\\\
0&0&1\\\ 1&2&0\\\ -2&-1&0\end{pmatrix},\;u=\begin{pmatrix}0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1/2}\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\ 1/2\\\
1\\\ 0\end{pmatrix}$

  * •

$F=C_{6v}$ ($c_{6,z}$, $\sigma_{x}$)

$zz;zz;zz=12,\;zz;xx;xx=zz;yy;yy=zz;xy;xy=6,$
$\sigma^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy},\;\sigma^{x;xz}=\sigma^{y;yz}$

    * $\circ$

$G=183-186$ [$\Gamma_{h}$] $\rightarrow$ Same as $G=177$

  * •

$F=D_{3h}$ ($s_{3,z}$, $c_{2,x}$)

$xx;xx;xx=-xx;xy;xy=3,\;$
$\sigma^{x;xx}=-\sigma^{x;yy}=-\sigma^{y;yx}_{I}=-\sigma^{y;yx}_{II}$

    * $\circ$

$G=187-190$ [$\Gamma_{h}$]

$U=\begin{pmatrix}1&0&0\\\ 0&0&-1\\\ 0&0&1\\\ 1&2&0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-2}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix},\;u=\begin{pmatrix}0\\\
0\\\ 1/2\\\ 1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix}$

  * •

$F=T$

$xx;yy;zz=xz;yx;zy=4$ $\sigma^{x;yz}=\sigma^{y;zx}=\sigma^{z;xy}$

    * $\circ$

$G=195$ & $198$ [$\Gamma_{c}$]

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}1/a&0&0\\\ 0&1/a&0\\\
0&0&1/a\end{pmatrix},\;U=\begin{pmatrix}0&1&0\\\ -1&0&1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}\\\
0&-1&1\\\ 1&0&0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\end{pmatrix},\;u=\begin{pmatrix}1/2\\\
0\\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
0\\\ 1/2\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix}$

    * $\circ$

$G=196$ [$\Gamma_{c}^{f}$]

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}-1/a&1/a&1/a\\\
1/a&-1/a&1/a\\\ 1/a&1/a&-1/a\end{pmatrix},\;U=\begin{pmatrix}1&1&1\\\
1&-1&1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\\\
1&0&-1\\\ -1&1&1\\\ 0&1&-1\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}3}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix},\;u=\begin{pmatrix}3/2\\\
1\\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
0\\\ 1\\\ 0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}3/2}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix}$

    * $\circ$

$G=197$ & $199$ [$\Gamma_{c}^{v}$]

$\begin{pmatrix}\bm{G}_{1}^{\text{ref}}\\\ \bm{G}_{2}^{\text{ref}}\\\
\bm{G}_{3}^{\text{ref}}\end{pmatrix}=\begin{pmatrix}0&1/a&1/a\\\ 1/a&0&1/a\\\
1/a&1/a&0\end{pmatrix},\;U=\begin{pmatrix}1&1&2\\\ 1&0&-1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
0&1&-1\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-2}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-2}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\end{pmatrix},\;u=\begin{pmatrix}1\\\
0\\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
0\\\ {\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix}$

  * •

$F=O$

$xx;yy;zz=xz;yx;zy=-xy;yx;zz=-xx;yz;zy=-xz;yy;zx=4$
$\imaginary\sigma^{x;yz}_{I}=\imaginary\sigma^{y;zx}_{I}=\imaginary\sigma^{z;xy}_{I}$

    * $\circ$

$G=207-208$ & $212-213$ [$\Gamma_{c}$]

$U=\begin{pmatrix}1&-1&0\\\ 0&1&0\\\ -1&0&1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\end{pmatrix},\;u=\begin{pmatrix}0\\\
1/2\\\ 0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix}$

    * $\circ$

$G=209-210$ [$\Gamma_{c}^{f}$]

$U=\begin{pmatrix}-1&1&0\\\ 1&1&1\\\ 1&-1&1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1}\\\
1&0&-1\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}3}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix},\;u=\begin{pmatrix}0\\\
3/2\\\ 1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\ 0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}3/2}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix}$

    * $\circ$

$G=211$ & $214$ [$\Gamma_{c}^{v}$]

$U=\begin{pmatrix}-1&1&0\\\ 1&1&2\\\ 1&0&-1\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-1}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-2}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}\end{pmatrix},\;u=\begin{pmatrix}0\\\
1\\\ 0\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{pmatrix}$

  * •

$F=T_{d}$

$xx;yy;zz=xz;yx;zy=xy;yx;zz=xx;yz;zy=xz;yy;zx=4$
$\sigma^{x;yz}_{I}=\sigma^{x;yz}_{II}=\sigma^{y;zx}_{I}=\sigma^{y;zx}_{II}=\sigma^{z;xy}_{I}=\sigma^{z;xy}_{II}$

    * $\circ$

$G=215$ & $218$ [$\Gamma_{c}$] $\rightarrow$ Same as $G=207$

    * $\circ$

$G=216$ & $219$ [$\Gamma_{c}^{f}$] $\rightarrow$ Same as $G=209$

    * $\circ$

$G=217$ & $220$ [$\Gamma_{c}^{v}$] $\rightarrow$ Same as $G=211$

The folded equation for magnetic point groups of type III, (20), can be
readily computed from the previous list. Since the evaluation is slightly more
involved, we next provide the explicit expressions for all non-vanishing
components $\sigma^{a;bc}_{\text{shift}}$ in the non-trivial cases, namely for
those groups $M$ where $H$ is non-centrosymmetric ($i\in
H\Rightarrow\sigma^{a;bc}_{\text{shift}}=0$) or there is
$\pazocal{P}\pazocal{T}$ symmetry ($i\in F-H\Rightarrow\mathbb{D}^{F}=0$), in
which case (20) reduces to

$\sigma^{a;bc}_{\text{shift}}=-\frac{i}{\absolutevalue{H}N^{\text{IBZ}}_{\bm{k}}}\sum_{\bm{k}\in\text{IBZ}}\sum_{a^{\prime},b^{\prime},c^{\prime}}\mathbb{D}^{H}_{a,a^{\prime};b,b^{\prime};c,c^{\prime}}\real\left(\tilde{\sigma}^{a^{\prime};b^{\prime}c^{\prime}}_{\text{shift}}(\bm{k})-\tilde{\sigma}^{a^{\prime};c^{\prime}b^{\prime}}_{\text{shift}}(\bm{k})\right)$

which can be evaluated straightforwardly as (18) or (19). The general form of
the shift conductivity tensor can then be immediately identified from these
expressions for any magnetic point group. The type III groups are labelled
with an arbitrary number and in the Shubnikov-Belov notation as in Reference
55. The orientation of each ordinary point group $H$ and $F$ is as in the
previous list, except when specified by new generators. For simplicity, we
omit the “shift” label, the $\bm{k}$ dependence and the ubiquitous
$\frac{1}{N^{\text{IBZ}}_{\bm{k}}}\sum_{\bm{k}\in\text{IBZ}}$ in the notation
of the conductivity. The general relations
$\sigma^{a;bc}_{\text{shift}}=(\sigma^{a;cb}_{\text{shift}})^{*}$ are also
omitted.

  *   * •

$M=2$ $(2^{\prime})$, $H=C_{1}$, $F=C_{2}$

$\sigma^{a;bc}=\left\\{\begin{aligned}
&\imaginary\left[\tilde{\sigma}^{a;bc}+\tilde{\sigma}^{a;cb}\right],&\text{ if
}{\\#x+\\#y}\text{ is even}\\\
&-i\real\left[\tilde{\sigma}^{a;bc}-\tilde{\sigma}^{a;cb}\right],&\text{ if
}{\\#x+\\#y}\text{ is odd}\end{aligned}\right.$

where $\\#x$ ($\\#y$) is the number of times that $x$ ($y$, resp.) appears in
the ($a,b,c$) triplet.

  * •

$M=3$ $(m^{\prime})$, $H=C_{1}$, $F=C_{s}$

$\sigma^{a;bc}=\left\\{\begin{aligned}
&\imaginary\left[\tilde{\sigma}^{a;bc}+\tilde{\sigma}^{a;cb}\right],&\text{ if
}{\\#x+\\#y}\text{ is odd}\\\
&-i\real\left[\tilde{\sigma}^{a;bc}-\tilde{\sigma}^{a;cb}\right],&\text{ if
}{\\#x+\\#y}\text{ is even}\end{aligned}\right.$

where $\\#x$ ($\\#y$) is the number of times that $x$ ($y$, resp.) appears in
the ($a,b,c$) triplet.

  * •

$M=7$ $(2^{\prime}2^{\prime}2)$, $H=C_{2}$, $F=D_{2}$

$\displaystyle\sigma^{x;xz}=-i\real\left[\tilde{\sigma}^{x;xz}-\tilde{\sigma}^{x;zx}\right],\;\sigma^{y;yz}=-i\real\left[\tilde{\sigma}^{y;yz}-\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=\imaginary\left[\tilde{\sigma}^{x;yz}+\tilde{\sigma}^{x;zy}\right],\;\sigma^{y;zx}=\imaginary\left[\tilde{\sigma}^{y;zx}+\tilde{\sigma}^{y;xz}\right],\;\sigma^{z;xy}=\imaginary\left[\tilde{\sigma}^{z;xy}+\tilde{\sigma}^{z;yx}\right]$

  * •

$M=8$ $(m^{\prime}m^{\prime}2)$, $H=C_{2}$, $F=C_{2v}$

$\displaystyle\sigma^{z;zz}=2\imaginary\tilde{\sigma}^{z;zz},\;\sigma^{z;xx}=2\imaginary\tilde{\sigma}^{z;xx},\;\sigma^{z;yy}=2\imaginary\tilde{\sigma}^{z;yy},\;$
$\displaystyle\sigma^{x;xz}=\imaginary\left[\tilde{\sigma}^{x;xz}+\tilde{\sigma}^{x;zx}\right],\;\sigma^{y;yz}=\imaginary\left[\tilde{\sigma}^{y;yz}+\tilde{\sigma}^{y;zy}\right],\;\sigma^{x;yz}=-i\real\left[\tilde{\sigma}^{x;yz}-\tilde{\sigma}^{x;zy}\right],\;$
$\displaystyle\sigma^{y;zx}=-i\real\left[\tilde{\sigma}^{y;zx}-\tilde{\sigma}^{y;xz}\right],\;\sigma^{z;xy}=-i\real\left[\tilde{\sigma}^{z;xy}-\tilde{\sigma}^{z;yx}\right]$

  * •

$M=9$ $(m^{\prime}m2^{\prime})$, $H=C_{s}$ ($\sigma_{y}$), $F=C_{2v}$

$\displaystyle\sigma^{z;zz}=2\imaginary\tilde{\sigma}^{z;zz},\;\sigma^{z;xx}=2\imaginary\tilde{\sigma}^{z;xx},\;\sigma^{z;yy}=2\imaginary\tilde{\sigma}^{z;yy},\;\sigma^{y;yx}=-i\real\left[\tilde{\sigma}^{y;yx}-\tilde{\sigma}^{y;xy}\right],\;$
$\displaystyle\sigma^{x;xz}=\imaginary\left[\tilde{\sigma}^{x;xz}+\tilde{\sigma}^{x;zx}\right],\;\sigma^{y;yz}=\imaginary\left[\tilde{\sigma}^{y;yz}+\tilde{\sigma}^{y;zy}\right],\;\sigma^{z;zx}=-i\real\left[\tilde{\sigma}^{z;zx}-\tilde{\sigma}^{z;xz}\right]$

  * •

$M=13$ $(4^{\prime})$, $H=C_{2}$, $F=C_{4}$

$\displaystyle\sigma^{z;zz}=2\imaginary\tilde{\sigma}^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy}=\imaginary\left[\tilde{\sigma}^{z;xx}+\tilde{\sigma}^{z;yy}\right],\;$
$\displaystyle\sigma^{x;xz}=(\sigma^{y;yz})^{*}=-\frac{i}{2}\left[\tilde{\sigma}^{x;xz}-(\tilde{\sigma}^{x;zx})^{*}-(\tilde{\sigma}^{y;yz})^{*}+\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=-\sigma^{y;zx}=-\frac{i}{2}\left[\tilde{\sigma}^{x;yz}-(\tilde{\sigma}^{x;zy})^{*}+(\tilde{\sigma}^{y;xz})^{*}-\tilde{\sigma}^{y;zx}\right]$

  * •

$M=14$ $(\overline{4}^{\prime})$, $H=C_{2}$, $F=S_{4}$

$\displaystyle\sigma^{z;xx}=-\sigma^{z;yy}=\imaginary\left[\tilde{\sigma}^{z;xx}-\tilde{\sigma}^{z;yy}\right],\;\sigma^{x;xz}=-(\sigma^{y;yz})^{*}=-\frac{i}{2}\left[\tilde{\sigma}^{x;xz}-(\tilde{\sigma}^{x;zx})^{*}+(\tilde{\sigma}^{y;yz})^{*}-\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=\sigma^{y;zx}=-\frac{i}{2}\left[\tilde{\sigma}^{x;yz}-(\tilde{\sigma}^{x;zy})^{*}-(\tilde{\sigma}^{y;xz})^{*}+\tilde{\sigma}^{y;zx}\right],\;\sigma^{z;xy}=-i\left[\tilde{\sigma}^{z;xy}-(\tilde{\sigma}^{z;yx})^{*}\right]$

  * •

$M=15$ $(42^{\prime}2^{\prime})$, $H=C_{4}$, $F=D_{4}$

$\displaystyle\sigma^{x;xz}=\sigma^{y;yz}=-\frac{i}{2}\real\left[\tilde{\sigma}^{x;xz}-\tilde{\sigma}^{x;zx}+\tilde{\sigma}^{y;yz}-\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=-\sigma^{y;zx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;yz}+\tilde{\sigma}^{x;zy}-\tilde{\sigma}^{y;xz}-\tilde{\sigma}^{y;zx}\right]$

  * •

$M=16$ $(4^{\prime}22^{\prime})$, $H=D_{2}$, $F=D_{4}$

$\sigma^{x;yz}=-\sigma^{y;zx}=-\frac{i}{2}\left[\tilde{\sigma}^{x;yz}-(\tilde{\sigma}^{x;zy})^{*}+(\tilde{\sigma}^{y;xz})^{*}-\tilde{\sigma}^{y;zx}\right]$

  * •

$M=20$ $(4m^{\prime}m^{\prime})$, $H=C_{4}$, $F=C_{4v}$

$\displaystyle\sigma^{z;zz}=2\imaginary\tilde{\sigma}^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy}=\imaginary\left[\tilde{\sigma}^{z;xx}+\tilde{\sigma}^{z;yy}\right],\;$
$\displaystyle\sigma^{x;xz}=\sigma^{y;yz}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;xz}+\tilde{\sigma}^{x;zx}+\tilde{\sigma}^{y;yz}+\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=\sigma^{y;zx}=-\frac{i}{2}\real\left[\tilde{\sigma}^{x;yz}-\tilde{\sigma}^{x;zy}-\tilde{\sigma}^{y;xz}+\tilde{\sigma}^{y;zx}\right],\;\sigma^{z;xy}=-i\real\left[\tilde{\sigma}^{z;xy}-\tilde{\sigma}^{z;yx}\right]$

  * •

$M=21$ $(4^{\prime}mm^{\prime})$, $H=C_{2v}$, $F=C_{4v}$

$\displaystyle\sigma^{z;zz}=2\imaginary\tilde{\sigma}^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy}=\imaginary\left[\tilde{\sigma}^{z;xx}+\tilde{\sigma}^{z;yy}\right],\;$
$\displaystyle\sigma^{x;xz}=(\sigma^{y;yz})^{*}=-\frac{i}{2}\left[\tilde{\sigma}^{x;xz}-(\tilde{\sigma}^{x;zx})^{*}-(\tilde{\sigma}^{y;yz})^{*}+\tilde{\sigma}^{y;zy}\right]$

  * •

$M=22$ $(\overline{4}2^{\prime}m^{\prime})$, $H=S_{4}$, $F=D_{2d}$

$\displaystyle\sigma^{x;xz}=-\sigma^{y;yz}=-\frac{i}{2}\real\left[\tilde{\sigma}^{x;xz}-\tilde{\sigma}^{x;zx}-\tilde{\sigma}^{y;yz}+\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=\sigma^{y;zx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;yz}+\tilde{\sigma}^{x;zy}+\tilde{\sigma}^{y;xz}+\tilde{\sigma}^{y;zx}\right],\;\sigma^{z;xy}=\imaginary\left[\tilde{\sigma}^{z;xy}+\tilde{\sigma}^{z;yx}\right]$

  * •

$M=23$ $(\overline{4}^{\prime}2m^{\prime})$, $H=D_{2}$, $F=D_{2d}$

$\sigma^{x;yz}=\sigma^{y;zx}=-\frac{i}{2}\left[\tilde{\sigma}^{x;yz}-(\tilde{\sigma}^{x;zy})^{*}-(\tilde{\sigma}^{y;xz})^{*}+\tilde{\sigma}^{y;zx}\right],\;\sigma^{z;xy}=-i\left[\tilde{\sigma}^{z;xy}-(\tilde{\sigma}^{z;yx})^{*}\right]$

  * •

$M=24$ $(\overline{4}^{\prime}m2^{\prime})$, $H=C_{2v}$, $F=D_{2d}$
($s_{4,z}$, $\sigma_{x}$)

$\sigma^{z;xx}=-\sigma^{z;yy}=\imaginary\left[\tilde{\sigma}^{z;xx}-\tilde{\sigma}^{z;yy}\right],\;\sigma^{x;xz}=-(\sigma^{y;yz})^{*}=-\frac{i}{2}\left[\tilde{\sigma}^{x;xz}-(\tilde{\sigma}^{x;zx})^{*}+(\tilde{\sigma}^{y;yz})^{*}-\tilde{\sigma}^{y;zy}\right]$

  * •

$M=30$ $(32^{\prime})$, $H=C_{3}$, $F=D_{3}$

$\displaystyle\sigma^{x;xx}=-\sigma^{x;yy}=-\sigma^{y;yx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;xx}-\tilde{\sigma}^{x;yy}-\tilde{\sigma}^{y;yx}-\tilde{\sigma}^{y;xy}\right],\;$
$\displaystyle\sigma^{x;xz}=\sigma^{y;yz}=-\frac{i}{2}\real\left[\tilde{\sigma}^{x;xz}-\tilde{\sigma}^{x;zx}+\tilde{\sigma}^{y;yz}-\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=-\sigma^{y;zx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;yz}+\tilde{\sigma}^{x;zy}-\tilde{\sigma}^{y;xz}-\tilde{\sigma}^{y;zx}\right]$

  * •

$M=31$ $(3m^{\prime})$, $H=C_{3}$, $F=C_{3v}$

$\displaystyle\sigma^{x;xx}=-\sigma^{x;yy}=-\sigma^{y;yx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;xx}-\tilde{\sigma}^{x;yy}-\tilde{\sigma}^{y;yx}-\tilde{\sigma}^{y;xy}\right],\;\sigma^{z;zz}=2\imaginary\tilde{\sigma}^{z;zz},\;$
$\displaystyle\sigma^{z;xx}=\sigma^{z;yy}=\imaginary\left[\tilde{\sigma}^{z;xx}+\tilde{\sigma}^{z;yy}\right],\;\sigma^{x;xz}=\sigma^{y;yz}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;xz}+\tilde{\sigma}^{x;zx}+\tilde{\sigma}^{y;yz}+\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=\sigma^{y;zx}=-\frac{i}{2}\real\left[\tilde{\sigma}^{x;yz}-\tilde{\sigma}^{x;zy}-\tilde{\sigma}^{y;xz}+\tilde{\sigma}^{y;zx}\right],\;\sigma^{z;xy}=-i\real\left[\tilde{\sigma}^{z;xy}-\tilde{\sigma}^{z;yx}\right]$

  * •

$M=32$ $(\overline{6}^{\prime})$, $H=C_{3}$, $F=C_{3h}$

$\displaystyle\sigma^{x;xx}=-\sigma^{x;yy}=-\sigma^{y;yx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;xx}-\tilde{\sigma}^{x;yy}-\tilde{\sigma}^{y;yx}-\tilde{\sigma}^{y;xy}\right],\;$
$\displaystyle\sigma^{y;yy}=-\sigma^{y;xx}=-\sigma^{x;xy}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{y;yy}-\tilde{\sigma}^{y;xx}-\tilde{\sigma}^{x;xy}-\tilde{\sigma}^{x;yx}\right],\;$
$\displaystyle\sigma^{x;xz}=\sigma^{y;yz}=-\frac{i}{2}\real\left[\tilde{\sigma}^{x;xz}-\tilde{\sigma}^{x;zx}+\tilde{\sigma}^{y;yz}-\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=\sigma^{y;zx}=-\frac{i}{2}\real\left[\tilde{\sigma}^{x;yz}-\tilde{\sigma}^{x;zy}-\tilde{\sigma}^{y;xz}+\tilde{\sigma}^{y;zx}\right],\;\sigma^{z;xy}=-i\real\left[\tilde{\sigma}^{z;xy}-\tilde{\sigma}^{z;yx}\right]$

  * •

$M=33$ $(\overline{6}m^{\prime}2^{\prime})$, $H=C_{3h}$, $F=D_{3h}$

$\sigma^{x;xx}=-\sigma^{x;yy}=-\sigma^{y;yx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;xx}-\tilde{\sigma}^{x;yy}-\tilde{\sigma}^{y;yx}-\tilde{\sigma}^{y;xy}\right]$

  * •

$M=34$ $(\overline{6}^{\prime}m2^{\prime})$, $H=C_{3v}$, $F=D_{3h}$

$\displaystyle\sigma^{x;xx}=-\sigma^{x;yy}=-\sigma^{y;yx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;xx}-\tilde{\sigma}^{x;yy}-\tilde{\sigma}^{y;yx}-\tilde{\sigma}^{y;xy}\right],\;$
$\displaystyle\sigma^{x;xz}=\sigma^{y;yz}=-\frac{i}{2}\real\left[\tilde{\sigma}^{x;xz}-\tilde{\sigma}^{x;zx}+\tilde{\sigma}^{y;yz}-\tilde{\sigma}^{y;zy}\right]$

  * •

$M=35$ $(\overline{6}^{\prime}m^{\prime}2)$, $H=D_{3}$, $F=D_{3h}$

$\displaystyle\sigma^{x;xx}=-\sigma^{x;yy}=-\sigma^{y;yx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;xx}-\tilde{\sigma}^{x;yy}-\tilde{\sigma}^{y;yx}-\tilde{\sigma}^{y;xy}\right],\;$
$\displaystyle\sigma^{x;yz}=\sigma^{y;zx}=-\frac{i}{2}\real\left[\tilde{\sigma}^{x;yz}-\tilde{\sigma}^{x;zy}-\tilde{\sigma}^{y;xz}+\tilde{\sigma}^{y;zx}\right],\;\sigma^{z;xy}=-i\real\left[\tilde{\sigma}^{z;xy}-\tilde{\sigma}^{z;yx}\right]$

  * •

$M=36$ $(6^{\prime})$, $H=C_{3}$, $F=C_{6}$

$\displaystyle\sigma^{z;zz}=2\imaginary\tilde{\sigma}^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy}=\imaginary\left[\tilde{\sigma}^{z;xx}+\tilde{\sigma}^{z;yy}\right],\;$
$\displaystyle\sigma^{x;xz}=\sigma^{y;yz}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;xz}+\tilde{\sigma}^{x;zx}+\tilde{\sigma}^{y;yz}+\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=-\sigma^{y;zx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;yz}+\tilde{\sigma}^{x;zy}-\tilde{\sigma}^{y;xz}-\tilde{\sigma}^{y;zx}\right]$

  * •

$M=41$ $(62^{\prime}2^{\prime})$, $H=C_{6}$, $F=D_{6}$

$\displaystyle\sigma^{x;xz}=\sigma^{y;yz}=-\frac{i}{2}\real\left[\tilde{\sigma}^{x;xz}-\tilde{\sigma}^{x;zx}+\tilde{\sigma}^{y;yz}-\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=-\sigma^{y;zx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;yz}+\tilde{\sigma}^{x;zy}-\tilde{\sigma}^{y;xz}-\tilde{\sigma}^{y;zx}\right]$

  * •

$M=42$ $(6^{\prime}2^{\prime}2)$, $H=D_{3}$, $F=D_{6}$

$\sigma^{x;yz}=-\sigma^{y;zx}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;yz}+\tilde{\sigma}^{x;zy}-\tilde{\sigma}^{y;xz}-\tilde{\sigma}^{y;zx}\right]$

  * •

$M=46$ $(6m^{\prime}m^{\prime})$, $H=C_{6}$, $F=C_{6v}$

$\displaystyle\sigma^{z;zz}=2\imaginary\tilde{\sigma}^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy}=\imaginary\left[\tilde{\sigma}^{z;xx}+\tilde{\sigma}^{z;yy}\right],\;$
$\displaystyle\sigma^{x;xz}=\sigma^{y;yz}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;xz}+\tilde{\sigma}^{x;zx}+\tilde{\sigma}^{y;yz}+\tilde{\sigma}^{y;zy}\right],\;$
$\displaystyle\sigma^{x;yz}=\sigma^{y;zx}=-\frac{i}{2}\real\left[\tilde{\sigma}^{x;yz}-\tilde{\sigma}^{x;zy}-\tilde{\sigma}^{y;xz}+\tilde{\sigma}^{y;zx}\right],\;\sigma^{z;xy}=-i\real\left[\tilde{\sigma}^{z;xy}-\tilde{\sigma}^{z;yx}\right]$

  * •

$M=47$ $(6^{\prime}m^{\prime}m)$, $H=C_{3v}$, $F=C_{6v}$

$\displaystyle\sigma^{z;zz}=2\imaginary\tilde{\sigma}^{z;zz},\;\sigma^{z;xx}=\sigma^{z;yy}=\imaginary\left[\tilde{\sigma}^{z;xx}+\tilde{\sigma}^{z;yy}\right],\;$
$\displaystyle\sigma^{x;xz}=\sigma^{y;yz}=\frac{1}{2}\imaginary\left[\tilde{\sigma}^{x;xz}+\tilde{\sigma}^{x;zx}+\tilde{\sigma}^{y;yz}+\tilde{\sigma}^{y;zy}\right]$

  * •

$M=54$ $(\overline{4}^{\prime}3m^{\prime})$, $H=T$, $F=T_{d}$

$\sigma^{x;yz}=\sigma^{y;zx}=\sigma^{z;xy}=-\frac{i}{3}\left[\tilde{\sigma}^{x;yz}-(\tilde{\sigma}^{x;zy})^{*}+\tilde{\sigma}^{y;zx}-(\tilde{\sigma}^{y;xz})^{*}+\tilde{\sigma}^{z;xy}-(\tilde{\sigma}^{z;yx})^{*}\right]$

  * •

$M=55$ $(4^{\prime}32^{\prime})$, $H=T$, $F=O$

$\sigma^{a;bc}=0,\;\forall a,b,c$

##  4 Conclusive remarks

The use of Gaussian basis sets to compute the shift conductivity has been
proved satisfactory. The analytical evaluation of real-space integrals
involving localized functions allows to readily compute the Berry connection
and velocity matrix elements, while the reduced dimension yields lighter
calculations, and the economical option of hybrid functionals allows to easily
reproduce the desired band gap with good precision in a wide variety of
systems. Furthermore, the (magnetic) space group symmetry is fully preserved
and it can be capitalised on to perform maximally-efficient reciprocal space
summations, in addition to immediately discerning the contributions to the
electrical current under any light polarization.

The numerical results for the chosen materials are in standard agreement with
the literature, and in all cases the length and velocity gauges have been
shown to yield nearly identical outcomes. One would be tempted to conclude
that the velocity gauge, in view of its comparative simplicity, should then be
the preferred option in general. However, the length gauge makes no use of
completeness relations and larger discrepancies may appear when employing
smaller bases. In addition, care should be taken when separating the shift and
injection contributions without time-reversal symmetry in the velocity gauge.
Nevertheless, the use of as-largest-as-possible basis sets is generally
advisable, typically between TZVP and QZVP and including diffuse exponents.
The use of ghost atoms, i.e., basis functions not located on atomic positions,
could be explored in order to facilitate the reproduction of particularly-
delocalized empty conduction states.

We note that all other single-particle contributions to the BPVE (in
particular, the injection current) and to the total second-order optical
response (in particular, the second-harmonic generation) can be computed from
this method since the corresponding expressions involve the same basic
quantities as the shift current. The transformation properties of these other
optical contributions are the same for spatial operations, but the role of
time-reversal may change. For example, the injection conductivity has the
opposite behaviour to $\sigma_{\text{shift}}$, in the sense that time-reversal
symmetry forces it to be imaginary (instead of real); and the results in
Section 3 can be adapted from that. The evaluation of metallic systems is also
possible in the length gauge, albeit it introduces additional terms with Fermi
surface derivatives 33.

## Supporting Information

Supporting Information available:

Input files for the self-consistent electronic structure calculations in
CRYSTAL23 for each material (.d12 files, in the terminology of the code). Band
structures along high-symmetry lines for each material. Explicit comparison
between the BZ and IBZ summations in BaTiO3.

The authors acknowledge financial support from Spanish MICINN (Grant Nos.
PID2019-109539GB-C43 & TED2021-131323B-I00 & PID2022-141712NB-C21), María de
Maeztu Program for Units of Excellence in R&D (Grant No. CEX2018-000805-M),
Comunidad Autónoma de Madrid through the Nanomag COST-CM Program (Grant No.
S2018/NMT-4321), Generalitat Valenciana through Programa Prometeo (2021/017),
Centro de Computación Científica of the Universidad Autónoma de Madrid, and
Red Española de Supercomputación.

The authors declare no competing financial interest.

## References

  * Sturman and Fridkin 1992 Sturman, B. I.; Fridkin, V. M. _The photovoltaic and photorefractive effects in noncentrosymmetric materials_ ; Gordon and Breach Science: Philadelphia, 1992
  * Chynoweth 1956 Chynoweth, A. G. Surface space-charge layers in barium titanate. _Phys. Rev._ 1956, _102_ , 705–714
  * Chen 1969 Chen, F. S. Optically induced change of refractive indices in LiNbO3 and LiTaO3. _J. Appl. Phys._ 1969, _40_ , 3389–3396
  * Glass et al. 1974 Glass, A. M.; Von der Linde, D.; Negran, T. J. High-voltage bulk photovoltaic effect and the photorefractive process in LiNb03. _Appl. Phys. Lett._ 1974, _25_ , 233–235
  * Koch et al. 1975 Koch, W. T. H.; Munser, R.; Ruppel, W.; Würfel, P. Bulk photovoltaic effect in BaTiO3. _Solid State Commun._ 1975, _17_ , 847–850
  * Fridkin et al. 1977 Fridkin, V. M.; Popov, B. N.; Verkhovskaya, K. A. Effect of anomalous bulk photovoltage in ferroelectrics. _Phys. Status Solidi A_ 1977, _39_ , 193–201
  * Kraut and von Baltz 1979 Kraut, W.; von Baltz, R. Anomalous bulk photovoltaic effect in ferroelectrics: a quadratic response theory. _Phys. Rev. B_ 1979, _19_ , 1548–1554
  * von Baltz and Kraut 1981 von Baltz, R.; Kraut, W. Theory of the bulk photovoltaic effect in pure crystals. _Phys. Rev. B_ 1981, _23_ , 5590–5596
  * Hornung et al. 1983 Hornung, D.; Von Baltz, R.; Rössler, U. Band structure investigation of the bulk photovoltaic effect in n-gap. _Solid State Commun._ 1983, _48_ , 225–229
  * Fridkin et al. 1993 Fridkin, V.; Dalba, G.; Fornasini, P.; Soldo, Y.; Rocca, F.; Burattini, E. The bulk photovoltaic effect in LiNbO3, crystals under x-ray synchrotron radiation. _Ferroelectrics Lett._ 1993, _16_ , 1–5
  * Batirov et al. 1997 Batirov, T.; Doubovik, E.; Djalalov, R.; Fridkin, V. M. The bulk photovoltaic effect in the piezoelectric crystal Pr3Ga5SiO14. _Ferroelectrics Lett._ 1997, _23_ , 95–98
  * Buse 1997 Buse, K. Light-induced charge transport processes in photorefractive crystals II: materials. _Appl. Phys. B_ 1997, _64_ , 391–407
  * Král et al. 2000 Král, P.; Mele, E. J.; Tománek, D. Photogalvanic effects in heteropolar nanotubes. _Phys. Rev. Lett._ 2000, _85_ , 1512–1515
  * Sipe and Shkrebtii 2000 Sipe, J. E.; Shkrebtii, A. I. Second-order optical response in semiconductors. _Phys. Rev. B_ 2000, _61_ , 5337–5352
  * Butler et al. 2015 Butler, K. T.; Frost, J. M.; Walsh, A. Ferroelectric materials for solar energy conversion: photoferroics revisited. _Energ. Environ. Sci._ 2015, _8_ , 838–848
  * Spanier et al. 2016 Spanier, J. E.; Fridkin, V. M.; Rappe, A. M.; Akbashev, A. R.; Polemi, A.; Qi, Y.; Gu, Z.; Young, S. M.; Hawley, C. J.; Imbrenda, D., et al. Power conversion efficiency exceeding the Shockley-Queisser limit in a ferroelectric insulator. _Nat. Photonics_ 2016, _10_ , 611–616
  * Dai and Rappe 2023 Dai, Z.; Rappe, A. M. Recent progress in the theory of bulk photovoltaic effect. _Chem. Phys. Rev._ 2023, _4_ , 011303.1–26
  * Cook et al. 2017 Cook, A. M.; Fregoso, B. M.; de Juan, F.; Coh, S.; Moore, J. E. Design principles for shift current photovoltaics. _Nat. Commun._ 2017, _8_ , 14176.1–9
  * Rangel et al. 2017 Rangel, T.; Fregoso, B. M.; Mendoza, B. S.; Morimoto, T.; Moore, J. E.; Neaton, J. B. Large bulk photovoltaic effect and spontaneous polarization of single-layer monochalcogenides. _Phys. Rev. Lett._ 2017, _119_ , 067402.1–6
  * Osterhoudt et al. 2019 Osterhoudt, G. B.; Diebel, L. K.; Gray, M. J.; Yang, X.; Stanco, J.; Huang, X.; Shen, B.; Ni, N.; Moll, P. J. W.; Ran, Y., et al. Colossal mid-infrared bulk photovoltaic effect in a type-I Weyl semimetal. _Nat. Mater._ 2019, _18_ , 471–475
  * Ahn et al. 2020 Ahn, J.; Guo, G. Y.; Nagaosa, N. Low-frequency divergence and quantum geometry of the bulk photovoltaic effect in topological semimetals. _Phys. Rev. X_ 2020, _10_ , 041041.1–28
  * Wang et al. 2020 Wang, M.; Wei, H.; Wu, Y.; Jia, J.; Yang, C.; Chen, Y.; Chen, X.; Cao, B. Polarization-enhanced bulk photovoltaic effect of BiFeO3 epitaxial film under standard solar illumination. _Phys. Lett. A_ 2020, _384_ , 126831.1–6
  * Xu et al. 2021 Xu, H.; Wang, H.; Zhou, J.; Li, J. Pure spin photocurrent in noncentrosymmetric crystals: bulk spin photovoltaic effect. _Nat. Commun._ 2021, _12_ , 4330.1–9
  * Blázquez-Martínez et al. 2022 Blázquez-Martínez, A.; Grysan, P.; Girod, S.; Glinsek, S.; Granzow, T. Direct evidence for bulk photovoltaic charge transport in a ferroelectric polycrystalline film. _Scripta Mater._ 2022, _211_ , 114498.1–5
  * Chaudhary et al. 2022 Chaudhary, S.; Lewandowski, C.; Refael, G. Shift-current response as a probe of quantum geometry and electron-electron interactions in twisted bilayer graphene. _Phys. Rev. Res._ 2022, _4_ , 013164.1–14
  * Zhang et al. 2022 Zhang, C.; Guo, P.; Zhou, J. Tailoring bulk photovoltaic effects in magnetic sliding ferroelectric materials. _Nano Lett._ 2022, _22_ , 9297–9305
  * Dai et al. 2021 Dai, Z.; Schankler, A. M.; Gao, L.; Tan, L. Z.; Rappe, A. M. Phonon-assisted ballistic current from first-principles calculations. _Phys. Rev. Lett._ 2021, _126_ , 177403.1–6
  * Xu et al. 2022 Xu, H.; Wang, H.; Li, J. Nonlinear nonreciprocal photocurrents under phonon dressing. _Phys. Rev. B_ 2022, _106_ , 035102.1–7
  * Dai and Rappe 2021 Dai, Z.; Rappe, A. M. First-principles calculation of ballistic current from electron-hole interaction. _Phys. Rev. B_ 2021, _104_ , 235203.1–6
  * Wang et al. 2022 Wang, H.; Tang, X.; Xu, H.; Li, J.; Qian, X. Generalized Wilson loop method for nonlinear light-matter interaction. _npj Quantum Mater._ 2022, _7_ , 61.1–7
  * Parker et al. 2019 Parker, D. E.; Morimoto, T.; Orenstein, J.; Moore, J. E. Diagrammatic approach to nonlinear optical response with application to Weyl semimetals. _Phys. Rev. B_ 2019, _99_ , 045121.1–20
  * Ventura et al. 2017 Ventura, G. B.; Passos, D. J.; Lopes Dos Santos, J. M. B.; Viana Parente Lopes, J. M.; Peres, N. M. R. Gauge covariances and nonlinear optical responses. _Phys. Rev. B_ 2017, _96_ , 035431.1–11
  * Taghizadeh et al. 2017 Taghizadeh, A.; Hipolito, F.; Pedersen, T. G. Linear and nonlinear optical response of crystals using length and velocity gauges: effect of basis truncation. _Phys. Rev. B_ 2017, _96_ , 195413.1–10
  * Passos et al. 2018 Passos, D. J.; Ventura, G. B.; Lopes, J. M., Viana Parente Lopes; Lopes dos Santos, J. M. B.; Peres, N. M. R. Nonlinear optical responses of crystalline systems: results from a velocity gauge analysis. _Phys. Rev. B_ 2018, _97_ , 235446.1–11
  * Young and Rappe 2012 Young, S. M.; Rappe, A. M. First principles calculation of the shift current photovoltaic effect in ferroelectrics. _Phys. Rev. Lett._ 2012, _109_ , 116601.1–5
  * Wang et al. 2017 Wang, C.; Liu, X.; Kang, L.; Gu, B. L.; Xu, Y.; Duan, W. First-principles calculation of nonlinear optical responses by Wannier interpolation. _Phys. Rev. B_ 2017, _96_ , 115147.1–9
  * Ibañez-Azpiroz et al. 2018 Ibañez-Azpiroz, J.; Tsirkin, S. S.; Souza, I. Ab initio calculation of the shift photocurrent by Wannier interpolation. _Phys. Rev. B_ 2018, _97_ , 245143.1–13
  * Aversa and Sipe 1995 Aversa, C.; Sipe, J. E. Nonlinear optical susceptibilities of semiconductors: results with a length-gauge analysis. _Phys. Rev. B_ 1995, _52_ , 14636–14645
  * Esteve-Paredes and Palacios 2023 Esteve-Paredes, J. J.; Palacios, J. J. A comprehensive study of the velocity, momentum and position matrix elements for Bloch states: application to a local orbital basis. _SciPost Phys. Core_ 2023, _6_ , 002.1–20
  * Doll et al. 2006 Doll, K.; Dovesi, R.; Orlando, R. Analytical Hartree-Fock gradients with respect to the cell parameter: systems periodic in one and two dimensions. _Theor. Chem. Acc._ 2006, _115_ , 354–360
  * 41 Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; Zicovich-Wilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; D’Arco, P.; Llunell, M.; Causà, M.; Noël, Y.; Maschio, L.; Erba, A.; Rerat, M.; Casassa, S. CRYSTAL17 user’s manual. https://www.crystal.unito.it/include/manuals/crystal17.pdf, Accesed on 07/2023
  * Helgaker et al. 2013 Helgaker, T.; Jorgensen, P.; Olsen, J. _Molecular electronic-structure theory_ ; John Wiley & Sons: New York, 2013
  * Gradshteyn and M. 2007 Gradshteyn, I. S.; M., R. I. _Table of integrals, series, and products_ , seventh ed.; Elsevier: New York, 2007; Page 365, section 3.462, formula 2
  * Bir and Pikus 1974 Bir, G. L.; Pikus, G. E. _Symmetry and strain-induced effects in semiconductors_ ; John Wiley & Sons: New York, 1974
  * Dovesi et al. 2018 Dovesi, R.; Erba, A.; Orlando, R.; Zicovich-Wilson, C. M.; Civalleri, B.; Maschio, L.; Rerat, M.; Casassa, S.; Baima, J.; Salustro, S.; Kirtman, B. Quantum-mechanical condensed matter simulations with CRYSTAL. _WIREs Comput. Mol. Sci._ 2018, _8_ , e1360.1–36
  * Erba et al. 2022 Erba, A.; Desmarais, J. K.; Casassa, S.; Civalleri, B.; Donà, L.; Bush, I. J.; Searle, B.; Maschio, L.; Daga, L.-E.; Cossard, A.; Ribaldone, C.; Ascrizzi, E.; Marana, N. L.; Flament, J.-P.; Kirtman, B.; Dovesi, R.; Erba, A.; Orlando, R.; Zicovich-Wilson, C. M.; Civalleri, B.; Maschio, L.; Rerat, M.; Casassa, S.; Baima, J.; Salustro, S.; Kirtman, B. CRYSTAL23: a program for computational solid state physics and chemistry. _J. Chem. Theory Comput._ 2022,
  * Pritchard et al. 2019 Pritchard, B. P.; Altarawy, D.; Didier, B.; Gibson, T. D.; Windus, T. L. New basis set exchange: an open, up-to-date resource for the molecular sciences community. _J. Chem. Inf. Model._ 2019, _59_ , 4814–4820
  * Laun and Bredow 2022 Laun, J.; Bredow, T. BSSE-corrected consistent Gaussian basis sets of triple-zeta valence with polarization quality of the fifth period for solid-state calculations. _J. Comput. Chem._ 2022, _43_ , 839–846
  * Heyd et al. 2005 Heyd, J.; Peralta, J. E.; Scuseria, G. E.; Martin, R. L. Energy band gaps and lattice parameters evaluated with the Heyd-Scuseria-Ernzerhof screened hybrid functional. _J. Chem. Phys._ 2005, _123_ , 174101.1–8
  * Laun and Bredow 2021 Laun, J.; Bredow, T. BSSE-corrected consistent Gaussian basis sets of triple-zeta valence with polarization quality of the sixth period for solid-state calculations. _J. Comput. Chem._ 2021, _42_ , 1064–1072
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. _Phys. Rev. Lett._ 1996, _77_ , 3865.1–4