Using the variance functional to optimize squeezing

QDYN also offers the possibility of selecting a variance of an operator as a target functional. We show how the system can be initialized using the QDYN-pylib and then how the target functional can be set inside the Fortran90 program.

As a model system we choose an anharmonic oscillator with Hamiltonian

\begin{equation} \hat{H} = \omega \hat{a}^\dagger \hat{a} + \kappa \hat{a}^\dagger \hat{a} \hat{a}^\dagger \hat{a} + \hat{H}_\text{dr}(t), \end{equation}

where \(\omega\) is the frequency of the harmonic oscillator and \(\kappa\) a small perturbation parameter. Further, we control the system via a time-dependent two-photon-drive

\begin{equation} \hat{H}_\text{dr}(t) = -i \left[u(t)e^{-2i\omega t}\left(\hat{a}^\dagger\right)^2 - u^*(t)e^{2i\omega t} \hat{a}^2\right]. \end{equation}

In the rotating frame, the Hamiltonian is

\begin{equation} \hat{H} = \kappa \hat{a}^\dagger \hat{a} \hat{a}^\dagger \hat{a} -i \left[u(t) \left(\hat{a}^\dagger\right)^2 - u^*(t) \hat{a}^2\right]. \end{equation}

The goal is now to drive the system into a squeezed state, i.e. a state for which the variance of one of the quadratures is decreased below the zero-point fluctuation achieved in the ground state.

For a mechanical oscillator the quadratures correspond to position and momentum, but they can also correspond to real and imaginary part of the complex amplitude given by \(\hat{a}\) for an electromagnetic field (cf. Wikipedia).

Thus, to keep it general we define the quadratures as

\begin{equation} \hat{X}_1 = \frac{1}{\sqrt{2}}\left(\hat{a}^\dagger + \hat{a}\right); \quad \hat{X}_2 = \frac{i}{\sqrt{2}} \left(\hat{a}^\dagger - \hat{a} \right) \end{equation}

The Heisenberg uncertainty principle then states:

\begin{equation} \Delta X_1^2 \Delta X_2^2 \geq \frac{1}{4} \end{equation}

In the following, we optimize towards \(\hat{X}_1\)-squeezed states (in the rotating frame). Therefore, we use real pulse amplitudes. A state is then considered squeezed if \(\Delta X_1^2 < \frac{1}{2}\).

We now set up the model in Python:

from pathlib import Path

import numpy as np
import qutip

import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable
plt.rcParams.update({'font.size': 15})


import qdyn.model
import qdyn.pulse
from qdyn.io import read_psi_amplitudes

!make clean #delete all compiled programs and all data from previous runs

Setting up the model

First, define a function which sets the Hamiltonian.

def get_Ham(
    omega=1,
    kappa=0.1,
    N_lvls=100,  #numerical number of levels
    u=None,      #functional form of guess pulse; if None: use constant
    E0=0.1       #constant amplitude of guess pulse
):

    Ham_1 = kappa * qutip.num(N_lvls)**2

    #if nothing given, implement constant guess pulse
    #implement real guess pulse since X1, the real part of the light field, is squeezed
    if not u:
        u = lambda t: E0

    #for a real pulse u = u*
    Ham_dr = -1j*(qutip.create(N_lvls)**2 - qutip.destroy(N_lvls)**2)

    return [Ham_1, [Ham_dr, np.vectorize(u)]]

Then, define a function which creates the runfolder containing a ‘config’ file with all the necessary data and files that contain the parts of the Hamiltonian, the observables of interest, the initial state and the time-dependent guess pulse. The runfolder is named “rf” as default.

All the optimization parameters are also set within this function. Feel free to play around with them!

