# Universal responses in nonmagnetic polar metals

Fabian Jäger  Materials Theory, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093
Zurich, Switzerland  Nicola A. Spaldin  Materials Theory, ETH Zurich,
Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland  Sayantika Bhowal
Materials Theory, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093 Zurich,
Switzerland

###### Abstract

We demonstrate that two phenomena, the kinetic magneto-electric effect and the
non-linear Hall effect, are universal to polar metals, as a consequence of
their coexisting and contraindicated polarization and metallicity. We show
that measurement of the effects provides a complete characterization of the
nature of the polar metal, in that the non-zero response components indicate
the direction of the polar axis, and the coefficients change sign on
polarization reversal and become zero in the non-polar phase. We illustrate
our findings for the case of electron-doped PbTiO3 using a combination of
density functional theory and model Hamiltonian-based calculations. Our model
Hamiltonian analysis provides crucial insight into the microscopic origin of
the effects, showing that they originate from inversion-symmetry-breaking-
induced inter-orbital hoppings, which cause an asymmetric charge density
quantified by odd-parity charge multipoles. Our work both heightens the
relevance of the kinetic magneto-electric and non-linear Hall effects, and
broadens the platform for investigating and detecting odd-parity charge
multipoles in metals.

##  I Introduction

The idea of combining electric polarization with metallicity, against the
common belief that polarization is screened by itinerant carriers, was first
conceived by Anderson and Blount more than fifty years ago Anderson and Blount
(1965). It has come to reality, however, rather recently with the practical
material realization of polar metals Shi et al. (2013); Fei et al. (2018);
Sharma et al. (2019). These have consequently opened up a new paradigm for
investigating numerous intriguing physical effects that result from the
coexistence of the seemingly mutually exclusive properties of polarity and
metallicity Zhou and Ariando (2020); Bhowal and Spaldin (2023); Hickox-Young
et al. (2023).

In the present work, we point out two such effects, the kinetic magneto-
electric effect (KME) and non-linear Hall effect (NHE), which are universal to
all polar metals. While these effects have been sporadically investigated in
some candidate polar metal systems Ma et al. (2019); Kang et al. (2019); Xiao
et al. (2020), a consensus in applying these effects to characterizing polar
metals is still missing. Here we show that both these effects carry
simultaneously the key signatures of the polar metal phase, that is the
direction of the polar axis, the switchability of the polarization, and the
ferroelectric-like nonpolar to polar structural transition, and so provide a
complete characterization of polar metals. Furthermore, we reveal the
microscopic origin of these two effects by analyzing asymmetries in the charge
density. While both effects are dominated by contributions from the electric
dipole moment i.e., the first-order asymmetry in the charge density, the
electric octupole moment, characterizing the third-order asymmetry in the
charge density, also plays an important role.

The kinetic magneto-electric effect is a linear effect, describing electric
field ($\cal E$) induced magnetization ${\cal M}_{j}={\cal K}_{ij}{\cal
E}_{i}$ in a nonmagnetic metal Levitov et al. (1984); Yoda et al. (2015);
Zhong et al. (2016); Şahin et al. (2018); Tsirkin et al. (2018). The resulting
magnetization, in turn, gives rise to a transverse Hall current ($J$) as a
second-order response to the applied electric field, $J_{i}=\chi_{ijk}{\cal
E}_{j}{\cal E}_{k}$, known as nonlinear Hall effect Sodemann and Fu (2015).
${\cal K}_{ij}$ and $\chi_{ijk}$ are the KME response and non-linear Hall
conductivity (NHC) tensor respectively with $i,j,k$ indicating the Cartesian
directions. Within the relaxation-time approximation for the nonequilibrium
electron distribution, both these responses can be elegantly recast in terms
of the equilibrium reciprocal-space magnetic (spin plus orbital) moment
$\vec{m}(\vec{k})$ and the Berry curvature dipole (BCD) ${\cal D}_{ij}$
respectively Zhong et al. (2016); Tsirkin et al. (2018); Sodemann and Fu
(2015):

$\displaystyle{\cal
K}_{ij}=-\frac{e\tau}{\hbar}\sum_{n}\frac{1}{(2\pi)^{3}}\int
d^{3}k~{}m^{n}_{j}(\vec{k})\partial_{k_{i}}\epsilon^{n}_{k}\Big{(}\frac{\partial
f_{0}}{\partial\epsilon^{n}_{k}}\Big{)}=\frac{e\tau}{\hbar}\tilde{\cal
K}_{ij}$
$\displaystyle\chi_{ijk}=-\varepsilon_{ilk}\frac{e^{3}\tau}{2(1+i\omega\tau)}{\cal
D}_{jl}$
$\displaystyle=\varepsilon_{ilk}\frac{e^{3}\tau}{2(1+i\omega\tau)}\sum_{n}\frac{1}{(2\pi)^{3}}\int
d^{3}k~{}\Omega^{n}_{l}(\vec{k})\partial_{k_{j}}\epsilon_{k}\Big{(}\frac{\partial
f_{0}}{\partial\epsilon^{n}_{k}}\Big{)}.$ (1)

Here $\tau$, $\varepsilon_{adc}$, $e,n$, $f_{0}$ and $\vec{\Omega}$ are
respectively the relaxation time-constant, Levi-Civita symbol, the electronic
charge, band index, equilibrium Fermi distribution function and Berry
curvature. Both the reduced KME response $\tilde{\cal K}_{ij}$ and the BCD are
intrinsic properties of a material and are given by,Zhong et al. (2016);
Tsirkin et al. (2018); Ma et al. (2019); Kang et al. (2019)

$\displaystyle\tilde{\cal K}_{ij}$ $\displaystyle=$
$\displaystyle-\sum_{n}\frac{1}{(2\pi)^{3}}\int
d^{3}k~{}m^{n}_{j}(\vec{k})\partial_{k_{i}}\epsilon^{n}_{k}\Big{(}\frac{\partial
f_{0}}{\partial\epsilon^{n}_{k}}\Big{)}$ (2) $\displaystyle=$
$\displaystyle\sum_{n}\frac{1}{(2\pi)^{3}}\int
d^{3}k~{}(\partial_{k_{i}}m^{n}_{j})f_{0}$

and,

$\displaystyle{\cal D}_{ij}$ $\displaystyle=$
$\displaystyle-\sum_{n}\frac{1}{(2\pi)^{3}}\int
d^{3}k~{}\Omega^{n}_{j}(\vec{k})\partial_{k_{i}}\epsilon_{k}\Big{(}\frac{\partial
f_{0}}{\partial\epsilon^{n}_{k}}\Big{)}$ (3) $\displaystyle=$
$\displaystyle\sum_{n}\frac{1}{(2\pi)^{3}}\int
d^{3}k~{}(\partial_{k_{i}}\Omega^{n}_{j})f_{0}.$

Both $\tilde{\cal K}_{ij}$ and ${\cal D}_{ij}$ are allowed in nonmagnetic
metals with gyrotropic point group symmetry Zhong et al. (2016); Tsirkin et
al. (2018); Ma et al. (2019); Kang et al. (2019). Since all polar point groups
are gyrotropic gyr ; He and Law (2020), both KME and NHE are allowed by
symmetry in all polar metals.

Interestingly, the components of the reciprocal-space magnetic moment
$\vec{m}(\vec{k})$, that contributes to the KME response, is determined by the
direction of the electric polarization Bhowal et al. (2022a). Similarly, the
antisymmetric component of the BCD, ${\cal D}^{-}=({\cal D}-{\cal D}^{T})/2$,
correlates with the orientation of the polar axis $\vec{d}$,
$d_{i}\equiv\varepsilon_{ijk}{\cal D}_{jk}^{-}/2$ Sodemann and Fu (2015),
suggesting a possible switching of both responses for a switchable orientation
of the polar distortion. Furthermore, since both effects are forbidden by
symmetry in an inversion symmetric structure, a structural transition from a
centrosymmetric to a noncentrosymmetric polar structure can be inferred from
the onset of these effects as the temperature is lowered.

We illustrate these concepts by explicitly considering the case of electron-
doped PbTiO3 (PTO) as an example material. Undoped PTO is a prototypical
conventional ferroelectric insulator Nelmes and Kuhs (1985). Interestingly,
even upon electron doping via replacing the Ti4+ ions by Nb5+ ions, the
resulting PbTi1-xNbxO3 was observed to sustain the electric polarization up to
$x=0.12$, at which point the system also becomes conducting Gu et al. (2017);
Iijima et al. (2000). In the present work, using both first-principles density
functional theory (DFT) and a model Hamiltonian-based approach we show that
the presence as well as the orientation of the polar axis in the polar metal
phase of doped PTO can be determined from the non-zero components of KME and
the NHE.

The remainder of this manuscript is organized as follows. We start by
describing the computational details in section II. This is followed by the
results and discussions in section III, where we present our computational
results for doped PTO, describing the existence and tuning of KME and NHE, the
momentum space distribution of the orbital moment and Berry curvature that
determine these effects, their microscopic origin within the model Hamiltonian
framework, and the role of odd-parity charge multipoles. Finally, we summarize
our results in section IV and give a proposal for measuring these effects.

##  II Computational details

