# Evaluating low-depth quantum algorithms for time evolution on fermion-boson
systems

Nathan Fitzpatrick<EMAIL_ADDRESS>Cambridge Quantum
Computing Ltd, 9a Bridge Street, Cambridge, United Kingdom  Harriet Apel
Cambridge Quantum Computing Ltd, 9a Bridge Street, Cambridge, United Kingdom
Department of Physics and Astronomy, University College London, Gower Street,
London, WC1E 6BT, United Kingdom  David Muñoz Ramo  Cambridge Quantum
Computing Ltd, 9a Bridge Street, Cambridge, United Kingdom

###### Abstract

Simulating time evolution of quantum systems is one of the most promising
applications of quantum computing and also appears as a subroutine in many
applications such as Green’s function methods. In the current era of NISQ
machines we assess the state of algorithms for simulating time dynamics with
limited resources. We propose the Jaynes-Cummings model and extensions to it
as useful toy models to investigate time evolution algorithms on near-term
quantum computers. Using these simple models, direct Trotterisation of the
time evolution operator produces deep circuits, requiring coherence times out
of reach on current NISQ hardware. Therefore we test two alternative responses
to this problem: variational compilation of the time evolution operator, and
variational quantum simulation of the wavefunction ansatz. We demonstrate
numerically to what extent these methods are successful in time evolving this
system. The costs in terms of circuit depth and number of measurements are
compared quantitatively, along with other drawbacks and advantages of each
method. We find that computational requirements for both methods make them
suitable for performing time evolution simulations of our models on NISQ
hardware. Our results also indicate that variational quantum compilation
produces more accurate results than variational quantum simulation, at the
cost of a larger number of measurements.

Quantum time evolution and Simulated time dynamics and Direct trotterisation
and Variational simulation and Incremental structured learning and Jaynes-
Cummings Model and Holstein–Primakoff encoding

##  I Introduction

Quantum computing has emerged as a promising technology to study quantum
dynamics problems, where the time evolution of a quantum system is essential
to understanding many complex physical processes. There has been much interest
in purely fermionic or bosonic dynamics on quantum computers, but little
interest has been shown to fermion-boson systems which are central to physical
processes such as electron-phonon interaction, quantum optics, and non-
adiabatic dynamics in chemical reactions Polyak  _et al._ (2019); Meyer and
Wang (2018); Spinlove  _et al._ (2018); Wang (2015); Richings  _et al._
(2015); Norell  _et al._ (2020). To our knowledge, only a scarce number of
works have been reported for the simulation of this interesting class of
systems Macridin  _et al._ (2018); Tacchino  _et al._ (2021).

The current quantum computing devices have small numbers of qubits with short
coherence times due to limited control Bharti  _et al._ (2021). Therefore,
Hamiltonian simulation frameworks that rely on many error corrected qubits
such as Phase Estimation and Quantum Signal Processing are not currently
possible Lee  _et al._ (2020); Kivlichan  _et al._ (2020); Dong  _et al._
(2020); Nam and Maslov (2019); Campbell (2019); Childs  _et al._ (2018);
Babbush  _et al._ (2018); Haah  _et al._ (2018). This situation has motivated
the development of many time evolution algorithms for near-term quantum
hardware. They can be classified into two broad families of algorithms:
variational compilation methods and wavefunction ansatz variational methods.
In the first case, the Trotterised time evolution operator is approximated by
a simplified circuit obtained via a variational optimisation procedure. In the
second case, time evolution is recast into an optimisation problem via an
appropriate variational principle, leading to the construction of the
wavefunction through an ansatz quantum circuit. Examples of the first family
include Variational Fast Forwarding Cîrstoiu  _et al._ (2020); Gibbs  _et al._
(2021), Incremental Structured Learning (ISL)Jaderberg  _et al._ , Variational
Time Dependent Phase Estimation Klymko  _et al._ (2021), circuit compilers for
time dependent applications Bassman  _et al._ (2020), Adaptive Product Formula
Zhang  _et al._ (2020), Quantum Imaginary Time Evolution Motta  _et al._
(2020) and truncated Dyson Series Kieferová  _et al._ (2019). Examples of the
second family of methods are Variational Quantum Simulation (VQS) Yuan  _et
al._ (2019); Mcardle  _et al._ (2019); McArdle  _et al._ (2019); Li and
Benjamin (2017), Hardware-Efficient Real Time Evolution Benedetti  _et al._
(2020), Quantum Assisted Simulation Lau  _et al._ (2021a) and Truncated Taylor
Series Lau  _et al._ (2021b).

The time evolution protocols mentioned above have been predominantly
demonstrated on purely fermionic systems like the hydrogen molecule or various
model Hamiltonians. However, their relative performance remains to be
understood, especially in the case of fermion-boson systems. As a particularly
appropriate instance of this class of problems, we propose the Jaynes-Cummings
Model (JCM) as a test case to benchmark the current NISQ quantum dynamics
methods. The JCM describes two-level oscillators interacting with a quantised
mode of an optical cavity in the semi-classical limit Jaynes and Cummings
(1963). This simple model can be represented with a two-qubit system and is a
convenient toy model since it produces periodic oscillations in time.
Extensions to the JCM that deal with an increasing number of oscillators
encoded in extra qubits can be easily considered. We shall test these model
systems with representative algorithms from the two families previously
described, namely the ISL and VQS methods. We will compare the relative merits
and flaws of the two methods at efficiently approximating simulations done
with direct Trotterisation of the time evolution operator.

The structure of this paper is as follows. Section II describes the theory of
the Jaynes-Cummings model and its extensions, laying out details of the
encoding used in the numerical simulations. Section III outlines the three
algorithms used for simulating time evolution. Section IV presents the results
for the three methods applied to the time evolution of our toy models,
including a detailed cost comparison. The final section provides our
conclusions regarding the relative performance of the time evolution
algorithms investigated in this work for our models.

##  II System Description

Dipolar coupling between a quantised electromagnetic field mode and a two-
level fermionic system is described by the quantum Rabi model Rabi (1936,
1937). The two-level system (hereafter referred to as ”the atom”) could indeed
be a natural atom, lending this model to atomic optics or instead an effective
system arising in a solid state device. In general the model describes the
simplest coherent coupling between a fermionic system and a bosonic mode or
photon. The quantum Rabi model Hamiltonian, $H_{R}$, is given by

$H_{R}=\frac{\hbar\omega_{c}}{2}\hat{\sigma}_{z}+\hbar\omega_{a}\hat{b}^{\dagger}\hat{b}+\hbar\Omega\hat{\sigma}_{x}(\hat{b}+\hat{b}^{\dagger}),$
(1)

where $\omega_{c}$, $\omega_{a}$ are resonance frequencies of the field mode
and atom respectively; $\Omega$ quantifies the strength of the coupling
between the field and the atom; $\hat{\sigma}_{z}$, $\hat{\sigma}_{x}$ are
Pauli operators acting on the atom and $\hat{b}^{\dagger}$, $\hat{b}$ are the
bosonic creation and annihilation operators.

Figure 1: The Jaynes-Cummings model of an atom coupled to a field mode in an
optical cavity. Excitation or relaxation between the two states is induced by
photon emission or absorption at frequency $\omega$ respectively.

###  II.1 Jaynes-Cummings Model

In order to make the quantum Rabi model Hamiltonian tractable for quantum time
evolution, various assumptions are then imposed. For ease, the system is
assumed to be at resonance, $\omega_{a}=\omega_{c}$ (although we briefly
explore off-resonant behaviour in Appendix A.1). The JCM arises by performing
the rotating wave approximation on the quantum Rabi model. A graphical
representation of this model is shown in Figure 1. Expressing the
$\hat{\sigma}_{x}$ in terms of Pauli ladder operators
($\hat{\sigma}_{x}=\hat{\sigma}_{+}+\hat{\sigma}_{-}$, where
$\hat{\sigma}_{\pm}=\sigma_{x}\pm i\sigma_{y}$), in this regime it is assumed
that there is strong coupling $\Omega/\omega<<1$ so the counter-rotating
terms, $\hat{a}^{\dagger}\hat{\sigma}_{+}$ and $\hat{a}\hat{\sigma}_{-}$,
contribute weakly to the dynamics and can be neglected. These approximations
leave the Hamiltonian as:

$H_{JC}=\hbar\omega(\hat{b}^{\dagger}\hat{b}+\frac{\hat{\sigma}_{z}}{2})+\frac{\hbar\Omega}{2}(\hat{b}\hat{\sigma}_{+}+\hat{b}^{\dagger}\hat{\sigma}_{-}).$
(2)

The interacting term conserves the system’s quanta through $U(1)$ symmetry and
only exchanges between the coupled subsystems are allowed; it is also
stationary in the interaction picture so the time dependence appears in the
state. Denoting the excited and ground state of the bosonic mode and the atom
as $\ket{e}/\ket{g}$ and $\ket{1}/\ket{0}$ respectively, consider a system
initially in state $\ket{\psi(t=0)}=\ket{g,1}$. At time $t$ the state will
evolve as:

$\ket{\psi(t)}=\cos\left(\frac{\Omega
t}{2}\right)\ket{g,1}-i\sin\left(\frac{\Omega t}{2}\right)\ket{e,0}.$ (3)

Since the system is closed, the total energy will be constant in time.
However, the probabilities of finding the system in a given eigenstate
oscillate in time e.g. $\mathbf{P}_{g,1}(t)=\cos^{2}\left(\frac{\Omega
t}{2}\right)$. This exchange of energy between the field and the atom can be
used to observe the system’s evolution in time by evaluating the expectation
value of the field and atom Hamiltonian at each time step.

###  II.2 Tavis-Cummings Model

An extension of the JCM comes from considering a single bosonic mode coupled
with multiple ($N>1$) atoms; this model is known as the Dicke or Tavis-
Cummings Model (TCM) Tavis and Cummings (1968). Physically this setup could
correspond to a higher dimensional spin or composite spin-$\frac{1}{2}$ system
coupled to the cavity. Increasing the number of quanta in the model grows the
Hilbert space exponentially, making it an ideal candidate for testing the
complexity of the quantum computational simulation achievable on NISQ devices.
The Hamiltonian for a system of $N$ atoms in the framework of this model is
given by:

$H_{TC}=\hbar\omega\hat{b}^{\dagger}\hat{b}+\hbar\omega\sum_{i=1}^{N}\frac{\hat{\sigma}_{z}^{i}}{2}-\frac{\hbar\Omega}{2}\sum_{i=1}^{N}(\hat{b}\hat{\sigma}_{+}^{i}+\hat{b}^{\dagger}\hat{\sigma}_{-}^{i}).$
(4)

As in the JCM case, this describes a closed system with a fixed number of
quanta.

###  II.3 Hamiltonian qubit encoding

General Hamiltonians contain both fermionic and bosonic operators. Encodings
translate Hamiltonians from the bosonic-fermionic representation into qubit
representations. The Jordan-Wigner encoding is a common method to transform
fermionic operators to quantum gates. However, as our numerical results
realise time evolution in a atomic-photon system, bosonic operators appear in
the Hamiltonian. Encoding bosonic operators for simulation on a quantum
computer is less common than the fermionic case, so the procedure followed in
order to obtain our qubit Hamiltonians is outlined here.

Using the Holstein–Primakoff transformation, we can relate spin operators to
bosonic creation and annihilation operators Holstein and Primakoff (1940):

$\hat{b}^{\dagger}\hat{b}=S+\hat{S}_{z},\qquad\hat{b}^{\dagger}=\frac{\hat{S}_{+}}{\sqrt{S-\hat{S}_{z}}},\qquad\hat{b}=\frac{\hat{S}_{-}}{\sqrt{S-\hat{S}_{z}}}.$
(5)

The non-trivial form of these expressions is greatly simplified by making the
semi-classical approximation, in which a truncated Taylor expansion in $1/S$
of these expressions is considered. This expansion is only valid for large
$S$, so this approximation would only be strict with a large number of qubits
representing the boson reservoir. However, this limit applied to the small
qubit register in our models captures a significant part of the essential
physics and is sufficient to demonstrate the time dynamics of the system. The
simplified transform is given by:

$\hat{b}^{\dagger}\approx\frac{\hat{S}_{+}}{\sqrt{2S}},\qquad\hat{b}\approx\frac{\hat{S}_{-}}{\sqrt{2S}}.$
(6)

In our models, we have one single qubit representing our boson reservoir, so
$S=1/2$ and the spin ladder operators can be easily translated into Pauli
terms.

This procedure is exemplified in the $N=1$ resonant JCM Hamiltonian with a
two-qubit register, where the 0 superscript in the operators corresponds to
the photon and the 1 superscript corresponds to the atom:

$\begin{split}H&=\omega(\hat{b}^{\dagger}\hat{b}+\hat{\sigma}_{z}/2)+\frac{\Omega}{2}(\hat{b}\hat{\sigma}_{+}+\hat{b}^{\dagger}\hat{\sigma}_{-})\\\
&=\omega(S+\frac{1}{2}\hat{\sigma}_{z}^{0}+\frac{1}{2}\hat{\sigma}_{z}^{1})+\frac{\Omega}{4\sqrt{2S}}[(\hat{\sigma}_{x}^{0}-i\hat{\sigma}_{y}^{0})(\hat{\sigma}_{x}^{1}+i\hat{\sigma}_{y}^{1})\\\
&+(\hat{\sigma}_{x}^{0}+i\hat{\sigma}_{y}^{0})(\hat{\sigma}_{x}^{1}-i\hat{\sigma}_{y}^{1})]\\\
&=\omega(S+\frac{1}{2}\hat{\sigma}_{z}^{0}+\frac{1}{2}\hat{\sigma}_{z}^{1})+\frac{\Omega}{4\sqrt{2S}}[2\hat{\sigma}_{x}^{0}\hat{\sigma}_{x}^{1}+2\hat{\sigma}_{y}^{0}\hat{\sigma}_{y}^{1}].\end{split}$
(7)

Encoding of the other Hamiltonians numerically simulated here follows
similarly.

##  III time evolution methods

The Hamiltonian determines the time evolution of quantum states. For a system
with a time independent Hamiltonian, $H$, an initial state, $\ket{\psi(t=0)}$,
can be propagated by the unitary time evolution operator $e^{-iHt}$. Unitary
operations can be performed on a quantum computer and, as previously
explained, different methods have been developed to implement the time
evolution operator on a quantum circuit. This section will briefly describe
three different representative methods, later used to obtain numerical
results.

###  III.1 Direct Trotterisation

Qubit Hamiltonians acting in $(\mathbb{C}^{2})^{\otimes n}$ can be represented
as a linear combination of tensor products of Pauli operators acting on
individual qubits:
$H=\sum_{k}c_{k}\hat{P}_{k}^{(1)}\otimes\hat{P}_{k}^{(2)}...\otimes\hat{P}_{k}^{(n)}$.
Direct exponentiation of this Hamiltonian is difficult, however quantum
simulation algorithms avoid this using the Trotter asymptotic approximation:

$\lim_{N\rightarrow\infty}\left(e^{iAt/N}e^{iBt/N}\right)^{N}=e^{i(A+B)t}.$
(8)

By performing a large but finite number of Trotter steps the unitary can be
well approximated to arbitrary precision at the cost of circuit depth. Using
this approximation the unitary time evolution operator becomes:

$\hat{U}(t)\approx\left(\sum_{k}e^{-ic_{k}\hat{P}_{k}^{(1)}\otimes...\otimes\hat{P}_{K}^{(n)}t/\hbar
N}\right)^{N}.$ (9)

Exponentiations of this kind are easily implemented in the circuit model.

###  III.2 Incremental structured learning

Linearly increasing circuit depth with time and accuracy is a significant
drawback of direct Trotterisation especially for near term applications.
Benedetti et al. Ostaszewski  _et al._ (2021) developed circuit structure
optimisation, later extended to the Incremental Structured Learning (ISL)
technique developed by Jaderberg et al.,Jaderberg  _et al._ which applies
variational compilation ideas inspired byJones and Benjamin (2018); Khatri
_et al._ (2019) to Trotterisation. In this technique, a so-called structural
ansatz is built up to approximate each unitary created by adding a Trotter
step to the time evolution operator. This optimisation is performed with
respect to minimising the cost function

$C(\vec{\theta})=1-|\bra{\bar{0}}\hat{Y}^{\dagger}(\vec{\theta})\hat{X}\ket{\bar{0}}|^{2},$
(10)

where the target state is engineered to always be the vacuum state
$\ket{\bar{0}}=\ket{0}^{\otimes n}$, $\hat{X}$ is the target unitary and
$\hat{Y}$ is the approximate structural ansatz unitary built up incrementally.
Adding the Trotter step to the previous structural ansatz creates the target
circuit for the next time step. This method circumvents the issue of linearly
increasing depth arising from direct Trotterisation by intelligently
restructuring a hardware efficient circuit at each time step. However, this
construction requires multiple evaluations of the cost function and introduces
a classical optimisation subroutine. This methodology is described in more
detail in Appendix B.

###  III.3 Variational Quantum Simulation

Wavefunction ansatz methods follow a different approach to the methods
previously discussed, and a good example is set out by the method of Li et al.
Li and Benjamin (2017), which uses a variational principle to simulate real
time evolution. Unlike the variational compilation methods, this method
requires a parameterised ansatz $\phi(\vec{\theta})$ to express the
wavefunction of the system. This is an initial disadvantage since postulating
a quantum circuit with sufficient freedom to represent the system at all times
is non-trivial and often requires prior knowledge of the quantum system. For
our numerical results, using a generalised or hardware-efficient ansatz was
sufficient but this constraint could pose a challenge for general use of this
technique. Instead of directly solving the Schrödinger equation, McLachlan’s
variational principle is applied to represent the evolution of the
wavefunction parameters with time. This evolution is given by the expression:

$\sum_{i}A^{R}_{ij}\dot{\theta}_{j}=C^{I}_{i},$ (11)

for time evolution governed by the Hamiltonian in its Pauli decomposition
$H=\sum_{k}c_{k}\hat{P}_{k}$. Given an ansatz with $n$ parameters,
$\mathbf{A}^{R}$ is a $(n\times n)$ matrix and $\mathbf{C}^{I}$ is a vector of
size $n$ defined as:

$A^{R}_{ij}=\Re\left(\frac{\partial\bra{\phi(t)}}{\partial\theta_{i}}\frac{\partial\ket{\phi(t)}}{\partial\theta_{j}}\right)$
(12)
$C^{I}_{i}=\Im\left(\sum_{k}c_{k}\frac{\partial\bra{\phi(t)}}{\partial\theta_{i}}\hat{P}_{k}\ket{\phi(t)}\right).$
(13)

We show the quantum circuits required to compute these objects for our models
in Appendix C. The nature of this method requires multiple circuit evaluations
to populate $\mathbf{A}^{R}$ and $\mathbf{C}^{I}$, whose size is heavily
dependent on having an efficient ansatz with few parameters. However, the
depth of these circuits is constant throughout the time evolution so trading
numerous measurements for low depth circuits may be particularly advantageous
in the NISQ era. Once these objects are computed, the parameters are advanced
using an update method:

$\vec{\theta}(t+\delta t)=\vec{\theta}(t)+\vec{\dot{\theta}}(t)\cdot\delta t.$
(14)

A more detailed discussion of this algorithm’s features is provided in
Appendix C.

##  IV Numerical results

We present here our results for simulations of time evolution of the different
models introduced previously with direct Trotterisation and the two
variational algorithms previously described. The simulations have been run
with our EUMEN molecular simulation platform,eum  interfaced with the tket
package for circuit optimisation Sivarajah  _et al._ (2020) and the QiskitAer
shot based simulator Abraham  _et al._ (2019). Parameters for our Hamiltonians
have been set to $S=1/2$, $\omega=2$ and $\Omega=10$.

###  IV.1 JCM results

The JCM on two qubits was simulated evolving in time using the three methods
described above. Figure 2 shows the energy evolution at each timestep for each
method, along with a representation of the error in the energy estimation. All
the methods agree qualitatively with the exact result calculated from the
model using $\mathbf{P}_{g,1}(t)=\cos^{2}\left(\frac{\Omega t}{2}\right)$,
validating the use of variational methods to study the JCM system. The error
at this stage is dominated by Trotterisation error from discretising the
system and so is similar for all three methods. This error can be decreased by
increasing the number of Trotter steps and decreasing $\delta t$ in VQS. Since
this discretising is causing a dephasing, it is intuitive that error is
periodic like the exact function and is minimised where the function’s
gradient is minimum (at extrema of energies). When averaged over a period, the
error is increasing with the number of time steps.

###  IV.2 TCM results

To demonstrate scaling, the TCM on three qubits (system with two quanta) was
simulated evolving in time using the three methods described above. Results
are displayed in Figure 3. The different methods now produce varying accuracy
and are compared with the analytical result calculated from the model using
$\mathbf{P}_{g,1}(t)=\cos^{2}\left(\frac{\sqrt{n}\Omega t}{2}\right)$. As in
the JCM case there is dephasing due to discretisation. However, VQS is now
underperforming with a generally higher error that the other two methods. This
could be due to the ansatz restricting the Hilbert space, demonstrating the
need of careful consideration of the ansatz when using these methods. In
contrast, ISL accuracy is comparable to direct Trotterisation. Increasing the
number of layers in the ansatz may be a good solution to increase VQS
accuracy. However, circuit depth and number of parameters to optimise also
increase significantly. As current NISQ hardware can only manipulate a limited
number of two-qubit gates before noise effects become untractable, this
solution is currently unrealistic. In addition, the extra number of parameters
may introduce additional complications to the ansatz optimisation process.

Figure 2: Left: Simulation of time evolution of the JCM system by three
methods: direct Trotterisation, VQS and ISL. Right: Absolute error in the
dynamics with respect to time, for both the atom and the whole system. The
methods simulate the time evolution to comparable accuracy, indicating error
is originating from discretisation for this model.

Figure 3: Left: Simulation of time evolution of the TCM system ($n=2$) by
three methods: direct Trotterisation, VQS and ISL. Right: Absolute error in
the dynamics with respect to time. Unlike the JCM case, in this more complex
system the methods’ performance is varied with VQS being less accurate.

###  IV.3 Cost comparison

#### Circuit depth

Due to the short decoherence times and high noise levels of current devices,
deep circuits cannot be meaningfully measured. When simulating time evolution
with all of the methods reported in this work, a discrete time step $\delta t$
is chosen, and a coarser discretisation results in a less accurate simulation.
The number of time/Trotter steps is given by $m=T/\delta t$ where $T$ is the
total evolution time. In all the results presented here $\delta t=0.01$.
Figure 4 demonstrates how the circuit depths vary for each method in the JCM
and TCM simulations. The linear increase with time steps for direct
Trotterisation highlights the main disadvantage of this method. VQS has a
constant depth determined by the choice of ansatz. ISL depth is also shown to
be approximately constant and to predominately outperform VQS in creating
shallow circuits with fewer two-qubit entangling gates. In this metric ISL and
VQS clearly outperform direct Trotterisation demonstrating the motivation for
their conception, however this advantage comes at other costs examined in the
next section.

Figure 4: Top row: Circuit depth comparison for the JCM simulations. Bottom
row: Circuit depth comparison for the TCM simulations. In both cases both
total circuit depth and 2-qubit gate depth are shown. Direct Trotterisation
displays a linear depth dependence whereas ISL and VQS show modest, almost
constant, circuit depths.

#### Number of circuit evaluations

Multiple shallow circuit evaluations, while feasible on NISQ hardware with
manageable noise, generate large computational overheads. The more shots
necessary to obtain the result, the more runtime is needed on the device, so
it is desirable to minimise the number of circuit evaluations. Consider the
scenario where one takes $k$ shots to produce an expectation value for a
circuit, the system Hamiltonian is encoded as a sum of $n_{p}$ Pauli strings
and the ansatz used for the VQS method contains $n_{pm}$ variable parameters.
For ISL, the number of layers used (nl) and the number of cost evaluations per
layer (nev) must be considered as well, taking into account that nev is
dependent on the efficiency of the chosen optimisation strategy. The number of
quantum circuit evaluations $c$ required per time step for each method scales
as:

$\displaystyle\text{DT:}\qquad c\approx kn_{p}$ (15)
$\displaystyle\text{ISL:}\qquad c\approx kn_{p}+n_{l}\times n_{ev}\times k$
(16) $\displaystyle\text{VQS:}\qquad c\approx
kn_{p}+\frac{1}{2}kn_{pm}(n_{pm}+1)+n_{pm}pk$ (17)

For the test case of the 3-qubit TCM Hamiltonian and a number of evaluations
$k=1000$, we have $n_{p}=7$, $n_{pm}=18$, average $n_{ev}=919$ and average
$n_{l}=5.5$. The different $c$ values for each method can be found in Table 1.
This test case gives an insight into the scale of measurements required for
the different methods, however this cost is strongly model dependent. It is
expected that the number of layers needed for ISL will increase as the number
of qubits required by the model increases. Similarly, as larger systems are
simulated, the ability to find an appropriate ansatz becomes less trivial and
the number of variational parameters increases. How these two models compare
for other larger models is uncertain, but in this test case VQS performed
better with fewer circuit evaluations. Variational simulation performs
multiple measurements but over 3 fixed circuit structures with variational
parameters, whereas ISL, by its nature, requires measurements of many
different circuit structures. At implementation level this is a disadvantage
of ISL.