def qdyn_model(
    Ham,                      #Hamiltonian
    qdyn_tlist,               #tuple of the form (T, Nt)
    psi0,                     #initial state
    N_lvls,                   #numerical system size
    runfolder="rf",
    prop_method="cheby",
    iter_stop=100,
    **user_kwargs,
):
    """ Create qdyn-model

    """
    #Create time grid
    dt = (qdyn_tlist[0]) / (qdyn_tlist[1] - 1)
    tgrid = np.linspace(
        float(dt / 2),
        float(qdyn_tlist[0] - dt / 2),
        qdyn_tlist[1] - 1,
        dtype=np.float64,
    )

    #Define quadrature operators
    X1 = (qutip.create(N_lvls) + qutip.destroy(N_lvls))/np.sqrt(2)
    X2 = 1j*(qutip.create(N_lvls) - qutip.destroy(N_lvls))/np.sqrt(2)

    #Initialize model
    model = qdyn.model.LevelModel()
    for H in Ham:
        if isinstance(H, list):
            model.add_ham(
                H[0],
                pulse=qdyn.pulse.Pulse(
                    tgrid,
                    amplitude=H[1](tgrid),
                    time_unit="iu",
                    ampl_unit="iu",
                    config_attribs={
                        'filename': 'pulse_initial.dat',
                        #Shape function by which the change in amplitude is multiplied
                        #after each iteration, const means no change
                        'oct_shape': 'const',
                        #Maximum amplitude increase factor per iteration
                        'oct_increase_factor': 50,
                        #Minimum allowed pulse amplitude, if limit_pulses is True
                        'oct_pulse_min': 0,
                        #Maximum allowed pulse amplitude, if limit_pulses is True
                        'oct_pulse_max': 1.0e10,
                        #scale for the step size in each iteration
                        #(the higher, the smaller the step size)
                        'oct_lambda_a': 5.0e0,
                        #file into which the optimized pulse will be written
                        'oct_outfile': 'pulse_oct.dat'
                    },
                ),
                op_unit="iu",
            )
        else:
            model.add_ham(H, op_unit="iu")

    model.set_propagation(
        qdyn_tlist[0], qdyn_tlist[1], time_unit="iu", prop_method=prop_method
    )
    model.set_oct(
        method="krotovpk",
        max_ram_mb=230000,
        #OCT is stopped when target functional reaches this value
        J_T_conv=1e-2,
        #Stop opt. when relative change in target functional is smaller than this value
        delta_J_T_conv=1e-5,
        iter_dat="oct_iters.dat",
        #If True, continue opt. from a previous run
        continue_=False,
        params_file="oct_params.dat",
        limit_pulses=False,
        keep_pulses="prev",
        #stop opt. independent of outcome after this number of steps
        iter_stop=iter_stop
    )

    #Add observables whose exp. value should be calculated during propagation,
    #here both quadratures and their squares
    model.add_observable(X1, "squeeze.dat", exp_unit="iu",
                         time_unit="iu", col_label="<X1>", square="<X1sq>")
    model.add_observable(X2, "squeeze.dat", exp_unit="iu",
                         time_unit="iu", col_label="<X2>", square="<X2sq>")

    model.add_state(psi0, 'initial')
    user_variables = {}
    user_variables["runfolder"] = str(Path(runfolder))
    user_variables.update(user_kwargs)

    model.user_data = user_variables
    model.write_to_runfolder(str(runfolder))  #write everything to runfolder


    return model

Define the relevant system and optimization parameters, time grid and initial state:

N_lvls = 150   #numerical system size
kappa  = 1e-3

T      = 5.0   #final time
Nt     = 501   #number of time steps
qdyn_tlist = (T, Nt)

N_iter = 20    #number of iterations the optimization should run

psi0 = qutip.basis(N_lvls, 0) #use ground state as initial state

Set up model and runfolder:

H = get_Ham(kappa=kappa, N_lvls=N_lvls, E0=0.01)

m = qdyn_model(
    H,
    qdyn_tlist,
    psi0,
    N_lvls=N_lvls,
    iter_stop=N_iter
)

Implementing the optimization in Fortran

Now that the runfolder is set up, we conduct the optimization inside Fortran. The Fortran file is structured as follows:

1. In the upper part all the necessary variables are defined, and functions imported.

2. The model is then initialized. This is done by reading all the parameters we set in Python using the read_para routine. The name of the runfolder is taken directly from the command line.

   Then the system is initialized with the init routine, in which also the Hamiltonian, pulses and observables are read from the files. The initial state is set inside the init_psi routine.

call get_command_argument(1, rf)
call read_para(para, rf, 'config')