The responses $\tilde{\cal K}_{ij}$ and ${\cal D}_{ij}$ are computed using the
QUANTUM ESPRESSOGiannozzi et al. (2009) and Wannier90 codes Marzari and
Vanderbilt (1997); Souza et al. (2001); Mostofi et al. (2014). We use fully
relativistic norm-conserving pseudo-potentials for all the atoms with the
following valence electron configurations: Pb ($6s^{2}6p^{2}$), Ti
($4s^{2}3d^{2}$), and O ($2s^{2}2p^{4}$). Self-consistency is achieved with a
12$\times$ 12$\times$10 $k$-point mesh and a convergence threshold of
$10^{-7}$ Ry. The ab-initio wave functions, thus obtained, are then projected
to maximally localized Wannier functions Marzari and Vanderbilt (1997); Souza
et al. (2001) using the Wannier90 code Mostofi et al. (2014). In the
disentanglement process, as initial projections, we choose 42 Wannier
functions per unit cell which include the $s$ and $p$ orbitals of Pb, $d$
orbitals of Ti and $s$ and $p$ orbitals of O atoms, excluding the rest. After
the disentanglement is achieved, the wannierisation process is converged to
$10^{-10}$ Å2. We then compute the $k$-space distribution of the orbital
moment and the Berry curvature as well as the reduced KME response,
$\tilde{\cal K}_{ij}$ and the BCD ${\cal D}_{ij}$ for a 150$\times$
150$\times$140 $k$-point mesh. To estimate the doped charge density, we also
compute the densities of states (DOS) for the same $k$-point mesh.

##  III Results and Discussion

###  III.1 NHE and KME and their tuning in polar metals

We start with the electronic structure of PTO, which crystallizes in the non-
centrosymmetric tetragonal ($P4mm$) structure with the polar $C_{4v}$ point
group symmetry Nelmes and Kuhs (1985). In tetragonal PTO, both Pb2+ ($6s^{2}$
lone pair) and Ti4+ ($3d^{0}$) ions off-center with respect to the surrounding
O2- ions, resulting in a net polarization along $\hat{z}$, which is switchable
to $-\hat{z}$ using an external electric field. We refer to these two
structures, schematically depicted in Fig. 1 (a), as $+P$ and $-P$
respectively. The electronic structure of the polar undoped PTO (corresponding
to $+P$) is shown in Fig. 1 (b), depicting the insulating band structure in
which the occupied O-$p$ states and the fromally empty Ti-$t_{2g}$ states form
the valence band maximum (VBM) and conduction band minimum (CBM) respectively.

Figure 1: (a) Schematic illustration of the crystal structure of PTO, showing
the off-centering of the Ti atom leading to a polarization $+P$ along
$\hat{z}$. The dashed circle indicates the displacement of the Ti atom in the
opposite direction, switching the direction of the polarization ($-P$)
indicated by the dashed arrow. (b) Comparison of the band structure of undoped
PTO, computed within DFT (dashed line) and that obtained from Wannier90 (solid
line), showing a good agreement between the two. (c) Computed reduced KME
response components $\tilde{\cal K}_{xy}$ (blue) and $\tilde{\cal K}_{yx}$
(red) for the two directions of polarization $+P$ (circle) and $-P$ (dashed
line), shown in (a), as a function of energy. (d) BCD components ${\cal
D}_{xy}$ (blue) and ${\cal D}_{yx}$ (red) for the polarization directions,
$+P$ (circles) and $-P$ (dashed line) as a function of energy. (e) Energy
variation of the spin and orbital contributions to the absolute value of the
antisymmetric component of the reduced KME response, $\tilde{\cal
K}^{-}_{xy}=\frac{1}{2}(\tilde{\cal K}_{xy}-\tilde{\cal K}_{yx}$) for
polarization $+P$. (f) Energy variation of the doped electron densities for
the two polarization directions, $+P$ (circles) and $-P$ (dashed line). The
vertical black dashed line corresponds to the experimentally achieved maximum
doped electron density, that maintains the polarity of the structure. The zero
of energy in (b)-(f) refers to the CBM of undoped PTO.

The doping electrons in doped PTO occupy the CBM, leading to a metallic band
structure within the rigid band approximation. In order to compute
$\tilde{\cal K}_{ij}$ and ${\cal D}_{ij}$, we first project the computed ab-
initio wave functions onto maximally localized Wannier functions, and then
disentangle the relevant bands (see section II for computational details) from
the rest using the Wannier90 code Mostofi et al. (2014). As depicted in Fig. 1
(b), the wannierised bands agree well with the full DFT band structure. The
central quantities $\tilde{\cal K}_{ij}$ and ${\cal D}_{ij}$ in determining
the magnitudes of the KME and NHE are then computed using Eqs. (2) and (3) as
implemented within the Wannier90 code Marzari and Vanderbilt (1997); Souza et
al. (2001); Mostofi et al. (2014).

Figure 2: (a) Comparison of the energy variation of the absolute value of
$\tilde{\cal K}^{-}_{xy}$ for two different displacements of the Ti ion, $P$
and $P_{1}$, with the former being larger than the latter. (b) Comparison of
the energy variation of the absolute value of the antisymmetric BCD component
${\cal D}^{-}_{xy}=\frac{1}{2}({\cal D}_{xy}-{\cal D}_{yx}$) for the same $P$
and $P_{1}$.

The computed non-zero components of the reduced KME response, $\tilde{\cal
K}_{xy}$ (blue circle), $\tilde{\cal K}_{yx}$ (red circle), and BCD ${\cal
D}_{xy}$ (blue circle) and ${\cal D}_{yx}$ (red circle) are shown as a
function of energy in Fig. 1 (c) and (d) for the $+P$ structure. To determine
whether the energy range used in the computation is experimentally achievable,
we further compute the doped electron density by integrating the corresponding
DOS and show the results in Fig. 1 (f). Note that the zero of the energy
corresponds to the CBM for the undoped case. The vertical dashed line in Fig.
1 (f) indicates the maximum doped electron density up to which the polarity of
the lattice persists in the experiments Gu et al. (2017); Iijima et al.
(2000), justifying the chosen energy range.

We note from Fig. 1 (c) and (d) that $\tilde{\cal K}_{xy}=-\tilde{\cal
K}_{yx}$ and ${\cal D}_{xy}=-{\cal D}_{yx}$, consistent with the $C_{4v}$
point group symmetry. Here $\tilde{\cal K}_{ij}$ has both spin and orbital
contributions. In order to understand the relative contributions of the two,
the individual spin and orbital contributions are also shown in Fig. 1 (e) for
the absolute value of the antisymmetric component of the reduced KME response,
$\tilde{\cal K}^{-}_{xy}=\frac{1}{2}(\tilde{\cal K}_{xy}-\tilde{\cal
K}_{yx})$. This clearly shows that the orbital contribution dominates over the
spin contribution. Such a current-induced orbital magnetization has also been
reported for other systems with broken inversion symmetry Yoda et al. (2015);
He et al. (2020); Bhowal and Satpathy (2020); Hara et al. (2020) and may have
important implications in the field of orbitronics.

To see the effect of the polarization direction, we reverse the direction of
the displacement of the ions, leading to the $-P$ structure (see Fig. 1 (a)).
The corresponding computed $\tilde{\cal K}_{xy}$ (blue dashed line),
$\tilde{\cal K}_{yx}$ (red dashed line), and ${\cal D}_{xy}$ (blue dashed
line), ${\cal D}_{yx}$ (red dashed line) are shown in Figs. 1 (c) and (d)
respectively. We note that in this case, all the computed quantities switch
sign compared to the $+P$ structure, while still maintaining the symmetry of
the $C_{4v}$ point group, as discussed above.

We further artificially decrease the amount of the Ti displacement to see the
effect of the magnitude of polarization. We refer to the corresponding
structure as $+P_{1}$. The computed absolute values of $\tilde{\cal
K}^{-}_{xy}$ and ${\cal D}^{-}_{xy}=\frac{1}{2}({\cal D}_{xy}-{\cal D}_{yx})$
for $+P_{1}$ are depicted in Figs. 2 (a) and (b) respectively, together with
the values for $+P$. We find that both $\tilde{\cal K}^{-}_{xy}$ and ${\cal
D}^{-}_{xy}$ have smaller magnitude for $+P_{1}$ compared to $+P$, suggesting
that both effects not only depend on the direction of polarization but also
depend on the magnitude of the polarization.

It is important to point out here that polarization is not the only factor
that contributes to the value of the responses. For example, both responses
also depend on the details of the electronic structure (see Eqs. 2 and 3). As
a result, the situation can be more complicated if there is a drastic change
in the band structure with the change in electric polarization. Nevertheless,
our analysis clearly shows that the polarization is an important factor and
that both KME and NHE are tunable by changing the direction or magnitude of
the electric polarization.

Figure 3: Comparison of the energy variation of the absolute value of (a)
$\tilde{\cal K}^{-}_{xy}$ and (b) ${\cal D}^{-}_{xy}$ in the absence and
presence of SOC.

Furthermore, to understand the dependence on the spin-orbit coupling (SOC), we
also perform additional calculations with the SOC turned off in our
computations. Comparisons of the computed $\tilde{\cal K}^{-}_{xy}$ and ${\cal
D}^{-}_{xy}$ both in the absence and presence of SOC are shown in Figs. 3 (a)
and (b). As seen from these figures, both $\tilde{\cal K}^{-}_{xy}$ and ${\cal
D}^{-}_{xy}$ exist even without the SOC. This suggests that both effects occur
due to the symmetry of the structure and the presence of SOC is not necessary.
Indeed, in the absence of SOC, the KME response is driven by the orbital
contribution. With the inclusion of SOC, the orbital degrees of freedom couple
to the spin degrees of freedom, and consequently, it leads to additional
current-induced spin magnetization in the system. The inclusion of SOC,
therefore, increases the magnitudes of both effects. It is important to point
out here that unlike these responses, the spin-splitting of the bands and the
resulting unconventional magnetic Compton scattering Bhowal et al. (2022b)
occur only in the presence of SOC.

###  III.2 $k$-space distribution of orbital moment and Berry curvature

Figure 4: Computed $k$-space distributions of the orbital moment components
(in units of eV-Å2) (a) $m^{\rm orb}_{x}$, (b) $m^{\rm orb}_{y}$, and the
components of the Berry curvature (in units of Å2) (c) $\Omega_{x}$, and (d)
$\Omega_{y}$ in PTO for the doped electron density $\sim 0.001$Å-3. The black
circle around the $\Gamma$ point indicates the Fermi surface corresponding to
this electron density. $\vec{m}^{\rm orb}(\vec{k})$ ($\vec{\Omega}(\vec{k})$)
is the sum of $\vec{m}^{n,\rm orb}(\vec{k})$ ($\Omega^{n}_{k}(\vec{k})$),
given in Eq. (4) (Eq. (5)), over the occupied bands. The high-symmetry
$k$-points are indicated on the $k_{z}=0$ plane for easy reference.

