# An Infinite Family of Integrable Sigma Models Using Auxiliary Fields

Christian Ferko<EMAIL_ADDRESS>Center for Quantum Mathematics and
Physics (QMAP), Department of Physics & Astronomy, University of California,
Davis, CA 95616, USA.  Liam Smith<EMAIL_ADDRESS>School of Mathematics
and Physics, University of Queensland, St Lucia, Brisbane, Queensland 4072,
Australia.

###### Abstract

We introduce a class of sigma models in two spacetime dimensions which are
parameterized by an interaction function of one real variable. In addition to
the physical group-valued field $g$, these models include an auxiliary vector
field $v_{\alpha}$ which mediates interactions in a prescribed way. We prove
that every model in this family is weakly integrable, in the sense that the
classical equations of motion are equivalent to flatness of a Lax connection
for any value of a spectral parameter. We also show that these models are
strongly integrable, in the sense that the Poisson bracket of the Lax
connection takes the Maillet form, which guarantees the existence of an
infinite set of conserved charges in involution. This class of theories
includes the principal chiral model (PCM) and all deformations of the PCM by
functions of the energy-momentum tensor, such as ${T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$ and root-${T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$.

##  I Introduction

Integrable models in two dimensions offer rare examples of interacting field
theories which can nonetheless be analyzed exactly, both at the classical and
quantum levels. The dynamics of such models are tightly constrained by the
presence of hidden symmetries, which can sometimes be leveraged to “solve” the
theories. It is therefore useful to identify many examples of such integrable
models, which has been a topic of great recent interest. These integrable $2d$
theories find applications in many areas of physics, such as in studies of
classical and quantum strings at finite tension [1, 2, 3] and in condensed
matter descriptions of spin chains (see e.g. [4, 5, 6]).

A fruitful method for generating new integrable theories is to deform existing
ones in a way which preserves integrability. One popular starting point for
this enterprise is the principal chiral model (PCM) [7, 8], which is an
integrable $2d$ model that shares certain properties with $4d$ Yang-Mills
theory [9, 10, 11, 12]. The PCM admits many integrable deformations, some of
which are reviewed in [13, 14, 15, 16, 17]; for instance, one can add a Wess-
Zumino term [18, 19, 20, 21], perform a Yang-Baxter deformation [22, 23],
implement a $\lambda$-deformation [24], and so on.

In this letter, we introduce a new, infinite class of integrable sigma models
by deforming the PCM via the activation of specific interactions that preserve
integrability. The key technical ingredient in our construction is the
coupling of the physical field in the PCM to an auxiliary field with algebraic
equations of motion. There are many examples in field theory where the
introduction of auxiliary fields can impose a desirable constraint, such as in
the Pasti-Sorokin-Tonin formalism [25, 26, 27] which uses an auxiliary vector
to implement a chirality condition on a tensor field in $6d$, or the Ivanov-
Zupnik auxiliary field formalism [28, 29] which imposes electric-magnetic
duality invariance in theories of $4d$ electrodynamics. In our case, the
inclusion of an auxiliary field imposes that the model remain integrable in a
particularly simple way, while nonetheless allowing for quite general
interactions.

An important step in the analysis of a classically integrable field theory is
the construction of an infinite set of Poisson-commuting conserved charges.
For the PCM and many of its deformations, this procedure can be performed in
two steps. The first step is to establish what is called “weak integrability,”
which means that the equations of motion for the model can be encoded in the
flatness of a one-form called the Lax connection. The second step is to show
“strong integrability” by computing a particular Poisson bracket involving the
spatial component of the Lax connection and demonstrating that it takes a
special form considered by Maillet [30, 31]. This Maillet form of the Poisson
bracket is parameterized by an $r$-matrix $r_{12}(z,z^{\prime})$, which in
many cases takes a standard form in terms of a “twist function” $\varphi(z)$.
Integrable deformations of the PCM often modify the form of this twist
function $\varphi(z)$; see [32] for examples.

In this work, we perform the two steps mentioned above for our new family of
models. We find that, unlike some other integrable deformations of the PCM,
the twist function is not modified for theories within our class, although the
dependence of the Lax connection on the fields is changed. This behavior is
similar to that of homogeneous Yang-Baxter deformations [33, 34, 35], which
can be interpreted via twisted boundary conditions (it would be interesting to
explore whether a similar interpretation applies to our models). As we will
see, the proofs of weak and strong integrability for our models are nearly as
simple as those for the PCM, even though a generic model in this family
exhibits highly non-linear interactions.

The structure of this letter is as follows. In Section II, we introduce the
family of sigma models of interest in this work and derive the equations of
motion for these theories. In Section III, we establish weak integrability by
proving that the equations of motion are equivalent to flatness of a Lax
connection. Section IV proves strong integrability by computing the Poisson
bracket of the spatial component of this Lax connection with itself, and
showing that it takes the Maillet form. In Section V, we show that our family
of models includes all deformations of the PCM by functions of the stress
tensor, such as ${T\mkern 1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern
1.5mu}$ and root-${T\mkern 1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern
1.5mu}$. Finally, Section VI summarizes and presents directions for future
research.

##  II Auxiliary Field Sigma Models

We now describe the class of models which are of interest in the present work.
The physical degree of freedom is a group-valued field $g:\Sigma\to G$ which
maps a flat two-dimensional spacetime $\Sigma$ with coordinates
$\sigma^{\alpha}=(\tau,\sigma)$ into a Lie group $G$. We write $\mathfrak{g}$
to denote the Lie algebra of $G$. The left-invariant Maurer-Cartan form is
defined by $j=g^{-1}dg$, and the pullback of this form to $\Sigma$ is

$\displaystyle j_{\alpha}=g^{-1}\partial_{\alpha}g\,.$ (1)

It will be convenient to use light-cone coordinates on $\Sigma$,

$\displaystyle\sigma^{\pm}=\frac{1}{2}\left(\tau\pm\sigma\right)\,,$ (2)

where light-cone indices are lowered or raised with the metric
$\eta_{+-}=\eta_{-+}=-2$ or inverse metric $\eta^{+-}=\eta^{-+}=-\frac{1}{2}$.
The field $j_{\alpha}$ satisfies the Maurer-Cartan identity, which is written
in light-cone coordinates as

$\displaystyle\partial_{+}j_{-}-\partial_{-}j_{+}+[j_{+},j_{-}]=0\,.$ (3)

In our notation, the Lagrangian of the PCM is written

$\displaystyle\mathcal{L}_{\text{PCM}}=\frac{1}{2}\eta^{\alpha\beta}\mathrm{tr}(j_{\alpha}j_{\beta})=-\frac{1}{2}\mathrm{tr}(j_{+}j_{-})\,.$
(4)

We now introduce an additional Lie algebra valued field $v_{\alpha}$, which is
not a physical degree of freedom, but merely an auxiliary field with algebraic
equations of motion. One can take traces to build scalars from $v_{\alpha}$,
such as

$\displaystyle\nu=\mathrm{tr}(v_{+}v_{+})\mathrm{tr}(v_{-}v_{-})\,.$ (5)

The family of models which we wish to study are described by Lagrangians of
the form

$\displaystyle\mathcal{L}$
$\displaystyle=\frac{1}{2}\mathrm{tr}(j_{+}j_{-})+\mathrm{tr}(v_{+}v_{-})+\mathrm{tr}(j_{+}v_{-}+j_{-}v_{+})+E(\nu)\,,$
(6)

which are parameterized by an arbitrary interaction function $E(\nu)$ that
depends on the single variable $\nu$. The structure of this family of
Lagrangians is inspired by the Ivanov-Zupnik formulation of theories of
duality-invariant electrodynamics in four dimensions [28, 29], which also
involves an interaction function that depends on a single variable constructed
from auxiliary fields.

Let us first consider the equations of motion associated with the Lagrangian
(6). Varying the auxiliary field $v_{\alpha}$ gives rise to the Euler-Lagrange
equation

$\displaystyle
j_{\pm}=-v_{\pm}-2E^{\prime}v_{\mp}\mathrm{tr}(v_{\pm}v_{\pm})\,,$ (7)

where we write $E^{\prime}=\frac{dE}{d\nu}$ for the derivative of $E$.

On the other hand, varying the group-valued field $g$ as $\delta g=\epsilon
g$, under which the Maurer-Cartan form varies as $\delta
j_{\pm}=\partial_{\pm}\epsilon+[j_{\pm},\epsilon]$, gives the $g$-field
equation of motion

$\displaystyle\partial_{+}j_{-}+\partial_{-}j_{+}=2\left([v_{-},j_{+}]+[v_{+},j_{-}]-\partial_{+}v_{-}-\partial_{-}v_{+}\right)\,.$
(8)

The interpretations of the two Euler-Lagrange equations (7) and (8) are rather
different. While (8) is a true dynamical condition, which implies that the
physical field $g$ is on-shell, equation (7) is merely a constraint imposed by
the auxiliary field. We introduce the symbol
$\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}$
to denote equality between two quantities which holds when equation (7) is
satisfied. For instance, (7) implies

$\displaystyle[v_{-},j_{+}]+[v_{+},j_{-}]\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}0\,,$
(9)

so the $g$-field equation of motion (8) can be written as

$\displaystyle\partial_{+}\left(j_{-}+2v-\right)+\partial_{-}\left(j_{+}+2v_{+}\right)\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}0\,,$
(10)

when the auxiliary field equation of motion (7) is satisfied. It is convenient
to define the new quantity

