Rabi Oscillator

In this example we set up a Rabi oscillator problem, and optimize for a Pauli-X gate.

The Rabi oscillator consists of a single qubit in the rotating frame, with the rotation frequency chosen perfectly so that the state vector is constant-in-time unless a control pulse is applied.

We have one control pulse available in the lab frame (and therefore two in the rotating frame), which will have a constant amplitude in the rotating frame.

Setting up the Problem

First we construct the Hamiltonians and initial conditions, and put them together in a SchrodingerProb. For a Rabi oscillator the system/drift Hamiltonian is zero, and the control Hamiltonian has real part $a+a^\dagger$ and imaginary part $a-a^\dagger$, where $a$ is the lowering/annihilation operator

\[a = \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix}.\]

Consequently, the dynamics of the problem are

\[\frac{dU}{dt} = \left(p_0\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} +iq_0 \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix}\right)U\]

(NOTE: not sure if there should be a $-i$ on the right-hand side when we work in the complex formulation).

using QuantumGateDesign

H0_sym = zeros(2,2)
H0_asym = zeros(2,2)
a = [0.0 1;
     0   0]
Hc_sym = a + a'
Hc_asym = a - a'

sym_ops = [Hc_sym]
asym_ops = [Hc_asym]

U0_complex = [1.0 0;
              0   1]

u0 = real(U0_complex)
v0 = imag(U0_complex)

Ω = 0.5 + 0.0im

# 5 Rabi Oscillations
tf = 10pi / (2*abs(Ω))
nsteps = 10
N_ess_levels = 2

prob = SchrodingerProb(
    H0_sym, H0_asym, sym_ops, asym_ops, u0, v0, tf, nsteps, N_ess_levels
)

Setting up the Control

Next we establish the controls. We use a GRAPEControl, which implements a piecewise constant control (like what is used in the GRAPE method, with one control coefficient (per real/imaginary part) in order to achieve control pulses that are constant throughout the duration of the gate.

N_control_coeff = 1
control = QuantumGateDesign.GRAPEControl(N_control_coeff, prob.tf)

pcof = [real(Ω), imag(Ω)]

Performing the Optimization

Finally, we set up our target and interface to IPOPT (using optimize_gate)to perform the optimization.

# Pauli-X gate, or 1-qubit SWAP
target_complex = [0.0 1;
                  1   0]
target = vcat(real(target_complex), imag(target_complex))

ret_dict = optimize_gate(prob, control, pcof, target, order=4)

2 Qudits in the Dispersive limit

QuantumGateDesign.jl provides convenience functions for setting up several SchrodingerProb's which arise frequently in quantum computing. The code below uses the DispersiveProblem constructor to make a SchrodingerProb representing two qubits in the dispersive limit.

# Set up the SchrodingerProb
subsystem_sizes = (4,4)
essential_subsystem_sizes = (2,2)
transition_freqs = [4.10595, 4.81526] .* 2pi
rotation_freqs = copy(transition_freqs)
# 0.1 for artificially large fast coupling. Actual value is 1e-6
kerr_coeffs = 2pi .* [2*0.1099   0.1
                      0.1        2*0.1126]
tf = 50.0
nsteps = 100
sparse_rep=true

prob = DispersiveProblem(
    subsystem_sizes, essential_subsystem_sizes,
    transition_freqs, rotation_freqs, kerr_coeffs,
    tf, nsteps;
    sparse_rep=sparse_rep
)

# Set up the Controls
carrier_wave_freqs = transition_freqs .- rotation_freqs
bspline_degree = 2
N_knots = 5

D1 = 10
fa = sum(transition_freqs) / length(transition_freqs)
om1 = 2pi .* [transition_freqs[1] - rotation_freqs[1],
              transition_freqs[1] - rotation_freqs[1] - kerr_coeffs[1,2],
             ]
om2 = 2pi .* [transition_freqs[2] - rotation_freqs[2],
              transition_freqs[2] - rotation_freqs[2] - kerr_coeffs[1,2],
             ]

control1 =  BSplineControl(tf, D1, om1)
control2 =  BSplineControl(tf, D1, om2)

controls = [control1, control2]

N_coeff = get_number_of_control_parameters(controls)

amax = 0.04 # Keeping pulse amplitudes below 40 MHz
pcof0 = rand(N_coeff) .* 0.04 / 10

# Set up target gate
CNOT_target = create_gate(
    subsystem_sizes, essential_subsystem_sizes,
    [(1,0) => (1,1),
     (1,1) => (1,0)]
)

# Transform target into rotating frame
rot_mat = rotation_matrix(subsystem_sizes, rotation_freqs, tf)
CNOT_target_complex = real_to_complex(CNOT_target)
CNOT_target_rotating_frame_complex = prod(rot_mat) * CNOT_target_complex
CNOT_target = complex_to_real(CNOT_target_rotating_frame_complex)

# Compare stepsize to relative error in solution history
run_convergence_test = false
if run_convergence_test
    convergence_dict = QuantumGateDesign.get_histories(prob, controls, pcof0, 15, base_nsteps=100, orders=(2,4))
    pl1 = QuantumGateDesign.plot_stepsize_convergence(convergence_dict)
    pl2 = QuantumGateDesign.plot_timing_convergence(convergence_dict)
end

# Use step sizes based on output of above tests
nsteps_order2 = trunc(Int, tf/10.0^(-3.5))
nsteps_order4 = trunc(Int, tf/10.0^(-0.5))

#prob.nsteps = nsteps_order2
prob.nsteps = 100
return_dict_order2 = optimize_gate(prob, controls, pcof0, CNOT_target, order=2, maxIter=50, pcof_L=-amax, pcof_U=amax)

prob.nsteps = nsteps_order4
return_dict_order4 = optimize_gate(prob, controls, pcof0, CNOT_target, order=4, maxIter=50, pcof_L=-amax, pcof_U=amax)