To better understand the responses, we further compute the $k$-space
distributions of the relevant orbital magnetic moment components $m^{\rm
orb}_{x}(\vec{k}),m^{\rm orb}_{y}(\vec{k})$ and Berry curvature components
$\Omega_{x}(\vec{k}),\Omega_{y}(\vec{k})$ in the $k_{x}$-$k_{y}$ plane. Since,
the $\tilde{\cal K}_{ij}$ response is dominated by the orbital contribution,
here for simplicity we only consider the orbital magnetic moment distribution.
The orbital magnetic moment is computed within the modern theory by evaluating
the expectation value of the orbital magnetization operator
$\frac{-e}{2}(\vec{r}\times\vec{v})$ Xiao et al. (2005); Thonhauser et al.
(2005); Ceresoli et al. (2006); Shi et al. (2007) with $-e<0$, as implemented
in the Wannier90 code Lopez et al. (2012),

$\displaystyle\vec{m}^{n,\rm orb}(\vec{k})$ $\displaystyle=$
$\displaystyle\frac{e}{2\hbar}{\rm Im}\langle\nabla_{k}u^{n}_{k}|\times[{\cal
H}(\vec{k})-\epsilon^{n}_{k}]|\nabla_{k}u^{n}_{k}\rangle$ (4) $\displaystyle+$
$\displaystyle\frac{e}{\hbar}{\rm
Im}\langle\nabla_{k}u^{n}_{k}|\times[\epsilon^{n}_{k}-\epsilon_{\rm
F}]|\nabla_{k}u^{n}_{k}\rangle.$

Here, $\epsilon^{n}_{k}$ and $u^{n}_{k}$ are the energy eigenvalues and
eigenfunctions of the Hamiltonian ${\cal H}(\vec{k})$ obtained from
Wannierization, and $\epsilon_{\rm F}$ is the Fermi energy. We note that since
the KME response is a Fermi surface property (see Eq. LABEL:formula), the
second term in Eq. 4 does not contribute to the KME response. This can be seen
easily by recognizing that $\frac{\partial
f_{0}}{\partial\epsilon^{n}_{k}}=-\delta(\epsilon^{n}_{k}-\epsilon_{\rm F})$,
and so has a non-zero value only if $\epsilon^{n}_{k}=\epsilon_{\rm F}$, in
which case the second term in Eq. (4) vanishes. The $k$-space distribution of
the Berry curvature is computed using the Kubo formula Thouless et al. (1982),

$\displaystyle\Omega^{n}_{k}(\vec{k})$ $\displaystyle=$
$\displaystyle-2\hbar^{2}\sum_{m\neq n}{\rm Im}\frac{\langle
u^{n}_{k}|v_{i}|u^{m}_{k}\rangle\langle
u^{m}_{k}|v_{j}|u^{n}_{k}\rangle}{(\epsilon_{k}^{n}-\epsilon_{k}^{m})^{2}},$
(5)

where $\vec{v}=\frac{1}{\hbar}\frac{\partial{\cal H}}{\partial\vec{k}}$ is the
velocity operator and ($i,j,k$) are cyclic permutations of the Cartesian
directions ($x,y,z$).

Both $\vec{m}^{\rm orb}(\vec{k})$ and $\vec{\Omega}(\vec{k})$ follow the same
symmetry relations: Under spatial inversion $\cal I$ symmetry both remain
invariant, with $\vec{m}^{\rm orb}(\vec{k})\xrightarrow{\cal I}\vec{m}^{\rm
orb}(-\vec{k})$, whereas under time-reversal ($\cal T$) symmetry they switch
sign, $\vec{m}^{\rm orb}(\vec{k})\xrightarrow{\cal T}-\vec{m}^{\rm
orb}(-\vec{k})$ (similarly for $\vec{\Omega}(\vec{k})$). Hence, for a non-zero
$\vec{m}^{\rm orb}(\vec{k})$ ($\vec{\Omega}(\vec{k})$), either of these two
symmetries must be broken. In the present case, the broken $\cal I$ symmetry
leads to non-zero values of $\vec{m}^{\rm orb}(\vec{k})$ and
$\vec{\Omega}(\vec{k})$. We plot our calculated $\vec{m}^{\rm orb}(\vec{k})$
and $\vec{\Omega}(\vec{k})$ in Fig. 4. Note that since $\cal T$ symmetry is
preserved, $\vec{m}^{\rm orb}$ ($\vec{\Omega}$) at $+\vec{k}$ has the opposite
sign to that at $-\vec{k}$, and as a result the sum of $\vec{m}^{\rm
orb}(\vec{k})$ ($\vec{\Omega}(\vec{k})$) over the occupied part of the
Brillouin zone (BZ) is zero, consistent with the overall nonmagnetic behavior
of PTO.

Figure 5: (a) Vertical ($\sigma_{v}$) and (b) diagonal ($\sigma_{d}$) mirror
planes in PTO.

The key features of the computed distributions in Figs. 4 (a)-(d) are the
following. First of all, $m^{\rm orb}_{x}$ ($\Omega_{x}$) is equal and
opposite at $\pm k_{y}$, while it has the same sign at $\pm k_{x}$, consistent
with the $\sigma_{v}$ mirror symmetries (see Fig. 5 (a)) that dictate:

$\displaystyle m^{\rm orb}_{x}(k_{x},k_{y},k_{z})\xrightarrow{M_{100}}m^{\rm
orb}_{x}(-k_{x},k_{y},k_{z})$ and $\displaystyle m^{\rm
orb}_{x}(k_{x},k_{y},k_{z})\xrightarrow{M_{010}}-m^{\rm
orb}_{x}(k_{x},-k_{y},k_{z}).$ (6)

In contrast, $m^{\rm orb}_{y}$ ($\Omega_{y}$) is equal and opposite at $\pm
k_{x}$, while having the same sign at $\pm k_{y}$ due to the same $\sigma_{v}$
symmetries, that is,

$\displaystyle m^{\rm orb}_{y}(k_{x},k_{y},k_{z})\xrightarrow{M_{100}}-m^{\rm
orb}_{y}(-k_{x},k_{y},k_{z})$ and $\displaystyle m^{\rm
orb}_{y}(k_{x},k_{y},k_{z})\xrightarrow{M_{010}}m^{\rm
orb}_{y}(k_{x},-k_{y},k_{z}).$ (7)

Furthermore, the $x$ and $y$ components of $\vec{m}^{\rm orb}(\vec{k})$
($\vec{\Omega}(\vec{k})$) are related to each other by the mirror
$M_{1\bar{1}0}$ symmetry (see Fig. 5 (b)), viz., $m^{\rm
orb}_{x}(k_{x},k_{y},k_{z})\xrightarrow{M_{1\bar{1}0}}-m^{\rm
orb}_{y}(k_{y},k_{x},k_{z})$. Moreover, since the velocity operator transforms
as $(v_{x},v_{y},v_{z})\xrightarrow{M_{1\bar{1}0}}(v_{y},v_{x},v_{z})$ under
the mirror $M_{1\bar{1}0}$ symmetry, Eq. (2) [Eq. (3)] leads to the constraint
$\tilde{\cal K}_{xy}=-\tilde{\cal K}_{yx}$ [${\cal D}_{xy}=-{\cal D}_{yx}$],
in agreement with our results in Fig. 1 (c) [Fig. 1 (d)].

###  III.3 Microscopic origin: role of odd-parity charge multipoles

Model Hamiltonian- To understand the microscopic origin of these effects, we
construct a minimal tight-binding (TB) model in the basis set of the
Ti-$t_{2g}$ orbitals, $\\{d_{xy},d_{yz},d_{xz}\\}$. For small doping, the
doped electrons occupy the bands around the $\Gamma$ point of the BZ that
correspond to the CBM for the undoped case, indicated by the black circles in
Fig. 4. We, therefore, expand the TB model around the $\Gamma$ point, and the
resulting low energy model Hamiltonian is given by

$\displaystyle{\cal H}(\vec{k})={\cal H}_{\rm inv}(\vec{k})+{\cal H}_{\rm
BI}(\vec{k}).$ (8)

Here ${\cal H}_{\rm inv}$ is the inversion symmetric part of the Hamiltonian,
and is given by,

$\cal{H}_{\rm inv}=\begin{pmatrix}h_{11}&h_{12}&h_{13}\\\
h_{12}&h_{22}&h_{23}\\\ h_{13}&h_{23}&h_{33}\end{pmatrix},$ (9)

with the explicit analytical forms of the elements $h_{ij}$ up to quadratic
order in $k$ given below:

$\displaystyle h_{11}$ $\displaystyle=$ $\displaystyle
t_{\text{eff}}^{1}-t_{\text{eff}}^{2}(k_{x}^{2}+k_{y}^{2})a^{2}-t_{\text{eff}}^{3}k_{z}^{2}c^{2}$
$\displaystyle h_{22}$ $\displaystyle=$ $\displaystyle
t_{\text{eff}}^{4}-t_{\text{eff}}^{5}k_{x}^{2}a^{2}-t_{\text{eff}}^{6}k_{y}^{2}a^{2}-t_{\text{eff}}^{7}k_{z}^{2}c^{2}$
$\displaystyle h_{33}$ $\displaystyle=$ $\displaystyle
t_{\text{eff}}^{4}-t_{\text{eff}}^{6}k_{x}^{2}a^{2}-t_{\text{eff}}^{5}k_{y}^{2}a^{2}-t_{\text{eff}}^{7}k_{z}^{2}c^{2}$
$\displaystyle h_{12}$ $\displaystyle=$ $\displaystyle t_{\rm
eff}^{8}k_{x}k_{z}ac$ $\displaystyle h_{13}$ $\displaystyle=$ $\displaystyle
t_{\rm eff}^{8}k_{y}k_{z}ac$ $\displaystyle h_{23}$ $\displaystyle=$
$\displaystyle t_{\rm eff}^{9}k_{x}k_{y}a^{2}.$ (10)

