# Classical benchmarking of zero noise extrapolation beyond the exactly-
verifiable regime

Sajant Anand  Dept. of Physics, University of California, Berkeley, CA 94720,
USA  Kristan Temme  IBM Quantum, IBM Thomas J. Watson Research Center,
Yorktown Heights, NY, USA  Abhinav Kandala  IBM Quantum, IBM Thomas J. Watson
Research Center, Yorktown Heights, NY, USA  Michael Zaletel  Dept. of Physics,
University of California, Berkeley, CA 94720, USA  Material Science Division,
Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA

###### Abstract

In a recent work a quantum error mitigation protocol was applied to the
expectation values obtained from circuits on the IBM Eagle quantum processor
with up $127$ \- qubits with up to $60\;-\;\mbox{CNOT}$ layers. To benchmark
the efficacy of this quantum protocol a physically motivated quantum circuit
family was considered that allowed access to exact solutions in different
regimes. The family interpolated between Clifford circuits and was
additionally evaluated at low depth where exact validation is practical. It
was observed that for highly entangling parameter regimes the circuits are
beyond the validation of matrix product state and isometric tensor network
state approximation methods. Here we compare the experimental results to
matrix product operator simulations of the Heisenberg evolution, find they
provide a closer approximation than these pure-state methods by exploiting the
closeness to Clifford circuits and limited operator growth. Recently other
approximation methods have been used to simulate the full circuit up to its
largest extent. We observe a discrepancy of up to $20\%$ among the different
classical approaches so far, an uncertainty comparable to the bootstrapped
error bars of the experiment. Based on the different approximation schemes we
propose modifications to the original circuit family that challenge the
particular classical methods discussed here.

Quantum error mitigation (QEM) has been proposed as a method for extending the
reach of near-term quantum hardware before the implementation of quantum error
correction [1, 2]. While quantum error correction is widely expected to be
necessary for implementing generic quantum algorithms, the overhead is
currently prohibitive. QEM instead compensates for the effect of errors by
engineering the noise in a manner which allows postprocessing to increase the
accuracy of expectation values.

Figure 1: Comparison of classical approximations for $\langle Z_{62}\rangle$
at Trotter depth $D=20$ against experimental ZNE results: (1) matrix-product-
state (MPS) representation of the pure state within a lightcone-reduced volume
[3]; (2) extrapolation of the MPS results with respect to the estimated
circuit fidelity (see Supp.); (3) Belief propagation tensor network states
[4]; (4) MPO representation of Heisenberg evolution (this work); (5)
simulation of a 31 qubit subset of the IBM Eagle device [5]; (6) Clifford
perturbation theory [6]. The latter four methods differ by $\sim 20\%$ amongst
themselves near $\theta_{h}\sim\pi/4$, an amount largely within the spread of
the ZNE error bars. Without further calculations it is not clear which of
these methods is most accurate.

In a recent work [3], we reported results benchmarking the efficacy of an
experimentally feasible QEM variant, “zero noise extrapolation” (ZNE) [1, 2,
7] for estimating observables using error-prone quantum hardware. To do so, we
considered a quantum circuit describing the discretized dynamics of the 2D
transverse field Ising model (TFI) (i.e., the “kicked” transverse Ising model
[8, 9]). Each round of the circuit takes the form

$\displaystyle U(\theta_{J},\theta_{h})$
$\displaystyle=\prod_{\braket{i,j}}R_{ZZ}^{ij}(\theta_{J})\prod_{i}R_{X}^{i}(\theta_{h})$
$\displaystyle=\prod_{\braket{i,j}}e^{-i\theta_{J}Z_{i}Z_{j}/2}\prod_{i}e^{-i\theta_{h}{X_{i}}/2}$
(1)

Here $\langle i,j\rangle$ runs over the bonds of IBM’s “Eagle” 127-qubit
heavy-hexagon lattice. The two-qubit angle $\theta_{J}=-\pi/2$ was chosen to
be Clifford so as to involve only a single CNOT gate, and we assume this value
unless otherwise specified. At $\theta_{h}=0,\pi/2$ the overall circuit
dynamics are Clifford and therefore exactly calculable, while at intermediate
$\theta_{h}$ they are expected to be ergodic.

In the reported experimental protocol, the system is first initialized in
state $\otimes^{127}_{i=1}\ket{\uparrow}$ and evolved under $D$-rounds of $U$,
at which point a variety of weight-1, 10 and 17 observables were estimated.
Results were compared with exact classical calculations where available, as
well as 1D and 2D tensor network simulations which approximate the evolution
of the pure quantum state. These classical approximations are expected to
break down at large depths due to the growth of quantum entanglement, which is
particularly rapid at $\theta_{h}=\pi/2$.

Experimentally, ZNE was found to provide a large improvement in the accuracy
of expectation values compared with unmitigated results. We found:

  1. 1.

In the regime where exact classical simulations can verify results (including
circuits with a depth of 15 CNOT layers, in which the measured operator spread
up to 68 sites), ZNE reproduces the exact results within the bootstrap error
bars (Fig. 1)

  2. 2.

Classical matrix-product state (MPS) and 2D isometric tensor-network (isoTNS)
methods which approximate the dynamics of the pure state, while accurate near
$\theta_{h}\sim 0$, struggle to accurately reproduce expectation values for
the full range of $\theta_{h}$ values.

  3. 3.

In regimes beyond exact verification, ZNE produces results which are correct
at both of the exactly solvable points $\theta_{h}=0,\pi/2$, while the
classical pure-state approaches fail badly near $\theta_{h}\sim\pi/2$.

The observation that ZNE could provide such reliable expectation values, at
this scale, was argued to be evidence for the utility of pre - fault tolerant
quantum computers. The argument for utility is that one can use noisy quantum
processors as a tool to reliably explore circuits and problems that will
ultimately pose difficult challenges for classical simulation methods. It is
therefore crucial to emphasize that the Ref. [3] did not claim “quantum
advantage” or any provable speedup, and it was left open whether other
approximate classical methods might perform better than the pure state
methods. Indeed, it was suggested that approximating the Heisenberg evolution
of the measured operators, rather than the quantum states themselves, was a
promising future direction.

In this work we report the results of such additional classical calculations.
In contrast to the pure-state methods, we find that in the verifiable regime,
matrix-product approximations of the Heisenberg evolution give excellent
agreement with the exact results, and therefor with ZNE. Going beyond the
verifiable regime, we find classical simulations remain in good agreement with
ZNE even up to 20 Trotter steps (60 CNOT layers), and reproduce the exact
results at $\theta_{h}=\pi/2$. These results provide further evidence that ZNE
produces accurate expectation values at scales orders of magnitude beyond both
exact verification and the unmitigated results.

At depth $D=20$, there remains a range of $\theta_{h}\sim\pi/4$ where
Heisenberg simulations are not fully converged, and they do not precisely
agree with results of the newly developed “BP-TNS” classical approximation
recently reported in [4]; Clifford perturbation theory (CPT) reported in [10,
6]; and simulations on a smaller 31-qubit geometry reported in [5]. The
various methods mutually disagree at about the 20% level, which happens to be
comparable to the ZNE error bars which sit between them - see Fig.1. The
Heisenberg evolution is a fully controlled approximation, meaning it is in
principle exact, but only in the limit of a bond dimension which scales
exponentially in $D^{2}$. BP-TNS, on the other hand, is fully converged at a
more favorable bond dimension of $2^{D}$, but makes uncontrolled
approximations which make it inexact even in this limit. CPT is a perturbative
expansion in $\tan(\theta_{h})$ (truncated at order $K=10$ in Ref. [6]) which
is no longer a small parameter at the $\theta_{h}\sim\pi/4$ point. Smaller-
size methods depend on the restricted growth of the operator, and ultimately
will require a careful finite-size scaling analysis. Without additional study,
it is not yet clear which of the methods is most accurate in this regime. We
anticipate further improvements in classical methods or quantum hardware will
prove fruitful here.