Table 1: Comparison of methods for time evolution on a quantum computer. The circuit evaluation is specific to the test case of TCM model with N=2, where $m$ is the number of time/Trotter steps. Method | Classical processing | Ansatz | Circuit Depth scaling | Circuit evaluations/time step
---|---|---|---|---
Trotter | None | No | $\mathcal{O}(m)$ | 7000
ISL | Optimisation | No | $\mathcal{O}(1)$ | 5061500
VQS | Tikhonov regularisation | Yes | $\mathcal{O}(1)$ | 304000

#### Ansatz

VQS is unique within our tested algorithms in that it requires a fixed ansatz
to express the wavefunction. The primary drawback of this feature is that the
additional work required is specific to each system you want to simulate. This
complication can be avoided by using hardware efficient ansätze. However, it
is not guaranteed that this type of ansatz will have sufficient degrees of
freedom to describe the system as it evolves. However, an informed choice of
ansatz that is physically motivated by the properties of the system may
increase the robustness of the simulation; for example by limiting the
associated Hilbert space to one which preserves a known symmetry.

#### Classical routines

Both ISL and VQS introduce classical subroutines. ISL requires classical
optimisation since each layer includes variational parameters which are
optimised to minimise the cost function. During our numerical simulations the
simplistic Nelder-Mead method was used and more sophisticated algorithms,
including conjugate gradient and BFGS, failed indicating barren plateaus on
the cost function surface. However since in ISL the optimisation is iterative,
there is freedom independent of system size to tailor the implementation of
the algorithm to keep the number of variational parameters manageable,
indicating that this method should scale efficiently. In contrast VQS avoids
optimisation but instead requires the inversion of the matrix $\mathbf{A}^{R}$
in order to solve the differential equation 11. For our simulation, Tikhonov
regularisation is used as this is robust when the matrix is close to singular.
Unfortunately, the matrix size scales with the number of variational
parameters in the chosen ansatz. Therefore, more complex systems requiring
more qubits for their description may strain this subroutine. The balance of
these costs will determine which method is most suited to a given task on a
given machine.

##  V Discussion and Conclusions

Time evolution in the JCM and TCM was successfully simulated using the
Holstein-Primakov tranform on a quantum computer, and the expected Rabi
oscillations were successfully reproduced. The numerical results demonstrate
that Trotterisation of the time evolution operator is an accurate method
provided an appropriate time step is considered, although it is not friendly
to current NISQ technology. Variational methods, either through circuit
compilation or ansatz optimisation, are viable for simulating time dynamics on
NISQ hardware. All of them can qualitatively, and in some cases
quantitatively, describe features of these models. Despite the problems
reported, these methods can be used to explore interesting physical phenomena
appearing in both the JCM and TCM. We show in Appendix A some examples that
illustrate this point.

Our benchmarking uncovered no overall superior time evolution simulation
algorithm: as expected, the best method will depend on the system of interest
and resources available. Direct Trotterisation is the simplest to implement
and for small systems and short times the circuit depth may be manageable. In
the current regime of NISQ machines, the alternative methods lend themselves
to longer time evolutions and provide another option when decoherence times
prevent the deep circuits generated by Trotterisation from being realised. In
this case the choice between VQS and ISL again lies in the particulars. If an
suitable ansatz can be physically motivated with few parameters, VQS is
efficient, requiring fewer circuit evaluations and with an attractively
prescriptive methodology. However, ISL prevails in situations where the system
is less understood and iteration to tolerance ensures the ansatz has
sufficient freedom to accurately describe the unitary. There is also freedom
within the ISL method regarding the optimisation protocol and construction of
the approximate unitary that could be tailored to be efficient for the given
problem.

Despite the simplifications introduced, computational requirements for the two
variational methods are still challenging for current quantum computers.
Further work is required to improve accuracy and reduce costs, including
number of measurements and circuit depth, in order to make simulations of more
complex systems feasible for NISQ devices.

## Acknowledgements

The authors thank Seyon Sivarajah for his helpful discussions on constructing
pytket quantum circuits.

## References

  * Polyak  _et al._ (2019) I. Polyak, G. W. Richings, S. Habershon, and P. J. Knowles, J. Chem. Phys 150, 41101 (2019).
  * Meyer and Wang (2018) H. D. Meyer and H. Wang, Journal of Chemical Physics 148, 124105 (2018).
  * Spinlove  _et al._ (2018) K. E. Spinlove, G. W. Richings, M. A. Robb, and G. A. Worth, Faraday Discussions 212, 191 (2018).
  * Wang (2015) H. Wang, Journal of Physical Chemistry A 119, 7951 (2015).
  * Richings  _et al._ (2015) G. W. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt, and B. Lasorne, International Reviews in Physical Chemistry 34, 269 (2015).
  * Norell  _et al._ (2020) J. Norell, M. Odelius, and M. Vacher, The Journal of Chemical Physics 7, 244115 (2020).
  * Macridin  _et al._ (2018) A. Macridin, P. Spentzouris, J. Amundson, and R. Harnik, Physical Review A 98, 042312 (2018).
  * Tacchino  _et al._ (2021) F. Tacchino, A. Chiesa, R. Sessoli, I. Tavernelli, and S. Carretta,  (2021), arXiv:2103.09706 .
  * Bharti  _et al._ (2021) K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik,  (2021), arXiv:2101.08448 .
  * Lee  _et al._ (2020) J. Lee, D. W. Berry, C. Gidney, W. J. Huggins, J. R. McClean, N. Wiebe, and R. Babbush,  (2020), arXiv:2011.03494 .
  * Kivlichan  _et al._ (2020) I. D. Kivlichan, C. Gidney, D. W. Berry, N. Wiebe, J. McClean, W. Sun, Z. Jiang, N. Rubin, A. Fowler, A. Aspuru-Guzik, H. Neven, and R. Babbush, Quantum 4, 296 (2020).
  * Dong  _et al._ (2020) Y. Dong, X. Meng, K. B. Whaley, and L. Lin,  _arXiv_ ,  (2020), arXiv:2002.11649 .
  * Nam and Maslov (2019) Y. Nam and D. Maslov, npj Quantum Information 5, 44 (2019).
  * Campbell (2019) E. Campbell, Physical Review Letters 123, 070503 (2019).
  * Childs  _et al._ (2018) A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, Proceedings of the National Academy of Sciences of the United States of America 115, 9456 (2018).
  * Babbush  _et al._ (2018) R. Babbush, N. Wiebe, J. Mcclean, J. Mcclain, H. Neven, G. Kin, and L. Chan, Physical Review X 8, 011044 (2018).
  * Haah  _et al._ (2018) J. Haah, M. B. Hastings, R. Kothari, and G. H. Low,  (2018), arXiv:1801.03922 .
  * Cîrstoiu  _et al._ (2020) C. Cîrstoiu, Z. Holmes, J. Iosue, L. Cincio, P. J. Coles, and A. Sornborger, npj Quantum Information 6, 82 (2020).
  * Gibbs  _et al._ (2021) J. Gibbs, K. Gili, Z. Holmes, B. Commeau, A. Arrasmith, L. Cincio, P. J. Coles, and A. Sornborger,  (2021), arXiv:2102.04313 .
  * (20) B. Jaderberg, A. Agarwal, K. Leonhardt, M. Kiffner, and D. Jaksch, Quantum Science and Technology 5, 3.
  * Klymko  _et al._ (2021) K. Klymko, C. Mejuto-Zaera, S. J. Cotton, F. Wudarski, M. Urbanek, D. Hait, M. Head-Gordon, K. B. Whaley, J. Moussa, N. Wiebe, W. A. de Jong, and T. N. M.,  (2021), arXiv:2103.08563 .
  * Bassman  _et al._ (2020) L. Bassman, S. Gulania, C. Powers, R. Li, T. Linker, K. Liu, T. K. Kumar, R. K. Kalia, A. Nakano, and P. Vashishta,  (2020), arXiv:2004.07418 .
  * Zhang  _et al._ (2020) Z. J. Zhang, J. Sun, X. Yuan, and M. H. Yung,  (2020), arXiv:2011.05283 .
  * Motta  _et al._ (2020) M. Motta, C. Sun, A. T. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. Brandão, and G. K. L. Chan, Nature Physics 16, 205 (2020).
  * Kieferová  _et al._ (2019) M. Kieferová, A. Scherer, and D. W. Berry, Physical Review A 99, 042314 (2019).
  * Yuan  _et al._ (2019) X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin, Quantum 3, 191 (2019).
  * Mcardle  _et al._ (2019) S. Mcardle, X. Yuan, and S. Benjamin,  (2019), arXiv:1807.02467v5 .
  * McArdle  _et al._ (2019) S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, npj Quantum Information 5, 75 (2019).
  * Li and Benjamin (2017) Y. Li and S. C. Benjamin, Physical Review X 7, 021050 (2017).
  * Benedetti  _et al._ (2020) M. Benedetti, M. Fiorentini, and M. Lubasch,  (2020), arXiv:2009.12361v2 .
  * Lau  _et al._ (2021a) J. W. Z. Lau, K. Bharti, T. Haug, and L. C. Kwek,  (2021a), arXiv:2101.07677 .
  * Lau  _et al._ (2021b) J. W. Z. Lau, T. Haug, L. C. Kwek, and K. Bharti,  (2021b), arXiv:2103.05500 .
  * Jaynes and Cummings (1963) E. Jaynes and F. Cummings, Proceedings of the IEEE 51, 89 (1963).
  * Rabi (1936) I. I. Rabi, Physical Review 49, 324 (1936).
  * Rabi (1937) I. I. Rabi, Physical Review 51, 652 (1937).
  * Tavis and Cummings (1968) M. Tavis and F. W. Cummings, Phys. Rev. 170, 379 (1968).
  * Holstein and Primakoff (1940) T. Holstein and H. Primakoff, Physical Review 58, 1098 (1940).
  * Ostaszewski  _et al._ (2021) M. Ostaszewski, E. Grant, and M. Benedetti, Quantum 5, 391 (2021).
  * Jones and Benjamin (2018) T. Jones and S. C. Benjamin,  (2018), arXiv:1811.03147 .
  * Khatri  _et al._ (2019) S. Khatri, R. LaRose, A. Poremba, L. Cincio, A. T. Sornborger, and P. J. Coles, Quantum 3, 140 (2019).
  * (41) “Eumen: Quantum chemistry and materials on quantum computers,” .
  * Sivarajah  _et al._ (2020) S. Sivarajah, S. Dilkes, A. Cowtan, W. Simmons, A. Edgington, and R. Duncan, Quantum Science and Technology 6, 014003 (2020).
  * Abraham  _et al._ (2019) H. Abraham, AduOffei, R. Agarwal, I. Y. Akhalwaya, G. Aleksandrowicz, T. Alexander, M. Amy, E. Arbel, Arijit02, A. Asfaw, A. Avkhadiev, C. Azaustre, AzizNgoueya, A. Banerjee, A. Bansal, P. Barkoutsos, A. Barnawal, G. Barron, G. S. Barron, L. Bello, Y. Ben-Haim, D. Bevenius, A. Bhobe, L. S. Bishop, C. Blank, S. Bolos, S. Bosch, Brandon, S. Bravyi, Bryce-Fuller, D. Bucher, A. Burov, F. Cabrera, P. Calpin, L. Capelluto, J. Carballo, G. Carrascal, A. Chen, C.-F. Chen, E. Chen, J. C. Chen, R. Chen, J. M. Chow, S. Churchill, C. Claus, C. Clauss, R. Cocking, F. Correa, A. J. Cross, A. W. Cross, S. Cross, J. Cruz-Benito, C. Culver, A. D. Córcoles-Gonzales, S. Dague, T. E. Dandachi, M. Daniels, M. Dartiailh, DavideFrr, A. R. Davila, A. Dekusar, D. Ding, J. Doi, E. Drechsler, Drew, E. Dumitrescu, K. Dumon, I. Duran, K. EL-Safty, E. Eastman, G. Eberle, P. Eendebak, D. Egger, M. Everitt, P. M. Fernández, A. H. Ferrera, R. Fouilland, FranckChevallier, A. Frisch, A. Fuhrer, B. Fuller, M. GEORGE, J. Gacon, B. G. Gago, C. Gambella, J. M. Gambetta, A. Gammanpila, L. Garcia, T. Garg, S. Garion, A. Gilliam, A. Giridharan, J. Gomez-Mosquera, Gonzalo, S. de la Puente González, J. Gorzinski, I. Gould, D. Greenberg, D. Grinko, W. Guan, J. A. Gunnels, M. Haglund, I. Haide, I. Hamamura, O. C. Hamido, F. Harkins, V. Havlicek, J. Hellmers, Ł. Herok, S. Hillmich, H. Horii, C. Howington, S. Hu, W. Hu, J. Huang, R. Huisman, H. Imai, T. Imamichi, K. Ishizaki, R. Iten, T. Itoko, JamesSeaward, A. Javadi, A. Javadi-Abhari, W. Javed, Jessica, M. Jivrajani, K. Johns, S. Johnstun, Jonathan-Shoemaker, V. K, T. Kachmann, A. Kale, N. Kanazawa, Kang-Bae, A. Karazeev, P. Kassebaum, J. Kelso, S. King, Knabberjoe, Y. Kobayashi, A. Kovyrshin, R. Krishnakumar, V. Krishnan, K. Krsulich, P. Kumkar, G. Kus, R. LaRose, E. Lacal, R. Lambert, J. Lapeyre, J. Latone, S. Lawrence, C. Lee, G. Li, D. Liu, P. Liu, Y. Maeng, K. Majmudar, A. Malyshev, J. Manela, J. Marecek, M. Marques, D. Maslov, D. Mathews, A. Matsuo, D. T. McClure, C. McGarry, D. McKay, D. McPherson, S. Meesala, T. Metcalfe, M. Mevissen, A. Meyer, A. Mezzacapo, R. Midha, Z. Minev, A. Mitchell, N. Moll, J. Montanez, G. Monteiro, M. D. Mooring, R. Morales, N. Moran, M. Motta, MrF, P. Murali, J. Müggenburg, D. Nadlinger, K. Nakanishi, G. Nannicini, P. Nation, E. Navarro, Y. Naveh, S. W. Neagle, P. Neuweiler, J. Nicander, P. Niroula, H. Norlen, NuoWenLei, L. J. O’Riordan, O. Ogunbayo, P. Ollitrault, R. Otaolea, S. Oud, D. Padilha, H. Paik, S. Pal, Y. Pang, V. R. Pascuzzi, S. Perriello, A. Phan, F. Piro, M. Pistoia, C. Piveteau, P. Pocreau, A. Pozas-Kerstjens, M. Prokop, V. Prutyanov, D. Puzzuoli, J. Pérez, Quintiii, R. I. Rahman, A. Raja, N. Ramagiri, A. Rao, R. Raymond, R. M.-C. Redondo, M. Reuter, J. Rice, M. Riedemann, M. L. Rocca, D. M. Rodríguez, RohithKarur, M. Rossmannek, M. Ryu, T. SAPV, SamFerracin, M. Sandberg, H. Sandesara, R. Sapra, H. Sargsyan, A. Sarkar, N. Sathaye, B. Schmitt, C. Schnabel, Z. Schoenfeld, T. L. Scholten, E. Schoute, J. Schwarm, I. F. Sertage, K. Setia, N. Shammah, Y. Shi, A. Silva, A. Simonetto, N. Singstock, Y. Siraichi, I. Sitdikov, S. Sivarajah, M. B. Sletfjerding, J. A. Smolin, M. Soeken, I. O. Sokolov, I. Sokolov, SooluThomas, Starfish, D. Steenken, M. Stypulkoski, S. Sun, K. J. Sung, H. Takahashi, T. Takawale, I. Tavernelli, C. Taylor, P. Taylour, S. Thomas, M. Tillet, M. Tod, M. Tomasik, E. de la Torre, K. Trabing, M. Treinish, TrishaPe, D. Tulsi, W. Turner, Y. Vaknin, C. R. Valcarce, F. Varchon, A. C. Vazquez, V. Villar, D. Vogt-Lee, C. Vuillot, J. Weaver, J. Weidenfeller, R. Wieczorek, J. A. Wildstrom, E. Winston, J. J. Woehr, S. Woerner, R. Woo, C. J. Wood, R. Wood, S. Wood, S. Wood, J. Wootton, D. Yeralin, D. Yonge-Mallo, R. Young, J. Yu, C. Zachow, L. Zdanski, H. Zhang, C. Zoufal, Zoufalc, a kapila, a matsuo, bcamorrison, brandhsn, nick bronn, brosand, chlorophyll zz, csseifms, dekel.meirom, dekelmeirom, dekool, dime10, drholmie, dtrenev, ehchen, elfrocampeador, faisaldebouni, fanizzamarco, gabrieleagl, gadial, galeinston, georgios ts, gruu, hhorii, hykavitha, jagunther, jliu45, jscott2, kanejess, klinvill, krutik2966, kurarrr, lerongil, ma5x, merav aharoni, michelle4654, ordmoj, sagar pahwa, rmoyard, saswati qiskit, scottkelso, sethmerkel, shaashwat, sternparky, strickroman, sumitpuri, tigerjack, toural, tsura crisaldo, vvilpas, welien, willhbang, yang.luh, yotamvakninibm, and M. Čepulkovskis, “Qiskit: An open-source framework for quantum computing,”  (2019).
  * Mclachlan (1964) A. D. Mclachlan, Molecular Physics 8, 39 (1964).

