M}^{1}}+\mathrm{i}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\phi_{\rm
M}^{2}})$ [12]. Then, to the order $\bm{p}^{2}$, we have the following
quadratic Lagrangian

$L=c_{1}\partial_{t}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}\partial_{t}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}-c_{2}\nabla{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}\cdot\nabla{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}+\mathrm{i}c_{3}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}^{*}}\partial_{t}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}}-c_{4}\nabla{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}^{*}}\cdot\nabla{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}}-Bc_{5}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}^{*}}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}}$ (69)

where $B$ denotes intensity of external magnetic field $\bm{B}=B\bm{T}_{3}$
which plays the role of chemical potential. The therm proportional to magnetic
field is
$({\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}})^{*}\langle\bm{\Sigma}(z),\bm{B}\rangle$, with $\Sigma_{i}(z)$ expanded up
to quadratic terms using (60). Also, $c_{1},c_{2},c_{3},c_{4}$ and $c_{5}$ are
unspecified constants whose values may be determined by matching with linear
spin-wave theory or numerical simulations.

Figure 1: Canted spin configuration on a square lattice – the symmetry of
ground state configuration is trivial, while emergent magnetization
$\bm{\Sigma}$ is invariant under the action of $\mathrm{SO}(2)$. Small
deviations from the ground state are described by NG fields: Each point in the
field space is a bundle so that massive NG fields
${\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\mathrm{M}}}$
take values on sphere $S^{2}$, while type A NG field
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}$
lives on $S^{1}$.

####  VI.3.2 Relativistic Bose-Einstein condensation

As a second example, consider a
$G=\mathrm{O}(4)\simeq\mathrm{SU}(2)\times\mathrm{SU}(2)$ invariant theory of
complex scalar doublet $\psi=[\psi^{1}\;\;\psi^{2}]^{\mathrm{T}}$ with
standard Lagrangian [27, 68]

$L=\partial_{\alpha}\psi^{\dagger}\partial^{\alpha}\psi-m^{2}\psi^{\dagger}\psi-\lambda\bigl{(}\psi^{\dagger}\psi\bigr{)}^{2}.$
(70)

Explicit symmetry breakdown in Lagrangian formalism is achieved by
substitution $\partial_{0}\to D_{0}=\partial_{0}-\mathrm{i}\mu$, where $\mu$
is the chemical potential, and this reduces the symmetry of Lagrangian to
$\widetilde{G}=\mathrm{SU}(2)\times\mathrm{U}(1)$. If $\mu>m$, this model
exhibits relativistic Bose-Einstein condensation and is used to study kaons in
dense quark matter. We shall take commutation relations between elements of
Lie algebra in this subsection to be
$[\bm{T}_{i},\bm{T}_{j}]=\mathrm{i}f_{ij}^{\;\;k}\bm{T}_{k}$.

The classical ground state of this model is $\psi_{0}=[0\;\;v]^{\mathrm{T}}$.
To find its symmetry group, we can choose the basis for
$\tilde{\mathfrak{g}}=\mathfrak{su}(2)\oplus\mathfrak{u}(1)$ as
$\\{\bm{T}_{i}=\bm{\sigma}_{i}/2,\bm{T}_{4}=I_{2\times 2}/2\\}$, where
$\bm{\sigma}_{i}$ are the Pauli matrices, and examine the action of an
arbitrary element of $\mathfrak{g}$

$\bm{A}=\sum_{i=1}^{4}A^{i}\bm{T}_{i}$ (71)

on $\psi_{0}$. The invariance condition $\bm{A}\psi_{0}=0$ yields
$A^{1}=A^{2}=0$ and $A^{3}=A^{4}$, leading to $\mathfrak{h}=\mathfrak{u}(1)$
and consequently $H=\mathrm{U}(1)$. Having determined $H$, we focus on the
order parameter. Giving the fact that the Lagrangian contains the term with
single time derivative
$\mathrm{i}\mu\left(\psi^{\dagger}\dot{\psi}-\dot{\psi}^{\dagger}\psi\right)$,
and that under the action of $\bm{A}$ the field transforms ae
$\delta\psi^{B}=\mathrm{i}A^{k}[\bm{T}_{k}]_{C}^{\;\;B}\psi^{C}\equiv\mathrm{i}A^{k}\mathcal{F}^{B}_{k}$,
where $C,B=1,2$ and $[\bm{T}_{k}]_{C}^{\;\;B}$ denoting matrix elements of
generator $\bm{T}_{k}$, we find the time component of conserved currents

$J^{0}_{k}(x)=-\mathrm{i}\frac{\partial\mathcal{L}}{\partial\dot{\psi}^{A}(x)}\mathcal{F}^{A}_{k}(x)=\mu\psi^{*}_{A}(x)\bigl{[}\bm{T}_{k}\bigr{]}_{B}^{\;\;A}\psi^{B}(x)+\dots$
(72)

which, for the ground state configuration $\psi_{0}$, gives
$J^{0}_{1}=J^{0}_{2}=0$, $J^{0}_{3}=-J^{0}_{4}=-\mu v^{2}$. Therefore, the
order parameter is given by

$\bm{\Sigma}=\mu v^{2}\biggl{(}-\bm{T}^{*3}+\bm{T}^{*4}\biggr{)}$ (73)

To identify $\widetilde{\mathfrak{stab}(\bm{\Sigma})}$, we start from the
definition (24) and find all $\bm{B}\in\mathfrak{su}(2)\oplus\mathfrak{u}(1)$
such that, for arbitrary $\bm{A}\in\mathfrak{su}(2)\oplus\mathfrak{u}(1)$ the
following holds

$0=\langle\mathrm{ad}^{*}(\bm{B})\bm{\Sigma},\bm{A}\rangle=\left\langle\bm{\Sigma},-\left[\bm{B},\bm{A}\right]\right\rangle=-\mathrm{i}f_{ij}^{\;\;k}B^{i}A^{j}\left\langle\bm{\Sigma},\bm{T}^{k}\right\rangle.$
(74)

which leads to $B^{1}=B^{2}=0$ and we get
$\widetilde{\mathrm{Stab}(\bm{\Sigma})}=\mathrm{U}(1)\times\mathrm{U}(1)$. By
using similar procedure we find
$\mathrm{Stab}(\bm{\Sigma})=\widetilde{\mathrm{Stab}(\bm{\Sigma})}$.
Therefore, spectrum of this model contains one massless NG field of type A,
two massless NG fields of type B (carrying one physical degree of freedom), as
well as two massive NG fields (carrying one physical degree of freedom as
well). This corresponds to one massless NG boson of type A, one of type B and
one massive NG boson in quantized theory. Indeed, this is found by explicit
calculations [68, 27, 69].

Effective Lagrangian for this model, which captures dynamics of non-
interacting NG fields up to $\bm{p}^{2}$, may be constructed in a similar
manner as (69). It is given by

$L=c_{1}\partial_{t}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}\partial_{t}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}-c_{2}\nabla{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}\cdot\nabla{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}+\mathrm{i}c_{3}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}^{*}}\partial_{t}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}}-c_{4}\nabla{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}^{*}}\cdot\nabla{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}}-c_{5}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}^{*}}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}}+\mathrm{i}c_{6}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\psi_{\rm
B}^{*}}\partial_{t}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\psi_{\rm
B}}-c_{7}\nabla{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\psi_{\rm
B}^{*}}\cdot\nabla{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\psi_{\rm
B}}$ (75)

where $c_{1}c_{2},\dots c_{7}$ are arbitrary constants with $c_{5}\propto\mu$
(As in previous example, we have absorbed $|\bm{\Sigma}|\propto\mu$ into the
definition of complex fields
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\psi_{\mathrm{B}}}$
and
${\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\psi_{\rm
M}}$). This term also stems from $\langle\bm{\Sigma}(z),\bm{\mu}\rangle$ and
can be obtained using (60) and structure constants of
$\mathfrak{su}(2)\oplus\mathfrak{u}(1)$, which vanish whenever $i,j$ or $k=4$.

####  VI.3.3 The $\mathrm{SO}(2N)$ vector model at fixed charge

