# Spin Statistics and Field Equations for any Spin

Feng Zixuan

###### Abstract

In this article we prove spin statistics theorem for arbitrary massive $(A,B)$
field in a representation theoretic manner. General Gamma matrices are
introduced, and explicit forms for low spin are calculated. Spin sums and
twisted spin sums are introduced to prove spin statistics and derive field
equations respectively, and their relations are discussed. Klein Gordan
equation is just the condition on the 4-momentum. Dirac equation (or massive
Weyl equations) and Proca equation are shown to be examples of our general
field equation.

###### Contents

  1. 1 Introduction
  2. 2 Gamma matrices
    1. 2.1 General spin
    2. 2.2 Examples
  3. 3 Spin statistics for massive particles
  4. 4 Field equations for massive particles
    1. 4.1 General equation
    2. 4.2 Examples
  5. 5 Conclusions

Keywords: Quantum field theory; Spin statistics; Gamma matrices; Dirac
equation; Proca equation;

##  1 Introduction

In this article we prove the spin statistics theorem in a representation
theoretic manner for arbitrary field spin $(A,B)$ for massive field. Also we
derive general field equation for massive field, and shows that the Dirac
equation and Proca equation are examples for it.

The logic chain of this article follows Weinberg’s path ([Wei64a]) of deriving
quantum field theory. We first have a unitary representation of Poincare group
on the state space (Fock space), then we define creation and annilation
operators explicitly, and then use them to construct free quantum fields. In
this approach the field equation does not serve as an equation that
characterizes the field uniquely, but is a property that the field possesses.
So field equation is not of great importance from certain point of view.

$\text{Fock space representation}\longrightarrow\text{Creation and
annilation}\longrightarrow\text{Free quantum field}\longrightarrow\text{Field
equation}$ (1)

First we recall the construction for free quantum field. We define
annihilation field $\psi_{l}^{+}(x)$ and creation field $\psi_{l}^{-}(x)$:

$\displaystyle\psi_{l}^{+}(x)=\sum_{\sigma n}\int
d^{3}pu_{l}(x;\mathbf{p},\sigma,n)a(\mathbf{p},\sigma,n),$ (2)
$\displaystyle\psi_{l}^{-}(x)=\sum_{\sigma n}\int
d^{3}pv_{l}(x;\mathbf{p},\sigma,n)a^{\dagger}(\mathbf{p},\sigma,n)$

We hope it to satisfy the Lorentz transformation of field: 111Our convention
is different from that of [Wei02].

######  Axiom 1.

(Lorentz transformation of field)

We have a finite dimensional representation $D$ of the universal cover
$SL(2,\mathbb{C})$ of the Lorentz group such that the ’Lorentz transformation
of the field’ are satisfied:

$\displaystyle
U_{0}(\Lambda,a)^{-1}\psi_{l}^{+}(\lambda(\Lambda)x+a)U_{0}(\Lambda,a)=D_{l\bar{l}}(\Lambda)\psi_{\bar{l}}^{+}(x)$
(3) $\displaystyle
U_{0}(\Lambda,a)^{-1}\psi_{l}^{-}(\lambda(\Lambda)x+a)U_{0}(\Lambda,a)=D_{l\bar{l}}(\Lambda)\psi_{\bar{l}}^{-}(x)$
$\displaystyle\forall\Lambda\in SL(2,\mathbb{C})$

In order for Axiom 1 to satisfy, it is equivalent that

$\displaystyle u(x;\mathbf{q})=(2\pi)^{-3/2}e^{ip\cdot x}u(\mathbf{q})$ (4)
$\displaystyle v(x;\mathbf{q})=(2\pi)^{-3/2}e^{-ip\cdot x}v(\mathbf{q})$

and:

$\displaystyle u(\mathbf{q})=(m/q^{0})^{1/2}D(L(q))u(\mathbf{0})$ (5)
$\displaystyle v(\mathbf{q})=(m/q^{0})^{1/2}D(L(q))v(\mathbf{0})$

and:

$\displaystyle u(\mathbf{0})D^{(j_{n})}(W)=D(W)u(\mathbf{0})$ (6)
$\displaystyle v(\mathbf{0})D^{(j_{n})^{*}}(W)=D(W)v(\mathbf{0})$

$\forall W\in\mathcal{W}$, where we have used the convention:
$u(q)=u(\mathbf{q})$ denotes the matrix
$\\{u_{l}(\mathbf{q},\sigma)\\}_{l|\sigma}$, and the multiplication is the
matrix multiplication.

Considering a irreducible field representation of spin $(A,B)$, we write

$\displaystyle\psi^{AB}_{ab}(x)=$ $\displaystyle(2\pi)^{-3/2}\sum_{\sigma}\int
d^{3}p[e^{ip\cdot
x}\kappa^{AB}u^{AB}_{ab}(\mathbf{p},\sigma)a(\mathbf{p},\sigma)+e^{-ip\cdot
x}\lambda^{AB}v^{AB}_{ab}(\mathbf{p},\sigma)a^{c\dagger}(\mathbf{p},\sigma)]$
(7)

with $\kappa$ and $\lambda$ arbitrary constants waiting for determination.

The next property we want the quantum field to possess is the microscopic
causality:

######  Axiom 2.

(Microscopic causality)

We should have

$[\psi^{AB}_{ab}(x),\psi^{CD\dagger}_{cd}(y)]_{\mp}=0,\quad(x-y)^{2}>0$ (8)

This requirement leads to the so called ’spin statistics theorem’, which is
discussed in Section 3. The central object here is the so called ’spin sum’:

$(2p^{0})^{-1}\pi^{AB,CD}(\mathbf{p})\equiv
u^{AB}(\mathbf{p})u^{CD\dagger}(\mathbf{p})=v^{AB}(\mathbf{p})v^{CD\dagger}(\mathbf{p})$
(9)

and we can show that it may be expressed as a covariant polynomial of
$p^{\mu}$ using generalized gamma matrices, which we will discuss in Section
2. As a result we can fix a normalization of $\kappa,\lambda$:

$\displaystyle\psi^{AB}_{ab}(x)=$ $\displaystyle(2\pi)^{-3/2}\sum_{\sigma}\int
d^{3}p[e^{ip\cdot
x}u^{AB}_{ab}(\mathbf{p},\sigma)a(\mathbf{p},\sigma)+e^{-ip\cdot
x}(-)^{2B}v^{AB}_{ab}(\mathbf{p},\sigma)a^{c\dagger}(\mathbf{p},\sigma)]$ (10)

where the factor $(-)^{2B}$ is essential for both spin statistics and the
existence for a field equation where we want to combine the two equations for
positive frequency part and the negative frequency part, as discussed in
Section 4. In writing down a field equation, we use a modified version of spin
sum, called ’twisted spin sum’ in this article.

##  2 Gamma matrices

###  2.1 General spin

We can prove that for any half integers $A,B,C,D,K$ satisfying

$\max\\{|A-D|,|B-C|\\}\leq K\leq\min\\{A+D,B+C\\}$ (11)

with $2K,2A+2D,2B+2C$ having the same parity, 222This parity condition will be
implicit when we write down the condition (11) afterwards.  there exists a set
of non-zero $(2A+1)(2B+1)\times(2C+1)(2D+1)$ matrices

${}^{ABCD}_{\quad
K}T^{\mu_{1}\mu_{2}\cdots\mu_{2K}},\quad\mu_{1},\mu_{2},\cdots,\mu_{2K}=0,1,2,3$
(12)

with the properties:

  * •

$T$ is symmetric in $\mu$’s.

  * •

$T$ is traceless in $\mu$’s, i.e.,

$g_{\mu_{1}\mu_{2}}T^{\mu_{1}\mu_{2}\cdots\mu_{2K}}=0$ (13)

  * •

$T$ is a tensor in $\mu$’s, i.e.,

$D^{AB}(\Lambda)T^{\mu_{1}\mu_{2}\cdots\mu_{2K}}D^{CD}(\Lambda)^{\dagger}=\Lambda_{\nu_{1}}{}^{\mu_{1}}\Lambda_{\nu_{2}}{}^{\mu_{2}}\cdots\Lambda_{\nu_{2K}}{}^{\mu_{2K}}T^{\nu_{1}\nu_{2}\cdots\nu_{2K}}$
(14)

where the ordinary matrix multiplication is understood on the LHS.

Consider the vector space $V$ consisting of all complex
$(2A+1)(2B+1)\times(2C+1)(2D+1)$-matrices. They furnish a representation of
$SL(2,\mathbb{C})$ by

$M\mapsto D^{AB}(\Lambda)MD^{CD\dagger}(\Lambda)$ (15)

One can see that this is isomorphic to the $(A,B)\otimes(D,C)$ representation,
which is because the complex conjugate of $D^{CD}$ is isomorphic (but not
directly equal) to $D^{DC}$. And we have 333From now on we use the convention.

$\displaystyle(A,B)\otimes(D,C)=(A\otimes D,B\otimes C)$
$\displaystyle=(\bigoplus_{|A-D|\leq K_{1}\leq A+D}K_{1},\bigoplus_{|B-C|\leq
K_{2}\leq B+C}K_{2})$ (16) $\displaystyle=\bigoplus_{|A-D|\leq K_{1}\leq
A+D,|B-C|\leq K_{2}\leq B+C}(K_{1},K_{2})$

So it contains a $(K,K)$-subrep. By the fact 444See [Wei64a] that
$(K,K)$-representation of $SL(2,\mathbb{C})$ consists of all symmetric
traceless tensors of rank $2K$, so we define $T$’s to be the standard basis in
this description and we are done.

###  2.2 Examples

Pauli matrices $\sigma^{\mu}$ are the $T$ matrices for
$(\frac{1}{2},0,\frac{1}{2},0)$, i.e. $\sigma^{\mu}=^{1/2,0,1/2,0}_{\quad
1/2}T^{\mu}$, where

$\sigma^{0}=I=\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\quad\sigma^{1}=X=\begin{pmatrix}0&1\\\
1&0\end{pmatrix}\quad\sigma^{2}=Y=\begin{pmatrix}0&-i\\\
i&0\end{pmatrix}\quad\sigma^{3}=Z=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}$
(17)

Indeed we see in the decomposition
$(\frac{1}{2},0)\otimes(0,\frac{1}{2})=(\frac{1}{2},\frac{1}{2})$ there is no
other subreps, so we only need to show they obey the property

$D^{1/2,0}(\Lambda)\sigma^{\mu}D^{1/2,0\dagger}(\Lambda)=\Lambda_{\nu}{}^{\mu}\sigma^{\nu}$
(18)

