Fundamentals
The full Hamiltonian of an \(n\)-qudit system with static coupling and drive terms is
\[H = H_0 + H_{\mathrm{int}} + H_{\mathrm{d}},\]
where
\[\begin{split}H_0 & = \sum_{j=1}^{n} \left[ \omega_j b_j^{\dagger} b_j + \frac{\Delta_j}{2} b_j^{\dagger} b_j (b_j^{\dagger} b_j - 1) \right]
= \sum_{j=1}^{n} \left[ \left( \omega_j - \frac{\Delta_j}{2} \right) N_j + \frac{\Delta_j}{2} N_j^2 \right], \\
H_{\mathrm{int}} & = \sum_{j<k} J_{jk} \left( b_j^{\dagger} b_k + b_j b_k^{\dagger} \right), \\
H_{\mathrm{d}} & = \sum_{jk} \alpha_{jk} \Omega_j \left( p_j(t) \cos (\nu_j t - \rho_{jk}) + q_j(t) \sin (\nu_j t - \rho_{jk}) \right)
\left( b_k^{\dagger} + b_k \right) \\
& = \sum_{jk} \alpha_{jk} \frac{\Omega_j}{2} \left( r_j(t) e^{-i(\nu_j t - \rho_{jk})} + \mathrm{c.c.} \right)
\left( b_k^{\dagger} + b_k \right)\end{split}\]
with \(b_j^{\dagger}\) and \(b_j\) the creation and annihilation operators for qudit \(j\) and
\(\omega_j\): Qubit frequency of qudit \(j\)
\(\Delta_j\): Anharmonicity of qudit \(j\)
\(J_{jk}\): Coupling between qudits \(j\) and \(k\)
\(\Omega_j\): Base amplitude of drive in channel \(j\)
\(p_j (t), q_j (t)\): I and Q components of the pulse envelope of drive in channel \(j\), and \(r_j (t) = p_j(t) + iq_j(t)\)
\(\nu_j\): Local oscillator frequency of drive in channel \(j\)
\(\alpha_{jk}\): Crosstalk attenuation factor of drive in channel \(j\) sensed by qudit \(k\)
\(\rho_{jk}\): Crosstalk phase shift of drive in channel \(j\) sensed by qudit \(k\).
When considering more than a single drive frequency per channel, it can be more convenient to express
the drive Hamiltonian in the frequency domain:
\[H_{\mathrm{d}} = \sum_{jk} \alpha_{jk} \frac{\Omega_j}{2} \int d\nu \left( \tilde{r}_j(\nu) e^{-i (\nu t - \rho_{jk})}
+ \mathrm{c.c.} \right) \left( b_k^{\dagger} + b_k \right)\]
Change of frame
Qudit frame
We move to the qudit frame through a transformation with \(U_q(t) := e^{i H_0 t}\):
\[\begin{split}\tilde{H} & := U_q H U_q^{\dagger} + i \dot{U_q} U_q^{\dagger} \\
& = U_q (H_{\mathrm{int}} + H_{\mathrm{d}}) U_q^{\dagger} =: \tilde{H}_{\mathrm{int}} + \tilde{H}_{\mathrm{d}}.\end{split}\]
\(\tilde{H}\) is the generator of time evolution for state \(|\tilde{\psi}(t)\rangle_q := U_q(t) |\psi(t)\rangle\):
\[\begin{split}i \frac{\partial}{\partial t} |\tilde{\psi}(t)\rangle_q & = (i \dot{U}_q(t) + U_q(t) H(t)) |\psi(t)\rangle \\
& = \tilde{H}(t) |\tilde{\psi}(t)\rangle_q.\end{split}\]
To write down \(\tilde{H}_{\mathrm{int}}\) and \(\tilde{H}_{\mathrm{d}}\) in terms of \(\{b_j\}_j\) and \(\{N_j\}_j\),
we first note that \(U_q\) can be factored into commuting subsystem unitaries:
\[U_q = \bigotimes_j \exp \left\{ i \left[\left( \omega_j - \frac{\Delta_j}{2} \right) N_j + \frac{\Delta_j}{2} N_j^2 \right] t \right\}
=: \bigotimes_j e^{i h_j t}.\]
Each \(h_j\) commutes with \(b_k\) and \(b_k^{\dagger}\) if \(k \neq j\), so
\[\begin{split}\tilde{H}_{\mathrm{int}} & = \sum_{j<k} J_{jk} \left( \tilde{b}_j^{\dagger} \tilde{b}_k + \tilde{b}_j \tilde{b}_k^{\dagger} \right) \\
\tilde{H}_{\mathrm{d}} & = \sum_{jk} \alpha_{jk} \frac{\Omega_j}{2} \left( r_j(t) e^{-i(\nu_j t - \rho_{jk})} + \mathrm{c.c.} \right)
\left( \tilde{b}_k^{\dagger} + \tilde{b}_k \right)\end{split}\]
where
\[\begin{split}\tilde{b}_{j} & = e^{i h_j t} b_j e^{-i h_j t}, \\
\tilde{b}_{j}^{\dagger} & = e^{i h_j t} b_j^{\dagger} e^{-i h_j t}.\end{split}\]
By definition \(b_j N_j = (N_j + 1) b_j\), which implies
\[b_j e^{-i h_j t} = \exp \left\{ -i \left[\left( \omega_j - \frac{\Delta_j}{2} \right) (N_j + 1)
+ \frac{\Delta_j}{2} (N_j + 1)^2) \right] t \right\} b_j\]
and therefore
\[\begin{split}\tilde{b}_{j} & = \exp \left\{ i \left[\left( \omega_j - \frac{\Delta_j}{2} \right) (N_j - (N_j + 1))
+ \frac{\Delta_j}{2} (N_j^2 - (N_j + 1)^2) \right] t \right\} b_j \\
& = e^{-i(\omega_j + \Delta_j N_j) t} b_j.\end{split}\]
Similarly, \(b_j^{\dagger} N_j = (N_j - 1) b_j^{\dagger}\) leads to
\[\begin{split}\tilde{b}_{j}^{\dagger} & = \exp \left\{ i \left[\left( \omega_j - \frac{\Delta_j}{2} \right) (N_j - (N_j - 1))
+ \frac{\Delta_j}{2} (N_j^2 - (N_j - 1)^2) \right] t \right\} b_j^{\dagger} \\
& = e^{i(\omega_j + \Delta_j (N_j - 1)) t} b_j^{\dagger}.\end{split}\]
The interaction Hamiltonian in the qudit frame is therefore
\[\begin{split}\tilde{H}_{\mathrm{int}} & = \sum_{j<k} J_{jk} \left( e^{i (\omega_j - \omega_k) t} e^{i [\Delta_j (N_j - 1) - \Delta_k N_k] t}
b_j^{\dagger} b_k + \mathrm{h.c.} \right) \\
& = \sum_{j<k} J_{jk} \left( e^{i (\omega_j - \omega_k) t} \sum_{lm} e^{i (\Delta_j l - \Delta_k m) t}
\sqrt{(l+1)(m+1)} | l + 1 \rangle_j \langle l |_j \otimes | m \rangle_k
\langle m + 1 |_k + \mathrm{h.c.} \right).\end{split}\]
In the last line, we used the expansion of the annihilation operator \(b_j = \sum_{l} \sqrt{l+1} | l \rangle_j \langle l + 1 |_j\)
and its Hermitian conjugate.
The drive Hamiltonian in the qudit frame is
\[\begin{split}\tilde{H}_{\mathrm{d}} & = \sum_{jk} \alpha_{jk} \frac{\Omega_j}{2} \left( r_j(t) e^{-i(\nu_j t - \rho_{jk})} + \mathrm{c.c.} \right)
\left( e^{i(\omega_k + \Delta_k (N_k - 1))t} b_k^{\dagger} + \mathrm{h.c.} \right) \\
& = \sum_{jk} \alpha_{jk} \frac{\Omega_j}{2} \left( r_j(t) e^{-i(\nu_j t - \rho_{jk})} + \mathrm{c.c.} \right)
\sum_l \left( e^{i \omega_k t} e^{i \Delta_k l t} \sqrt{l+1} | l + 1 \rangle_k \langle l |_k + \mathrm{h.c.} \right).\end{split}\]
Dressed frame
Even in the absense of a drive, \(\tilde{H}_{\mathrm{int}}\) above actually causes slow phase drifts in the qudit frame. As it is difficult
to see this from a time-dependent \(\tilde{H}_{\mathrm{int}}\), we move back once again to the lab frame and diagonalize
\(H_{\mathrm{stat}} = H_0 + H_{\mathrm{int}}\) as
\[H_{\mathrm{stat}} = V E V^{\dagger}.\]
The unitary \(V\) is chosen to be
\[V = I + \eta\]
that minimizes \(|\eta|\) while satisfying the diagonalization condition above. Given that \(H_{\mathrm{int}}\) is off-diagonal and
\(|H_{\mathrm{int}}| \ll |H_0|\), this results in
\[E = H_0 + \delta\]
for some small diagonal \(\delta\).
The time evolution by \(H_{\mathrm{stat}}\) is
\[\begin{split}e^{-i H_{\mathrm{stat}} t} & = V e^{-iEt} V^{\dagger} \\
& = e^{-iEt} + \eta e^{-iEt} + e^{-iEt} \eta^{\dagger} + \eta e^{-iEt} \eta^{\dagger}.\end{split}\]
Because \(\tilde{H}_{\mathrm{int}}\) is the generator of time evolution in the qudit frame
\[\begin{split}& T \left[ \exp \left(-i \int_{0}^{t} dt' \tilde{H}_{\mathrm{int}} (t') \right) \right] U_q(0) |\psi(0)\rangle
= U_q(t) |\psi(t)\rangle \\
& = U_q(t) e^{-i H_{\mathrm{stat}} t} |\psi(0)\rangle.\end{split}\]
Therefore, for any free-Hamiltonian eigenstate \(|l\rangle\),
\[T \left[ \exp \left(-i \int_{0}^{t} dt' \tilde{H}_{\mathrm{int}} (t') \right) \right] |l\rangle
= e^{-i \delta_{l} t} |l\rangle + e^{i H_0 t} \left(\eta e^{-iEt} + e^{-iEt} \eta^{\dagger} + \eta e^{-iEt} \eta^{\dagger} \right) |l\rangle,\]
where \(\delta_{l}\) is the \(l\)-th element of \(\delta\).
To eliminate these phase drifts, we would like to work in the frame defined by \(U_E = e^{i E t}\). However, this mathematically trivial change
of frame is not physically practical, because \(E\) does not necessarily render itself to a sum of single-qudit operators, while all drive and
readout are performed in terms of individual qudits. Therefore we fall back to a frame defined by \(U_d = e^{i D t}\), where
\[D = \sum_{j=1}^{n} \sum_{l} \left( \langle l |_j \otimes \langle 0 |^{\otimes n-1} E | l \rangle_j \otimes | 0 \rangle^{\otimes n-1} \right)
| l \rangle_j \langle l |_j,\]
which eliminates the phase drifts of single-qudit excitations. This frame rotates at “dressed” frequencies, i.e., free-qudit frequencies
shifted by the effects of inter-qudit interactions.
General frame
We can also move to an arbitrary frame specifying the frequency and phase offset for each level gap of each qudit. Let for qudit \(j\) the frequency
and the phase offset between level \(l\) and \(l+1\) be \(\xi_{j}^{l}\) and \(\phi_{j}^{l}\), and \(\Xi_{j}^{l} := \sum_{m<l} \xi_{j}^{m}\),
\(\Phi_{j}^{l} := \sum_{m<l} \phi_{j}^{m}\). Then the transformation unitary is
\[\begin{split}U_f & := \bigotimes_j \exp \left[ i \sum_l \left( \Xi_j^l t + \Phi_j^{l} \right) |l\rangle_j \langle l |_j \right] \\
& = \bigotimes_j \sum_l \exp \left[ i \left( \Xi_j^l t + \Phi_j^{l} \right) \right] |l\rangle_j \langle l |_j.\end{split}\]
\(U_f\) commutes with the free Hamiltonian \(H_0\) but \(i \dot{U}_f U_f^{\dagger} \neq -H_0\) in general, so
\[\tilde{H} = U_f H U_f^{\dagger} + i \dot{U_f} U_f^{\dagger} = H_{\mathrm{diag}} + \tilde{H}_{\mathrm{int}} + \tilde{H}_{\mathrm{d}}.\]
The three terms can be expressed in terms of individual qudit levels as
\[\begin{split}H_{\mathrm{diag}} & = \sum_{j} \sum_{l} \left[ \left( \omega_j - \frac{\Delta_j}{2} \right) l + \frac{\Delta_j}{2} l^2 - \Xi_j^{l} \right] |l\rangle_j \langle l|_j, \\
\tilde{H}_{\mathrm{int}} & = \sum_{j<k} J_{jk} \sum_{lm} \left( e^{i [(\xi_j^{l} - \xi_{k}^{m}) t + (\phi_j^{l} - \phi_k^{m})]} \sqrt{(l+1)(m+1)} |l+1\rangle_j \langle l|_j \otimes |m\rangle_k \langle m+1|_k + \mathrm{h.c.} \right), \\
\tilde{H}_{\mathrm{d}} & = \sum_{jk} \alpha_{jk} \frac{\Omega_j}{2} \left( r_j(t) e^{-i (\nu_j t - \rho_{jk})} + \mathrm{c.c.} \right) \sum_l \left( e^{i (\xi_k^{l} t + \phi_k^{l})} \sqrt{l+1} |l + 1 \rangle_k \langle l |_k + \mathrm{h.c.} \right).\end{split}\]
Frame change operations
Frame change between \(U_f(t)\) and \(U_g(t)\) is effected by the change-of-frame unitary
\(V_{gf}(t) := U_g(t) U_f^{\dagger}(t)\):
\[V_{gf}(t) |\tilde{\psi}(t)\rangle_f = V_{gf}(t) U_f(t) |\psi(t)\rangle = U_g(t) |\psi(t)\rangle = |\tilde{\psi}(t)\rangle_g.\]
The Schrodinger equation in integral form in frame \(f\) is
\[\begin{split}|\tilde{\psi}(t_1)\rangle_f = U_f(t_1) K(t_1; t_0) |\psi(t_0)\rangle & = U_f(t_1) K(t_1; t_0) U_f^{\dagger}(t_0) |\tilde{\psi}(t_0)\rangle_f \\
& =: \tilde{K}_f(t_1; t_0) |\tilde{\psi}(t_0)\rangle_f\end{split}\]
where \(K(t_1; t_0) := T \left[ \exp \left(-i \int_{t_0}^{t_1}dt H(t) \right)\right]\) is the lab-frame time evolution operator. Time evolution operator in frame \(g\) is
\[\tilde{K}_g(t_1; t_0) = V_{gf}(t_1) \tilde{K}_f(t_1; t_0) V_{gf}^{\dagger}(t_0).\]
As noted above, \(\tilde{H}_f := U_f H U_f^{\dagger} + i \dot{U}_f U_f^{\dagger}\) (ommitting the time dependence) is the Hamiltonian in frame \(f\). The Hamiltonian in frame \(g\) is obtained from \(\tilde{H}_f\) by
\[\begin{split}\tilde{H}_g & = U_g U_f^{\dagger} \tilde{H}_f U_f U_g^{\dagger} - i U_g U_f^{\dagger} \dot{U}_f U_g^{\dagger} + i \dot{U}_g U_g^{\dagger} \\
& = V_{gf} \tilde{H}_f V_{gf}^{\dagger} + i \dot{V}_{gf} V_{gf}^{\dagger}.\end{split}\]
Finally, given the frame transformation of the states, the expectation value of an observable \(O\) in frame \(f\) changes its value in frame \(g\) to
\[\langle O \rangle_f = \langle \tilde{\psi}|_g V_{gf} O V_{gf}^{\dagger} |\tilde{\psi}\rangle_g = \langle V_{gf} O V_{gf}^{\dagger} \rangle_g.\]
Rotating-wave approximation
When \(|\nu_j + \xi_k^l| \gg |\nu_j - \xi_k^l|\) for all \(j, k, l\), we can apply the rotating-wave approximation (RWA) to the drive Hamiltonian and ignore
the fast-oscillating terms:
\[\bar{H}_{\mathrm{d}} = \sum_{jk} \alpha_{jk} \frac{\Omega_j}{2} \left( r_j(t) e^{i \rho_{jk}} \sum_l e^{-i (\epsilon_{jk}^l t - \phi_{k}^l)} \sqrt{l+1} |l+1\rangle_k \langle l |_k + \mathrm{h.c.} \right),\]
where \(\epsilon_{jk}^l := \nu_j - \xi_k^l\).
The RWA drive Hamiltonian in the frequency domain is (assuming \(\tilde{r}_j\) has support only around the frame frequencies)
\[\bar{H}_{\mathrm{d}} = \sum_{jk} \alpha_{jk} \frac{\Omega_j}{2} \int d\nu \left( \tilde{r}_j(\nu) e^{i \rho_{jk}} \sum_l e^{-i [(\nu - \xi_k^l) t - \phi_{k}^l]} |l+1\rangle_k \langle l |_k + \mathrm{h.c.} \right).\]
Drive Hamiltonian
The drive Hamiltonian for a given channel \(j\), qudit \(k\), level \(l\) is
\[\tilde{H}_{\mathrm{d}}|_{jk}^{l} = \alpha_{jk} \frac{\Omega_j}{2} \left( r_j(t) e^{-i (\nu_j t - \rho_{jk})} + \mathrm{c.c.} \right) \left( e^{i (\xi_k^{l} t + \phi_k^{l})} \sqrt{l+1} |l + 1 \rangle_k \langle l |_k + \mathrm{h.c.} \right).\]
Let
\[R_{jk}(t) = \alpha_{jk} e^{i\rho_{jk}} \frac{\Omega_j}{2} r_j(t)\]
and
\[\begin{split}A^{l}_{k} & = e^{-i \phi^{l}_{k}} \sqrt{l + 1} | l \rangle_k \langle l + 1 |_k, \\
X^{l}_{k} & = A^{l\dagger}_{k} + A^{l}_{k} \\
Y^{l}_{k} & = i(A^{l\dagger}_{k} - A^{l}_{k}).\end{split}\]
Then the drive term above is
\[\begin{split}\tilde{H}_{\mathrm{d}}|_{jk}^{l} = & R_{jk}(t) e^{-i (\nu_j - \xi_k^{l}) t} A^{l\dagger}_{k} + R^{*}_{jk}(t) e^{i (\nu_j - \xi_k^{l}) t} A^{l}_{k}
+ R^{*}_{jk}(t) e^{i (\nu_j + \xi_k^{l}) t} A^{l\dagger}_{k} + R_{jk}(t) e^{-i (\nu_j + \xi_k^{l}) t} A^{l}_{k} \\
= & \mathrm{Re}[R_{jk}(t) e^{-i (\nu_j - \xi_k^{l}) t} + R_{jk}(t) e^{-i (\nu_j + \xi_k^{l}) t}] X^{l}_{k}
+ \mathrm{Im}[R_{jk}(t) e^{-i (\nu_j - \xi_k^{l}) t} - R_{jk}(t) e^{-i (\nu_j + \xi_k^{l}) t}] Y^{l}_{k} \\
= & 2 \mathrm{Re}[R_{jk}(t) e^{-i \nu_j t}] [\cos(\xi_k^{l} t) X^{l}_{k} + \sin(\xi_k^{l} t) Y^{l}_{k}].\end{split}\]
With the rotating-wave approximation we instead have
\[\begin{split}\bar{H}_{\mathrm{d}}|_{jk}^{l} = & R_{jk}(t) e^{-i (\nu_j - \xi_k^{l}) t} A^{l\dagger}_{k} + R^{*}_{jk}(t) e^{i (\nu_j - \xi_k^{l}) t} A^{l}_{k} \\
= & \mathrm{Re}[R_{jk}(t) e^{-i (\nu_j - \xi_k^{l}) t}] X^{l}_{k} + \mathrm{Im}[R_{jk}(t) e^{-i (\nu_j - \xi_k^{l}) t}] Y^{l}_{k}.\end{split}\]
The representation in terms of \(X^{l}_{k}\) and \(Y^{l}_{k}\) operators has several advantages over using \(A^{l}_{k}\) and \(A^{l\dagger}_{k}\):
When \(r_j(t)\) is a callable, QuTiP sesolve seems to run slightly faster when \(X\) and \(Y\) with real coefficients are passed as Hamiltonian terms.
The Hamiltonian is manifestly Hermitian.
For a pure real or imaginary \(R_{jk}(t)\), on-resonant (\(\nu_j = \xi_k^{l}\)) drive, the RWA Hamiltonian reduces to a single term.