##  Appendix A Exploring model features

In our simulations we have focused on the basics of the model in order to
compare the different methods for simulating time evolution. However, the
Jaynes-Cummings model is interesting in its own right and by varying
parameters, the aspects of the system’s physics can be demonstrated using the
time evolution methods described in this work. Obtaining the expected results
is also an important sanity check that our computation is successfully
simulating the physical system.

###  A.1 Off resonant behaviour

The main numerical results of this paper are developed for resonant Rabi
oscillations. The cavity and oscillators can be detuned by adjusting the
parameters of the toy model which generates an additional term in the
Hamiltonian:

$H=\hbar\omega_{c}(S+\frac{1}{2}\hat{\sigma}^{0}_{z})+\frac{\hbar\omega_{a}}{2}\hat{\sigma}^{1}_{z}+\frac{\hbar(\omega_{a}-\omega_{c})}{2}\hat{\sigma}^{1}_{z}+\frac{\hbar\Omega}{4\sqrt{2S}}\left[2\hat{\sigma}^{0}_{x}\hat{\sigma}^{1}_{x}+2\hat{\sigma}^{0}_{y}\hat{\sigma}^{1}_{y}\right].$
(18)

The system, when evolved in time under this modified Hamiltonian, performed
Rabi oscillations at a slightly different frequency and reduced amplitude.
Figure 5 demonstrates how detuning changes the state evolution by comparing
the on and off resonance behaviour of the field energy. This behaviour is
expected theoretically as detuning prevents the subsystem from transitioning
between the two sublevels.

Figure 5: Field energy evolution with time for the JCM at on and off resonance
regimes.

###  A.2 Investigating number of quanta

Our results focused on small systems of 2 or 3 qubits but the time evolution
methods are all general to systems of any size. The toy model can also be
adapted to demonstrate larger systems by adding quanta in the TCM which
corresponds to adding qubits to the register. Theoretically, the Rabi
oscillation associated with an atom-field system of $n$ quanta will have
frequency:

$\Omega_{n}=\frac{\Omega\sqrt{n}}{2}.$ (19)

A manifestation of this relationship is seen numerically in Figure 6, where
the ratio of gradients agree approximately with what is expected from equation
19:

$\displaystyle\frac{0.465}{0.332}$
$\displaystyle=1.40\approx\frac{\sqrt{2}}{\sqrt{1}}=1.41$ (20)
$\displaystyle\frac{0.549}{0.332}$
$\displaystyle=1.66\approx\frac{\sqrt{3}}{\sqrt{1}}=1.73.$ (21)

Figure 6: Evolution of oscillation frequency with Rabi frequency calculated
with trotterised time evolution.

##  Appendix B ISL method

We provide here a detailed overview of the ISL method, starting with the flow
of the algorithm, then providing information about the optimisation steps and
the qubit selection process implemented here. For a more in-depth explanation
of the algorithm, we refer the reader to the work of Jaderberg et al.
Jaderberg  _et al._ whose method we follow closely.

###  B.1 Algorithm

We start with the state of the system at time $t$; this state could be just
the initialisation of the qubits or the circuit produced from the previous
time step. The state is propagated in time by appending the usual Trotter step
giving the target unitary $\hat{X}$. Instead of directly measuring properties
of $\hat{X}$, an approximate unitary is built that describes with arbitrary
accuracy the state of the system at time $t+\delta t$.

This approximate unitary, $\hat{Y}$, is constructed by optimising a circuit
with variational parameters $\vec{\theta}$. The optimisation is done with
respect to minimising the cost function
$C(\vec{\theta})=1-|\bra{\bar{0}}\hat{Y}^{\dagger}(\vec{\theta})\hat{X}\ket{\bar{0}}|^{2}$,
which has a minimum when $\hat{X}\hat{Y}^{\dagger}=\mathbb{I}$. The unitary
$\hat{Y}$ is then obtained from $\hat{Y}^{\dagger}$. The low depth $\hat{Y}$
circuit is used to evaluate any quantities for the model (i.e. in this case
the system, atom and field energies) before another Trotter step is appended
and the system is repeated to continue propagating the state in time whilst
avoiding the linear increase in depth.

###  B.2 Optimisation of the structural ansatz

There is freedom in how $\hat{Y}^{\dagger}$ is constructed, which allows for
flexibility depending on the available resources. A suitable approach is
choosing a structural ansatz that is then optimised; this would greatly reduce
the number of measurements required while introducing other issues. Choosing a
suitable structural ansatz with enough freedom to describe the system is non-
trivial for large systems, and the state’s evolution quickly deviates from the
true system dynamics if the approximated unitary is unable to drive the cost
function below a strict tolerance value. To avoid these problems the ISL
method constructs $\hat{Y}^{\dagger}$ by incrementally adding small circuit
blocks. For our implementation the circuit block acts on two qubits and is a
CNOT gate surrounded by single qubit rotation gates (see Figure 7).

Figure 7: Nested CNOT circuit block used to incrementally grow the ansatz to
optimse $\hat{B}^{\dagger}$.

In order to minimise the circuit depth each layer is optimised carefully:

  1. 1.

two qubits are selected to which a new circuit block is appended;

  2. 2.

in this circuit block the rotation angle and axis of three rotations is fixed
while the fourth is optimised. This process is done sequentially for all four
rotations in the block before the new cost function is evaluated. This step is
then repeated until the improvement in cost between cycles is below certain
threshold (1% in this work);

  3. 3.

fixing all rotation axes, all angles in the structural ansatz are optimised
together before a final tket circuit optimisation is applied to remove
redundancies;

  4. 4.

if the cost is below the pre-defined tolerance level, the cycle terminates. If
it is above the tolerance level, the whole process is repeated.

A scheme of the compilation process is shown in Figure 8. For the numerical
results a tolerance of $1\times 10^{-4}$ was used.

Figure 8: A circuit description of the ISL update procedure evaluating the
cost function $|\bra{0}\hat{Y}^{\dagger}(\vec{\theta})\hat{X}\ket{0}|^{2}$ for
an individual trotter step denoted by the circuit $U_{a+1}$, where the
intermediate compiled circuit is given by $A_{a+1:0}$. Nested CNOTs will
continue to be appended until a pre-defined tolerance is reached.

###  B.3 Choosing qubits

Since the quantum system being simulated will, in general, require $n>2$
qubits to span the Hilbert space, careful placement of the 2 qubit circuit
block at each augmentation of $\hat{Y}^{\dagger}$ has the potential to
increase method performance. Selecting the pair of qubits with the highest
entanglement outlined is the approached employed here. The entropy of
formation of each pair of qubits is calculated by tracing out all other qubits
and then calculating the entropy of formation, $E$. The qubit pair with the
highest $E$ is selected conditional on it being a different pair to the
previous circuit block.

To calculate the entropy of formation of 2 qubits in a larger circuit the
procedure in Ostaszewski  _et al._ (2021) and Jaderberg  _et al._ was
followed. In summary:

  1. 1.

the full density operator is obtained from the circuit;

  2. 2.

all but the two qubits of interest are traced out leaving a ($4\times 4$)
reduced density operator $\rho$;

  3. 3.

the concurrence $C(\rho)$ is calculated

$C(\rho)=\text{max}\\{0,2\lambda_{\text{max}}-\text{Tr}(R)\\},$ (22)

where $\lambda_{\text{max}}$ is the largest eigenvalue of $\rho$.
$R=\sqrt{\sqrt{\rho}\tilde{\rho}\sqrt{\rho}}$ and $\tilde{\rho}$ is the
Hermitian conjugate of $\rho$ expressed in the magic basis:

$\begin{split}\ket{e_{1}}&=\frac{1}{2}(\ket{00}+\ket{11}),\qquad\ket{e_{2}}=\frac{1}{2}(\ket{00}-\ket{11}),\\\
\ket{e_{3}}&=\frac{1}{2}(\ket{01}+\ket{10}),\qquad\ket{e_{4}}=\frac{1}{2}(\ket{01}-\ket{10});\end{split}$
(23)

  4. 4.

