4.6. Floquet optimization: non-equilibrium steady states (NESSs)

This tutorial demonstrates a Floquet optimization on a simple three-level quantum model of the Nitrogen-vancancy center in diamond, treated as an open system. The goal is to find a periodic driving that leads to a non-equilibrium steady state (NESS) that is optimal in some given way: in this case, we will define this optimaility in terms of the time-averaged expectation value of some operator.

The model is defined by the three-level Hamiltonian:

\begin{eqnarray} H(u, t) &=& H_0 + V(u, t), \\ H_0 &=& -B_s S_z + N_z S_z^2 + N_{xy}(S_x^2-S_y^2), \\ \label{eq:tdpart} V(u, t) &=& - g_x(u, t) B_d S_x - g_y(u, t) B_d S_y. \end{eqnarray}

This model is taken from [Ikeda et al, Science Advances 6, eabb4019 (2020)]. \(S_x, S_y\) and \(S_z\) are the spin operators, whereas \(N_z\), \(N_{xy}\), \(B_s\), and \(B_d\) are constants. The shape of the real time-dependent control functions \(g_x\) and \(g_y\), dependent on the control parameters \(u\), is explained below.

The system is governed by Lindblad’s equation:

\begin{align} \dot{\rho}(t) = -i\left[H(u, t), \rho(t)\right] + \sum_{ij} \gamma_{ij} \left( V_{ij}\rho(t)V^\dagger_{ij} - \frac{1}{2} \lbrace V_{ij}^\dagger V_{ij}, \rho(t) \rbrace\right)\,. \label{eq:lindblad0-eq} \end{align}

The transition operators are \(V_{ij} = \vert E_i\rangle\langle E_j\vert\), and the dissipative constants are \(\gamma_{ij} = \gamma e^{-\beta E_i} / (e^{-\beta E_i}+e^{-\beta E_j})\) and \(\gamma_{ii}=0\), where \(\beta = 1/(k_{\rm B}T)\) is the inverse of the temperature, and \(\gamma\) is a rate constant. Notice that this dissipation model ensures the detailed balance condition, \(\gamma_{ij}e^{-\beta E_j} = \gamma_{ji}e^{-\beta E_i}\).

import os
import numpy as np
import qutip as qt
import nlopt
import matplotlib as mpl
import matplotlib.pyplot as plt

import qocttools
import qocttools.pulses as pulses
import qocttools.qoct as qoct
import qocttools.hamiltonians as hamiltonians
import qocttools.floquet as floquet
import qocttools.target as target

It is good practice to print the precise version of the software that you are using.

qocttools.about()

No git dir found

Running qocttools version 0.0.14

data = []

Now, we build the static Hamiltonian \(H_0\) (stored into the Qobj object H0), and the two coupling operators \(V_1 = -B_d S_x\) and \(V_2 = -B_d S_y\) (stored into the Qobj objects V[0] and V[1]). The function system_definition also returns some Lindblad operators into the A list, as defined above. The field-free eigenvalues are stored stored in the array e and the eigenfunctions in psi.

Sx = qt.jmat(1, "x")
Sy = qt.jmat(1, "y")
Sz = qt.jmat(1, "z")
Bs = 0.3
Nz = 1.00
Nxy = 0.05
Bd = 0.1
omega = 1.00
gamma = 0.2
beta = 3.0
d = 3
dim = d**2

def system_definition():
    H0 = -Bs * Sz + Nz * Sz**2 + Nxy * (Sx**2 - Sy**2)
    Vx = -Bd * Sx
    Vy = -Bd * Sy
    A = []
    e, psi = H0.eigenstates()
    for i in range(d):
        for j in range(d):
            if j == i:
                continue
            gammaij = gamma * np.exp(-beta*e[j]) / (np.exp(-beta*e[i])+np.exp(-beta*e[j]))
            A.append( np.sqrt(gammaij) * psi[j] * psi[i].dag())
    return H0, [Vx, Vy], A, e, psi

H0, V, A, e, psi = system_definition()

print("Field-free eigenvalues = {}".format(e))

Field-free eigenvalues = [0.         0.69586187 1.30413813]

We must set a value for the period of the driving, and the corresponding frequency. That is of course problem dependent; for this tutorial we will arbitrarily use \(\omega = 0.5\).

omega0 = 0.5
T = (2.0*np.pi/omega0)
nts = 100
times = np.linspace(0, T, nts + 1)

We now build the hamiltonian object.

Ham = hamiltonians.hamiltonian(H0, V, A)

The target operator is \(S_z\). Thus, the goal would be to find the control parameters that lead to a NESS that maximizes the time-average of \(\langle S_z\rangle\).

target_operator = Sz

These functions are used to to create the pulses.

def pulse_definition(T, p, bound = 4.0, seed = 0):
    if seed >= 0:
        np.random.seed(seed)
        u = (bound-(-bound)) * np.random.random_sample(p) + (-bound)
        g1 = pulses.pulse("fourier", T, u)
        u = (bound-(-bound)) * np.random.random_sample(p) + (-bound)
        g2 = pulses.pulse("fourier", T, u)
    else:
        #M = p
        K = 1
        u = np.zeros(p)
        u[K] = bound #np.sqrt(T)/2
        g1 = pulses.pulse("fourier", T, u)
        u = np.zeros(p)
        u[K+1] = bound #np.sqrt(T)/2
        g2 = pulses.pulse("fourier", T, u)
    return [g1, g2]


def pulse_set_new(g, bound = 4.0, seed = 0):
    np.random.seed(seed)
    p = g[0].u.shape[0]
    u = (bound-(-bound)) * np.random.random_sample(p) + (-bound)
    g[0].set_parameters(u)
    u = (bound-(-bound)) * np.random.random_sample(p) + (-bound)
    g[1].set_parameters(u)

M = 4
bound = 4.0
g = pulse_definition(T, 2*M+1, bound = bound, seed = -1)
gref = pulse_definition(T, 2*M+1, bound = bound, seed = -1)
u = pulses.pulse_collection_get_parameters(g)
pulses.pulse_collection_set_parameters(gref, u)

tg = target.Target('floquet', operator = target_operator, T = T)
opt = qoct.Qoct(Ham, T, nts, tg, g, None, floquet_mode = 'ness')

print("G(u) = {}".format(opt.gfunc(u)))
print("G(u=0) = {}".format(opt.gfunc(np.zeros_like(u))))

G(u) = 0.0907525770787369
G(u=0) = 0.08966865961492401

We will not do the optimization, because we want the tutorial to run quickly. But let us compute the gradient of the target function, and check it against a finite-difference formula.

check_gradient = True
if check_gradient:
    derqoct, dernum, error, elapsed_time = opt.check_grad(u)
    print("QOCT calculation: \t{}".format(derqoct))
    print("Ridders calculation: \t{} +- {}".format(dernum, error))

QOCT calculation:       0.0021752268449753524
Ridders calculation:    0.0021752268447795702 +- 8.762131645245752e-14

data.append(derqoct)

# This file is used by the testing script of the code.
with open("data", "w") as f:
    for i in data:
        f.write("{:.14e}\n".format(i))