Universal Quantum Computing#

Hide code cell content
import mmf_setup; mmf_setup.nbinit()

This cell adds /home/docs/checkouts/readthedocs.org/user_builds/physics-555-quantum-technologies/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:

  • Choose "Trust Notebook" from the "File" menu.
  • Re-execute this cell.
  • Reload the notebook.

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 [Nielsen and Chuang, 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*}\]
  5. Describe how to implement an arbitrary single-qubit gate \(U\), given the matrix \(\mat{U}\).

  6. 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.

  1. (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)\):

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)  
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[2], line 6
      3 from scipy.integrate import solve_ivp
      4 sp = scipy
----> 6 from phys_555_2022.utils import sigmas
      9 def get_t_max_E0_B(t, mu=1.0, **kw):
     10     """Return (t_max, E0, (Bx, By, Bz)).
     11     
     12     Arguments
   (...)
     30         Magnetic field components at time t.
     31     """

ModuleNotFoundError: No module named 'phys_555_2022'
# 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)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 16
     13 mu, hbar = rng.random(size=2)
     14 theta = rng.normal(size=3)
---> 16 R_exact = sp.linalg.expm(-1j*np.einsum('a,a...', theta, sigmas/2))
     17 U = get_U(B_R, theta=theta, mu=mu, hbar=hbar)
     18 assert np.allclose(R_exact, U)

NameError: name 'sigmas' is not defined
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))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 44
     41     return t_max, E0, B
     43 mu, hbar = rng.random(size=2)
---> 44 X, Y, Z = sigmas
     45 H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
     46 S = np.diag([1, 1j])

NameError: name 'sigmas' is not defined
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))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 13
     10 U = sp.linalg.expm(A - A.T.conj())
     11 assert np.allclose(U.T.conj() @ U, np.eye(2))
---> 13 assert np.allclose(U, get_U(B_U, U=U, mu=mu, hbar=hbar))

NameError: name 'get_U' is not defined