the entropy function is calculated as

$E(p)=p\log_{2}(p)+(1-p)\log_{2}(1-p)$ (24)

where $p(\rho)=\frac{1}{2}+\frac{1}{2}\sqrt{1-C(\rho)}$.

##  Appendix C VQS method

In this appendix, we provide some additional background on the variational
quantum simulation method. For full details, we refer to the original work in
Refs. Yuan  _et al._ (2019); Mcardle  _et al._ (2019); McArdle  _et al._
(2019); Li and Benjamin (2017). We first describe the method to recast the
time evolution operator into a variational problem, and then we provide
information about the wavefunction ansatz used for this method.

###  C.1 Variational Principle

We express the system’s wavefunction as a series of unitary operators acting
on the $\ket{\overline{0}}$ state:

$\ket{\psi}=U_{N}...U_{j}U_{j-1}...U_{1}\ket{\overline{0}}=V\ket{\overline{0}}$
(25)

If we allow $U_{j}$ to be a parametrised gate represented in terms of a Pauli
generators in our state preparation ansatz

$U_{j}=e^{-i\sum_{a}\theta_{j}\sigma_{a}/2}$ (26)

Using McLachlan’s variational principle, Mclachlan (1964) the real time
evolution of the wavefunction is given by:

$\sum_{i}A^{R}_{ij}\dot{\theta}_{j}=C^{I}_{i}$ (27)

After calculating $\mathbf{A}^{R}$ and $\mathbf{c}^{I}$ the parameters are
advanced using an update method such as the Euler method:

$\vec{\theta}(\tau+\delta\tau)=\vec{\theta}(\tau)+\delta\tau\times\vec{\dot{\theta}}(\tau)$
(28)

With this information we can then derive a form of the real time evolution
equations in terms of parameterised unitaries. Given an ansatz with $n$
parameters $\vec{\theta}$, $\mathbf{A}^{R}$ is a $(n\times n)$ matrix and
$\mathbf{C}^{I}$ is a vector size $n$ defined as:

$\begin{split}A^{R}_{ij}=\Re\left(\frac{\partial\bra{\psi}}{\partial\theta_{i}}\frac{\partial\ket{\psi}}{\partial\theta_{j}}\right)&=\Re\left(\sum_{a,b}-\frac{i}{2}\cdot-\frac{i}{2}\bra{\overline{0}}U_{1}^{\dagger}...U_{i-1}^{\dagger}\sigma_{a,i}^{\dagger}U_{i}^{\dagger}...U^{\dagger}_{N}U_{N}...\sigma_{b,j}U_{j}U_{j-1}...U_{1}\ket{\overline{0}}\right)\\\
&=\Re\left(\sum_{k,l}-\frac{1}{4}\bra{\overline{0}}\tilde{V}_{k,i}^{\dagger}\tilde{V}_{l,j}\ket{\overline{0}}\right)\end{split}$
(29)

$\begin{split}C^{I}_{i}=\Im\left(-\sum_{k}c_{k}\frac{\partial\bra{\phi}}{\partial\theta_{i}}\hat{P}_{k}\ket{\phi}\right)&=\Im\left(\sum_{k,a}-\frac{i}{2}c_{k}\bra{\overline{0}}U_{1}^{\dagger}...U_{i-1}^{\dagger}\sigma_{a,i}^{\dagger}U_{i}^{\dagger}...U^{\dagger}_{N}\hat{P}_{k}U_{N}...U_{j}U_{j-1}...U_{1}\ket{\overline{0}}\right)\\\
&=\Im\left(\sum_{k,a}f^{*}_{a,i}c_{k}\bra{\overline{0}}\tilde{V}_{a,i}^{\dagger}\hat{P}_{k}V\ket{\overline{0}}\right)\end{split}$
(30)

Where $\mathbf{C}^{I}$ is the imaginary part of the gradient vector with
respect to wavefuction parameters and $\mathbf{A}^{R}$ is the real part of the
metric tensor, which is related to the geometry of the parameter space. With
$H=\sum_{k}c_{k}\hat{P}_{k}$ and $f_{a,i}$ as a phase factor. In order to
calculate the $\mathbf{A}^{R}$ and $\mathbf{C}^{I}$ objects, we show the
quantum circuits to calculate their elements $A^{R}_{i,j}$ and $C^{I}_{i}$ in
Figures 9 and 10. The Hadamard test was chosen as the circuit primitive of
choice but one could equally use another gradient protocol such as the phase
shift method.

$\ket{0}$ & H 1 1 H
$\ket{\bar{0}}$ U_1 U_i-1 σ_a,i U_i U_j-1 σ_b,j

Figure 9: Quantum circuit for the measurement of the $A^{R}_{ij}$ term in the
JCM calculation.

$\ket{0}$ & H S 1 1 H
$\ket{\bar{0}}$ U_1 U_i-1 σ_a,i U_i U_j-1 ^P_k

Figure 10: Quantum circuit for the measurement of the $C^{I}_{i}$ term in the
JCM calculation.

###  C.2 Ansatz for numerical VQS simulations

The VQS method requires an ansatz to express the model’s wavefunction. Chosing
an appropriate ansatz is non-trivial as it must have sufficient freedom
through the variational parameters to access the system Hilbert space, while
keeping a reasonable circuit depth. For the simulations presented in this
work, hardware-efficient style ansatze where sufficient, with initialising $X$
gates so that the system begins with the oscillators in excited states and the
field in the vacuum state at $t=0$. Figures 11 and 12 show the quantum
circuits corresponding to the ansatze used in this work for the JCM and TCM,
respectively.

$|0_{1}\rangle$ & R_x(θ) R_y(θ) R_z(θ) 1 R_x(θ) R_y(θ) R_z(θ) 1
$|0_{0}\rangle$ X R_x(θ) R_y(θ) R_z(θ) R_x(θ) R_y(θ) R_z(θ)

Figure 11: Hardware-efficient ansatz used for the JCM in this work. General
ansatz form with two layers. Two CNOT gates are needed to perform real time
evolution with ansatz initialised to act as $U(\vec{\theta})=\mathbb{I}$. An
$X$ gate is applied on the atom qubit to start the time evolution from the
excited state.

$|0_{2}\rangle$ & R_x(θ) R_y(θ) R_z(θ) 1 2 R_x(θ) R_y(θ) R_z(θ) 1 2
$|0_{1}\rangle$ X R_x(θ) R_y(θ) R_z(θ) 1 R_x(θ) R_y(θ) R_z(θ) 1
$|0_{0}\rangle$ X R_x(θ) R_y(θ) R_z(θ) R_x(θ) R_y(θ) R_z(θ)

Figure 12: Hardware-efficient ansatz used for TCM with $N=2$ in this work.
General ansatz form with two layers. The ordering of the CNOTs in the
entangling layers is specific to this setup (first qubit assigned to the
field, second and third qubits assigned to the atom). The system is
initialised with $X$ gates on the two atom qubits in order to introduce two
quanta in the system.


# Thin film notch filters as platforms for biological image processing

Shaban B. Sulejman<EMAIL_ADDRESS>ARC Centre of Excellence for
Transformative Meta-Optical Systems, School of Physics, The University of
Melbourne, Victoria 3010, Australia.  Niken Priscilla
<EMAIL_ADDRESS>ARC Centre of Excellence for Transformative
Meta-Optical Systems, School of Physics, The University of Melbourne, Victoria
3010, Australia.  Lukas Wesemann  ARC Centre of Excellence for Transformative
Meta-Optical Systems, School of Physics, The University of Melbourne, Victoria
3010, Australia.  Wendy S. L. Lee  ARC Centre of Excellence for Transformative
Meta-Optical Systems, School of Physics, The University of Melbourne, Victoria
3010, Australia.  ARC Centre of Excellence for Transformative Meta-Optical
Systems, Department of Electrical and Electronic Engineering, The University
of Melbourne, Victoria 3010, Australia.  Jieqiong Lou  School of Physics, The
University of Melbourne, Victoria 3010, Australia.  Elizabeth Hinde  School of
Physics, The University of Melbourne, Victoria 3010, Australia.  Timothy J.
Davis  ARC Centre of Excellence for Transformative Meta-Optical Systems,
School of Physics, The University of Melbourne, Victoria 3010, Australia.  Ann
Roberts  ARC Centre of Excellence for Transformative Meta-Optical Systems,
School of Physics, The University of Melbourne, Victoria 3010, Australia.

###### Abstract

Many image processing operations involve the modification of the spatial
frequency content of images. Here we demonstrate object-plane spatial
frequency filtering utilizing the angular sensitivity of a commercial spectral
bandstop filter. This approach to all-optical image processing is shown to
generate real-time pseudo-3D images of transparent biological and other
samples, such as human cervical cancer cells. This work demonstrates the
potential of non-local, non-interferometric approaches to image processing for
uses in label-free biological cell imaging and dynamical monitoring.

##  I Introduction

Transparent objects, including most biological cells, interact weakly with
light resulting in little contrast in conventional bright field microscopy.
However, spatial variations in their morphology and optical properties
introduce local spatial phase variations onto light transmitted through them.
In the simplest case, this can be characterized by a transmission function
$O(x,y)\approx O_{0}e^{i\varphi(x,y)}$. While the amplitude $O_{0}$ is
approximately spatially invariant, the shape and refractive index information
is contained in the phase function $\varphi(x,y)$. Such phase variations
cannot be directly sensed by conventional cameras and requires indirect
detection. Popular optical phase visualization methods include Zernike Zernike
(1942) and differential interference contrast microscopy Lang (1982). However,
these can require expensive components or Fourier plane access that increases
system complexity and size. Digital methods include ptychography Balaur  _et
al._ (2021); Horstmeyer  _et al._ (2016); Giewekemeyer  _et al._ (2010), the
use of the transport of intensity equation Engay  _et al._ (2021); Bostan  _et
al._ (2016); Li  _et al._ (2016) or phase retrieval algorithms such as the
Gerchberg-Saxton Gerchberg and Saxton (1972) and Fienup algorithms Fienup
(1978). However, these are limited by their extensive computational
requirements.

All-optical, object-plane image processing offers a non-interferometric
compact alternative for phase visualization. It is enabled by 2D space-
invariant linear optical systems, such as thin films Wesemann  _et al._
(2019); Zhu  _et al._ (2017), with angular responsivities that directly filter
the spatial frequency of wavefields Wesemann, Davis, and Roberts (2021).
Unlike common computational or all-optical methods utilizing the classical
$4f$-configuration Goodman (1996), it avoids optical phase information losses,
energy consuming post processing and bulky configurations associated with
accessing Fourier planes. The importance of compact optical systems for all-
optical, object-plane image processing is motivated by the potential for
integration into portable devices. This can have applications as diverse as
mobile diagnostics, environmental monitoring and remote sensing.

To explain how a device exhibiting angular dispersion can perform image
processing, we ignore any polarization effects for simplicity. In this case,
the impact of object-plane Fourier filtering on the spatial frequency spectrum
of the field can be described by an optical transfer function
$\mathcal{H}(k_{x},k_{y})$ Davis  _et al._ (2019). By taking the $z$-axis as
the optical axis, $k_{x}$ and $k_{y}$ denote the transverse spatial frequency
components of the wave-vector $\vec{k}=(k_{x},k_{y},k_{z})$ and
$k_{z}=\sqrt{|\vec{k}|^{2}-k_{x}^{2}-k_{y}^{2}}$. The transfer function
relates the processed output to the input by the convolution theorem,

$E_{\text{out}}(x,y,z)=\mathscr{F}^{-1}\left\\{\mathcal{H}(k_{x},k_{y})\tilde{E}_{\text{in}}(k_{x},k_{y};z)\right\\}(x,y)\
,$ (1)

where $\mathscr{F}$ denotes the Fourier transform, $E$ represents any
component of the electric field and
$\tilde{E}_{\text{in}}\equiv\mathscr{F}\left\\{E_{\text{in}}\right\\}$. For
example, high-pass filters block low spatial frequencies to eliminate
unscattered field components for edge detection Roberts, Gomez, and Davis
(2018), which is fundamental to data compression Hui  _et al._ (2019) and
machine vision Brosnan and Sun (2004); Park, Chae, and Cho (2021). A notable
sub-class is that of linear optical transfer functions, i.e.
$\mathcal{H}\propto k_{x}$ or $\mathcal{H}\propto k_{y}$. These can directly
translate phase gradients into intensity variations through spatial
differentiation in the case of transparent objects to permit phase
visualization Wesemann  _et al._ (2019). The influence of polarization can be
incorporated into this approach by utilizing a $2\times 2$ transfer function
dyadic tensor.

