# Magnetic response of nematic superconductors: skyrmion stripes and their
signatures in muon spin relaxation experiments

Martin Speight  School of Mathematics, University of Leeds, Leeds LS2 9JT,
United Kingdom  Thomas Winyard  School of Mathematics, University of Leeds,
Leeds LS2 9JT, United Kingdom  School of Mathematics, Statistics, and
Actuarial Science, University of Kent, Canterbury, United Kingdom  Egor Babaev
Department of Physics, KTH-Royal Institute of Technology, Stockholm, SE-10691
Sweden

###### Abstract

We investigate the magnetic response of nematic superconductors, presenting a
new approach to find vortex and skyrmion structures beyond symmetry-
constraining anzatses. Using this approach we show that nematic
superconductors form distinctive skyrmion stripes. Our approach lends itself
to accurately determining the field distribution for muon spin rotation
probes. We use this to show that the skyrmion structure manifests as a double
peak in the field distribution, markedly different from the signal of standard
vortex lattices.

Nematic states are currently one of the topics of central interest in
superconductivity Fu and Berg (2010); Fu (2014); Venderbos et al. (2016a);
Zyuzin et al. (2017); Wu and Martin (2017); Chirolli (2018); Barkman et al.
(2019); Uematsu et al. (2019); Dentelski et al. (2020); How and Yip (2020a);
Chichinadze et al. (2020). One particular material that has been the target of
intense experimental study is the doped topological insulator
$\mathrm{Bi_{2}Se_{3}}$ Hor et al. (2010); Wray et al. (2010); Kriener et al.
(2011). Superconducting candidate materials of interest include
$\mathrm{Cu_{x}Bi_{2}Se_{3}}$, $\mathrm{Sr_{x}Bi_{2}Se_{3}}$ Shruti et al.
(2015); Liu et al. (2015); Pan et al. (2016); Matano et al. (2016); Tao et al.
(2018) and $\mathrm{Nb_{x}Bi_{2}Se_{3}}$ Asaba et al. (2017). One of the
intriguing properties of nematic superconductors is that the flux carrying
objects are skyrmions Zyuzin et al. (2017); How and Yip (2020b); Speight et
al. (2022), which are fundamentally distinct from the vortices that normally
appear in superconductivity.

The appearance of skyrmions raises two questions: how do the skyrmions affect
the magnetic response of the material, and what unique experimental signatures
does this produce? The first question is technically nontrivial, since the
inter-skyrmion interactions are complicated and anisotropic, driven by mode
mixing in general multicomponent anistropic models Silaev et al. (2018);
Winyard et al. (2019a, b); Speight et al. (2019, 2022), which dictates the
appearance of multiple scales in the inter-skyrmion forces. This anisotropic
behaviour cannot be captured by a simple ansatz, thus it is challenging to
approximate the resulting lattice with any certainty. Using brute-force
numerics, on the other hand, would require simulating complicated Ginzburg-
Landau models at large system sizes, with boundaries that may affect structure
formation.

In this work and in the companion extended paper Speight et al. (2022) we
study this magnetic response by developing and applying a new method that can
determine vortex lattice structure for any model. We find that nematic
superconductors exhibit an unusual magnetic response in the form of skyrmion
stripes.

Since surface probes visualize vortex structures at limited scales and
structure formation may be affected by pinning, strain,inhomogeneities and
surface effects, the most reliable method to detect these structures is muon
spin rotation ($\mu SR$) experiments. These fire muons into the sample,
detecting the orientation of positrons emitted when the muons decay. This data
can then be used to calculate the statistical distribution of the magnetic
field in the bulk of a superconductor Sonier et al. (2000).

We will demonstrate that for nematic superconductors, this magnetic field
distribution exhibits a double peak structure. We posit that the number of
peaks in the distribution correlates with the number of length scales in the
lattice. For a nematic superconductor this is two, the separation between
chains and the inter-skyrmion separation along the chain, which we will call
the chain link length.

The model we consider in this paper is derived from the proposal that certain
topological insulators, such as electron-doped nematic superconductors (e.g.
$Cu_{x}Bi_{2}Se_{3}$) are described by a two-component order parameter Hor et
al. (2010); Wray et al. (2010); Kriener et al. (2011); Fu and Berg (2010).
This order-parameter breaks rotational symmetry in the basal plane to a mixed
symmetry, combining the phase difference of the order-parameter and spatial
rotations Speight et al. (2022). This matches the observed two fold rotational
symmetry of specific heat and upper critical field Venderbos et al. (2016b).
The associated Ginzburg-Landau Gibbs free energy for such a nematic
superconductor was microscopically derived in Zyuzin et al. (2017). In
dimensionless units, it reads

$G=\int_{\Omega}\left\\{\frac{1}{2}Q^{\alpha\beta}_{ij}\overline{D_{i}\psi_{\alpha}}D_{j}\psi_{\beta}+\frac{1}{2}|B-H|^{2}+F_{p}(|\psi_{\alpha}|,\theta_{\alpha\beta})\right\\},$
(1)

where $\psi_{\alpha}=\rho_{\alpha}e^{i\varphi_{\alpha}}$ are two complex
superconducting order parameters $\alpha\in\\{1,2\\}$. Note that Greek indices
will always be used for superconducting components, Latin indices will be used
for spatial directions and repeated indices imply summation. $F_{p}$ denotes
the collected potential terms and varies between models, but we will use the
standard form,

$F_{p}=\frac{\eta^{2}}{2}\left(-|\psi_{1}|^{2}-|\psi_{2}|^{2}+\frac{1}{2}(|\psi_{1}|^{4}+|\psi_{2}|^{4})+\gamma|\psi_{1}|^{2}|\psi_{2}|^{2}\right).$
(2)

The constant matrices $Q^{\alpha\beta}$ define the anisotropy in the system
and must obey the symmetry
$Q^{\alpha\beta}_{ij}=\overline{Q}^{\beta\alpha}_{ji}$, ensuring that the
energy is real. $D_{i}=\partial_{i}-iA_{i}$ is the covariant derivative, where
$A$ is the vector potential corresponding to a $U(1)$ gauge field, with
corresponding magnetic field $B=\nabla\times A$. For a nematic superconductor
the anisotropy matrices take the form

$\displaystyle Q^{11}=Q^{22}=\left(\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\
0&0&\beta_{z}\end{array}\right),\quad
Q^{12}=\beta_{\perp}\left(\begin{array}[]{ccc}1&i&0\\\ i&-1&0\\\
0&0&0\end{array}\right).$ (9)

The field equations in the bulk of the superconductor can then be found by
varying the Gibbs free energy w.r.t. the fields $\psi_{\alpha}$ and $A_{i}$.

Solutions for isolated flux-carrying objects show that in nematic
superconductors these are skyrmions, that have the form of spatially separated
half-quantum vortices Zyuzin et al. (2017); How and Yip (2020b); Speight et
al. (2022). These solutions were found by simulating the field equations on a
plane ($\Omega=\mathbb{R}^{2}$), however to find a lattice solution, we must
solve them over the unit cell of a lattice with periodic boundary conditions
(up to gauge). While it is often assumed that the unit cell has hexagonal
symmetry, there is no reason to expect this for skyrmions in an anisotropic
model. Hence we take a new approach, considering a general unit cell, which
takes the form of a parallelogram that is periodic (up to gauge), as shown in
fig. 1. A lattice solution is then given by the unit cell and field
configuration that minimises $G/\mathcal{A}$, where $\mathcal{A}$ is the area
of the unit cell. Hence numerically, we seek minimisers of $G/\mathcal{A}$
w.r.t. the fields and w.r.t. the geometry of the unit cell itself, represented
as the two vectors $V_{1}$ and $V_{2}$ in fig. 1 (note that this will also
control the area $\mathcal{A}$). The full technical detail of the procedure we
propose, is given in a companion paper Speight et al. (2022) (along with a
generalisation that is not of use in this paper). A brief account is also
presented in the appendix.

Figure 1:  Diagram of the geometry of a general unit cell for a vortex
lattice, defined by the two vectors $V_{1}$ and $V_{2}$ with angle $\alpha$.

By tessellating the solution on the periodic unit cell, we find a solution
that covers the plane with a vortex lattice. Hence, in applying the above
method we find the field configuration that minimises $G/\mathcal{A}$ in the
bulk of a superconductor with no boundary effects. In addition, as we know the
shape of the unit cell, we know the symmetry of the lattice.

Applying the above approach for the basal plane of a nematic system gives the
results in fig. 2, where we see that for low and medium strength external
field $H$, the unit cell is rectangular ($\alpha=\pi/2$), containing four
half-quantum vortices (totalling winding two) in each cell, forming a Skyrmion
chain. The lattice vectors lie parallel to the chain and perpendicular to the
chain. Hence, we can understand such a nematic unit cell by considering two
lengths:

  * •

$|V_{1}|$ \- chain link length,

  * •

$|V_{2}|$ \- the distance between the chains.