Figure 2: Evolution of operators using matrix product operators (MPO). (a) The
1D MPO representation of an operator can be viewed as a vectorized state in a
larger Hilbert space of onsite dimension $k=4$. (b) An operator expectation
value is found by evolving the vectorized operator by $U^{\dagger}\otimes U$,
in this work represented by the product of 13 MPOs, each of bond dimension 4.

##  I Heisenberg-MPO Method

Figure 3: Comparison of Heisenberg evolution against experimental ZNE results
for (a) weight-10, (b) weight-17, and (c) weight-17 with modified circuit, all
at depth 5. For the first two, exact answer is available either by brute force
or lightcone-reduced MPS simulations, and Heisenberg evolution accurately
matches this. For (c), where the exact answer is not available, Heisenberg and
ZNE agree.

We begin by describing the numerical method used to approximate the Heisenberg
evolution of an observable,
$\mathcal{O}(D)=(U^{\dagger})^{D}\,\mathcal{O}\,U^{D}$. We take the standard
approach of approximating $\mathcal{O}(D)$ via a 1D matrix-product operator
representation (MPO) [11, 12]. The MPO can be viewed as a vectorized MPS in a
doubled Hilbert space of local bond dimension 4, as shown in Fig. 2(a). To map
the 2D heavy-hex lattice to a 1D chain, sites are ordered according to the
“snake” ordering of Ref. [3]. This particular ordering minimizes the number of
long-range connections in the resulting MPS. As in the standard method for
time-evolving 1D operators, we interleave unitary conjugation with variational
matrix-product compression [13].

However, mapping the 2D lattice to the 1D chain introduces long-range
couplings which require some care to implement efficiently. If a full round of
$U$ is applied, the MPO dimension would increase by a factor of 4096, making
subsequent truncation intractable111We note that the bra and ket legs of the
operator must be evolved together, without truncation in between. If not,
evolution of a Pauli string under a non-trivial Clifford circuit is not
guaranteed to yield a single Pauli string.. Instead, we decompose $U$ into 13
layers, $U=\prod_{r=1}^{13}U_{r}$, chosen such that the exact application
$O\to U_{r}^{\dagger}OU_{r}$ increases $\chi\to 4\chi$ on the evolved bonds.
As the two-qubit gates commute, the layers $U_{r}$ are chosen to evenly
distribute the gates among the layers. To prevent blowup in $\chi$, two-site
variational matrix product compression [13] is then applied between each
layer. Simulations are conducted at fixed matrix-product bond dimension
$\chi$, leading to errors at long times. After $D$ steps, we then compute
$\langle\mathcal{O}(D)\rangle=\braket{\psi}{\mathcal{O}(D)}{\psi}$ for
$\ket{\psi}=\otimes_{i=1}^{N}\ket{\uparrow}$, $N=127$. Note that the evolved
operator can be measured with any initial state, pure or otherwise. This is
shown in Fig. 2(b) in the vectorized picture. Each Trotter step takes time
$N\chi^{3}$, for total complexity $DN\chi^{3}$.

##  II Benchmarking the Kicked Ising Circuits

We now consider the circuits of Ref. [3]. In Fig. 3, we show results for the
weight-10 operator
$X_{\\{3\\}}Y_{\\{2\\}}Z_{\\{5\\}}=U(\pi/2)^{5}Z_{13}\left[U^{\dagger}(\pi/2)\right]$
(a) and the weight-17 operator
$X_{\\{8\\}}Y_{\\{1\\}}Z_{\\{8\\}}=U(\pi/2)^{5}Z_{58}\left[U^{\dagger}(\pi/2)\right]^{5}$
(b), both measured at circuit depth 5. Both of these operators trivially have
expectation value 1 (0) at the Clifford $\theta_{h}=\pi/2$ ($0$) point. We
find that the results from Heisenberg evolution of the measured operator agree
to good precision with the exact results, and thus with the ZNE-experimental
results, across all ranges of $\theta_{h}$. In Appendix A, we quantitatively
compare to the exact results and demonstrate agreement better than $10^{-4}$
(often significantly) for all $\theta_{h}$. In fact, for the weight-10
operator in Fig. 3(a), a bond dimension of $\chi=384$ is sufficient for exact
simulation. Note that unlike for pure state methods, which were considered
previously in Ref. [3], the bond dimension needed at the non-trivial Clifford
point $\theta_{h}=\pi/2$ is 1, despite the circuit generating a volume law
entangled stabilizer state. This is because a Pauli operator remains a single
Pauli string, albeit expanded in extent, throughout the evolution. Moving away
from the either Clifford point $\theta_{h}=0,\,\pi/2$ has a perturbative
effect on the Heisenberg evolution, as the evolved operator becomes a
superposition of Pauli strings concentrated in the vicinity of the Pauli
string resulting from Clifford evolution. Thus near either Clifford point,
Heisenberg evolution performs well and needs minimal bond dimension to capture
the dynamics.

Figure 4: Operator spreading in the heavy-hex Kicked Ising dynamics. Growth in
operator entanglement entropy (a) consistent with $O_{EE}(D)\propto D^{2}$ is
consistent with an operator spreading over a disk. The plateau is an artifact
of finite-$\chi$, which causes the circuit fidelity (b) to decay. Data is for
the Kicked TFI with $\theta_{J}=-\pi/2$, a non-Clifford version
$\theta_{J}=-\pi/4$, and modified $\theta_{J}=-\pi/2$ Kicked TFI model with
additional layers of $R_{X}(\theta_{h})$ gates between each of the three
layers of $R_{ZZ}$ gates that constitute a Trotter step so that the two-qubit
gate layers become non-commuting. All models have $\theta_{h}=0.7$. (c)
Spatial profile of the out-of-time order correlator $C(7,x)=\langle
Z_{x}Z_{62}(7)Z_{x}Z_{62}(7)\rangle$ for the experimental model, showing slow
spread of the OTOC with circuit depth. The 54 sites in the lightcone at depth
$D=7$ are shown by the crosses.

Having built confidence in Heisenberg MPO evolution to match the exact results
where available, we move to circuits where an exact result is not available.
To clarify, by exact we mean either a full contraction of the light cone
reduced circuit, or MPS dynamics with an analytically known bond dimension
that entails no truncation. In Fig. 3(c), we report results of the weight-17
operator
$X_{\\{8\\}}Y_{\\{8\\}}Z_{\\{1\\}}=\prod_{i}e^{-i\theta_{h}{X_{i}}/2}U(\pi/2)^{5}Z_{58}\left[U^{\dagger}(\pi/2)\right]^{5}\prod_{i}e^{i\theta_{h}{X_{i}}/2}$
measured at circuit depth 5 with an additional layer of
$\prod_{i}e^{-i\theta_{h}{X_{i}}/2}$ gates at the end. The expectation value
of this operator is entirely equivalent to the expectation value of
$U(\pi/2)^{6}Z_{58}\left[U^{\dagger}(\pi/2)\right]^{6}$ measured at circuit
depth 6. We find good agreement between experiment and Heisenberg results and
additionally find the known answers at the two Clifford angles. Given the
convergence in bond dimension and reasonable circuit fidelities shown in
Appendix A, this demonstrates the applicability of ZNE beyond the exactly
verifiable regime.

We next turn to single-site magnetization $\langle Z_{62}\rangle$ at depth 20,
already shown in Fig. 1. Before discussing the results, let us describe the 6
different numerical methods: (1) Heisenberg MPO evolution described in Sec. I,
(2) pure state MPS evolution shown in [3], (3) extrapolation of pure state MPS
results with respect to the estimated fidelity $F_{D}\to 1$ (see Appendix A),
(4) the recently developed heavy-hexagon belief propagation TNS (BP-TNS)
method [4], (5) Clifford Perturbation Theory [10, 6], and (6) simulations of a
31 qubit subset of the device [5]. Methods (1), (2), and (4) are run at a
fixed bond dimension $\chi$, and confidence in the results stems from apparent
convergence in bond dimension, truncation error, and for (1) and (2), the
circuit fidelity. Method (5) was run at order $K=10$.