Here $a$ and $c$ are the lattice constants for the tetragonal unit cell. Note
that since $\cal{H}_{\rm inv}$ is inversion symmetric, it contains only terms
that are even in $k$. The effective hopping parameters $t_{\text{eff}}^{i}$,
$i=1,9$ are linear combinations of the different effective $t_{2g}$-$t_{2g}$
electronic hopping parameters that we extract using the Nth order muffin-tin
orbital (NMTO) downfolding technique Andersen and Saha-Dasgupta (2000a). The
computed parameters for one direction of polarization ($+P$) are listed in
Table 1. We considered up to fourth nearest neighbor (NN) interactions. It is
important to consider further neighbor interactions which are needed to
capture the physics of the two effects of interest, as we discuss later.

Table 1: Effective hopping parameters (in units of $10^{-2}$ Ry) in Eq. (III.3), derived from the computed TB hopping parameters and onsite energies for PTO using the NMTO downfolding technique. $t^{1}_{\rm eff}$ | $t^{2}_{\rm eff}$ | $t^{3}_{\rm eff}$ | $t^{4}_{\rm eff}$ | $t^{5}_{\rm eff}$ | $t^{6}_{\rm eff}$ | $t^{7}_{\rm eff}$ | $t^{8}_{\rm eff}$ | $t^{9}_{\rm eff}$
---|---|---|---|---|---|---|---|---
4.95 | -2.26 | 0.4 | 10.09 | -0.19 | -0.97 | -1.59 | -0.48 | -0.28

On the other hand, ${\cal H}_{\rm BI}$ includes the hopping parameters that
are induced by the broken ${\cal I}$ symmetry. It can be written in terms of
the components of the orbital angular momentum operator $\hat{\vec{L}}$,

$\displaystyle{\cal H}_{\rm BI}$ $\displaystyle=$ $\displaystyle\frac{\alpha
a}{\hbar}(k_{x}\hat{L}_{y}-k_{y}\hat{L}_{x})-\frac{\alpha
a^{3}}{6\hbar}(k_{x}^{3}\hat{L}_{y}-k_{y}^{3}\hat{L}_{x})$ $\displaystyle-$
$\displaystyle\frac{\beta
ac^{2}}{\hbar}k_{z}^{2}(k_{x}\hat{L}_{y}-k_{y}\hat{L}_{x})-\frac{\gamma
a^{3}}{\hbar}k_{x}k_{y}(k_{y}\hat{L}_{y}-k_{x}\hat{L}_{x}).$

The parameters $\alpha,\beta,\gamma$ are determined by the broken ${\cal
I}$-symmetry-induced hopping parameters and have opposite signs for $+P$ and
$-P$. In centrosymmetric PTO, $\alpha,\beta,\gamma$ are zero so that
$\cal{H}=\cal{H}_{\rm inv}$. In addition, $t_{\rm eff}^{8}=-2(t^{x}-t^{y})=0$
in the centrosymmetric structure, where $t^{x}$ and $t^{y}$ are the fourth
nearest neighbor inter-orbital hopping integrals, which we discuss in detail
later.

The components of the orbital angular momentum operator in Eq. (III.3) in the
$t_{2g}$ orbital basis $\\{d_{xy},d_{yz},d_{xz}\\}$ are given by,

$\displaystyle L_{x}^{(t_{2g})}=$ $\displaystyle\hbar\begin{pmatrix}0&0&-i\\\
0&0&0\\\ i&0&0\end{pmatrix},\quad
L_{y}^{(t_{2g})}=\hbar\begin{pmatrix}0&i&0\\\ -i&0&0\\\ 0&0&0\end{pmatrix},$
(12) $\displaystyle L_{z}^{(t_{2g})}=\hbar\begin{pmatrix}0&0&0\\\ 0&0&i\\\
0&-i&0\end{pmatrix}.$

The advantage of writing ${\cal H}_{\rm BI}$ in terms of the $\hat{\vec{L}}$
operators is that we can readily identify the resulting orbital texture in
momentum space. For example, the first term in Eq. (III.3), which is linear in
$\vec{k}$, depicts a toroidal arrangement of orbital magnetic moment in
reciprocal space (see the inset of Fig. 6 (a)). Such a toroidal arrangement of
the orbital moment in $k$ space is also in agreement with our DFT results (see
Fig. 4 (a) and (b)) and the symmetry analysis presented in section III.2. We
note that the first term in Eq. (III.3) has a form
$\sim(\vec{k}\times\vec{L})$, which is an orbital counterpart of the (spin)
Rashba effect $\sim(\vec{k}\times\vec{\sigma})$ and, hence, is often referred
to as an orbital Rashba effect Go et al. (2017, 2021). In the presence of SOC,
the orbital texture in the orbital Rashba effect couples to the spin,
additionally leading to spin texture and the Rashba effect in PTO Arras et al.
(2019). The Rashba spin-splitting $\Delta\varepsilon_{s}(\vec{k})$ is
antisymmetric in $\vec{k}$, corresponding to $p$-wave symmetry, due to the
presence of time-reversal symmetry, which means that
$\Delta\varepsilon_{s}(\vec{k})=\varepsilon_{\uparrow}(\vec{k})-\varepsilon_{\downarrow}(\vec{k})=-\Delta\varepsilon_{s}(-\vec{k})$.
Here, for simplicity, we do not include SOC in our model Hamiltonian in Eq.
(8), since both KME and BCD exist even in its absence (See Fig. 3).

Role of odd-parity charge multipoles- Interestingly, each term of different
order in $\vec{k}$ in the Hamiltonian ${\cal H}_{\rm BI}$ of Eq. (III.3) has a
direct correlation to a corresponding odd-parity charge multipole. Recently,
we showed that the $k$-space orbital and spin textures in ferroelectrics
result from the $k$-space magnetoelectric multipoles that are reciprocal to
the real-space odd-parity charge multipoles Bhowal et al. (2022b). The odd-
parity charge multipoles characterize the asymmetries in the charge density
that are present due to the broken ${\cal I}$ symmetry. For example, the
electric dipole dictates the first-order asymmetry in the charge density,
while the electric octupole corresponds to the third-order asymmetry, and so
on. The first term within the parentheses in Eq. (III.3), which is linear in
$\vec{k}$, corresponds to the $k$-space representation of the electric dipole
moment ($p_{10}$) whereas the remaining terms, which are all cubic in
$\vec{k}$, correspond to the electric octupole moment (${\cal O}_{30}$).

Figure 6: Atomic-site charge dipole moment component $p_{10}$ and octupole
moment component ${\cal O}_{30}$ on the Ti4+ ions as a function of the
displacement (in units of out-of-plane lattice constant $c$) of the Ti ion
from the center of the unit cell in PTO. The inset shows the schematic for the
toroidal arrangement of the orbital (spin) moment (indicated in thick arrows)
in the $k_{x}$-$k_{y}$ plane due to the first term in the Hamiltonian (III.3)
driven by the charge dipole. (b) Fourth nearest neighbor Ti atoms (connected
by the brown straight lines) along $(\pm a,0,\pm c)$ and $(0,\pm a,\pm c)$.
Note that in the cubic high-symmetry structure with $c=a$, these are second
nearest neighbors.

To verify the existence of the local electric dipoles and octupoles in PTO, we
decompose the $\cal T$ symmetric density matrix
$\rho_{lm,l^{\prime}m^{\prime}}$, computed within the DFT framework, into
parity-odd tensor moments and explicitly compute the atomic-site electric
dipole and octupole moments, for which only the odd $l-l^{\prime}$ terms
contribute Spaldin et al. (2013). The computed odd-parity charge multipoles on
the Ti${{}^{4}+}$ ions are non-zero in the polar structure, as shown in Fig. 6
(a), and confirm the presence of a ferrotype ordering of electric dipole
component $p_{10}$ and octupole component ${\cal O}_{30}$ at the Ti site. Here
the indices at the suffix of the multipole components represent the $l$ and
$m$ indices of the spherical harmonics that are used to build these charge
multipoles. The electric dipole moment $\vec{p}$ is a tensor of rank 1
(vector), with $p_{10}$ indicating its $z$ component. Similarly, the octupole
moment ${\cal O}_{ijk}$ is a totally symmetric tensor of rank 3 with seven
components. The ${\cal O}_{30}$ component has the representation
$\frac{1}{2}z(5z^{2}-r^{2})$.

Results and discussion- Now that we have correlated the individual terms of
the Hamiltonian to the charge multipoles, we diagonalize the Hamiltonian
${\cal H}(\vec{k})$ in Eq. (8) for the realistic parameters listed in Table 1,
extracted using the NMTO downfolding technique Löwdin (2004); Andersen and
Saha-Dasgupta (2000b). We, then, use the computed eigenvalues
$\epsilon_{k}^{n}$ and eigenfunctions $u_{k}^{n}$ to obtain the $k$-space
distribution of the orbital moment and the Berry curvature using Eqs. (4) and
(5) for the lowest energy band of the Hamiltonian in Eq. (8).

Note that the second term in Eq. (4) does not contribute to the KME response,
as stated before, and hence, we ignore this term for the computation of the
orbital moment. We then compute the BCD density
$d_{ij}(\vec{k})=\partial_{k_{i}}\Omega_{j}(\vec{k})$ and the reduced KME
density ${\kappa}_{ij}(\vec{k})=\partial_{k_{i}}m^{\rm orb}_{j}(\vec{k})$ for
$i,j=x,y$, the integrals of which over the occupied part of the BZ determine
the magnitude of ${\cal D}_{ij}$ and $\tilde{\cal K}_{ij}$ respectively [see
Eqs. (2) and (3)]. The computed densities show that they have the same sign (+
or -) over $k$-space only if $i\neq j$ and hence when integrated over the
occupied part of the BZ, only the $xy$ and $yx$ components of ${\cal D}$ and
$\tilde{\cal K}$ have non-zero values. The variations of these components
along a specific momentum direction are shown in Fig. 7 (see the solid lines).