$\displaystyle\mathfrak{J}_{\pm}=-\left(j_{\pm}+2v_{\pm}\right)\,,$ (11)

which allows us to express the Euler-Lagrange equation for $g$ as a
conservation equation
$\partial_{\alpha}\mathfrak{J}^{\alpha}\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}0$.
Here $j_{\alpha}$ is flat but not conserved, while $\mathfrak{J}_{\alpha}$ is
conserved but not flat.

Note that, when the interaction function $E(\nu)$ vanishes, the auxiliary
field equation of motion reduces to $j_{\pm}=-v_{\pm}$. In this case, the
Lagrangian (6) reduces to the PCM Lagrangian (4), and the quantity
$\mathfrak{J}_{\alpha}$ simply becomes $j_{\alpha}$, which is indeed conserved
in the PCM.

More generally, when the auxiliary field $v_{\alpha}$ has been eliminated
using its equation of motion, the quantity $\mathfrak{J}_{\alpha}$ becomes the
Noether current associated with left-multiplication of the group-valued field
$g$ by an arbitrary element of $G$. In the theory without auxiliary fields,
conservation of this Noether current is equivalent to the equation of motion
for the model.

##  III Weak Integrability

Ordinarily, we say that a model is weakly integrable if its equations of
motion admit a zero-curvature representation, i.e. if they are equivalent to
the flatness condition

$\displaystyle d\mathfrak{L}=0\,,$ (12)

for a Lax connection $\mathfrak{L}(z)$, at any value of a spectral parameter
$z\in\mathbb{C}$. The curvature of the Lax connection, $d\mathfrak{L}$, can be
expressed in light-cone coordinates as

$\displaystyle
d\mathfrak{L}=\partial_{+}\mathfrak{L}_{-}-\partial_{-}\mathfrak{L}_{+}+[\mathfrak{L}_{+},\mathfrak{L}_{-}]\,.$
(13)

Let us first define what we mean by integrability of an auxiliary field model,
for which we have both the $v_{\alpha}$-field equation of motion (7) and the
$g$-field equation of motion (8). Given a Lagrangian
$\mathcal{L}(j_{\alpha},v_{\alpha})$, one could eliminate the auxiliary field
$v_{\alpha}$ using its equation of motion to write

$\displaystyle\mathcal{L}(j_{\alpha},v_{\alpha})\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}\mathcal{L}(j_{\alpha})\,.$
(14)

Likewise, one could eliminate the auxiliary fields in an expression for a Lax
connection $\mathfrak{L}(j_{\alpha},v_{\alpha})$ to obtain

$\displaystyle\mathfrak{L}(j_{\alpha},v_{\alpha})\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}\mathfrak{L}(j_{\alpha})\,.$
(15)

If the flatness of the resulting Lax connection $\mathfrak{L}(j_{\alpha})$ for
any $z\in\mathbb{C}$ is equivalent to the $g$-field equation of motion arising
from $\mathcal{L}(j_{\alpha})$, we say that the original auxiliary field model
$\mathcal{L}(j_{\alpha},v_{\alpha})$ is weakly integrable with Lax connection
$\mathfrak{L}(j_{\alpha},v_{\alpha})$. An equivalent, but more succinct,
version of this definition is as follows. If we define

$\displaystyle\mathcal{E}_{g}=\partial_{\alpha}\left(\frac{\partial\mathcal{L}}{\partial(\partial_{\alpha}g)}\right)-\frac{\partial\mathcal{L}}{\partial
g}\,,$ (16)

so that $\mathcal{E}_{g}=0$ is the $g$-field equation of motion, then weak
integrability of an auxiliary field model is the statement

$\displaystyle\left(d\mathfrak{L}\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}0\,,\,\forall
z\in\mathbb{C}\,\right)\,\iff\,\left(\mathcal{E}_{g}\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}0\right)\,.$
(17)

We will now demonstrate that every Lagrangian of the form (6) satisfies this
weak integrability condition (17), where the Lax connection is

$\displaystyle\mathfrak{L}_{\pm}=\frac{j_{\pm}\pm
z\mathfrak{J}_{\pm}}{1-z^{2}}\,.$ (18)

First note that the auxiliary field equation (7) implies

$\displaystyle[\mathfrak{J}_{+},\mathfrak{J}_{-}]\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}[j_{+},j_{-}]\,,\qquad[\mathfrak{J}_{+},j_{-}]\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}[j_{+},\mathfrak{J}_{-}]\,,$
(19)

which agree with relations that were found in the study of the root-${T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$-deformed PCM without
auxiliary fields [36]. These formulas allow us to simplify the commutator

$\displaystyle[\mathfrak{L}_{+},\mathfrak{L}_{-}]$
$\displaystyle=\frac{[j_{+},j_{-}]-z\left([j_{+},\mathfrak{J}_{-}]-[\mathfrak{J}_{+},j_{-}]\right)-z^{2}[\mathfrak{J}_{+},\mathfrak{J}_{-}]}{\left(1-z^{2}\right)^{2}}$
$\displaystyle\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}\,\frac{[j_{+},j_{-}]}{1-z^{2}}\,.$
(20)

The curvature (13) of the Lax connection (18) is then

$\displaystyle
d\mathfrak{L}\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}\frac{\partial_{+}j_{-}-\partial_{-}j_{+}+[j_{+},j_{-}]-z\left(\partial_{+}\mathfrak{J}_{-}+\partial_{-}\mathfrak{J}_{+}\right)}{1-z^{2}}\,.$
(21)

The first three terms in the numerator of (21) vanish due to the Maurer-Cartan
identity (3), and the final two terms vanish if and only if
$\partial_{\alpha}\mathfrak{J}^{\alpha}\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}0$,
which is the equation of motion for $g$. We conclude that every model in this
family obeys the weak integrability condition (17).

##  IV Strong Integrability

We say that a field theory is strongly integrable if the theory possesses an
infinite set of conserved charges which are in involution, or mutually
Poisson-commuting. Demonstrating weak integrability, i.e. that the equations
of motion for a model admit a zero-curvature representation, is a useful step
in this direction, since one can then construct an infinite set of conserved
charges using the monodromy matrix. However, the Poisson-commutativity of
these charges is not guaranteed unless the Poisson bracket of the Lax
connection,

$\displaystyle\left\\{\mathfrak{L}_{\sigma,1}(\sigma,z),\mathfrak{L}_{\sigma,2}(\sigma^{\prime},z^{\prime})\right\\}\,,$
(22)

takes a special form. Here we have introduced, for any Lie algebra valued
quantity $X$, the subscript notation

$\displaystyle X_{1}=X\otimes 1\,,\qquad X_{2}=1\otimes X\,,$ (23)

which tensors $X$ with the identity on either side. We do not carefully
distinguish between $\mathfrak{g}$ and $U(\mathfrak{g})$, the universal
enveloping algebra of $\mathfrak{g}$, but strictly speaking $1\in
U(\mathfrak{g})$ and $X_{1},X_{2}\in U(\mathfrak{g})\otimes U(\mathfrak{g})$.
If (22) takes a Sklyanin [37] or Maillet [30, 31] form, then one can appeal to
general theorems to conclude that the theory possesses an infinite set of
conserved charges in involution.

We will now prove that, for every model in the family (6), the Poisson bracket
(22) takes the non-ultralocal Maillet form, which establishes strong
integrability. Remarkably, this argument proceeds almost exactly as in the
analogous proof for the PCM [31], except replacing every instance of
$j_{\tau}$ with $\mathfrak{J}_{\tau}$; see [38] for a review.

To study the canonical structure of models with group-valued fields [39], it
is convenient to introduce local coordinates $\phi^{\mu}$ on the Lie group
$G$, not to be confused with $\varphi(z)$, so that
$g=g(\phi^{\mu}(\sigma^{\alpha}))$. We use early Greek letters like $\alpha$,
$\beta$ for indices on $\Sigma$ and middle Greek letters like $\mu$, $\nu$ for
coordinates on $G$. We also introduce capital early Latin letters (e.g. $A$,
$B$) which label the generators $T_{A}$ of the Lie algebra $\mathfrak{g}$. For
any Lie algebra valued quantity $X$, one can expand in generators as
$X=X^{A}T_{A}$. Furthermore, we may use the pull-back map
$\partial_{\alpha}\phi^{\mu}$ to convert between indices on $\Sigma$ and
indices on $G$. For instance,

$\displaystyle
j_{\alpha}=j_{\mu}^{A}\frac{\partial\phi^{\mu}}{\partial\sigma^{\alpha}}T_{A}\,.$
(24)

We write $\gamma_{AB}$ for the Killing form on $\mathfrak{g}$, which we assume
to be non-degenerate with inverse $\gamma^{AB}$, and we denote the structure
constants by $\tensor{f}{{}_{A}{}_{B}^{C}}$. In our conventions, these two
objects are defined by the relations

$\displaystyle\gamma_{AB}=\text{Tr}\left[{T_{A}T_{B}}\right]\,,\qquad[T_{A},T_{B}]=\tensor{f}{{}_{A}{}_{B}^{C}}T_{C}\,.$
(25)

In terms of these quantities, the canonical momentum $\pi_{\mu}$ which is
conjugate to the coordinates $\phi^{\mu}$ on $G$ is

$\displaystyle\pi_{\mu}=\frac{\partial\mathcal{L}}{\partial(\partial_{\tau}\phi^{\mu})}=j_{\mu}^{A}\left(j_{\tau}^{B}+2v_{\tau}^{B}\right)\gamma_{AB}\,.$
(26)