For our final example, consider a Lorentz-invariant model with additional
$G=\mathrm{SO}(2N)$ internal symmetry. In what follows, we shall use the fact
that coadjoint orbits for compact groups are equivalent with adjoint ones,
where the identification between elements of $\mathfrak{g}$ and
$\mathfrak{g}^{*}$ can be done using an Ad-invariant inner product on
$\mathfrak{g}$ [42]. Suppose that the Lorentz symmetry, together with
$\mathrm{SO}(2N)$ internal symmetry, is explicitly broken by adding
$-\mu\left(Q_{1}+Q_{2}+\dots+Q_{k}\right)$, $k\leq N$ to Hamiltonian, with
$Q_{i}$ may now be viewed as basis elements of $\mathfrak{so}(2N)$. Since the
same chemical potential is coupled to all $k$ generators, we find [43]

$\widetilde{G}=\mathrm{Stab}\biggl{(}Q_{1}+Q_{2}+\dots+Q_{k}\biggr{)}=\mathrm{SO}(2N-2k)\times\mathrm{U}(k).$
(76)

The $k$ fields transform under $U(k)$ and their configurations take part in
subsequent spontaneous symmetry breakdown while the remaining $2N-2k$ real
field are spectators. Following [70] we take the symmetry of ground state
configuration to be $H=\mathrm{SO}(2N-2k)\times\mathrm{U}(k-1)$. With Lorentz
symmetry being lost, one of the $k$ generators may develop a vacuum
expectation value $\Sigma\propto\langle Q_{1}\rangle_{\mathrm{g.s.}}\neq 0$
and $\bm{\Sigma}=\langle Q_{1}\rangle_{\mathrm{g.s.}}\bm{T}^{*1}$. Since
$\mathfrak{g}^{*}$ and $\mathfrak{g}$ have been identified, we have [43]

$\mathrm{Stab}(\bm{\Sigma})=\mathrm{Stab}\bigl{(}Q_{1}\bigr{)}=\mathrm{SO}(2N-2)\times\mathrm{U}(1)$
(77)

under the action of $G$. Likewise, the order parameter stabilizer for the
action of $\mathrm{U}(k)$ is $\mathrm{U}(k-1)\times\mathrm{U}(1)$ [39], thus

$\widetilde{\mathrm{Stab}(\bm{\Sigma})}=\mathrm{SO}(2N-2k)\times\mathrm{U}(k-1)\times\mathrm{U}(1).$
(78)

Having specified all relevant symmetry groups, we can analyze the spectrum of
NG fields in this model for explicit $G\to\widetilde{G}$ and subsequent
spontaneous symmetry breakdown $\widetilde{G}\to H$.

First, there is
$\mathrm{dim}\left(\widetilde{\mathrm{Stab}(\bm{\Sigma})}/H\right)=1$ NG field
of type A. Second, the number of physical degrees of freedom corresponding to
type B NG fields is given by (62) and this yields
$\widetilde{N}_{\mathrm{B}}=k-1$. Finally, the number of degrees of freedom
carried by massive fields can be read off from (63), which gives $2N-k-1$.
However, $2N-2k$ of these massive modes are, by construction, spectator fields
and thus ”frozen” under the action of full group $G$. Therefore, the number of
massive NG fields is

$N_{\mathrm{mNGF}}=k-1$ (79)

In other words, the NG sector of corresponding quantum theory contains one NG
boson of type A, $k-1$ NG bosons of type B, $k-1$ massive NG bosons and the
full spectrum also contains $2N-2k$ massive particles corresponding to
spectator fields. The results presented here are in complete agreement with
detailed calculations performed in [70] (See also [71, 72]). As in the
previous example, a direct calculation reveals the presence of a one more
massive degree of freedom corresponding to the Higgs mode which is
inaccessible to this construction. The effective Lagrangian which describes
dynamics of non-interacting NG fields to $\bm{p}^{2}$ in this case is obtained
as a direct generalization of (75)

$\displaystyle L$ $\displaystyle=$ $\displaystyle
c_{1}\partial_{\mu}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}\partial^{\mu}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}+\frac{c_{2}}{2}\Sigma_{1}f_{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\alpha}}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\beta}}}^{1}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\phi_{\mathrm{B}}^{\alpha}}}\partial_{t}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\phi_{\mathrm{B}}^{\beta}}}-c_{3}\nabla{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\phi_{\rm
B}^{\alpha}}\cdot\nabla{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\phi_{\rm
B}^{\alpha}}$ (80) $\displaystyle+$
$\displaystyle\frac{c_{4}}{2}\Sigma_{1}f_{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\alpha\beta}^{1}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}{\phi_{\mathrm{M}}^{\alpha}}}\partial_{t}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}{\phi_{\mathrm{M}}^{\beta}}}-c_{5}\nabla{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\phi_{\rm
M}^{\alpha}}\cdot\nabla{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\phi_{\rm
M}^{\alpha}}+c_{6}\frac{\mu^{l}}{2}\Sigma_{1}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}{\phi_{\mathrm{M}}^{\alpha}}}{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}{\phi_{\mathrm{M}}^{\beta}}}f_{{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\alpha}l}^{\;\;j}f_{{\color[rgb]{0,0.9,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.9,0.8}\beta}j}^{\;\;1}$

where $\alpha,\beta=1,2,\dots k-1$ for all fields and
$\mu_{1}=\mu_{2}\dots=\mu_{k}\equiv\mu$. Note that $G/\widetilde{G}$ is a
subspace of $G/\widetilde{\mathrm{Stab}(\bm{\Sigma})}$ and thus inherits the
symplectic structure defined by $\bm{\Sigma}$. Hadn’t we chosen $2N-2k$
spectator fields, the spectrum would contain one NG boson of type A, $k-1$ NG
bosons of type B, $2N-k-1$ massive NG bosons and one Higgs boson. The total
number of particles is the same in both cases and equals $2N$.

##  VII Sumarry

Nambu-Goldstone fields, the low energy degrees of freedom, play prominent role
in physical systems exhibiting spontaneous symmetry breakdown. Their dynamics
follows from the effective action which is the most general $G$ invariant
functional defined on an associated bundle whose typical fiber is $G/H$. While
the form of effective Lagrangian, which involves local coordinates on $G/H$,
is needed to derive equations of motion and detailed prediction of theory
based on perturbative calculations, some important characteristics of NG
fields follow directly from the symmetry breaking pattern and the symmetry of
the order parameter $\bm{\Sigma}\in\mathfrak{g}^{*}$. A novel feature of the
present work is precisely the coordinate-free formulation of the first-order
effective field theory which allows for a direct derivation of the number of
physical degrees of freedom assigned to type A and type B Nambu-Goldstone
fields, including so-called massive NG fields.

The construction of classical field theory for type A and B NG fields is
implemented in several steps. By assuming symmetry breaking pattern $G\to H$,
we have shown that classical field theory in which NG fields are maps
$\phi:\mathrm{M}\to G/H$ could be cast in the language of jets where NG fields
are identified as sections on associated bundle
$(\mathrm{E},\pi_{\mathrm{E}},\mathrm{M})$ whose fibers are $G/H$. Next, with
nonzero Lorentz-symmetry breaking order parameter
$\bm{\Sigma}\in\mathfrak{g}^{*}$ and $H\subset\mathrm{Stab}(\bm{\Sigma})$, the
coset space $G/H$ becomes a bundle, the base space of which is the symplectic
manifold $G/\mathrm{Stab}(\bm{\Sigma})$ and the fibers are
$\mathrm{Stab}(\bm{\Sigma})/H$. By using the structure of associated bundle,
we demonstrated that to each field configuration $\phi(\mathrm{p})\in G/H$
corresponds a pair of points:
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\phi_{\mathrm{B}}}(\mathrm{p})\in
G/\mathrm{Stab}(\bm{\Sigma})$ and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}(\mathrm{p})\in\mathrm{Stab}(\bm{\Sigma})/H$.
The maps
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\phi_{\mathrm{B}}}:\mathrm{M}\to
G/\mathrm{Stab}(\bm{\Sigma})$ and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}:\mathrm{M}\to\mathrm{Stab}(\bm{\Sigma})/H$
are local representatives of classical type B and type A NG fields and the
number of physical degrees of freedom carried by both types of fields is
determined by $G,H$ and $\mathrm{Stab}(\bm{\Sigma})$. Besides directly
yielding correct number of physical degrees of freedom, this approach also
provides a path to Lagrangian formulated in terms of
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi_{\mathrm{A}}}$
and
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\phi_{\mathrm{B}}}$.
We also examine the case of explicit symmetry breaking by finite charge
density in the ground state which reduces the symmetry of model from $G$ to
$\widetilde{G}$. By adapting arguments from first four sections to this case,
we demonstrated how effective field theory in geometric formulation may be
used to predict the number of degrees of freedom corresponding to type A, type
B and massive NG fields, latter of which are defined to be field
configurations obtained by the the action of elements of $G$ not contained in
$\widetilde{G}$. As the target space for massive NG fields is
$G/\widetilde{G}$, the base space of bundle $G/\mathrm{Stab}(\bm{\Sigma})$,
massive NG fields also couple through (the pullback of) KK symplectic form
generated by $\bm{\Sigma}$. In this sense, they turn out be quite similar to
type B NG fields. It is also possible that, under the action of
$\widetilde{G}$, the isotropy group $\mathrm{Stab}(\bm{\Sigma})$ reduces to
$\widetilde{\mathrm{Stab}(\bm{\Sigma})}$. The effect of this change is the
reduction in the number of massive NG fields by
$\mathrm{Stab}(\bm{\Sigma})/\widetilde{\mathrm{Stab}(\bm{\Sigma})}$. Finally,
we illustrate general ideas presented here with several examples.

