<jupyter_start><jupyter_text>Landau-Zener transitionsThis notebook illustrates using the `toy_systems` package to simulate the classic Landau-Zener system that illustrates the behaviour of quantum states at avoided energy crossings. The LZ system is a 2-level system with a coupling of constant strength $\Omega$ between the states, and an energy difference that varies linearly over time from $\Delta E\rightarrow-\infty$ as $t\rightarrow-\infty$ thorugh $\Delta E = 0$ at $t = 0$ to $\Delta E\rightarrow\infty$ as $t\rightarrow\infty$. If we take the initial state of the system at $t \rightarrow -\infty$ to be $|\psi\rangle = |g\rangle$, the Landau-Zener formula gives the probability of an ImportsStart by importing the necessary packages<jupyter_code>%load_ext autoreload
%autoreload 2

import joblib
import matplotlib.pyplot as plt
plt.style.use("ggplot")
import numpy as np
import qutip
from sympy import Symbol

from toy_systems.couplings import ToyCoupling, ToyEnergy
from toy_systems.decays import ToyDecay
from toy_systems.hamiltonian import Hamiltonian
from toy_systems.quantum_system import QuantumSystem
from toy_systems.states import Basis, BasisState, ToyQuantumNumbers
from toy_systems.visualization import Visualizer<jupyter_output>The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload<jupyter_text>Define statesWe start by defining the two states of the system: the ground state $|g\rangle$ and the excited $|e\rangle$, which we will later set to have a decay path to $|g\rangle$:<jupyter_code>g = BasisState(qn=ToyQuantumNumbers(label="g"))
e = BasisState(qn=ToyQuantumNumbers(label="e"))

# Print the states as a check
print(f"|g> = {g}")
print(f"|e> = {e}")

# Check that the states are orthogonal by taking inner product
print(f"<e|g> = {e @ g}")<jupyter_output>|g> = |g>
|e> = |e>
<e|g> = 0<jupyter_text>Define basisWe use the two orthogonal states to define a basis as shown below. The indexing of the matrix and vector representations of operators is determined by the order of the states when generating the state<jupyter_code>basis = Basis((g,e))
basis.print()<jupyter_output>|0> = |g>
|1> = |e><jupyter_text>Define energiesI'm going to define the system in the rotating frame so that the Hamiltonian doesn't have any quickly rotating terms of the form $e^{i\omega t}$. The detuning between the driving field and the natural frequency of the system is taken to be $\Delta(t)$ (time-varying), so I can set the energies of the two states to be $E_g = -\Delta/2$ and $E_e = \Delta/2$.Note: the time dependence will be multiplied by the value of $\Delta$ when time-evolving the Hamiltonian.<jupyter_code>Δ = Symbol('Delta') # Using sympy to generate a symbolic Hamiltonian first (to make a "Δ" type "\Delta", then press tab)
E_g = ToyEnergy([g], -Δ/2, time_dep = "t") 
E_e = ToyEnergy([e], +Δ/2, time_dep = "t") # Note the linear time-dependence
E_g.time_dep<jupyter_output><empty_output><jupyter_text>Define couplingI'm taking the strength of the driving field to vary as a Gaussian over time. This time dependence of the drive field may result in the probability of state transfer durin RAP being slightly different from<jupyter_code>Ω = Symbol('Omega') 
coupling = ToyCoupling(g, e, Ω)<jupyter_output><empty_output><jupyter_text>Define HamiltonianNow that we have the states, energies and couplings, we can define the Hamiltonian<jupyter_code>H = Hamiltonian((E_g, E_e, coupling), basis)
display(H)<jupyter_output><empty_output><jupyter_text>Define decay matrix from |e> to |g>Can add a spontaneous decay from the excited state to the ground state as shown below:<jupyter_code>decay = ToyDecay(e, g, Symbol("Gamma"))<jupyter_output><empty_output><jupyter_text>Define a QuantumSystemThe QuantumSystem object combines the basis, Hamiltonian and decays to make setting parameters for time evolution using QuTiP more convenient.<jupyter_code># Define the system
system = QuantumSystem(
    basis=basis,
    couplings=[E_g, E_e, coupling],
    decays=[decay],
)
visualizer = Visualizer(system, vertical={"label":10}, horizontal={"label":0})

# Get representations of the Hamiltonian and the decays that will be accepted by qutip
Hqobj, c_qobj = system.get_qobjs()
Hqobj.ops[2].coeff<jupyter_output><empty_output><jupyter_text>Time-evolution using `QuTiP`The `toy_systems` package is designed to take care of generating the Hamiltonian and collapse operators for the quantum system, and spit out a `qutip.QobjEVo`-object. The actual time-evolution is then taken care of by `qutip` (or something else if the user so desires).<jupyter_code># Get a pointer to the time-evolution arguments
args = Hqobj.args
print("Keys for setting arguments:")
print(f"args = {args}")
# Generate a Qobj representing the initial state
psi0 = (1*g).qobj(basis)

# Make operators for getting the probability of being in each state
P_g_op = qutip.Qobj((1*g).density_matrix(basis), type = "oper")
P_e_op = qutip.Qobj((1*e).density_matrix(basis), type = "oper")

# Set the parameters for the system
args["Delta"] = 1e-2 # This is actually the rate at which delta varies, i.e. dD/dt 
args["sigma"] = 2/args["Delta"] # Width of the Gaussian describing time-dependence of drive field
args["Omega"] = 1 # The maximum value of the drive field Rabi rate
args["Gamma"] = 0 