We find that all numerical methods agree broadly with the experimental ZNE
results, yet differences are present across the range of single qubit angles
$\theta_{h}$. We find very good agreement between all methods and experiment
for $\theta_{h}\leq\pi/8$. For $\theta_{h}>\pi/8$, the experimental error bars
fall in between the results from Heisenberg MPO evolution and BP-TNS
evolution. As the MPO evolution is exact as $\chi\to\infty$, a natural
improvement would be to extrapolate the MPO results with respect to the
estimated circuit fidelity, as we did for pure state MPS. However, we often
find the results are not monotonic in $\chi$, as shown in Appendix A, making
such extrapolation difficult. Further improvements in both the classical
methods and experiment will be necessary to resolve the remaining discrepancy.

Why do MPO calculations prove sufficient at this scale, despite the naive
scaling $\chi\sim k^{D^{2}}$? To investigate we characterize the operator
growth in the Kicked Ising circuit, which governs the difficulty of MPO
simulations. We first consider the operator entanglement entropy (OEE) of the
Heisenberg-evolved $Z_{62}(D)$. The OEE is defined by interpreting an operator
as a wavefunction in a doubled Hilbert space, [15] and computing the resulting
entanglement entropy of a bipartition; we focus on the largest OEE across all
bipartitions of the 1D ordering. Focusing on $\theta_{h}=0.7$, we observe a
quadratic growth with circuit depth, $\sqrt{O_{EE}}\sim 0.11D$ (the plateau at
large times is a numerical artifact of finite-$\chi$ truncation; see Fig.
4(a)). This scaling is consistent with an operator growing over a disk of area
$\propto D^{2}$, but with a prefactor far below the maximal amount allowed by
the lightcone (e.g., $O_{EE}\leq 0.7D^{2}$ at $D=7$). A priori this may either
because the “butterfly velocity” $v_{B}$, which governs the rate of operator
spread, is less than than that of the lightcone, or because the nearby
Clifford point delays the growth in operator entanglement.

In Fig. 4(c), we plot the out of time ordered correlator (OTOC) defined as
$C(D,x)=\langle Z_{x}Z_{62}(D)Z_{x}Z_{62}(D)\rangle$, which is straightforward
to evaluate once $Z_{62}(D)$ is obtained in MPO form222We simply evaluate the
expectation value of the doubled operator $Z\otimes Z$ using the vectorized
MPO as an MPS.. At $D=7$, we see that the OTOC signal is confined well within
the lightcone, indicating $v_{B}$ quite a bit below 1. This was also recently
observed in Ref. [5], where it was exploited to efficiently simulate the
Kicked Ising circuits with reduced qubit count.

##  III Comparison of operator methods with BP-TNS

Belief propagation tensor network states (BP-TNS) [17, 4] on the heavy-hexagon
lattice were recently used to simulate the Kicked Ising experiments and were
found to match the exact results where available and provide good agreement
with the ZNE results more generally [4]. The BP-TNS method shows great promise
for approximating the local observables of moderate-depth dynamics, so we take
a moment to consider some strengths and weaknesses relative to the operator
based approaches.

BP-TNS is two-dimensional (2D) tensor network (PEPs) ansatz in which the
“environment” of each region (an input for the calculation of local
expectations values) is approximated by assuming it takes the form of an
uncorrelated tensor-product across ancilla bonds. Under this approximation -
which is uncontrolled - the environment can be self-consistently calculated in
a procedure analogous to statistical “belief propagation” [17]. In contrast to
1D MPS / MPO methods, BP-TNS need not reproduce the exact results even as the
bond dimension $\chi\to\infty$, a point we will return to.

Figure 5: Expectation value of the depth-$D$ Clifford stabilizer
$\mathcal{S}_{D}=U(\pi/2)^{D}Z_{i}\left[U^{\dagger}(\pi/2)\right]^{D}$ with
respect to the depth-$D$ circuit $U(\pi/2)^{D}$. BP-TNS calculations are
conducted at $\chi=2^{7}=128$. We see that once $\chi<2^{D}$ there is a
catastrophic loss in accuracy, as expected for any pure-state TNS method.
Heisenberg simulations remain highly accurate due to the proximity to the
Clifford point, though still suffer from an exponential blowup more generally.

Like other 2D TNS methods, BP-TNS naturally adapts to the 2D geometry of the
heavy hexagon lattice, avoiding the need of MPS/MPO algorithms to turn short-
range 2D interactions into long-range 1D interactions. Additionally,
approximating the environments as uncorrelated reduces the computational cost
compared to more expressive 2D tensor network algorithms (e.g. PEPs / isoTNS)
[18, 19], allowing BP-TNS to scale to much larger bond dimensions ($\chi\sim
200$) than is possible with standard 2D TNS ($\chi\leq 12$). This is crucial
for the kicked TFI circuits near the $\theta_{h}=\pi/2$ Clifford point, where
the required bond dimension for an accurate 2D TNS representation scales as
$2^{D}$. In particular, at $\theta_{h}=\pi/2$ the entanglement spectrum on any
bond of the network, given by the diagonal $\Lambda$ matrices of BP-TNS, has
rank $2^{D}$ and is exactly flat. Thus any truncation of the bond dimension
will immediately lead to noticeable errors in expectation values.

The effect of BP-TNS bond truncation can be seen in Figure 5, where we use BP-
TNS to measure a high-weight stabilizer $\mathcal{S}_{D}\equiv
U(\pi/2)^{D}Z_{i}\left[U^{\dagger}(\pi/2)\right]^{D}$ of the
$\theta_{h}=\pi/2$ circuit as a function of depth $D$. In order to compare
with brute-force exact results we consider the “two-hexagon” geometry with 21
qubits, shown in Figure 5(a), which retains the salient features of the 127
qubit geometry. Using the “extended time-evolution” method proposed in Ref.
[4], we measure $\mathcal{S}_{D}$ by first evolving a circuit forward $D$
steps with single-qubit angle $\theta_{h}$ and then backward $D$ steps with
angle $\pi/2$:
$\ket{\psi}=\left[U^{\dagger}(\pi/2)\right]^{D}U(\theta_{h})^{D}\ket{\uparrow}^{N}$.
This allows for high-weight $\mathcal{S}_{D}$ to be measured as a single-site
observable:
$\braket{\uparrow}{\left[U^{\dagger}(\theta_{h})\right]^{D}\mathcal{S}_{D}U(\theta_{h})^{D}}{\uparrow}=\braket{\psi}{Z_{i}}{\psi}$,
where $i$ is shown as a red qubit in Figure 5(a). We consider up to depth
$D=10$ with fixed bond dimension $\chi=128$ for both BP-TNS and Heisenberg
evolution. In Figure 5(b), we demonstrate that as one approaches the Clifford
point $\theta_{h}=\pi/2$, using a BP-TNS bond dimension less than $2^{D}$
leads to a catastrophic reduction in the expectation value. The difficulties
of BP-TNS for large $\theta_{h}$ occur for the same reason MPS and isoTNS
pure-state simulations struggled in the original simulations of Ref. [3]: the
evolved state is highly entangled. Heisenberg evolution, on the other hand,
reproduce the dynamics of the $\pi/2$ Clifford point, regardless of circuit
depth $D$, with bond dimension $\chi=1$ (or CPT order $K=0$) as the evolved
operator is a single Pauli string. Of course away from the Clifford points
these methods also suffer from an eventual exponential blowup, but the
difficulty is pushed out to higher depths.

Figure 6: Kicked-Ising ($\theta_{J}=-\pi/4,\,\theta_{h}=\pi/2$) realization of
a many-body double slit experiment. An initial $x$-polarized state is prepared
with a spin flip on site 0 of a double-hexagon geometry (Fig.5a), which then
propagates under kicked-Ising dynamics as revealed by the site-resolved
magnetization. In the exact result (bottom row), the magnetization density
exhibits constructive / deconstructive interference on site 6, $D\sim 10$
depending on whether $\phi=0/\pi$ flux threads the left hexagon. Within the
BP-TNS approximation (top row), there is no sign of interference: the result
is entirely independent of the flux through the hexagon. This is because the
BP-TNS approximation cannot capture the loop-like entanglement required by the
double-slit experiment.