For the opposite polarization direction ($-P$), the parameters
$\alpha,\beta,\gamma$ switch sign and, consequently, as shown in Fig. 7, the
$xy$ and $yx$ components of $d$ and $\kappa$ switch signs, keeping their
magnitudes unaltered. In an ${\cal I}$-symmetric system, on the other hand,
$\alpha=\beta=\gamma=0$, and consequently, we find that $d_{ij}$,
$\kappa_{ij}$ become zero as shown in the insets of Fig. 7, emphasizing the
important role of ${\cal I}$ symmetry breaking.

Figure 7: Results of the tight-binding analysis. Computed variation (solid
line with circles) of the reduced KME density components (a) $\kappa_{xy}$,
(b) $\kappa_{yx}$, and the BCD density components (c) $d_{xy}$, (d) $d_{yx}$
around the $\Gamma$ point for the $+P$ polarization. The variation is shown
along $k_{y}$ for (a) and (c), and along $k_{x}$ for (b) and (d). The same
variation (indicated in dashed lines with diamonds) for $-\kappa_{xy}$,
$-\kappa_{yx}$, $-d_{xy}$, and $-d_{yx}$ for the polarization $-P$ are also
shown in (a)-(d). The same variation of (e) $\kappa_{xy}$, (f) $\kappa_{yx}$,
(g) $d_{xy}$ and (h) $d_{yx}$ in the presence of inversion symmetry (black
solid line), in absence of fourth NN inter-orbital hopping parameters $t^{x}$
and $t^{y}$ (dashed brown line), and in absence of the first term (linear in
$\vec{k}$) in Eq. (III.3) (green line with circles). The parameters used for
the plots are listed in Table 1, and $\alpha=0.22,\beta=0.02,$ and
$\gamma=-0.10$ (in units of 10-2 Ry) for $+P$ polarization.

Further to gain insight into the origin of these two effects, we switch off
the linear term in Eq. (III.3), which originates from the electric dipole
moment. Interestingly, in this case, we find that while all the considered
components of $d$, $\kappa$ still survive, their values reduce drastically by
an order of magnitude. This suggests that the linear terms in $k$ in Eq.
(III.3), originating from the electric dipole moment, play an important role
in determining the magnitudes of both these effects, although the importance
of the electric octupole-driven $k^{3}$ terms can not be ignored. Our findings
are consistent with the multipole description of the KME response, proposed by
Hayami et. al. based on symmetry analysis Hayami et al. (2018). Indeed, we
find that the antisymmetric part of the KME response ${\cal K}^{-}_{ij}$ in
PTO can be described by the existence of an electric dipole moment component,
$\tilde{\cal K}^{-}_{ij}=\frac{1}{2}({\cal K}_{ij}-{\cal
K}_{ji})=\varepsilon_{ijk}p_{k}$. It is important to point out here that the
KME, although universal to all polar metals, can also occur in
noncentrosymmetric but non-polar systems, e.g., chiral materials, in which
case other multipoles such as the monopole of the electric toroidal dipole
moment will dictate the symmetric part (with the trace) of the KME response
Hayami et al. (2018).

We further note that the fourth NN (see Fig. 6 (b)), inter-orbital
($d_{xy}-d_{xz}$ and $d_{xy}-d_{yz}$) hopping integrals, $t^{x}$ and $t^{y}$,
induced by the broken $\cal I$ symmetry, are the key ingredients for both
these effects. While both these hopping integrals contribute to the parameters
$\alpha$ and $\beta$, $\beta$ is solely determined by $t^{x}$ and $t^{y}$
while $\alpha$ has additional contributions. As a result, in the absence of
these hoppings, $\beta$ and the effective hopping, $t_{\rm eff}^{8}$, in
$\cal{H}_{\rm inv}$ vanish. In this case of $t^{x}=t^{y}=0$, we find that the
components of both $d$ and $\kappa$ also vanish, as shown in the insets of
Fig. 7 (see the dashed brown line), emphasizing the importance of the further
neighbor interactions.

To understand why the fourth NN hopping parameters are crucial, we first note
that the non-zero $\beta$ and $t_{\rm eff}^{8}$ resulting from the fourth NN
hopping parameters appear in the third term of Eq. (III.3) and the off-
diagonal elements $h_{12}$ and $h_{13}$ of Eq. (9) respectively.
Interestingly, these are the only inter-orbital contributions in our minimal
model that are also responsible for the band dispersion along the out-of-plane
$k_{z}$ direction. Since the inter-orbital hopping parameters drive the non-
zero Berry curvature Bhowal and Satpathy (2019) and since the dispersion along
$k_{z}$ is crucial for the existence of the in-plane components of both
orbital moment and Berry curvature [see Eqs. (4) and (5)], we see that both
quantities vanish in the absence of fourth NN hopping. This, in turn, also
leads to an absence of $xy$ and $yx$ components of $d$ and $\kappa$,
explaining the crucial role of the fourth NN hopping integrals in driving the
KME and NHE in doped PTO.

##  IV Summary and Outlook

To summarize, taking the example of doped PTO, we have shown that both the KME
and the NHE, are universal to all polar metals and can be used for a complete
characterization of this class of materials. Our work paves the way for the
broad applicability of these two effects in polar metals in general, going
beyond their earlier investigation in topological systems Zhong et al. (2016);
Johansson et al. (2018); Tsirkin et al. (2018); Roy and Narayan (2022). Our
detailed tight-binding analysis reveals the importance of the broken-symmetry-
induced inter-orbital hopping parameters, correlated to the odd-parity charge
multipoles, in mediating these effects. In particular, we have identified the
broken-inversion-induced fourth NN inter-orbital hopping parameters as being
essential in driving these effects in doped PTO.

Proposal for experiments. Before concluding, here we briefly discuss possible
routes to detecting the two effects. The second-order NHE in polar metals can
be detected by measuring the second harmonic current $J_{2\omega}$ at a
frequency $2\omega$ for an applied ac electric field $\vec{E}$ of frequency
$\omega$, Sodemann and Fu (2015)

$\vec{j}_{2\omega}=\frac{e^{3}\tau}{2(1+i\omega\tau)}{\vec{E}_{\omega}}\times({\vec{p}}\times{\vec{E}_{\omega}}).$
(13)