or its Lie algebra level (where $J^{\mu\nu}$ denotes the image of Lie algebra
basis under representation):

$i[J^{\mu\nu},\sigma^{\rho}]=-g^{\rho\mu}\sigma^{\nu}+g^{\rho\nu}\sigma^{\mu}$
(19)

This is by the definition of $(\frac{1}{2},0)$ representation:

$\displaystyle J_{1}$ $\displaystyle\mapsto\frac{1}{2}X$ (20) $\displaystyle
J_{2}$ $\displaystyle\mapsto\frac{1}{2}Y$ $\displaystyle J_{3}$
$\displaystyle\mapsto\frac{1}{2}Z$ $\displaystyle K_{1}$
$\displaystyle\mapsto-\frac{i}{2}X$ $\displaystyle K_{2}$
$\displaystyle\mapsto-\frac{i}{2}Y$ $\displaystyle K_{3}$
$\displaystyle\mapsto-\frac{i}{2}Z$

Similarly we can show that $\bar{\sigma}^{\mu}\equiv(I,-\mathbf{\sigma})$ are
the $T$ matrices for $(0,\frac{1}{2},0,\frac{1}{2})$, i.e.
$\bar{\sigma}^{\mu}=^{0,1/2,0,1/2}_{\quad 1/2}T^{\mu}$.

##  3 Spin statistics for massive particles

The results in this subsection apply only when the field representation is
irreducible and the particle is massive.

Let’s compute $[\psi^{AB}_{ab}(x),\psi^{CD\dagger}_{cd}(y)]_{\mp}$. Firstly we
have the equality

$\displaystyle{}_{\mp}=$ $\displaystyle(2\pi)^{-3/2}\int
d^{3}p(2p^{0})^{-1}\pi^{AB,CD}_{ab,cd}(\mathbf{p})$ (21)
$\displaystyle\times[\kappa^{AB}\kappa^{CD*}e^{ip\cdot(x-y)}\mp\lambda^{AB}\lambda^{CD*}e^{-ip\cdot(x-y)}]$

where

######  Definition 3.1.

(Spin sum)

Define the spin sum as

$(2p^{0})^{-1}\pi^{AB,CD}(\mathbf{p})\equiv
u^{AB}(\mathbf{p})u^{CD\dagger}(\mathbf{p})=v^{AB}(\mathbf{p})v^{CD\dagger}(\mathbf{p})$
(22)

It is a $(2A+1)(2B+1)\times(2C+1)(2D+1)$-matrix.

######  Remark 3.1.

We use the convention that $u^{AB}(p)$ is a matrix
$(u_{ab}^{AB}(p,\sigma))_{ab,\sigma}$, having $(2A+1)(2B+1)$ rows and $2j+1$
columns. The above is understood as a matrix multiplication. Written
explicitly in components, it reads

$(2p^{0})^{-1}\pi^{AB,CD}_{ab,cd}(\mathbf{p})\equiv\sum_{\sigma}u^{AB}_{ab}(\mathbf{p},\sigma)u^{CD*}_{cd}(\mathbf{p},\sigma)=\sum_{\sigma}v^{AB}_{ab}(\mathbf{p},\sigma)v^{CD*}_{cd}(\mathbf{p},\sigma)$
(23)

But for convenience, we adapt the more concise notation involving matrix
multiplication.

We immediately see its Lorentz transformation property:

######  Proposition 3.1.

$\pi^{AB,CD}(\Lambda
p)=D^{AB}(\Lambda)\pi^{AB,CD}(p)D^{CD}(\Lambda)^{\dagger}$ (24)

###### Proof.

$\displaystyle D^{AB}(\Lambda)\pi^{AB,CD}(p)D^{CD\dagger}(\Lambda)$
$\displaystyle=2p^{0}D^{AB}(\Lambda)u^{AB}(p)u^{CD\dagger}(p)D^{CD\dagger}(\Lambda)$
(25) $\displaystyle=2(\Lambda p)^{0}u^{AB}(\Lambda
p)D^{j}(W(\Lambda,p))D^{j\dagger}(W(\Lambda,p))u^{CD\dagger}(\Lambda p)$
$\displaystyle=2(\Lambda p)^{0}u^{AB}(\Lambda p)u^{CD\dagger}(\Lambda p)$
$\displaystyle=\pi^{AB,CD}(\Lambda p)$

∎

The major part of the proof of spin statistics is that we can express the spin
sum as a polynomial with a good parity property, and it is this polynomial’s
parity that determine the boson/fermion of the particle. The parity is used in
absorbing the minus sign in $e^{-ip\cdot x}$ of $\psi^{(-)}$. 555Until not we
do not assume the particle to be massive.

######  Theorem 3.1.

666Weinberg in [Wei02] proved the special case when $\mathbf{p}=(0,0,p^{3})$.
But it does not directly imply the general case.

In the massive case there is a 4-variables polynomial $P$ of coefficients
$(2A+1)(2B+1)\times(2C+1)(2D+1)$-matrices, such that on the mass shell its
value always equal the spin sum

$\pi^{AB,CD}(\mathbf{p},\sqrt{\mathbf{p}^{2}+m^{2}})=P(\mathbf{p},\sqrt{\mathbf{p}^{2}+m^{2}})$
(26)

and it is either an odd function or an even function determined by the parity
of $2A+2D$:

$P(-\mathbf{p},-p^{0})=(-1)^{2A+2D}P(\mathbf{p},p^{0})$ (27)

for all 4-vectors $p$.

Let’s prove the theorem.

First we can prove that the initial value $\pi^{AB,CD}(\mathbf{0})$ can be
expressed as a linear combination of above matrices

$\pi^{AB,CD}(\mathbf{0})=\sum_{K\text{satisfying
}\eqref{ABCDK}}\xi_{K}^{ABCD}{}^{ABCD}_{\quad K}T^{0_{1}0_{2}\cdots 0_{2K}}$
(28)

This is because by (24), $\pi^{AB,CD}(\mathbf{0})$ is invariant under the
representation restricted to $SU(2)$, so it must lie in the direct sum of
$0$-subrepresentations of $SU(2)$ when we decompose

$\displaystyle(A,B)\otimes(D,C)|_{SU(2)}$ $\displaystyle=\bigoplus_{|A-D|\leq
K_{1}\leq A+D,|B-C|\leq K_{2}\leq B+C}(K_{1},K_{2})|_{SU(2)}$ (29)
$\displaystyle=\bigoplus_{|A-D|\leq K_{1}\leq A+D,|B-C|\leq K_{2}\leq
B+C}\bigoplus_{|K_{1}-K_{2}|\leq K_{3}\leq K_{1}+K_{2}}K_{3}$

We see that each $0$-subrep of $SU(2)$ comes from a $(K,K)$-subrep of
$SL(2,\mathbb{C})$ where $K$ satisfies (11). Also, by (14), $T^{00\cdots 0}$
is non-zero and invariant under $SU(2)$, so they form a basis for the
$0$-subspace of $SU(2)$.

Finally, combine (24), (14) and (28), we have the equality on mass shell

$\pi^{AB,CD}(\mathbf{p})=\sum_{K\text{satisfying
}\eqref{ABCDK}}\xi_{K}^{ABCD}{}^{ABCD}_{\quad
K}T^{\mu_{1}\mu_{2}\cdots\mu_{2K}}m^{-2K}p_{\mu_{1}}p_{\mu_{2}}\cdots
p_{\mu_{2K}}$ (30)

So we define the polynomial as the RHS.

The proof is done.

With this theorem we can finally begin the proof for spin statistics. First we
show that we need only verify a simpler case:

######  Lemma 3.1.

777Weinberg in [Wei02] mentioned and used this result.

For $[\psi^{AB}_{ab}(x),\psi^{CD\dagger}_{cd}(y)]_{\mp}$ vanishes for space-
like $x-y$, it is sufficient that it vanishes when $x^{0}=y^{0}$.

###### Proof.

We can choose a Lorentz transformation $\Lambda$ such that $(\Lambda
x)^{0}=(\Lambda y)^{0}$. Under Lorentz transformation we have

$[\psi^{AB}(x),\psi^{CD\dagger}(y)]_{\mp}=U_{0}(\Lambda)^{-1}D^{AB}(\Lambda^{-1})[\psi^{AB}(\Lambda
x),\psi^{CD\dagger}(\Lambda
y)]_{\mp}D^{CD\dagger}(\Lambda^{-1})U_{0}(\Lambda)$ (31)

So it follows. ∎

Now we rewrite the above polynomial: substitute all exponentials of $p^{0}$ of
degree larger than 1 by substituting $(p^{0})^{2}=\mathbf{p}^{2}+m^{2}$. Then
the value on mass shell do not change (but the polynomial changes) and can be
written as

$\pi^{AB,CD}(\mathbf{p})=P(\mathbf{p})+2\sqrt{\mathbf{p}^{2}+m^{2}}Q(\mathbf{p})$
(32)

where $P$ and $Q$ are 3-variables polynomials in $\mathbf{p}$ alone, with

$\displaystyle P(-\mathbf{p})=(-)^{2A+2D}P(\mathbf{p})$ (33) $\displaystyle
Q(-\mathbf{p})=-(-)^{2A+2D}Q(\mathbf{p})$

For $x-y$ space-like, by the above lemma we can adopt a Lorentz frame in which
$x^{0}=y^{0}$, and thus

$\displaystyle{}_{\mp}=$
$\displaystyle[\kappa^{AB}\kappa^{CD*}\mp(-)^{2A+2D}\lambda^{AB}\lambda^{CD*}]P(-i\nabla)\Delta_{+}(\mathbf{x}-\mathbf{y},0)$
(34)
$\displaystyle+[\kappa^{AB}\kappa^{CD}\pm(-)^{2A+2D}\lambda^{AB}\lambda^{CD*}]Q(-i\nabla)\delta^{3}(\mathbf{x}-\mathbf{y})$

In order that this should vanish when $\mathbf{x}\neq\mathbf{y}$, we must have

$\kappa^{AB}\kappa^{CD*}=\pm(-1)^{2A+2D}\lambda^{AB}\lambda^{CD*}$ (35)

For the case where $A=C$ and $B=D$.

$|\kappa^{AB}|^{2}=\pm(-)^{2A+2B}|\lambda^{AB}|^{2}$ (36)

This is possible if and only if

$\pm(-1)^{2A+2B}=+1$ (37)

and

$|\kappa^{AB}|^{2}=|\lambda^{AB}|^{2}$ (38)

Return to the general case we have

