---
jupytext:
  formats: ipynb,md:myst
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.14.0
kernelspec:
  display_name: Python 3 (phys-555)
  language: python
  name: phys-555
---

# Universal Quantum Computing

```{code-cell} ipython3
:tags: [hide-cell]

import mmf_setup; mmf_setup.nbinit()
```

## Single Qubit Operations

Here you will demonstrate how to physically implement any single qubit gate by
manipulating an external magnetic field $\vec{B}(t)$, using a spin-½ particle as the
qubit.  For details about the gates, and various relationships, see §4.2 in the text
{cite:p}`Nielsen:2010`.

Consider a spin-½ qubit.  The Hamiltonian for such a particle in a magnetic field
$\vec{B}$ is:

\begin{gather*}
    \op{H} = -\mu \vec{B}\cdot \vec{\op{S}}, \qquad
    \op{S}_{i} = \frac{\hbar}{2}\mat{σ}_{i},
\end{gather*}

where $\mu$ is the magnetic dipole moment, and $\vec{\op{S}}$ are the spin operators,
which can be expressed in terms of the Pauli matrices

\begin{gather*}
  \mat{σ}_x = \begin{pmatrix}
    0 & 1\\
    1 & 0
  \end{pmatrix}, \qquad
  \mat{σ}_x = \begin{pmatrix}
    0 & -\I\\
    \I & 0
  \end{pmatrix}, \qquad
  \mat{σ}_z = \begin{pmatrix}
    1 & 0\\
    0 & -1
  \end{pmatrix}
\end{gather*}

Intuitively, think of $\vec{S}$ as the classical dipole moment of the spin-½ particle.
The Hamiltonian $\op{H}$ represents the "energy" of the particle in the external field,
which "wants" to align the spin.  Thus, the lowest energy state should have $\vec{B}$
parallel to the spin $\vec{S}$, maximizing $\vec{B} \cdot \vec{S}$, hence the minus
sign.  The operators $\vec{\op{S}}$ provide the quantum mechanical realization of the
classical spin, and have the given form in the standard basis for a spin-½ particle.

Explicitly, in the standard basis:

\begin{gather*}
  \mat{H} = \frac{-\mu\hbar}{2}\begin{pmatrix}
    B_z & B_x - \I B_y\\
    B_x + \I B_y & -B_z
  \end{pmatrix}.
\end{gather*}

For the purposes of most of this assignment, we will express the magnetic field $\vec{b} = -\mu
  \vec{B} /2$ and set $\hbar = 1$ so that

\begin{gather*}
  \mat{H} \equiv \frac{\mat{H}}{\hbar} = \vec{b}\cdot\vec{\mat{σ}}
   = b_x(t)\mat{σ}_x + b_y(t)\mat{σ}_y + b_z(t)\mat{σ}_z\\
   = \begin{pmatrix}
    b_z & b_x - \I b_y\\
    b_x + \I b_y & -b_z
  \end{pmatrix}.
\end{gather*}

1. Describe how to implement the rotation gate $R_{\vec{n}}(\theta)$ by angle $\theta$
   about the $\uvec{n} = (n_x, n_y, n_z)$ axis (Eq. (4.8) in the text):
   
   \begin{gather*}
     R_{\uvec{n}}(\theta) = \exp\left(
       -\I \theta \uvec{n}\cdot \frac{\vec{\mat{σ}}}{2}
     \right)
   \end{gather*}
   
   *Hint: Recall that the propagator $\mat{U} = \exp(\mat{H} t /\I\hbar)$.*
   
   I.e. Specify how $\vec{b}(t)$ should depend on time.  Assume that the magnitude $b =
   \norm{\vec{b}(t)}$ is fixed: how long we evolve to implement such a gate. 
   
2. Describe how to use an external magnetic field $\vec{b}(t)$ to implement the Pauli gates:

   \begin{gather*}
     X = \mat{σ}_x, \qquad
     Y = \mat{σ}_y, \qquad
     Z = \mat{σ}_z.
   \end{gather*}
   