Here $\vec{p}$ is the direction of the electric dipole moment, which is along
$\hat{z}$ for doped PTO. This suggests that for $\vec{E}$ along $\hat{z}$
(i.e., with polar angle $\theta=0$), the Hall current vanishes as we found
also from our explicit calculations discussed above. Furthermore, for a
general form of the field, $\vec{E}=Ee^{i\omega
t}(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$, it is also easy to see
from Eq. (13) that the Hall current does not depend on the azimuthal angle
$\phi$ made by $\vec{E}$ with $\hat{x}$ for an in-plane $\vec{E}$ (i.e.,
$\theta=\pi/2$). This means that rotation of $\vec{E}$ within the $x$-$y$
plane will leave the Hall current invariant.

The current-induced magnetization in the KME should be detectable using the
magneto-optical Kerr effect. In doped PTO the generated magnetization is
dominated by the orbital moment for a reasonable doping concentration (see the
inset of Fig. 1 (c)) and has a magnitude of $~{}1.8\times 10^{-4}\mu_{B}/$atom
at the experimentally observed maximum doping concentration
($n_{x=0.12}=1.9\times 10^{21}$ cm-3) up to which the system retains the
ferroelectricity, for an applied field of $10^{5}$ V/m and a typical
relaxation time constant $\tau\simeq 1$ ps. The computed orbital magnetization
is about four orders of magnitude larger than that reported in Te Tsirkin et
al. (2018), while it is about an order of magnitude smaller than the orbital
magnetization in BCC iron Lopez et al. (2012). The computed total (spin plus
orbital) magnetization of $\sim 1.0\times 10^{-3}\mu_{B}$ per unit cell is
also comparable to the magnetization of the Rashba system Bi/Ag(111), the
(001) surface of the topological insulator $\alpha$-Sn, and the Weyl semimetal
TaAs Johansson et al. (2018) and, hence, likely to be discernible in
measurements.

In the present work, we considered a rigid band approximation to describe the
doped PTO case. While we expect this to provide a good description of the NHE
and KME for the small doping concentration achievable in the measurements,
future work should investigate computationally how electron doping affects the
electronic structure of PTO. The dominance of the orbital magnetization in the
KME response of doped PTO that emerges from our work, opens the door for the
application of polar metals in orbitronics with the additional advantage of
switchable orbital texture by reversal of the electric polarization. We hope
that our work will motivate both theoretical and experimental work in these
directions in the near future.

## Acknowledgements

The authors thank Awadhesh Narayan and Dominic Varghese for stimulating
discussions. NAS and SB were supported by the ERC under the EU’s Horizon 2020
Research and Innovation Programme grant No 810451 and by the ETH Zurich.
Computational resources were provided by ETH Zurich’s Euler cluster, and the
Swiss National Supercomputing Centre, project ID eth3.

## References

  * Anderson and Blount (1965) P. W. Anderson and E. I. Blount, Phys. Rev. Lett. 14, 217 (1965), URL https://link.aps.org/doi/10.1103/PhysRevLett.14.217.
  * Shi et al. (2013) Y. Shi, Y. Guo, X. Wang, A. J. Princep, D. Khalyavin, P. Manuel, Y. Michiue, A. Sato, K. Tsuda, S. Yu, et al., Nat. Mater. 12, 1024 (2013), ISSN 1476-1122, URL www.nature.com/naturematerialshttp://www.nature.com/articles/nmat3754.
  * Fei et al. (2018) Z. Fei, W. Zhao, T. A. Palomaki, B. Sun, M. K. Miller, Z. Zhao, J. Yan, X. Xu, and D. H. Cobden, Nature 560, 336 (2018), URL https://doi.org/10.1038/s41586-018-0336-3.
  * Sharma et al. (2019) P. Sharma, F.-X. Xiang, D.-F. Shao, D. Zhang, E. Y. Tsymbal, A. R. Hamilton, and J. Seidel, Sci. Adv. 5, eaax5080 (2019), eprint https://www.science.org/doi/pdf/10.1126/sciadv.aax5080, URL https://www.science.org/doi/abs/10.1126/sciadv.aax5080.
  * Zhou and Ariando (2020) W. X. Zhou and A. Ariando, Jpn. J. Appl. Phys. 59, SI0802 (2020), URL https://doi.org/10.35848/1347-4065/ab8bbf.
  * Bhowal and Spaldin (2023) S. Bhowal and N. A. Spaldin, Annual Review of Materials Research 53, null (2023), eprint https://doi.org/10.1146/annurev-matsci-080921-105501, URL https://doi.org/10.1146/annurev-matsci-080921-105501.
  * Hickox-Young et al. (2023) D. Hickox-Young, D. Puggioni, and J. M. Rondinelli, Phys. Rev. Mater. 7, 010301 (2023), URL https://link.aps.org/doi/10.1103/PhysRevMaterials.7.010301.
  * Ma et al. (2019) Q. Ma, S.-Y. Xu, H. Shen, D. MacNeill, V. Fatemi, T.-R. Chang, A. M. Mier Valdivia, S. Wu, Z. Du, C.-H. Hsu, et al., Nature 565, 337 (2019), URL https://doi.org/10.1038/s41586-018-0807-6.
  * Kang et al. (2019) K. Kang, T. Li, E. Sohn, J. Shan, and K. F. Mak, Nat. Mater. 18, 324 (2019), URL https://doi.org/10.1038/s41563-019-0294-7.
  * Xiao et al. (2020) R.-C. Xiao, D.-F. Shao, W. Huang, and H. Jiang, Phys. Rev. B 102, 024109 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.102.024109.
  * Levitov et al. (1984) L. S. Levitov, Y. V. Nazarov, and G. M. Éliashberg, Zh. Eksp. Teor. Fiz 88, 229 (1984), URL http://jetp.ras.ru/cgi-bin/dn/e_061_01_0133.pdf.
  * Yoda et al. (2015) T. Yoda, T. Yokoyama, and S. Murakami, Sci. Rep. 5, 12024 (2015), URL https://doi.org/10.1038/srep12024.
  * Zhong et al. (2016) S. Zhong, J. E. Moore, and I. Souza, Phys. Rev. Lett. 116, 077201 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.077201.
  * Şahin et al. (2018) C. Şahin, J. Rou, J. Ma, and D. A. Pesin, Phys. Rev. B 97, 205206 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.205206.
  * Tsirkin et al. (2018) S. S. Tsirkin, P. A. Puente, and I. Souza, Phys. Rev. B 97, 035158 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.035158.
  * Sodemann and Fu (2015) I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.115.216806.
  * (17) The gyrotropic point groups are those crystal point groups that allow the presence of both polar ($\vec{P}$) and axial ($\vec{L}$) vector components transforming according to equivalent representations of the point group. The 18 gyrotropic point groups, therefore, allow for a second-rank tensor, connecting physical quantities $P_{i}$ and $L_{j}$, that describes physical effects such as natural optical activity.
  * He and Law (2020) W.-Y. He and K. T. Law, Phys. Rev. Res. 2, 012073 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.012073.
  * Bhowal et al. (2022a) S. Bhowal, S. P. Collins, and N. A. Spaldin, Phys. Rev. Lett. 128, 116402 (2022a), URL https://link.aps.org/doi/10.1103/PhysRevLett.128.116402.
  * Nelmes and Kuhs (1985) R. Nelmes and W. Kuhs, Solid State Commun. 54, 721 (1985), ISSN 0038-1098, URL https://www.sciencedirect.com/science/article/pii/0038109885905952.
  * Gu et al. (2017) J.-x. Gu, K.-j. Jin, C. Ma, Q.-h. Zhang, L. Gu, C. Ge, J.-s. Wang, C. Wang, H.-z. Guo, and G.-z. Yang, Phys. Rev. B 96, 165206 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.165206.
  * Iijima et al. (2000) T. Iijima, H. Näfe, and F. Aldinger, Integr. Ferroelectr. 30, 9 (2000), eprint https://doi.org/10.1080/10584580008222248, URL https://doi.org/10.1080/10584580008222248.
  * Giannozzi et al. (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al., J. Phys. Condens. Matter. 21, 395502 (2009), URL https://dx.doi.org/10.1088/0953-8984/21/39/395502.
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997), URL https://link.aps.org/doi/10.1103/PhysRevB.56.12847.
  * Souza et al. (2001) I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001), URL https://link.aps.org/doi/10.1103/PhysRevB.65.035109.
  * Mostofi et al. (2014) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 185, 2309 (2014), ISSN 0010-4655, URL https://www.sciencedirect.com/science/article/pii/S001046551400157X.
  * He et al. (2020) W.-Y. He, D. Goldhaber-Gordon, and K. T. Law, Nat. Commun. 11, 1650 (2020), URL https://doi.org/10.1038/s41467-020-15473-9.
  * Bhowal and Satpathy (2020) S. Bhowal and S. Satpathy, Phys. Rev. B 102, 201403 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.102.201403.
  * Hara et al. (2020) D. Hara, M. S. Bahramy, and S. Murakami, Phys. Rev. B 102, 184404 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.102.184404.
  * Bhowal et al. (2022b) S. Bhowal, S. P. Collins, and N. A. Spaldin, Phys. Rev. Lett. 128, 116402 (2022b), URL https://link.aps.org/doi/10.1103/PhysRevLett.128.116402.
  * Xiao et al. (2005) D. Xiao, J. Shi, and Q. Niu, Phys. Rev. Lett. 95, 137204 (2005), URL https://link.aps.org/doi/10.1103/PhysRevLett.95.137204.
  * Thonhauser et al. (2005) T. Thonhauser, D. Ceresoli, D. Vanderbilt, and R. Resta, Phys. Rev. Lett. 95, 137205 (2005), URL https://link.aps.org/doi/10.1103/PhysRevLett.95.137205.
  * Ceresoli et al. (2006) D. Ceresoli, T. Thonhauser, D. Vanderbilt, and R. Resta, Phys. Rev. B 74, 024408 (2006), URL https://link.aps.org/doi/10.1103/PhysRevB.74.024408.
  * Shi et al. (2007) J. Shi, G. Vignale, D. Xiao, and Q. Niu, Phys. Rev. Lett. 99, 197202 (2007), URL https://link.aps.org/doi/10.1103/PhysRevLett.99.197202.
  * Lopez et al. (2012) M. G. Lopez, D. Vanderbilt, T. Thonhauser, and I. Souza, Phys. Rev. B 85, 014435 (2012), URL https://link.aps.org/doi/10.1103/PhysRevB.85.014435.
  * Thouless et al. (1982) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982), URL https://link.aps.org/doi/10.1103/PhysRevLett.49.405.
  * Andersen and Saha-Dasgupta (2000a) O. K. Andersen and T. Saha-Dasgupta, Phys. Rev. B 62, R16219 (2000a), URL https://link.aps.org/doi/10.1103/PhysRevB.62.R16219.
  * Go et al. (2017) D. Go, J.-P. Hanke, P. M. Buhl, F. Freimuth, G. Bihlmayer, H.-W. Lee, Y. Mokrousov, and S. Blügel, Sci. Rep. 7, 46742 (2017), URL https://doi.org/10.1038/srep46742.
  * Go et al. (2021) D. Go, D. Jo, T. Gao, K. Ando, S. Blügel, H.-W. Lee, and Y. Mokrousov, Phys. Rev. B 103, L121113 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.L121113.
  * Arras et al. (2019) R. Arras, J. Gosteau, H. J. Zhao, C. Paillard, Y. Yang, and L. Bellaiche, Phys. Rev. B 100, 174415 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.174415.
  * Spaldin et al. (2013) N. A. Spaldin, M. Fechner, E. Bousquet, A. Balatsky, and L. Nordström, Phys. Rev. B 88, 094429 (2013), URL https://link.aps.org/doi/10.1103/PhysRevB.88.094429.
  * Löwdin (2004) P. Löwdin, The Journal of Chemical Physics 19, 1396 (2004), ISSN 0021-9606, eprint https://pubs.aip.org/aip/jcp/article-pdf/19/11/1396/11158843/1396_1_online.pdf, URL https://doi.org/10.1063/1.1748067.
  * Andersen and Saha-Dasgupta (2000b) O. K. Andersen and T. Saha-Dasgupta, Phys. Rev. B 62, R16219 (2000b), URL https://link.aps.org/doi/10.1103/PhysRevB.62.R16219.
  * Hayami et al. (2018) S. Hayami, M. Yatsushiro, Y. Yanagi, and H. Kusunose, Phys. Rev. B 98, 165110 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.165110.
  * Bhowal and Satpathy (2019) S. Bhowal and S. Satpathy, npj Computational Materials 5, 61 (2019), URL https://doi.org/10.1038/s41524-019-0198-8.
  * Johansson et al. (2018) A. Johansson, J. Henk, and I. Mertig, Phys. Rev. B 97, 085417 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.085417.
  * Roy and Narayan (2022) S. Roy and A. Narayan, J. Condens. Matter Phys. 34, 385301 (2022), URL https://dx.doi.org/10.1088/1361-648X/ac8091.


Thus, after considering the case in which both classes pick out the most
possible subsets of $D$ simultaneously, we have
$s(\widetilde{\mathcal{C}_{1}},n)\leq 4((n-1)^{d-1}+1)^{2}.$