$\frac{\kappa^{AB}}{\kappa^{CD}}=(-1)^{2A-2C}\frac{\lambda^{AB}}{\lambda^{CD}}$
(39)

To summarize what we have proved above:

######  Theorem 3.2.

(Spin statistics, massive irreducible case) 888This result is stated in
[Wei02], but there is no rigorous proof. 999The factor $(-)^{2B}$ is
irrelevant when we just consider one field, but is essential when considering
different fields.

In order for all fields $\psi^{AB}(x)$ constructed from a massive particle
with spin $j$ to satisfy Axiom 1, it is equivalent that whether it is boson or
fermion depends on whether $2j$ is even or odd, and
$\lambda^{AB}=(-)^{2B}\kappa^{AB}c$, where $c$ is a constant independent of
$A,B$.

######  Example 3.1.

Massive scalar field has $(A,B)=(0,0)$ and $j=0$, so it describes bosons.

Massive vector field has $(A,B)=(\frac{1}{2},\frac{1}{2})$ and $j=0,1$, so it
describes bosons.

Massive Weyl field has $(\frac{1}{2},0)$ or $(0,\frac{1}{2})$ and
$j=\frac{1}{2}$, so it describes fermions.

So we can choose a normalization for each $(A,B)$ field separately, and for
$c$ we can adjust the phase in the definition of creation and annilation
operators and $c=1$. As the result we have:

$\displaystyle\psi^{AB}_{ab}(x)=$ $\displaystyle(2\pi)^{-3/2}\sum_{\sigma}\int
d^{3}p[e^{ip\cdot
x}u^{AB}_{ab}(\mathbf{p},\sigma)a(\mathbf{p},\sigma)+e^{-ip\cdot
x}(-)^{2B}v^{AB}_{ab}(\mathbf{p},\sigma)a^{c\dagger}(\mathbf{p},\sigma)]$ (40)

##  4 Field equations for massive particles

###  4.1 General equation

For a quantum field $\psi_{l}^{+}(x)=(2\pi)^{-3/2}\int d^{3}pe^{ip\cdot
x}u_{l}(\mathbf{p},\sigma)a(\mathbf{p},\sigma)$, denote by

$\psi^{(+)AB}_{l}(p)=(2\pi)^{-3/2}u^{AB}_{l}(\mathbf{p},\sigma)a(\mathbf{p},\sigma)$
(41)

its Fourier transformation, and the similarly

$\psi^{(-)AB}_{l}(p)=(2\pi)^{-3/2}(-)^{2B}v^{AB}_{l}(\mathbf{p},\sigma)a^{\dagger}(\mathbf{p},\sigma)$
(42)

Clearly by unitarity of $D^{CD}|_{SU(2)}$, we have
$u^{CD\dagger}(0)u^{CD}(0)=I$, so we have

$\displaystyle\psi^{(+)AB}(p)$ $\displaystyle=\Pi^{AB,CD}(p)\psi^{(+)CD}(p)$
(43) $\displaystyle\psi^{(-)AB}(p)$
$\displaystyle=(-)^{2B+2D}\Pi^{AB,CD}(p)\psi^{(-)CD}(p)$

where we have defined:

######  Definition 4.1.

(Twisted spin sum)

Define the twisted spin sum

$\displaystyle\Pi^{AB,CD}(p)$ $\displaystyle\equiv
D^{AB}(Lp)u^{AB}(\mathbf{0})u^{CD\dagger}(\mathbf{0})D^{CD-1}(Lp)$ (44)
$\displaystyle=D^{AB}(Lp)v^{AB}(\mathbf{0})v^{CD\dagger}(\mathbf{0})D^{CD-1}(Lp)$

We have its transformation rule just like for spin sum:

######  Proposition 4.1.

$\Pi^{AB,CD}(\Lambda p)=D^{AB}(\Lambda)\Pi^{AB,CD}(p)D^{CD-1}(\Lambda)$ (45)

###### Proof.

$\displaystyle
D^{AB}(\Lambda)\Pi^{AB,CD}(p)D^{CD-1}(\Lambda)=D^{AB}(\Lambda)D^{AB}(Lp)u^{AB}(0)u^{CD\dagger}(0)D^{CD-1}(Lp)D^{CD-1}(\Lambda)$
(46) $\displaystyle=D^{AB}(L(\Lambda
p))D^{AB}(W(\Lambda,p))u^{AB}(0)u^{CD\dagger}(0)D^{CD-1}(W(\Lambda,p))D^{CD-1}(L(\Lambda
p))$ $\displaystyle=D^{AB}(L(\Lambda
p))D^{AB}(W(\Lambda,p))u^{AB}(0)u^{CD\dagger}(0)D^{CD\dagger}(W(\Lambda,p))D^{CD-1}(L(\Lambda
p))$ $\displaystyle=D^{AB}(L(\Lambda
p))u^{AB}(0)D^{j}(W(\Lambda,p))D^{j\dagger}(W(\Lambda,p))u^{CD\dagger}(0)D^{CD-1}(L(\Lambda
p))$ $\displaystyle=D^{AB}(L(\Lambda
p))u^{AB}(0)u^{CD\dagger}(0)D^{CD-1}(L(\Lambda p))$
$\displaystyle=\Pi^{AB,CD}(\Lambda p)$

where we have used the unitarity of $D^{CD}|_{SU(2)}$ in the third line. ∎

So we have the similar result for twisted spin sum:

######  Theorem 4.1.

In the massive case there is a 4-variables polynomial $P$ of coefficients
$(2A+1)(2B+1)\times(2C+1)(2D+1)$-matrices, such that on the mass shell its
value always equal the twisted spin sum

$\Pi^{AB,CD}(\mathbf{p},\sqrt{\mathbf{p}^{2}+m^{2}})=P(\mathbf{p},\sqrt{\mathbf{p}^{2}+m^{2}})$
(47)

and it is either an odd function or an even function determined by the parity
of $2A+2C$:

$P(-\mathbf{p},-p^{0})=(-1)^{2A+2C}P(\mathbf{p},p^{0})$ (48)

for all 4-vectors $p$. And the degrees $d$ of monomials in $P$ satisfies

$\max\\{|A-C|,|B-D|\\}\leq d/2\leq\min\\{A+C,B+D\\}$ (49)

For the proof we just modify the representation on $V$ by $M\mapsto
D^{AB}(\Lambda)MD^{CD-1}(\Lambda)$, and the LHS of the transformation rule for
$T$ matrices by $D^{AB}(\Lambda)T^{\cdots}D^{CD-1}$.

Note that $C,D$ are interchanged, because $\Lambda\mapsto
D^{CDT}(\Lambda^{-1})$ is isomorphic to $D^{CD}$, not $D^{DC}$ as in the spin
sum case.

At this point we can derive from (43) two field equations:

$\displaystyle\psi^{(+)AB}(x)$
$\displaystyle=\Pi^{AB,CD}(-i\partial)\psi^{(+)CD}(x)$ (50)
$\displaystyle\psi^{(-)AB}(x)$
$\displaystyle=(-)^{2A+2C}(-)^{2B+2D}\Pi^{AB,CD}(-i\partial)\psi^{(-)CD}(x)$

where we have abused the notation to let $\Pi^{AB,CD}(-i\partial)$ instead of
$P(-i\partial)$ denote replacing $p_{\mu}$ by $-i\partial_{\mu}$ in the
polynomial $P$. The two phases in the second line cancel, so we have the
following one field equation:

######  Theorem 4.2.

(Field equation of any spin, massive case) 101010[Tun67] neglect the phase of
negative energy part.  111111In the combination of positive and negative part,
the phase $(-)^{2B}$ plays an essential role.

The following field equation is satisfied:

$\psi^{AB}(x)=\Pi^{AB,CD}(-i\partial)\psi^{CD}(x)$ (51)

So clearly the field equation is not unique: we can adjust the polynomial $P$
by replacing $(p_{0})^{2}$ with $\mathbf{p}^{2}+m^{2}$ in any monomial, which
do not change the value on the mass shell. But this is just adding a term of
Klein-Gordan operator connected by several partial derivative operators.

The degrees of monomials in $P$ satisfy:

$\max\\{|A-C|,|B-D|\\}\leq d/2\leq\min\\{A+C,B+D\\}$ (52)

As a special case, the degrees of $\Pi^{j,0;j,0}$ and $\Pi^{0,j;0,j}$ are all
0. So we do not have a nontrivial field equation for $(j,0)$ or $(0,j)$ field
with itself from the above procedure. 121212This coincides with Weinberg. And
$\Pi^{j,0;0,j}$ and $\Pi^{0,j;j,0}$ must be homogeneous of degree $2j$. This
is the case that we have a homogeneous (except one term) field equation.
131313This is the result in Weinberg.

But in the general case, the polynomial may not be homogeneous.

Let’s talk about the relation between ordinary spin sum and twisted spin sum.
As is used in the proof of the field equation, $D^{CD-T}$ is isomorphic (but
not generally equal) $D^{CD}$, and $D^{CD*}$ is isomorphic to $D^{DC}$. So we
may expect that $\pi^{AB,CD}(p)$ is almost the same as $\Pi^{AB,DC}(p)$,
because their Lorentz transformation rule is under the same representation.

Define

$\Omega_{CD}:\begin{pmatrix}\mathbb{C}^{(2C+1)(2D+1)}\to\mathbb{C}^{(2C+1)(2D+1)}\\\
e_{(c-1)(2D+1)+d}\mapsto e_{(d-1)(2C+1)+c}\end{pmatrix}$ (53)

This is just a basis transformation

$\Omega_{CD}:\begin{pmatrix}V^{C}\otimes V^{D}\to V^{D}\otimes V^{C}\\\
v\otimes w\mapsto w\otimes v\end{pmatrix}$ (54)

So clearly

$\Omega_{CD}u^{CD}(\mathbf{0})=u^{DC}(\mathbf{0})$ (55)

Also we can show that $\Omega_{CD}$ is a homomorphism of representations:

$D^{DC-\dagger}(\Lambda)\Omega_{CD}=\Omega_{CD}D^{CD}(\Lambda)$ (56)

This is because $\exp(-ig)^{-\dagger}=\exp(-ig^{\dagger})$ implies that the
Lie algebra representation of $D^{DC-\dagger}$ is:

$\displaystyle\mathbf{J}\mapsto$ $\displaystyle\mathbf{J}^{D}\otimes
I^{C}+I^{D}\otimes\mathbf{J}^{C}$ (57) $\displaystyle\mathbf{K}\mapsto$
$\displaystyle+i(\mathbf{J}^{D}\otimes I^{C}-I^{D}\otimes\mathbf{J}^{C})$

