Using the MCWF method

Introduction

An easy example for the use of the MCWF method [13, 14, 15, 16] in QDYN is to consider a two-qubit system with spontaneous decay and dephasing on both qubits.\(\newcommand{\bra}[1]{\langle #1 \vert}\newcommand{\ket}[1]{\vert #1 \rangle}\),

The qdyn_prop_traj utility that comes as part of the QDYN installation can propagate the system both as trajectories in the MCWF method, as well as the full density matrix, giving the “exact” solution for comparison. We assume the QDYN has been compiled with MPI support, so that the trajectories can be calculated in parallel.

! qdyn_prop_traj -v

QDYN 2.0dev revision e5a08b2e2c6c21846bec4b7a4f7cafde09ba9497 (master)
  features: no-check-cheby, no-check-newton, no-parallel-ham, use-mkl=sequential, use-mpi=intel, parallel-oct, backtraces, no-debug, no-no-ipo
  compiled with ifort on Tue Sep 19 09:17:49 2017 on host mlhpc2
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=1
:
system msg for write_line failure : Bad file descriptor

Two-Qubit Hilbert space

For basic linear algebra support, as well as input/output, we use the standard numpy library,

import numpy as np

We can now define the basic building blocks of the Hilbert space, including the vector representation of \(\ket{0}\), \(\ket{1}\), a basis for arbitrary operators through ket-bras, as well as the identity.

def ket(i):
    k = np.matrix(np.zeros((2,1), dtype=np.complex128))
    k[int(i)] = 1.0
    return k

def ketbra(i, j):
    return ket(i) * ket(j).H

Id = np.matrix(np.eye(2))

Parameters

The system system energies are chosen somewhat arbitrarily in the GHZ regime, with decay term on the order of magnitude of MHz. This is sufficiently strong to see strong decay on the order of a few hundred nanosecond.

w1 = 1.0 # GHz
w2 = 1.2 # GHz
gamma_decay_1 = 2.3 # MHz
gamma_decay_2 = 2.1 # MHz
gamma_dephase_1 = 1.1 # MHz
gamma_dephase_2 = 1.2 # MHz

Note that a factor of \(2 \pi\) is implicit (QDYN knows how to work with GHz/MHz as an energy unit)

Hamiltonian, Lindblad operators, and Observables

We new construct the Hamiltonian, as well as the four Lindblad operators for decay and dephasing on each qubit.

H = w1 * np.kron(ketbra(1,1), Id) + w2 * np.kron(Id, ketbra(1,1))

Ls = [
    np.sqrt(gamma_decay_1) * np.kron(ketbra(0,1), Id),
    np.sqrt(gamma_decay_2) * np.kron(Id, ketbra(0,1)),
    np.sqrt(gamma_dephase_1) * np.kron(ketbra(1,1), Id),
    np.sqrt(gamma_dephase_2) * np.kron(Id, ketbra(1,1))
]

As observables, we would like to show the population on each of the four basis levels

Os = []
for i in range(4):
    O = np.matrix(np.zeros((4,4), dtype=np.complex128))
    O[i,i] = 1.0
    Os.append(O)

The initial state is \(\ket{11}\), which will decay to \(\ket{01}\) and \(\ket{10}\), and eventually to \(\ket{11}\)

psi0 = np.kron(ket(1), ket(1))

Setting up the MCWF Model

To set up the model (i.e., generate the config file and all the depdendent data), we will use the QDYN Python library.

import qdyn

mcwf_model = qdyn.model.LevelModel()
mcwf_model.add_ham(H, specrad_method='diag', op_unit='GHz')
for L in Ls:
    mcwf_model.add_lindblad_op(L, op_unit='sqrt_MHz')
nt = 1000
mcwf_model.set_propagation(
    T=500, nt=nt, time_unit='ns', use_mcwf=True, prop_method='newton', initial_state=psi0)
for i, O in enumerate(Os):
    mcwf_model.add_observable(O, 'pop.dat', 'unitless', 'ns', '<P_%d>'%(i+1))
mcwf_model.write_to_runfolder('./rf_2q_decay')

The resulting config file will have all the necessary settings for the MCWF propagation

! cat ./rf_2q_decay/config

tgrid: t_start = 0, t_stop = 500_ns, nt = 1000

prop: method = newton, use_mcwf = T, mcwf_order = 2

ham:
* type = matrix, filename = H0.dat, real_op = F, n_surf = 4, &
  specrad_method = diag, op_unit = GHz, sparsity_model = banded, op_type = pot

dissipator: type = lindblad_ops, real_op = F, op_unit = sqrt_MHz, &
  add_to_H_jump = banded, conv_to_superop = F, sparsity_model = banded
* filename = L1.dat
* filename = L2.dat
* filename = L3.dat
* filename = L4.dat

observables: type = matrix, real_op = F, n_surf = 4, outfile = pop.dat, &
  exp_unit = unitless, is_real = T, time_unit = ns, op_unit = unitless, &
  sparsity_model = banded
* filename = O1.dat, column_label = <P_1>
* filename = O2.dat, column_label = <P_2>
* filename = O3.dat, column_label = <P_3>
* filename = O4.dat, column_label = <P_4>

psi:
* type = file, filename = psi.dat

When running an MCWF propagation of this runfolder, the averaged population dynamics will be written to the file pop.dat inside the runfolder. By using qdyn_prop_traj’s --write-statistical-squares option, we will also get a file pop.dat.statsq that contains the square roots of the average of the squares of the population. This data can be used to estimate the statistical variance of the results, from which the statistical error can be derived.

Density Matrix Propagation

We want to see how well the MCWF propagation does in terms of the absolute error. Thus, we first propagate the corresponding density matrix, which gives us the “exact” solution of the system

rho_model = qdyn.model.LevelModel()
rho_model.add_ham(H, specrad_method='diag', op_unit='GHz')
for L in Ls:
    rho_model.add_lindblad_op(L, op_unit='sqrt_MHz')
rho_model.set_propagation(T=500, nt=nt, time_unit='ns', prop_method='newton', initial_state=psi0)
for i, O in enumerate(Os):
    rho_model.add_observable(O, 'pop.dat', 'unitless', 'ns', '<P_%d>'%(i+1))
rho_model.write_to_runfolder('./rf_2q_decay_rho')

! qdyn_prop_traj --rho ./rf_2q_decay_rho

***** START OF PROGRAM qdyn_prop_traj ******
Tue Sep 19 10:19:55 -0700 2017

*** Read config file ./rf_2q_decay_rho/config ***
*** Done reading config file ***
*** Initializing system ***
  *** Initializing grid ***
  No explicit grid found in config file! Initializing grid as one dimensional 1x1 fake grid
    Initializing grid as 1D cartesian grid without mapping
      Number of grid points was set to nr = 1 for dimension dim = 1! Initializing dimension 1 as a 1x1 fake-grid...
  *** Initializing pulses ***
  *** Initializing dynamical generator ***
*** Done with initialization ***

Completed traj 1/1 in proc 1/1
Writing avg in proc 1/1

Tue Sep 19 10:19:55 -0700 2017
***** END OF PROGRAM qdyn_prop_traj ******

tgrid, pop00_rho, pop01_rho, pop10_rho, pop11_rho = np.genfromtxt('./rf_2q_decay_rho/pop.dat', unpack=True)

We keep these population dynamics for future reference:

EXACT = (pop00_rho, pop01_rho, pop10_rho, pop11_rho)

Visualization

In the following, we prepare some visualizations for the population dynamcis and the errors, using the matplotlib library.

%matplotlib inline

import matplotlib.pylab as plt

The statistical error for an observable \(P\) (the population, in our case) when averaging over \(N\) MCWF trajectories is given by \(\sqrt{\text{var}(P) / N}\), where the variance is \(\text{var}(P) = E(P^2) - E(P)^2\). See Mølmer et al, J. Opt. Soc. Am. B 10, 524 (1993)] for a more detailed discussion. The mean

\[E(P^2) = \frac{1}{N} \sum_k \langle P \rangle_k^2\]

can be obtained from pop.dat.statsq. Note that pop.dat.statsq stores \(\sqrt{E(P^2)}\), not \(E(P^2)\), just so that the columns have the same units as the expectation values \(E(P)\) stored in pop.dat.

from os.path import join
def variances(rf):
    tgrid, pop00_mcwf, pop01_mcwf, pop10_mcwf, pop11_mcwf = np.genfromtxt(
        join(rf, 'pop.dat'), unpack=True)
    tgrid, pop00_s, pop01_s, pop10_s, pop11_s = np.genfromtxt(
        join(rf,'pop.dat.statsq'), unpack=True)
    var_pop00 = pop00_s**2 - pop00_mcwf**2
    var_pop01 = pop01_s**2 - pop01_mcwf**2
    var_pop10 = pop10_s**2 - pop10_mcwf**2
    var_pop11 = pop11_s**2 - pop11_mcwf**2
    for var in (var_pop00, var_pop01, var_pop10, var_pop11):
        # avoid negative variance due to numerical error
        assert np.min(var) > -1e-14
        var[var < 0] = 0
    return var_pop00, var_pop01, var_pop10, var_pop11

Again, the statistical error is now simply the normalized \(\sqrt{\text{var}(P) / N}\):

def statistical_errors(rf, n_trajs):
    var_pop00, var_pop01, var_pop10, var_pop11 = variances(rf)
    delta_pop00 = np.sqrt(var_pop00) / np.sqrt(n_trajs)
    delta_pop01 = np.sqrt(var_pop01) / np.sqrt(n_trajs)
    delta_pop10 = np.sqrt(var_pop10) / np.sqrt(n_trajs)
    delta_pop11 = np.sqrt(var_pop11) / np.sqrt(n_trajs)
    return delta_pop00, delta_pop01, delta_pop10, delta_pop11

The plot of the population dynamics will include the averaged population, the exact population, and the statistical error \(\delta P\) (as “error bars” \(\pm 3 \delta P\))

def show_popdyn(rf, n_trajs, exact=EXACT, ax=None, figsize=None):
    tgrid, pop00_mcwf, pop01_mcwf, pop10_mcwf, pop11_mcwf = np.genfromtxt(
        join(rf, 'pop.dat'), unpack=True)
    delta_pop00, delta_pop01, delta_pop10, delta_pop11 = statistical_errors(rf, n_trajs)
    pop00_rho, pop01_rho, pop10_rho, pop11_rho = exact
    show_fig = False
    if ax is None:
        fig = plt.figure(figsize=figsize)
        ax = fig.add_subplot(111)
        show_fig = True
    civ = 3  # confidence interval to plot (sigmas)
    ax.fill_between(
        tgrid, pop11_mcwf-civ*delta_pop11, pop11_mcwf+civ*delta_pop11,
        color='red', alpha=0.2)
    ax.fill_between(
        tgrid, pop10_mcwf-civ*delta_pop10, pop10_mcwf+civ*delta_pop10,
        color='blue', alpha=0.2)
    ax.fill_between(
        tgrid, pop01_mcwf-civ*delta_pop01, pop01_mcwf+civ*delta_pop01,
        color='purple', alpha=0.2)
    ax.fill_between(
        tgrid, pop00_mcwf-civ*delta_pop00, pop00_mcwf+civ*delta_pop00,
        color='magenta', alpha=0.2)
    ax.plot(tgrid, pop11_mcwf, label='11', color='red')
    ax.plot(tgrid, pop10_mcwf, label='10', color='blue')
    ax.plot(tgrid, pop01_mcwf, label='01', color='purple')
    ax.plot(tgrid, pop00_mcwf, label='00', color='magenta')
    ax.plot(tgrid, pop11_rho, ls='dotted', label='', color='red')
    ax.plot(tgrid, pop10_rho, ls='dotted', label='', color='blue')
    ax.plot(tgrid, pop01_rho, ls='dotted', label='', color='purple')
    ax.plot(tgrid, pop00_rho, ls='dotted', label='', color='magenta')
    ax.set_ylabel("population")
    ax.set_xlabel("time (ns)")
    ax.legend()
    if show_fig:
        plt.show(fig)

We also prepare a plot of the absolute error (absolute value of average population minus exact population)

def show_error(rf, exact=EXACT, ax=None, figsize=None):
    tgrid, pop00_mcwf, pop01_mcwf, pop10_mcwf, pop11_mcwf = np.genfromtxt(join(rf, 'pop.dat'), unpack=True)
    pop00_rho, pop01_rho, pop10_rho, pop11_rho = exact
    show_fig = False
    if ax is None:
        fig = plt.figure(figsize=figsize)
        ax = fig.add_subplot(111)
        show_fig = True
    ax.plot(tgrid, np.abs(pop11_mcwf - pop11_rho), label='err_11', color='red')
    ax.plot(tgrid, np.abs(pop10_mcwf - pop10_rho), label='err_10', color='blue')
    ax.plot(tgrid, np.abs(pop01_mcwf - pop01_rho), label='err_01', color='purple')
    ax.plot(tgrid, np.abs(pop00_mcwf - pop00_rho), label='err_00', color='magenta')
    ax.set_ylabel("population error")
    ax.set_xlabel("time (ns)")
    ax.legend()
    if show_fig:
        plt.show(fig)

Both of these plots will be shown side-by-side in a combined plot:

def show_combined(rf, n_trajs, exact=EXACT, figsize=None):
    fig, axs = plt.subplots(ncols=2, figsize=figsize)
    show_popdyn(rf, n_trajs, exact, ax=axs[0])
    show_error(rf, exact, ax=axs[1])

All plots will have the same (arbitrary) figure size (inches, at 72 dpi)

FIGSIZE = (16,4)

MCWF Propagation with varying number of trajectories

We will run multiple trajectories in parallel using MPI.

N_PROCS = 78 # the number of CPUs

The above value must be less than than the number of CPUs available on the machine — otherwise, MPI will get extremely slow.

We can now propagate for varying numbers of trajectories, with the results shown below

1e1 trajectories

mcwf_model.write_to_runfolder('./rf_2q_decay_1e1')

! mpirun -n 1 qdyn_prop_traj --write-statistical-squares --seed=0 --n-trajs=10 ./rf_2q_decay_1e1 > out.log

show_combined('./rf_2q_decay_1e1', 1e1, figsize=FIGSIZE)

1e2 trajectories

mcwf_model.write_to_runfolder('./rf_2q_decay_1e2')

! mpirun -n {N_PROCS} qdyn_prop_traj --write-statistical-squares --seed=0 --n-trajs=100 ./rf_2q_decay_1e2 > out.log

show_combined('./rf_2q_decay_1e2', 1e2, figsize=FIGSIZE)

1e3 trajectories

mcwf_model.write_to_runfolder('./rf_2q_decay')

! mpirun -n {N_PROCS} qdyn_prop_traj --write-statistical-squares --seed=0 --n-trajs=1000 ./rf_2q_decay > out.log

show_combined('./rf_2q_decay', 1e3, figsize=(16,4))

1e4 trajectories

mcwf_model.write_to_runfolder('./rf_2q_decay_1e4')

! mpirun -n {N_PROCS} qdyn_prop_traj --write-statistical-squares --seed=0 --n-trajs=10000 ./rf_2q_decay_1e4 > out.log

show_combined('./rf_2q_decay_1e4', 1e4, figsize=FIGSIZE)

1e5 trajectories

mcwf_model.write_to_runfolder('./rf_2q_decay_1e5')

! mpirun -n {N_PROCS} qdyn_prop_traj --write-statistical-squares --seed=0 --n-trajs=100000 ./rf_2q_decay_1e5 > out.log

show_combined('./rf_2q_decay_1e5', 1e5, figsize=FIGSIZE)

1e6 trajectories

mcwf_model.write_to_runfolder('./rf_2q_decay_1e6')

! mpirun -n {N_PROCS} qdyn_prop_traj --write-statistical-squares --seed=0 --n-trajs=1000000 ./rf_2q_decay_1e6 > out.log

show_combined('./rf_2q_decay_1e6', 1e6, figsize=FIGSIZE)

Convergence Analysis

As discussed above, for an observable \(P\) averaged over \(N\) MCWF trajectories the statistical errr is \(\delta P = \sqrt{\text{var}(P) / N}\). To get a single number that we can compare across the varying number of trajectories, we average over the time grid, as well as over the four expectation values \(P_{00}\), \(P_{01}\), \(P_{10}\) and \(P_{11}\).

def avg_variance(rf):
    var_pop00, var_pop01, var_pop10, var_pop11 = variances(rf)
    var_total = (var_pop00 + var_pop01 + var_pop10 + var_pop11) / 4.0
    return np.mean(var_total)

Likewise, we average the actual absolute error, \(\Delta E\)

def avg_error(rf, exact=EXACT):
    tgrid, pop00_mcwf, pop01_mcwf, pop10_mcwf, pop11_mcwf = np.genfromtxt(join(rf, 'pop.dat'), unpack=True)
    pop00_rho, pop01_rho, pop10_rho, pop11_rho = exact
    err_total = abs(pop00_mcwf - pop00_rho)
    err_total += abs(pop01_mcwf - pop01_rho)
    err_total += abs(pop10_mcwf - pop10_rho)
    err_total += abs(pop11_mcwf - pop11_rho)
    err_total *= 0.25
    return np.mean(err_total)

We have derived the variances using the values in pop.dat.statsq, which was generated by giving the --write-statistical-squares flag. However, we can get an upper bound of the statistical variance from the quantum variance \(\langle P^2 \rangle - \langle P \rangle^2\), with

\[\langle P^2 \rangle = \frac{1}{N} \sum_k \langle P^2 \rangle_k \quad\geq\quad E(P^2) = \frac{1}{N} \sum_k \langle P \rangle^2_k\]

Thus, we find:

\[\text{var}(P) \leq \langle P^2 \rangle - \langle P \rangle^2\]

.

For our specific case, \(P\) is a projector, i.e. \(P^2 = P\). Moreover, in every trajectory, the system is always in an eigenstate of \(P\), in which case the inequality becomes exact [cf. Mølmer et al, J. Opt. Soc. Am. B 10, 524 (1993), Section 7]

We thus define the following alternative for the variance:

def avg_variance2(rf):
    tgrid, pop00_mcwf, pop01_mcwf, pop10_mcwf, pop11_mcwf = np.genfromtxt(
        join(rf, 'pop.dat'), unpack=True)
    var_total = pop00_mcwf - pop00_mcwf**2
    var_total += pop01_mcwf - pop01_mcwf**2
    var_total += pop10_mcwf - pop10_mcwf**2
    var_total += pop11_mcwf - pop11_mcwf**2
    var_total *= 0.25
    return np.mean(var_total)

We now collect the variance, the estimated statistical error, and the actual absolute error in a table:

import pandas as pd
from collections import OrderedDict
rfs = OrderedDict([
    (1e1, './rf_2q_decay_1e1'), (1e2, './rf_2q_decay_1e2'),
    (1e3, './rf_2q_decay'), (1e4, './rf_2q_decay_1e4'),
    (1e5, './rf_2q_decay_1e5'), (1e6, './rf_2q_decay_1e6')])
convergence_table = pd.DataFrame(
    data=OrderedDict([
        ('N', list(rfs.keys())),
        ('variance', [avg_variance(rfs[N]) for N in rfs.keys()]),
        ('<P²> - <P>²', [avg_variance2(rfs[N]) for N in rfs.keys()]),
        ('δP', np.sqrt([avg_variance2(rfs[N])/float(N) for N in rfs.keys()])),
        ('ΔP', [avg_error(rfs[N]) for N in rfs.keys()]),
    ])
)
convergence_table['ΔP/δP'] = convergence_table['ΔP']  / convergence_table['δP']

convergence_table

	N	variance	<P²> - <P>²	δP	ΔP	ΔP/δP
0	10.0	0.050750	0.050750	0.071239	0.048418	0.679649
1	100.0	0.057399	0.057399	0.023958	0.012505	0.521956
2	1000.0	0.059202	0.059202	0.007694	0.004177	0.542915
3	10000.0	0.059975	0.059975	0.002449	0.001308	0.534225
4	100000.0	0.060375	0.060375	0.000777	0.000458	0.589062
5	1000000.0	0.060316	0.060316	0.000246	0.000176	0.717991

The variance converges quite rapidly towards the the exact value:

variance_exact = avg_variance2('./rf_2q_decay_rho');
variance_exact

0.060243964718507251

To a good approximation, we can assume the variance as independent of \(N\), in which case the statistical error scales exacly as \(1/\sqrt{N}\). For an order of magnitude improvement in the statistical error, a two order of magnitude increase in the number of trajectories is required. Note that we can also see this in the y-axis of the plots of the absolute error over time in the previous section.

If you cannot reproduce an absolute error that is compatible with the \(\delta P\) in the above table, the culprit is most likely an insufficient random number generator.