This same argument can be repeated to show
$s(\widetilde{\mathcal{C}_{2}},n)\leq 4((n-1)^{d-1}+1)^{2},$ where

$\displaystyle\widetilde{\mathcal{C}}_{2}$
$\displaystyle=\\{\\{(\boldsymbol{x},y)\in\mathbb{R}^{d+1}:(1,\boldsymbol{x}^{\top})\boldsymbol{u}<0,\;\textup{sign}(y-(1,\boldsymbol{x}^{\top})\boldsymbol{\theta})\leq
0\\}:\boldsymbol{u}\neq\boldsymbol{0},\;\boldsymbol{\theta}\in\Theta\\}.$

Since each element of $\mathcal{C}_{R}$ is given by
$\widetilde{C}_{1}\cup\widetilde{C}_{2},$ where
$\widetilde{C}_{1}\in\widetilde{\mathcal{C}}_{1}$ and
$\widetilde{\mathcal{C}}_{1},$ a similar argument that was made when defining
classes using intersections results in the bound, $s(\mathcal{C}_{R},n)\leq
16((n-1)^{d-1}+1)^{4}.$ ∎

Note that this bound on the shattering coefficient also implies
$\sup_{\boldsymbol{\theta}}\lvert\textup{rdepth}(\boldsymbol{\theta},D)/n-\inf_{\boldsymbol{u}\neq\boldsymbol{0}}G(\boldsymbol{\theta},\boldsymbol{u})\rvert$
converges almost surely to zero at a rate of $O(1/\sqrt{n}),$ which is already
a known result in the literature; see for example, Gao, (2020). Lemma 13 also
implies our required probability bound, which is provided below.

Figure 4. The mean squared error at the 25% and 75% percentile of $x_{i}$ of
each DP estimator, for each DGP $x_{i},v_{i}\sim N(\mu=0,\sigma^{2}=1)$ (top),
$x_{i},v_{i}\sim\textup{Laplace}(\mu=0,b=1)$ (middle), $x_{i},v_{i}\sim
t(\mu=0,\nu=3)$ (bottom), and for $\epsilon=4$ (left) and $\epsilon=12$
(right). For NoisyStats$(\cdot)$ and the DPMedTS_exp$(\cdot)$ mechanisms,
$\psi$ was set to 0.95. Both NoisyStats$(\cdot)$ and DPMedTS_exp$(\cdot)$
satisfy $\epsilon-$DP; DPGradDesc$(\cdot),$ DPMedsweep$(\cdot),$ and
DPDeepestReg$(\cdot)$ satisfy $(\epsilon,10^{-6})-$DP.

Figure 5. The median absolute error at the 25% and 75% percentile of $x_{i}$
of each DP estimator, for each DGP $x_{i},v_{i}\sim N(\mu=0,\sigma^{2}=1)$
(top), $x_{i},v_{i}\sim\textup{Laplace}(\mu=0,b=1)$ (middle), $x_{i},v_{i}\sim
t(\mu=0,\nu=3)$ (bottom), and for $\epsilon=4$ (left) and $\epsilon=12$
(right). For NoisyStats$(\cdot)$ and the DPMedTS_exp$(\cdot)$ mechanisms,
$\psi$ was set to 0.95. Both NoisyStats$(\cdot)$ and DPMedTS_exp$(\cdot)$
satisfy $\epsilon-$DP; DPGradDesc$(\cdot),$ DPMedsweep$(\cdot),$ and
DPDeepestReg$(\cdot)$ satisfy $(\epsilon,10^{-6})-$DP.

######  Lemma 14.

Suppose $D$ is composed of observations that were sampled independently from
the population distribution $T.$ Then, for $\kappa>0,$ the following bound
holds

(39) $\displaystyle
P\left(\sup_{\boldsymbol{\theta}\in\Theta,\;\boldsymbol{u}\neq\boldsymbol{0}}\lvert
G(\boldsymbol{\theta},\boldsymbol{u})-G_{n}(\boldsymbol{\theta},\boldsymbol{u})\rvert\geq\kappa\right)\leq
64\left((n^{2}-1)^{d-1}+1\right)^{4}\exp\left(2\kappa(2+(2-n)\kappa)\right).$

###### Proof.

This a direct consequence of Lemma 13 and the main result of (Devroye,, 1982);
see Theorem 12 for more detail. ∎

######  Remark 5.

The remark following Theorem 12, regarding the rate of convergence of the
probability bounds provided above being suboptimal, are also applicable in
this case.

##  Appendix C Simulation Results For Different Choices of $\psi$

As described previously, we set the bounds on each $x_{i}$ and $y_{i}$ in our
simulations for both the NoisyStats$(\cdot)$ and the DPMedTS_exp$(\cdot)$
mechanisms by finding the smallest bounding box of the form $[-c,c]^{2}$ such
that a proportion $\psi\in(0,1)$ of the data points are in $[-c,c]^{2}.$ The
simulations provided above used $\psi=0.98,$ and, since this choice did have a
significant impact on the simulation results, this Appendix C provides the
simulation results for the alternative choices of $\psi=0.95$ and $\psi=1.$
Specifically, Figures 4 and 5 provide the mean squared error and median
absolute error when $\psi=0.95,$ and Figures 6 and 7 provide these error
metrics for the case in which $\psi=1.$

For the case in which $\psi=0.95,$ all of the comments regarding the relative
performance of the DP estimators from Section 5 appear to still be the case;
however, in almost all cases the performance appears to improve when
$\psi=0.95,$ relative to when $\psi=0.98,$ by between 50 and 200%, with the
largest improvements for cases in which $\epsilon=12$ and for the DGPs with
thicker tails. In contrast, for the cases in which $\psi=1,$ the statistical
efficiency appears to change more noticeably; unsurprisingly, this is most
apparent when $\epsilon$ is small and for both of the DGPs with thicker tails.

Figure 6. The mean squared error at the 25% and 75% percentile of $x_{i}$ of
each DP estimator, for each DGP $x_{i},v_{i}\sim N(\mu=0,\sigma^{2}=1)$ (top),
$x_{i},v_{i}\sim\textup{Laplace}(\mu=0,b=1)$ (middle), $x_{i},v_{i}\sim
t(\mu=0,\nu=3)$ (bottom), and for $\epsilon=4$ (left) and $\epsilon=12$
(right). For NoisyStats$(\cdot)$ and the DPMedTS_exp$(\cdot)$ mechanisms,
$\psi$ was set to 1. Both NoisyStats$(\cdot)$ and DPMedTS_exp$(\cdot)$ satisfy
$\epsilon-$DP; DPGradDesc$(\cdot),$ DPMedsweep$(\cdot),$ and
DPDeepestReg$(\cdot)$ satisfy $(\epsilon,10^{-6})-$DP.

Figure 7. The median absolute error at the 25% and 75% percentile of $x_{i}$
of each DP estimator, for each DGP $x_{i},v_{i}\sim N(\mu=0,\sigma^{2}=1)$
(top), $x_{i},v_{i}\sim\textup{Laplace}(\mu=0,b=1)$ (middle), $x_{i},v_{i}\sim
t(\mu=0,\nu=3)$ (bottom), and for $\epsilon=4$ (left) and $\epsilon=12$
(right). For NoisyStats$(\cdot)$ and the DPMedTS_exp$(\cdot)$ mechanisms,
$\psi$ was set to 1. Both NoisyStats$(\cdot)$ and DPMedTS_exp$(\cdot)$ satisfy
$\epsilon-$DP; DPGradDesc$(\cdot),$ DPMedsweep$(\cdot),$ and
DPDeepestReg$(\cdot)$ satisfy $(\epsilon,10^{-6})-$DP.