The existence of Lorentz-symmetry breaking order parameter
$\bm{\Sigma}\in\mathfrak{g}^{*}$, which induces symplectic structure on
$G/\mathrm{Stab}(\bm{\Sigma})$, is a question which can not be completely
resolved within classical field theory discussed here and needs to be
addressed in corresponding quantum or statistical theory. The transition to
quantized theory would allow for detailed perturbative studies and geometric
methods presented here could also be of use when analyzing the effects of
interactions between NG bosons. All results presented in this paper are
unrelated to a choice of local coordinates on $G/H$, and, by extension, the
Lagrangian function, as well as to a choice of basis in $\mathfrak{g}$, which
makes this formulation general and applicable to various physical situations
ranging from high energy particle physics all the way to condensed matter and
statistical systems.

###### Acknowledgements.

The author would like to thank professor Julio Guerrero for useful
discussions. The author also acknowledges the generous hospitality of the
Department of Mathematics at the University of Jaén during the spring semester
of 2024.

## References

  * S. Weinberg [2010] S. Weinberg,  _The Quantum Theory of Fields, Vol. II_ (Cambridge University Press, 2010).
  * Burgess [2000] C. Burgess, Goldstone and pseudo-goldstone bosons in nuclear, particle and condensed-matter physics, Phys. Repts. 330, 193 (2000).
  * Burgess [2020] C. P. Burgess,  _Introduction to effective field theory_ (Cambridge University Press, 2020).
  * Brauner [2024] T. Brauner,  _Effective Field Theory for Spontaneously Broken Symmetry_ (Springer, 2024).
  * Hofmann [1999a] C. P. Hofmann, Spin-wave scattering in the effective Lagrangian perspective, Phys. Rev. B 60, 388 (1999a).
  * Pavaskar  _et al._ [2022] S. Pavaskar, R. Penco, and I. Z. Rothstein, An effective field theory of magneto-elasticity, SciPost Phys. 12, 155 (2022).
  * Gasser and Leutwyler [1985] J. Gasser and H. Leutwyler, Chiral perturbation theory: Expansions in the mass of the strange quark, Nucl. Phys. B 250, 465 (1985).
  * Gasser and Leutwyler [1984] J. Gasser and H. Leutwyler, Chiral perturbation theory to one loop, Ann. Phys. 158, 142 (1984).
  * Gasser and Leutwyler [1988] J. Gasser and H. Leutwyler, Spontaneously broken symmetries: Effective Lagrangians at finite volume, Nucl. Phys. B 307, 763 (1988).
  * Hofmann [2011] C. P. Hofmann, Spontaneous magnetization of an ideal ferromagnet: Beyond Dyson’s analysis, Phys. Rev. B 84, 064414 (2011).
  * Weinberg [2016] S. Weinberg, Effective field theory, past and future, International Journal of Modern Physics A 31, 1630007 (2016).
  * Radošević [2015] S. M. Radošević, Magnon–magnon interactions in O(3) ferromagnets and equations of motion for spin operators, Annals of Physics 362, 336 (2015).
  * Nielsen and Chadha [1976] H. Nielsen and S. Chadha, On how to count Goldstone bosons, Nuclear Physics B 105, 445 (1976).
  * Brauner [2010] T. B. Brauner, Spontaneous symmetry breaking and Nambu–Goldstone bosons in quantum many-body systems, Symmetry 2, 609 (2010).
  * Watanabe [2020] H. Watanabe, Counting rules of Nambu–Goldstone modes, Annual Review of Condensed Matter Physics 11, 169 (2020).
  * Leutwyler [1994a] H. Leutwyler, Nonrelativistic effective Lagrangians, Phys. Rev. D 49, 3033 (1994a).
  * Nambu [2004] Y. Nambu, Spontaneous breaking of lie and current algebras, J. Stat. Phys. 115, 7 (2004).
  * Schäfer  _et al._ [2001] T. Schäfer, D. Son, M. Stephanov, D. Toublan, and J. Verbaarschot, Kaon condensation and Goldstone’s theorem, Phys. Lett. B 522, 67 (2001).
  * Watanabe and Murayama [2012] H. Watanabe and H. Murayama, Unified Description of Nambu-Goldstone Bosons without Lorentz Invariance, Phys. Rev. Lett. 108, 251602 (2012).
  * Watanabe and Murayama [2014] H. Watanabe and H. Murayama, Effective Lagrangian for nonrelativistic systems, Phys. Rev. X 4, 031057 (2014).
  * Watanabe and Brauner [2011] H. Watanabe and T. Brauner, Number of Nambu-Goldstone bosons and its relation to charge densities, Phys. Rev. D 84, 125013 (2011).
  * Alonso  _et al._ [2016a] R. Alonso, E. E. Jenkins, and A. V. Manohar, Sigma models with negative curvature, Physics Letters B 756, 358 (2016a).
  * Alonso  _et al._ [2016b] R. Alonso, E. E. Jenkins, and A. V. Manohar, Geometry of the scalar sector, Journal of High Energy Physics 2016, 1 (2016b).
  * Dobado and Espriu [2020] A. Dobado and D. Espriu, Strongly coupled theories beyond the Standard Model, Progress in Particle and Nuclear Physics 115, 103813 (2020).
  * Hidaka  _et al._ [2015] Y. Hidaka, T. Noumi, and G. Shiu, Effective field theory for spacetime symmetry breaking, Phys. Rev. D 92, 045020 (2015).
  * Tanii [2014] Y. Tanii,  _Introduction to supergravity_ (Springer, 2014).
  * Watanabe  _et al._ [2013] H. Watanabe, T. Brauner, and H. Murayama, Massive Nambu-Goldstone bosons, Phys. Rev. Lett. 111, 021601 (2013).
  * Nicolis  _et al._ [2013] A. Nicolis, R. Penco, F. Piazza, and R. A. Rosen, More on gapped Goldstones at finite density: More gapped Goldstones, Journal of High Energy Physics 2013, 1 (2013).
  * Hayata and Hidaka [2015] T. Hayata and Y. Hidaka, Dispersion relations of Nambu-Goldstone modes at finite temperature and density, Phys. Rev. D 91, 056006 (2015).
  * Bleecker [2005] D. Bleecker,  _Gauge theory and variational principles_ (Dover, 2005).
  * Binz  _et al._ [1988] E. Binz, J. Sniatycki, and H. R. Fisher,  _Geometry of classical fields_ (North-Hoolland, 1988).
  * Watanabe and Murayama [2013] H. Watanabe and H. Murayama, Redundancies in Nambu-Goldstone Bosons, Phys. Rev. Lett. 110, 181601 (2013).
  * Bogers and Brauner [2018] M. P. Bogers and T. Brauner, Lie-algebraic classification of effective theories with enhanced soft limits, Journal of High Energy Physics 2018, 1 (2018).
  * Saunders [1989] D. J. Saunders,  _The Geometry of Jet Bundles_, London Mathematical Society Lecture Note Series (Cambridge University Press, 1989).
  * Aldaya [2021] V. Aldaya, Structure of gauge theories, The European Physical Journal Plus 136, 304 (2021).
  * Gotay  _et al._ [1998] M. J. Gotay, J. Isenberg, J. E. Marsden, and R. Montgomery,  _Momentum maps and classical relativistic fields. Part I: Covariant field theory_ (1998).
  * De Leon  _et al._ [2004] M. De Leon, D. Martin de Diego, and A. Santamaria-Merino, Symmetries in classical field theory, International Journal of Geometric Methods in Modern Physics 1, 651 (2004).
  * Aldaya and de Azcárraga [1978] V. Aldaya and J. A. de Azcárraga, Variational principles on rth order jets of fibre bundles in field theory, Journal of Mathematical Physics 19, 1869 (1978).
  * Arvanitoyeorgos [2003] A. Arvanitoyeorgos,  _An introduction to Lie groups and the geometry of homogeneous spaces_, Vol. 22 (American Mathematical Soc., 2003).
  * Kowalski and Szenthe [2000] O. Kowalski and J. Szenthe, On the existence of homogeneous geodesics in homogeneous Riemannian manifolds, Geometriae Dedicata 81, 209 (2000).
  * Kolár  _et al._ [2013] I. Kolár, P. W. Michor, and J. Slovák,  _Natural Operations in Differential Geometry_ (Springer, 2013).
  * Kirillov [2004] A. A. Kirillov,  _Lectures on the orbit method_, Vol. 64 (American Mathematical Soc., 2004).
  * Besse [2007] A. L. Besse,  _Einstein manifolds_ (Springer, 2007).
  * Guerrero  _et al._ [2004] J. Guerrero, J. L. Jaramillo, and V. Aldaya, Group-cohomology refinement to classify $G$-symplectic manifolds, Journal of Mathematical Physics 45, 2051 (2004).
  * Delacrétaz  _et al._ [2022] L. V. Delacrétaz, Y.-H. Du, U. Mehta, and D. T. Son, Nonlinear bosonization of Fermi surfaces: The method of coadjoint orbits, Phys. Rev. Res. 4, 033131 (2022).
  * Craig and Lee [2024] N. Craig and Y.-T. Lee, Effective field theories on the jet bundle, Phys. Rev. Lett. 132, 061602 (2024).
  * Alminawi  _et al._ [2023] M. Alminawi, I. Brivio, and J. Davighi, Jet bundle geometry of scalar field theories, arXiv preprint 10.48550/arXiv.2308.00017 (2023).
  * Leutwyler [1994b] H. Leutwyler, On the foundations of chiral perturbation theory, Ann. Phys. 235, 165 (1994b).
  * Hidaka [2013] Y. Hidaka, Counting rule for Nambu-Goldstone modes in nonrelativistic systems, Phys. Rev. Lett. 110, 091601 (2013).
  * Hofmann [1999b] C. P. Hofmann, Effective analysis of the $\mathrm{O}(N)$ antiferromagnet: Low-temperature expansion of the order parameter, Phys. Rev. B 60, 406 (1999b).
  * Radošević  _et al._ [2013] S. M. Radošević, M. R. Pantić, M. V. Pavkov-Hrvojević, and D. V. Kapor, Magnon energy renormalization and low-temperature thermodynamics of O(3) Heisenberg ferromagnets, Ann. Phys. 339, 382 (2013).
  * A. Auerbach [1994] A. Auerbach,  _Interacting Electrons and Quantum Magnetism_ (Springer-Verlag, 1994).
  * Bajpai  _et al._ [2021] U. Bajpai, A. Suresh, and B. K. Nikolić, Quantum many-body states and Green’s functions of nonequilibrium electron-magnon systems: Localized spin operators versus their mapping to Holstein-Primakoff bosons, Phys. Rev. B 104, 184425 (2021).
  * Akyuz  _et al._ [2024] C. O. Akyuz, G. Goon, and R. Penco, The schwinger-keldysh coset construction, Journal of High Energy Physics 2024, 1 (2024).
  * Marsden and Ratiu [2013] J. E. Marsden and T. S. Ratiu,  _Introduction to mechanics and symmetry: a basic exposition of classical mechanical systems_, Vol. 17 (Springer Science & Business Media, 2013).
  * Coleman  _et al._ [1969] S. Coleman, J. Wess, and B. Zumino, Structure of Phenomenological Lagrangians. I, Phys. Rev. 177, 2239 (1969).
  * Callan  _et al._ [1969] C. G. Callan, S. Coleman, J. Wess, and B. Zumino, Structure of Phenomenological Lagrangians. II, Phys. Rev. 177, 2247 (1969).
  * Andersen  _et al._ [2014] J. O. Andersen, T. Brauner, C. P. Hofmann, and A. Vuorinen, Effective Lagrangians for quantum many-body systems, Journal of High Energy Physics 2014, 1 (2014).
  * Helgason [2001] S. Helgason,  _Differential geometry, Lie groups, and symmetric spaces_ (American Mathematical Society, 2001).
  * Casas and Murua [2009] F. Casas and A. Murua, An efficient algorithm for computing the Baker–Campbell–Hausdorff series and some of its applications, Journal of Mathematical Physics 50, 033513 (2009).
  * Araújo [2010] F. Araújo, Some Einstein homogeneous Riemannian fibrations, Differential Geometry and its Applications 28, 241 (2010).
  * Nicolis and Piazza [2012] A. Nicolis and F. Piazza, Spontaneous symmetry probing, Journal of High Energy Physics 2012, 1 (2012).
  * Cuomo  _et al._ [2021] G. Cuomo, A. Esposito, E. Gendy, A. Khmelnitsky, A. Monin, and R. Rattazzi, Gapped goldstones at the cut-off scale: a non-relativistic eft, Journal of High Energy Physics 2021, 1 (2021).
  * Akhiezer  _et al._ [1961] A. I. Akhiezer, V. G. Bar’yakhtar, and M. I. Kaganov, Spin waves in ferromagnets and antiferromagnets. I, Soviet Physics Uspekhi 3, 567 (1961).
  * Yung-Li and Callen [1964] W. Yung-Li and H. B. Callen, Spin waves in the spin-flop phase of an antiferromagnet, and metastability of the spin-flop transition, Journal of Physics and Chemistry of Solids 25, 1459 (1964).
  * Román and Soto [2000] J. M. Román and J. Soto, Spin waves in canted phases: An application to doped manganites, Phys. Rev. B 62, 3300 (2000).
  * Brauner and Hofmann [2017] T. Brauner and C. P. Hofmann, Two-loop free energy of three-dimensional antiferromagnets in external magnetic and staggered fields, Annals of Physics 386, 178 (2017).
  * Brauner and Jakobsen [2018] T. Brauner and M. F. Jakobsen, Scattering amplitudes of massive Nambu-Goldstone bosons, Phys. Rev. D 97, 025021 (2018).
  * Argurio  _et al._ [2015] R. Argurio, A. Marzolla, A. Mezzalira, and D. Naegels, Note on holographic nonrelativistic Goldstone bosons, Phys. Rev. D 92, 066009 (2015).
  * Alvarez-Gaume  _et al._ [2017] L. Alvarez-Gaume, O. Loukas, D. Orlando, and S. Reffert, Compensating strong coupling with large charge, Journal of High Energy Physics 2017, 1 (2017).
  * Alvarez-Gaume  _et al._ [2021] L. Alvarez-Gaume, D. Orlando, and S. Reffert, Selected topics in the large quantum number expansion, Physics Reports 933, 1 (2021).
  * Antipin  _et al._ [2020] O. Antipin, J. Bersini, F. Sannino, Z.-W. Wang, and C. Zhang, Charging the $\mathrm{O}(N)$ model, Phys. Rev. D 102, 045011 (2020).


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Simulating quantum circuit expectation values by Clifford perturbation
theory