# Times at which result is requested
times = np.linspace(-40/args["Delta"], 40/args["Delta"], 40001) 

# The solver used by qutip sometimes doesn't seem to notice that adiabatic passage is supposed to be happening, so reduce max step size to 
# mitigate this 
options = qutip.solver.Options(method = 'adams', nsteps=1000000, max_step=1e2/args["Delta"],)

# Run the time-evolution
result = qutip.mesolve(Hqobj, psi0, times, c_ops = c_qobj, e_ops = [P_g_op, P_e_op], options = options)<jupyter_output><empty_output><jupyter_text>Plot the result:<jupyter_code>fig, (ax,ax2) = plt.subplots(2,1, figsize = (16,9))
ax.plot(times, result.expect[0], label = "P_g")
ax.plot(times, result.expect[1], label = "P_e")
ax.legend()
ax.set_xlabel("Time / (1/Ω)")
ax.set_ylabel("Population in each state")

# Plot the coupling strength and detuning also
ax2c = ax2.twinx()
# ax2.grid(False)
# ax2c.grid(False)
ln1 = ax2.plot(times, E_e.eval(times, args), label = "Δ")
ln2 = coupling.plot_time_dep(times, args, ax = ax2c, label = 'Ω', c = 'tab:blue')
ax2c.set_ylabel('')
ax2.legend(ln1+ln2, [l.get_label() for l in (ln1 + ln2)])<jupyter_output><empty_output><jupyter_text>Testing if Landau-Zener formula appliesThe famous Landau-Zener formula for the probability of adiabatic transitions assumes a constant coupling strength and a linear detuning ramp. Under these conditions the probability of staying in the initial state can be solved for exactly and is given by $P_{ini} = 1-\exp\left(-2\pi\Omega^2/\left(\frac{dD}{dt}\right)\right)$.To investigate, I will do a parameter scan with the detuning rate. This will also illustrate using QuTiP and parallelization.<jupyter_code># Define the function we want to run in parallel to perform the parameter scan
def parallel_func(Δ: float):
        
    # Generate a Qobj representing the initial state
    psi0 = (1*g).qobj(basis)

    # Make operators for getting the probability of being in each state
    P_g_op = qutip.Qobj((1*g).density_matrix(basis), type = "oper")
    P_e_op = qutip.Qobj((1*e).density_matrix(basis), type = "oper")
    
    args = {}
    
    # Set the parameters for the system
    args["Delta"] = Δ # This is actually the rate at which delta varies, i.e. dD/dt 
    args["sigma"] = 2/Δ # Width of the Gaussian describing time-dependence of drive field
    args["Omega"] = 1 # The maximum value of the drive field Rabi rate
    args["Gamma"] = 0 


    # Times at which result is requested
    times = np.linspace(-500/args["Delta"], 500/args["Delta"], 10001) 

    # The solver used by qutip sometimes doesn't seem to notice that adiabatic passage is supposed to be happening, so reduce max step size to 
    # mitigate this 
    options = qutip.solver.Options(method = 'adams', nsteps=1000000, max_step=1e2/Δ, rhs_reuse=True)

    # Run the time-evolution
    result = qutip.mesolve(Hqobj, psi0, times, args = args, c_ops = c_qobj, e_ops = [P_g_op, P_e_op], options = options)
    
    # Plot the result
    fig, ax = plt.subplots(figsize = (16,9))
    ax.plot(times, result.expect[0], label = "P_g")
    ax.plot(times, result.expect[1], label = "P_e")
    ax.legend()
    ax.set_xlabel("Time / (1/Ω)")
    ax.set_ylabel("Population in each state")
    ax.set_title(f"Δ = {Δ}")
        
    return np.mean(result.expect[0][-1000:])

# Define the values of dD/dt we want to try
Δs = 10**np.linspace(-0.5,3,30)

# Define a loop that will be excecuted in parallel using joblib
probs = (joblib.Parallel(n_jobs = 10, backend='loky', verbose = True)(joblib.delayed(parallel_func)(Δ) for Δ in Δs))
# Plot the result
fig, (ax, ax2) = plt.subplots(2,1, figsize = (16,9))
ax.plot(np.log10(Δs), probs, 'x')
ax.set_xlabel("dΔ/dt")
ax.set_ylabel("P_g")
ax2.plot(1/Δs, np.log(np.array(probs)), 'x')
ax2.set_xlabel("1/(dΔ/dt)")
ax2.set_ylabel("ln(P_g)")

# Fit a line to the log plot
from scipy.stats import linregress
x_fit = 1/Δs[3:]
y_fit = np.log(np.array(probs[3:]))
fit_result = linregress(x_fit, y_fit)
ax2.plot(x_fit, fit_result.slope*x_fit+fit_result.intercept)

print(f"slope: ({fit_result.slope/(2*np.pi):.4f})*2*pi, intercept: {fit_result.intercept:.4f}")<jupyter_output>slope: (-0.9990)*2*pi, intercept: -0.0008
{'fraction_non_alphanumeric': 0.11476276022403044, 'fraction_numerical': 0.015745535242523513, 'mean_word_length': 5.0127064803049555, 'pattern_counts': {'":': 2, '<': 43, '<?xml version=': 0, '>': 55, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '622258', 'license': 'MIT', 'n_tokens_mistral': 3291, 'n_tokens_neox': 3101, 'n_words': 1058, 'path': 'examples/Landau-Zener transitions.ipynb', 'repo_name': 'otimgren/toy-systems', 'size': '9463'}