_The effect of loop-like entanglement._ A more interesting aspect of BP-TNS
arises from its approximate treatment of the TNS “environment.” Even if
$\chi=2^{D}$, the BP-TNS environment approximation becomes inexact once the
lightcone of any site encompasses a loop within the lattice. As pointed out in
Ref. [4], BP-TNS is thus well suited to the heavy-hex lattice, whose shortest
loop has length 12. For the kicked TFI circuit $U(\theta_{h})$, BP-TNS is then
exact up to depth $D=6$ if bond dimension $2^{6}=64$ is used, which happens to
be depth of the circuits amenable to brute-force calculation where
verification was conducted. Beyond $D=6$, the BP-TNS makes uncontrolled
approximations, and it is an interesting question whether the Kicked Ising
circuit ever features the “loop-like” correlations which evades the BP-TNS
approximation. At $D=20$, it is at least quantitatively suggestive that near
the $\theta_{h}\sim\pi/4$ point, there is a $20\%$ discrepancy in $\langle
Z_{62}\rangle$ between BP-TNS and the CPT / Heisenberg methods _even once
extrapolating_ the BP-TNS $\chi\to\infty$. This is even more clear near
$\theta_{h}=3\pi/8$, where ZNE, CPT, and Heisenberg all give $\langle
Z_{62}\rangle\sim 0$, while BP-TNS gives $\langle Z_{62}\rangle\sim 0.05$. CPT
and Heisenberg methods are expected to be quite accurate here.

One fanciful way to illustrate the generation of loop-like entanglement in
kicked-Ising dynamics is to conduct a many-body version of the two-slit
experiment. To do so we consider a modified circuit in which the $ZZ$-gates
are made non-Clifford:
$V=U(\theta_{J}=-\pi/4,\,\theta_{h}=\pi/2)=R_{ZZ}(-\pi/4)R_{X}(\pi/2)$. We
then consider a state polarized along $X=1$, flip a spin at site $i$, and
measure the change in $X$-magnetization at site $j$:

$\displaystyle
C_{ij}(D)=\bra{\rightarrow^{N}}Z_{i}\left[V^{\dagger}\right]^{D}X_{j}V^{D}Z_{i}\ket{\rightarrow^{N}}.$
(2)

The dynamics cause the spin flip to propagate from $i\to j$, albeit on top of
a rapidly thermalizing background. Choosing $i/j$ to lie on the green / red
sites of Fig. 5(a), the dynamics of the spin flip effectively realize the two-
slit experiment as the spin flip enters into a superposition of the top and
bottom arms of the hexagon. In Fig.6, we show the dynamics of the
$X$-magnetization in the left 12-site hexagon for both the exact and BP-TNS
approximation ($\chi=128$, re-gauging with 15 “message passing” iterations).
In the exact result, we see a wave front which propagates and constructively
interferes when the top and bottom fronts collide on the “red” site. BP-TNS
reproduce the exact results beautifully before the collision but fail to
capture the interference phenomena as the fronts overlap. This is because a
spin-flip in such a superposition is exactly the sort of entanglement
neglected in the BP-TNS approximation: the influence of the top and bottom
arms on the “red” site is assumed to factorize. This failure can be quantified
by threading $\pi$-flux (e.g. flipping the sign of one $ZZ$-bond) through the
hexagon and repeating the experiment. In the exact simulation, we see the
wavefronts now interfere _destructively_ ; BP-TNS, in contrast, gives
identical results for both fluxes.

We do not expect such a simple picture to apply to deeper circuits; the
double-split example required a coherence length comparable to the length of
the minimal loop. As the loop develops more connections to the surrounding
geometry, or as the effective energy density increases due to the driving, we
expect additional many-body degrees of freedom will obtain “which-path”
information which destroys the interference, and the BP-TNS approximation may
become more accurate.

##  IV Circuit extensions

Having highlighted some strengths and weaknesses of both pure state and
operator evolution methods, we now propose extensions to experiments that may
lead to more challenging simulations, both classically and for ZNE.

_Non-Clifford two qubit gates._ Reducing the two-qubit angle away from the
Clifford point, say $\theta_{J}\to-\pi/4$ will increase the difficulty of the
CPT and MPO-Heisenberg methods as it increases the rate at which new Pauli
strings are generated. In Fig. 4(a,b), we compute the $O_{EE}$ and MPO
fidelity at $\theta_{J}=-\pi/4,\theta_{h}=0.7$. We see $\sqrt{O_{EE}}\sim
0.18D$, almost 2x the rate of the Clifford 2Q case, and correspondingly the
MPO fidelity drops much faster with gate depth. Note, however, that a non-
Clifford two-qubit gate may make pure-state simulations such as BP-TNS easier,
as such a gate is less entangling than one with $\theta_{J}=-\pi/2$.

_Non-commuting two qubit gates._ $R_{ZZ}$ is composed of three layers of
disjoint $ZZ$ gates. While naively each layer could advance the lightcone,
because they commute, the lightcone grows by only one site in each direction
per $R_{ZZ}$. The simplification is obstructed if the two-qubit gates are non-
commuting; or equivalently, by sprinkling the $R_{X}$ gates amongst the three
layers. In Fig. 4(a,b), we compute the $O_{EE}$ and fidelity at
$\theta_{J}=-\pi/2,\theta_{h}=0.7$, performing a round of $R_{X}$ gates after
each layer of $ZZ$ gates, thus tripling the number of single qubit gates. In
addition to increasing the size of the lightcone, we find that the OEE grows
much more rapidly than for the commuting kicked TFI models considered
($\sqrt{O_{EE}}\sim 0.33D$), and the plateau imposed by finite bond dimension
is much higher than for the commuting version also with $\theta_{J}=-\pi/2$.
However, the region with non-unital OTOC is still a fraction of the circuit
lightcone, indicating $v_{B}<1$ but higher than the commuting case.
Experimentally, this modification has the advantage of an identical two-qubit
gate count.

Figure 7: OTOC growth for Kicked TFI variants. The lightcone at depth $D=7$ is
depicted by the crosses. (a) $C(7,x)$ for the
$(\theta_{J},\theta_{h})=(-\pi/4,0.7)$ model. Compared to the version with
$\theta_{J}=-\pi/2$, shown in Fig. 4(c), we see that the OTOC grows both in
size and spatial extent by depth 7. (b) $C(4,x)$ for the
$(\theta_{J},\theta_{h})=(-\pi/2,0.7)$ non-commuting model introduced in the
text. We see the lightcone at depth $D=4$ is larger (69 sites) compared to
that for commuting models (54 sites) at a greater depth and that the OTOC has
spread more than it has in other models.

Of course at fixed $\theta_{h}$, the final signal $\langle Z\rangle$ depends
on these modifications as well, and for ZNE to be performant it must be
measurable. In the non-trotterized dynamics ($\theta_{h},\theta_{J}\to 0$
holding $\theta_{h}/\theta_{J}$ fixed), we expect a continuous phase
transition tuned by $\theta_{h}/\theta_{J}$, corresponding to finite-$T$
equilibrium symmetry breaking of the Ising model. Even at finite Trotter step,
the moderate-time dynamics will be affected by the transition before ultimate
ergodicity, making this an interesting parameter regime for comparison of ZNE
and classical approaches.

_Echo experiments._ From the perspective of benchmarking ZNE, another future
direction is to measure the analog of high-weight operators designed to ensure
a large signal. Away from the Clifford point the high-weight operators are
complex, so can instead be measured according to