Tomislav Begušić  Kasra Hejazi  Garnet Kin-Lic Chan<EMAIL_ADDRESS>Division of Chemistry and Chemical Engineering, California Institute of
Technology, Pasadena, California 91125, USA

###### Abstract

The classical simulation of quantum circuits is of central importance for
benchmarking near-term quantum devices. The fact that gates belonging to the
Clifford group can be simulated efficiently on classical computers has
motivated a range of methods that scale exponentially only in the number of
non-Clifford gates. Here, we consider the expectation value problem for
circuits composed of Clifford gates and non-Clifford Pauli rotations, and
introduce a heuristic perturbative approach based on the truncation of the
exponentially growing sum of Pauli terms in the Heisenberg picture. Numerical
results are shown on a Quantum Approximate Optimization Algorithm (QAOA)
benchmark for the E3LIN2 problem and we also demonstrate how this method can
be used to quantify coherent and incoherent errors of local observables in
Clifford circuits. Our results indicate that this systematically improvable
perturbative method offers a viable alternative to exact methods for
approximating expectation values of large near-Clifford circuits.

##  1 Introduction

Validating near-term quantum devices with classical simulations is paramount
for their future development [1]. In general, the classical simulation of
quantum circuits is hard and thus limited to a small number of qubits or
gates. One exception is Clifford circuits, which can be simulated efficiently
using the stabilizer formalism [2, 3, 4]. From this starting point, circuit
simulation methods have been developed which scale exponentially only with the
number of non-Clifford gates [5, 6, 7, 8, 9, 10]. Currently, the low-rank
stabilizer method proposed by Bravyi et al. [8] is one of the best-scaling
methods for sampling the output distribution of a quantum circuit with a small
number of non-Clifford gates.