3. As discussed in §4.2, exercise 4.8, Eq. (4.9), an arbitrary unitary operation $U$ can
   be expressed as
   
   \begin{gather*}
     U = e^{\I\alpha}R_{\uvec{n}}(\theta).
   \end{gather*}
   
   Prove this, and show that the phase $\alpha$ can be implemented by providing a
   constant energy offset to the Hamiltonian
   
   \begin{gather*}
     \op{H} = E_0\mat{1} + \vec{b} \cdot \vec{\op{σ}}.
   \end{gather*}

4. Describe how to use $E_0$ and an external magnetic field to implement the Hadamard
   gate $H$, phase gate $S$, and $\pi/8$ gate $T$:
   
   \begin{align*}
     H &= \frac{1}{\sqrt{2}} 
     \begin{pmatrix}
       1 & 1\\
       1 & -1
     \end{pmatrix}, \\
     S &= 
     \begin{pmatrix}
       1 & 0\\
       0 & \I
     \end{pmatrix}, \\
     T &= 
     \begin{pmatrix}
       1 & 0\\
       0 & \exp(\I\pi/4)
     \end{pmatrix}.
   \end{align*}

4. Describe how to implement an arbitrary single-qubit gate $U$, given the matrix $\mat{U}$.

5. As implied by Eq. (4.9), all of these gates can be implemented by applying a
   constant magnetic field $\vec{b}$ and energy offset $E_0$ for a fixed period of time
   $t$.  Suppose you can only generate a magnetic field of maximum strength $b \leq 
   b_{\max}$.  This means that implementing the various gates this way will take a
   certain minimum time.  Which gate takes the longest, and how long does this take?

   Restore units, and figure out what this timescale is for current experiments, asking
   members of the department (Brian Saam and Peter Engels are good sources) what types
   of magnetic fields $B_{\max}$ they can produce, and what typical magnetic moments
   $\mu$ are for the species they work with.

:::{margin}
A solution to this problem is described in {cite:p}`Barnes:2013`.
:::
7. (Hard): Can you implement these gates faster by allowing $E_0(t)$ and $\vec{b}(t)$ to
   depend on time, respecting the condition $\norm{\vec{b}} \leq b_\max$?  *(For
   simplicity, ignore the phase $\exp(\I\alpha)$ and just implement
   $R_{\uvec{n}}(\theta)$ for one of the most difficult cases.)*  This is sometimes
   referred to as the **quantum speed limit** and represents an engineering challenge
   for real systems.

Test your answers using the following function, which numerically integrates the
Schrödinger equation using your specified function $B(t)$:

```{code-cell} ipython3
import numpy as np
import scipy.linalg
from scipy.integrate import solve_ivp
sp = scipy

from phys_555_2022.utils import sigmas


def get_t_max_E0_B(t, mu=1.0, **kw):
    """Return (t_max, E0, (Bx, By, Bz)).
    
    Arguments
    ---------
    t : float
        Time 0 <= t <= t_max.
    mu : float
        Constant magnetic moment.
    hbar : float
        Constant hbar.        
    kw : dict
        Additional arguments relevant for the problem.
    
    Returns
    -------
    t_max : float
        Total evolution time.
    E0 : float
        Constant energy offset to define overall phase.
    Bx, By, Bz : float
        Magnetic field components at time t.
    """
    t_max = 1.0
    E0 = 0.0
    B = (0.0, 0.0, 0.0)
    return (t_max, E0, B)


def get_U(get_B, mu=1.0, hbar=1.0, **kw):
    """Return the propagator U by solving the SEQ.
    
    This implementation directly integrates the SEQ using
    the identity matrix as the initial state.
    """
    t_max, E0, B0 = get_B(t=0, mu=mu, hbar=hbar, **kw)
    
    def dy_dt(t, y):
        """Return dy/dt"""
        U = y.reshape(U0.shape)
        t_max, E0, B = get_B(t=t, mu=mu, hbar=hbar, **kw)
        Ss = sigmas * hbar / 2
        H = -mu*np.einsum('a,a...', B, Ss)
        dU_dt = (H @ U + E0 * U) / 1j / hbar
        dy_dt = dU_dt.ravel()
        return dy_dt
    
    U0 = np.eye(2, dtype=complex)
    res = solve_ivp(dy_dt, t_span=(0, t_max), y0=U0.ravel(),
                    rtol=1e-8, atol=1e-8)
    if not res.success:
        raise ValueError(res.message)
    return res.y.T[-1].reshape(U0.shape)  
```