$\displaystyle Z_{i}(D|\theta,\theta^{\prime})\equiv$
$\displaystyle\bra{\textrm{CPS}}\left[U^{\dagger}(\theta^{\prime})\right]^{D}U(\theta)^{D}$
$\displaystyle
Z_{i}\left[U^{\dagger}(\theta)\right]^{D}U(\theta^{\prime})^{D}\ket{\textrm{CPS}}$
(3)

Such an experiment is essentially a Loschmidt echo, but restricting the return
probability to a single site rather than the global fidelity.

##  V Conclusion

We have simulated the circuits considered in Ref. [3], which we previously
studied with pure state methods, via matrix product compression Heisenberg
evolution of the operators. We match the exact answer where available and find
general agreement with ZNE experimental results. For the largest circuit
depths considered, we find 10 - 20% disagreements between other recently
reported classical methods [4, 6, 5]. The classical uncertainty appears to be
within the bounds of the ZNE uncertainty, and likely can be further reduced
with additional resources. The Heisenberg approach reveals a relatively slow
growth of operators in the studied parameter regime, and modifications to the
circuits which would make similar calculations more difficult are discussed.

###### Acknowledgements.

SA and MZ were supported by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences, under Early Career Award No. DE-SC0022716. We thank Y.
Kim and A. Eddins for insightful experimental data and simulations. We thank
S. Bravyi for suggestions on additional simulations. We thank J. Tindall, M.
Fishmann, M. Stoudenmire, K. Siva, T. Soejima, and S. Garratt for helpful
discussions. Computing resources were provided by National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility located at Lawrence Berkeley National Laboratory,
operated under Contract No. DE-AC02-05CH11231 using NERSC Award No. BES-
ERCAP0020043.

##  Appendix A Additional numerics

Here we present additional numerics for the simulations discussed in the main
text.

###  A.1 Exactly verifiable operators

We begin in Fig. 8 with the error in simulations of the exactly solvable
weight-10 and weight-17 observables shown in Fig. 3(a,b). We compare results
with Heisenberg evolution to the exact answer and find good agreement across a
range of $\theta_{h}$ values. Note that $\chi=384$ is sufficient to obtain the
exact answer for the weight-10 operator, as no truncation occurs during the
depth-5 simulations. For the weight-17 operator, we find that increasing the
bond dimension does not always reduce the error, as is not uncommon for time
evolution which is not variational. Individual simulations for the weight-10
(17) operator take less than 2 (90) minutes at the maximum bond dimensions
considered.

Figure 8: Error for Heisenberg simulation of the (a) weight-10 and (b)
weight-17 operators where the exact solution is available. We generally find
accurate quantitative agreement with the exact result across a wide range of
$\theta_{h}$, building confidence when moving to more difficult circuits.

###  A.2 MPO Circuit Fidelity

Before moving on to more difficult circuits, let us introduce a metric to
bound the accuracy of matrix-product operator based time evolution. Recall
that at each Trotter step, i.e. each application of the 13 MPOs that apply a
layer of two-qubit gates on every bond on the lattice, the bond dimension
$\chi$ grows and must be truncated via variational compression. Upon
compressing the evolved an enlarged MPO
$\tilde{O}_{t}=U_{r}^{\dagger}O_{t-1}U_{r}\to O_{t}$ back to fixed $\chi$, a
truncation error $\epsilon_{t}=|\tilde{O}_{t}-O_{t}|$ is incurred. Here we use
the Frobenius norm, scaled so that all Pauli strings have unit norm. We can
then upper-bound the total error after $D$ steps by
$\epsilon(D,\chi)\leq\sqrt{\sum_{t}\epsilon_{t}^{2}}$. We record this as a
fidelity, $f_{t}=|\braket{O_{t}}{\tilde{O}_{t}}|^{2}$, so that
$F_{D}\equiv\prod_{t=1}^{D}f_{t}\leq e^{-\sum_{t=1}^{D}\epsilon_{t}^{2}}\leq
e^{-\epsilon^{2}(D,\chi)}$. As $\chi\to\infty$, $F_{D}\to 1$ as each step of
the evolution will be done exactly. We use the fidelity to both extrapolate
results from finite bond dimension and qualitatively gauge the confidence in
our simulation. An equivalent estimate of the final fidelity is also
applicable in the context of MPS pure states, as discussed for example in
Ref.[20].

###  A.3 Modified Weight-17 Operator

Next we turn to the weight-17 operator for the modified circuit: 5 layers of
the circuit defined by Eq. 1, followed by an additional layer of $R_{X}$
rotation gates. The depth-5 stabilizer of this circuit is equivalent to the
depth-6 stabilizer of the original circuit and was chosen to access this
larger operator with the same circuit resources. In Fig. 9(a), we show the
circuit fidelity $F_{D}$ as a function of $\theta_{h}$ for several bond
dimenions. We see that at the two Clifford points $\theta_{h}=0,\,\pi/2$ that
the fidelity is exactly 1, as a bond dimension of $\chi=1$ is sufficient to
track the evolution of a Pauli string. We see that the fidelity increases
demonstrably with increasing chi, which as been argued to indicated that the
exact bond dimension for this circuit isn’t exponentially far removed from
what is achievable in simulations [20]. In Fig. 9(b), we investigate the
convergence of the expectation values at different $\theta_{h}$ as a function
of bond dimension. Note that the expectation value is not changing
monotonically with bond dimension.

Figure 9: (a) Fidelity as a function of single qubit angle $\theta_{h}$ for
different $\chi$ for Heisenberg simulation. Near the Clifford points, the
simulations have fidelity $1$ with small bond dimension, but there is a
significant dip in the intermediate regime. (b) Convergence of
$\langle\mathcal{O}\rangle$ as a function of $\chi$. As the exact answer is
not available, we look at differences from the value at $\chi=1024$, the
largest bond dimension considered. Non-monotonic behavior in $\chi$ is
observed.

###  A.4 $\langle Z_{62}\rangle$ \- MPS Extrapolation

Finally we turn to the depth-20 circuits for $\langle Z_{62}\rangle$. In Fig.
10 we demonstrate how we extrapolate the pure state, Schrödinger picture MPS
simulations from finite bond dimension to $\chi=\infty$. A clear trend in
$\langle Z\rangle$ in $1/\chi$ is not seen in Fig. 10(a), so we instead use
our circuit fidelity. We fit the results for the largest three bond dimensions
as a linear function in $\log(F_{D})$ and extrapolate to $F_{D}\to 1$. As
shown in Fig. 1, extrapolating MPS data to $F_{D}=1$ (achieved at
$\chi=\infty$ where each step of the time evolution can be done exactly as no
truncation is ever done) brings the results closer to the ZNE experiment as
well as the other computational methods, but it is not sufficient to find
agreement.

Figure 10: Extrapolation of $\langle Z_{62}\rangle$ using circuit fidelity.
(a) Expectation value as a function of $1/\chi$ does not display a clear power
law for extrapolation. (b) Extrapolation performed at $\theta_{h}=0.7$ with
circuit fidelity, $\langle Z\rangle=a\cdot\log(F_{D})+b$.

###  A.5 $\langle Z_{62}\rangle$ \- MPO Extrapolation

A natural next step is to try extrapolation of the Heisenberg results in
circuit fidelity. We show the results of this in Fig. 11(a), where the
extrapolation for $\theta_{h}=0.7$ is shown in Fig. 11(c). Note that unlike
for the MPS simulations, the expectation value is often non-monotonic in
$\chi$ (or equivalently $F_{D}$), which can be seen in Fig. 11(b), which makes
extrapolation less justified. To demonstrate, we extrapolate the same way as
with MPS, using the three largest bond dimension results. We see that the
extrapolated results are reasonable (and lie amongst the other classical
methods and ZNE experimental results) for $\theta_{h}\leq\pi/4$, but for
$\pi/4<\theta_{h}\leq 3\pi/8$ are generally inconsistent with other methods.
The expectation value in this region is decreasing with $\chi$, which leads to
the negative extrapolated value. If trends at smaller $\theta_{h}$ are to
continue, the expectation value will begin increasing with $\chi$ after some
point, but the necessary $\chi$ to see this may be exponentially large [20].