This means that the quantity
$\mathfrak{J}_{\tau}^{A}=-(j_{\tau}^{B}+2v_{\tau}^{A})$ is simply related to
the canonical momentum as

$\displaystyle\mathfrak{J}_{\tau}^{A}=-\gamma^{AB}\pi_{\mu}j^{\mu}_{B}\,,$
(27)

where we have defined the inverse field $j^{\mu}_{B}$ which satisfies
$j^{A}_{\mu}j^{\mu}_{B}=\tensor{\delta}{{}^{A}_{B}}$. The fundamental Poisson
brackets are then

$\displaystyle\left\\{\pi_{\mu}(\sigma),\phi^{\nu}(\sigma^{\prime})\right\\}=\tensor{\delta}{{}_{\mu}^{\nu}}\delta(\sigma-\sigma^{\prime})\,,$
(28)

which hold at equal times $\tau$. Because no time derivatives of the auxiliary
field $v_{\alpha}$ appear in the Lagrangian, the momentum
$\mathfrak{p}^{\alpha}$ which is conjugate to $v_{\alpha}$ vanishes:

$\displaystyle\mathfrak{p}^{\alpha}=\frac{\partial\mathcal{L}}{\partial(\partial_{\tau}v_{\alpha})}=0\,.$
(29)

In the canonical formulation, our model (6) is therefore a constrained
Hamiltonian system, for which (29) is a primary constraint. This structure is
reminiscent of the Hamiltonian formulation of Maxwell theory, where the
temporal gauge field $A^{0}$ also has vanishing conjugate momentum. A
systematic treatment of the Hamiltonian structure of our models would proceed
using the method of Dirac [40], which implements the primary constraint in the
Hamiltonian via a Lagrange multiplier, and then imposes that the primary
constraint is preserved under time evolution. However, we will not undertake a
detailed discussion of this procedure here, since it is not necessary for our
present purposes. Because we are only interested in computing Poisson brackets
involving $j_{\alpha}$ and $\mathfrak{J}_{\alpha}$, where the latter is
related to the canonical momentum by (27), it suffices to compute brackets
using the fundamental relations (28) and ignore the fact that
$\mathfrak{J}_{\alpha}$ is itself composed from a constrained auxiliary field.
In fact, the auxiliary field equation of motion will play no role whatsoever
in the following calculations.

Using the fundamental brackets (28), one can show

$\displaystyle\\{\mathfrak{J}_{\tau}^{A}(\sigma),\mathfrak{J}_{\tau}^{B}(\sigma^{\prime})\\}$
$\displaystyle=\tensor{f}{{}^{A}{}^{B}_{C}}\mathfrak{J}_{\tau}^{C}\delta(\sigma-\sigma^{\prime})\,,$
$\displaystyle\left\\{\mathfrak{J}_{\tau}^{A}(\sigma),j_{\sigma}^{B}(\sigma^{\prime})\right\\}$
$\displaystyle=\tensor{f}{{}^{A}{}^{B}_{C}}j_{\sigma}^{C}(\sigma)\delta(\sigma-\sigma^{\prime})-\gamma^{AB}\delta^{\prime}(\sigma-\sigma^{\prime})\,,$
$\displaystyle\left\\{j_{\sigma}^{A}(\sigma),j_{\sigma}^{B}(\sigma^{\prime})\right\\}$
$\displaystyle=0\,,$ (30)

where
$\delta^{\prime}(\sigma-\sigma^{\prime})=\partial^{(\sigma)}\delta(\sigma-\sigma^{\prime})$.
To obtain these results, we have used the Maurer-Cartan identity (3) in the
form

$\displaystyle\partial_{\mu}j_{\nu}^{A}-\partial_{\nu}j_{\mu}^{A}+j_{\mu}^{B}j_{\nu}^{C}\tensor{f}{{}_{B}{}_{C}^{A}}=0\,,$
(31)

along with $\delta$-function identities such as

$\displaystyle f(y)\partial^{(x)}\delta(x-y)$
$\displaystyle=f(x)\partial^{(x)}\delta(x-y)+f^{\prime}(x)\delta(x-y)\,.$ (32)

One can then contract each of the relations (IV) with the tensor product
$T_{A}\otimes T_{B}$ to find

$\displaystyle\\{\mathfrak{J}_{\tau,1}(\sigma),\mathfrak{J}_{\tau,2}(\sigma^{\prime})\\}$
$\displaystyle=[\mathfrak{J}_{\tau,2},C_{12}]\delta(\sigma-\sigma^{\prime})\,,$
$\displaystyle\left\\{\mathfrak{J}_{\tau,1}(\sigma),j_{\sigma,2}(\sigma^{\prime})\right\\}$
$\displaystyle=[j_{\sigma,2},C_{12}]\delta(\sigma-\sigma^{\prime})-C_{12}\delta^{\prime}(\sigma-\sigma^{\prime})\,,$
$\displaystyle\left\\{j_{\sigma,1}(\sigma),j_{\sigma,2}(\sigma^{\prime})\right\\}$
$\displaystyle=0\,,$ (33)

where we have defined the Casimir $C_{12}=\gamma^{AB}T_{A}\otimes T_{B}$.

From (IV) one can show that

$\displaystyle\left\\{\mathfrak{L}_{\sigma,1}(\sigma,z),\mathfrak{L}_{\sigma,2}(\sigma^{\prime},z^{\prime})\right\\}$
$\displaystyle=\frac{1}{(1-z^{2})(1-z^{\prime
2})}\Big{(}\left[zz^{\prime}\mathfrak{J}_{\tau,2}+(z+z^{\prime})j_{\sigma,2},C_{12}\right]\delta(\sigma-\sigma^{\prime})-(z+z^{\prime})C_{12}\delta^{\prime}(\sigma-\sigma^{\prime})\Big{)}$
$\displaystyle=\left[r_{12}(z,z^{\prime}),\mathfrak{L}_{\sigma,1}(\sigma,z)\right]\delta(\sigma-\sigma^{\prime})-\left[r_{21}(z^{\prime},z),\mathfrak{L}_{\sigma,2}(\sigma,z^{\prime})\right]\delta(\sigma-\sigma^{\prime})-s_{12}(z,z^{\prime})\delta^{\prime}(\sigma-\sigma^{\prime})\,,$
(34)

where $s_{12}(z,z^{\prime})=r_{12}(z,z^{\prime})+r_{21}(z^{\prime},z)$, the
$r$-matrix is

$\displaystyle r_{12}(z,z^{\prime})$
$\displaystyle=\frac{C_{12}}{z-z^{\prime}}\varphi^{-1}(z^{\prime})\,,$ (35)

which solves the classical Yang-Baxter equation, and

$\displaystyle\varphi(z)$ $\displaystyle=\frac{z^{2}-1}{z^{2}}$ (36)

is the same twist function as for the undeformed PCM. The structure in the
second line of (IV) is precisely the desired non-ultralocal Maillet form of
the Poisson bracket.

Perhaps surprisingly, the result (IV) takes exactly the same form as in the
PCM. This calculation has been completely blind to the form of the interaction
function $E(\nu)$, and the auxiliary field has entered only through the
combination $\mathfrak{J}_{\tau}$ of equation (27). We stress that, although
the auxiliary field equation of motion was needed in the proof of weak
integrability, this equation is not used anywhere in the computation of the
bracket (IV). We therefore conclude that every theory within the class (6) is
strongly integrable, in the sense that it possesses an infinite collection of
Poisson-commuting conserved charges.

##  V Stress Tensor Flows

The family of sigma models (6) is classically equivalent to the set of all
deformations of the PCM by functions of the stress tensor, which we define by

$\displaystyle T_{\alpha\beta}=-\frac{2}{\sqrt{-g}}\frac{\delta S}{\delta
g^{\alpha\beta}}\,.$ (37)

In components, $\tensor{T}{{}^{\alpha}_{\beta}}$ is a $2\times 2$ matrix,
which means that there are two functionally independent Lorentz scalars that
can be constructed from it. One is the trace,

$\displaystyle\tensor{T}{{}^{\alpha}_{\alpha}}\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}2\left(E+2E^{\prime}\nu\right)\,,$
(38)

where we have simplified using the auxiliary field equation of motion (7). The
other is $\tensor{T}{{}^{\alpha}_{\beta}}\tensor{T}{{}^{\beta}_{\alpha}}$,
which satisfies

$\displaystyle\tensor{T}{{}^{\alpha}_{\beta}}\tensor{T}{{}^{\beta}_{\alpha}}-\frac{1}{2}\left(\tensor{T}{{}^{\alpha}_{\alpha}}\right)^{2}\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}\frac{1}{2}\nu\left(1-4(E^{\prime})^{2}\nu\right)^{2}\,.$
(39)

Importantly, both of these contractions can be expressed entirely in terms of
$\nu$, $E(\nu)$, and $E^{\prime}(\nu)$, when the auxiliary field equation of
motion is satisfied. This means that any classical flow equation of the form