call init(para, grid=grid, gen=gen, pulses=pulses)

call init_psi(psi0, grid, gen%ham%nsurf, 1, para, psi0_state_label)

3. Now that the system is set up, we can start with the optimization. First we need to define our target functional. Therefore we read the \(\hat{X}_1\) - quadrature from the file in which it is defined as an observable. Then we pass it to the set_oct_target routine, since we need it during the optimization.

   Then we can already start with the optimization by calling the optimize_pulses routine. For that we have to pass a routine that calculates the target functional, here J_matrix_variance, and a routine that calculates the co-states, here chis_matrix_variance.

   Note that in this case, we use a single Hermitian operator. If you want to use multiple or non-Hermitian operators, see the documentation of J_matrix_variance.

call read_op_matrix(X1, join_path(rf, 'O1.dat'), gen%ham%nsurf,              &
                    op_type="potential", pulse=0, n_photons=0,               &
                    conjg_pulse=.false., expand_hermitian=.false.,           &
                    sparsity_model = "indexed")

call set_oct_target(targets, 1, 1, psi0, gen, grid, matrix_op=(/X1/))

call optimize_pulses(targets, fw_states_T, pulses, para, J_matrix_variance,  &
&                    chis_matrix_variance, g_b_zero)

4. After the optimization stopped, we finally propagate the system to be able to analyze the system dynamics under the influence of the optimized pulse. Therefore, we initialize the propagation using the init_prop routine and propagate by calling prop.

   During the propagation, the expectation values of the observables we set in the beginning for each time step is calculated and written inside the corresponding file. If additional data is needed, additional observables can be set inside the qdyn_model function above.

   Alternatively, one can define a prop_info_hook routine inside Fortran and pass it to prop. The routine is then called after each time step.

   As a last step we write out the final state.

call init_prop(gen, grid, prop_work, para, pulses)
call prop(psi0, grid, gen, prop_work, para, pulses)

call write_state(psi0, join_path(rf, 'psi_final_after_oct.dat'))

Run the optimization

By putting an exclamation mark in front of a command we can pass it directly to the shell. make now compiles the Fortran code and creates a program from it, which we then execute with ./squeeze. The name of the runfolder, “rf”, we pass as an argument after the command.

!make

make: 'squeeze' is up to date.

!./squeeze rf

QDYN 23.04+dev revision bc6d8d9659b96baa7c1b7b5ac1f461283e45fbec (110-add-example-for-minimizing-expectation-value-variance-of-an-)
  features: no-check-cheby, no-check-newton, no-parallel-ham, use-mkl=sequential, parallel-oct, backtraces, no-debug, no-no-ipo, no-shared-lib
  compiled with ifort on Wed Jul 19 17:01:17 2023 on host talbiya
*** Read config file rf/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 ***
    1 pulses in config file
    Initializing pulse 1
      setting pulse shape as constant   1.00 in [   0.000E+00,   5.000E+00]
  *** Initializing dynamical generator ***
*** Done with initialization ***

Initializing Krotov
Storage used for target     1:
  bw prop (old):         4 MB (RAM);          0 MB (disk)
Total: 4 MB (RAM); 0 MB (disk)
Number of processes: 1
OMP number of threads per process: 16
Process 1, Thread 1/8: Initial forward propagation of target 1
   Completed forward propagation of system   1 in        0.0 seconds