In this work, we are concerned not with sampling the full output distribution
of a quantum circuit, but only with evaluating the expectation value (mean
value) of Pauli operators. Quantum mean values [11] are not only a central
aspect of the quantum output, but also form the cost function in variational
quantum algorithms [12], such as the variational quantum eigensolver (VQE)
[13] or QAOA [14, 15]. Here, we introduce a technique to obtain expectation
values of quantum circuits based on a Clifford-based perturbation method. The
method is most efficient when the gates are close to Clifford. Specifically,
in the proposed heuristic inspired by time-dependent perturbation theory, the
general Heisenberg evolution of a Pauli observable is computed as a polynomial
in non-Clifford parameters, and the polynomial is truncated to control the
cost. We use numerical examples to show that accurate results can be obtained
even with low orders of perturbation, and we demonstrate on a QAOA benchmark
that the method is faster than computing expectation values by sampling
methods by orders of magnitude. Finally, we show how the method can be used to
model coherent and incoherent noise in Clifford circuits, which is relevant to
benchmarking stabilizer error-correcting codes with many qubits.

##  2 Method

We consider the expectation value

$\langle O\rangle=\langle 0^{\otimes
n}|C_{N}^{{\dagger}}U_{N}^{{\dagger}}\cdots
C_{1}^{{\dagger}}U_{1}^{{\dagger}}OU_{1}C_{1}\cdots U_{N}C_{N}|0^{\otimes
n}\rangle$ (1)

of a Pauli operator $O$ evolved under a set of general unitary operators
(quantum gates) $U_{i}$ and Clifford gates $C_{i}$ in the Heisenberg picture.
For convenience, we label the gates in the order in which they are applied to
the operator. Note that since Clifford gates both include the identity and are
subsets of general unitary gates, the above expression can represent any
quantum circuit with $N$ gates.

Without loss of generality, through circuit compilation we can consider
circuits (1) containing only Pauli rotations $U_{i}\equiv
U_{i}(\theta_{i})=\exp(-i\theta_{i}P_{i}/2)$, where $P_{i}$ are Pauli
operators, and, further, that the angles satisfy $|\theta_{i}|\leq\pi/4$. The
latter is because for any $|\theta_{i}|>\pi/4$, we can use the angle
transformation $U_{i}(\theta_{i})=U_{i}(\tilde{\theta_{i}}+k\pi/2)C_{U_{i}}$,
where $C_{U_{i}}=\exp(-ik\pi P_{i}/4)$ ($k\in\mathbb{Z}$) is again a Clifford
gate. Similarly, because Clifford gates map the Pauli group onto itself [2],
every gate $C_{i}$ can be applied to the observable $O$ and all $P_{j<i}$
without increasing the number of terms, leaving only Pauli rotations with
modified $P_{i}$. Applying these two transformations is efficient and yields
the Clifford interaction picture, whereby Eq. (1) becomes

$\displaystyle\langle O\rangle=\langle 0^{\otimes
n}|{U}_{N}(\theta_{N})^{{\dagger}}\cdots{U}_{1}(\theta_{1})^{{\dagger}}O{U}_{1}(\theta_{1})\cdots{U}_{N}(\theta_{N})|0^{\otimes
n}\rangle$ (2)

with all $|\theta_{i}|\leq\pi/4$, and the $U_{i}$ have been transformed by the
Cliffords.

To reduce the cost of evaluating Eq. (2), we first note that many circuits
possess a certain structure that reduces the number of gates that affect the
final expectation value. Specifically, the application of a Pauli rotation to
a Pauli operator yields:

$e^{i\theta P/2}Oe^{-i\theta P/2}=\begin{cases}O,&[P,O]=0,\\\
\cos(\theta)O+i\sin(\theta)PO&\\{P,O\\}=0.\end{cases}$ (3)

Therefore, only gates that anticommute with the observable increase the number
of Pauli terms needed to represent the Heisenberg-evolved observable and thus
contribute to the computational cost. The full circuit is applied by iterating
Eq. (3), that is, at step $i+1$, the evolved observable is
$O_{i+1}=U_{i+1}(\theta_{i+1})^{{\dagger}}O_{i}U_{i+1}(\theta_{i+1})=\sum_{j}U_{i+1}(\theta_{i+1})^{{\dagger}}\sigma_{i,j}U_{i+1}(\theta_{i+1})$,
where $O_{i}=\sum_{j}\sigma_{i,j}$ and $\sigma_{i,j}$ are Pauli operators. The
worst-case scaling of this method is $2^{N}$, which is attained only if all
$P_{i}$ anticommute with $O$ and commute with each other. For random circuits,
we can expect that $P_{i+1}$ will commute on average with half of the Pauli
terms in the evolved observable and generate $M_{i+1}=3/2M_{i}$ Pauli terms,
leading to an average scaling of $(3/2)^{n}$.

In this work, we do not aim to propose an exact or rigorously
$\delta$-approximate method with an asymptotic exponential scaling better than
that already available in the literature. Rather, we are interested in
computing expectation values approximately, but in a systematically improvable
way, via a perturbation expansion. For this purpose, we note that the ratio
between coefficients of the $O$ and $PO$ branches in Eq. (3) is
$\tan(\theta)$. Since $|\theta|\leq\pi/4$, we know that $|\tan(\theta)|\leq
1$, i.e., the two branches will be weighted equally in the worst case, but
otherwise, the term involving the unmodified Pauli operator $O$ will be
larger. This motivates an ordering of terms by the number of multiplications
by $\sin(\theta_{i})$, which we call a perturbation order $k$. We can rewrite
the mean value (1) as

$\displaystyle\langle O\rangle$ $\displaystyle=\sum_{k=0}^{\tilde{N}}E^{(k)},$
(4) $\displaystyle E^{(k)}$ $\displaystyle=i^{k}\sum_{1\leq
j_{1}<j_{2}<\cdots<j_{k}\leq\tilde{N}}c_{j}\sin(\theta_{j_{1}})\cdots\sin(\theta_{j_{k}})\langle
0^{\otimes n}|P_{j_{k}}\cdots P_{j_{1}}O|0^{\otimes n}\rangle,$ (5)

where $\tilde{N}\leq N$ is the highest order of perturbation and $c_{j}$ are
products of cosine functions,

$c_{j}=\prod_{l\notin\\{j_{1},\dots,j_{k}\\}}^{N}\cos^{\xi_{l}}(\theta_{l}).$
(6)

In Eq. (6), $\xi_{l}=0$ if the Pauli operator $P_{l}$, associated with the
Pauli rotation gate $U_{l}(\theta_{l})$, commutes with $P_{j_{m}}\cdots
P_{j_{1}}O$ for $j_{m}<l<j_{m+1}$ and $\xi_{l}=1$ otherwise. The goal then is
to explore the truncated series

$\langle O\rangle^{(K)}=\sum_{k=0}^{K}E^{(k)},\quad K\leq\tilde{N}$ (7)

as an approximation for $\langle O\rangle$. This series expansion resembles
the series expansion of time-dependent perturbation theory, thus we refer to
this method as Clifford perturbation theory. The angles $|\tan(\theta)|$ serve
as a measure of the non-Clifford nature of the gates and are analogous to
perturbation parameters, and the truncated series is thus most accurate for
near-Clifford circuits where $|\tan(\theta)|\sim 0$. Note that because the
circuit is finite and discrete, unlike in perturbation theories associated
with continuous Hamiltonian evolution, the series always truncates at
$K=\tilde{N}$ and cannot diverge.