So putting these together we have

$\displaystyle u^{DC\dagger}(\mathbf{0})D^{DC-1}(Lp)$
$\displaystyle=u^{CD\dagger}(\mathbf{0})\Omega_{CD}^{\dagger}\Omega_{CD}^{-\dagger}D^{CD\dagger}(Lp)\Omega_{CD}^{\dagger}$
(58)
$\displaystyle=u^{CD\dagger}(\mathbf{0})D^{CD\dagger}(Lp)\Omega_{CD}^{\dagger}$

Comparing with what we called the spin sum

$\pi^{AB,CD}(\mathbf{p})\equiv(2p^{0})u^{AB}(\mathbf{p})u^{CD\dagger}(\mathbf{p})=(2m)D^{AB}(Lp)u^{AB}(\mathbf{0})u^{CD\dagger}(\mathbf{0})D^{CD\dagger}(Lp)$
(59)

we derive:

######  Proposition 4.2.

(Relation between spin sum and twisted spin sum)

$\Pi^{AB,DC}(p)=\frac{1}{2m}\pi^{AB,CD}(p)\Omega_{CD}^{\dagger}$ (60)

Now we can formulate the field equation in terms of ordinary spin sum:

######  Theorem 4.3.

(Field equation of any spin, massive case)

The following field equation is satisfied:

$\psi^{AB}(x)=\frac{1}{2m}\pi^{AB,DC}(-i\partial)\Omega_{DC}^{\dagger}\psi^{CD}(x)$
(61)

In the special case where one of $C,D$ is 0, $\Omega_{CD}$ is the identity
matrix, so the twisted spin sum is the same as the spin sum. This is the case
in [Wei64a]. This can be seen more directly by noticing

$\displaystyle D^{0,j}(\Lambda)=D^{j,0.-\dagger}(\Lambda)$ (62) $\displaystyle
D^{j,0}(\Lambda)=D^{0,j.-\dagger}(\Lambda)$

is not only isomorphic as representations, but directly equal in their
standard basis.

###  4.2 Examples

The first example is massive spin-$\frac{1}{2}$ particles with field rep
$(\frac{1}{2},0)$ or $(0,\frac{1}{2})$. The corresponding field equation is
what we called Weyl equations or Dirac equation.

We can easily see that the zero momentum value coefficients and spin sums are:

$\begin{pmatrix}u^{1/2,0}(\mathbf{0})\\\
u^{0,1/2}(\mathbf{0})\end{pmatrix}=\begin{pmatrix}1&0\\\ 0&1\\\ 1&0\\\
0&1\end{pmatrix}$ (63)
$\begin{pmatrix}\pi^{1/2,0,1/2,0}(\mathbf{0})&\pi^{1/2,0,0,1/2}(\mathbf{0})\\\
\pi^{0,1/2,1/2,0}(\mathbf{0})&\pi^{0,1/2,0,1/2}(\mathbf{0})\end{pmatrix}=2m\begin{pmatrix}I_{2}&I_{2}\\\
I_{2}&I_{2}\end{pmatrix}$ (64)

By (11), the the upper-left term and the lower-right term when the momentum
varies equal a polynomial of degree $1$, while lower-left term and the upper-
right term when the momentum varies equal polynomials of degree 0. And by the
$T$ matrices for $(1/2,0,1/2,0)$ and $(0,1/2,0,1/2)$ computed in section 2,
the expansion of zero momentum spin sums are:

$\pi^{1/2,0,1/2,0}(\mathbf{0})=2m\sigma^{0},\quad\pi^{0,1/2,0,1/2}(\mathbf{0})=2m\bar{\sigma}^{0}$
(65)

so we have

$\begin{pmatrix}\pi^{1/2,0,1/2,0}(p)&\pi^{1/2,0,0,1/2}(p)\\\
\pi^{0,1/2,1/2,0}(p)&\pi^{0,1/2,0,1/2}(p)\end{pmatrix}=2\begin{pmatrix}-p_{\mu}\sigma^{\mu}&mI_{2}\\\
mI_{2}&-p_{\mu}\bar{\sigma}^{\mu}\end{pmatrix}$ (66)

We abbreviate $\varphi=\psi^{1/2,0},\chi=\psi^{0,1/2}$, and the field
equations (61) write:

$\displaystyle\varphi=\frac{1}{2m}\pi^{1/2,0,0,1/2}(-i\partial)\varphi$ (67)
$\displaystyle\varphi=\frac{1}{2m}\pi^{1/2,0,1/2,0}(-i\partial)\chi$
$\displaystyle\chi=\frac{1}{2m}\pi^{0,1/2,1/2,0}(-i\partial)\chi$
$\displaystyle\chi=\frac{1}{2m}\pi^{0,1/2,0,1/2}(-i\partial)\varphi$

The first and the third line are trivial:

$\displaystyle\varphi$ $\displaystyle=\varphi$ (68) $\displaystyle\chi$
$\displaystyle=\chi$

and the other two lines are:

$\displaystyle m\varphi$ $\displaystyle=i\sigma^{\mu}\partial_{\mu}\chi$ (69)
$\displaystyle m\chi$ $\displaystyle=i\bar{\sigma}^{\mu}\partial_{\mu}\varphi$

They are exactly Weyl equations for massive particles.

The second example is massive spin-$1$ particles with field rep
$(\frac{1}{2},\frac{1}{2})$. The corresponding field equation is what we
called Proca equation, also called the massive Maxwell equations.

The field rep $(\frac{1}{2},\frac{1}{2})$ is not directly equal, but is
equivalent to the vector representation (or called canonical representation)
$D(\Lambda)=\Lambda$ of $SO^{+}(1,3)$. We take this form.

The spin-$1$ rep of $su(2)$ is defined as:

$J_{z}\mapsto\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{pmatrix}\quad\quad\quad
J_{+}\mapsto\begin{pmatrix}0&\frac{1}{\sqrt{2}}&0\\\ 0&0&\frac{1}{\sqrt{2}}\\\
0&0&0\end{pmatrix}\quad\quad\quad J_{-}\mapsto\begin{pmatrix}0&0&0\\\
\frac{1}{\sqrt{2}}&0&0\\\ 0&\frac{1}{\sqrt{2}}&0\end{pmatrix}$ (70)

and the vector representation maps everything to itself:

$J_{z}=-i\begin{pmatrix}0&0&0&0\\\ 0&0&1&0\\\ 0&-1&0&0\\\
0&0&0&0\end{pmatrix}\quad\quad\quad
J_{\pm}=\frac{-i}{\sqrt{2}}\begin{pmatrix}0&0&0&0\\\ 0&0&0&\mp i\\\ 0&0&0&1\\\
0&\pm i&-1&0\end{pmatrix}$ (71)

Under the vector representation, the eigenvectors of z-axis angular-momentum
is:

$e_{0}=\begin{pmatrix}0\\\ 0\\\ 0\\\ 1\end{pmatrix}\quad\quad\quad
e_{+}=J_{+}e_{0}=-\frac{1}{\sqrt{2}}\begin{pmatrix}0\\\ 1\\\ i\\\
0\end{pmatrix}\quad\quad\quad
e_{-}=J_{-}e_{0}=\frac{1}{\sqrt{2}}\begin{pmatrix}0\\\ 1\\\ -i\\\
0\end{pmatrix}$ (72)

So the initial value of the coefficients is:

$u(0)=\begin{pmatrix}0&0&0\\\ -\frac{1}{\sqrt{2}}&0&\frac{1}{\sqrt{2}}\\\
-\frac{i}{\sqrt{2}}&0&-\frac{i}{\sqrt{2}}\\\ 0&1&0\end{pmatrix}$ (73)

and the initial value of the spin sum is:

$\pi(0)=2m\begin{pmatrix}0&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{pmatrix}$ (74)

We can show that

${}^{vec,vec}_{1}(T^{\mu\rho})^{\nu\sigma}=g^{\mu\sigma}g^{\nu\rho}$ (75)

Indeed we should verify:

$\Lambda^{\mu}{}_{\nu}(T^{ab})^{\nu}{}_{\rho}(\Lambda^{-1})^{\rho}{}_{\sigma}=\Lambda_{\rho}{}^{a}\Lambda_{\nu}{}^{b}(T^{\rho\nu})^{\mu}{}_{\sigma}$
(76)

This is equivalent to

$\Lambda^{\mu}{}_{\nu}g^{a\nu}g^{b}{}_{\rho}\Lambda_{\sigma}{}^{\rho}=\Lambda_{\rho}{}^{a}\Lambda_{\nu}{}^{b}g^{\mu\rho}g^{\nu}{}_{\sigma}$
(77)

which is obvious. Also clearly we have

${}^{vec,vec}_{0}T=I_{4}$ (78)

So the initial value of the spin sum is decomposed as:

$\pi(0)=2m[I_{4}-\begin{pmatrix}1&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{pmatrix}]=2m[{}^{vec,vec}_{0}(T)^{\circ\circ}-{}^{vec,vec}_{1}(T^{\mu\rho})^{\circ\circ}]$
(79)

So the field equation is

$(I_{4}-\frac{\partial_{\mu}\partial_{\rho}}{m^{2}}{}^{vec,vec}_{1}(T^{\mu\rho})^{\circ\circ})B_{\circ}=B_{\circ}$
(80)

It is equivalent to:

$\partial_{\mu}\partial_{\rho}g^{\mu\sigma}g^{\nu\rho}B_{\sigma}=0$ (81)

So the field equation becomes:

$\partial^{\nu}\partial^{\mu}B_{\mu}=0$ (82)

Taking another $\partial_{\nu}$ and using Klein-Gordan equation, we have:

$\partial^{\mu}B_{\mu}=0$ (83)

This, combined with Klein-Gordan equation, is equivalent to the famous Proca
equation:

$\partial_{\mu}(\partial^{\mu}B^{\nu}-\partial^{\nu}B^{\mu})+m^{2}B^{\nu}=0$
(84)

######  Remark 4.1.

Notice that for vector field of massless particles with helicity-$1$, (83) is
a gauge fixing condition we artificially add to the field, while in the
massive case it is a field equation that must be satisfied. So the vector
field for massive spin-$1$ particles do not have gauge invariance.

##  5 Conclusions

In this article we proved spin statistics for arbitrary massive $(A,B)$ field.
The quantity ’spin sum’ $\pi^{AB}(p)$ is also used to derive the massive field
equations, after a little modification. Klein-Gordan equation is just a
condition of 4-momentum in this context. Explicit calculation shows that
(Massive) Weyl equation, Dirac equation and Proca equation are all examples of
our general field equation.