Done with Krotov Initialization
Iter. |         J_T |     g_a_int |     g_b_int |           J |      Delta J |
    0 | 4.09372E-01 | 0.00000E+00 | 0.00000E+00 | 4.09372E-01 |  4.09372E-01 |
    1 | 4.92125E-02 | 3.47937E-01 | 0.00000E+00 | 3.97149E-01 | -1.22230E-02 |
    2 | 4.17514E-02 | 5.48660E-03 | 0.00000E+00 | 4.72380E-02 | -1.97450E-03 |
    3 | 3.95118E-02 | 1.55089E-03 | 0.00000E+00 | 4.10627E-02 | -6.88726E-04 |
    4 | 3.81555E-02 | 8.24652E-04 | 0.00000E+00 | 3.89801E-02 | -5.31671E-04 |
    5 | 3.70142E-02 | 6.57760E-04 | 0.00000E+00 | 3.76719E-02 | -4.83569E-04 |
    6 | 3.59928E-02 | 5.77573E-04 | 0.00000E+00 | 3.65703E-02 | -4.43822E-04 |
    7 | 3.50651E-02 | 5.19378E-04 | 0.00000E+00 | 3.55844E-02 | -4.08319E-04 |
    8 | 3.42165E-02 | 4.71784E-04 | 0.00000E+00 | 3.46882E-02 | -3.76823E-04 |
    9 | 3.34363E-02 | 4.31357E-04 | 0.00000E+00 | 3.38676E-02 | -3.48825E-04 |
   10 | 3.27161E-02 | 3.96411E-04 | 0.00000E+00 | 3.31125E-02 | -3.23791E-04 |
   11 | 3.20490E-02 | 3.65845E-04 | 0.00000E+00 | 3.24148E-02 | -3.01250E-04 |
   12 | 3.14293E-02 | 3.38845E-04 | 0.00000E+00 | 3.17682E-02 | -2.80817E-04 |
   13 | 3.08524E-02 | 3.14776E-04 | 0.00000E+00 | 3.11671E-02 | -2.62176E-04 |
   14 | 3.03142E-02 | 2.93140E-04 | 0.00000E+00 | 3.06073E-02 | -2.45072E-04 |
   15 | 2.98113E-02 | 2.73535E-04 | 0.00000E+00 | 3.00849E-02 | -2.29299E-04 |
   16 | 2.93410E-02 | 2.55644E-04 | 0.00000E+00 | 2.95966E-02 | -2.14694E-04 |
   17 | 2.89006E-02 | 2.39214E-04 | 0.00000E+00 | 2.91399E-02 | -2.01129E-04 |
   18 | 2.84881E-02 | 2.24047E-04 | 0.00000E+00 | 2.87121E-02 | -1.88501E-04 |
   19 | 2.81014E-02 | 2.09988E-04 | 0.00000E+00 | 2.83114E-02 | -1.76731E-04 |
   20 | 2.77387E-02 | 1.96915E-04 | 0.00000E+00 | 2.79356E-02 | -1.65753E-04 |
*** Propagate with optimized pulse... ***
Finished on Fri Jul 21 09:54:20 +0200 2023

Analyze results

Now that the optimization is done, we can analyze the result. We first compare the initial and optimized pulse. Therefore we load them from their files, whose names we set to ‘pulse_initial.dat’ and ‘pulse_oct.dat’ in the qdyn_model function. As one can see, the amplitude was increased during the optimization.

We also set the quadratures and their squares as observables. During the final propagation, their expectation values were then calculated for each time step and stored in the corresponding file whose name we specified as ‘squeeze.dat’. From these we calculate the variances and plot them against time. One can see that the variance of \(\hat{X}_1\) decreases almost exponentially while the variance of \(\hat{X}_2\) is increased by the same amount.

These results are consistent with the theory since without the anharmonicity, the evolving unitary is