For simplicity, let us consider an example in which all Pauli operators
$P_{i}$ anticommute with $O$ and commute with each other. In this case,
$\tilde{N}=N$ and

$E^{(k)}=i^{k}\left(\prod_{i=1}^{N}\cos(\theta_{i})\right)\sum_{1\leq
j_{1}<j_{2}<\cdots<j_{k}\leq
N}\tan(\theta_{j_{1}})\cdots\tan(\theta_{j_{k}})\langle 0^{\otimes
n}|P_{j_{k}}\cdots P_{j_{1}}O|0^{\otimes n}\rangle.$ (8)

On one hand, we can rely on the fact that each term at order $k$ will have an
additional factor of $\tan(\theta_{j_{k}})$ compared to a term at order $k-1$.
On the other hand, there are $M_{k}=\binom{N}{k}$ terms at order $k$. In
practice, however, many terms in Eq. (8) might be zero regardless of the value
of $\theta$. For example, for a random Pauli operator $P$ with a weight $w$
the probability of $\langle 0^{\otimes n}|P|0^{\otimes n}\rangle\neq 0$ is
$3^{-w}$. Since this possibility is not known before the Pauli operator
$P_{j_{k}}\cdots P_{j_{1}}O$ in (8) is computed, we choose to study this
perturbative expansion numerically. In Appendix A, we present a more detailed
analysis of the perturbative treatment as applied to random quantum circuits.

We note that a work closely related to ours recently appeared, where the
authors consider the truncated Fourier series as an approximation to the
quantum mean value problem [16]. However, their expansion orders terms in
Fourier levels that are defined differently from the perturbation series
considered in this work. In addition, another perturbative method [17] was
proposed for approximately optimizing VQE parameters, in this case, the angles
of Pauli rotation gates, in a circuit composed of alternating Pauli and
Clifford gates. There, the authors expanded the mean value to second order in
the parameters, thereby effectively reducing the problem to a classically
efficient Clifford-circuit simulation of the expectation value, its gradient,
and Hessian with respect to the parameters.

##  3 Numerical examples

The Clifford perturbation method was implemented using Qiskit [18] and studied
on two different problems. First, we applied the method to evaluate the QAOA
cost function of the combinatorial Max E3LIN2 problem [19, 8]. Second, we
simulated mean values of a Clifford circuit that was subject to coherent
noise, a common source of errors in modern-day implementations of quantum
circuits [20].

###  3.1 QAOA applied to the Max E3LIN2 problem

The Max E3LIN2 problem is a combinatorial optimization task of finding
bitstrings $z_{1},z_{2},...,z_{n}$ that maximize the cost function

$C_{z}=\frac{1}{2}\sum_{1\leq u<v<w\leq n}d_{uvw}z_{u}z_{v}z_{w},$ (9)

determined by coefficients $d_{uvw}\in\\{0,\pm 1\\}$ and a parameter $D$ that
defines the number of occurrences of a given $z_{i}$ in the sum (9). QAOA
solves this task by constructing an ansatz

$|\mathbf{\beta},\mathbf{\gamma}\rangle=\prod_{j=1}^{p}e^{-\beta_{j}\sum_{i}^{n}X_{i}}e^{-i\gamma_{j}C}H^{\otimes
n}|0^{\otimes n}\rangle,$ (10)

and optimizing the expectation value $\langle
C\rangle=\langle\mathbf{\beta},\mathbf{\gamma}|C|\mathbf{\beta},\mathbf{\gamma}\rangle$,
where $C$ is obtained from $C_{z}$ by replacing the bits $z_{u}$ in Eq. (9) by
Pauli operators $Z_{u}$. In what follows, we evaluate the expectation value

$\langle C\rangle=\frac{1}{2}\sum_{1\leq u<v<w\leq
n}d_{uvw}\langle\mathbf{\beta},\mathbf{\gamma}|Z_{u}Z_{v}Z_{w}|\mathbf{\beta},\mathbf{\gamma}\rangle$
(11)

by applying the algorithm introduced earlier to individual terms on the right-
hand side of (11).

Figure 1: (a) QAOA cost function (11) for the E3LIN2 problem with $n=50$
qubits and $D=4$. Only positive values of $\gamma$ are shown because of the
symmetry relation $C(-\gamma)=-C(\gamma)$. (b) Histogram over 100 runs of the
maximum order $K$ with a nonzero contribution to the cost function (11) with
different values of $D$ and $n=50$ qubits.

Following Bravyi et al. [8], we consider an example with $n=50$ qubits, $p=1$
(one-layer ansatz), and $\beta=\pi/4$. The expectation value $\langle
C\rangle$ is shown as a function of $\gamma$ in Fig. 1 (top), where each point
took only a few seconds to compute on a laptop. In comparison, in [8],
generating the same data as in Fig. 1a with the low-rank stabilizer method was
reported to take about 3 days. This reflects the fact that whereas sampling
from the probability distribution $p(x)=|\langle x|\beta,\gamma\rangle|^{2}$
is generically hard, the expectation values (11) can be computed in polynomial
time using a Monte Carlo algorithm introduced by Van Den Nest [8, 21].
Therefore, this model illustrates the significant simplification afforded if
one targets the computation of expectation values, due to the light-cone
structure and the shallow depth. Indeed, even though the whole circuit is
composed of $N=\lfloor n\times D/3\rfloor=66$ non-Clifford gates, only a
fraction of those, namely $\tilde{N}\leq 3(D-1)+1$, anticommute with the
individual Pauli terms in the observable [22].

However, in Fig. 1a we also show that in this case there is no error if we
truncate the perturbation order to $K=1$, regardless of the value of $\gamma$.
This arises from the algebraic structure of the gates and observable, which
mean that not all operators in the observable lightcone yield a non-zero
expectation value, giving additional savings. Prompted by this result, we
evaluated randomly generated circuits for $1\leq D\leq 5$ and for each
instance evaluated the number of instances with contributions from higher
orders to the expectation value. As seen in Fig. 1b, almost all instances
require only the first-order ($K=1$) perturbation contribution, whereas at
higher $D$ there are examples with nonzero third- and fifth-order
contributions. These contributions will additionally be weighted by factors of
$\tan(\theta)$ which will usually make such terms small. Thus in certain
shallow QAOA circuits [23], the perturbation expansion can take advantage of
additional truncations and approximations to the observable lightcone,
speeding up the classical simulation.

###  3.2 Coherent error in Clifford circuits

Clifford gates are essential for the implementation of quantum error
correction [24] and for the validation of new quantum hardware [1]. Quantum
coherent errors introduce a bias in the rotation angles, which implies that
simulating such noisy circuits can become difficult because the gates are no
longer Clifford. Here, we show that our method can be of use in such cases,
especially in the limit of small errors in the rotation angles. Moreover, we
show in Appendix B that our approach is applicable to incoherent errors as
well, including those represented by Pauli and amplitude- or phase-damping
channels.

We construct a Clifford circuit

$U=\prod_{j=1}^{p}B_{j}A_{j}$ (12)

from alternating layers of one-qubit

$A=e^{-\frac{i}{2}\sum_{i=1}^{n}\theta_{i}\sigma_{i}}$ (13)

and two-qubit

$B=e^{-\frac{i}{2}\sum_{i=1}^{n/2}\theta_{i}P_{i}},\quad
P_{i}=\sigma_{i_{1}}^{\prime}\sigma_{i_{2}}^{\prime\prime}$ (14)

Pauli rotation gates, where
$\sigma,\sigma^{\prime},\sigma^{\prime\prime}\in\\{X,Y,Z\\}$ are selected
randomly. In (14), all $i_{k}$ in the pairs are distinct and chosen randomly,
i.e., exactly $n/2$ random disjoint pairs appear within a single two-qubit
layer. The rotation angles are set to $\theta\in\\{\pm\pi/2,\pm\pi\\}$, so
that all gates are Clifford. Then, all rotations are distorted
($\theta\rightarrow\theta+\delta\theta$) by $\delta\theta$, which represents
the error. We use $n=50$ qubits and measure the expectation value of
$Z_{1}Z_{26}$.