In the case of $(A,B,C,D)=(j,0,0,j)$ or $(A,B,C,D)=(0,j,j,0)$, the polynomial
of the twisted spin sum is homogeneous, so the field equation contains a
constant term plus a homogeneous part, thanks to the condition on $K$. But in
the general case it may not be homogeneous.

In this article we just proved the ’existence’ part of the polynomial of the
spin sum, what remains undiscovered is an explicit formula for the
coefficients $\xi_{K}^{ABCD}$. Weinberg in [Wei69] discussed its explicit
formula.

Also in the massless case the proof in this article cannot be directly used.
One solution is to derive and explicit expansion for the spin sum
$\pi^{j,0}(p)$ and $\pi^{0,j}(p)$, as in [Wei64b], and then pass to the $m\to
0$ limit. Another is to notice that $(A,B)$ fields can be constructed out of
$(j,0)$ and $(0,j)$ fields.

Also, Rarita-Schwinger equation remains undiscovered in this article.

## References

  * [Fen23] Zixuan Feng.  Representation theory in the construction of free quantum field.  arXiv preprint arXiv:2302.13808, 2023.
  * [Tun67] Wu-Ki Tung.  Relativistic wave equations and field theory for arbitrary spin.  Physical Review, 156(5):1385, 1967.
  * [Wei64a] Steven Weinberg.  Feynman rules for any spin.  Physical Review, 133(5B):B1318, 1964.
  * [Wei64b] Steven Weinberg.  Feynman rules for any spin. ii. massless particles.  Physical Review, 134(4B):B882, 1964.
  * [Wei69] Steven Weinberg.  Feynman rules for any spin. iii.  Physical Review, 181(5):1893, 1969.
  * [Wei02] Steven Weinberg.  The quantum theory of fields: Foundations.  Cambridge University Press, 2002.


# Quantum Amplitude Estimation for Probabilistic Methods in Power Systems