\begin{equation} \hat{U}(t) = \exp\left[\int_0^t u(t') \hat{a}^2 \text{d}t' - \text{H.c.}\right], \end{equation}

which is equivalent to the squeezing operator with squeezing parameter \(r(t)=\int_0^t u(t') \text{d}t'\) where \(\Delta \hat{X}_1^2(r) = \Delta \hat{X}_1^2(0) e^{-r}.\) For more information, see the Wikipedia page of Squeezed states.

Without the anharmonicity (and a sufficient amount of numerical levels to cover the phase space), one could in theory produce an infinitely squeezed state. However, the anharmonicity poses a bound to the amount of squeezing that can be acquired. This is the reason why it would not be sensible to further increase the amplitude and probably also why for the optimized pulse the amplitude is increased towards the end, since the perturbation is quadratic and thus gets more relevant for later times in which the oscillator’s higher levels get populated.

Note that depending on your guess pulse, you might achieve a different outcome. Feel free to play around with the anharmonicity and the guess pulse and see how the number of numerical levels limits the amount of squeezing that can be achieved!

fig, ax = plt.subplots(1, 2, figsize=(12,5))

### pulse plot ###
t_pulse, amp_oct = np.split(np.loadtxt("rf/pulse_oct.dat", usecols=(0,1)), 2, axis=1)
t_pulse, amp_init = np.split(np.loadtxt("rf/pulse_initial.dat", usecols=(0,1)), 2, axis=1)

ax[0].plot(t_pulse, amp_oct,  label='optimized pulse', linewidth=2.5)
ax[0].plot(t_pulse, amp_init, label='guess pulse', linestyle='--', linewidth=2.5)
ax[0].set_xlabel("time (i.u.)")
ax[0].set_ylabel("amplitude")
ax[0].legend()
ax[0].grid();


### variance plot ###
t, X1, X1sq, X2, X2sq = np.split(np.loadtxt("rf/squeeze.dat"), 5, axis=1)
Delta_X1sq = X1sq - X1**2
Delta_X2sq = X2sq - X2**2

ax[1].plot(t, Delta_X1sq, label="$\Delta X_1^2$", linestyle='-', linewidth=2.5)
ax[1].plot(t, Delta_X2sq, label="$\Delta X_2^2$", linestyle=':', linewidth=2.5, color='k')
ax[1].plot(t, Delta_X1sq*Delta_X2sq, label="$\Delta X_1^2\Delta X_2^2$",
           linestyle='--', linewidth=2.5)
ax[1].legend()
ax[1].set_xlabel("time (i.u.)")
ax[1].set_ylabel("magnitude")
ax[1].set_title("Optimized variances")
ax[1].set_yscale('log')
ax[1].grid();

Finally, we visualize the final state. Therefore, we load it from ‘psi_final_after_oct.dat’ and plot the population and phase of the different levels as well as the Wigner quasiprobability distribution in phase space using the plot_wigner function from qutip.

For a squeezed state, only even levels should be populated, which is already ensured by the two-photon driving. The phase of the coefficients then determines the angle of the squeezed state in phase space, i.e., which one of the quadratures is squeezed or a combination of both. For an \(\hat{X}_1\)-squeezed state, the coefficients should all be real with an alternating sign.

As one can see, this is the case for the first levels. However, the anharmonicity introduces an extra level-dependent phase , which is why the higher levels deviate from this pattern. This expresses itself also in the probability density in phase space. As required, the probability distribution of \(\hat{X}_1\) is very localized, i.e., squeezed while the other is very delocalized, which means that the variance of \(\hat{X}_2\) is very high. Due to the anharmonicity, the density is somewhat curved on the outer edges.

psi_final = read_psi_amplitudes("rf/psi_final_after_oct.dat", N_lvls, normalize=False)
levels    = np.arange(N_lvls)
pops      = np.abs(psi_final)
phases    = np.angle(psi_final)/(2*np.pi) + 1/2



fig, ax = plt.subplots(1, 2, figsize=(16, 6), width_ratios=[2,1])
fig.tight_layout()

### bar plot ###
my_cmap = plt.get_cmap('twilight')
colors = my_cmap(phases)
bars = ax[0].bar(levels, pops, color=colors, edgecolor='k')
ax[0].set_facecolor("white")

sm = ScalarMappable(cmap=my_cmap, norm=plt.Normalize(-1, 1))
sm.set_array([])

cbar = plt.colorbar(sm, ax=ax[0], pad=0.01)
cbar.set_label('Phase in units of $\pi$', rotation=270, labelpad=15)

ax[0].set_xlabel("level")
ax[0].set_ylabel("population")
ax[0].set_yscale('log')
ax[0].set_ylim(1e-3,1)
ax[0].set_xlim(-2, 61)
ax[0].set_title("Final state populations and phase");


### Wigner plot ###
qutip.plot_wigner(qutip.Qobj(psi_final), colorbar=True, fig=fig, ax=ax[1], )
ax[1].set_xlabel(r'$\hat{X}_1$', fontsize=15)
ax[1].set_ylabel(r'$\hat{X}_2$', fontsize=15)
ax[1].set_title("Wigner function");

When you have become familiar with this problem, try to implement the optimization of a squeezed state for an anharmonic oscillator, but with linear driving.