Figure 2: $\langle Z_{1}Z_{26}\rangle$ evaluated for a Clifford circuit
defined by Eqs. (12)–(14) with $n=50$ qubits, $p=4$, and a coherent error
defined by $\delta\theta$ (see text). (a) Expectation value as a function of
$\delta\theta$ for several perturbation orders $K$; (b, c) error in the
expectation value as a function of the perturbation order (on linear (b) and
logarithmic (c) scales); (d) cumulative (dashed, right-hand axis) and per-
order (solid, left-hand axis) number of Pauli operators generated in the
evaluation of $\langle Z_{1}Z_{26}\rangle$.

Figure 2a plots the expectation value as a function of $\delta\theta$ for
different values of perturbation order $K$. As $\delta\theta$ increases, the
expectation value under coherent noise diverges from its noiseless result.
Unsurprisingly, the low-order perturbation approach works better for small
$\delta\theta$, whereas for larger $\delta\theta$ higher-order perturbation
terms are needed. This is shown clearly in Figs. 2b, c. In this example with
$p=4$ ($300$ gates), we can simulate the expectation values accurately, e.g.,
with error $<10^{-2}$, already with $K=4$ or even $K=2$ for the smallest
values of $\delta\theta$. For comparison, the total number of Pauli operators
generated by evolving the observable as $U^{{\dagger}}Z_{1}Z_{26}U$ is
$135930$, whereas only $2313$ operators are generated for orders up to $K=4$
(see Fig. 2d).

Figure 3: Mean absolute values of errors in $\langle Z_{1}Z_{26}\rangle$ for
1000 random circuits generated according to Eqs. (12)–(14) with $n=50$ qubits,
$p=3$, and $\delta\theta=0.2$. The results are shown for two different
calculations: one in which the perturbation theory is applied directly to the
Pauli rotation gates (“no angle transform”) and the other in which the
rotation gates are first transformed so that $|\theta|<\pi/4$ (“angle
transform”).

Next, we illustrate the advantage of carrying out perturbation theory after
using the angle transformation on the circuit so that $|\theta|\leq\pi/4$, as
outlined in Sec. 2. Figure 3 compares the result from Eq. 2 with and without
the angle transformation. Here, the mean absolute error is computed from 1000
examples constructed with $n=50$ qubits and $p=3$ layers. Low perturbation
orders contribute little to the final expectation values if the rotation gates
are not transformed, whereas an opposite trend is observed within the angle
transformed picture.

Finally, let us consider an example where performing the exact computation
would require extensive computational resources, beyond what can be done with
a basic implementation on a personal computer. We set $n=100$, $p=3$ (450
gates), and $\delta\theta=0.2$. The calculation up to perturbation order $K=7$
took approximately 2 hours on a laptop and generated around $10^{7}$ Pauli
terms. Given that the expectation values for all $K\geq 4$ are similar, we can
estimate $\langle Z_{1}Z_{50}\rangle\approx\langle
Z_{1}Z_{50}\rangle^{(7)}=0.0023$.

Figure 4: $\langle Z_{1}Z_{50}\rangle$ expectation value (solid, left-hand
$y$-axis) for a Clifford circuit ($n=100$, $p=3$) with coherent error
($\delta\theta=0.2$) computed up to perturbation order $K=7$. The dashed line
(right-hand $y$-axis) shows the total number of Pauli operators generated up
to order $K$.

##  4 Conclusion

To conclude, we have introduced a Clifford-based perturbation theory that
enables the efficient computation of expectation values of near-Clifford
quantum circuits. We studied the merits of this approximate scheme in the
context of QAOA and Clifford circuits subject to coherent noise. These
numerical examples helped illustrate the two sources of speed-up in our
method: First, the fact that in some circuits, for example shallow QAOA, only
low-order terms in the perturbation expansion are found to be non-zero, and
second, the fact that sufficiently small angles of Pauli rotation gates
suppress contributions from higher-order terms. In conjunction, these
simplifications enable the practical determination of expectation values of
observables of large near-Clifford circuits, for example, with over a hundred
qubits and hundreds of gates, using only laptop computational resources.

##  5 Acknowledgments

We acknowledge informative discussions with Steve Flammia, Ali Lavasani, David
Gosset, Sergey Bravyi, and Alex Dalzell. TB and GKC were supported by the US
Department of Energy, Office of Science, Office of Advanced Scientific
Computing Research and Office of Basic Energy Sciences, Scientific Discovery
through Advanced Computing (SciDAC) program under Award Number DE-SC0022088.
TB acknowledges financial support from the Swiss National Science Foundation
through the Postdoc Mobility Fellowship (grant number P500PN-214214). GKC is a
Simons Investigator in Physics.

##  Appendix A Random quantum circuits

Let us consider a random quantum circuit composed of Pauli rotations with a
fixed small angle $\theta_{i}=\theta$, such that $\cos(\theta)\approx 1$.
Then,

$E^{(k)}=i^{k}\sin^{k}(\theta)\sum_{1\leq
j_{1}<\cdots<j_{k}\leq\tilde{N}}\langle 0^{\otimes n}|P_{j_{k}}\cdots
P_{j_{1}}O|0^{\otimes n}\rangle.$ (15)

In the remainder, we will not discuss the value of $i^{k}\langle 0^{\otimes
n}|P_{j_{k}}\cdots P_{j_{1}}O|0^{\otimes n}\rangle\in\\{0,\pm 1\\}$, but
rather focus on the average value of the sum of absolute values of the
coefficients

$|E|^{(k)}=M^{(k)}|\sin(\theta)|^{k},$ (16)

where $M^{(k)}$ denotes the average number of Pauli operators at perturbation
order $k$. We can deduce $M^{(k)}$ by assuming that the Pauli rotation gate
$U_{i+1}(\theta)$ commutes with each term in the evolved observable $O_{i}$
with a probability of 0.5. Let us label one Pauli operator in $O_{i}$ as
$o_{i}^{(k)}$, where $k$ indicates the number of times it has been multiplied
by $\sin(\theta)$, i.e., its order. If $[P_{i+1},o_{i}^{(k)}]=0$, according to
Eq. (3), $U_{i+1}(\theta)^{{\dagger}}o_{i}^{(k)}U_{i+1}(\theta)=o_{i}^{(k)}$,
i.e., the number of terms remains the same. Otherwise, if
$\\{P_{i+1},o_{i}^{(k)}\\}=0$,
$U_{i+1}(\theta)^{{\dagger}}o_{i}^{(k)}U_{i+1}(\theta)=\cos(\theta)o_{i}^{(k)}+i\sin(\theta)P_{i+1}o_{i}^{(k)}$.
Here, the number of Pauli operators of order $k$ remains the same, but the
number of terms of order $k+1$ is increased by 1. This results in the
following recursive formula

$M_{i}^{(k)}=M_{i-1}^{(k)}+\frac{1}{2}M_{i-1}^{(k-1)},$ (17)

where $M_{i}^{(k)}$ denotes the number of terms of order $k$ in the observable
$O_{i}$ evolved up to step $i$ and $M_{0}^{(k)}=\delta_{k0}$. By induction, we
can prove the following closed-form expression

$M^{(k)}=M_{N}^{(k)}=2^{-k}\binom{N}{k},$ (18)

which agrees with the average total number of Pauli operators
$M=\sum_{k=0}^{N}M^{(k)}=(3/2)^{N}$. Then,

$|E|^{(k)}=|\sin(\theta)/2|^{k}\binom{N}{k}$ (19)

and

$|O|^{(K)}=\sum_{k=0}^{K}|E|^{(k)}=\sum_{k=0}^{K}|\sin(\theta)/2|^{k}\binom{N}{k}.$
(20)

In particular,

$|O|=|O|^{(N)}=(1+|\sin(\theta)/2|)^{N}.$ (21)

Let us now consider an upper bound on the relative error

$\frac{|O|-|O|^{(K)}}{|O|}<\delta.$ (22)

Specifically, we wish to find the smallest value of perturbation order $K$
that satisfies the equation above for given $N$, $\theta$, and $\delta$. As an
example, Fig. 5 shows $K$ as a function of $N$ for $\delta=0.01,0.05$, and
$\theta=0.2$.

Figure 5: Minimal value of the perturbation order $K$ that satisfies Eq. (22)
for different numbers of gates $N$, angle $\theta=0.2$, and $\delta=0.05$ (red
markers) or $\delta=0.01$ (blue markers). Solid lines (right-hand $y$-axis)
correspond to $\sum_{k=0}^{K}M^{(k)}$ for each pair of $N$ and $K$.

##  Appendix B A note on incoherent errors