Figure 11: (a) Attempted extrapolation of Heisenberg results with circuit
fidelity show behavior inconsistent with other classical methods and ZNE
experimental results at intermediate $\theta_{h}$. (b) Examples of non-
monotonic behavior found at various $\theta_{h}$. (c) Example of extrapolation
performed at $\theta_{h}=0.7$.

## References

  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Phys. Rev. Lett. 119, 180509 (2017).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Phys. Rev. X 7, 021050 (2017).
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. van den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, and A. Kandala, Nature 618, 500 (2023).
  * Tindall  _et al._ [2023] J. Tindall, M. Fishman, M. Stoudenmire, and D. Sels, “Efficient tensor network simulation of ibm’s kicked ising experiment,” (2023), arXiv:2306.14887 [quant-ph] .
  * Kechedzhi  _et al._ [2023] K. Kechedzhi, S. V. Isakov, S. Mandrà, B. Villalonga, X. Mi, S. Boixo, and V. Smelyanskiy, “Effective quantum volume, fidelity and computational cost of noisy quantum processing experiments,” (2023), arXiv:2306.15970 [quant-ph] .
  * Begušić and Chan [2023] T. Begušić and G. K.-L. Chan, “Fast classical simulation of evidence for the utility of quantum computing before fault tolerance,” (2023), arXiv:2306.16372 [quant-ph] .
  * Kandala  _et al._ [2019] A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Nature 567, 491 (2019).
  * Prosen [2000] T. Prosen, Progress of Theoretical Physics Supplement 139, 191 (2000).
  * Pineda  _et al._ [2014] C. Pineda, T. Prosen, and E. Villaseñor, New Journal of Physics 16, 123044 (2014).
  * Begušić  _et al._ [2023] T. Begušić, K. Hejazi, and G. K.-L. Chan, “Simulating quantum circuit expectation values by clifford perturbation theory,” (2023), arXiv:2306.04797 [quant-ph] .
  * Verstraete  _et al._ [2004] F. Verstraete, J. J. Garcia-Ripoll, and J. I. Cirac, Physical review letters 93, 207204 (2004).
  * Zwolak and Vidal [2004] M. Zwolak and G. Vidal, Physical review letters 93, 207205 (2004).
  * Schollwöck [2011] U. Schollwöck, Annals of physics 326, 96 (2011).
  * Note [1] We note that the bra and ket legs of the operator must be evolved together, without truncation in between. If not, evolution of a Pauli string under a non-trivial Clifford circuit is not guaranteed to yield a single Pauli string.
  * Prosen and Pižorn [2007] T. Prosen and I. Pižorn, Physical Review A 76, 032316 (2007).
  * Note [2] We simply evaluate the expectation value of the doubled operator $Z\otimes Z$ using the vectorized MPO as an MPS.
  * Sahu and Swingle [2022] S. Sahu and B. Swingle, “Efficient tensor network simulation of quantum many-body physics on sparse graphs,” (2022), arXiv:2206.04701 [quant-ph] .
  * Verstraete and Cirac [2004] F. Verstraete and J. I. Cirac, arXiv preprint cond-mat/0407066 (2004).
  * Zaletel and Pollmann [2020] M. P. Zaletel and F. Pollmann, Physical Review Letters 124 (2020), 10.1103/physrevlett.124.037201.
  * Zhou  _et al._ [2020] Y. Zhou, E. M. Stoudenmire, and X. Waintal, Physical Review X 10, 041038 (2020).


# Multiple feedback based wavefront shaping method
to retrieve hidden signal

Nazifa Rumman  Electrical, Computer, and Systems Engineering Department,
Rensselaer Polytechnic Institute,Troy, NY, USA  Tianhong Wang  Kaitlin
Jennings  Pascal Bassène  Finn Buldt  Department of Physics, Applied Physics
and Astronomy, Rensselaer Polytechnic Institute, Troy, NY, USA  Moussa N’Gom
Department of Physics, Applied Physics and Astronomy, Rensselaer Polytechnic
Institute, Troy, NY, USA  Center for Ultrafast Optical Sciences, University of
Michigan, Ann Arbor, Michigan 48109<EMAIL_ADDRESS>
###### Abstract

We present an optical wavefront shaping approach that allows tracking and
localization of signal hidden inside or behind a scattering medium. The method
combines traditional feedback based wavefront shaping together with a switch
function, controlled by two different signals. A simple, in transmission
imaging system is used with two detectors: one monitors the speckle signature
and the other tracks the fully hidden signal (e.g., fluorescent beads). The
algorithm initially finds the optimal incident wavefront to maximize light
transmission to generate a focus in the scattering medium. This modulation
process redirects the scattered input signal inducing instantaneous changes in
both monitored signals, which in turn locates the hidden objects. Once the
response from the hidden target becomes distinct, the algorithm switches to
use this signal as the feedback. We provide experimental demonstrations as a
proof of concept of our approach. Potential applications of our method include
extracting information from biological samples and developing non-invasive
diagnosis methods.

††preprint: AIP/123-QED

##  I Introduction

Turbid materials including biological tissue or human skin scatter light and
create complex interference patterns, known as speckle. This prevents light
from being focused inside biological materials or retrieving information from
objects behind scattering media, posing a fundamental limitation to many
information processing and visualization technologies. Biomedical applications
such as photodynamic therapy (PDT) [1] and optogenetics applications [2]
require directing light through skin or extracting information from target
cells. These applications are either constraint to the outer layer of skin or
become invasive through the insertion of an optical fiber [3].
Wavefront shaping techniques using spatial light modulators (SLM) have emerged
as an attractive proposition to address these issues [4, 5]. The SLM modulates
the incident wavefront in order to focus light through a scattering medium
[6], using methodologies based on phase conjugation [7, 8, 9], feedback based
optimization [10, 11], the transmission matrix, [12, 13] or semi-definite
programming [14] approach. The objective of most studies is to demonstrate
focusing through scattering medium by utilizing speckle feedback e.g., speckle
correlations feedback [15] rather than recovering signals/information hidden
by a scattering medium. Nonlinear fluorescent signal has also been used as
feedback to obtain diffraction-limited focusing through turbid samples [16]. A
shaped focus can be utilized to image masked objects with the help of memory
effect, however, this does not work for thick scattering samples [17, 18, 19].
Approaches that only employ fluorescence feedback are often limited by the low
intensity of the feedback signal. Many of these experiments must have the
fluorescent feedback above a threshold signal to noise ratio (SNR) before the
optimization process is initiated [20].

In this paper, we present a new approach that allows tracking, localization,
and optimization of fluorescent signal emitted by objects hidden within and/or
behind a scattering medium. We combine the traditional principle of feedback
scheme based wavefront shaping together with a switch function that utilizes
both the speckle intensity and the masked signal. The feedback signals are
characterized by two different detectors: the first measures the speckle
intensity from the highly scattering medium and the second detector monitors
the signal generated by the hidden targets. The algorithm exploits the fact
that modulating the incident wavefront to generate a focus in the speckle,
redirects the scattered input fields. These redirected light fields will
simultaneously induce a response from objects hidden within or behind the
scattering medium, thereby locating targets. Once the target signal is located
or reaches a threshold intensity, the feedback mechanism now switches to use
the signal from the hidden object, then proceeds to enhance or optimize it.
Our method is not dependent on a high SNR, or memory effect. This technique is
a stride towards identifying and imaging hidden objects within highly
scattering media. This approach will continue to improve biomedical imaging.
More specifically, it can be an ideal biomedical tool for non–destructive,
non–invasive diagnosis.

##  II principle