## References

  * Agarwal, (1990) Agarwal, P. K. (1990).  Partitioning arrangements of lines I: An efficient deterministic algorithm.  Discrete & Computational Geometry, 5(5):449–483.
  * Alabi et al., (2020) Alabi, D., McMillan, A., Sarathy, J., Smith, A., and Vadhan, S. (2020).  Differentially private simple linear regression.  arXiv preprint arXiv:2007.05157.
  * Amenta et al., (2000) Amenta, N., Bern, M., Eppstein, D., and H Teng, S. (2000).  Regression depth and center points.  Discrete & Computational Geometry, 23(3):305–323.
  * Awan and Slavković, (2021) Awan, J. and Slavković, A. (2021).  Structure and sensitivity in differential privacy: Comparing k-norm mechanisms.  Journal of the American Statistical Association, 116(534):935–954.
  * Bai and He, (2008) Bai, Z.-D. and He, X. (2008).  Asymptotic distributions of the maximal depth estimators for regression and mulitvariate location.  In Advances In Statistics, pages 241–262. World Scientific.
  * Barrientos et al., (2019) Barrientos, A. F., Reiter, J. P., Machanavajjhala, A., and Chen, Y. (2019).  Differentially private significance tests for regression coefficients.  Journal of Computational and Graphical Statistics, 28(2):440–453.
  * Bassily et al., (2014) Bassily, R., Smith, A., and Thakurta, A. (2014).  Private empirical risk minimization: Efficient algorithms and tight error bounds.  In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pages 464–473. IEEE.
  * Bun and Steinke, (2016) Bun, M. and Steinke, T. (2016).  Concentrated differential privacy: Simplifications, extensions, and lower bounds.  In Theory of Cryptography Conference, pages 635–658. Springer.
  * Bun and Steinke, (2019) Bun, M. and Steinke, T. (2019).  Average-case averages: Private algorithms for smooth sensitivity and mean estimation.  Advances in Neural Information Processing Systems, 32.
  * Chen et al., (2016) Chen, Y., Machanavajjhala, A., Reiter, J. P., and Barrientos, A. F. (2016).  Differentially private regression diagnostics.  In ICDM, pages 81–90.
  * Chetty et al., (2018) Chetty, R., Friedman, J. N., Hendren, N., Jones, M. R., and Porter, S. R. (2018).  The opportunity atlas: Mapping the childhood roots of social mobility.  Technical report, National Bureau of Economic Research.
  * Devroye, (1982) Devroye, L. (1982).  Bounds for the uniform deviation of empirical measures.  Journal of Multivariate Analysis, 12(1):72–79.
  * Devroye et al., (2013) Devroye, L., Györfi, L., and Lugosi, G. (2013).  A probabilistic theory of pattern recognition, volume 31.  Springer Science & Business Media.
  * Donoho and Gasko, (1992) Donoho, D. L. and Gasko, M. (1992).  Breakdown properties of location estimates based on halfspace depth and projected outlyingness.  The Annals of Statistics, pages 1803–1827.
  * Dvoretzky et al., (1956) Dvoretzky, A., Kiefer, J., and Wolfowitz, J. (1956).  Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator.  The Annals of Mathematical Statistics, pages 642–669.
  * (16) Dwork, C., Kenthapadi, K., McSherry, F., Mironov, I., and Naor, M. (2006a).  Our data, ourselves: Privacy via distributed noise generation.  In Annual International Conference on the Theory and Applications of Cryptographic Techniques, pages 486–503. Springer.
  * Dwork and Lei, (2009) Dwork, C. and Lei, J. (2009).  Differential privacy and robust statistics.  In Proceedings of the forty-first annual ACM symposium on Theory of computing, pages 371–380.
  * (18) Dwork, C., McSherry, F., Nissim, K., and Smith, A. (2006b).  Calibrating noise to sensitivity in private data analysis.  In Theory of cryptography conference, pages 265–284. Springer.
  * Dwork and Roth, (2014) Dwork, C. and Roth, A. (2014).  The algorithmic foundations of differential privacy.  Theoretical Computer Science, 9(3-4):211–407.
  * Dwork et al., (2014) Dwork, C., Talwar, K., Thakurta, A., and Zhang, L. (2014).  Analyze gauss: optimal bounds for privacy-preserving principal component analysis.  In Proceedings of the forty-sixth annual ACM symposium on Theory of computing, pages 11–20.
  * Edelsbrunner, (1987) Edelsbrunner, H. (1987).  Algorithms in combinatorial geometry, volume 10.  Springer Science & Business Media.
  * Foulds et al., (2016) Foulds, J., Geumlek, J., Welling, M., and Chaudhuri, K. (2016).  On the theory and practice of privacy-preserving bayesian data analysis.  In Proceedings of the Thirty-Second Conference on Uncertainty in Artificial Intelligence, pages 192–201.
  * Gao, (2020) Gao, C. (2020).  Robust regression via mutivariate regression depth.  Bernoulli, 26(2):1139–1170.
  * Hall et al., (2013) Hall, R., Wasserman, L., and Rinaldo, A. (2013).  Random differential privacy.  Journal of Privacy and Confidentiality, 4(2).
  * He and Portnoy, (1998) He, X. and Portnoy, S. (1998).  Asymptotics of the deepest line.  In Applied Statistical Science III: Nonparametric Statistics and related Topics. Citeseer.
  * Kamath et al., (2022) Kamath, G., Mouzakis, A., Singhal, V., Steinke, T., and Ullman, J. (2022).  A private and computationally-efficient estimator for unbounded gaussians.  In Conference on Learning Theory, pages 544–572. PMLR.
  * Karwa and Vadhan, (2018) Karwa, V. and Vadhan, S. (2018).  Finite sample differentially private confidence intervals.  9th Innovations in Theoretical Computer Science Conference (ITCS 2018).
  * Kifer and Machanavajjhala, (2014) Kifer, D. and Machanavajjhala, A. (2014).  Pufferfish: A framework for mathematical privacy definitions.  ACM Transactions on Database Systems (TODS), 39(1):1–36.
  * Kifer et al., (2012) Kifer, D., Smith, A., and Thakurta, A. (2012).  Private convex empirical risk minimization and high-dimensional regression.  In Conference on Learning Theory, pages 25–1. JMLR Workshop and Conference Proceedings.
  * Koenker, (2000) Koenker, R. (2000).  Galton, Edgeworth, Frisch, and prospects for quantile regression in econometrics.  Journal of Econometrics, 95(2):347–374.
  * Koenker and Bassett Jr, (1978) Koenker, R. and Bassett Jr, G. (1978).  Regression quantiles.  Econometrica: journal of the Econometric Society, pages 33–50.
  * (32) Langerman, S. and Steiger, W. (2003a).  The complexity of hyperplane depth in the plane.  Discrete & computational geometry, 30(2):299–309.
  * (33) Langerman, S. and Steiger, W. (2003b).  Optimization in arrangements.  In Annual Symposium on Theoretical Aspects of Computer Science, pages 50–61. Springer.
  * Langerman and Steiger, (2000) Langerman, S. and Steiger, W. L. (2000).  An optimal algorithm for hyperplane depth in the plane.  In SODA, pages 54–59.
  * Masse, (2002) Masse, J.-C. (2002).  Asymptotics for the tukey median.  Journal of Multivariate Analysis, 81(2):286–300.
  * McSherry and Mironov, (2009) McSherry, F. and Mironov, I. (2009).  Differentially private recommender systems: Building privacy into the netflix prize contenders.  In Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 627–636.
  * McSherry and Talwar, (2007) McSherry, F. and Talwar, K. (2007).  Mechanism design via differential privacy.  In 48th Annual IEEE Symposium on Foundations of Computer Science (FOCS’07), pages 94–103. IEEE.
  * Miller et al., (2001) Miller, K., Ramaswami, S., Rousseeuw, P., Sellares, T., Souvaine, D., Streinu, I., and Struyf, A. (2001).  Fast implementation of depth contours using topological sweep.  In SODA, pages 690–699.
  * Mizera, (2002) Mizera, I. (2002).  On depth and deep points: a calculus.  The Annals of Statistics, 30(6):1681–1736.
  * Nissim et al., (2007) Nissim, K., Raskhodnikova, S., and Smith, A. (2007).  Smooth sensitivity and sampling in private data analysis.  In Proceedings of the thirty-ninth annual ACM symposium on Theory of computing, pages 75–84.
  * Ramsay and Chenouri, (2021) Ramsay, K. and Chenouri, S. (2021).  Differentially private depth functions and their associated medians.  arXiv preprint arXiv:2101.02800.
  * Rousseeuw and Hubert, (1999) Rousseeuw, P. J. and Hubert, M. (1999).  Regression depth.  Journal of the American Statistical Association, 94(446):388–402.
  * Rousseeuw and Struyf, (1998) Rousseeuw, P. J. and Struyf, A. (1998).  Computing location depth and regression depth in higher dimensions.  Statistics and Computing, 8(3):193–203.
  * Segaert et al., (2017) Segaert, P., Hubert, M., and Rousseeuw, P. (2017).  mrfdepth: Depth measures in multivariate, regression and functional settings.  https://cran.r-project.org/package=mrfDepth.
  * Sen, (1968) Sen, P. K. (1968).  Estimates of the regression coefficient based on kendall’s tau.  Journal of the American statistical association, 63(324):1379–1389.
  * Sheffet, (2017) Sheffet, O. (2017).  Differentially private ordinary least squares.  In International Conference on Machine Learning, pages 3105–3114. PMLR.
  * Small, (1990) Small, C. G. (1990).  A survey of multidimensional medians.  International Statistical Review/Revue Internationale de Statistique, pages 263–277.
  * Struyf and Rousseeuw, (2000) Struyf, A. and Rousseeuw, P. J. (2000).  High-dimensional computation of the deepest location.  Computational Statistics & Data Analysis, 34(4):415–426.
  * Talagrand, (1994) Talagrand, M. (1994).  Sharper bounds for gaussian and empirical processes.  The Annals of Probability, pages 28–76.
  * Theil, (1950) Theil, H. (1950).  A rank-invariant method of linear and polynomial regression analysis.  Indagationes mathematicae, 12(85):173.
  * Topsøe, (2004) Topsøe, F. (2004).  Some bounds for the logarithmic function.  RGMIA Res. Rep. Collection, 7(2):1–20.
  * Tukey, (1975) Tukey, J. W. (1975).  Mathematics and the picturing of data.  In Proceedings of the International Congress of Mathematicians, Vancouver, 1975, volume 2, pages 523–531.
  * Van Aelst et al., (2002) Van Aelst, S., Rousseeuw, P. J., Hubert, M., and Struyf, A. (2002).  The deepest regression method.  Journal of Multivariate Analysis, 81(1):138–166.
  * van Kreveld et al., (2008) van Kreveld, M., Mitchell, J. S., Rousseeuw, P., Sharir, M., Snoeyink, J., and Speckmann, B. (2008).  Efficient algorithms for maximum regression depth.  Discrete & computational geometry, 39(4):656–677.
  * Vapnik and Chervonenkis, (2015) Vapnik, V. N. and Chervonenkis, A. Y. (2015).  On the uniform convergence of relative frequencies of events to their probabilities.  In Measures of complexity, pages 11–30. Springer.
  * Vu and Slavkovic, (2009) Vu, D. and Slavkovic, A. (2009).  Differential privacy for clinical trial data: Preliminary evaluations.  In 2009 IEEE International Conference on Data Mining Workshops, pages 138–143. IEEE.
  * Wang, (2018) Wang, Y.-X. (2018).  Revisiting differentially private linear regression: optimal and adaptive prediction & estimation in unbounded domain.  arXiv preprint arXiv:1803.02596.
  * Wasserman and Zhou, (2010) Wasserman, L. and Zhou, S. (2010).  A statistical framework for differential privacy.  Journal of the American Statistical Association, 105(489):375–389.
  * Wei and Dudley, (2012) Wei, F. and Dudley, R. M. (2012).  Two-sample Dvoretzky–Kiefer–Wolfowitz inequalities.  Statistics & Probability Letters, 82(3):636–644.