Several basic models of errors [25] can be represented as a Pauli channel

$\tilde{\rho}=(1-s_{x}-s_{y}-s_{z})\rho+s_{x}X\rho X+s_{y}Y\rho Y+s_{z}Z\rho
Z,$ (23)

including bit flip errors

$\tilde{\rho}=(1-s)\rho+sX\rho X,$ (24)

phase flip errors

$\tilde{\rho}=(1-s)\rho+sZ\rho Z,$ (25)

and a depolarizing channel

$\tilde{\rho}=(1-s)\rho+\frac{s}{3}(X\rho X+Y\rho Y+Z\rho Z).$ (26)

Pauli channels can be implemented with our method deterministically and at
almost no additional computational cost. First, we note that the Pauli
operators and the scaling factor in Eq. (23) can be applied directly to the
observable (i.e. in the Heisenberg picture) instead of on the density
operator. Second, the conjugation of one Pauli operator by another results
only in a potential sign change:

$\sigma_{\text{Error}}\sigma_{\text{Observable}}\sigma_{\text{Error}}=\pm\sigma_{\text{Observable}},$
(27)

where the sign depends on whether the two Pauli operators commute. Therefore,
a one-qubit Pauli channel applied to a Pauli operator will only introduce a
scaling factor

$\eta=1-s_{x}-s_{y}-s_{z}\pm s_{x}\pm s_{y}\pm s_{z}.$ (28)

Another example is the amplitude-damping error [25], which can be represented
as:

$\tilde{\rho}=E_{0}^{\dagger}\rho E_{0}+E_{1}^{\dagger}\rho E_{1},$ (29)

where

$\displaystyle E_{0}$ $\displaystyle=\begin{pmatrix}1&0\\\
0&\sqrt{1-\lambda}\end{pmatrix},$ (30) $\displaystyle E_{1}$
$\displaystyle=\begin{pmatrix}0&\sqrt{\lambda}\\\ 0&0\end{pmatrix}.$ (31)

Equation (29) applied to single-qubit Pauli operators yields:

$\displaystyle\tilde{I}$ $\displaystyle=I$ (32) $\displaystyle\tilde{X}$
$\displaystyle=\sqrt{1-\lambda}X$ (33) $\displaystyle\tilde{Y}$
$\displaystyle=\sqrt{1-\lambda}Y$ (34) $\displaystyle\tilde{Z}$
$\displaystyle=(1-\lambda)Z+\lambda I.$ (35)

$I$, $X$, and $Y$ can be efficiently simulated, but any occurrence of $Z$ will
produce twice as many Pauli terms. Fortunately, the two branches in the last
equation have different weights if $\lambda$ is small, so the problem can be
again treated perturbatively by keeping track of powers of $\lambda$.

Finally, the phase-damping channel can be simulated efficiently. It is similar
to the amplitude damping channel [Eq. (29)] but with

$\displaystyle E_{0}$ $\displaystyle=\begin{pmatrix}1&0\\\
0&\sqrt{1-\lambda}\end{pmatrix},$ (36) $\displaystyle E_{1}$
$\displaystyle=\begin{pmatrix}0&0\\\ 0&\sqrt{\lambda}\end{pmatrix}.$ (37)

The application to $Z$ now leaves it unchanged ($\tilde{Z}=Z$), while the
action on all other operators is the same as in the amplitude-damping channel
[Eqs. (32)–(34)]. Therefore, the phase-damping channel can be simulated at no
additional cost.

## References

  * [1] Andrew W. Cross, Lev S. Bishop, Sarah Sheldon, Paul D. Nation, and Jay M. Gambetta.  “Validating quantum computers using randomized model circuits”.  Phys. Rev. A 100, 032328 (2019).
  * [2] Daniel Gottesman.  “The Heisenberg Representation of Quantum Computers” (1998).  arXiv:9807006.
  * [3] Scott Aaronson and Daniel Gottesman.  “Improved simulation of stabilizer circuits”.  Phys. Rev. A 70, 052328 (2004).
  * [4] Maarten Van Den Nes.  “Classical simulation of quantum computation, the gottesman-knill theorem, and slightly beyond”.  Quantum Info. Comput. 10, 258–271 (2010).
  * [5] Sergey Bravyi and David Gosset.  “Improved Classical Simulation of Quantum Circuits Dominated by Clifford Gates”.  Phys. Rev. Lett. 116, 250501 (2016).
  * [6] Ryan S. Bennink, Erik M. Ferragut, Travis S. Humble, Jason A. Laska, James J. Nutaro, Mark G. Pleszkoch, and Raphael C. Pooser.  “Unbiased simulation of near-Clifford quantum circuits”.  Phys. Rev. A 95, 062337 (2017).
  * [7] Hammam Qassim, Joel J. Wallman, and Joseph Emerson.  “Clifford recompilation for faster classical simulation of quantum circuits”.  Quantum 3, 170 (2019).
  * [8] Sergey Bravyi, Dan Browne, Padraic Calpin, Earl Campbell, David Gosset, and Mark Howard.  “Simulation of quantum circuits by low-rank stabilizer decompositions”.  Quantum 3, 181 (2019).
  * [9] Yifei Huang and Peter Love.  “Feynman-path-type simulation using stabilizer projector decomposition of unitaries”.  Phys. Rev. A 103, 022428 (2021).
  * [10] Aleks Kissinger and John van de Wetering.  “Simulating quantum circuits with zx-calculus reduced stabiliser decompositions”.  Quantum Sci. Technol. 7, 044001 (2022).
  * [11] Sergey Bravyi, David Gosset, and Ramis Movassagh.  “Classical algorithms for quantum mean values”.  Nat. Phys. 17, 337–341 (2021).
  * [12] M. Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C. Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R. McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, and Patrick J. Coles.  “Variational quantum algorithms”.  Nat. Rev. Phys. 3, 625–644 (2021).
  * [13] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O’Brien.  “A variational eigenvalue solver on a photonic quantum processor”.  Nat. Commun. 5, 4213 (2014).
  * [14] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  “A quantum approximate optimization algorithm” (2014).  arXiv:1411.4028.
  * [15] Edward Farhi and Aram W Harrow.  “Quantum Supremacy through the Quantum Approximate Optimization Algorithm” (2016).  arXiv:1602.07674.
  * [16] Nikita A. Nemkov, Evgeniy O. Kiktenko, and Aleksey K. Fedorov.  “Fourier expansion in variational quantum algorithms” (2023).  arXiv:2304.03787.
  * [17] Kosuke Mitarai, Yasunari Suzuki, Wataru Mizukami, Yuya O. Nakagawa, and Keisuke Fujii.  “Quadratic Clifford expansion for efficient benchmarking and initialization of variational quantum algorithms”.  Phys. Rev. Res. 4, 033012 (2022).
  * [18] Qiskit contributors.  “Qiskit: An open-source framework for quantum computing” (2023).
  * [19] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  “A Quantum Approximate Optimization Algorithm Applied to a Bounded Occurrence Constraint Problem” (2015).  arXiv:1412.6062.
  * [20] Yingkai Ouyang.  “Avoiding coherent errors with rotated concatenated stabilizer codes”.  npj Quantum Inf. 7, 87 (2021).
  * [21] Maarten Van Den Nest.  “Simulating quantum computers with probabilistic methods”.  Quantum Info. Comput. 11, 784–812 (2011).
  * [22] Patrick Rall, Daniel Liang, Jeremy Cook, and William Kretschmer.  “Simulation of qubit quantum circuits via Pauli propagation”.  Phys. Rev. A 99, 062337 (2019).  arXiv:1901.09070.
  * [23] Wei-Feng Zhuang, Ya-Nan Pu, Hong-Ze Xu, Xudan Chai, Yanwu Gu, Yunheng Ma, Shahid Qamar, Chen Qian, Peng Qian, Xiao Xiao, Meng-Jun Hu, and Dong E. Liu.  “Efficient Classical Computation of Quantum Mean Values for Shallow QAOA Circuits” (2021).  arXiv:2112.11151.
  * [24] Sergey Bravyi, Joseph A. Latone, and Dmitri Maslov.  “6-qubit optimal Clifford circuits”.  npj Quantum Inf. 8, 79 (2022).
  * [25] Michael A. Nielsen and Isaac L. Chuang.  “Quantum computation and quantum information”.  Cambridge University Press. Cambridge (2010).