The angular sensitivity provides the mechanism for image processing through
the correspondence between angles of incidence and spatial frequencies. This
is given by representing the spatial frequency components in spherical
coordinates Wesemann, Davis, and Roberts (2021),

$\begin{gathered}k_{x}=k_{0}\sin{\theta}\cos{\phi}\\\
k_{y}=k_{0}\sin{\theta}\sin{\phi}\\\ k_{z}=k_{0}\cos{\theta}\ ,\end{gathered}$
(2)

where $k_{0}=|\vec{k}|$ is the wavenumber, while $(\theta,\phi)$ are the polar
and azimuthal propagation angles of plane waves with respect to the $z$-axis.
Given that optical transfer functions represent plane wave responses in
$k$-space, and that light can be decomposed into weighted plane waves by the
spatial Fourier transform Wesemann, Davis, and Roberts (2021); Goodman (2000),
then devices exhibiting angular dispersive transmission are capable of object-
plane image processing.

Recently, meta-optical devices have attracted considerable attention as ultra-
compact image processors Silva  _et al._ (2014). For example, Zhou et al. Zhou
_et al._ (2020) employed photonic crystals for edge detection of organic
samples, while Wesemann et al. Wesemann  _et al._ (2021) obtained phase
contrast images of human cancer cells using a resonant waveguide grating.
Other approaches have involved Mie Wan  _et al._ (2020); Komar  _et al._
(2021) or Fano Cordaro  _et al._ (2019) resonances, photonic spin-orbit
coupling effects Zhu  _et al._ (2019); Wesemann  _et al._ (2022) and bound
states in the continuum Pan  _et al._ (2021). With this rapidly growing
interest in meta-optics, it is timely to consider other optical elements
capable of performing an equivalent role. Earlier works have investigated
various structures for edge detection such as volume hologram filters Case
(1979), Fabry-Pérot etalons Indebetouw (1980), detuned interference filters
Molesini and Quercioli (1981); Molesini, Bertani, and Cetica (1982); Cetica,
Bertani, and Molesini (1982), acousto-optic modulators Indebetouw and Poon
(1984) and gratings Angell (1985); Peri and Ritter (1985); Márquez  _et al._
(2003). With the exception of Fourier plane phase contrast methods, the
emergence of successful digital methods for image processing held back further
progress in all-optical techniques.

Here we demonstrate the use of a commercially available thin film spectral
notch filter applied to phase contrast imaging of weakly absorbing phase
objects. Notch filters are band-stop filters commonly employed in various
types of spectroscopy to remove temporal frequencies over a specific range
Marques  _et al._ (2012); Antonacci, Elsayad, and Polli (2022). We show that
the angular dispersion of the filter’s rejection band produces a high-pass
spatial frequency filter at the operating wavelength. Introducing a phase bias
by tilting the filter with respect to the optical axis can produce pseudo-
three dimensional images similar to those obtained in differential
interference contrast microscopy. Moreover, the contrast generated is
determined by the rotation axis and angle. We demonstrate enhanced contrast
imaging of wavefields introduced by a spatial light modulator and unstained
biological samples, including human cervical cancer cells. Our results confirm
real-time phase contrast imaging without post-processing, permitting direct
imaging with either a camera or the eye. This method provides a novel all-
optical approach for biological and other image processing using off-the-shelf
band-stop devices. It has the potential for developments in machine vision,
biological imaging and dynamical monitoring.

##  II Experimental Methods & Results

###  II.1 Device Performance

The device investigated here is a commercially obtained notch filter (Thorlabs
NF633-25), which has a specified central operating wavelength of
$633\text{\,}\mathrm{nm}$ and band-width of $25\text{\,}\mathrm{nm}$ at normal
incidence. Transmission spectra as a function of angle of incidence were
experimentally measured using the configuration in Fig. 1, with details
provided in the supplementary information (§S2A). The results (Fig. 2(a))
obtained using circularly polarized light are consistent with the
manufacturer’s specifications. Meanwhile, the results for $p$\- and
$s$-polarized light are given in the supplementary material (§2A). Noting the
correspondence between $(\theta,\phi)$ and $(k_{x},k_{y})$, it is found that
the notch filter suppresses low spatial frequencies associated with near-
normal incidence angles at the band-stop wavelength.

Figure 1: The experimental configuration used to capture transmission spectra
of the notch filter. Here, LP and QWP denote linear polarizer and quarter
wave-plate, respectively.

Modulation transfer functions $|\mathcal{H}(k_{x},k_{y})|$ were mapped from
the measured transmission spectra for $p$-, $s$\- and circular polarizations.
Line profiles along $k_{y}=0$ (Fig. 2(b)) at the band-stop wavelength exhibit
approximately polarization-insensitive, high-pass behaviour with a suppression
zone of numerical aperture (NA) range $\sim 0.2$. This is produced by the
blue-shifting of the band-stop wavelength with increasing angle of incidence.

Figure 2: The experimental transmission spectra obtained by incrementally
rotating the filter is given in (a). The modulation transfer function along
$k_{y}=0$ produces the plot shown in (b) for various polarizations, which is
compared to data provided by the manufacturer. Linear fitting of the
modulation transfer function in the contrast zone is given in (c).

A significant feature of the device is the region of approximate linear
dependence on $k_{x}$ from zero to near-unity transmission beyond the
suppression zone. This is supported by linear fitting (Fig. 2(c)) from
$k_{x}/k_{0}\sim 0.2$ to $0.3$, represented by the relation $Ak_{x}/k_{0}+B$.
The polyfit function and curve fitting toolbox in MATLAB were used to obtain
the fitting parameters and their error intervals as $A=5.9\pm 0.83$ and
$B=-0.97\pm 0.22$. Rotating the filter by
$12\text{\,}\mathrm{\SIUnitSymbolDegree}$ about a line perpendicular to the
optical axis shifts operation to $k_{x}/k_{0}\approx\pm 0.2$ to access this
region, referred to here as the contrast zone. Such a linear optical transfer
function substituted into Eq. (1) approximately produces the spatial
derivative around an intensity offset along the axis of rotation. This can
generate intensity images directly related to phase gradients Wesemann, Davis,
and Roberts (2021), given by

$\nabla O(x,y)\approx iO_{0}\left(\nabla\varphi(x,y)\right)e^{i\phi(x,y)}\ .$
(3)

Consequently, intensity contrast is created in regions where the phase is
varying along the relevant direction by the direct relationship in Eq. (3).
Operating near the edge of the contrast zone removes unscattered field
components leaving only relatively large phase gradients that enhance edges.
However, operating within the contrast zone at a rotation angle of
$14\text{\,}\mathrm{\SIUnitSymbolDegree}$, corresponding to
$k_{x}/k_{0}\approx 0.24$, preserves some unscattered components and enables
discrimination between positive and negative phase gradients. These manifest
as different grayscale levels above or below the shifted $k$-space origin,
respectively. Therefore, operation within the contrast zone enables phase
visualization by converting phase gradients to intensity. This includes
information about the sign of the gradient to produce pseudo-3D images. Phase
gradients along different directions can also be visualized by changing the
axis of rotation. These form the key underpinnings of this paper to describe
the capacity of notch filters for phase visualization.

###  II.2 Phase Contrast Imaging

All-optical image processing was experimentally performed on various samples
with the notch filter using circularly polarized, $635\text{\,}\mathrm{nm}$
laser light. Edge detection of amplitude objects was first demonstrated on a
USAF resolution test target, with details and results provided in the
supplementary information (§S2B). Phase imaging was then performed using the
configuration (Fig. 3) detailed in the supplementary information (§S2C). A
computer-controlled spatial light modulator emulated phase profiles of human
red blood cells (Fig. 4(a)), which were modelled using the properties obtained
from Ref. 41. Collimated light reflected from the spatial light modulator
passed through the notch filter, which was rotated within the focal plane
between paired microscope objectives. A camera then captured the transmitted
images.

Figure 3: The configuration for image processing with a spatial light
modulator (SLM), where L and MO denote lenses and microscope objectives,
respectively.

The results are presented in Fig. 4(b)-(e), where simulated and experimental
images are compared. An experimental control image without the filter (Fig.
4(b)) displays poor contrast, rendering the object virtually invisible as
expected for weakly absorbing phase objects. Simulated and experimental phase
images (Fig. 4(c) and (d)) are obtained by introducing the filter at a
rotation angle of $14\text{\,}\mathrm{\SIUnitSymbolDegree}$. These preserved
some unscattered contributions leading to the appearance of pseudo-3D phase
contrast. Moreover, the simulated phase image (Fig. 4(c)) is consistent with
experiment (Fig. 4(d)). In both cases, regions where a change in phase is
present, i.e. $\nabla\varphi(x,y)\neq 0$, are transformed to intensity
variations arising from otherwise invisible phase modulations. Moreover,
ringing artifacts associated with the Gibbs phenomenon can also be seen. The
resultant images additionally possess the capacity to distinguish between
positive and negative phase gradients. Furthermore, line profiles (Fig. 4(e))
exhibit intensity variations associated with phase variation introduced into
the field by the sample. Some intensity artifacts can also be seen, such as
those arising from amplitude variations in the field.

Figure 4:  A red blood cell (a) was emulated by a spatial light modulator with
a phase excursion of $6\pi/5$. Simulated and experimental images normalized to
their brightest pixels are given in (b)-(d), comprising control (b), simulated
(c) and experimental (d) phase contrast images. These are supported by line
profiles in (e) along the dashed lines shown in (b)-(d).

To illustrate potential microscopy applications, phase contrast microscopy
images were obtained of biological samples with weak amplitude contrast. This
was performed using an inverted microscope configuration depicted in Fig.
5(a). Moreover, human cervical cancer (HeLa) cells were used as the sample and
the preparation steps are outlined in the supplementary information (§S2E).
Their bright field image (Fig. 5(b)) displays little to no contrast. However,
placing the notch filter immediately beneath the sample enabled visualization
of the phase variations when illuminated with $635\text{\,}\mathrm{nm}$ laser
light. The phase contrast image (Fig. 5(c)) obtained within the contrast zone
contains morphological detail absent in the corresponding bright field image,
such as those highlighted within the green curves. The contrast produced
exhibits significantly enhanced phase visualization compared to the relatively
low contrast bright field image. Cellular thickness and local refractive index
deviations are accented to discriminate the cells from their background. A
conventional differential interference contrast image was also taken of the
same regions, which is given in the supplementary information (§S2D) as a
baseline for comparison.

Figure 5:  An experimental schematic for biological phase imaging is given in
(a) and outlined in the supplementary information (§S2D). A bright field image
of the HeLa cells obtained without the filter is given in (b). The
corresponding phase contrast image obtained using the filter is given in (c).

##  III Discussion

Inspired by the recent surge in interest towards meta-optical imaging, these
results represent the first demonstration of biological image processing using
a commercially available filter. The spatial frequency content of images was
directly altered by the notch filter in the object plane without needing to
access the Fourier plane. The transfer function exhibited the required
behaviour for high-pass filtering enabling edge detection, in addition to
approximately linear regions allowing phase visualization through spatial
differentiation. Rotating the filter enabled access to these regimes to offset
the Fourier origin in $k$-space. Phase gradients were converted into intensity
variations to generate images with significantly enhanced contrast. As a
result, notch filters offer a new approach to non-interferometric phase
visualization through image processing without bulky or expensive components.

Commercial notch filters offer a relatively cost efficient and readily
available option for phase visualization compared to alternative optics.
Although the band-stop region was limited within the visible band, various
filters have been designed with coverage across the electromagnetic spectrum.
For example, Yuan et al. Yuan, Lou, and Dong (2020) presented a photonic
filter employable at the absorption band of acetylene gas. Others include
terahertz Bark  _et al._ (2021), infrared Kwon and Shin (2004), microwave Long
_et al._ (2017); Gertler  _et al._ (2022) and commercial ultra-violet filters.
These hold potential for phase imaging across the spectrum using the methods
outlined in this paper, with applications in machine vision and biological
imaging beyond the visible range. Live biological monitoring is possible with
the notch filter dynamically capturing phase information of specimens. Further
applications include integration into conventional CMOS sensors to form
spatial frequency and phase-sensitive detectors. Incorporating notch filters
into their detection planes would introduce spatial frequency selectivity and
phase variations in otherwise indiscriminate, intensity-based devices.

Although phase contrast was produced by the notch filter, its behaviour in the
contrast zone was not perfectly linear. Furthermore, its restricted numerical
aperture limits the range of samples to which it could be applied. The large
angular range of the suppression zone also restricts contrast to sufficiently
sharp features. Not only is the range of objects limited, but it necessitates
rotation of the notch filter. Finally, our theoretical and numerical modelling
neglected the non-negligible notch filter thickness, which can introduce
perturbations into the beam path. This can lead to the appearance of
aberrations such as those in the edge enhanced images in the supplementary
information (§S2B). However, the success of experiments performed here provide
confidence that custom thin film devices could be developed with reduced
thicknesses to minimize aberrations and tailor the transfer function.