The algorithm is designed to find and probe a target deep inside a highly
scattering medium. We start with the traditional iterative method utilizing
continuous sequential algorithm [21] but a switch function is added. The
algorithm starts with finding the optimal wavefront for creating a focus at
the detector that monitors the speckle intensity. The pixels of the SLM are
subdivided into N groups of pixels, known as superpixels. The phase of each
superpixel is iterated from $0$ to $2\pi$ to find the phase that maximizes the
feedback signal. These initial optimization processes are taken to generate a
focus. The continuous variation due to the phase modulation can increase the
signal from the hidden objects. After a predetermined number of iterations
($n$), e.g. at minimum the first row of superpixels must be optimized; the
algorithm starts to ‘look’ for this hidden signal and locates the buried
target using the second detector. Once the signal from the target reaches a
threshold intensity, the feedback signal is switched to come from the second
detector. The algorithm now proceeds to enhance the signal from the target.
The localization and enhancement of the target signal can be categorized into
two cases:
$i$) We predetermine the target position where fluorescence could be expected.
The iterative algorithm is then started to generate a focus at the chosen
location using the speckle as feedback. The algorithm switches feedback when
the fluorescence detected is above a defined threshold and then proceeds to
enhance the signal.
$ii$) The algorithm finds the signal: the iterative optimization is started to
generate a focus at a random speckle point. The algorithm independently
monitors fluorescence signal from different targets. After modulating $n$
number of superpixels with speckle feedback the current fluorescence
intensity($I_{n}$) is compared to the initial fluorescence intensity
($I_{int}$). The hidden fluorescent source that is emitting maximum intensity
at this point of the optimization process is located and chosen as the target.
The algorithm now switches to use the fluorescence as feedback to enhance its
signal (see Figure 1).

Figure 1: A block diagram showing the steps of the switching feedback
algorithm: The optimization process is initiated with speckle feedback. If the
fluorescent target is predefined (case $i$), only that target is monitored,
whereas in case $ii$, the algorithm finds the target after a preset number of
iterations. In both cases, the switching condition has to be met to utilize
fluorescence as the feedback. After switching feedback, the subsequent
superpixels are modulated to optimize the fluorescent signal.

For both of these cases, a seamless switch from the speckle to the
fluorescence feedback is permitted, once the hidden signal’s intensity at the
target position ($I_{fluo}(target)$) attains a value of at least $\eta$ times
higher than the initial intensity ($I_{int}(target)$):
$I_{fluo}(target)\geq\eta\times I_{int}(target)$, where $1.1\leq\eta\leq 2$.
The switching parameter, $\eta$, is both sample and case dependent. If the
value of $\eta$ is too small the algorithm will switch feedback before the
target intensity, $I_{fluo}(target)$ is sufficiently higher than the noise. In
such a case the algorithm will not converge. If $\eta$ is too high the
algorithm will not switch its feedback and the algorithm will proceed to
generate a focus at the chosen speckle point[22]. Once the switching
conditions are fulfilled, the algorithm proceeds to adjust the phases of the
remaining superpixels to enhance the masked signal. The speckle is not used as
the feedback for the remainder of the process and is only monitored. A focus
spot in the speckle will also be obtained as a result of the optimization
process as shown in Figure 2.

Figure 2: a) Fluorescent beads hidden by the scattering medium are excited by
the speckle and optimized with wavefront shaping. b) The fluorescence
intensity detected initially is very low, so the phase of each SLM superpixel
is modulated to maximize the speckle feedback, which induces variations in the
monitored fluorescent signal at the same time. c) Once the fluorescent target
has become distinct or gained enough intensity to perform the optimization
process, the algorithm switches to update the remaining superpixels with
fluorescence feedback d) Optimized fluorescence along with a focus in the
speckle (marked by blue rectangle) can be obtained.

##  III Experimental setup

The experimental optical setup is illustrated in Figure 3. The light source is
a continuous wave Helium–Neon laser (HeNe, $\lambda=632.8$ nm). The beam is
expanded and incident on a phase-only SLM (Santec 200). The polarization of
the incident beam is controlled by a half-wave plate (HWP) and a polarizing
beam splitter (PBS). The reflected beam is transmitted through a 4$f$ imaging
system and focused on the sample by a microscope objective (MO1). A second
microscopic objective (MO2) is used to collect the speckle and fluorescence
from the sample. A dichroic mirror (DM) splits the signals emanating from the
sample. The speckle is detected by CAM1 (Thorlabs CS2100M) and the
fluorescence by CAM2 (Andor Sona sCMOS). We use an additional longpass filter
($\lambda_{c}=650$ nm) to ensure CAM2 only detects fluorescence.

Figure 3: The experimental setup. Incident beam ($\lambda=632.8$ nm) passes
through a half-wave plate (HWP) and a polarizing beam splitter (PBS) to the
spatial light modulator (SLM). From the SLM, the beam is focused onto the
scattering sample by a microscope objective (MO1). Speckle and fluorescence
emission are collected by MO2, split by a dichroic mirror (DM), and monitored
by CAM1 and CAM2, respectively. A longpass filter is added before CAM2.

Plain yogurt and pig skin are used as the scattering media and the hidden
objects are fluorescent microspheres (carboxylate-modified polystyrene beads,
0.04 $\mu$m diameter). The beads are excited at $\lambda=632.8$ nm and the
emission peaks at $\lambda=720$ nm. To prepare a 750 $\mu$m thick yogurt
sample, we spread plain white yogurt on two microscope slides and sandwich the
fluorescent beads between them. The pig skin sample is 2300 $\mu$m thick and
is placed directly on a slide with the fluorescent beads randomly dispersed on
the surface.

##  IV Results and discussion

Starting with case $i$), we choose a point $(x_{speck},y_{speck})$ in the
speckle (CAM1) and define a target point $(x_{fluo},y_{fluo})$ in the
fluorescence image (CAM2). The algorithm begins with optimizing the intensity
at ($x_{speck},y_{speck}$) using the feedback from CAM1 and keeps track of the
fluorescence intensity at ($x_{fluo},y_{fluo}$). The speckle pattern created
by the yogurt sample (Figure 4a) has dark patches due to the absorption of the
beads. A point near the darker region in the speckle is chosen as
$(x_{speck},y_{speck})$, which is marked by the white rectangle.
Prior to the optimization process, the fluorescence emitted by the beads is
negligible, as seen in Figure 4b. Once the switching condition
($I_{fluo}(target)\geq\eta\times I_{int}(target)$) is met during speckle
optimization, the algorithm now switches to acquire its feedback signal from
CAM2.
For this experiment, we set the switching parameter $\eta=1.16$. The switching
condition was met after optimizing only 16 superpixels (Figure 4c). As the
feedback switched to CAM2, the intensity variations at ($x_{speck},y_{speck}$)
are no longer important as they do not affect the the remainder of the
optimization process. But the speckle is still monitored by CAM1 and the
optimization process generates a focus spot in the speckle which is indicated
by the blue rectangle in Figure 4d. As seen in Figure 4e and 4f, fluorescence
intensity at ($x_{fluo},y_{fluo}$) increased significantly after optimization.

Figure 4: Optimization of fluorescent signal via switching feedback (yogurt
sample). (a) Initial speckle pattern, (b) fluorescence prior to optimization;
(c) speckle feedback is being used initially for wavefront shaping (d) speckle
and (e) fluorescence after optimization; and (f) intensity enhancement after
switching to fluorescence feedback. The chosen points ($x_{speck},y_{speck}$)
or ($x_{fluo},y_{fluo}$) are at the center of the white boxes and the
resulting focus in the speckle is indicated by the blue rectangle. Figure 5:
Fluorescence optimization at target chosen by the algorithm. The top and
bottom row have results with yogurt and pig skin sample, respectively.
Fluorescence image (a) before and (b) after optimization, (c) fluorescence
intensity enhancement after switching (yogurt sample); Fluorescence emitted by
beads hidden by pig skin (d) before and (e) after optimization, (f)
fluorescence intensity enhancement after choosing the target position and
switching to fluorescence feedback; The targets picked by the algorithm are at
the center of the white boxes.