Emilie Jong, Brynjar Sævarsson, Hjörtur Jóhannsson, Spyros Chatzivasileiadis
Department of Wind and Energy Systems
Technical University of Denmark (DTU)
Kgs. Lyngby, Denmark
{emijo, brysa, hjjo<EMAIL_ADDRESS>
###### Abstract

This paper introduces quantum computing methods for Monte Carlo simulations in
power systems which are expected to be exponentially faster than their
classical computing counterparts. Monte Carlo simulations is a fundamental
method, widely used in power systems to estimate key parameters of unknown
probability distributions, such as the mean value, the standard deviation, or
the value at risk. It is, however, very computationally intensive. Approaches
based on Quantum Amplitude Estimation can offer a quadratic speedup, requiring
orders of magnitude less samples to achieve the same accuracy. This paper
explains three Quantum Amplitude Estimation methods to replace the Classical
Monte Carlo method, namely the Iterative Quantum Amplitude Estimation (IQAE),
Maximum Likelihood Amplitude Estimation (MLAE), and Faster Amplitude
Estimation (FAE), and compares their performance for three different types of
probability distributions for power systems.

###### Index Terms:

Monte Carlo Simulation, Quantum Amplitude Estimation, IBM Qiskit

Submitted to the 23rd Power Systems Computation Conference (PSCC 2024).

##  I Introduction

The introduction of renewable energy sources in the power system together with
an increasingly fluctuating net electricity demand drastically increases the
degree of stochasticity in power system operation. To ensure a secure and
efficient operation of power systems, we need to estimate both accurately and
efficiently probability distributions for wind forecast errors, line loadings,
electricity demand patterns, voltage profiles, and many others. Key
parameters, such as the mean value, standard deviation, and value at risk are
important in order to extract information that will help assess the degree of
risk and determine the optimal power system operation. Traditionally, Monte
Carlo simulations have been widely used to estimate key parameters of
probability distributions, such as the mean value or the value at risk.
Considering, however, the high computational burden of such approaches, new
computational algorithms are extremely important to help drastically reduce
the computational time for these tasks.
Quantum Amplitude Estimation (QAE) algorithms, as first introduced in [1],
have proven to show a quadratic speed-up compared to Classical Monte Carlo
simulations (CMC). Although the initial algorithm [1] requires extensive
operations which make it unfeasible for current quantum computers due to
Quantum Phase Estimation (QPE) (i.e. quantum computers belonging to the Noisy
Intermediate-Scale Quantum era), in the previous years, variants on QAE have
been developed requiring less gates and qubits and thus making them more
suitable for implementation on near-term quantum computers. Among them, the
most promising are the Iterative Quantum Amplitude Estimation (IQAE) [2], the
Maximum Likelihood Amplitude Estimation (MLAE) [3], and the Faster Amplitude
Estimation (FAE) [4]. All algorithms require less samples to achieve the same
accuracy and confidence level as Classical Monte-Carlo Simulations (CMC) and
do not require the Quantum Phase Estimation. FAE was introduced as a faster
variant that required even less samples than the first version of IQAE. Two
years after its first introduction in 2019, IQAE was updated in 2021 [2] with
improved theoretical bounds, giving it a renewed theoretical advantage over
FAE. When it comes to MLAE, Ref. [5] compares the performance of IQAE and MLAE
in terms of oracle queries and circuit depth for a given degree of accuracy.
In this paper it is concluded that MLAE and IQAE are comparable in accuracy
and number of oracle calls and predominate classical Monte Carlo integration.
The goal of this paper is to introduce the three most promising Quantum
Amplitude Estimation algorithms for power systems applications by providing a
condensed overview of the research conducted in [6], and compare them
theoretically and empirically with regards to their efficiency and estimation
error. We implement them using Qiskit [7], the most widely used Python-based
toolbox for Quantum Computing, and explore their quantum advantage compared to
the Classical Monte Carlo Simulations. We apply the algorithms to different
probability distributions for power system applications to estimate parameters
such as the mean, the value at risk and conditional value at risk.
From a classical point of view, the required amount of samples for estimating
the mean with a certain accuracy and confidence level is largely influenced by
the type of probability distribution that it is applied on. Depending on the
power system application, stochastic values follow different types of
probability distributions. For example, wind forecasts typically follow a
Weibull distribution, while wind forecast errors are often considered to
follow a normal distribution with zero mean. Uniform distributions are also
widely used, e.g. for the distribution of voltage magnitudes across system
nodes. This paper introduces the three most promising quantum computing
variants which have the potential to replace classical Monte Carlo simulations
for power systems, and assesses them against Monte Carlo simulations for three
different types of probability distributions: normal, Weibull, and uniform. We
evaluate the performance of these algorithms based on the number of samples
they need, their accuracy, estimation errors, and circuit depth.
The remainder of the paper is structured as follows: In Section II, we delve
into the fundamental principles of quantum computing. Section III introduces
the Quantum Amplitude Estimation (QAE) algorithms and provides a theoretical
comparison. Section IV describes the practical application and implementation
of the algorithms. This is followed by the presentation of the results
obtained from simulations on quantum computers in Section V, including a
comparative analysis of the QAE algorithms. Moving forward to Section VI, we
discuss the limitations in this research. Finally, we conclude in Section VII,
summarizing the key findings and implications.

##  II Quantum Computing Fundamentals

Classical bits process information using bits, which are either 0 or 1, while
quantum bits (qubits) can be in a superposition state. The qubit’s state
remains uncertain until measured, at which point the qubit collapses into a
definite state. By making use of this property, quantum computers can evaluate
multiple states at the same time, allowing them to be more efficient than
their classical counterparts. A qubit is represented by the Dirac bra-ket
notation. A vector is denoted as $\ket{v}$, i.e. the ket. The bra, i.e
$\bra{v}$, is the conjugate transpose of $\ket{v}$. We can write the state of
qubit $\ket{\psi}$ as:

$\ket{\psi}=\alpha\ket{0}+\beta\ket{1}$ (1)

where $\alpha,\beta\in\mathbb{C}$, i.e they are complex numbers.
$|\alpha|^{2}$ and $|\beta|^{2}$ are the probability of the qubit being in
state $\ket{0}$ and $\ket{1}$ respectively, so it must hold that
$|\alpha|^{2}+|\beta|^{2}=1$. Quantum circuits make use of gates, which are
operations that are applied to qubits to change their states. Important gates
are the Hadamard gate, $H=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\\
1&-1\end{bmatrix}$, and the Pauli matrices:

$\displaystyle X=\begin{bmatrix}0&1\\\ 1&0\end{bmatrix},$ $\displaystyle
Y=\begin{bmatrix}0&-i\\\ i&0\end{bmatrix},$ $\displaystyle
Z=\begin{bmatrix}1&0\\\ 0&-1\end{bmatrix}.$

The more gates are used, the more the depth of the circuit increases. Although
Quantum Computing has been rapidly evolving and the prospects for quantum
computing applications look bright, we still find ourselves in the Noisy
Intermediate-Scale Quantum (NISQ) era, where the circuit depth matters
significantly due to the amount of noise generated by the gates. Meanwhile,
publicly available quantum computers during this research were limited to only
7 qubits, while in October 2023 quantum computers with up to 127 qubits became
publicly available for a time-limited use.

##  III Quantum Amplitude Estimation

Quantum Amplitude Estimation is the quantum alternative of the Monte Carlo
method. Theoretically, QAE shows a quadratic speed-up over classical
algorithms [8]. This means that instead of requiring 100’000 samples to
estimate the mean value of a distribution, we would only need approx. 320
samples to estimate the mean with the same accuracy. In 2000, the initial
canonical QAE version was invented, however it required a form of Quantum
Phase Estimation (QPE) resulting in a very deep circuit. In the previous
years, several other algorithms have been introduced, which are more suitable
for quantum computers in the Noisy Intermediate-Scale Quantum (NISQ) era. This
section introduces first the canonical version, and then the IQAE, MLAE and
FAE algorithms.

###  III-A Canonical Amplitude Estimation – Proof of the Quadratic Speedup

QAE was first introduced in [1] using Grover’s search algorithm and QPE. It
works as follows. Suppose we have an operator $\mathcal{A}$ acting on $n+1$
qubits resulting in:

$\ket{\psi}=\mathcal{A}\ket{0}_{n+1}=\sqrt{1-a}\ket{\psi_{0}}_{n}\ket{0}+\sqrt{a}\ket{\psi_{1}}_{n}\ket{1}$
(2)

The application of $\mathcal{A}$ corresponds to an oracle query (quantum
sample). Now, we can find an estimate for the amplitude of $a$ for state
$\ket{\psi_{1}}$ using the Grover operator:
$\mathcal{Q}=\mathcal{A}\mathcal{S}_{0}\mathcal{A}^{\dagger}\mathcal{S}_{\psi_{0}}$.
$\mathcal{S}_{0}$ and $\mathcal{S}_{\psi_{0}}$ are reflections, defined as
$\mathcal{S}_{0}=\mathbb{I}-2\ket{0}_{n+1}\bra{0}_{n+1}$ and
$\mathcal{S}_{\psi_{0}}=\mathbb{I}-2\ket{\psi_{0}}_{n}\bra{\psi_{0}}_{n}\otimes\ket{0}\bra{0}$
with $\mathbb{I}$ being the identity operator. For canonical QAE $m$ ancilla
bits111Ancilla bits are auxiliary bits which are used as “flags” in quantum
operations. are needed to represent the final result, which are initialized in
a superposition state created by Hadamard gates. The number of quantum samples
is defined as $M=2^{m}$. The ancilla bits control the increasing powers of
$\mathcal{Q}$. Then a Quantum Fourier Transform (QFT) is applied to the
ancilla bits, after which they are measured. The state of the ancilla bits
corresponds to an integer $y\in\\{0,...,M-1\\}$ which is mapped to an angle
$\theta_{a}=y\pi/M$. The estimate for $a$,
$\tilde{a}=\sin^{2}{\theta_{a}}\in[0,1]$. Using this, (2) can be rewritten as
with $0\leq\theta_{a}\leq\pi/2$:

$\ket{\psi}=\mathcal{A}\ket{0}_{n+1}=\cos{\theta_{a}}\ket{\psi_{0}}_{n}\ket{0}+\sin{\theta_{a}}\ket{\psi_{1}}_{n}\ket{1}$
(3)

The estimate satisfies the bound error with a probability of
$\frac{8}{\pi^{2}}\approx 81\%$:

$\frac{2\sqrt{a(1-a)}\pi}{M}+\frac{\pi^{2}}{M^{2}}=\mathcal{O}\left(\frac{1}{M}\right)$
(4)

This shows a quadratic speedup compared to classical Monte Carlo methods which
converge with $\mathcal{O}\left(\frac{1}{\sqrt{M}}\right)$ [8]. In [1], it is
shown that after $k$ applications of $\mathcal{Q}$ on $\ket{\psi}$, we get:

$\mathcal{Q}^{k}\ket{\psi}=\cos{((2k+1)\theta_{a})}\ket{\psi_{0}}\ket{0}+\sin{((2k+1)\theta_{a})}\ket{\psi_{1}}\ket{1}$
(5)

The canonical algorithm then uses QPE which requires many controlled
$\mathcal{Q}$ operation and additional ancilla qubits. IQAE, MLAE and FAE
differ from the canonical algorithm by eliminating the need for quantum phase
estimation through a combination of post-processing techniques and by applying
different powers of the Grover operators.

###  III-B Iterative Quantum Amplitude Estimation

In 2019, a new quantum amplitude estimation algorithm without phase estimation
is proposed in [2], called Iterative Quantum Amplitude Estimation (IQAE). It
uses Grover iterations to find an estimate for the target amplitude. IQAE
starts with a small $\mathcal{Q}^{0}$ circuit which has one solution. This
solution is then used to evaluate correct solutions from deeper circuits with
more precise results (with higher powers of the Grover operator). It
iteratively applies rotations to find the correct solution of the deeper
circuits by using the outcomes of smaller circuits. It does this until it
obtains an estimate with a certain confidence level. The advantage of IQAE is
that it is possible to specify the desired confidence level. However, the
circuit depth of IQAE cannot be controlled.

###  III-C Maximum Likelihood Amplitude Estimation

The Maximum Likelihood Amplitude Estimation algorithm implements QAE without
QPE, by using maximum likelihood as a post-processing method. Maximum
likelihood is a statistical method used to estimate the parameters of a
probability distribution by obtaining the values that would be the most
probable when assuming that the observed data follows a certain distribution.
In Ref. [3] it has been proven that the MLAE algorithm achieves nearly the
optimal quantum speedup with a reasonable circuit length allowing MLAE to be
executed on NISQ devices. MLAE runs experiments with multiple
$\mathcal{Q}^{k}$ circuits. MLAE takes the likelihood functions of the results
of these circuits and then applies the classical maximum likelihood estimation
to retrieve an estimate for $\theta$. MLAE can thus be parallelized. There are
two ways of finding the $k$ values for the power of the Grover operator:
linearly and exponentially. The exponential approach results in a deeper
circuit, however, it provides a more accurate estimation. The amount of query
calls to achieve estimation error $\epsilon$ for MLAE is implied to be
$\mathcal{O}\left(\frac{1}{\epsilon^{3/4}}\right)$ and
$\mathcal{O}\left(\frac{1}{\epsilon}\right)$ for linearly and exponentially
incremental sequence respectively compared to
$\mathcal{O}\left(\frac{1}{\epsilon^{2}}\right)$ for the required samples of
the classical approach in order to achieve the same precision. Because of the
more accurate estimation of the exponentially incremental sequence, we will
consider just the exponentially incremental sequence.

###  III-D Faster Amplitude Estimation

Faster Amplitude Estimation (FAE) is proposed in [4]. This algorithm is
similar to IQAE, however, it solves an ambiguity differently in the post-
processing. For FAE, it is possible to fix the amount of iterations to control
the circuit depth as opposed to IQAE.

###  III-E Theoretical Quantum Advantage

In quantum, a sample is defined as the amount of times $\mathcal{A}$ is
applied and is called an oracle query or oracle call. The estimation accuracy
depends on the confidence level $\alpha$ and the estimation error $\epsilon$.
The estimation error is the difference between the actual value and the
estimated value of the parameter. The confidence level indicates the
probability that the estimate has a maximal estimation error of $\epsilon$ and
is defined as: 100(1-$\alpha$)%. For example, if $\alpha$ is set to 0.05,
$P[|\theta-\hat{\theta}|\leq\epsilon]\geq 95\%$, where $\hat{\theta}$ is the
estimate of $\theta$.

For MLAE, the lower bound for the amount of samples $N$ is given by (6).

$N_{min}^{MLAE}\geq\frac{\sqrt{\alpha(1-\alpha)}}{\epsilon}$ (6)

For the other algorithms, upper bounds are provided. For IQAE, different upper
bounds are presented in [2]. The loose upper bound presented in the paper is:

$N_{max}^{IQAE}<\frac{50}{\epsilon}\log{\left(\frac{2}{\alpha}\log_{2}{\left(\frac{\pi}{4\epsilon}\right)}\right)}$
(7)

Another smaller upper bound is given for IQAE, which is based on an empirical
complexity analysis: Clopper-Pearson (CP). The Clopper-Pearson bound is:

$N_{max}^{CP-
avg}\leq\frac{0.8}{\epsilon}\log{\left(\frac{2}{\alpha}\log_{2}{\left(\frac{\pi}{4\epsilon}\right)}\right)}$
(8)

Finally, the upper bound for FAE is:

$N_{max}^{FAE}<\frac{4.1\cdot
10^{3}}{\epsilon}\log{\left(\frac{2}{\alpha}\log_{2}{\left(\frac{2\pi}{3\epsilon}\right)}\right)}$
(9)

For classical Monte Carlo, the relationship between the accuracy and the
amount of samples follows the central limit theorem:

$N_{min}^{CMC}\sim\frac{z_{c}(\alpha)^{2}s_{n}^{2}}{\epsilon^{2}},$ (10)

where the value of $z_{c}$ depends on the confidence level $\alpha$ and
$s_{n}$ is the sample standard deviation of the distribution.

Fig. 1 shows a graphical comparison of the algorithms. Here, for CMC,
$z_{c}=1.96$ for a confidence interval of 95% ($\alpha=0.05$) and
$s_{n}^{2}=1$, but in reality $s_{n}$ depends on the sampling and on the type
of probability distribution, e.g. the more complex and asymmetric the
distribution, the larger $s_{n}$. The graph depicted for CMC is thus more an
indication for its trend, since this graph shifts upwards or downwards
depending on the $s_{n}$. It should be noted that for all algorithms, except
MLAE and CMC, upper bounds are shown, while for MLAE a lower bound is
visualized.

Figure 1: Comparison of the theoretical speedup of QAE algorithms and the
classical Monte Carlo approach ($\alpha=0.05$) [6]

##  IV Implementation

###  IV-A Method for estimating the mean, VaR and CVaR

The mean value, the Value-at-Risk (VaR) and the Conditional-Value-at-Risk
(CVaR) are the statistical parameters of interest for the investigated
probability distributions, because they provide a measure for risk assessment.
The VaR is a measure that estimates the maximum value in a system for a
specified confidence level, while the CVaR measures the expected value if the
VaR threshold is exceeded. The methods as explained in [9] are used to
evaluate the mean, VaR and CVaR. In order to be able to estimate these
parameters, the qubits need to be initialized so that they represent the
probability distribution. Suppose we have a random variable $X$ that follows a
certain probability distribution. The probability distribution has to be
discretized into $2^{n}$ bins so that it can be represented by $n$ qubits in
the interval [0, …, N - 1], where $N=2^{n}$. The discrete probability
distribution has to be loaded into the amplitudes of the qubits. The quantum
state for $X$ is then described as:

$\ket{\psi}_{n}=\sum_{x=0}^{N-1}\sqrt{p_{x}}\ket{x}_{n},$ (11)

where $p_{x}$ is the probability of the state $\ket{x}_{n}$ and the sum of the
probabilities of all states should be 1. This way, the sampling probability,
which is the square of the amplitude, equals the probability of the
distribution. The linear amplitude function $F$ implements piecewise linear
functions on qubit amplitudes constructed by linearly controlled Pauli-Y
rotations. The operator $F$ acts on an ancilla bit and maps the function
$f\in[0,...,2^{n-1}]$ to:

$F|x\rangle|0\rangle=\sqrt{1-f(x)}|x\rangle|0\rangle+\sqrt{f(x)}|x\rangle|1\rangle.$
(12)

Applying $F$ to $\ket{\psi_{n}}\ket{0}$ gives:

$\sum_{x=0}^{N-1}\sqrt{1-f(x)}\sqrt{p_{x}}\ket{x}_{n}\ket{0}+\sum_{x=0}^{N-1}\sqrt{f(x)}\sqrt{p_{x}}\ket{x}_{n}\ket{1}$
(13)

The probability of measuring $\ket{1}$ in the rightmost qubit is then equal to
the expected value of $f(x)$.
To estimate the value at risk, an operator, $\mathcal{F}_{l}$, is constructed
which forms a cumulative distribution function. $F_{l}$ implements the
function $f_{l}$ where $f_{l}(x)=1$ for $x\leq l$ and $f_{l}(x)=0$ otherwise.
Applying $\mathcal{F}_{l}$ to $\ket{\psi}_{n}\ket{0}$ results in:

$\sum_{x=l+1}^{N-1}\sqrt{p_{x}}\ket{x}\ket{0}+\sum_{x=0}^{l}\sqrt{p_{x}}\ket{x}\ket{1}$
(14)

By using QAE, we can estimate the probability of obtaining $\ket{1}$ which is
equal to $\sum_{x=0}^{l}p_{x}=P[X\leq l]$. Then, classically, we perform a
bisection search. The bisection search finds the smallest value $l_{\alpha}$
for which $P[X\leq l_{\alpha}]\geq 1-\alpha$. Quantum amplitude estimation is
specifically suitable for estimating tail probabilities of distributions [9].
This holds, because for the VaR, we can replace $a$ in (4) for $\alpha$, since
the estimated probability is larger than or equal to $\alpha$. If $\alpha$ is
close to either 0 or 1, the error bound becomes very small. The error bound is
then independent of other properties of the problem.

For estimating the CVaR, we identify the function
$f(x)=\frac{x}{l_{\alpha}}f_{l}$. Applying the corresponding operator
$\mathcal{F}$ yields:

$F\ket{\psi}_{n}\ket{0}=\left(\sum_{x=l_{\alpha}+1}^{N-1}\sqrt{p_{x}}\ket{x}+\sum_{x=0}^{l_{\alpha}}\sqrt{1-\frac{x}{l_{\alpha}}}\sqrt{p_{x}}\ket{x}\right)\ket{0}\\\
+\sum_{x=0}^{l_{\alpha}}\sqrt{\frac{x}{l_{\alpha}}}\sqrt{p_{x}}\ket{x}\ket{1}$
(15)

Then, we apply QAE to estimate the CVaR. Applying QAE gives
$\sum_{x=0}^{l_{\alpha}}\frac{x}{l_{\alpha}}p_{x}$, where
$\sum_{x=0}^{l_{\alpha}}p_{x}$ is $P[X\leq l_{\alpha}]$. Then, it follows
that:

$CVaR(X)=\frac{l_{\alpha}}{P[X\leq
l_{\alpha}]}\sum_{x=0}^{l_{\alpha}}\frac{x}{l_{\alpha}}p_{x}$ (16)

###  IV-B Simulation procedure

The QAE algorithms (IQAE, MLAE and FAE) are applied to the normal, Weibull and
uniform distributions to find an estimate for the mean, VaR and CVaR. The
parameters of the distributions are often based on historical data [10]. Here,
we will follow a more general approach for determining the parameters. Since
the probability distribution is loaded on the amplitude of the qubits, the
distribution can be scaled as desired by multiplying the qubit states
accordingly. The normal distribution is denoted as $N(\mu,\sigma)$ with $\mu$
being the mean of the distribution and $\sigma$ the standard deviation. We
investigate two normal distributions with the same mean ($\mu=0.1$) and
different standard deviations i.e. $\sigma=0.01$ and $\sigma=0.05$. Normal
distributions are applied in power systems to model uncertainty. We set the
standard deviation to 10% and 50% as percentage of the mean value. In [11],
for the uncertainty modelling of wind farms, the standard deviation can be up
to 50% of the mean value, while for load modelling the standard deviation is
assumed to be 5% of the mean value. Another study fixes the standard deviation
to 7% of the mean value [12]. The Weibull distribution is referred to as
$W(\beta)$, where $\beta$ is the shape parameter of the Weibull slope. We
follow a general approach and just consider the shape parameter to model the
shape of the distribution, while the scale parameter is set to 1. The shape
parameter depends on the wind speed data, namely the average wind speed and
the standard deviation [13]. We choose $\beta=1.8$ for the Weibull
distribution, since this parameter is typically between 1.63 to 2.97 [14],
determining the shape of the distribution. The uniform distribution is defined
as $U(a,b)$ where $U$ is uniform on the interval $[a,b)$. The interval that is
investigated is $[0,1)$. The probability distribution is loaded onto 4 qubits.
There is a trade-off between accuracy and runtime, while taking the limited
amount of qubits in publicly available quantum computers into account.
For the QAE algorithms, the functions ‘IterativeAmplitudeEstimation’,
‘MaximumLikelihoodEstimation’ and ‘FasterAmplitudeEstimation’ are used, which
are available in Qiskit. To account for the varying scaling of the different
probability distributions, we assess the relative errors of the application of
different QAE algorithms applied to the different probability distributions in
comparison to the actual distribution parameters. All algorithms are simulated
10 times on the quantum computer simulator to find the mean, VaR and CVaR for
each QAE algorithm for all probability distributions. Furthermore, some tests
are ran on a real quantum computer of IBM (IBM Guadalupe). Then, a more
detailed comparison of the QAE algorithms versus CMC is performed for the
estimate of the mean. For the comparison of QAE algorithms when estimating the
mean, the average relative estimation error is computed for a range of oracles
from $\sim 10^{2}$ to $\sim 2\cdot 10^{5}$. In [15], it is suggested that
10000 iterations would provide sufficiently accurate results for probabilistic
(optimal) power flow problems, however, is stated to be computationally
intensive. In [16], it is mentioned that tens of thousands of deterministic
power flow calculations are needed for convergence with random samples
generated from the probability distributions of the input variables. For this
reason, the range that is investigated, is from order of magnitude $10^{2}$ to
$10^{5}$. Thereafter, the practical confidence interval is computed to
evaluate the spread of the estimates.

##  V Simulation

In Table I, the average amount of oracle queries is shown with the
corresponding relative error for IQAE, MLAE and FAE when estimating the mean
of a normal distribution $N(0.1,0.01)$. Here, we can observe that the order of
magnitude for the average error differs, as well as the number of oracle
queries.

TABLE I: Oracle queries and relative error (%) for estimating the mean of N($\mu=0.1,\sigma=0.01$) ($\alpha=0.05$, $N_{shots}=100$, $\epsilon=10^{-3}$, $q=4$, $max\\_iter=3$) QAE algorithm | Average oracle queries |  $\Delta_{avg}$ (%)
---|---|---
IQAE | 15960 | $0.03132$
MLAE | 700 | $0.70745$
FAE | 136238 | $1.55281$

For a more consistent comparison of the algorithms, we assess the relative
error for estimating the mean of the four probability distributions for each
of the QAE algorithms and plot them using the same range of oracle queries.
This is shown in Fig. 2. For a 1% relative error when estimating the mean, the
QAE algorithms do not seem to have any advantage over classical Monte Carlo in
terms of oracle queries for the normal distribution N(0.1, 0.01). This could
be the case because the standard deviation of the normal distribution is
relatively small. Starting from a relative estimation error of 0.3%, the IQAE
and MLAE algorithms overtake CMC. Both algorithms show similar trends for all
probability distributions regarding the average relative error and the amount
of oracle queries. For N(0.1, 0.05) the quantum advantage for IQAE and MLAE is
already clear for a relative error of around 10%. The graph indicates that a
smaller amount of oracles are needed when compared to the amount of samples
needed for CMC to achieve the same accuracy. The advantage increases
monotonically when we require smaller relative errors. Remarkably, the
performance of FAE for N(0.1, 0.05) is much closer to CMC than for N(0.1,
0.01).

Figure 2: The relative error vs the amount of oracles plotted for estimating
the mean of different types of probability distributions

Tables II and III show the results for estimating the VaR and the CVaR for the
different distributions. The VaR levels deviate a few percentage points, since
only 4 qubits are used for loading the probability distribution. Increasing
the amount of qubits from 4 to 7 significantly decreases the relative error.
For the uniform distribution, 4 qubits are not sufficient to find a VaR of
95%. The average relative error for the CVaR is around 10 times larger than
for the VaR. Again, increasing the amount of qubits for the loading of the
probability distribution to 7, significantly decreases the average relative
error.
From Fig. 2, it can be concluded that the IQAE and MLAE algorithms show
similar performance in terms of accuracy and required amount of oracle calls.
For further comparison, we investigate the confidence interval for the
estimates that IQAE, MLAE and CMC provide. Fig. 3 shows the average confidence
interval for the estimates given by IQAE, MLAE and CMC. FAE has been left out
of the graph for readibility reasons due to its significantly larger
confidence interval compared to MLAE, IQAE and CMC. For samples in the order
of magnitude of $10^{2}$, the confidence interval of the CMC is smaller than
for the QAE algorithms. However, from $\sim 10^{3}$ oracles, IQAE and MLAE
overtake CMC. The confidence interval of IQAE is the smallest, and thus gives
the most accurate estimate.

TABLE II: Results for estimating the VaR for different distributions ($\alpha=0.05$, $N_{shots}=100$, $\epsilon=10^{-3}$, $q=4$, $max\\_iter=3$) including the actual VaR’s | IQAE | MLAE | FAE
---|---|---|---
N (0.1, 0.01), $VaR_{95\%}=0.1164$
Oracles | 26840 | 700 | 121907
VaR estimate | 0.1233 | 0.1233 | 0.1233
VaR level (%) | 97.92 | 98.77 | 97.92
Actual VaR for VaR level | 0.1204 | 0.1225 | 0.1204
$\Delta_{VaR}$ (%) | 2.41 | 0.65 | 2.41
N (0.1, 0.05), $VaR_{95\%}=0.1822$
Oracles | 23890 | 700 | 121907
VaR estimate | 0.2100 | 0.2080 | 0.2040
VaR level (%) | 97.86 | 97.01 | 97.01
Actual VaR for VaR level | 0.2013 | 0.1941 | 0.1941
$\Delta_{VaR}$ (%) | 4.14 | 7.16 | 5.10
Weibull ($\beta=1.8$), $VaR_{95\%}=1.8396$
Oracles | 18030 | 700 | 121907
VaR estimate | 2.0667 | 2.0893 | 2.0667
VaR level (%) | 96.52 | 96.66 | 96.50
Actual VaR for VaR level | 1.9601 | 1.9734 | 1.9582
$\Delta_{VaR}$ (%) | 5.44 | 5.87 | 5.54
Uniform [0, 1), $VaR_{95\%}=0.95$
Oracles | 19890 | 700 | 121907
VaR estimate | 1.0000 | 1.0000 | 1.0000
VaR level (%) | 100 | 100 | 100
TABLE III: Results for estimating the CVaR for different distributions ($\alpha=0.05$, $N_{shots}=100$, $\epsilon=10^{-3}$, $q=4$, $max\\_iter=3$) including actual CVaRs | IQAE | MLAE | FAE
---|---|---|---
N (0.1, 0.01)
Oracles | 18390 | 700 | 268924
CVaR estimate | 0.1407 | 0.1390 | 0.1472
VaR level (%) | 97.92 | 98.77 | 97.92
Actual CVaR for VaR level | 0.1241 | 0.1259 | 0.1241
$\Delta_{CVaR}$ (%) | 13.40 | 16.85 | 12.08
N (0.1, 0.05)
Oracles | 19420 | 700 | 280744
CVaR estimate | 0.1875 | 0.1877 | 0.1873
VaR level (%) | 97.86 | 97.01 | 97.01
Actual CVaR for VaR level | 0.2198 | 0.2135 | 0.2125
$\Delta_{CVaR}$ (%) | 14.69 | 12.06 | 11.86
Weibull ($\beta=1.8$)
Oracles | 19470 | 700 | 277788
CVaR estimate | 2.7305 | 2.7165 | 2.604
VaR level (%) | 96.52 | 96.66 | 96.50
Actual CVaR for VaR level | 2.2528 | 2.651 | 2.2515
$\Delta_{CVaR}$ (%) | 21.21 | 19.93 | 15.64

(a) The relative deviation of the relative estimate from the actual mean and
confidence interval for estimating the mean of N(0.1, 0.05)

(b) Zoom-in on the relative deviation of the relative estimate from the actual
mean and confidence interval for estimating the mean of N(0.1, 0.05) for range
$10^{3}-10^{6}$

Figure 3: A zoom-in on the relative deviation of the relative estimate from
the actual mean and confidence interval for estimating the mean of N(0.1,
0.05)

(a) Estimating the mean using QAE with 1 iteration (k=1)

(b) Estimating the mean using QAE with 5 iterations (k=5)

Figure 4: Circuits with different iterations for QAE for loading the
probability distribution onto 4 qubits, $\mathcal{Q}$ is the Grover operator

The QAE algorithms have also been ran on a real IBM quantum computer (IBM
Guadalupe) in order to estimate the mean value. However, due to the noise
already generated in the gates of the circuit for the construction of the
probability distribution, there are discrepancies in the representation of the
probability distributions. This is why the CVaR and VaR have not been
estimated on a real quantum computer. The construction of the probability
distribution can already take up to 62 gates (in case of the Weibull
distribution) and applying IQAE requires around 300 more gates. This does not
yet allow for accurate results on a real quantum computer.

##  VI Discussion

The results show potential for the application of QAE to estimate the mean,
VaR and CVaR for probability distributions relevant to power systems. The
empirical comparison shows that MLAE and IQAE provide a speed-up compared to
classical Monte Carlo simulation. For the tested probability distributions,
the IQAE and MLAE algorithms were able to provide a more accurate estimate for
the mean using less samples for all distributions. In case of N(0.1, 0.01),
MLAE and IQAE gained an advantage over CMC around $\sim 3\cdot 10^{3}$
samples, while for the other distributions, this happened from $\sim 3\cdot
10^{2}$ samples. FAE, however, did not show any advantage for the investigated
range of samples.

For the normal distribution and the uniform distribution similar trends are
visible for the empirical and theoretical comparison while for the Weibull
distribution at some point, the quantum advantage decreases. One of the
reasons for this might be that the mean can not be estimated more accurately
when more oracles are used due to the limited amount of qubits with which the
probability distribution is represented. It could be that for the other
distributions, the same happens when investigating an even larger range of
oracles. It is assumed that increasing the amount of qubits would solve this
problem.

The accepted average relative error depends on the type of application for
which they are used. Although IQAE seems to be the preferred algorithm in
terms of number of oracle queries and accuracy, it should be noted that its
circuit depth cannot be controlled, while for MLAE and FAE this is possible.
For current quantum computers, the depth of a circuit still matters
significantly due to the amount of noise generated by the quantum computers.
The circuit depth increases drastically when the amount of Grover operations
increases. For example, estimating the mean of $N(0.1,0.01)$ using 1
iteration, the amount of required gates would be 299, while for 5 iterations,
this number incresases to 1183. Fig. 4 shows the circuits with 1 and 5
iterations respectively. Additionally, as the number of qubits within the
circuit grows, the circuit depth increases as well.

Although, IQAE and MLAE have proven to provide a speed-up compared to
classical CMC empirically, it has to be noted that this advantage is still not
visible in terms of runtime. While it takes several minutes to do one
iteration of a QAE algorithm on the quantum computer simulator, it takes
around 10-20 seconds to perform it on a quantum computer. However, CMC just
needs several seconds in classical computation time. The advantage is thus
purely seen in computation complexity.

While in this paper, we have showed that we can apply QAE to known probability
distributions, in a typical Monte Carlo simulation the probability
distribution is not known a priori. The next step involves loading a truly
unknown distribution on the qubits.

##  VII Conclusion

In this paper, we have introduced quantum computing algorithms for Monte Carlo
simulations which have the potential to find wide use in power system
applications. Monte Carlo simulations is a standard method employed in both
power system literature and practice to measure the mean value or the value at
risk of unknown distributions such as the wind forecast error or the line
limit violations (among several others). Quantum Computing algorithms, based
on the Quantum Amplitude Estimation method, have been theoretically proven to
have a quadratic advantage versus classical Monte Carlo methods. For example,
instead of requiring 100’000 samples to estimate the mean value of a
distribution, they can use only approx. 320 samples to estimate it with the
same level of accuracy.

In this paper, we have described and tested the three most promising quantum
algorithms that could potentially replace the Classical Monte Carlo
simulations: the IQAE (Iterative Quantum Amplitude Estimation), the MLAE
(Maximum Likelihood Amplitude Estimation), and the FAE (Faster Amplitude
Estimation).

Besides their theoretical advantage, we have shown that all quantum algorithms
have demonstrated a practical quantum advantage compared to the Classical
Monte Carlo method, when tested on a simulated Quantum Computer. Among them,
IQAE demonstrated the most accurate estimates for a desired level of accuracy.

The contribution of this work is to introduce the quantum algorithms which
have the potential to drastically improve the Classical Monte Carlo
simulations in the power systems community and test them for distributions
that we often encounter in power system applications. Still, there is
significant work necessary both for power system engineers and for quantum
computing engineers to bring these methods into practice. Current quantum
computers are at the Noisy Intermediate Scale Era. Although we attempted tests
in real quantum computers, for the comparatively larger size of circuits used
in the methods we presented (300 gates or more) the level of noise in real
quantum computers does not allow for obtaining useful results. For the same
reason, although we have shown that the number of ”oracles”, i.e. iterations,
is indeed substantially lower in the quantum algorithms, quantum algorithms
still do not run faster than Classical Monte Carlo simulations. Quantum
scientists and engineers are working intensively to drastically reduce the
noise level within the next few years, which will have a positive impact both
in the number of gates that a quantum computer can handle and in reducing the
runtime. Considering that real Quantum Computers are expected to deliver
significant practical advantages very soon, we think that power system
engineers shall already start developing the required methods that can address
effectively power system challenges.

## References

  * [1] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, “Quantum amplitude amplification and estimation,” _AMS Contemporary Mathematics Series_ , vol. 305, no. 4, June 2000.
  * [2] D. Grinko, J. Gacon, C. Zoufal, and S. Woerner, “Iterative quantum amplitude estimation,” _npj Quantum Information_ , vol. 7, no. 1, March 2021.
  * [3] Y. Suzuki, S. Uno, R. Raymond, T. Tanaka, T. Onodera, and N. Yamamoto, “Amplitude estimation without phase estimation,” _Quantum Information Processing_ , vol. 19, January 2020.
  * [4] K. Nakaji, “Faster amplitude estimation,” _Quantum Information and Computation_ , vol. 20, no. 13&14, pp. 1109–1123, November 2020.
  * [5] K. Yu, H. Lim, P. Rao, and D. Jin, “Comparison of amplitude estimation algorithms by implementation,” 2020, (unpublished).
  * [6] E. Jong, “Quantum computing for probabilistic power flow in power systems,” Master’s thesis, Technical University of Denmark, Lyngby, July 2023.
  * [7] Qiskit, “An open-source framework for Quantum Computing,” Online. [Online]. Available: https://qiskit.org/
  * [8] S. Woerner and D. Egger, “Quantum risk analysis,” _npj Quantum Information_ , vol. 5, 12 2019.
  * [9] D. J. Egger, R. García Gutiérrez, J. C. Mestre, and S. Woerner, “Credit risk analysis using quantum computers,” _IEEE Transactions on Computers_ , vol. 70, no. 12, pp. 2136–2145, 2021.
  * [10] U. H. Ramadhani, M. Shepero, J. Munkhammar, J. Widén, and N. Etherden, “Review of probabilistic load flow approaches for power distribution systems with photovoltaic generation and electric vehicle charging,” _International Journal of Electrical Power & Energy Systems_, vol. 120, p. 106003, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0142061519341730
  * [11] M. Aien, M. Fotuhi-Firuzabad, and F. Aminifar, “Probabilistic load flow in correlated uncertain environment using unscented transformation,” _IEEE Transactions on Power Systems_ , vol. 27, no. 4, pp. 2233–2241, 2012.
  * [12] N. Bhat, B. R. Prusty, and D. Jena, “Cumulant-based correlated probabilistic load flow considering photovoltaic generation and electric vehicle charging demand,” _Frontiers in Energy_ , vol. 11, pp. 184–196, 05 2017.
  * [13] S.-E. Gryning, R. Floors, A. Peña, E. Batchvarova, and B. Brümmer, “Weibull wind-speed distribution parameters derived from a combination of wind-lidar and tall-mast measurements over land, coastal and marine sites,” _Boundary-Layer Meteorology_ , vol. 159, no. 2, pp. 329–348, 2016, .
  * [14] Z. Shu and M. Jesson, “Estimation of weibull parameters for wind energy analysis across the uk,” _Journal of Renewable and Sustainable Energy_ , vol. 13, p. 023303, 03 2021.
  * [15] W. Sun, M. Zamani, H.-T. Zhang, and Y. Li, “Probabilistic optimal power flow considering correlation of wind farms via markov chain quasi-monte carlo sampling,” 2018.
  * [16] Z. Ren, W. Li, R. Billinton, and W. Yan, “Probabilistic power flow analysis based on the stochastic response surface method,” _IEEE Transactions on Power Systems_ , vol. 31, no. 3, pp. 2307–2315, 2016.