##  IV Conclusion

In conclusion, this paper has demonstrated real-time, all-optical, object-
plane image processing using a commercial spectral notch filter. Spectroscopic
measurements verified angle dispersive band-rejection necessary for high-pass
spatial frequency filtering and phase contrast imaging. Edge detection was
realizable through the suppression zone where unscattered field components
were removed. Meanwhile, offsetting to approximately linear regions within the
contrast zone produced phase contrast images. Unstained biological samples,
including human cervical cancer cells, imposing otherwise invisible phase
modulations were visualized by the filter. Given its real-time capabilities,
this has significant implications in label-free biomedical imaging Park,
Depeursinge, and Popescu (2018) including medical diagnostics, non-invasive
micro-organism growth and dynamical monitoring. The results open developmental
possibilities in extending to beyond the visible band and constructing
monolithic spatial frequency-sensitive cameras for commercialization.

## Funding

This research was funded by the Australian Government through the Australian
Research Council Centre of Excellence grant (CE200100010). S.B.S also
gratefully acknowledges the support of the Ernst & Grace Matthaei Scholarship
and the Australian Government Research Training Program Scholarship.

## Supplementary Material

Supporting information is available in the supplementary material, including
the underlying theory and detailed experimental methods.

## References

  * Zernike (1942) F. Zernike, “Phase contrast, a new method for the microscopic observation of transparent objects,” Physica 9, 686–698 (1942).
  * Lang (1982) W. Lang,  _Nomarski Differential Interference-Contrast Microscopy_ (Carl Zeiss, 1982).
  * Balaur  _et al._ (2021) E. Balaur, G. A. Cadenazzi, N. Anthony, A. Spurling, E. Hanssen, J. Orian, K. A. Nugent, B. S. Parker, and B. Abbey, “Plasmon-induced enhancement of ptychographic phase microscopy via sub-surface nanoaperture arrays,” Nature Photonics 15, 222–229 (2021).
  * Horstmeyer  _et al._ (2016) R. Horstmeyer, J. Chung, X. Ou, G. Zheng, and C. Yang, “Diffraction tomography with fourier ptychography,” Optica 3, 827–835 (2016).
  * Giewekemeyer  _et al._ (2010) K. Giewekemeyer, P. Thibault, S. Kalbfleisch, A. Beerlink, C. M. Kewish, M. Dierolf, F. Pfeiffer, and T. Salditt, “Quantitative biological imaging by ptychographic x-ray diffraction microscopy,” Proceedings of the National Academy of Sciences 107, 529–534 (2010).
  * Engay  _et al._ (2021) E. Engay, D. Huo, R. Malureanu, A.-I. Bunea, and A. Lavrinenko, “Polarization-dependent all-dielectric metasurface for single-shot quantitative phase imaging,” Nano Letters 21, 3820–3826 (2021).
  * Bostan  _et al._ (2016) E. Bostan, E. Froustey, M. Nilchian, D. Sage, and M. Unser, “Variational phase imaging using the transport-of-intensity equation,” IEEE Transactions on Image Processing 25, 807–817 (2016).
  * Li  _et al._ (2016) J. Li, Q. Chen, J. Sun, J. Zhang, and C. Zuo, “Multimodal computational microscopy based on transport of intensity equation,” Journal of Biomedical Optics 21, 126003 (2016).
  * Gerchberg and Saxton (1972) R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35, 237–246 (1972).
  * Fienup (1978) J. R. Fienup, “Reconstruction of an object from the modulus of its fourier transform,” Optics Letters 3, 27–29 (1978).
  * Wesemann  _et al._ (2019) L. Wesemann, E. Panchenko, K. Singh, E. D. Gaspera, D. E. Gómez, T. J. Davis, and A. Roberts, “Selective near-perfect absorbing mirror as a spatial frequency filter for optical image processing,” APL Photonics 4, 100801 (2019).
  * Zhu  _et al._ (2017) T. Zhu, Y. Zhou, Y. Lou, H. Ye, M. Qiu, Z. Ruan, and S. Fan, “Plasmonic computing of spatial differentiation,” Nature Communications 8, 15391 (2017).
  * Wesemann, Davis, and Roberts (2021) L. Wesemann, T. J. Davis, and A. Roberts, “Meta-optical and thin film devices for all-optical information processing,” Applied Physics Reviews 8, 031309 (2021).
  * Goodman (1996) J. W. Goodman,  _Introduction to Fourier Optics_ (McGraw-Hill, 1996).
  * Davis  _et al._ (2019) T. J. Davis, F. Eftekhari, D. E. Gómez, and A. Roberts, “Metasurfaces with asymmetric optical transfer functions for optical signal processing,” Physical Review Letters 123, 013901 (2019).
  * Roberts, Gomez, and Davis (2018) A. Roberts, D. E. Gomez, and T. J. Davis, “Optical image processing with metasurface dark modes,” Journal of the Optical Society of America A 35, 1575–1584 (2018).
  * Hui  _et al._ (2019) G. Hui, H. Ruyong, W. Chaopeng, L. Zequn, W. Le, and Z. Shengmei, “Compressed ghost edge imaging,” Chinese Optics Letters 17, 071101 (2019).
  * Brosnan and Sun (2004) T. Brosnan and D.-W. Sun, “Improving quality inspection of food products by computer vision––a review,” Journal of Food Engineering 61, 3–16 (2004), applications of computer vision in the food industry.
  * Park, Chae, and Cho (2021) K. Park, M. Chae, and J. H. Cho, “Image pre-processing method of machine learning for edge detection with image signal processor enhancement,” Micromachines 12 (2021), 10.3390/mi12010073.
  * Goodman (2000) J. W. Goodman,  _Statistical Optics_ (John Wiley & Sons, 2000).
  * Silva  _et al._ (2014) A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343, 160–163 (2014).
  * Zhou  _et al._ (2020) Y. Zhou, H. Zheng, I. I. Kravchenko, and J. Valentine, “Flat optics for image differentiation,” Nature Photonics 14, 316–323 (2020).
  * Wesemann  _et al._ (2021) L. Wesemann, J. Rickett, J. Song, J. Lou, E. Hinde, T. J. Davis, and A. Roberts, “Nanophotonics enhanced coverslip for phase imaging in biology,” Light: Science & Applications 10, 98 (2021).
  * Wan  _et al._ (2020) L. Wan, D. Pan, S. Yang, W. Zhang, A. A. Potapov, X. Wu, W. Liu, T. Feng, and Z. Li, “Optical analog computing of spatial differentiation and edge detection with dielectric metasurfaces,” Optics Letters 45, 2070–2073 (2020).
  * Komar  _et al._ (2021) A. Komar, R. A. Aoni, L. Xu, M. Rahmani, A. E. Miroshnichenko, and D. N. Neshev, “Edge detection with mie-resonant dielectric metasurfaces,” ACS Photonics 8, 864–871 (2021).
  * Cordaro  _et al._ (2019) A. Cordaro, H. Kwon, D. Sounas, A. F. Koenderink, A. Alù, and A. Polman, “High-index dielectric metasurfaces performing mathematical operations,” Nano Letters 19, 8418–8423 (2019).
  * Zhu  _et al._ (2019) T. Zhu, Y. Lou, Y. Zhou, J. Zhang, J. Huang, Y. Li, H. Luo, S. Wen, S. Zhu, Q. Gong, M. Qiu, and Z. Ruan, “Generalized spatial differentiation from the spin hall effect of light and its application in image processing of edge detection,” Phys. Rev. Applied 11, 034043 (2019).
  * Wesemann  _et al._ (2022) L. Wesemann, J. Rickett, T. J. Davis, and A. Roberts, “Real-time phase imaging with an asymmetric transfer function metasurface,” ACS Photonics 9, 1803–1807 (2022).
  * Pan  _et al._ (2021) D. Pan, L. Wan, M. Ouyang, W. Zhang, A. A. Potapov, W. Liu, Z. Liang, T. Feng, and Z. Li, “Laplace metasurfaces for optical analog computing based on quasi-bound states in the continuum,” Photonics Research 9, 1758–1766 (2021).
  * Case (1979) S. K. Case, “Fourier processing in the object plane,” Opt. Lett. 4, 286–288 (1979).
  * Indebetouw (1980) G. Indebetouw, “Tunable spatial filtering with a fabry-perot etalon,” Applied Optics 19, 761–764 (1980).
  * Molesini and Quercioli (1981) G. Molesini and F. Quercioli, “Interference filters as fourier processors,” Applied Optics 20, 1893–1893 (1981).
  * Molesini, Bertani, and Cetica (1982) G. Molesini, D. Bertani, and M. Cetica, “Dark ground microscopy with detuned interference filters,” Optical Engineering 21, 1061 – 1063 (1982).
  * Cetica, Bertani, and Molesini (1982) M. Cetica, D. Bertani, and G. Molesini, “Interference filters for fourier filtering and for frequency coding,” Journal of Non-Crystalline Solids 47, 279–281 (1982).
  * Indebetouw and Poon (1984) G. Indebetouw and T.-C. Poon, “Incoherent spatial filtering with a scanning heterodyne system,” Appl. Opt. 23, 4571–4574 (1984).
  * Angell (1985) D. K. Angell, “Incoherent spatial filtering with grating interferometers,” Appl. Opt. 24, 2903–2906 (1985).
  * Peri and Ritter (1985) D. Peri and D. Ritter, “Spatial filtering with volume gratings,” Appl. Opt. 24, 1535–1540 (1985).
  * Márquez  _et al._ (2003) A. Márquez, C. Neipp, A. Beléndez, S. Gallego, M. O. no, and I. Pascual, “Edge-enhanced imaging with polyvinyl alcohol /acrylamide photopolymer gratings,” Opt. Lett. 28, 1510–1512 (2003).
  * Marques  _et al._ (2012) C. A. F. Marques, V. de Oliveira, H. J. Kalinowski, and R. N. Nogueira, “Production of optical notch filters with fine parameter control using regenerated fiber bragg gratings,” Opt. Lett. 37, 1697–1699 (2012).
  * Antonacci, Elsayad, and Polli (2022) G. Antonacci, K. Elsayad, and D. Polli, “On-chip notch filter on a silicon nitride ring resonator for brillouin spectroscopy,” ACS Photonics 9, 772–777 (2022).
  * Evans and Fung (1972) E. Evans and Y.-C. Fung, “Improved measurements of the erythrocyte geometry,” Microvascular Research 4, 335–347 (1972).
  * Yuan, Lou, and Dong (2020) Z. Yuan, X. Lou, and Y. Dong, “Microwave photonic notch filter with a tunable frequency and a bandwidth based on gas absorption,” Optics Letters 45, 248–251 (2020).
  * Bark  _et al._ (2021) H. S. Bark, K.-H. Jang, K. Lee, Y. U. Jeong, and T.-I. Jeon, “Thz guided-mode resonance notch filter with variable filtering strength,” Scientific Reports 11, 1307 (2021).
  * Kwon and Shin (2004) M.-S. Kwon and S.-Y. Shin, “Tunable notch filter using a thermooptic long-period grating,” J. Lightwave Technol. 22, 1968 (2004).
  * Long  _et al._ (2017) Y. Long, J. Xia, Y. Zhang, J. Dong, and J. Wang, “Photonic crystal nanocavity assisted rejection ratio tunable notch microwave photonic filter,” Scientific Reports 7, 40223 (2017).
  * Gertler  _et al._ (2022) S. Gertler, N. T. Otterstrom, M. Gehl, A. L. Starbuck, C. M. Dallo, A. T. Pomerene, D. C. Trotter, A. L. Lentine, and P. T. Rakich, “Narrowband microwave-photonic notch filters using brillouin-based signal transduction in silicon,” Nature Communications 13, 1947 (2022).
  * Park, Depeursinge, and Popescu (2018) Y. Park, C. Depeursinge, and G. Popescu, “Quantitative phase imaging in biomedicine,” Nature Photonics 12, 578–589 (2018).

Supplementary Material

##  S1 Theoretical Modelling

###  S1.1 Phase Visualization via Spatial Frequency Filtering

The core idea underlying this work is to directly modify spatial frequencies
of wavefields without accessing Fourier planes to enable phase extraction. The
essential feature of suitable devices is a sensitivity to angles of incidence.
As explained in the main article, this associates with a selective
transmission to spatial frequencies necessary for image processing.