$\displaystyle\frac{\partial\mathcal{L}}{\partial\lambda}\stackrel{{\scriptstyle\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{{=}}f\left(\tensor{T}{{}^{\alpha}_{\alpha}},\tensor{T}{{}^{\alpha}_{\beta}}\tensor{T}{{}^{\beta}_{\alpha}},\lambda\right)\,,$
(40)

gives rise to an ordinary differential equation for the interaction function
$E(\nu)$, whose solution is another model within the same class of theories.
This structure is similar to that of classical stress tensor deformations of
theories of duality-invariant electrodynamics in $4d$ [41], chiral tensor
theories in $6d$ [42], or chiral bosons in $2d$ [43].

Examples of deformations (40) include the ${T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$ [44, 45] and
root-${T\mkern 1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$ [46,
47, 48] flows, driven by the functions

$\displaystyle\mathcal{O}_{{T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}}$
$\displaystyle=\tensor{T}{{}^{\alpha}_{\beta}}\tensor{T}{{}^{\beta}_{\alpha}}-\left(\tensor{T}{{}^{\alpha}_{\alpha}}\right)^{2}\,,$
(41) $\displaystyle\mathcal{R}$
$\displaystyle=\frac{1}{\sqrt{2}}\sqrt{\tensor{T}{{}^{\alpha}_{\beta}}\tensor{T}{{}^{\beta}_{\alpha}}-\frac{1}{2}\left(\tensor{T}{{}^{\alpha}_{\alpha}}\right)^{2}}\,.$
(42)

The solution for the interaction function $E(\nu)$ deformed by the
root-${T\mkern 1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$ flow
driven by $\mathcal{R}$, subject to the initial condition
$E(\lambda=0,\nu)=0$, is

$\displaystyle
E(\lambda,\nu)=\tanh\left(\frac{\lambda}{2}\right)\sqrt{\nu}\,,$ (43)

which is the same as the interaction function for the Modified Maxwell theory
[49] in the $4d$ Ivanov-Zupnik formulation. The solution to the ${T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$ flow equation driven
by (41) with the same initial condition can be written in terms of a
hypergeometric function [42].

It is well-known that applying the ${T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$ deformation to a
classically integrable theory produces another classically integrable theory
[50], and the Lax connection for the ${T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$-deformed PCM was
obtained in [51]. Likewise, the Lax connection for the PCM deformed by both
${T\mkern 1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$ and
root-${T\mkern 1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$
appeared in [36]. Both of these results can be viewed as special cases of the
general Lax connection (18), which is appropriate for a deformation of the PCM
by an arbitrary function of the energy-momentum tensor.

##  VI Conclusions and Outlook

In this work, we have introduced a new auxiliary field technique for proving
integrability of $2d$ sigma models. We used this machinery to establish both
weak and strong integrability for an infinite class of interacting PCM-like
theories, which includes all deformations of the PCM by functions of the
energy-momentum tensor.

We emphasize that, by virtue of the auxiliary fields in our model, the proofs
of weak and strong integrability for this entire infinite class are hardly
more difficult than those for the PCM. It would be very interesting to see
whether similar techniques could be used to streamline proofs of integrability
for other theories, or to construct larger families of integrable models. For
instance, it seems likely that stress tensor deformations of symmetric and
semi-symmetric space sigma models, with or without Wess-Zumino term, may also
be addressed within our framework; weak integrability of ${T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$ and root-${T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$ deformations of these
models was established in [36].

Our discussion has been entirely classical. It would be interesting to see
whether our auxiliary field framework aids in the study of integrability for
models at the quantum level. At least in the case of the ${T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$-deformed PCM, it is
believed that the theory is quantum-mechanically well-defined due to the
point-splitting definition of the ${T\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$ operator [44]. Perhaps
other models within this class also admit quantum definitions, and auxiliary
fields may be useful for studying quantum integrability.

A third avenue for future research is to investigate whether there is a
connection between these deformed integrable sigma models and $4d$ Chern-
Simons theory [52, 53, 54]; see [55] for a recent review. The conventional
${T\mkern 1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu}$ deformation
admits an interpretation via $4d$ Chern-Simons, at least to leading order in
the deformation parameter [56], and it might be that a similar narrative
applies to the general family of models constructed here.

###### Acknowledgements.

We thank Gabriele Tartaglino-Mazzucchelli and Alessandro Sfondrini for helpful
discussions and comments on a draft of this letter. We also acknowledge kind
hospitality and financial support at the MATRIX Program “New Deformations of
Quantum Field and Gravity Theories,” where this work was initiated, and thank
the participants of this meeting for productive conversations on related
topics. C. F. is supported by U.S. Department of Energy grant DE-SC0009999 and
funds from the University of California. L. S. is supported by a postgraduate
scholarship at the University of Queensland and was partially supported by the
Australian Research Council (ARC) grants FT180100353 and DP240101409. This
research was supported in part by grant NSF PHY-2309135 to the Kavli Institute
for Theoretical Physics (KITP).

## References

  * Arutyunov and Frolov [2009] G. Arutyunov and S. Frolov, J. Phys. A A42, 254003 (2009), arXiv:0901.4937 [hep-th] .
  * Beisert  _et al._ [2012] N. Beisert, C. Ahn, L. F. Alday, Z. Bajnok, J. M. Drummond, L. Freyhult, N. Gromov, R. A. Janik, V. Kazakov, T. Klose, G. P. Korchemsky, C. Kristjansen, M. Magro, T. McLoughlin, J. A. Minahan, R. I. Nepomechie, A. Rej, R. Roiban, S. Schäfer-Nameki, C. Sieg, M. Staudacher, A. Torrielli, A. A. Tseytlin, P. Vieira, D. Volin, and K. Zoubos, Lett. Math. Phys. 99, 3 (2012), arXiv:1012.3982 [hep-th] .
  * Demulder  _et al._ [2023] S. Demulder, S. Driezen, B. Knighton, G. Oling, A. L. Retore, F. K. Seibold, A. Sfondrini, and Z. Yan, arXiv:2312.12930 [hep-th] (2023).
  * Haldane [1983] F. Haldane, Physics Letters A 93, 464 (1983).
  * Faddeev [1996] L. D. Faddeev, in  _Les Houches School of Physics: Astrophysical Sources of Gravitational Radiation_ (1996) pp. pp. 149–219, arXiv:hep-th/9605187 .
  * Essler  _et al._ [2005] F. H. L. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin,  _The One-Dimensional Hubbard Model_ (Cambridge University Press, 2005).
  * Gross and Wilczek [1973] D. J. Gross and F. Wilczek, Phys. Rev. Lett. 30, 1343 (1973).
  * Politzer [1973] H. D. Politzer, Phys. Rev. Lett. 30, 1346 (1973).
  * Polyakov [1975] A. Polyakov, Physics Letters B 59, 79 (1975).
  * Polyakov and Belavin [1975] A. M. Polyakov and A. A. Belavin, JETP Lett. 22, 245 (1975).
  * Migdal [1975] A. A. Migdal, Sov. Phys. JETP 42, 413 (1975).
  * Balog  _et al._ [1992] J. Balog, S. Naik, F. Niedermayer, and P. Weisz, Phys. Rev. Lett. 69, 873 (1992).
  * Zarembo [2017] K. Zarembo, arXiv:1712.07725 [hep-th] (2017).
  * Orlando  _et al._ [2020] D. Orlando, S. Reffert, J.-i. Sakamoto, Y. Sekiguchi, and K. Yoshida, J. Phys. A 53, 443001 (2020), arXiv:1912.02553 [hep-th] .
  * Seibold [2020] F. K. Seibold,  _Integrable deformations of sigma models and superstrings_ , Ph.D. thesis, ETH, Zurich (2020).
  * Klimčík [2021] C. Klimčík, Rev. Math. Phys. 33, 2130004 (2021), arXiv:2101.05230 [hep-th] .
  * Hoare [2022] B. Hoare, J. Phys. A 55, 093001 (2022), arXiv:2109.14284 [hep-th] .
  * Wess and Zumino [1971] J. Wess and B. Zumino, Physics Letters B 37, 95 (1971).
  * Novikov [1982] S. P. Novikov, Usp. Mat. Nauk 37N5, 3 (1982).
  * Witten [1984] E. Witten, Commun. Math. Phys. 92, 455 (1984).
  * Abdalla  _et al._ [1982] E. Abdalla, M. Gomes, and M. Forger, Nuclear Physics B 210, 181 (1982).
  * Klimcik [2002] C. Klimcik, JHEP 12, 051, arXiv:hep-th/0210095 .
  * Klimcik [2009] C. Klimcik, J. Math. Phys. 50, 043508 (2009), arXiv:0802.3518 [hep-th] .
  * Sfetsos [2014] K. Sfetsos, Nucl. Phys. B 880, 225 (2014), arXiv:1312.4560 [hep-th] .
  * Pasti  _et al._ [1995] P. Pasti, D. P. Sorokin, and M. Tonin, Phys. Rev. D 52, R4277 (1995), arXiv:hep-th/9506109 .
  * Pasti  _et al._ [1997a] P. Pasti, D. P. Sorokin, and M. Tonin, Phys. Rev. D 55, 6292 (1997a), arXiv:hep-th/9611100 .
  * Pasti  _et al._ [1997b] P. Pasti, D. P. Sorokin, and M. Tonin, Phys. Lett. B 398, 41 (1997b), arXiv:hep-th/9701037 .
  * Ivanov and Zupnik [2002] E. A. Ivanov and B. M. Zupnik, in  _Supersymmetries and Quantum Symmetries_ (2002) pp. 235–250, arXiv:hep-th/0202203 [hep-th] .
  * Ivanov and Zupnik [2004] E. A. Ivanov and B. M. Zupnik, Phys. Atom. Nucl. 67, 2188 (2004), arXiv:hep-th/0303192 .
  * Maillet [1986a] J.-M. Maillet, Nuclear Physics B 269, 54 (1986a).
  * Maillet [1986b] J.-M. Maillet, Physics Letters B 167, 401 (1986b).
  * Lacroix [2018] S. Lacroix,  _Integrable models with twist function and affine Gaudin models_ , Ph.D. thesis, Lyon, Ecole Normale Superieure (2018), arXiv:1809.06811 [hep-th] .
  * Kawaguchi  _et al._ [2014] I. Kawaguchi, T. Matsumoto, and K. Yoshida, JHEP 1404, 153, arXiv:1401.4855 [hep-th] .
  * van Tongeren [2015] S. J. van Tongeren, JHEP 06, 048, arXiv:1504.05516 [hep-th] .
  * Borsato  _et al._ [2022] R. Borsato, S. Driezen, and J. L. Miramontes, JHEP 04, 053, arXiv:2112.12025 [hep-th] .
  * Borsato  _et al._ [2023] R. Borsato, C. Ferko, and A. Sfondrini, Phys. Rev. D 107, 086011 (2023), arXiv:2209.14274 [hep-th] .
  * Sklyanin [1980] E. K. Sklyanin, Zap. Nauchn. Semin. 95, 55 (1980).
  * Driezen [2022] S. Driezen, PoS Modave2021, 002 (2022), arXiv:2112.14628 [hep-th] .
  * Bowcock [1989] P. Bowcock, Nuclear Physics B 316, 80 (1989).
  * Dirac [1964] P. A. M. Dirac,  _Lectures on quantum mechanics_ , Belfer Graduate School of Science Monographs Series, Vol. 2 (Belfer Graduate School of Science, New York, 1964).
  * Ferko  _et al._ [2023] C. Ferko, S. M. Kuzenko, L. Smith, and G. Tartaglino-Mazzucchelli, Phys. Rev. D 108, 106021 (2023), arXiv:2309.04253 [hep-th] .
  * Ferko  _et al._ [2024] C. Ferko, S. M. Kuzenko, K. Lechner, D. P. Sorokin, and G. Tartaglino-Mazzucchelli, arXiv:2402.06947 [hep-th] (2024).
  * Ebert  _et al._ [2024] S. Ebert, C. Ferko, C. L. Martin, and G. Tartaglino-Mazzucchelli, arXiv:2403.18242 [hep-th] (2024).
  * Zamolodchikov [2004] A. B. Zamolodchikov, arXiv:hep-th/0401146 [hep-th] (2004).
  * Cavaglià  _et al._ [2016] A. Cavaglià, S. Negro, I. M. Szécsényi, and R. Tateo, JHEP 10, 112, arXiv:1608.05534 [hep-th] .
  * Ferko  _et al._ [2022] C. Ferko, A. Sfondrini, L. Smith, and G. Tartaglino-Mazzucchelli, Phys. Rev. Lett. 129, 201604 (2022), arXiv:2206.10515 [hep-th] .
  * Babaei-Aghbolagh  _et al._ [2022a] H. Babaei-Aghbolagh, K. Babaei Velni, D. Mahdavian Yekta, and H. Mohammadzadeh, Phys. Rev. D 106, 086022 (2022a), arXiv:2206.12677 [hep-th] .
  * Babaei-Aghbolagh  _et al._ [2022b] H. Babaei-Aghbolagh, K. B. Velni, D. M. Yekta, and H. Mohammadzadeh, Phys. Lett. B 829, 137079 (2022b), arXiv:2202.11156 [hep-th] .
  * Bandos  _et al._ [2020] I. Bandos, K. Lechner, D. Sorokin, and P. K. Townsend, Phys. Rev. D 102, 121703 (2020), arXiv:2007.09092 [hep-th] .
  * Smirnov and Zamolodchikov [2017] F. A. Smirnov and A. B. Zamolodchikov, Nucl. Phys. B915, 363 (2017), arXiv:1608.05499 [hep-th] .
  * Chen  _et al._ [2021] B. Chen, J. Hou, and J. Tian, Chin. Phys. C 45, 093112 (2021), arXiv:2102.01470 [hep-th] .
  * Costello and Yamazaki [2019] K. Costello and M. Yamazaki, arXiv:1908.02289 [hep-th] (2019).
  * Nekrasov [1996] N. Nekrasov,  _Four Dimensional Holomorphic Theories_ , Ph.D. thesis, Princeton University (1996).
  * Costello [2013] K. Costello, arXiv:1303.2632 [hep-th] (2013).
  * Lacroix [2022] S. Lacroix, J. Phys. A 55, 083001 (2022), arXiv:2109.14278 [hep-th] .
  * Py [2022] V. Py, JHEP 08, 101, arXiv:2202.08841 [hep-th] .


# Spin polarisation and non-isotropic effective mass in the conduction band of
GdN

W. F. Holmes-Hewett1,3,∗ ,E. X. M. Trewick2,3,, H. J. Trodahl2, R. G.
Buckley1,3 and B. J. Ruck2,3 1Robinson Research Institute, Victoria University
of Wellington, P.O. Box 33436, Petone 5046, New Zealand  2School of Chemical
and Physical Sciences, Victoria University of Wellington, P.O. Box 600,
Wellington 6140, New Zealand  3MacDiarmid Institute for Advanced Materials and
Nanotechnology, P.O. Box 600, Wellington 6140, New Zealand  *Corresponding
author: W. F. Holmes-Hewett<EMAIL_ADDRESS>
###### Abstract

GdN is a ferromagnetic semiconductor which has seen increasing interest in the
preceding decades particularly in the areas of spin- and superconducting-
based electronics. Here we report a detailed computational study and optical
spectroscopy study of the electronic structure of stoichiometric and nitrogen
vacancy doped GdN. Based on our calculations we provide the effective mass
tensor for undoped GdN, and some indicative values for electron doped GdN.
Such a property is valuable as it can directly affect device design, and be
directly measured experimentally to validate the existing computation results.

###### pacs:

71.27.+a, 75.50.Pp

##  I Introduction

The increasing global demand for computational efficiency Andrae (2015) has
pushed the development of conventional CMOS based logic and memory beyond the
imagination of the materials scientists and physicists who pioneered the first
generation of transistor based devices decades ago. With the limits of the
present family of conventional-materials based devices (perpetually) on the
horizon there has been an ongoing effort to design more exotic logic and
memory systems, systems which can surpass the conventional technologies
leading to orders of magnitude increases in operating efficiency and speed IRD
(2018a, b).

Among the most attractive of these technologies are spintronics Hirohata  _et
al._ (2020); Bhatti  _et al._ (2017) and superconducting electronics Soloviev
_et al._ (2017); Likharev and Semenov (1991). Both of these operate at
cryogenic temperatures, temperatures which limit the use of conventional
materials. As such, a search for materials with the functionalities required
for these technologies and the ability to operate at the required temperatures
is an active area of materials science Li and Yang (2016); Shamiul Alam
(2023). In this context the rare-earth nitride (LN - L a lanthanide) series of
ferromagnetic semiconductors Natali  _et al._ (2010) has seen a burst of
interest in recent years. First studied in the 1950s Eick  _et al._ (1956) the
LN were found to be mostly ferromagnetic with Curie temperatures on the order
of 50 K Hulliger (1978). The simple rock-salt structure and limited quenching
of the orbital angular momentum results in rich and varied magnetic behaviours
across the series, which is driven by the 4f occupation ranging from 0 in LaN
to 14 in LuN. The series was initially thought to be half-metallic, similar to
the other rare-earth pnicitdes Hulliger (1979). However, improved experimental
procedures allowing the growth of high quality thin films, along with improved
computational methods, has in recent decades revealed the semiconducting
nature of at least several members of the series Granville  _et al._ (2006);
Aerts  _et al._ (2004); Anton  _et al._ (2016); Holmes-Hewett  _et al._ (2020,
2023); Devese  _et al._ (2022).

Figure 1: Calculated band structure of stoichiometric GdN (left - majority
spin, right - minority spin) along the line $\Gamma$-X. The wavefunctions are
projected onto the atomic orbitals, Gd 5d \- green, Gd 4f \- red N 2p \-
purple.

GdN is the most studied of the rare-earth nitrites, with a Curie temperature
reported between 50 K and 70 K Ludbrook  _et al._ (2009). The half filled 4f
shell yields a strong magnetic moment of 7 $\mu_{B}$ per Gd ion, or 32
$\mu_{B}$ per nm3. GdN is insulating in the stoichiometric state with a
valence band formed from N 2p states and a conduction band minimum formed from
Gd 5d states. The stoichiometric band structure calculated in the present
study is shown in Figure 1. The valence band maximum is at $\Gamma$ and there
is a $\sim$ 0.9 eV optical band gap between the majority spin states at the X
point Trodahl  _et al._ (2007), to which we have tuned our calculation (see
section II). The minority spin gap is significantly larger in the calculation
($\sim$ 2.5 eV) between the minority spin N 2p and Gd 5d states at X, a result
of a large spin-splitting in the conduction ($1.3$ eV) and the valence ($0.2$
eV) bands both resulting from exchange with the Gd 4f states Larson and
Lambrecht (2006). Previous computational studies Larson  _et al._ (2007);
Preston  _et al._ (2010), also tuned to the optical gap in the ferromagnetic
phase, find a significantly smaller minority spin gap of $\sim~{}1.7$ eV.

The onset of optical absorption in the ferromagnetic phase represents the
majority spin gap, as both the conduction band minimum and valence band
maximum are majority spin. There is largely agreement in the existing
experimental studies regarding the 0.9 eV majority spin gap Trodahl  _et al._
(2007); Vilela  _et al._ (2024); Azeem (2016); Yoshitomi  _et al._ (2011).
There is however much less agreement regarding the minority spin gap and
associated spin-splitting, the latter of which is reported variously as 0.16
eV, 0.4 eV and 0.6 eV respectively in Refs Vidyasagar  _et al._ (2012); Azeem
(2016) and Vilela  _et al._ (2024). In the context of majority-minority spin
gaps the paramagnetic gap can be estimated as an average of the two. The
present calculation shown in Figure 1 yields 1.75 eV, somewhat above the
experimental paramagnetic gap of 1.5 eV.

As with many of the other rare earth nitrides GdN can be doped with electrons
though the inclusion of nitrogen vacancies Maity  _et al._ (2018); Holmes-
Hewett  _et al._ (2020); Devese  _et al._ (2022). These have a low activation
energy, dependent on growth conditions Punya  _et al._ (2011), so can be
readily formed in thin films to the order of a few % by either reducing the N2
partial pressure during deposition, or controlling the substrate temperature.
Experimental results on intentionally doped GdN films support the simple
argument of the nitrogen vacancy donors raising the Fermi level into the
(initially) majority spin only Gd 5d bands Trodahl  _et al._ (2017).

Figure 2: Calculated band structure for electron doped GdN1-δ highlighting the
dispersion around the conduction band minimum at X: left - majority spin,
right - minority spin. The wave functions are projected onto the atomic
orbitals, Gd 5d \- green, Gd 4f \- red N 2p \- purple. The Gd 5d states from
the stochometric crystal are shown in the left hand panel in black dashed
lines. The inset shows the Fermi surface highlighted on the first Brillouin
zone.

The rare confluence of a combine semiconducting-ferromagnetic ground state has
seen GdN included in many prototype cryogenic devices. The material is
particularly attractive in the area of superconducting spintronics. The
inclusion of ferromagnetic materials as the weak link in a Josephson junction
structure adds functionality based on the hysteric nature of the magnetisation
Birge and Satchell (2024), however the low resistance of the ferromagnetic-
metals commonly used results in switching speeds orders of magnitude lower
than in conventional superconducting logic, rendering conventional
ferromagnetic-metal based Josephson junctions incompatible Vernik  _et al._
(2013). In this context GdN has been used as the ferromagnetic insulator of
choice in the superconducting spintronics community, with various devices
already demonstrated including superconducting- spin-filters Senapati  _et
al._ (2011); Massarotti  _et al._ (2015); Ahmad  _et al._ (2020), diodes
Sharma  _et al._ (2023), primitive memory devices Cascales  _et al._ (2019);
Pot  _et al._ (2023), and even the recent proposal of a GdN based ferro-
transmon qubit Ahmad  _et al._ (2022).

Given the roughly a decade of intense interest in the application of GdN in
such devices it is surprising that the most recent computational studies
regarding the electronic structure of GdN are themselves over a decade old
Larson  _et al._ (2007); Punya  _et al._ (2011); Mitra and Lambrecht (2008).
In the present manuscript we attempt to fill this void by undertaking a
computational study of GdN, informed by experiment, where we aim to provide
details regarding electron transport in GdN and realistic material parameters
which will aid experimentalists and engineers in device design.

##  II Methods

###  II.1 Computational

Density functional theory (DFT) based calculations were undertaken using
Quantum Espresso Giannozzi  _et al._ (2009); Cococcioni and de Gironcoli
(2005) and rare earth pseudo-potentials developed using the rare-earth
nitrides Topsakal and Wentzcovitch (2014). Self-consistent calculations on the
primitive cell were completed using a $k$-mesh with $10\times 10\times 10$
divisions, while super-cell calculations were on a $4\times 4\times 4$
division k-mesh. The wave function and charge density cut-off energies were 50
Ry and 200 Ry respectively for all calculations. Following our DFT
calculations the output from Quantum Espresso was used to generate maximally
localised Wannier functions using Wannier 90 Mostofi  _et al._ (2014); Marzari
and Vanderbilt (1997); Souza  _et al._ (2001). The resulting Wannier functions
were then used to calculate the electronic properties on denser k meshes, a
$25\times 25\times 25$ grid for calculation of the density of states and Fermi
level) and a $125\times 125\times 125$ grid for the optical conductivity. The
effective mass was determined using conventional DFT based on 400 k along each
high symmetry direction.

The 4$f$ electrons of the LnN series are strongly correlated and thus require
careful treatment beyond the traditional DFT methods Larson and Lambrecht
(2006); Larson  _et al._ (2007). In the basic DFT the 4$f$ states are found at
or near the Fermi level for most of the stoichiometric rare earth nitrides,
when in reality this is not the case. The strongly correlated nature of these
electrons pushes the filled states below and unfilled states above the Fermi
level. This physics can be approximated using the DFT+$U$ method where the
behaviour of the correlated orbitals is determined by an adjustable parameter
$U$. In the present study two $U$ parameters are used, as described first in
reference Larson and Lambrecht (2006). One to account for the strongly
correlated 4$f$ states ($U_{f}$), and a second applied to the 5$d$ states
($U_{d}$) to adjust to the experimental optical band gap in the ferromagnetic
phase. The process is discussed further in references Holmes-Hewett (2021) and
Holmes-Hewett  _et al._ (2023).

###  II.2 film growth and optical spectroscopy

GdN thin films on the order of 100 nm were prepared in a Thermionics ultrahigh
vacuum chamber with a base pressure of $10^{-9}$ mbar. Gd metal was evaporated
by an electron gun in a N2 partial pressure of $1\times 10^{-4}$ mbar at a
rate of $\sim 1$ Å s-1. Varying the N2 growth pressure acts to control the
concentration of nitrogen vacancies (VN) in the films. The films were capped
with $\sim 50$ nm insulating AlN to minimise deterioration when exposed
atmospheric water vapour and oxygen. The films were simultaneously deposited
at ambient temperature on Si and Al2O3 substrates to suit the required
measurements.

Samples for electrical transport measurements were deposited on
$10~{}\times~{}10$ mm2 Al2O3 substrates predeposited with Cr/Au contacts in
the van der Pauw configuration. Electrical measurements were conducted in a
Janis closed-cycle helium cryostat and a Quantum Design Physical Properties
System at temperatures from 300 to 4 K. Reflection and transmission
measurements were performed with a Bruker Vertex 80v Fourier transform
spectrometer in the 0.012-3 eV spectral range using films deposited on both Si
and Al2O3 substrates. The reflectivity was measured relative to a 200 nm Al
film and then corrected for the finite reflectivity of Al using data from
Ehrenreich Ehrenreich  _et al._ (1963). The reflectivity and transmission
spectra were simultaneously modelled with refFIT software Kuzmenko (2005)
using a Kramers-Kronig consistent consistent sum of Drude-Lorentzians to
calculate the energy-dependent dielectric function from which the presented
optical conductivity data was generated.

##  III Results

###  III.1 Density function theory

As is commonly reported, GdN is grown with varying nitrogen vacancy doping to
control the carrier concentration. In the doped crystal the Gd ions in the
lattice retain their 3+ oxidation state, so to maintain charge neutrality the
three electrons formerly on the nitrogen ion must be accounted for. Two
electrons occupy defect states which localise to the vacancy site, while the
third electron finds a higher energy defect state, drawing the Fermi level up
into the extended state Gd 5d conduction band. This one mobile electron per
nitrogen vacancy is responsible for electron transport. The Fermi surface is
then a prolate spheroid with its centre at the X point, major axis pointing
along X$\Gamma$ and minor axes in the XW and XU plane. The first Brillouin
zone with the Fermi surface highlighted is shown in the inset to Figure 2. It
is then the dispersion along these three high symmetry directions which will
dictate the transport properties of the material.

To represent the above description Figure 2 shows the calculated band
structure of nitrogen vacancy doped GdN along these high symmetry directions.
It is clear that the curvature of the bands contrasts significantly, resulting
in a differing effective mass, and ultimately transport properties along
different high symmetry directions. The figure shows the calculated wave-
functions projected onto the atomic orbitals (colours noted in caption), in
addition we have calculated the dispersion for the stoichiometric Gd 5d
conduction band along the same high symmetry directions, shown for the
majority spin as black dashed lines in the left hand panel of Figure 2. It is
significant that the stoichiometric and nitrogen vacancy doped conduction
bands are essentially superimposed on one another, fitting the picture that
the electrons donated by nitrogen vacancy doping are simply lifting the Fermi
level into the intrinsic conduction band. Although the shape of the bands is
retained between the stoichiometric and doped calculations the gap between the
Gd 5d and N 2p states at the X point has reduced from $\sim$ 0.9 eV in the
stoichiometric case to $\sim$ 0.6 eV in the doped case, the renormalzation
likely a result of screening from the delocalised electrons. The onset of
optical absorption is still near $\sim$ 0.9 eV, however now between $\Gamma$
and X, and expected to be much weaker as the initial states are the occupied
defect states, here with a concentration of only 1/27 to that of the N 2p. The
intrinsic majority spin optical band gap (i.e. not involving defect states)
increases to $\sim$ 1.2 eV ($\sim$ 0.9 eV in the stoichiometric crystal) in a
Moss-Berstein like shift.

The minority spin states are shown in the right hand panel of Figure 2. The Gd
5d states are again very similar to the stoichiometric crystal. However, the
presence of the occupied defect states in the gap, below the Fermi level,
effectively reduces minority spin optical gap from 2.5 eV in the
stoichiometric to 1.8 eV in the doped crystal. The average absorption onset of
the majority and minority spin gaps in the doped crystal (both with initial
impurity states and Gd 5d final states) is now $\sim$ 1.35 eV (1.75 eV in the
stoichiometric crystal). This is now less than is commonly reported for the
paramagnetic band gap.

Using the present results from the doped crystal the spin splitting can be
reported in various ways. The most appropriate for electrical transport (i.e.
as in a spin filter) is the minimum energy difference between available (hole
or electron) majority and minority spin states. As Fermi level lies in the
majority spin Gd 5d band, for electron carriers this is the energy difference
between the Fermi level and the bottom of the minority spin conduction band,
roughly 1.1 eV, while for hole carriers this is the energy difference between
the Fermi level and the top of the minority spin valance band, roughly 0.6 eV.
For optical studies the most relevant spin splitting is the difference between
the majority spin and minority spin optical gaps, which now $\sim$0.9 eV. This
is closer to the experimental reports, yet still beyond the maximum 0.6 eV
reported in reference Vilela  _et al._ (2024).

###  III.2 Optical spectroscopy

Figure 3: Calculated optical conductivity for pristine (top panel) and
nitrogen vacency doped (bottom panel) GdN. Both majority (black) and minority
(red) spin contributions are shown. The insets show the majority spin
contribution and in addition experimental data collected at 8 K.

As an initial comparison of our computational results to experiment we have
calculated the optical conductivity for both the pristine and doped cells
discussed above. These are displayed in Figure 3. The top panel of Figure 3
shows the calculated optical conductivity for stoichiometric GdN. The onset of
optical absorption is in the majority spin states at $\sim$ 0.9 eV and
initially shows an increase with a near $A+(\epsilon-E_{g})^{1/2}$ form, which
is expected from the Van Hove singularity at the conduction band minimum.
Above 2 eV the gradient increases relating to the flattening of both the 5d
and the N 2p bands near +2 eV and -2 eV respectively. At roughly 2.5 eV there
is a similar onset of absorption in the minority spin states. The inset show
the spectra between 0.8 eV and 3 eV now with experimental data (blue) measured
at 8 K. Although the onset of absorption is a roughly the same energy in both
spectra the form is very different, with the experiment not taking the typical
shape expected for the onset of optical transitions in a parabolic band.

The lower panel of Figure 3 shows the calculated optical conductivity for the
nitrogen vacancy doped super-cell. Here, as discussed in section III.1, the
absorption onset is at roughly the same energy (0.9 eV) as the pristine cell.
For the doped material the form of the optical conductivity differs greatly
from the simple $\epsilon^{1/2}$ shape as the Fermi level now lies within the
conduction band. The inset again shows the comparison between the calculation
and experiment with a clear qualitative similarity between the two (the
calculation has been scaled by a factor of 0.2 to best match the experiment).
The film used for the measurement was moderately conductive with a resistivity
of 0.1 $\Omega$cm at room temperature, increasing to 10 $\Omega$cm at 4 K.
Given this resistivity, and based on Hall measurements on other films Lee  _et
al._ (2015), one could assume a carrier concentration on the order 10 17cm-3
at 4 K, however the striking resemblance between the optical spectra for this
film and the super-cell calculation (carrier concentration $\sim 1~{}\times
10^{21}$cm-3) may rather imply that in this is incorrect. This film and
possibly other nominally undoped rare-earth nitrides may harbour a significant
population of localised carriers. These do not contribute to electron
transport in the usual extended state manner; however, they do occupy states
at the bottom of the conduction band, affecting the optical spectra and other
related electrical properties. Similar optical spectra have been seen in other
rare-earth nitrides Holmes-Hewett  _et al._ (2019, 2020, 2023) and along with
electrical transport measurements were discussed in the context of localised
hopping type conductivity and a mobility edge scenario.

###  III.3 Effective mass

Although the details of the spin-splitting and related exchange energy remain
to be determined the present data, calculations and existing experimental
reports Trodahl  _et al._ (2017) seem to agree that when GdN is doped with
electrons via nitrogen vacancies the Fermi level is raised into the Gd 5d
band, which is majority spin only at dilute to moderate doping levels.
Motivated by these observations, and the similarity between the stoichiometric
and nitrogen vacancy doped Gd 5d bands we have used the former to calculate
the effective mass along the X$\Gamma$, XW and XU directions, which may be of
use in the context of device design. The results are shown in Figure 4 with
the dispersion along the X$\Gamma$, XW and XU directions in the top panel, and
the effective mass in the bottom panel, both as a function of wave-vector
measured relative to the X point at $k=0$.

At the X point the first Brillouin zone of the FCC lattice has four-fold
symmetry (see inset of Figure 2), which dictates that along any two directions
in the plane of the X-face the effective mass is degenerate at $k=$X. Without
loss of generality we choose directions along XU and XW to form our basis, the
final perpendicular direction is then along X$\Gamma$ and will have a distinct
effective mass at the X point. The $k=$X effective mass tensor is diagonal in
this coordinate system and given by

$\frac{m^{*}}{m_{e}}=\begin{pmatrix}0.18&0&0\\\ 0&0.18&0\\\
0&0&1.8\end{pmatrix}.$

The enhanced curvature along the XU and XW directions results in an order of
magnitude effective mass difference between the longitudinal and transverse
directions of the Fermi surface. Using these values we can estimate the
conductivity effective mass as $m^{*}_{c}=0.25~{}m_{e}$, which is dominated by
the low mass components. Given the relative orientation between the equivalent
X points in the first Brillouin zone an electrical field along any direction
will result in contributions from both the transverse the longitudinal
directions. The order of magnitude lower mass in the transverse direction will
result in this dominating in any electrical transport measurement.

The effective mass precisely at the X point is useful in the context of
tunnelling experiments where the GdN is likely intended to be insulating and
close to stoichiometric, in addition the geometry of these experiments is well
controlled and for example the common case of NbN/GdN/NbN multi-layers
Senapati  _et al._ (2011) results in transport through the GdN along a 001
direction. As discussed above the degeneracy of the X point results in
contributions from both the transverse and longitudinal directions; at the X
point the combination of these contributions can in principle be easily
treated.

As one dopes the crystal with electrons the Fermi surface expands from the X
point and the effective mass is a function of wavevector as shown in Figure 4.
Along the X$\Gamma$ direction the effective mass is relatively stable with a
slight drop before an increase near $k=\pi/a$. As one continues along this
direction there is ultimately a sign inversion near the first Brillouin zone
boundary as the curvature of the band itself changes sign.

Due to the rapidly changing mass along the two transverse directions it may be
difficult to predict the transport properties of even moderately electron
doped GdN. The super-cell used in the calculation of Figure 2 was constructed
from $3\times 3\times 3$ primitive unit cells with a N atom removed. The
single mobile electron in a volume of $\sim 915\mathrm{\mathring{A}}^{3}$
results in a carrier concentration of $\sim$1$\times 10^{21}$cm-3 and is
representative of the heavily doped GdN from experimental reports Maity  _et
al._ (2018). The presence of this electron raises the Fermi level roughly 0.6
eV into the conduction band, the resulting Fermi wave vectors are indicated on
Figure 4. Far from the X point the effective mass tensor will no longer take
the simple diagonal form presented above, but even so the wave vector
dependent effective mass shown in Figure 4 gives a qualitative indication of
the evolution of the effective mass with electron doping in GdN.

Figure 4: Calculated dispersion (top) and effective mass (bottom) for GdN as a
function of wave vector along various high symmetry lines with $k=$X$=0$. The
Fermi wave vector for each dispersion, with reference to Figure 2, is
indicated with the dashed vertical lines.

##  IV Conclusions

We have undertaken a computational study of the band structure and defect
states in stoichiometric and nitrogen vacancy doped GdN. Our most direct
result is the effective mass tensor at the X point, and calculation of the
effective mass in the high symmetry directions moving away from the centre of
this pocket. We have provided updated band structure calculations of GdN and
compared these to conventional optical spectroscopy. The comparison between
conventional spectroscopy and the calculation indicates that even nominally
dilutly doped GdN has evidence of a large population of localised carriers in
occupying states at the conduction band minimum.

##  V Data Availability

The data used during this study are available from the corresponding author
upon reasonable request.

##  VI ACKNOWLEDGMENTS

This research was supported by Quantum Technologies Aotearoa, a research
programme of Te Whai Ao – the Dodd Walls Centre, funded by the New Zealand
Ministry of Business Innovation and Employment through International Science
Partnerships, contract number UOO2347. The computations were performed on the
Rāpoi high performance computing facility of Victoria University of
Wellington.

## References

  * Andrae (2015) T. Andrae, A.S.G.; Edler, Challenges 6, 117 (2015).
  * IRD (2018a) _International roadmap for devices and systems - Beyond CMOS_ (2018).
  * IRD (2018b) _International roadmap for devices and systems - More Moore_ (2018).
  * Hirohata  _et al._ (2020) A. Hirohata, K. Yamada, Y. Nakatani, I.-L. Prejbeanu, B. Diény, P. Pirro, and B. Hillebrands, Journal of Magnetism and Magnetic Materials 509, 166711 (2020).
  * Bhatti  _et al._ (2017) S. Bhatti, R. Sbiaa, A. Hirohata, H. Ohno, S. Fukami, and S. Piramanayagam, Materials Today 20, 530 (2017).
  * Soloviev  _et al._ (2017) I. I. Soloviev, N. V. Klenov, S. V. Bakurskiy, M. Y. Kupriyanov, A. L. Gudkov, and A. S. Sidorenko, Beilstein Journal of Nanotechnology 8, 2689 (2017).
  * Likharev and Semenov (1991) K. Likharev and V. Semenov, IEEE Transactions on Applied Superconductivity 1, 3 (1991).
  * Li and Yang (2016) X. Li and J. Yang, National Science Review 3, 365 (2016).
  * Shamiul Alam (2023) S. R. S. . A. A. Shamiul Alam, Md Shafayat Hossain, Nature Electronics 6, 185–198 (2023).
  * Natali  _et al._ (2010) F. Natali, N. O. Plank, J. Galipaud, B. J. Ruck, H. J. Trodahl, F. Semond, S. Sorieul, and L. Hirsch, Journal of Crystal Growth 312, 3583 (2010).
  * Eick  _et al._ (1956) H. A. Eick, N. C. Baenziger, and L. Eyring, Journal of the American Chemical Society 78, 5987 (1956).
  * Hulliger (1978) F. Hulliger, Journal of Magnetism and Magnetic Materials 8, 183 (1978).
  * Hulliger (1979) F. Hulliger, in  _Non-Metallic Compounds - II_, Handbook on the Physics and Chemistry of Rare Earths, Vol. 4 (Elsevier, 1979) pp. 153 – 236.
  * Granville  _et al._ (2006) S. Granville, B. J. Ruck, F. Budde, A. Koo, D. J. Pringle, F. Kuchler, A. R. H. Preston, D. H. Housden, N. Lund, A. Bittar, G. V. M. Williams, and H. J. Trodahl, Phys. Rev. B 73, 235335 (2006).
  * Aerts  _et al._ (2004) C. M. Aerts, P. Strange, M. Horne, W. M. Temmerman, Z. Szotek, and A. Svane, Phys. Rev. B 69, 045115 (2004).
  * Anton  _et al._ (2016) E.-M. Anton, J. F. McNulty, B. J. Ruck, M. Suzuki, M. Mizumaki, V. N. Antonov, J. W. Quilty, N. Strickland, and H. J. Trodahl, Phys. Rev. B 93, 064431 (2016).
  * Holmes-Hewett  _et al._ (2020) W. F. Holmes-Hewett, C. Pot, R. G. Buckley, A. Koo, B. J. Ruck, F. Natali, A. Shaib, J. D. Miller, and H. J. Trodahl, Appl. Phys. Lett. 117, 222409 (2020).
  * Holmes-Hewett  _et al._ (2023) W. F. Holmes-Hewett, K. V. Koughnet, J. D. Miller, E. X. M. Trewick, B. J. Ruck, H. J. Trodahl, and R. G. Buckley, Sci Rep 13 (2023).
  * Devese  _et al._ (2022) S. Devese, K. Van Koughnet, R. G. Buckley, F. Natali, P. P. Murmu, E.-M. Anton, B. J. Ruck, and W. F. Holmes-Hewett, AIP Advances 12, 035108 (2022).
  * Ludbrook  _et al._ (2009) B. Ludbrook, I. Farrell, M. Kübel, B. Ruck, A. Preston, J. Trodahl, L. Ranno, R. Reeves, and S. Durbin, Journal of Applied Physics 106, 063910 (2009).
  * Trodahl  _et al._ (2007) H. J. Trodahl, A. R. H. Preston, J. Zhong, B. J. Ruck, N. M. Strickland, C. Mitra, and W. R. L. Lambrecht, Phys. Rev. B 76, 085211 (2007).
  * Larson and Lambrecht (2006) P. Larson and W. R. L. Lambrecht, Phys. Rev. B 74, 085108 (2006).
  * Larson  _et al._ (2007) P. Larson, W. R. L. Lambrecht, A. Chantis, and M. vanSchilfgaarde, Phys. Rev. B 75, 045114 (2007).
  * Preston  _et al._ (2010) A. R. H. Preston, B. J. Ruck, W. R. L. Lambrecht, L. F. J. Piper, J. E. Downes, K. E. Smith, and H. J. Trodahl, Appl. Phys. Lett. 96, 032101 (2010).
  * Vilela  _et al._ (2024) G. L. S. Vilela, G. M. Stephen, X. Gratens, G. D. Galgano, Y. Hou, Y. Takamura, D. Heiman, A. B. Henriques, G. Berera, and J. S. Moodera, Phys. Rev. B 109, L060401 (2024).
  * Azeem (2016) M. Azeem, Chin. Phys. Lett. 33, 027501 (2016).
  * Yoshitomi  _et al._ (2011) H. Yoshitomi, S. Kitayama, T. Kita, O. Wada, M. Fujisawa, H. Ohta, and T. Sakurai, Phys. Rev. B 83, 155202 (2011).
  * Vidyasagar  _et al._ (2012) R. Vidyasagar, S. Kitayama, H. Yoshitomi, T. Kita, T. Sakurai, and H. Ohta, Appl. Phys. Lett. 100, 232410 (2012).
  * Maity  _et al._ (2018) T. Maity, H. J. Trodahl, F. Natali, B. J. Ruck, and S. Vézian, Phys. Rev. Materials 2, 014405 (2018).
  * Punya  _et al._ (2011) A. Punya, T. Cheiwchanchamnangij, A. Thiess, and W. Lambrecht, MRS Proceedings 1290 (2011).
  * Trodahl  _et al._ (2017) H. J. Trodahl, F. Natali, B. J. Ruck, and W. R. L. Lambrecht, Phys. Rev. B 96, 115309 (2017).
  * Birge and Satchell (2024) N. O. Birge and N. Satchell, APL Materials 12, 041105 (2024).
  * Vernik  _et al._ (2013) I. V. Vernik, V. V. Bol’ginov, S. V. Bakurskiy, A. A. Golubov, M. Y. Kupriyanov, V. V. Ryazanov, and O. A. Mukhanov, IEEE Transactions on Applied Superconductivity 23, 1701208 (2013).
  * Senapati  _et al._ (2011) K. Senapati, M. G. Blamire, and Z. H. Barber, Nat. Mater. 10, 849 (2011).
  * Massarotti  _et al._ (2015) D. Massarotti, A.Pal, G. Rotoli, L. Longobardi, M. Blamire, and F. Tafuri, Nature Communications 6, 7376 (2015).
  * Ahmad  _et al._ (2020) H. Ahmad, R. Caruso, A. Pal, G. Rotoli, G. Pepe, M. Blamire, F. Tafuri, and D. Massarotti, Phys. Rev. Appl. 13, 014017 (2020).
  * Sharma  _et al._ (2023) P. K. Sharma, S. Banerjee, B. Dutta, V. Singhal, P. Banerjee, H. K. Pal, and A. Pal, (2023), arXiv:2312.04650 [cond-mat.supr-con] .
  * Cascales  _et al._ (2019) J. P. Cascales, Y. Takamura, G. M. Stephen, D. Heiman, F. S. Bergeret, and J. S. Moodera, Applied Physics Letters 114, 022601 (2019).
  * Pot  _et al._ (2023) C. Pot, W. F. Holmes-Hewett, E.-M. Anton, J. D. Miller, B. J. Ruck, and H. J. Trodahl, Appl. Phys. Lett. 13 (2023).
  * Ahmad  _et al._ (2022) H. G. Ahmad, V. Brosco, A. Miano, L. Di Palma, M. Arzeo, D. Montemurro, P. Lucignano, G. P. Pepe, F. Tafuri, R. Fazio, and D. Massarotti, Phys. Rev. B 105, 214522 (2022).
  * Mitra and Lambrecht (2008) C. Mitra and W. R. L. Lambrecht, Phys. Rev. B 78, 195203 (2008).
  * Giannozzi  _et al._ (2009) P. Giannozzi _et al._ , J. Phys.: Condens. Matter 21, 395502 (2009).
  * Cococcioni and de Gironcoli (2005) M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 035105 (2005).
  * Topsakal and Wentzcovitch (2014) M. Topsakal and R. Wentzcovitch, Comput. Mater. Sci. 95, 263 (2014).
  * Mostofi  _et al._ (2014) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 185, 2309 (2014).
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  * Souza  _et al._ (2001) I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001).
  * Holmes-Hewett (2021) W. F. Holmes-Hewett, Phys. Rev. B 104, 075124 (2021).
  * Ehrenreich  _et al._ (1963) H. Ehrenreich, H. R. Philipp, and B. Segall, Phys. Rev. 132, 1918 (1963).
  * Kuzmenko (2005) A. B. Kuzmenko, Rev. Sci. Instrum. 76, 083108 (2005).
  * Lee  _et al._ (2015) C. M. Lee, H. Warring, S. Vézian, B. Damilano, S. Granville, M. Al Khalfioui, Y. Cordier, H. J. Trodahl, B. J. Ruck, and F. Natali, Appl. Phys. Lett. 106, 022401 (2015).
  * Holmes-Hewett  _et al._ (2019) W. F. Holmes-Hewett, R. G. Buckley, B. J. Ruck, F. Natali, and H. J. Trodahl, Phys. Rev. B 99, 205131 (2019).