We can see how these values change in relation to an increasing external field
strength in fig. 2. As the strength of the external field $H$ is initially
increased, the density of chains increases and $|V_{2}|$ decreases. Then we
see that $|V_{1}|$ also decreases as the chain segments are squashed.
Eventually when the external field $H$ becomes large and approaches
$H_{c_{2}}$, the vortices are forced together and form the familiar tightly
packed triangular lattice.

Figure 2:  Contour plots of the lattice solutions (skyrmion chains) for the
model eq. 1 in the Basal plane, with $\eta=3$ and increasing external field
$H$. The fields plotted are the order parameter magnitudes $|\psi_{1}|$ and
$|\psi_{2}|$, the phase difference between the two complex order parameters
$\theta_{12}=\varphi_{1}-\varphi_{2}$ and the local magnetic field strength
$B$. The lattice consists of skyrmion chains, formed by rectangular unit cells
(marked with black lines) containing four half-quantum vortices or winding
two. As the external field increases the skyrmion chains squash together. Then
near $H_{c_{2}}$ a transition occurs and the vortices are forced into the
conventional triangular lattice with unit cells of winding one. Note that
there is no simple symmetry between $|\psi_{1}|$ and $|\psi_{2}|$ as the
symmetry is broken to a mixed symmetry with spatial rotations Speight et al.
(2022).

Theoretically in clean samples these skyrmion structures should be visible on
the surface using scanning SQUID or scanning Hall probes. However, flux
structure near the surface can be affected by surface physics and it is
difficult to probe large areas. Hence, one of the most powerful tool to
identify skyrmions is the muon spin relaxation technique, which gives
statistical information regarding the flux distribution in the bulk of the
system. To compare the lattice configurations in fig. 2 with $\mu$SR
experiment, we must convert the numerical field configuration for a given unit
cell to a probability density distribution $p(B)$. Our method is ideal for
this as it produces the field configuration on a single unit cell, with no
boundary effects and no deformations. We will use a standard kernel density
estimation Parzen (1962) to approximate a continuous distribution from our
discrete field configuration. Namely, our numerical algorithm produces a
discrete configuration, where the field distribution is a set of normalised
delta functions. We choose to smooth these delta functions by replacing them
with localised Gaussians and taking the normalised sum. Hence, the probability
density of field strength $B$ in the unit cell is,

$p(B)=\frac{1}{n_{1}n_{2}h}\sum^{n_{1}}_{i=1}\sum_{j=1}^{n_{2}}\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{B-B_{ij}}{h})^{2}},$
(10)

where $B_{ij}$ is the local magnetic field at lattice site $(i,j)$ with
$i\in[1,n_{1}]$ and $j\in[1,n_{2}]$. The parameter $h$ determines the amount
of smoothing and should be set relative to the number of grid points used to
simulate the lattice. We used $h=0.01$ throughout.

We used the above method to construct the field distribution $p(B)$ for a
nematic superconductor and a standard single component superconductor for
comparison. The latter’s parameters were chosen so that $H_{c_{1}}$ and
$H_{c_{2}}$ match those of the nematic model for $\eta=3$. This process and
the parameters are described in the appendix. We have plotted the nematic
field distribution for field strengths nearer $H_{c1}$ than $H_{c2}$ in fig.
3, and also compared the results to the standard GL model. There is a clear
difference, with the nematic superconductor exhibiting a double peak structure
as opposed to the well known single peak of the standard GL model. It is also
different from the $\mu SR$ signature of a system that phase separates into
vortex clusters and domains of Meissner states Biswas et al. (2020); Ray et
al. (2014).

Figure 3:  Magnetic field distribution for a skyrmion state in nematic
superconductor (top) and vortex lattice in single component superconductor
(bottom) for external fields closer to $H_{c1}$ than $H_{c2}$. The nematic
superconductor is described in eq. 1 with $\eta=3$. The structure of the field
distribution is markedly different for the skyrmion state in the nematic
superconductor with a double peak, rather than the familiar single peak for
the standard GL model.

Figure 4:  Magnetic field distribution for a nematic superconductor (top) and
single component superconductor (bottom) for external fields $H_{c2}\approx
H\gg H_{c1}$. The nematic superconductor is described in eq. 1 with $\eta=3$.
As one approaches $H_{c2}$, the structure of the field distribution for the
nematic superconductor becomes similar to the single peak structure of the
familiar single peak for the standard GL model. This corresponds with the
lattice becoming a triangular lattice for external fields near $H_{c2}$.

Consider the single component Ginzburg-Landau model. Here, the neighbouring
vortex separation is determined by a single parameter, the Ginzburg-Landau
parameter $\kappa$. While this parameter is all that is needed to describe the
lattice solution, it is actually the ratio of two length scales. The form of
the magnetic field distribution is determined by how the magnetic field decays
between neighbouring vortices, hence as there is a single defining parameter,
there is a single corresponding peak in the field distribution. In contrast,
multicomponent anisotropic models have a continuous set of coupled length
scales that are directionally dependent Speight et al. (2019, 2022). However,
given the rectangular nematic unit cell shown in fig. 2, we can describe the
vortex separation in the unit cell using just two length scales, $|V_{1}|$ and
$|V_{2}|$. Hence, the decay of the magnetic field between two neighbouring
vortices depends on whether they are in the same chain or neighbouring chains.
If the vortices are in the same chain the decay is determined by $|V_{1}|$, if
they are in different chains the decay is determined by $|V_{2}|$. Hence, each
of these neighbouring decay possibilities, or length scales, contributes a
different peak to the field distribution $p(B)$. This can be seen in fig. 5,
where we have compared the probability density distribution $p(B)$ with the
corresponding magnetic field plot $B(x)$ for $H=2$. It can be seen in this
figure that the first sharp peak corresponds to chain separation, and the
second peak corresponds to the chain link length.

The two peaks for low fields start far apart, due to the large chain
separation $|V_{2}|$. Then as the field strength increases, the chains are
packed in much tighter and $|V_{2}|$ decreases, causing the two peaks to get
closer. Then as the external field strength is increased further, both
$|V_{1}|$ and $|V_{2}|$ shrink. This explains why the peaks begin far apart
for particularly low fields. Hence, we ultimately end up with a very
distinctive double peak structure for nematic skyrmion lattices. We expect
that in the presence of weak pinning the statistics of the magnetic fields in
the bulk will still retain these features.

Figure 5:  Comparing the field distribution $p(B)$ of a nematic superconductor
in the lower plot with the areas of the local magnetic field configuration
$B(x)$ each peak corresponds to for external field strength $H=2$.

While the above analysis holds for low and medium strength external fields,
once the external field is increased to be near $H_{c2}$, the vortices are
forced into a tightly packed triangular lattice. This triangular structure
leads to the most efficient packing for the vortices, hence it is natural that
the increased pressure will eventually force a triangular lattice to be lower
in energy. The vortex separation for the nematic superconductor is now
described by a single length scale and forms a single peak, similar to that of
the single component model. This can be seen for strong external fields in
fig. 4, where we have again compared the nematic and single component results.

In conclusion we have shown that the magnetic response of nematic
superconductors takes the form of skyrmion stripes, which is markedly
different from the standard vortex lattice. As the field is increased, these
skyrmion stripes decrease their separation. Eventually the magnetization
process compresses the skyrmion stripes themselves, forming the conventional
triangular vortex lattice close to upper critical field $H_{c2}$. We also find
that skyrmion stripes give qualitatively distinct signatures in muon spin
rotation experiments. In particular, we observe a distinctive two peak
structure, where the second peak changes relative to the first as the external
field strength is increased. These findings can be used to experimentally
confirm (i) nematic superconducting states and (ii) that flux-carrying objects
in superconductors take the form of skyrmions.

## Acknowledgements

We thank Vadim Grinenko for discussions. EB was supported by the Swedish
Research Council Grants 2016-06122, 2018-03659, and Olle Engkvists Stiftelse.
MS and TW were supported by the UK Engineering and Physical Sciences Research
Council, through grant EP/P024688/1. TW was also supported by an academic
development fellowship, awarded by the University of Leeds. The numerical work
of this paper was performed using the code library Soliton Solver, developed
by one of the authors, and was undertaken on ARC4, part of the High
Performance Computing facilities at the University of Leeds.

##  Appendix A Finding vortex lattice solutions

In this section we give a brief account of the new method used to find vortex
lattice solutions of arbitrary complexity in the presence of an external field
$H$. Detailed discussion of the method is presented in the companion paper
Speight et al. (2022). Namely, we find the unit cell of the periodic vortex
solution in the bulk, with external field strength
$H_{c_{1}}<H_{0}<H_{c_{2}}$. The standard approach for a Ginzburg-Landau model
is to consider a unit cell of degree $n=1$ with either triangular
($\alpha=\pi/3$) or square ($\alpha=\pi/2$) symmetry. However, as the nematic
system we are considering is anisotropic, there is no reason to expect that a
lattice with such high symmetry will be the global minimizer. Hence, the
correct approach is to minimize energy, not just w.r.t. the periodic fields,
but also w.r.t. the geometry of the unit cell itself. We present here a
general method to find the optimal unit cell, without assuming the symmetry of
the underlying lattice, based upon the method first presented in Speight et
al. (2022).

