How to define continuous analytical Hamiltonians

import matplotlib.pyplot as plt
import numpy as np
import qctrlvisualizer as qv
import boulderopal as bo

plt.style.use(qv.get_qctrl_style())

# Create the data flow graph describing the system.
graph = bo.Graph()

# Define the duration of the signals.
duration = 2e-6  # s

# Define the number of samples for the discretizations.
sample_count = 512
sample_times = np.linspace(0.0, duration, sample_count)

# Define maximum values for the signals.
omega_max = 1.5e6  # rad/s
delta_max = 2.1e6  # rad/s

# Define the signals.
t = graph.identity_stf()
omega = omega_max * graph.sin(np.pi * t / duration) ** 2
delta = -delta_max * graph.cos(np.pi * t / duration)

# Construct the Hamiltonian.
hamiltonian = omega * graph.pauli_matrix("X") + delta * graph.pauli_matrix("Z")

# Compute the time-evolution unitaries.
unitaries = graph.time_evolution_operators_stf(
    hamiltonian=hamiltonian, sample_times=sample_times
)

# Define the initial state.
initial_state = graph.fock_state(2, 0)[:, None]

# Evolve the initial state with the unitary to get the final state.
evolved_states = unitaries @ initial_state
evolved_states.name = "states"

# Discretize the signals as STFs can't be directly exported.
graph.discretize_stf(omega, duration, sample_count, name=r"$\Omega$")
graph.discretize_stf(delta, duration, sample_count, name=r"$\Delta$")

# Run simulation.
result = bo.execute_graph(
    graph=graph, output_node_names=["states", r"$\Omega$", r"$\Delta$"]
)

# Plot the analytical signals used.
qv.plot_controls(
    {name: result["output"][name] for name in [r"$\Omega$", r"$\Delta$"]}, smooth=True
)
plt.suptitle("Landau–Zener signals")

# Plot the population dynamics.
populations = np.abs(result["output"]["states"]["value"].squeeze()) ** 2
qv.plot_population_dynamics(
    sample_times, {rf"$|{k}\rangle$": populations[:, k] for k in range(2)}
)
plt.suptitle("Landau–Zener transition dynamics")
plt.show()

graph = bo.Graph()

# Define the duration of the signals.
duration = 3e-6  # s
sample_count = 64

# Define α(t) as an arctan ramp.
t = graph.identity_stf()
alpha_ramp = duration / 8
alpha_zero = graph.optimizable_scalar(
    lower_bound=-3e7, upper_bound=3e7, name="alpha_amplitude"
)
alpha_center = graph.optimizable_scalar(
    lower_bound=0.0, upper_bound=duration, name="alpha_center"
)
alpha_stf = alpha_zero * graph.arctan((t - alpha_center) / alpha_ramp)

# Define γ(t) as a Gaussian pulse.
gamma_max = 2.0 * np.pi * 1e6
gamma_width = duration / 8
gamma_modulus = graph.optimizable_scalar(lower_bound=-gamma_max, upper_bound=gamma_max)
gamma_phase = graph.optimizable_scalar(
    lower_bound=0,
    upper_bound=2.0 * np.pi,
    is_lower_unbounded=False,
    is_upper_unbounded=False,
)
gamma_amplitude = gamma_modulus * graph.exp(1j * gamma_phase)
gamma_amplitude.name = "gamma_amplitude"

gamma_center_time = graph.optimizable_scalar(
    lower_bound=0.0, upper_bound=duration, name="gamma_center"
)
gamma_stf = graph.signals.gaussian_pulse_stf(
    gamma_amplitude, gamma_width, gamma_center_time
)

# Discretize the signals.
alpha = graph.discretize_stf(alpha_stf, duration, sample_count, name=r"$\alpha$")
gamma = graph.discretize_stf(gamma_stf, duration, sample_count, name=r"$\gamma$")

# Create Hamiltonian.
hamiltonian = alpha * graph.pauli_matrix("Z") + graph.hermitian_part(
    gamma * graph.pauli_matrix("M")
)

# Add detuning term.
delta = 2 * np.pi * 0.25e6  # rad/s
hamiltonian += delta * graph.pauli_matrix("Z")

# Define target gate.
target = graph.target(operator=graph.pauli_matrix("X"))

# Add dephasing noise term.
zeta = 2 * np.pi * 20e3  # rad/s
dephasing = zeta * graph.pauli_matrix("Z")

# Robust infidelity.
robust_infidelity = graph.infidelity_pwc(
    hamiltonian=hamiltonian,
    noise_operators=[dephasing],
    target=target,
    name="robust_infidelity",
)

result = bo.run_optimization(
    graph=graph,
    cost_node_name="robust_infidelity",
    output_node_names=[
        "alpha_amplitude",
        "alpha_center",
        "gamma_amplitude",
        "gamma_center",
        r"$\alpha$",
        r"$\gamma$",
    ],
    optimization_count=8,
)

print(f"Optimized cost: {result['cost']:.3e}")

# Print out the converged signal parameters.
print("Optimal signal parameters")
print(f"\tα(t) amplitude:\t{result['output']['alpha_amplitude']['value']:.3e} rad/s")
print(f"\tα(t) center:\t{result['output']['alpha_center']['value']:.3e} s")
print(f"\tγ(t) amplitude:\t{result['output']['gamma_amplitude']['value']:.3e} rad/s")
print(f"\tγ(t) center:\t{result['output']['gamma_center']['value']:.3e} s")

# Plot the optimal analytical signals found.
qv.plot_controls({name: result["output"][name] for name in [r"$\alpha$", r"$\gamma$"]})
plt.suptitle("Optimal signals")
plt.show()