```{code-cell} ipython3
# Sample solution for R_theta gate.
# Here we let theta == theta*n be the full vector.
rng = np.random.default_rng(seed=1)

def B_R(t, theta, mu=1.0, hbar=1.0):
    """Return B(t) implementing R(theta)."""
    t_max = 1 / mu
    E0 = 0.0
    B = -theta
    return (t_max, E0, B)
    

mu, hbar = rng.random(size=2)
theta = rng.normal(size=3)

R_exact = sp.linalg.expm(-1j*np.einsum('a,a...', theta, sigmas/2))
U = get_U(B_R, theta=theta, mu=mu, hbar=hbar)
assert np.allclose(R_exact, U)
```

```{code-cell} ipython3
def B_X(t, mu=1.0, hbar=1.0):
    """Return (t_max, E0, B(t)) implementing X."""
    t_max = 1.0
    E0 = 0
    B = (0, 0, 0)
    return t_max, E0, B

def B_Y(t, mu=1.0, hbar=1.0):
    """Return (t_max, E0, B(t)) implementing Y."""
    t_max = 1.0
    E0 = 0
    B = (0, 0, 0)
    return t_max, E0, B

def B_Z(t, mu=1.0, hbar=1.0):
    """Return (t_max, E0, B(t)) implementing Z."""
    t_max = 1.0
    E0 = 0
    B = (0, 0, 0)
    return t_max, E0, B

def B_H(t, mu=1.0, hbar=1.0):
    """Return (t_max, E0, B(t)) implementing H."""
    t_max = 1.0
    E0 = 0
    B = (0, 0, 0)
    return t_max, E0, B

def B_S(t, mu=1.0, hbar=1.0):
    """Return (t_max, E0, B(t)) implementing S."""
    t_max = 1.0
    E0 = 0
    B = (0, 0, 0)
    return t_max, E0, B

def B_T(t, mu=1.0, hbar=1.0):
    """Return (t_max, E0, B(t)) implementing T."""
    t_max = 1.0
    E0 = 0
    B = (0, 0, 0)
    return t_max, E0, B

mu, hbar = rng.random(size=2)
X, Y, Z = sigmas
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
S = np.diag([1, 1j])
T = np.diag([1, np.exp(np.pi * 0.25j)])

assert np.allclose(X, get_U(B_X, mu=mu, hbar=hbar))
assert np.allclose(Y, get_U(B_Y, mu=mu, hbar=hbar))
assert np.allclose(Z, get_U(B_Z, mu=mu, hbar=hbar))
assert np.allclose(H, get_U(B_H, mu=mu, hbar=hbar))
assert np.allclose(S, get_U(B_S, mu=mu, hbar=hbar))
assert np.allclose(T, get_U(B_T, mu=mu, hbar=hbar))
```

```{code-cell} ipython3
def B_U(t, U, mu=1.0, hbar=1.0):
    """Return (t_max, E0, B(t)) implementing U."""
    t_max = 1.0
    E0 = 0
    B = (0, 0, 0)
    return t_max, E0, B
    
# Random unitary matrix as an exponential of a random anti-hermitian matrix
A = rng.normal(size=(2, 2*2)).view(dtype=complex)
U = sp.linalg.expm(A - A.T.conj())
assert np.allclose(U.T.conj() @ U, np.eye(2))

assert np.allclose(U, get_U(B_U, U=U, mu=mu, hbar=hbar))
```

```{code-cell} ipython3

```