We first assume that far from the boundary of the system (deep in the bulk),
the gauge invariant quantities $(\rho_{\alpha},\theta_{12},B)$ are periodic on
the unit cell. We also note that we can represent a general periodic structure
as a tessellation of a general parallelogram (unit cell), as seen in figure 1,
formed by two vectors $V_{1}$ and $V_{2}$ with angle $\alpha$.

We then seek local minimisers of the Gibbs free energy in eq. 1 w.r.t. the
fields $\psi_{\alpha}(x),A_{i}(x)$, now defined on a general flat periodic
2-torus $T_{\Lambda}^{2}=\mathbb{R}^{2}/\Lambda$, where the period lattice
$\Lambda$ is spanned by vectors $V_{1}$, $V_{2}$ in fig. 1,

$\Lambda=\left\\{n_{1}V_{1}+n_{2}V_{2}\,|\,n_{1},n_{2}\in\mathbb{Z},V_{1},V_{2}\in\mathbb{R}^{2}\right\\}.$
(11)

The unit cell has periodic boundary conditions, up to gauge, so that the
fields have boundary conditions,

$\displaystyle\psi_{\alpha}(x+V_{i})$
$\displaystyle=\psi_{\alpha}(x)e^{if_{i}},$ $\displaystyle A_{j}(x+V_{i})$
$\displaystyle=A_{j}(x)+\partial_{j}f_{i}$ (12)

Note that while the fields are not purely periodic, the above boundary
conditions do leave all physical quantities periodic as desired,

$\displaystyle B(x+V_{i})=B(x),\quad\,J(x+V_{i})=J(x),$
$\displaystyle|\psi(x+V_{i})|=|\psi(x)|.\quad\quad\quad\quad$ (13)

To simulate the fields on the unit cell, we will simplify the above
formulation by transforming to a more convenient coordinate system in the
$x,y$ plane. Let us define $X_{1},X_{2}$ so that
$(x,y)=X_{1}V_{1}+X_{2}V_{2}$. The unit cell spanned by $V_{1},V_{2}$ is now
covered by $(X_{1},X_{2})\in[0,1]\times[0,1]$. Let $L$ be the matrix with
columns $V_{1},V_{2}$. Note that $\det L=\mathcal{A}$, the area of the unit
cell. Now define the unimodular matrix $M=\sqrt{\mathcal{A}}L^{-1}$ such that
$X_{i}=\frac{M_{ij}}{\sqrt{\mathcal{A}}}x_{j}$. Then

$\displaystyle F(M,\mathcal{A},\psi_{\alpha},A)$
$\displaystyle=\int_{[0,1]^{2}}\left\\{\frac{1}{2}(MQ^{\alpha\beta}M^{T})_{ij}\overline{D_{i}\psi_{\alpha}}D_{j}\psi_{\beta}+\frac{1}{2\mathcal{A}}\left((\partial_{1}A_{2}-\partial_{2}A_{1})^{2}\right)+\mathcal{A}F_{p}(\psi_{\alpha})\right\\}dX_{1}\,dX_{2},$

leading in turn to the Gibbs free energy,

$\displaystyle G(M,\mathcal{A},\psi_{\alpha},A)$
$\displaystyle=F-\int_{[0,1]^{2}}H_{0}(\partial_{1}A_{2}-\partial_{2}A_{1})\,dX_{1}\,dX_{2}$
$\displaystyle=F-2n\pi H_{0}.$ (14)

where $n$ is the winding number of the field configuration. Note that this is
only valid in the basal plane, as otherwise one must consider non-zero $A_{3}$
in the free energy. This would allow in-plane spontaneous fields and magnetic
field twisting, which occurs away from the basal plane due to the coupled
length scalesSpeight et al. (2022). Models that exhibit such fields have been
studied in Speight et al. (2021); Benfenati et al. (2020); Speight et al.
(2022).

If we exclude boundary effects, the optimal lattice solution is the one that
minimises the the total Gibbs free energy of the system. If the system area is
$\mathcal{A}_{sys}$, this can be calculated from the unit cell as
$G_{sys}=\mathcal{A}_{sys}G/\mathcal{A}$. Hence we seek minimisers of
$G/\mathcal{A}$ or Gibbs free energy per unit area, with respect to the
fields, geometry of the unit cell $M$ and area $\mathcal{A}$. Note that the
degree or winding number of a given unit cell $n$ is fixed and strictly
speaking we must the find the global minimiser of all $G_{n}/\mathcal{A}$,
where $G_{n}$ is the Gibbs free energy of a unit cell of degree $n$ and then
find the infimum of this set.

In practice, it is sufficient to find minima of $G_{n}/\mathcal{A}$ until we
get a repeated minimiser, that is, until we find a minimiser of
$G_{n}/\mathcal{A}$ whose cell and field configuration is two cells of the
$G_{n/2}/\mathcal{A}$ minimiser joined together.

###  A.1 Numerical method

To find minimisers of $G(M,\mathcal{A},\psi_{\alpha},A)$ we discretise the
standard square unit cell on an $n_{1}\times n_{2}$ regular grid. We utilize
an arrested Newton flow algorithm subject to the periodic boundary conditions,

$\displaystyle\psi_{\alpha}(x+(1,0))$ $\displaystyle=\psi_{\alpha}(x)e^{in2\pi
nx_{2}},$ $\displaystyle\psi_{\alpha}(x+(0,1))$
$\displaystyle=\psi_{\alpha}(x),$ $\displaystyle A_{2}(x+(1,0))$
$\displaystyle=A_{1}(x)+2\pi n$ $\displaystyle A_{i}(x+(1,0))$
$\displaystyle=A_{2}(x)\quad i\neq 2$ $\displaystyle A_{i}(x+(0,1))$
$\displaystyle=A_{i}(x)$ (15)

where $n$ is the winding number of the unit cell. We set the fields
$\phi=(\psi,A)$ on the unit square torus and the geometry of the unit cell
$(M,\mathcal{A})$ to be some initial condition. Fixing the unit cell
$(M,\mathcal{A})$, we find a local minimum w.r.t. the collected fields $\phi$,
using arrested Newton flow for a notional particle subject to the potential
$\ddot{\phi}=-grad\,G_{dis}(M,\mathcal{A},\phi)$. A more detailed description
of arrested newton flow is given in Speight and Winyard (2020). This is
continued until the discretised function is minimised with respect to some
tolerance. We then fix the field configuration $\phi$ and area $\mathcal{A}$,
and minimise $G_{dis}(M,\mathcal{A},\phi)/\mathcal{A}$ with respect to $M\in
SL(2,\mathbb{R})$, where we note that only the gradient term of the energy is
dependent on this matrix. This can be done exactly using linear algebra,
described in the next subsection. Finally, we fix the fields $\phi$ and the
lattice shape $M$ and minimise the free energy w.r.t. $\mathcal{A}$, by
solving the algebraic equation in $\mathcal{A}$ with coefficients given by the
integrals of the different energy terms.

###  A.2 Finding the minimal geometry

We can find the minimal geometry of a unit cell with a given configuration
$(\psi_{\alpha},A)$ explicitly. We first note that the only term that is
dependent on the geometry $M$ is the gradient term of the free energy. Hence
we write the gradient energy as a $4\times 4$ complex matrix,

$P^{c}_{ik,jl}:=\int_{[0,1]^{2}}Q^{\alpha\beta}_{kl}\overline{D_{i}\psi_{\alpha}}D_{j}\psi_{\beta},$
(16)

where we take each index pair $(i,k),(j,l)$ to range over the values
$\left\\{(1,1),(1,2),(2,1),(2,2)\right\\}$. Due to the symmetry
$Q^{\alpha\beta}_{ij}=\overline{Q^{\beta\alpha}_{ji}}$ we know $P^{c}$ is
hermitian. Hence if we split it into its real and imaginary parts
$P^{c}=P+iP_{I}$ we know that the real part $P$ is symmetric and the imaginary
part $P_{I}$ is skew. If we then consider the vector,

$X:=\left(M_{11}M_{12}M_{21}M_{22}\right)^{T},$ (17)

we can rewrite the gradient part of the free energy as,

$F_{grad}(X)=\frac{1}{2}X^{T}P^{c}X=\frac{1}{2}X^{T}PX.$ (18)