For case $ii$), we only define the speckle point, $(x_{speck},y_{speck})$.
After optimizing the first row of superpixels with speckle feedback, the
algorithm compares the initial and current fluorescence intensity for all the
beads and chooses the one that exhibited the largest increase in the emitted
fluorescence thus far as the target $(x_{fluo},y_{fluo})$. Once the target is
chosen the speckle optimization continues until the feedback switching
condition is met.
For the data shown in Figure 5, the top row (a$-$c) was obtained with plain
yogurt as the scattering medium, and the bottom row (d$-$f) with pig skin. For
both samples $\eta$ was set to 1.5. In Figure 5a and 5d, it can be seen that
before the optimization process the fluorescence intensity is very low. After
optimization, one can observe in Figure 5b and 5e that the intensity at the
targets chosen by the algorithm have increased significantly. Despite the pig
skin being thicker by an order of magnitude, the algorithm switched feedbacks
after optimizing approximately the same number of superpixels (Figure 5c and
5f).
One of the benefits of case $ii$) over case $i$) is that the fluorescence
optimization is to be successful even with higher values of $\eta$. Unlike
case $i$) no prior knowledge of expected fluorescence position is needed for
case $ii$), however, case $i$) allows one to control the target position.

Contrary to other wavefront shaping methods based on only fluorescence
feedback that require higher SNR to begin the optimization process, this
method is applicable even when there is no fluorescence detected initially.
Additionally, this technique doesn’t require any memory effect. While in most
instances, the enhancement is proportional to the number of iterations[23],
it’s not true for fluorescence optimization as increasing number of iterations
will require more time and eventually result in photobleaching. In our
experimental demonstration the fluorescent emission was collected in
transmission geometry, however, we expect the method to work with collecting
fluorescence in reflection.

##  V Conclusion

In summary, we have proposed a new wavefront shaping method that allows for
the tracking and enhancement of initially weak signal emitted by sources
hidden behind highly scattering media. We have introduced in the algorithm a
smart switching function that toggles between feedback signals. We have
demonstrated two main scenarios to show the flexibility of our approach. In
one case the hidden target position is known, and in the other the algorithm
locates the position of the hidden object and identifies the target. Our
optical method is simple, easy to implement, and flexible enough to be used in
either transmission or in reflection. We believe, this technique will assist
in retrieving information from hidden objects inside or behind turbid media
and is a stride towards the development of non-invasive optical tool for
medical diagnosis.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors declare no conflicts of interest.

### Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

## References

  * Agostinis  _et al._ [2011] P. Agostinis, K. Berg, K. A. Cengel, T. H. Foster, A. W. Girotti, S. O. Gollnick, S. M. Hahn, M. R. Hamblin, A. Juzeniene, D. Kessel, M. Korbelik, J. Moan, P. Mroz, D. Nowis, J. Piette, B. C. Wilson, and J. Golab, “Photodynamic therapy of cancer: an update,” CA: a cancer journal for clinicians 61, 250–281 (2011).
  * Boyden [2011] E. S. Boyden, “A history of optogenetics: the development of tools for controlling brain circuits with light,” F1000 biology reports 3 (2011), 10.3410/B3-11.
  * Pashaie and Falk [2013] R. Pashaie and R. Falk, “Single optical fiber probe for fluorescence detection and optogenetic stimulation,” IEEE Transactions on Biomedical Engineering 60, 268–280 (2013).
  * Park  _et al._ [2018] J. H. Park, Z. Yu, K. Lee, P. Lai, and Y. Park, “Perspective: Wavefront shaping techniques for controlling multiple light scattering in biological tissues: Toward in vivo applications,” APL Photonics 3, 100901 (2018).
  * Yu  _et al._ [2015] H. Yu, J. Park, K. Lee, J. Yoon, K. Kim, S. Lee, and Y. Park, “Recent advances in wavefront shaping techniques for biomedical applications,” Current Applied Physics 15, 632–641 (2015).
  * Vellekoop and Mosk [2007] I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Opt. Lett. 32, 2309–2311 (2007).
  * Yaqoob  _et al._ [2008] Z. Yaqoob, D. Psaltis, M. Feld, and C. Yang, “Optical phase conjugation for turbidity suppresion in biological samples,” Nat. Photonics 116, 110–115 (2008).
  * Papadopoulos  _et al._ [2012] I. N. Papadopoulos, S. Farahi, C. Moser, and D. Psaltis, “Focusing and scanning light through a multimode optical fiber using digital phase conjugation,” Opt. Express 20, 10583–10590 (2012).
  * Liu  _et al._ [2017] Y. Liu, C. Ma, Y. Shen, J. Shi, and L. V. Wang, “Focusing light inside dynamic scattering media with millisecond digital optical phase conjugation,” Optica 4, 280–288 (2017).
  * Stockbridge  _et al._ [2012] C. Stockbridge, Y. Lu, J. Moore, S. Hoffman, R. Paxman, K. Toussaint, and T. Bifano, “Focusing through dynamic scattering media,” Opt. Express 20, 15086–15092 (2012).
  * Conkey  _et al._ [2012] D. B. Conkey, A. N. Brown, A. M. Caravaca-Aguirre, and R. Piestun, “Genetic algorithm optimization for focusing through turbid media in noisy environments,” Opt. Express 20, 4840–4849 (2012).
  * Popoff  _et al._ [2010] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, “Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media,” Phys. Rev. Lett. 104, 100601 (2010).
  * N’Gom  _et al._ [2018] M. N’Gom, T. B. Norris, E. Michielssen, and R. R. Nadakuditi, “Mode control in a multimode fiber through acquiring its transmission matrix from a reference-less optical system,” Opt. Lett. 43, 419–422 (2018).
  * N’Gom  _et al._ [2017] M. N’Gom, M.-B. Lien, N. M. Estakhri, T. B. Norris, E. Michielssen, and R. R. Nadakuditi, “Controlling light transmission through highly scattering media using semi-definite programming as a phase retrieval computation method,” Scientific Reports 7, 2045–2322 (2017).
  * Stern and Katz [2019] G. Stern and O. Katz, “Noninvasive focusing through scattering layers using speckle correlations,” Opt. Lett. 44, 143–146 (2019).
  * Katz  _et al._ [2014] O. Katz, E. Small, Y. Guan, and Y. Silberberg, “Noninvasive nonlinear focusing and imaging through strongly scattering turbid layers,” Optica 1, 170–174 (2014).
  * Hsieh  _et al._ [2010] C.-L. Hsieh, Y. Pu, R. Grange, G. Laporte, and D. Psaltis, “Imaging through turbid layers by scanning the phase conjugated second harmonic radiation from a nanoparticle,” Opt. Express 18, 20723–20731 (2010).
  * Daniel, Oron, and Silberberg [2019] A. Daniel, D. Oron, and Y. Silberberg, “Light focusing through scattering media via linear fluorescence variance maximization, and its application for fluorescence imaging,” Opt. Express 27, 21778–21786 (2019).
  * Yang, Pu, and Psaltis [2014] X. Yang, Y. Pu, and D. Psaltis, “Imaging blood cells through scattering biological tissue using speckle scanning microscopy,” Opt. Express 22, 3405–3413 (2014).
  * Boniface  _et al._ [2019] A. Boniface, B. Blochet, J. Dong, and S. Gigan, “Noninvasive light focusing in scattering media using speckle variance optimization,” Optica 6, 1381–1385 (2019).
  * Vellekoop and Mosk [2008] I. Vellekoop and A. Mosk, “Phase control algorithms for focusing light through turbid media,” Opt. Communications 281, 3071–3080 (2008).
  * Vellekoop [2015] I. M. Vellekoop, “Feedback-based wavefront shaping,” Opt. Express 23, 12189–12206 (2015).
  * Yang  _et al._ [2021] J. Yang, Q. He, L. Liu, Y. Qu, R. Shao, B. Song, and Y. Zhao, “Anti-scattering light focusing by fast wavefront shaping based on multi-pixel encoded digital-micromirror device,” Light: Science & Applications 10, 2047–7538 (2021).