Image processing for phase visualization can be explained using classical
scalar electromagnetic theory Roberts, Gomez, and Davis (2018). Polarization
effects can be accounted through tensorial analysis Wesemann, Davis, and
Roberts (2021), but are ignored here for simplicity. Suppose that the source
electric field can be modelled as monochromatic plane waves $E_{s}(x,y,z)$ of
wavelength $\lambda$ and wavevector $\vec{k}=(k_{x},k_{y},k_{z})$. With the
$z$-axis as the optical axis, $k_{x}$ and $k_{y}$ denote the transverse
spatial frequency components and $k_{0}\equiv|\vec{k}|=\frac{2\pi}{\lambda}$
is the wavenumber. The source illuminates a weakly absorbing object located in
the $z=0$ plane, which is assumed as optically thin in the Kirchhoff
approximation. Its transmitted light can be modelled by

$E_{\text{obj}}(x,y,0)=O(x,y)E_{s}(x,y,0)\ ,$ (4)

where the object transmission function is approximately given by

$O(x,y)\approx O_{0}e^{i\varphi(x,y)}\ .$ (5)

Its intensity image $|O(x,y)|^{2}=|O_{0}|^{2}$ has little to no contrast by
construction and is therefore referred to as a phase object. Therefore, its
visualization requires access to the phase function $\varphi(x,y)$ that
contains its shape information. This is made possible through optical systems,
such as thin films Zhu  _et al._ (2017); Wesemann  _et al._ (2019), described
by specific optical transfer functions that can extract $\varphi(x,y)$ from
the complex exponential of Eq. (5). To shed light on this, we note the
following property for the Fourier transform of derivatives,

$\mathscr{F}\left\\{\frac{\partial^{m}E}{\partial
x^{m}}+\frac{\partial^{m}E}{\partial
y^{m}}\right\\}(k_{x},k_{y})=i^{m}(k_{x}^{m}+k_{y}^{m})\mathscr{F}\\{E(x,y)\\}(k_{x},k_{y})\
,$ (6)

where $\mathscr{F}$ denotes the Fourier transform and $m\in\mathbb{N}$. As in
Eq. (1) of the main article, the object light filtered by the notch filter is
related to its optical transfer function $\mathcal{H}(k_{x},k_{y})$,

$E_{\text{out}}(x,y,0)=\mathscr{F}^{-1}\left\\{\mathcal{H}(k_{x},k_{y})\mathscr{F}\left\\{E_{\text{obj}}(x^{\prime},y^{\prime},0)\right\\}(k_{x},k_{y})\right\\}(x,y)\
,$ (7)

where the notch filter is taken to immediately succeed the object. By
comparing Eqs. (6) and (7), the $m^{th}$-order derivatives of the input field
are produced if the optical transfer function is given by

$\mathcal{H}(k_{x},k_{y})=\begin{cases}i^{m}k_{x}^{m},\ &\text{for
x-differentiation}\\\ i^{m}k_{y}^{m},\ &\text{for
y-differentiation}\end{cases}\ ,$ (8)

up to a multiplicative constant. It was shown in the main article that the
transfer function exhibited approximately linear variation within the contrast
zone. The notch filter is therefore capable of approximate first-order spatial
differentiation in this region required for phase contrast imaging.
Substituting a linear optical transfer function into Eq. (7) and taking
$E_{s}(x,y,0)=\text{const.}$ produces Eq. (3) of the main article, given by

$\nabla O(x,y)\approx iO_{0}\left(\nabla\varphi(x,y)\right)e^{i\phi(x,y)}\ .$
(9)

Phase gradients of the object are translated into measurable intensity
variations described by the amplitude term in Eq. (9). Thus, its intensity
image is now given by

$|\nabla O(x,y)|^{2}\approx|\nabla\varphi(x,y)|^{2}\ ,$ (10)

with contrast generated where $\nabla\varphi\neq 0$. This holds for
approximately pure phase objects, with the gradient taken along the axis of
rotation. However, non-negligible amplitude variations in $O(x,y)$ would
produce additional amplitude derivatives in Eq. (9). Although, it is possible
to subtract these by using measured intensities in the absence of the filter
to isolate the phase contribution.

###  S1.2 Simulations

Numerical simulations were performed in Python by implementing the theoretical
model presented in §S2.1. Our experimentally measured modulation transfer
function and the phase response provided the manufacturer were combined to
model the optical response of the notch filter.

##  S2 Experimental Methods

###  S2.1 Transfer Function Determination

Spectroscopic measurements of the notch filter (Thorlabs NF633-25) were made
using the experimental configuration depicted in Fig. 1 of the main article.
The filter was $3.5\text{\,}\mathrm{mm}$ thick and comprised tantalum
pentoxide (Ta2O5) and silicon dioxide (SiO2) thin films on a fused quartz
substrate. White light from a fibre-coupled (Thorlabs SM600 single-mode fibre)
halogen lamp (Ocean Insight HL-2000-FHSA) was collimated by a microscope
objective (Nikon UPlanFL 20x/0.15NA). This light was polarized by a linear
polarizer (Thorlabs LPVIS050-MP) to produce $p$ or $s$-polarized light.
Circularly polarized light was produced by introducing an additional quarter
waveplate (Thorlabs AQWP05M-600) before illuminating the notch filter. A
plano-convex lens (Thorlabs LA1068-A, $f=$75\text{\,}\mathrm{mm}$$) focused
the transmitted light onto a port that was fibre-coupled (Thorlabs M15L01) to
a spectrometer (Ocean Insight QE6500). Measurements of transmission spectra
were made for incident angles ranging from
$-30\text{\,}\mathrm{\SIUnitSymbolDegree}$ to
$30\text{\,}\mathrm{\SIUnitSymbolDegree}$ by incrementally rotating the filter
by $2\text{\,}\mathrm{\SIUnitSymbolDegree}$.

Figure 6: Transmission spectra of the notch filter obtained with $p$\- and
$s$-polarized light are given in (a) and (b), respectively. These exhibit
similar behaviour with suppression around $633\text{\,}\mathrm{nm}$ (dashed
lines) that blue-shifts with increasing incident angle.

Spectra were collected for $p$-, $s$\- and circularly polarized light. The
former two are given in Fig. 6, while the latter is given in Fig. 2(a) of the
main article. Each of these confirm suppression about the band-stop wavelength
and its angular dispersion.

###  S2.2 Edge Detection

All-optical amplitude image processing was experimentally performed to
demonstrate edge detection with the notch filter. As depicted in Fig. 7,
fibre-coupled (Thorlabs SM600) laser light (Thorlabs S1FC635
$635\text{\,}\mathrm{nm}$) was collimated by a microscope objective (Nikon LU
Plan 5x/0.15NA) and two lenses (Thorlabs LB1901-B $f=$75\text{\,}\mathrm{mm}$$
& Thorlabs LB1761-B $f=$25.4\text{\,}\mathrm{mm}$$). Light reflected from a
mirror was projected onto the the largest regions of a USAF resolution test
target (Thorlabs R2L2S1N1 NBS 1963A) using a second microscope objective
(Olympus UPlanFI 20x/0.50NA) and lens (Thorlabs LA1509-A
$f=$100\text{\,}\mathrm{mm}$$). Upon reflection from a second mirror, a fourth
lens (Thorlabs LA1433-A $f=$150\text{\,}\mathrm{mm}$$) allowed image filtering
with the notch filter between paired microscope objectives (Nikon UPlanFI
20x/0.4NA & Nikon LU Plan 50x/0.55NA). Finally, a fifth lens (Thorlabs
LA1131-A $f=$50\text{\,}\mathrm{mm}$$) projected the filtered light onto a
camera (Thorlabs DCC1645C).

Figure 7: A schematic of the experimental configuration for edge detection
with a notch filter is shown. Here L, MO and M each denote lenses, microscope
objectives and mirrors, respectively.

Given in Figs. 8(a) and (b), the simulated and experimental results
demonstrate edge enhancement produced by the notch filter. The uniform
intensity regions associated with unscattered field components were removed
while retaining only the edges. This is supported by line profiles given in
Fig. 8(c). Some slight asymmetries are present in the experimental edge-
enhanced image owing to small rotations of the filter.

Figure 8: Edge-enhanced images of USAF test target bars overlayed on their
unfiltered images are given in (a) and (b). Line profiles in (c) depict
suppression within the uniform regions of the bars.

###  S2.3 Phase Imaging

All-optical image processing for phase visualization was performed using the
configuration given in Fig. 3 of the main article. Fibre-coupled (Thorlabs
SM600) laser light (Thorlabs S1FC635, $635\text{\,}\mathrm{nm}$) was
collimated by a microscope objective (Nikon LU Plan 5x/0.15NA) and two plano-
convex lenses (Thorlabs LA1027-A $f=$75\text{\,}\mathrm{mm}$$ and Thorlabs
LA1509-A $f=$100\text{\,}\mathrm{mm}$$). Required for its operation, a linear
polarizer (Thorlabs LPVIS050-MP) ensured that linearly polarized light
illuminated a reflective spatial light modulator (Holoeye Pluto-2-VIS-014
LCOS-SLM). Comprising 1920 $\times$ 1080 liquid crystal pixels of
$8\text{\,}\mathrm{\SIUnitSymbolMicro m}$ pitch, it prescribed computer-
generated phase profiles of human red blood cells. The reflected light from
the SLM passed through a quarter wave-plate (Thorlabs AQWP05M-600) to
circularly polarize it. A third plano-convex lens (Thorlabs LA1433-A
$f=$150\text{\,}\mathrm{mm}$$) and microscope objective (Olympus PlanN
20x/0.4NA) together formed a telescope that de-magnified the SLM image onto
the notch filter in its focal plane. To access the contrast zone, the latter
was longitudinally rotated by $14\text{\,}\mathrm{\SIUnitSymbolDegree}$ using
a rotation mount. A second, identical microscope objective (Olympus PlanN
20x/0.4NA) and a fourth plano-convex lens (Thorlabs LA1131-A
$f=$50\text{\,}\mathrm{mm}$$) re-magnified and focused the filtered image onto
a camera (Thorlabs DCC1645C).

###  S2.4 Phase Microscopy

Bright field and phase contrast microscopy of the HeLa cells were performed
using the configuration given in Fig. 5(a) of the main article. First, bright
field images were obtained using white light from an inverted microscope
(Nikon Ti-80i) without the notch filter. The source was then replaced by
fibre-coupled (Thorlabs SM600) laser light (Thorlabs S1FC635
$635\text{\,}\mathrm{nm}$) collimated by an objective (Olympus A4 4x/0.1NA),
and circularly polarized by a linear polarizer (Thorlabs LPVISC100-MP2) and
quarter waveplate (Thorlabs AQWP05M-600). These were contained within an
optical cage mounted on an XYZ-stage to vary incident angles. Moreover, the
notch filter was place immediately beneath the sample to perform phase
imaging. A second microscope objective (Nikon LU Plan 50x/LWD) collected the
filtered light onto a camera (Andor Zyla sCMOS 4.2P). Given in Fig. 9(b), a
differential interference contrast microscopy image was taken on a second
microscope (Olympus BX60) using a separate microscope objective (Olympus PlanN
20x/0.4NA).

Figure 9: A differential interference contrast microscopy image of the HeLa
cells is given in (b). For comparison purposes, its bright field image (a) and
experimental phase contrast image obtained with the notch filter (c) are
reproduced from the main article.

###  S2.5 Biological Sample Preparation

The HeLa cells were grown in Dulbecco’s modified Eagle’s medium (DMEM) (Lonza)
supplemented with 10% heat inactivated bovine growth serum (Gibco), 1x Pen-
Strep (Lonza) at $37\text{\,}\mathrm{\SIUnitSymbolCelsius}$ with 5% carbon
dioxide. The cells were plated $24\text{\,}\mathrm{h}$ before fixation onto
$35\text{\,}\mathrm{mm}$ glass bottom dishes, before being fixed with 4%
paraformaldehyde for $15\text{\,}\mathrm{min}$ at room temperature and washed
3 times with phosphate-buffered saline (PBS).

## References

  * Roberts, Gomez, and Davis (2018) A. Roberts, D. E. Gomez, and T. J. Davis, “Optical image processing with metasurface dark modes,” Journal of the Optical Society of America A 35, 1575–1584 (2018).
  * Wesemann, Davis, and Roberts (2021) L. Wesemann, T. J. Davis, and A. Roberts, “Meta-optical and thin film devices for all-optical information processing,” Applied Physics Reviews 8, 031309 (2021).
  * Zhu  _et al._ (2017) T. Zhu, Y. Zhou, Y. Lou, H. Ye, M. Qiu, Z. Ruan, and S. Fan, “Plasmonic computing of spatial differentiation,” Nature Communications 8, 15391 (2017).
  * Wesemann  _et al._ (2019) L. Wesemann, E. Panchenko, K. Singh, E. D. Gaspera, D. E. Gómez, T. J. Davis, and A. Roberts, “Selective near-perfect absorbing mirror as a spatial frequency filter for optical image processing,” APL Photonics 4, 100801 (2019).