Hence we must minimise the quadratic form $X^{T}PX$ subject to the constraint
that $\det M=1$ (or $M\in SL(2,\mathbb{R}$), which can be written,

$\frac{1}{2}X^{T}JX=1,\quad\quad J:=\left(\begin{array}[]{cccc}0&0&0&1\\\
0&0&-1&0\\\ 0&-1&0&0\\\ 1&0&0&0\end{array}\right).$ (19)

Hence, we seek critical points of
$F_{\text{grad}}:SL(2,\mathbb{R})\rightarrow\mathbb{R}$, which are
$X\in\mathbb{R}^{4}$ with $X^{T}JX=1$ and,

$PX=\lambda\left(X_{4},-X_{3},-X_{2},X_{1}\right)^{T}=\lambda JX,$ (20)

for some $\lambda\in\mathbb{R}$. We can see that $J^{2}=\mathbb{I}_{4}$ is the
identity matrix, hence eq. 20 becomes,

$JPX=\lambda X.$ (21)

In other words $X$ is an eigenvector of $JP$.

This allows us to minimise $F_{\text{grad}}$ w.r.t. $M$ explicitly by
constructing $JP$ and finding it’s 4 eigenvectors $Y$. We discard any complex
solutions and any s.t. $Y^{T}JY\leq 0$. We compute the solution we seek
$X=Y/\sqrt{Y^{T}JY}$ such that,

$F_{\text{grad}}=\frac{1}{2}X^{T}PX=\frac{1}{2}X^{T}\lambda JX=\lambda.$ (22)

We can then read off the new period lattice as,

$L=\sqrt{\mathcal{A}}M^{-1},\quad\quad
M=\left(\begin{array}[]{cc}X_{1}&X_{2}\\\ X_{3}&X_{4}\end{array}\right).$ (23)

The period lattice is spanned by the columns of $L$, named $V_{1}$, $V_{2}$.

##  Appendix B Single component field distribution

It is instructive to compare the results for nematic systems to a standard
single component Ginzburg Landau model. We will consider the following model,

$G=\int_{\mathbb{R}^{2}}\left\\{\frac{1}{2}\overline{D_{i}\psi}D_{i}\psi+\frac{1}{2}|B-H|^{2}+F_{p}(|\psi|)\right\\},$
(24)

where there is a single order parameter $\psi$, and

$F_{p}(|\psi|)=-a|\psi|^{2}+\frac{b}{2}|\psi|^{4}.$ (25)

To effectively compare this model to the results for the nematic model, we
will select the parameters $a,b$ so that $H_{c_{1}}$ and $H_{c_{2}}$ equal
those of the nematic model. By a standard computation, for example see Tinkham
page 134 Tinkham (2004),

$H_{c_{2}}=2a,$ (26)

while, by definition, $H_{c1}$ is the value of external field strength $H$ for
which the Gibbs energy of a vortex solution equals that of the superconducting
ground state, yielding

$H_{c1}=\frac{F_{\scriptscriptstyle vortex}-F_{\scriptscriptstyle
ground}}{2\pi}.$ (27)

The right hand side of eq. 27 depends implicitly on $a,b$. We find that
$a=5.70$, $b=2.21$ produces $H_{c1}=1.18$ and $H_{c2}=11.41$, matching the
nematic model with $\eta=3$.

## References

  * Fu and Berg (2010) L. Fu and E. Berg, Physical review letters 105, 097001 (2010).
  * Fu (2014) L. Fu, Physical Review B 90, 100509 (2014).
  * Venderbos et al. (2016a) J. W. Venderbos, V. Kozii, and L. Fu, Physical Review B 94, 180504 (2016a).
  * Zyuzin et al. (2017) A. Zyuzin, J. Garaud, and E. Babaev, Physical review letters 119, 167001 (2017).
  * Wu and Martin (2017) F. Wu and I. Martin, Physical Review B 95, 224503 (2017).
  * Chirolli (2018) L. Chirolli, Physical Review B 98, 014505 (2018).
  * Barkman et al. (2019) M. Barkman, A. A. Zyuzin, and E. Babaev, Physical Review B 99, 220508 (2019).
  * Uematsu et al. (2019) H. Uematsu, T. Mizushima, A. Tsuruta, S. Fujimoto, and J. Sauls, Physical review letters 123, 237001 (2019).
  * Dentelski et al. (2020) D. Dentelski, V. Kozii, and J. Ruhman, Physical Review Research 2, 033302 (2020).
  * How and Yip (2020a) P. T. How and S.-K. Yip, Phys. Rev. Research 2, 043192 (2020a), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.043192.
  * Chichinadze et al. (2020) D. V. Chichinadze, L. Classen, and A. V. Chubukov, Phys. Rev. B 102, 125120 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.102.125120.
  * Hor et al. (2010) Y. S. Hor, A. J. Williams, J. G. Checkelsky, P. Roushan, J. Seo, Q. Xu, H. W. Zandbergen, A. Yazdani, N. P. Ong, and R. J. Cava, Physical review letters 104, 057001 (2010).
  * Wray et al. (2010) L. A. Wray, S.-Y. Xu, Y. Xia, Y. San Hor, D. Qian, A. V. Fedorov, H. Lin, A. Bansil, R. J. Cava, and M. Z. Hasan, Nature Physics 6, 855 (2010).
  * Kriener et al. (2011) M. Kriener, K. Segawa, Z. Ren, S. Sasaki, and Y. Ando, Physical Review Letters 106, 127004 (2011).
  * Shruti et al. (2015) Shruti, V. K. Maurya, P. Neha, P. Srivastava, and S. Patnaik, Phys. Rev. B 92, 020506 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.92.020506.
  * Liu et al. (2015) Z. Liu, X. Yao, J. Shao, M. Zuo, L. Pi, S. Tan, C. Zhang, and Y. Zhang, Journal of the American Chemical Society 137, 10512 (2015), URL http://dx.doi.org/10.1021/jacs.5b06815.
  * Pan et al. (2016) Y. Pan, A. M. Nikitin, G. K. Araizi, Y. K. Huang, Y. Matsushita, T. Naka, and A. de Visser, Scientific Reports 6, 28632 (2016), URL http://dx.doi.org/10.1038/srep28632.
  * Matano et al. (2016) K. Matano, M. Kriener, K. Segawa, Y. Ando, and G.-q. Zheng, Nat Phys 12, 852 (2016), ISSN 1745-2473, URL http://dx.doi.org/10.1038/nphys3781.
  * Tao et al. (2018) R. Tao, Y.-J. Yan, X. Liu, Z.-W. Wang, Y. Ando, Q.-H. Wang, T. Zhang, and D.-L. Feng, Physical Review X 8, 041024 (2018).
  * Asaba et al. (2017) T. Asaba, B. J. Lawson, C. Tinsman, L. Chen, P. Corbae, G. Li, Y. Qiu, Y. S. Hor, L. Fu, and L. Li, Phys. Rev. X 7, 011009 (2017), URL https://link.aps.org/doi/10.1103/PhysRevX.7.011009.
  * How and Yip (2020b) P. T. How and S.-K. Yip, Physical Review Research 2, 043192 (2020b).
  * Speight et al. (2022) M. Speight, T. Winyard, and E. Babaev, _Symmetries, length scales, magnetic response and skyrmion chains in nematic superconductors_ (2022), eprint 2202.13674.
  * Silaev et al. (2018) M. Silaev, T. Winyard, and E. Babaev, Physical Review B 97, 174504 (2018).
  * Winyard et al. (2019a) T. Winyard, M. Silaev, and E. Babaev, Physical Review B 99, 064509 (2019a).
  * Winyard et al. (2019b) T. Winyard, M. Silaev, and E. Babaev, Physical Review B 99, 024501 (2019b).
  * Speight et al. (2019) M. Speight, T. Winyard, and E. Babaev, Physical Review B 100, 174514 (2019).
  * Sonier et al. (2000) J. E. Sonier, J. H. Brewer, and R. F. Kiefl, Reviews of Modern Physics 72, 769 (2000).
  * Venderbos et al. (2016b) J. W. Venderbos, V. Kozii, and L. Fu, Physical Review B 94, 094522 (2016b).
  * Parzen (1962) E. Parzen, The annals of mathematical statistics 33, 1065 (1962).
  * Biswas et al. (2020) P. Biswas, F. N. Rybakov, R. Singh, S. Mukherjee, N. Parzyk, G. Balakrishnan, M. Lees, C. Dewhurst, E. Babaev, A. Hillier, et al., Physical Review B 102, 144523 (2020).
  * Ray et al. (2014) S. Ray, A. Gibbs, S. Bending, P. Curran, E. Babaev, C. Baines, A. Mackenzie, and S. Lee, Physical Review B 89, 094504 (2014).
  * Speight et al. (2021) M. Speight, T. Winyard, A. Wormald, and E. Babaev, Physical Review B 104, 174515 (2021).
  * Benfenati et al. (2020) A. Benfenati, M. Barkman, T. Winyard, A. Wormald, M. Speight, and E. Babaev, Physical Review B 101, 054507 (2020).
  * Speight and Winyard (2020) M. Speight and T. Winyard, Physical Review B 101, 134420 (2020).
  * Tinkham (2004) M. Tinkham, _Introduction to superconductivity_ (Courier Corporation, 2004).


#  Magnetic Compton profile in non-magnetic ferroelectrics

Sayantika Bhowal  Materials Theory, ETH Zurich, Wolfgang-Pauli-Strasse 27,
8093 Zurich, Switzerland  Stephen P. Collins  Diamond Light Source Ltd,
Diamond House, Harwell Science & Innovation Campus, Didcot, Oxfordshire, OX11
0DE  Nicola A. Spaldin  Materials Theory, ETH Zurich, Wolfgang-Pauli-Strasse
27, 8093 Zurich, Switzerland

###### Abstract

Magnetic Compton scattering is an established tool for probing magnetism in
ferromagnetic or ferrimagnetic materials with a net spin polarization. Here we
show that, counterintuitively, non-magnetic systems can also have a non-zero
magnetic Compton profile, provided that space-inversion symmetry is broken.
The magnetic Compton profile is antisymmetric in momentum and, if the
inversion symmetry is broken by an electric-field switchable ferroelectric
distortion, can be reversed using an electric field. We show that the
underlying physics of the magnetic Compton profile and its electrical control
are conveniently described in terms of $k$-space magnetoelectric multipoles,
which are reciprocal to the real-space charge dipoles associated with the
broken inversion symmetry. Using the prototypical ferroelectric lead titanate,
PbTiO3, as an example, we show that the ferroelectric polarization introduces
a spin asymmetry in momentum space that corresponds to a pure $k$-space
magnetoelectric toroidal moment. This in turn manifests in an antisymmetric
magnetic Compton profile which can be reversed using an electric field. Our
work suggests an experimental route to directly measuring and tuning hidden
$k$-space magnetoelectric multipoles via their magnetic Compton profile.

The inelastic Compton scattering of x-ray photons by electrons was an early
confirmation of quantum mechanical behavior, in its revolutionary step of
assigning momentum to electromagnetic waves Compton (1923). The effect is
widely used today in fields as diverse as radio-biology Redler  _et al._
(2018), astrophysics Wil (1985), and condensed matter, where it is used to
measure the distribution of electron density in momentum space Cooper (1971).
The extension to measuring the spin-dependent momentum distribution, now known
as magnetic Compton scattering, was proposed as early as 1970 Platzman and
Tzoar (1970), and dramatic advances in synchrotron light sources mean that the
weak scattering cross section in the spin channel of materials can now be
captured.

Over the past decades, magnetic Compton scattering Cooper  _et al._ (2007);
Duffy (2013); Ahuja (2013) has been used effectively to investigate the
magnetic properties of ferro/ferrimagnetic materials. The measured quantity,
the magnetic Compton profile (MCP), $J_{mag}(p_{z})$, is the one-dimensional
projection of the spin-polarized electron momentum density which, for
backscattering, lies along the direction of the scattering vector $p_{z}$,

$J_{mag}(p_{z})=\int\int[\rho^{\uparrow}(\vec{p})-\rho^{\downarrow}(\vec{p})]dp_{x}dp_{y}\quad.$
(1)

Here $\rho^{\uparrow}(\vec{p})$ $[\rho^{\downarrow}(\vec{p})]$ is the momentum
density of the majority [minority] spin bands. The MCP is insensitive to
orbital magnetism due to the instantaneous interaction of the photon with the
electrons Timms  _et al._ (1993); Cooper  _et al._ (1992), and therefore
allows separation of the spin and orbital contributions to the total magnetic
moment Cooper  _et al._ (1991); Duffy  _et al._ (2010); Itou  _et al._ (2013),
as well as the magnetic contributions from localized and itinerant electrons
Zukowski  _et al._ (1993); Duffy  _et al._ (1998, 2000); Banfield  _et al._
(2005); Shenton-Taylor  _et al._ (2007). It is also used to measure the spin
polarization at the Fermi surface Duffy (2013); Mijnarends  _et al._ (2007);
Mizoroki  _et al._ (2011), which is at the heart of spintronic applications.
Since a non-zero MCP was believed to require broken time-reversal ($\cal T$)
symmetry, it has to date only been studied in magnetic materials.

Here we show that a non-zero MCP can exist in non-magnetic (NM) materials
provided that space-inversion ($\cal I$) symmetry is broken. The MCP in this
case is antisymmetric, since the $\cal T$ symmetry dictates
$J_{mag}(p_{z})\overset{\cal T}{\rightarrow}-J_{mag}(-p_{z})$. This behavior
is in striking contrast to the centrosymmetric, NM case with both $\cal I$ and
$\cal T$ symmetries, for which the MCP is zero since
$J_{mag}(p_{z})\overset{\cal I}{\rightarrow}J_{mag}(-p_{z})$.

These simple symmetry relations indicate the possibility of a non-zero MCP for
a NM system, and the underlying physics can be elegantly described in terms of
the odd-parity electric multipoles induced by the broken $\cal I$ symmetry,
which determine both the occurrence of a non-zero $J_{mag}$ and its direction
in momentum space. The key physics is illustrated in Fig. 1, showing the
duality between real and $k$ space Watanabe and Yanase (2018); Bhowal and
Spaldin (2021) in odd parity ($\cal I$ asymmetric) materials, between the
real- and momentum-space multipole representations in magnetic and NM species.
As seen from the top panel of Fig. 1 (a), the real-space ME multipole (MEM)
(antisymmetric under both $\cal I$ and $\cal T$) Ederer and Spaldin (2007);
Spaldin  _et al._ (2008, 2013) corresponds to an electric dipole (no
magnetization dependence) in reciprocal space, which gives rise to an
antisymmetric part in the regular Compton profile Bhowal and Spaldin (2021).
Here we point out that an exactly opposite situation exists for a real-space
electric dipole (bottom panel of Fig. 1 (a)), which corresponds to a MEM in
momentum space and, in turn, gives rise to an antisymmetric MCP.

The transformation of the electric dipole in real space into a MEM in
$k$-space follows from the fact that both $r$ and $k$ change sign upon
inversion, i.e. $\vec{r}\overset{\cal I}{\rightarrow}-\vec{r}$, and
$\vec{k}\overset{\cal I}{\rightarrow}-\vec{k}$, but behave differently under
time-reversal, with $\vec{r}\overset{\cal T}{\rightarrow}\vec{r}$, but
$\vec{k}\overset{\cal T}{\rightarrow}-\vec{k}$. As a result, broken inversion
symmetry in real space alone is sufficient to break both inversion and time-
reversal symmetries in momentum space, leading to $k$-space MEMs.

Figure 1: (a) Cross-links between real (left) and $k$-space (right) multipoles
in $\cal T$ asymmetric (top) and symmetric (bottom) odd-parity systems. (b)
Magnetic Compton scattering with a circularly polarized photon beam in PTO.
The solid-outlined Ti-atom position, $\mathfrak{t}^{(k)}_{z}$ spins and MCP
profile indicate one polarization orientation. With reversed E (dashed
outlines), the Ti ion displacement reverses, the toroidal moment
$\mathfrak{t}^{(k)}_{z}\rightarrow-\mathfrak{t}^{(k)}_{z}$, and the sign of
the MCP reverses.

The following implications of these $k$-space MEMs are key to our work: First,
since the odd-parity charge multipoles [bottom panel of Fig. 1 (a)], which
measure the asymmetry of the electron density in real space (left), create an
asymmetry in the magnetization density in $k$ space (right), they can be
probed using magnetic measurements in momentum space, such as magnetic Compton
scattering. Second, as the $k$-space MEMs are the measure of this momentum-
space asymmetry in the magnetization density, the MCP provides a direct
signature of the specific $k$-space MEMs. Third, since the $k$-space MEMs
originate from broken $\cal I$ symmetry, structural tuning that changes these
multipoles will modify the MCP. Such structural changes can be easily
accessible in ferroelectric materials, where the lattice distortions (electric
polarization) can be controlled by an external electric field. Such an
electrical tuning of the MCP typically allows faster switching than a
conventional magnetic field, leading to reduced systematic errors in
experiments.

Here, we illustrate the above ideas by explicitly computing the MCP and the
odd parity multipoles for the prototype tetragonal ferroelectric PbTiO3 (PTO),
using first-principles methods based on density functional theory (DFT) as
implemented in extended versions Ernsting  _et al._ (2014); Spaldin  _et al._
(2013) of the Elk code cod . The odd-parity charge multipoles are extracted by
decomposing the $\cal T$ symmetric density matrix
$\rho_{lm,l^{\prime}m^{\prime}}$ into parity-odd tensor moments, where only
the odd $l-l^{\prime}$ terms contribute to the desired multipoles Spaldin  _et
al._ (2013). To illustrate how different multipoles result in a different MCP,
we also introduce the case of the rhombohedral ferroelectric GeTe.

The key finding of our work is the antisymmetric MCP for NM ferroelectrics,
and its switching using an applied electric field (Fig. 1 (b)). We show that
odd-parity charge multipoles and spin-orbit interaction (SOI) are the key
ingredients for MCP in broken $\cal I$-symmetry systems. We reveal how the
structural asymmetry, governed by the odd-parity charge multipoles, introduces
additional inter-orbital hopping parameters. These cause the magnetic-moment
dependence of the charge multipoles in $k$ space, and the MEMs in momentum
space. An interesting outcome is the realization of a pure ME toroidal moment,
$\mathfrak{t}^{(k)}_{z}$, in the momentum space of PTO; such pure toroidal
moments are hard to realize in real space, where they are usually accompanied
by ME quadrupole moments Spaldin  _et al._ (2008, 2013); Thöle  _et al._
(2020); Spaldin (2021); Bhowal and Spaldin (2021). The magnetic asymmetry,
governed by these $k$-space MEMs, combined with SOI, results in the
antisymmetric MCP.

We start by analyzing the symmetry of ferroelectric PTO, which crystallizes in
the tetragonal $P4mm$ structure Nelmes and Kuhs (1985), with no ${\cal I}$
symmetry [middle panel of Fig. 1 (b)]. The absence of ${\cal I}$ in the
$C_{4v}$ point group allows for several odd-parity charge multipoles, of which
the $A_{1}^{+}$ irreducible representation (IR), that corresponds to the PTO
structure, allows for the odd parity $Q_{10}$ dipole as its lowest order
charge multipole (note that higher order multipoles of odd parity, such as the
$Q_{30}$ octupole, etc. are also allowed) as shown in Table I not ; Watanabe
and Yanase (2018); Inui  _et al._ (1990); Perez-Mato  _et al._ (2015). As
expected from the duality and cross-links discussed above, the corresponding
$k$-space representation has explicit dependence on the magnetic moment
($\vec{\mu}$) and represents a pure $k$-space toroidal moment
$\mathfrak{t}_{z}^{(k)}=k_{x}\mu_{y}-k_{y}\mu_{x}$ (this can be seen simply by
replacing $\vec{r}\rightarrow\vec{k}$), without any ME quadrupole moment
components. In fact, the $k$-space ME quadrupole moments belong to different
IRs (${\cal Q}_{xy}^{(k)}$ to $B_{1}^{+}$, etc). By definition, the existence
of $\mathfrak{t}_{z}^{(k)}$ indicates the presence of magnetically polarized
bands with magnetic orientation along $y$ ($x$) in the $k_{x}$ ($k_{y}$)
direction of momentum space, which, in turn, implies the presence of a MCP
along $p_{x}$ ($p_{y}$). Here, $x,y,z$ denote the Cartesian axes.

Table 1: The basis functions of the odd parity charge multipoles (lowest order in $k$) corresponding to the IR representations of PTO and GeTe. The $k$-space ME monopole ${\cal A}^{(k)}$, toroidal moment $\vec{\mathfrak{t}}^{(k)}$, and quadrupole moment ${\cal Q}_{ij}^{(k)}$ are also indicated within the parentheses. The magnetic moment $\vec{\mu}\equiv\mu_{B}\big{(}\frac{2\vec{l}}{l+1}+2\vec{s}\big{)}$ includes contributions from both spin $\vec{s}$ and orbital $\vec{l}$ angular momenta. IR | Multipole | Real space |  $k$ space
---|---|---|---
PTO ($C_{4v}$)
$A_{1}^{+}$ | $Q_{10}$ | $z$ |  $k_{x}\mu_{y}-k_{y}\mu_{x}$ ($\mathfrak{t}_{z}^{(k)}$)
GeTe ($C_{3v}$)
$A_{1}^{+}$ | $Q_{10}$ | $z$ |  $k_{x}\mu_{y}-k_{y}\mu_{x}$($\mathfrak{t}_{z}^{(k)}$)
$E^{+}$ | $\\{Q_{11}^{+},Q_{11}^{-}\\}$ | $\\{x,y\\}$ | $\\{k_{y}\mu_{z}+k_{z}\mu_{y},k_{z}\mu_{x}+k_{x}\mu_{z}\\}$
|  |  | ($\\{{\cal Q}_{yz}^{(k)},{\cal Q}_{xz}^{(k)}\\}$)
| $\\{Q_{31}^{+},Q_{31}^{-}\\}$ |  | $\\{k_{y}\mu_{z}-k_{z}\mu_{y},k_{z}\mu_{x}-k_{x}\mu_{z}\\}$
|  |  | ($\\{\mathfrak{t}_{x}^{(k)},\mathfrak{t}_{y}^{(k)}\\}$)
| $\\{Q_{32}^{+},Q_{32}^{-}\\}$ | $\\{(x^{2}-y^{2})z,$ | $\\{k_{x}\mu_{y}+k_{y}\mu_{x},k_{x}\mu_{x}-k_{y}\mu_{y}\\}$
|  | $xyz\\}$ | $(\\{{\cal Q}_{xy}^{(k)},{\cal Q}_{x^{2}-y^{2}}^{(k)}\\})$

Next, we explicitly compute the MCP of PTO Com  and show our result in Fig. 2
(a). We note first that the computed MCPs satisfy
$m=\int_{-\infty}^{\infty}J_{mag}(p_{x})dp_{x}=\int_{-\infty}^{\infty}J_{mag}(p_{y})dp_{y}=0$,
as required for a NM system, since they are purely antisymmetric functions.
The computed MCP is about three orders of magnitude smaller than the
corresponding total Compton profile and is about an order of magnitude smaller
than the MCP for ferromagnetic Ni Dixon  _et al._ (1998). The signs of the
MCPs along $p_{x}$ and $p_{y}$ are opposite, directly reflecting the opposite
signs of $k_{x}$ and $k_{y}$ in
$\mathfrak{t}_{z}^{(k)}=k_{x}\mu_{y}-k_{y}\mu_{x}$. These asymmetries are also
reflected in the spin asymmetry in the calculated bandstructure along the $\pm
k_{x}$ (or $\pm k_{y}$) direction of the momentum space, shown in Fig. 2 (b).
This spin asymmetry is also consistent with the presence of $\cal T$ symmetry,
which dictates that the $\uparrow$-spin band at $+\vec{k}$ has the same energy
as the $\downarrow$-spin band at $-\vec{k}$, that is
$\varepsilon_{n}(\vec{k}\uparrow)=\varepsilon_{n}(-\vec{k}\downarrow)$,
creating the left-right asymmetry of the spin-polarized bands in Fig. 2 (b).

Figure 2: (a) Calculated MCPs of PTO along the $p_{x}$ and $p_{y}$ directions
in momentum space. (b) Bandstructure of PTO, showing the same spin asymmetry
in the $S_{x}$ spin component, shown in color map, along $\pm k_{y}$. The
Fermi energy is set at zero eV. (c) Illustration of the antisymmetric hopping
$\pm\beta$ between $p_{x}$-$p_{z}$ and $p_{y}$-$p_{z}$ orbitals of the O atoms
along $\pm\hat{x}$ and $\pm\hat{y}$ directions respectively, induced by broken
$\cal I$ symmetry. Only hoppings between apical O atoms are indicated. In-
plane O atoms are shown in dashed circles. (d) Dependence of the MCP on the
enforced SOI strength $\lambda=\lambda_{c}\times\lambda_{r}$, with
$\lambda_{r}$ and $\lambda_{c}$ the actual SOI strength and the scaling factor
respectively.

Since the spin asymmetry that leads to the MCP is largest at the top of the
valence bands, which are predominantly of O-$p$ character, we next analyze the
role of the O $p$ orbitals in the chemical bonding, by analyzing the hopping
parameters between neighboring O $p$ orbitals along $\hat{x}$ and $\hat{y}$
directions. We note that the off-centering of the atoms in the non-
centrosymmetric ferroelectric structure induces additional interorbital
hoppings between the $p_{x}$-$p_{z}$ and $p_{y}$-$p_{z}$ orbitals of the O
atoms, which are forbidden in the Slater-Koster Slater and Koster (1954);
Harrison (1989) parameters in the presence of $\cal I$ symmetry [see Fig. 2
(c)]. Interestingly, the computed effective hopping parameters for the broken
symmetry structure, extracted using the NMTO downfolding method Andersen and
Saha-Dasgupta (2000), have left-right asymmetry (Fig. 2 (c)). Such
antisymmetric hopping parameters lead to sine terms in the tight-binding (TB)
model, in contrast to the usual cosine terms associated with symmetric
hopping, e.g. $\langle p_{x}|{\cal H_{A}}|p_{z}\rangle=-2i\beta\sin(k_{x}a)$.
In the small $k$ limit, $\sin(k_{x}a)\rightarrow k_{x}a$, and the
corresponding $p$-orbital TB Hamiltonian is ${\cal
H_{A}}(k)=\gamma(k_{x}\hat{L}_{y}-k_{y}\hat{L}_{x})$, where the constant
$\gamma\propto\beta$ is a measure of $\cal I$ asymmetry and
$\hat{L}_{x},\hat{L}_{y}$ are the Cartesian components of the orbital angular
momentum operator for the $p$ orbitals Sugano (1970). We recognize that ${\cal
H_{A}}$ represents the orbital part of the toroidal moment
$\mathfrak{t}^{(k)}_{z}$. In the presence of SOI, the orbitals couple with the
spins, giving rise to a spin component of $\mathfrak{t}^{(k)}_{z}$. Since the
MCP is insensitive to orbital magnetism, we see here the importance of SOI in
the MCP of PTO. To further clarify the role of SOI in MCP, we next
artificially change the strength of the coupling $\lambda$ in our
calculations. As seen from Fig. 2 (d), this results in a drastic change in the
magnitude of the MCP, with vanishing MCP at $\lambda=0$, consistent with the
hopping analysis.

Figure 3: (a) Variation of MCP with displacement $\delta$ (in units of lattice
vector $c$) of the Ti atom from the centrosymmetric position in PTO. The MCP
increases with $\delta$, and switches sign as the sign of $\delta$ changes.
(b) Variation of the odd-parity charge multipole $Q_{10}$ at Pb, Ti and O
sites [apical O(1) and in-plane O(2) and O(3)] for the same displacements
$\delta$ as in (a). The correlation between $Q_{10}$ and MCP is apparent. The
vertical dashed line corresponds to the centrosymmetric case with $\delta=0$,
at which both multipoles and MCP vanish. The odd-parity multipoles at Pb and
Ti sites are obtained by adding the $p$-$d$ and $p$-$s$ contributions, while
for the O atoms only the $p$-$s$ contribution is considered.

We now turn to the tuning of the MCP by manipulating the odd-parity MEMs
through changing the charge multipoles. For this, we artificially change the
ferroelectric displacements of the atoms in PTO along the polarization
direction $[001]$. The resulting changes in MCP are depicted in Fig. 3 (a). We
also study the corresponding variation in the charge multipoles at Pb, Ti and
the three O atoms, shown in Fig. 3 (b). As expected, the multipole at the Pb
atom dominates due to the $6s^{2}$ lone pair that causes the broken $\cal I$
symmetry. The magnitudes of the multipoles at the Ti and O atoms are rather
weak and have opposite sign to that of Pb atom. It is also apparent from Figs.
3 (a) and (b) that the MCP follows the trend in odd parity multipoles: The
magnitudes of both increase with increasing ferroelectric displacement and
both reverse sign as the displacement is switched. Since the switching of
ferroelectric displacements in PTO can be achieved using an external electric
field, our results show the possibility of switching MCP using electrical
means, which may be of practical importance in reducing the experimental
uncertainties in a measurement of the effect. Reciprocally, the MCP provides a
possible probe for the detection of odd parity electric multipoles.

To emphasize the correlation between MCP and odd-parity multipoles, we now
briefly analyze ferroelectric GeTe, with structural polarization along the
[111] direction of the rhombohedral unit cell (Fig. 4 (a)). The ferroelectric
distortion of GeTe corresponds to $A_{1}$ and $E$ IRs of its $C_{3v}$ point
group symmetry Chatterji  _et al._ (2018); Kagdada  _et al._ (2018). These
allow the $\mathfrak{t}_{x}^{(k)},\mathfrak{t}_{y}^{(k)}$ components of the
$k$-space toroidal moment, as well as quadrupole moment components ${\cal
Q}_{xy}^{(k)},{\cal Q}_{yz}^{(k)},{\cal Q}_{xz}^{(k)},$ and ${\cal
Q}_{x^{2}-y^{2}}^{(k)},$ in addition to $\mathfrak{t}_{z}^{(k)}$ [see Table
1], and the corresponding MCP has more components than in the case of PTO.
First, in contrast to PTO, the MCPs in GeTe have simultaneous contributions
from toroidal and quadrupole moments rather than pure toroidal contribution.
For example, $\\{\mathfrak{t}_{x}^{(k)},{\cal Q}_{yz}^{(k)}\\}$ contribute to
the MCP along $p_{y}$ with [001] spin quantization axis. Secondly, even along
the same momentum direction, the MCP can appear for different spin
quantization axes. For example, along the $p_{y}$ direction in momentum space,
an MCP exists for $\hat{x},\hat{y},$ and $\hat{z}$ spin quantization axes.
Knowledge of the odd-parity multipoles is invaluable in interpreting such a
complex MCP. We see that the MCP can again be switched by reversing the
direction of polarization [see Fig. 4 (b)]. While the magnitude of the MCP is
slightly smaller than that of PTO, in GeTe MCP components with parallel spin
and momentum directions [e.g., along $\hat{x}$ and $\hat{y}$ driven by ${\cal
Q}_{x^{2}-y^{2}}$ (see Table 1)] exist, allowing for experimental detection
within the simple back-scattering geometry.

Figure 4: Switching of MCP in GeTe. (a) The rhombohedral unit cell of GeTe.
The dashed atoms indicate the centrosymmetric Ge (0.25,0.25,0.25) and Te
(0.75,0.75,0.75) positions; the solid atoms the positions with polarization P.
(b) Calculated MCP for structures with polarization P, -P and 0 (for which the
MCP vanishes). Inset: odd-parity multipoles Q${}_{11}^{+}$ and Q${}_{31}^{+}$
(in a.u.), which contribute to $J_{mag}(p_{y})$, for the antisymmetric spin
density with spin quantization along [001]. Multipoles at Ge and Te have
opposite sign. Q${}_{11}^{+}$ is summed from its $p$-$d$ and $p$-$s$
contributions; Q${}_{31}^{+}$ has only $p$-$d$ contribution.

Finally, we discuss the experimental set up for detecting the MCP in PTO. The
experiment will be similar to conventional magnetic Compton scattering, with a
circularly polarized beam Cooper  _et al._ (2007), except that no magnetic
field is required. The antisymmetric spin profile can be obtained by reversing
the circular polarization of the incoming photon and subtracting the two
signals, similar to the case of antisymmetric Compton measurements Bhowal  _et
al._ , or, more conveniently, the electric polarization can be flipped rapidly
with fixed photon helicity by reversing the electric field. Sensitivity to the
MCP arises from a relativistic correction to the scattering cross section
(details are given in Ref. Cooper  _et al._ (2007)). For the ideal case of
backscattering, the relative sensitivity to the MCP, compared to to the NM
profile, is reduced by $\sim\frac{2E}{m_{e}c^{2}}$, where $E$ is the incident
photon energy. For 100 keV photons, this amounts to a relatively modest
reduction of $\sim$0.4.

It is important to note that we need a single ferroelectric domain of PTO or
GeTe for the magnetic Compton scattering measurements, which may be obtained
at room temperature (the ferroelectric transition temperatures for PTO and
GeTe being T${}_{\rm c}\sim 765$ K Samara (1971) and around 650 KChattopadhyay
_et al._ (1987); Wdowik  _et al._ (2014) respectively) using an electric
field. The switching of the MCP can be observed by reversing the direction of
the electric field, which switches the ferroelectric domain, and, hence, the
odd parity multipoles.

The odd parity real-space charge multipoles and associated $k$-space MEMs also
have implications in other physical effects beyond the MCP Hayami  _et al._
(2018); Watanabe and Yanase (2020); Onimaru  _et al._ (2011); Fu (2015);
Yatsushiro and Hayami (2020). Examples include an exotic superconducting state
driven by odd-parity multipole fluctuations Ishizuka and Yanase (2018); Sumita
and Yanase (2020), and the predicted Rashba effect in PTO Arras  _et al._
(2019), which is a direct consequence of $\mathfrak{t}_{z}^{(k)}$. Higher
order multipoles in PTO can give rise to a Berry curvature dipole which should
lead to a nonlinear Hall effect. Furthermore, the $k$-space MEMs could
describe the recently observed polar skyrmions Das  _et al._ (2019), in the
same spirit as their real-space counterpart forms a basis for describing
magnetic skyrmions Göbel  _et al._ (2019). We hope that our proposal
stimulates experimental efforts to measure the MCP, or these related MEM-
driven behaviors, in NM ferroelectrics.

## Acknowledgements

We thank Jon Duffy, Stephen Dugdale and Urs Staub 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

  * Compton (1923) A. H. Compton, Phys. Rev. 21, 483 (1923).
  * Redler  _et al._ (2018) G. Redler, K. C. Jones, A. Templeton, D. Bernard, J. Turian, and J. C. H. Chu, Med Phys. 45, 1233 (2018).
  * Wil (1985) _Radiative Processes in Astrophysics_ (John Wiley & Sons, Ltd, 1985) Chap. 7, pp. 195–223.
  * Cooper (1971) M. Cooper, Adv. Phys 20, 453 (1971).
  * Platzman and Tzoar (1970) P. M. Platzman and N. Tzoar, Phys. Rev. B 2, 3556 (1970).
  * Cooper  _et al._ (2007) M. Cooper, C. Shenton-Taylor, J. Duffy, C. Steer, and L. Blaauw, Nucl. Instrum. Methods Phys. Res. A: Accel. Spectrom. Detect. Assoc. Equip. 580, 1 (2007), proceedings of the 10 th International Symposium on Radiation Physics.
  * Duffy (2013) J. A. Duffy, J. Phys. Conf. Ser. 443, 012011 (2013).
  * Ahuja (2013) B. L. Ahuja, AIP Conf Proc 1512, 26 (2013).
  * Timms  _et al._ (1993) D. N. Timms, E. Zukowski, M. J. Cooper, D. Laundy, S. P. Collins, F. Itoh, H. SakuraI, T. Iwazumi, H. Kawata, M. Ito, N. Sakai, and Y. Tanaka, J. Phys. Soc. Japan 62, 1716 (1993).
  * Cooper  _et al._ (1992) M. J. Cooper, E. Zukowski, S. P. Collins, D. N. Timms, F. Itoh, and H. Sakurai, J. Condens. Matter Phys. 4, L399 (1992).
  * Cooper  _et al._ (1991) M. J. Cooper, S. P. Collins, S. W. Lovesey, D. Laundy, and D. N. Timms, Phys. Scr. T35, 103 (1991).
  * Duffy  _et al._ (2010) J. A. Duffy, J. W. Taylor, S. B. Dugdale, C. Shenton-Taylor, M. W. Butchers, S. R. Giblin, M. J. Cooper, Y. Sakurai, and M. Itou, Phys. Rev. B 81, 134424 (2010).
  * Itou  _et al._ (2013) M. Itou, A. Koizumi, and Y. Sakurai, Appl. Phys. Lett. 102, 082403 (2013).
  * Zukowski  _et al._ (1993) E. Zukowski, S. P. Collins, M. J. Cooper, D. N. Timms, F. Itoh, H. Sakurai, H. Kawata, Y. Tanaka, and A. Malinowski, J. Condens. Matter Phys. 5, 4077 (1993).
  * Duffy  _et al._ (1998) J. A. Duffy, J. E. McCarthy, S. B. Dugdale, V. Honkimäki, M. J. Cooper, M. A. Alam, T. Jarlborg, and S. B. Palmer, J. Condens. Matter Phys. 10, 10391 (1998).
  * Duffy  _et al._ (2000) J. A. Duffy, S. B. Dugdale, J. E. McCarthy, M. A. Alam, M. J. Cooper, S. B. Palmer, and T. Jarlborg, Phys. Rev. B 61, 14331 (2000).
  * Banfield  _et al._ (2005) Z. F. Banfield, J. A. Duffy, J. W. Taylor, C. A. Steer, A. Bebb, M. J. Cooper, L. Blaauw, C. Shenton-Taylor, and R. Ruiz-Bustos, J. Condens. Matter Phys. 17, 5533 (2005).
  * Shenton-Taylor  _et al._ (2007) C. Shenton-Taylor, J. A. Duffy, J. W. Taylor, C. A. Steer, D. N. Timms, M. J. Cooper, and L. V. Blaauw, J. Condens. Matter Phys. 19, 186208 (2007).
  * Mijnarends  _et al._ (2007) P. E. Mijnarends, S. Kaprzyk, B. Barbiellini, Y. Li, J. F. Mitchell, P. A. Montano, and A. Bansil, Phys. Rev. B 75, 014428 (2007).
  * Mizoroki  _et al._ (2011) T. Mizoroki, M. Itou, Y. Taguchi, T. Iwazumi, and Y. Sakurai, Appl. Phys. Lett. 98, 052107 (2011).
  * Watanabe and Yanase (2018) H. Watanabe and Y. Yanase, Phys. Rev. B 98, 245129 (2018).
  * Bhowal and Spaldin (2021) S. Bhowal and N. A. Spaldin, Phys. Rev. Research 3, 033185 (2021).
  * Ederer and Spaldin (2007) C. Ederer and N. A. Spaldin, Phys. Rev. B 76, 214404 (2007).
  * Spaldin  _et al._ (2008) N. A. Spaldin, M. Fiebig, and M. Mostovoy, J. Condens. Matter Phys. 20, 434203 (2008).
  * Spaldin  _et al._ (2013) N. A. Spaldin, M. Fechner, E. Bousquet, A. Balatsky, and L. Nordström, Phys. Rev. B 88, 094429 (2013).
  * Ernsting  _et al._ (2014) D. Ernsting, D. Billington, T. D. Haynes, T. E. Millichamp, J. W. Taylor, J. A. Duffy, S. R. Giblin, J. K. Dewhurst, and S. B. Dugdale, J. Condens. Matter Phys. 26, 495501 (2014).
  * (27) “The Elk Code,” http://elk.sourceforge.net/.
  * Thöle  _et al._ (2020) F. Thöle, A. Keliri, and N. A. Spaldin, J. Appl. Phys. 127, 213905 (2020).
  * Spaldin (2021) N. A. Spaldin, J. Exp. Theor. Phys. 159, 594 (2021).
  * Nelmes and Kuhs (1985) R. Nelmes and W. Kuhs, Solid State Commun. 54, 721 (1985).
  * (31) The real and momentum space representations are obtained following Ref. Watanabe and Yanase (2018) and using the compatibility relations Inui  _et al._ (1990); Perez-Mato  _et al._ (2015).
  * Inui  _et al._ (1990) T. Inui, Y. Tanabe, and Y. Onodera, Springer, Berlin 78 (1990).
  * Perez-Mato  _et al._ (2015) J. Perez-Mato, S. Gallego, E. Tasci, L. Elcoro, G. de la Flor, and M. Aroyo, Annu. Rev. Mater. Res. 45, 217 (2015).
  * (34) We used the LDA+SOC formalism and a $6\times 6\times 5$ k-point sampling of the Brillouin Zone (with no reduction of $k$-point set, i.e., reducek =0) to achieve the self-consistency. The MCP is, then, calculated following the implementations, discussed in Ref. Ernsting  _et al._ (2014). In order to get access to the MCPs along $p_{y}$ ($p_{x}$) for the $S_{x}$ ($S_{y}$) components we rotate the structure about the Cartesian $\hat{y}$ ($\hat{x}$) directions by $90^{\circ}$ ($270^{\circ}$). The computed MCP is then scaled to the factor that normalizes the valence contribution of the total Compton profile to the total number of valence electrons in the calculation, which is 44 for PTO.
  * Dixon  _et al._ (1998) M. A. G. Dixon, J. A. Duffy, S. Gardelis, J. E. McCarthy, M. J. Cooper, S. B. Dugdale, T. Jarlborg, and D. N. Timms, J. Condens. Matter Phys. 10, 2759 (1998).
  * Slater and Koster (1954) J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).
  * Harrison (1989) W. A. Harrison,  _Electronic Structure and the Properties of Solids: The Physics of the Chemical Bond_ (Dover, New York, 1989).
  * Andersen and Saha-Dasgupta (2000) O. K. Andersen and T. Saha-Dasgupta, Phys. Rev. B 62, R16219 (2000).
  * Sugano (1970) S. Sugano,  _Multiplets of Transition-Metal Ions in Crystals_ (Academic Press, London, 1970).
  * Chatterji  _et al._ (2018) T. Chatterji, S. Rols, and U. D. Wdowik, Front. Phys. 14, 23601 (2018).
  * Kagdada  _et al._ (2018) H. L. Kagdada, P. K. Jha, P. Śpiewak, and K. J. Kurzydłowski, Phys. Rev. B 97, 134105 (2018).
  * (42) S. Bhowal, D. O’Neill, N. A. Spaldin, U. Staub, J. Duffy, and S. P. Collins, “Anti-symmetric Compton scattering in LiNiPO4: Towards a direct probe of the magneto-electric multipole moment,” arXiv:2106.02431(2021).
  * Samara (1971) G. A. Samara, Ferroelectrics 2, 277 (1971).
  * Chattopadhyay  _et al._ (1987) T. Chattopadhyay, J. X. Boucherle, and H. G. vonSchnering, J. phys., C, Solid state phys. 20, 1431 (1987).
  * Wdowik  _et al._ (2014) U. D. Wdowik, K. Parlinski, S. Rols, and T. Chatterji, Phys. Rev. B 89, 224306 (2014).
  * Hayami  _et al._ (2018) S. Hayami, M. Yatsushiro, Y. Yanagi, and H. Kusunose, Phys. Rev. B 98, 165110 (2018).
  * Watanabe and Yanase (2020) H. Watanabe and Y. Yanase, Phys. Rev. Research 2, 043081 (2020).
  * Onimaru  _et al._ (2011) T. Onimaru, K. T. Matsumoto, Y. F. Inoue, K. Umeo, T. Sakakibara, Y. Karaki, M. Kubota, and T. Takabatake, Phys. Rev. Lett. 106, 177001 (2011).
  * Fu (2015) L. Fu, Phys. Rev. Lett. 115, 026401 (2015).
  * Yatsushiro and Hayami (2020) M. Yatsushiro and S. Hayami, Phys. Rev. B 102, 195147 (2020).
  * Ishizuka and Yanase (2018) J. Ishizuka and Y. Yanase, Phys. Rev. B 98, 224510 (2018).
  * Sumita and Yanase (2020) S. Sumita and Y. Yanase, Phys. Rev. Research 2, 033225 (2020).
  * Arras  _et al._ (2019) R. Arras, J. Gosteau, H. J. Zhao, C. Paillard, Y. Yang, and L. Bellaiche, Phys. Rev. B 100, 174415 (2019).
  * Das  _et al._ (2019) S. Das, Y. L. Tang, Z. Hong, M. A. P. Gonçalves, M. R. McCarter, C. Klewe, K. X. Nguyen, F. Gómez-Ortiz, P. Shafer, E. Arenholz, V. A. Stoica, S. L. Hsu, B. Wang, C. Ophus, J. F. Liu, C. T. Nelson, S. Saremi, B. Prasad, A. B. Mei, D. G. Schlom, J. Íñiguez, P. García-Fernández, D. A. Muller, L. Q. Chen, J. Junquera, L. W. Martin, and R. Ramesh, Nature 568, 368 (2019).
  * Göbel  _et al._ (2019) B. Göbel, A. Mook, J. Henk, and I. Mertig, Phys. Rev. B 99, 060406 (2019).