# Heisenberg-limited Hamiltonian learning for interacting bosons

Haoya Li 111These authors contributed equally to this work. Department of
Mathematics, Stanford University, Stanford, CA 94305  Yu Tong 11footnotemark:
1 Institute for Quantum Information and Matter, California Institute of
Technology, Pasadena, CA 91125  Hongkang Ni  Institute for Computational and
Mathematical Engineering, Stanford University, Stanford, CA 94305  Tuvia Gefen
Institute for Quantum Information and Matter, California Institute of
Technology, Pasadena, CA 91125  Lexing Ying  Department of Mathematics,
Stanford University, Stanford, CA 94305  Institute for Computational and
Mathematical Engineering, Stanford University, Stanford, CA 94305

###### Abstract

We develop a protocol for learning a class of interacting bosonic Hamiltonians
from dynamics with Heisenberg-limited scaling. For Hamiltonians with an
underlying bounded-degree graph structure, we can learn all parameters with
root mean squared error $\epsilon$ using $\mathcal{O}(1/\epsilon)$ total
evolution time, which is independent of the system size, in a way that is
robust against state-preparation and measurement error. In the protocol, we
only use bosonic coherent states, beam splitters, phase shifters, and homodyne
measurements, which are easy to implement on many experimental platforms. A
key technique we develop is to apply random unitaries to enforce symmetry in
the effective Hamiltonian, which may be of independent interest.

##  1 Introduction

Many tasks in quantum metrology and quantum sensing can be reduced to the task
of learning the Hamiltonian $H$ of a quantum system, whose evolution is
described by the operator $e^{-iHt}$ [1, 2, 3, 4, 5, 6, 7, 8, 9]. We call this
task _Hamiltonian learning_ , a name that is commonly used in the literature
[10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. Besides quantum
metrology and quantum sensing, Hamiltonian learning is also useful for quantum
device engineering [24, 25, 26, 27, 28, 29], and quantum many-body physics
[10, 11, 30, 31, 32, 33, 34, 35].

Previous works on Hamiltonian learning for many-body quantum systems are
generally subject to the standard quantum limit (SQL), where to estimate the
parameters in the Hamiltonian to precision $\epsilon$,
$\mathcal{O}(\epsilon^{-2})$ samples are required [12, 13, 14, 15, 16, 17, 18,
19, 20, 21, 22, 23, 36]. On the other hand, for simple systems such as those
consisting of a single spin, the Heisenberg limit can be achieved, where to
obtain $\epsilon$ precision, only $\mathcal{O}(\epsilon^{-1})$ total amount of
resources is needed. Achieving the Heisenberg limit requires using quantum-
enhanced protocols that either use $\mathcal{O}(\epsilon^{-1})$ entangled
probes [5, 4, 3] or coherent evolution for $\mathcal{O}(\epsilon^{-1})$ time
[1, 37, 38]. The resources consumed are the number of probes and the length of
time evolution, respectively.

The natural question is, can we achieve the Heisenberg limit for many-body
quantum systems? When applying the existing quantum-enhanced protocols to the
many-body setting, one quickly encounters difficulties. When many entangled
probes are used, one needs many copies of the quantum system with the same
parameters that can evolve simultaneously without interacting with each other.
It is often unclear how one can create these copies, except for certain
scenarios, such as when many probes undergo evolution under the same field
strength. For long coherent time-evolution, the many-body nature of the
quantum systems becomes problematic as subsystems undergo open-system
dynamics, and phenomena such as thermalization prevent local observables from
having enough sensitivity to achieve the Heisenberg limit. One can consider
performing entangled measurements across all parts of the many-body system.
Still, the difficulty in simulating the system makes finding a good
measurement strategy extremely difficult.

Recently, a method was proposed in [39] to perform Hamiltonian learning for
many-body spin systems with Heisenberg-limited scaling. The main technique is
to apply quantum control in the form of random Pauli operators during time
evolution so that the system evolves with an effective Hamiltonian that is
easy to learn and, at the same time, preserves the parameters that one wants
to learn. Another recent work proved that some form of quantum control is
necessary for achieving the Heisenberg limit in this task [40].

The above works are all focused on multi-qubit systems, and Heisenberg-limited
Hamiltonian learning for bosonic systems is relatively less studied. Bosonic
systems, such as superconducting circuits [41, 42], integrated photonic
circuits [43] and optomechanical platforms [44, 45] are widely used for
quantum sensing, communication, and computing [46, 47, 48, 49]. These quantum
applications require efficient calibration [16], and it is thus highly
desirable to develop optimal algorithms for characterizing bosonic
Hamiltonians. For example, quantum computing and sensing with transmons
require learning the energy levels and interactions between the transmons and
microwave resonators.

For bosonic systems, there is a different set of “easy” quantum states,
unitaries, and measurements than for spins. This work assumes that one can
prepare coherent states, apply phase shifters and beam splitters, and perform
the homodyne measurement. We note that although we may use terms from quantum
optics, such as “phase shifters”, we do not constrain our discussion to the
optical setting. Additionally, in our protocol, we do not require any
squeezing, which can be experimentally difficult to implement [50, 51, 52,
53]. Using these resources, we present a protocol to learn a class of
interacting bosonic Hamiltonians with Heisenberg-limited scaling. These
Hamiltonians involve terms that are quadratic or quartic in the creation and
annihilation operators, and are particle-number preserving. The specific form
of the Hamiltonians is given in (1). Our protocol can also tolerate a constant
amount of noise in the state preparation and measurement (SPAM) procedures and
has a small classical post-processing cost.

In our method, we apply random unitaries during time evolution to reshape the
Hamiltonian into an effective Hamiltonian that is easier to learn. This
follows the same high-level idea as [39] but is specifically tailored to the
bosonic setting. Moreover, we can interpret the procedure as enforcing a
target symmetry in the effective Hamiltonian, thus putting constraints on the
dynamics. We believe this technique may be useful for other problems in
quantum simulation as well [54, 55]. In analyzing the deviation from the
effective dynamics, the unboundedness of the bosonic Hamiltonian terms poses a
challenge, as the analysis in [39] requires Hamiltonian terms to be bounded.
We use more involved techniques to overcome this difficulty in Section E.

##  2 Results

In this work, we focus on quantum systems on $N$ bosonic modes forming a
$d$-dimensional lattice, with the Hamiltonian of the form

$H=\sum_{\braket{i,j}}h_{ij}b_{i}^{\dagger}b_{j}+\sum_{i}\omega_{i}b_{i}^{\dagger}b_{i}+\frac{\xi_{i}}{2}\sum_{i}n_{i}(n_{i}-1),$
(1)

where $b_{i}$ ($b_{i}^{\dagger}$) are bosonic annihilation (creation)
operators, and $n_{i}=b_{i}^{\dagger}b_{i}$ are the number opreators.
$\braket{i,j}$ means that the summation is over sites $i,j$ that are adjacent
to each other. $h_{ij}=h_{ji}^{*}$, and each $\xi_{i}$ and $\omega_{i}$ is a
real number. We also assume that $|h_{ij}|,|\omega_{i}|,|\xi_{i}|\leq 1$. This
class of Hamiltonians is relevant for superconducting quantum processors [16],
arrays of coupled cavities [56], and phonon dynamics in ion crystals [57, 58].
We will present a protocol that generates estimates $\hat{h}_{ij}$,
$\hat{\omega}_{i}$, and $\hat{\xi}_{i}$ such that

$\mathbb{E}[|\hat{h}_{ij}-h_{ij}|^{2}],\
\mathbb{E}[|\hat{\omega}_{i}-\omega_{i}|^{2}],\
\mathbb{E}[|\hat{\xi}_{i}-\xi_{i}|^{2}]\leq\epsilon^{2},$ (2)

for all $i$ and $j$. The protocol has the following properties:

  1. 1.

The _total evolution time_ is $\mathcal{O}(\epsilon^{-1})$;

  2. 2.

The number of experiments is
$\mathcal{O}(\operatorname{polylog}(\epsilon^{-1}))$;

  3. 3.

A constant amount of SPAM error can be tolerated.

More precisely, our protocol consists of
$N_{\exp}=\mathcal{O}(\operatorname{polylog}(\epsilon^{-1}))$ experiments,
which we number by $1,2,\cdots,N_{\exp}$. In the $j$th experiment, we will
initialize each bosonic mode in the system in a coherent state, let the system
evolve for time $t_{j}>0$, and perform homodyne measurement for the bosonic
modes. During time evolution, we will apply random beam splitters (on two
modes) or phase shifters (on one mode). The total evolution time is defined to
be $\sum_{j=1}^{N_{\exp}}t_{j}$, which is the amount of time required to run
all the experiments. We assume that after we prepare the initial state and
before we perform the measurement, the system goes through error channels
$\mathcal{E}_{1}$ and $\mathcal{E}_{2}$, which model the SPAM error. If
$\|\mathcal{E}_{1}-\mathcal{I}\|_{\diamond}+\|\mathcal{E}_{2}-\mathcal{I}\|_{\diamond}$
is upper-bounded by a small constant, then our protocol will still be able to
reach arbitrary precision $\epsilon$. Here $\|\cdot\|_{\diamond}$ is the
diamond norm [59], and $\mathcal{I}$ is the identity channel. The precision is
measured by the mean squared error (MSE). We are using the big-$\mathcal{O}$
notation to hide the constants for simplicity, and we note that these
constants never depend on the system size. Our protocol generates
$\mathcal{O}(NN_{\exp})=\mathcal{O}(N\operatorname{polylog}(\epsilon^{-1}))$
classical data and it takes a similar amount of time to process these data to
compute the estimates.

Below we will describe the protocol in detail. We will start with a protocol
to learn a single anharmonic oscillator, which forms the basic building block
for more complex situations.

###  2.1 Learning an anharmonic oscillator

We first consider the simple case in which

$H_{\mathrm{AHO}}=\omega n+\frac{\xi}{2}n(n-1),$ (3)

where $n=b^{\dagger}b$. We want to estimate the coefficients $\omega$ and
$\xi$ with root mean squared error (RMSE) at most $\epsilon$.

This is a quantum sensing problem with two parameters to be estimated. In
quantum sensing, one usually calculates the quantum Cramér-Rao bound (QCRB)
that provides a lower bound on the MSE of unbiased estimators. Because the two
parameters correspond to Hamiltonian terms that commute with each other, the
QCRB scales inverse quadratically with time, allowing us to achieve the
Heisenberg-limited scaling. This bound, however, is valid only for local
estimation where the prior distribution of the estimators is already
concentrated around the exact value. Here we provide an estimation protocol
that achieves this scaling without any prior knowledge of the parameters.

Our protocol builds upon a robust frequency estimation algorithm similar to
the robust phase estimation algorithm proposed in [38] as well as the
alternative version in [60]. In the robust phase estimation algorithm, we
assume that through performing certain experiments that we will specify when
introducing our protocol, we have access to a random variable $Z_{\delta}(t)$
from measurement results, such that $|Z_{\delta}(t)-e^{-i\omega t}|\leq 1$
with probability at least $1-\delta$, and generating such a random variable
requires evolution time $\mathcal{O}(t\log(\delta^{-1}))$. With multiple
samples of this variable for different values of $t$ and $\delta$, we can
generate an estimate of $\omega$ with RMSE at most $\epsilon$ using
$\mathcal{O}(\epsilon^{-1})$ total evolution time. The algorithm proceeds by
iteratively obtaining estimates with increasing accuracy through longer time
evolution until the target precision is achieved. A detailed description of
the algorithm and proof of its correctness can be found in Section A.

We initialize the system in a coherent state
$\ket{\alpha}=e^{-|\alpha|^{2}/2}\sum_{k}(\alpha^{k}/\sqrt{k!})\ket{k}$, and
let the system evolve under the Hamiltonian $H_{\mathrm{AHO}}$. In the end we
perform homodyne measurements with quadrature operators
$X=(b+b^{\dagger})/\sqrt{2}$ and $P=i(b^{\dagger}-b)/\sqrt{2}$ in separate
experiments. With these measurement results we will be able to estimate
$\braket{b}_{\alpha,t}=\braket{\alpha}{e^{iH_{\mathrm{AHO}}t}be^{-iH_{\mathrm{AHO}}t}}{\alpha}$,
which can be exactly computed to be

$\braket{b}_{\alpha,t}=\alpha e^{-|\alpha|^{2}}e^{-i\omega
t}e^{|\alpha|^{2}e^{-i\xi t}}.$ (4)

We perform this calculation in Section B.

Using (4), we can extract the values of $\omega$ and $\xi$ from
$\braket{b}_{\alpha,t}$. For $\omega$, note that
$\braket{b}_{\alpha,t}/\alpha=e^{-i\omega t}+\mathcal{O}(|\alpha|^{2})$, and
therefore we can choose $|\alpha|$ to be below a small constant so that an
estimate for $\braket{b}_{\alpha,t}/\alpha$ will be close to $e^{-i\omega t}$
within some small constant distance, which enables us to apply the robust
frequency estimation algorithm to estimate $\omega$ with RMSE at most
$\epsilon$ using total evolution time $\mathcal{O}(\epsilon^{-1})$.

For $\xi$, we can extract its value by constructing a periodically oscillating
signal through

$e^{-i\xi
t}=\frac{1}{|\alpha_{1}|^{2}-|\alpha_{2}|^{2}}\log\left(\frac{\alpha_{2}\braket{b}_{\alpha_{1},t}}{\alpha_{1}\braket{b}_{\alpha_{2},t}}\right)+1.$
(5)

This enables us to estimate $\xi$ using the robust frequency estimation
algorithm. Note that, once again, $\braket{b}_{\alpha_{1},t}$ and
$\braket{b}_{\alpha_{2},t}$ only need to be estimated to constant precision,
rather than $\epsilon$ precision which would result in an
$\mathcal{O}(\epsilon^{-2})$ scaling that would destroy the Heisenberg-limited
scaling.

In the above procedure, we need to estimate the expectation of $X$ and $P$
operators, which are unbounded operators that can infinitely amplify any error
in the quantum state. Fortunately, we found that we can replace them with
operators $X\mathbf{1}_{\\{|X|\leq M\\}}$ and $P\mathbf{1}_{\\{|P|\leq M\\}}$,
where $\mathbf{1}_{\\{|X|\leq M\\}}=\int_{|x|\leq M}\ket{x}\bra{x}\mathrm{d}x$
and $\mathbf{1}_{\\{|P|\leq M\\}}$ is similarly defined. This means truncating
the eigenvalues of these operators at a threshold $M=\mathcal{O}(1)$. In
practice, we can simply discard any $X$ and $P$ samples that are above the
threshold $M$ to implement the measurement associated with these truncated
operators. This fact, together with the error tolerance in the robust
frequency estimation algorithm, enables us to tolerate a constant amount of
error from SPAM and time evolution. The combined error from all sources should
be below a small constant, which is sufficient for achieving arbitrarily high
precision.

###  2.2 Learning two coupled anharmonic oscillators

Next, we consider a system consisting of two coupled anharmonic oscillators,
where the Hamiltonian is of the following form:

$H=\omega_{1}b_{1}^{\dagger}b_{1}+\omega_{2}b_{2}^{\dagger}b_{2}+h_{12}b_{1}^{\dagger}b_{2}+h_{21}b_{2}^{\dagger}b_{1}+\frac{\xi_{1}}{2}n_{1}(n_{1}-1)+\frac{\xi_{2}}{2}n_{2}(n_{2}-1)$
(6)

The goal is to learn all the coefficients $\omega_{1}$, $\omega_{2}$,
$\xi_{1}$, $\xi_{2}$, and $h_{12}$ ($h_{21}=h^{*}_{12}$).

We first focus on learning the single-mode coefficients $\omega_{1}$,
$\omega_{2}$, $\xi_{1}$, and $\xi_{2}$. To do this, we will insert random
unitaries during time evolution to decouple the bosonic modes from each other.
In other words, the time evolution operator undergoes the following
transformation

$e^{-iHt}\mapsto\prod_{j=1}^{r}U_{j}^{\dagger}e^{-iH\tau}U_{j}=\prod_{j=1}^{r}e^{-iU_{j}^{\dagger}HU_{j}\tau},$
(7)

where the $U_{j}$, $j=1,2,\cdots,r$, are the random beam splitters or phase
shifters that we insert, $r=t/\tau$, and the product goes from right to left.
Each $U_{j}$ is independently drawn from a distribution that we denote by
$\mathcal{D}$. In the limit of $\tau\to 0$, the dynamics can be described by
an effective Hamiltonian

$H_{\mathrm{effective}}=\mathbb{E}_{U\sim\mathcal{D}}U^{\dagger}HU.$ (8)

This can be seen by considering the Taylor expansion of the time-evolved state
in a small time step:

$\displaystyle\mathbb{E}_{U\sim\mathcal{D}}[e^{-iU^{\dagger}HU\tau}\rho
e^{iU^{\dagger}HU\tau}]$
$\displaystyle=\rho-i\tau\mathbb{E}_{U\sim\mathcal{D}}[[U^{\dagger}HU,\rho]]+\mathcal{O}(\tau^{2})$
(9) $\displaystyle=e^{-i\mathbb{E}_{U\sim\mathcal{D}}[U^{\dagger}HU]\tau}\rho
e^{i\mathbb{E}_{U\sim\mathcal{D}}[U^{\dagger}HU]\tau}+\mathcal{O}(\tau^{2}).$

Note that the above is not a rigorous proof, because the
$\mathcal{O}(\tau^{2})$ residue is an unbounded operator. We provide a
rigorous bound of how far the actual dynamics deviate from the limiting
effective dynamics with finite $\tau>0$ in Section E. The above procedure
introduces additional randomness to our protocol, but it does not introduce
any sample complexity overhead, because we only need the final quantum states
to be close in terms of the trace distance.

To learn all the single mode coefficients, we let the unitary $U$ drawn from
the distribution $\mathcal{D}$ be

$U=e^{-i\theta n_{1}},\quad\theta\sim\mathcal{U}([0,2\pi]).$ (10)

Here $\mathcal{U}([0,2\pi])$ is the uniform distribution over $[0,2\pi]$. We
can then compute the effective Hamiltonian

$H_{\mathrm{effective}}=\frac{1}{2\pi}\int_{0}^{2\pi}e^{i\theta
n_{1}}He^{-i\theta
n_{1}}\mathrm{d}\theta=\omega_{1}n_{1}+\omega_{2}n_{2}+\frac{\xi_{1}}{2}n_{1}(n_{1}-1)+\frac{\xi_{2}}{2}n_{2}(n_{2}-1).$
(11)

In other words, the coupling term
$h_{12}b_{1}^{\dagger}b_{2}+h_{21}b_{2}^{\dagger}b_{1}$ is cancelled in the
process, due to the equality $e^{i\theta n_{1}}b_{1}e^{-i\theta
n_{1}}=e^{-i\theta}b_{1}$. We can interpret this procedure as enforcing a
particle number conservation on the first bosonic mode. In the effective
Hamiltonian, the two bosonic modes are no longer coupled together, and
therefore we can apply the learning algorithm described in Section 2.1 to
learn the parameters of the two modes separately. For a more detailed
description of the protocol see Section C.1.

Next, we will learn the coupling coefficient $h_{12}$. We will use the
following unitaries

$U_{x}(\theta)=e^{i\theta(b_{1}^{\dagger}b_{2}+b_{2}^{\dagger}b_{1})},\quad
U_{y}(\theta)=e^{\theta(b_{1}^{\dagger}b_{2}-b_{2}^{\dagger}b_{1})}.$ (12)

Our protocol is based on the observation that under a single-particle basis
rotation, $h_{12}$ can be estimated from the new single-mode coefficients.
More precisely, we let
$\tilde{b}_{1}=U_{y}(\pi/4)b_{1}U_{y}^{\dagger}(\pi/4)$,
$\tilde{b}_{2}=U_{y}(\pi/4)b_{2}U_{y}^{\dagger}(\pi/4)$, and the new bosonic
modes will be related to the old ones through

$\begin{pmatrix}\tilde{b}_{1}\\\
\tilde{b}_{2}\end{pmatrix}=\begin{pmatrix}\cos(\pi/4)&\sin(\pi/4)\\\
-\sin(\pi/4)&\cos(\pi/4)\end{pmatrix}\begin{pmatrix}b_{1}\\\
b_{2}\end{pmatrix}.$ (13)

We will then rewrite the Hamiltonian (6) in terms of $\tilde{b}_{1}$ and
$\tilde{b}_{2}$. The quadratic part of $H$ can be written as
$\tilde{\omega}_{1}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}+\tilde{\omega}_{2}\tilde{b}_{2}^{\dagger}\tilde{b}_{2}+\tilde{h}_{12}\tilde{b}_{1}^{\dagger}\tilde{b}_{2}+\tilde{h}_{21}\tilde{b}_{2}^{\dagger}\tilde{b}_{1}$,
where

$\tilde{\omega}_{1}=\frac{\omega_{1}+\omega_{2}}{2}+\operatorname{Re}h_{12}.$
(14)

Therefore, $\operatorname{Re}h_{12}$ can be estimated if we can learn
$\tilde{\omega}_{1}$. The quartic part becomes more complicated, but the
procedure we describe next will yield an effective Hamiltonian of a simpler
form.

In our protocol for learning $\operatorname{Re}h_{12}$, we will let the random
unitaries $U_{j}$ in (7) be

$U_{j}=U_{x}(-\theta/2),\quad\theta\sim\mathcal{U}([0,2\pi]),$ (15)

where $\mathcal{U}([0,2\pi])$ denotes the uniform distribution on $[0,2\pi]$.
Note that
$e^{-i\theta\tilde{n}_{1}}=e^{-i(\theta/2)(n_{1}+n_{2})}U_{x}(-\theta/2)$
where $\tilde{n}_{1}=\tilde{b}_{1}^{\dagger}\tilde{b}_{1}$, and because the
total particle number $n_{1}+n_{2}$ is conserved, the random unitary
$U_{x}(-\theta/2)$ is equivalent to $e^{-i\theta\tilde{n}_{1}}$ up to a global
phase. This random unitary, as in (11), results in an effective Hamiltonian in
which $\tilde{n}_{1}$ is conserved. The effective Hamiltonian can be written
as the following

$H_{\mathrm{effective}}=\tilde{\omega}_{1}\tilde{n}_{1}+\tilde{\omega}_{2}\tilde{n}_{2}+\frac{\tilde{\xi}_{11}}{2}\tilde{n}_{1}(\tilde{n}_{1}-1)+\frac{\tilde{\xi}_{22}}{2}\tilde{n}_{2}(\tilde{n}_{2}-1)+\tilde{\xi}_{12}\tilde{n}_{1}\tilde{n}_{2}.$
(16)

In this effective Hamiltonian, the two bosonic modes $\tilde{b}_{1}$ and
$\tilde{b}_{2}$ are still coupled through the term
$\tilde{\xi}_{12}\tilde{n}_{1}\tilde{n}_{2}$. However, because the particle
numbers on both modes are conserved, we can simply initialize the system with
no particle on the mode $\tilde{b}_{2}$, and the coupling term will have no
effect. More specifically, the initial state we use is
$U_{y}(\pi/4)\ket{\alpha}\ket{0}$, which is an $\alpha$-eigenstate for
$\tilde{b}_{1}$ and a $0$-eigenstate for $\tilde{b}_{2}$. The effective
Hamiltonian can then be further reduced to

$H_{\mathrm{effective}}^{\prime}=\tilde{\omega}_{1}\tilde{n}_{1}+\frac{\tilde{\xi}_{11}}{2}\tilde{n}_{1}(\tilde{n}_{1}-1).$
(17)

This enables us to learn $\tilde{\omega}_{1}$ using the single-mode protocol
in Section 2.1, which then gives us $\operatorname{Re}h_{12}$ through (14).
When performing homodyne measurement in the end, we also need to apply
$U_{y}(-\pi/4)$ to rotate back to the original single-particle basis. We write
down the quantum state we get right before measurement to summarize the whole
procedure:

$U_{y}\left(-\frac{\pi}{4}\right)\prod_{j=1}^{r}\left(U_{x}\left(\frac{\theta_{j}}{2}\right)e^{-iH\tau}U_{x}\left(-\frac{\theta_{j}}{2}\right)\right)U_{y}\left(\frac{\pi}{4}\right)\ket{\alpha}\ket{0},$
(18)

where all $\theta_{j}$ are independently drawn from the uniform distribution
over $[0,2\pi]$.

The above procedure yields $\operatorname{Re}h_{12}$. For
$\operatorname{Im}h_{12}$, we only need to switch the roles of $U_{x}(\theta)$
and $U_{y}(\theta)$ and go through the same procedure. For a more detailed
discussion, see Section C.2.

###  2.3 Learning an $N$-mode system

So far, we have concerned ourselves with learning small systems with one or
two modes, but the protocol we develop can be easily generalized to $N$-mode
systems. This section will focus on $N$ bosonic modes arranged on a 1D chain.
For the more general situation with a bounded degree graph, e.g.,
$D$-dimensional square lattice, Kagome lattice, etc., see Section D. The
Hamiltonian is described by (1), where the bosonic modes are labeled
$1,2,\cdots,N$, and $i$ and $j$ are adjacent only when $j=i\pm 1$.

For this $N$-mode system, we consider a divide-and-conquer approach. We will
apply random unitaries so that in the effective dynamics, the system is
divided into clusters of one or two modes, each of which does not interact
with the rest of the system. In this way, we can learn the parameters
associated with each cluster independently and in parallel using our protocol
in Section 2.2.

More specifically, we apply random unitaries in the same way as described in
(7). The random unitary $U_{j}$ is first chosen to be

$U_{j}=\prod_{k=1}^{\lfloor N/3\rfloor}e^{-i\theta_{3k}n_{3k}},$ (19)

where the random variables $\theta_{3k}$ are independently drawn from
$\mathcal{U}([0,2\pi])$, the uniform distribution over $[0,2\pi]$. Randomly
applying the unitaries from this distribution enforces particle number
conservation on sites with indices that are integer multiples of $3$.
Therefore, any Hamiltonian term $b_{i}^{\dagger}b_{j}$ that involves sites
$3,6,9,\cdots$ are canceled. The effective Hamiltonian is

$\displaystyle H$
$\displaystyle=\omega_{1}n_{1}+\omega_{2}n_{2}+h_{12}b_{1}^{\dagger}b_{2}+h_{21}b_{2}^{\dagger}b_{1}$
(20)
$\displaystyle+\omega_{4}n_{4}+\omega_{5}n_{5}+h_{45}b_{4}^{\dagger}b_{5}+h_{54}b_{5}^{\dagger}b_{4}$
$\displaystyle+\cdots$
$\displaystyle+\sum_{i}\frac{\xi_{i}}{2}n_{i}(n_{i}-1),$

where we did not include the terms $\omega_{3}n_{3}$, $\omega_{6}n_{6}$, etc.,
because they only contribute a global phase. In this Hamiltonian, the two
modes $1$ and $2$ form a cluster: they only interact with each other but not
with the rest of the system. The same is true for modes 5 and 6, 7 and 8, etc.
We can then apply the two-mode protocol in Section 2.2 to learn all
coefficients associated with modes 1, 2, 5, 6, … Note that coefficients
associated with different clusters can be learned in parallel in the same
experiment.

Other coefficients remain to learn, such as $\omega_{3}$, $h_{23}$, and
$h_{34}$. We can adopt the same strategy but choose the random unitary
$U_{j}=\prod_{k=0}^{\lfloor N/3\rfloor-1}e^{-i\theta_{3k+1}n_{3k+1}}$ so that
modes 2 and 3, 5 and 6, etc. now form clusters. Similarly, we can let modes 3
and 4, 6 and 7, etc., form clusters. In this way, we can learn all the
coefficients in the Hamiltonian using three different clustering schemes. The
total evolution time required for carrying out all experiments will only be
three times the cost of a two-mode protocol because different clusters can be
learned in parallel.

More generally, we consider a system whose interaction can be described by a
bounded-degree graph. We can design similar clustering schemes based on an
appropriate coloring of its link graph, i.e., the graph whose vertices are the
edges of the original graph. The overhead introduced will be quadratic in the
degree of the original graph and independent of the system size $N$. This is
discussed in more detail in Section D.

##  3 Discussion

In this work, we propose a protocol to learn a class of interacting bosonic
Hamiltonians with Heisenberg-limited scaling. Our protocol uses only elements
of linear optics that that can be implemented on various experimental
platforms. Besides achieving the Heisenberg-limited scaling, our protocol can
also tolerate a constant amount of SPAM noise thanks to the robust frequency
estimation subroutine discussed in Section 2.1. As a part of the protocol, we
also propose a method to enforce symmetry on the effective Hamiltonian
governing the system’s evolution as discussed in more detail in Section 4.1.

To our knowledge, our work is the first to propose a method that learns
interacting bosonic Hamiltonians with Heisenberg-limited scaling in a scalable
way. However, many open problems remain to be solved in this research
direction. In this work, we only consider the particle-number preserving
Hamiltonian in (1), but realistic Hamiltonians may contain terms that do not
preserve the particle number, such as the coupling term in the Jaynes–Cummings
model [61] and capacitive and inductive couplings between superconducting
circuits [41]. Also, higher-order anharmonic effects beyond the fourth order
may be non-negligible in certain quantum systems.

In our protocol, we need to apply random unitaries with a frequency that
depends on the target precision. For higher precision, the speed of applying
these unitaries will also need to be faster, which may be a problem for
experimental implementation. A possible solution is to use some form of
continuous control as considered in [62, 63, 64, 65, 66, 40]. Moreover, since
our protocol requires letting the system evolve coherently for
$\mathcal{O}(\epsilon^{-1})$ times to reach $\epsilon$ precision, the
achievable precision will be limited by quantum noise such as dephasing and
photon losses that limit the coherence time of most experimental Bosonic
systems. It would be therefore interesting to explore whether noise
suppression techniques such as dynamical decoupling [67, 68, 69, 70] and
quantum error correction [71, 72, 73, 74, 75, 76] can mitigate this limitation
and whether they can be incorporated into our protocol in a useful and
scalable way.

Random Clifford unitaries played a crucial role in the classical shadow
formalism [35, 77] as well as Hamiltonian learning [39]. Similarly, one may
wonder whether the random gaussian unitaries used in this work can be useful
for other quantum information tasks for bosonic systems, such as classical
shadow tomography for continuous-variable systems [78, 79, 80].

##  4 Methods

###  4.1 Enforcing symmetry using random unitaries

This section will describe how to enforce symmetry using random unitaries.
This strategy is similar in spirit to the symmetry protection strategies in
[54, 55], but is easier to scale to an $N$-mode system in the current setting.

Let us first consider the general case where we have a compact Lie group $G$
that describes the symmetry we want in the quantum system. Our quantum system
is evolving under a Hamiltonian $H$ that does not necessarily satisfy this
symmetry, i.e., there may exist $g\in G$ such that $gHg^{-1}\neq H$ (here we
equate an element of the Lie group with its matrix representation). We want to
have the system evolve under an effective Hamiltonian $H_{\mathrm{effective}}$
that satisfies the symmetry, i.e.,

$gH_{\mathrm{effective}}g^{-1}=H_{\mathrm{effective}}.$ (21)

We achieve this by inserting random unitaries in the same way as in (7), which
gives us an effective Hamiltonian according to (8). The distribution from
which we draw the random unitaries is the Haar measure on $G$, which we denote
by $\mu$. The effective Hamiltonian can be computed as

$H_{\mathrm{effective}}=\int gHg^{-1}\mu(\mathrm{d}g).$ (22)

When the Hamiltonian $H$ is unbounded, the above equality may only hold in a
weak sense. We can verify that this effective Hamiltonian satisfies the
desired symmetry because

$g^{\prime}H_{\mathrm{effective}}g^{\prime-1}=\int
g^{\prime}gH(g^{\prime}g)^{-1}\mu(\mathrm{d}g)=\int
g^{\prime}gH(g^{\prime}g)^{-1}\mu(\mathrm{d}(g^{\prime}g))=H_{\mathrm{effective}}.$
(23)

Here we have used the property of the Haar measure that
$\mu(\mathrm{d}(g^{\prime}g))=\mu(\mathrm{d}g)$.

It may not be easy to randomly apply elements from the symmetry group $G$.
Still, in our learning protocol, we will only enforce symmetries that are
either U(1) or
$\mathrm{U}(1)\times\mathrm{U}(1)\times\cdots\times\mathrm{U}(1)=U(1)^{\times
N}$, where sampling can easily be done for each U(1) group separately.

##  5 Acknowledgements

The authors thank Matthias Caro for helpful discussions. Y.T. acknowledges
funding from the U.S. Department of Energy Office of Science, Office of
Advanced Scientific Computing Research (DE-NA0003525, and DE-SC0020290). Work
supported by DE-SC0020290 is supported by the DOE QuantISED program through
the theory consortium “Intersections of QIS and Theoretical Particle Physics”
at Fermilab. The Institute for Quantum Information and Matter is an NSF
Physics Frontiers Center. The work of H.L. and L.Y. is partially supported by
National Science Foundation under awards DMS-2011699 and DMS-2208163. T.G.
acknowledges funding provided by the Institute for Quantum Information and
Matter and the Quantum Science and Technology Scholarship of the Israel
Council for Higher Education.

## References

  * de Burgh and Bartlett [2005] Mark de Burgh and Stephen D Bartlett.  Quantum methods for clock synchronization: Beating the standard quantum limit without entanglement.  _Physical Review A_ , 72(4):042301, 2005.
  * Valencia et al. [2004] Alejandra Valencia, Giuliano Scarcelli, and Yanhua Shih.  Distant clock synchronization using entangled photon pairs.  _Applied Physics Letters_ , 85(13):2655–2657, 2004.
  * Leibfried et al. [2004] Dietrich Leibfried, Murray D Barrett, T Schaetz, Joseph Britton, J Chiaverini, Wayne M Itano, John D Jost, Christopher Langer, and David J Wineland.  Toward heisenberg-limited spectroscopy with multiparticle entangled states.  _Science_ , 304(5676):1476–1478, 2004.
  * Bollinger et al. [1996] John J Bollinger, Wayne M Itano, David J Wineland, and Daniel J Heinzen.  Optimal frequency measurements with maximally correlated states.  _Physical Review A_ , 54(6):R4649, 1996.
  * Lee et al. [2002] Hwang Lee, Pieter Kok, and Jonathan P Dowling.  A quantum Rosetta stone for interferometry.  _Journal of Modern Optics_ , 49(14-15):2325–2338, 2002.
  * McKenzie et al. [2002] Kirk McKenzie, Daniel A Shaddock, David E McClelland, Ben C Buchler, and Ping Koy Lam.  Experimental demonstration of a squeezing-enhanced power-recycled michelson interferometer for gravitational wave detection.  _Physical review letters_ , 88(23):231102, 2002\.
  * Holland and Burnett [1993] MJ Holland and K Burnett.  Interferometric detection of optical phase shifts at the heisenberg limit.  _Physical review letters_ , 71(9):1355, 1993.
  * Wineland et al. [1992] David J Wineland, John J Bollinger, Wayne M Itano, FL Moore, and Daniel J Heinzen.  Spin squeezing and reduced quantum noise in spectroscopy.  _Physical Review A_ , 46(11):R6797, 1992.
  * Caves [1981] Carlton M Caves.  Quantum-mechanical noise in an interferometer.  _Physical Review D_ , 23(8):1693, 1981.
  * Wiebe et al. [2014a] Nathan Wiebe, Christopher Granade, Christopher Ferrie, and David Cory.  Quantum hamiltonian learning using imperfect quantum resources.  _Physical Review A_ , 89(4), April 2014a.  ISSN 1094-1622.  doi:10.1103/physreva.89.042314.
  * Wiebe et al. [2014b] Nathan Wiebe, Christopher Granade, Christopher Ferrie, and D. G. Cory.  Hamiltonian learning and certification using quantum resources.  _Physical Review Letters_ , 112(19), May 2014b.  doi:10.1103/physrevlett.112.190501.
  * Li et al. [2020] Zhi Li, Liujun Zou, and Timothy H Hsieh.  Hamiltonian tomography via quantum quench.  _Physical review letters_ , 124(16):160502, 2020\.
  * Che et al. [2021] Liangyu Che, Chao Wei, Yulei Huang, Dafa Zhao, Shunzhong Xue, Xinfang Nie, Jun Li, Dawei Lu, and Tao Xin.  Learning quantum hamiltonians from single-qubit measurements.  _Physical Review Research_ , 3(2):023246, 2021\.
  * Haah et al. [2021] Jeongwan Haah, Robin Kothari, and Ewin Tang.  Optimal learning of quantum hamiltonians from high-temperature gibbs states.  _arXiv preprint arXiv:2108.04842_ , 2021.
  * Yu et al. [2022] Wenjun Yu, Jinzhao Sun, Zeyao Han, and Xiao Yuan.  Practical and efficient hamiltonian learning, 2022.
  * Hangleiter et al. [2021] Dominik Hangleiter, Ingo Roth, Jens Eisert, and Pedram Roushan.  Precise hamiltonian identification of a superconducting quantum processor.  _arXiv preprint arXiv:2108.08319_ , 2021.
  * Franca et al. [2022] Daniel Stilck Franca, Liubov A Markovich, VV Dobrovitski, Albert H Werner, and Johannes Borregaard.  Efficient and robust estimation of many-qubit hamiltonians.  _arXiv preprint arXiv:2205.09567_ , 2022.
  * Zubida et al. [2021] Assaf Zubida, Elad Yitzhaki, Netanel H Lindner, and Eyal Bairey.  Optimal short-time measurements for hamiltonian learning.  _arXiv preprint arXiv:2108.08824_ , 2021.
  * Bairey et al. [2019] Eyal Bairey, Itai Arad, and Netanel H Lindner.  Learning a local hamiltonian from local measurements.  _Physical review letters_ , 122(2):020504, 2019\.
  * Granade et al. [2012] Christopher E Granade, Christopher Ferrie, Nathan Wiebe, and David G Cory.  Robust online hamiltonian learning.  _New Journal of Physics_ , 14(10):103013, 2012\.
  * Gu et al. [2022a] Andi Gu, Lukasz Cincio, and Patrick J Coles.  Practical black box hamiltonian learning.  _arXiv preprint arXiv:2206.15464_ , 2022a.
  * Wilde et al. [2022] Frederik Wilde, Augustine Kshetrimayum, Ingo Roth, Dominik Hangleiter, Ryan Sweke, and Jens Eisert.  Scalably learning quantum many-body hamiltonians from dynamical data, 2022\.  URL https://arxiv.org/abs/2209.14328.
  * Krastanov et al. [2019] Stefan Krastanov, Sisi Zhou, Steven T Flammia, and Liang Jiang.  Stochastic estimation of dynamical variables.  _Quantum Science and Technology_ , 4(3):035003, 2019.
  * Boulant et al. [2003] Nicolas Boulant, Timothy F. Havel, Marco A. Pravia, and David G. Cory.  Robust method for estimating the Lindblad operators of a dissipative quantum process from measurements of the density operator at multiple time points.  _Physical Review A_ , 67(4), April 2003.  ISSN 1050-2947, 1094-1622.  doi:10.1103/PhysRevA.67.042322.
  * Innocenti et al. [2020] Luca Innocenti, Leonardo Banchi, Alessandro Ferraro, Sougato Bose, and Mauro Paternostro.  Supervised learning of time-independent Hamiltonians for gate design.  _New Journal of Physics_ , 22(6), June 2020.  ISSN 1367-2630.  doi:10.1088/1367-2630/ab8aaf.
  * Ben Av et al. [2020] Eitan Ben Av, Yotam Shapira, Nitzan Akerman, and Roee Ozeri.  Direct reconstruction of the quantum-master-equation dynamics of a trapped-ion qubit.  _Physical Review A_ , 101(6), June 2020.  ISSN 2469-9926, 2469-9934.  doi:10.1103/PhysRevA.101.062305.
  * Shulman et al. [2014] Michael D. Shulman, S. P. Harvey, J. M. Nichol, S. D. Bartlett, A. C. Doherty, V. Umansky, and A. Yacoby.  Suppressing qubit dephasing using real-time Hamiltonian estimation.  _Nature Communications_ , 5(1), December 2014.  ISSN 2041-1723.  doi:10.1038/ncomms6156.
  * Sheldon et al. [2016] Sarah Sheldon, Easwar Magesan, Jerry M. Chow, and Jay M. Gambetta.  Procedure for systematically tuning up cross-talk in the cross-resonance gate.  _Physical Review A_ , 93(6), June 2016.  ISSN 2469-9926, 2469-9934.  doi:10.1103/PhysRevA.93.060302.
  * Sundaresan et al. [2020] Neereja Sundaresan, Isaac Lauer, Emily Pritchett, Easwar Magesan, Petar Jurcevic, and Jay M. Gambetta.  Reducing Unitary and Spectator Errors in Cross Resonance with Optimized Rotary Echoes.  _PRX Quantum_ , 1(2), December 2020.  ISSN 2691-3399.  doi:10.1103/PRXQuantum.1.020318.
  * Verdon et al. [2019] Guillaume Verdon, Jacob Marks, Sasha Nanda, Stefan Leichenauer, and Jack Hidary.  Quantum hamiltonian-based models and the variational quantum thermalizer algorithm, 2019.
  * Burgarth and Ajoy [2017] Daniel Burgarth and Ashok Ajoy.  Evolution-Free Hamiltonian Parameter Estimation through Zeeman Markers.  _Physical Review Letters_ , 119(3), July 2017.  ISSN 0031-9007, 1079-7114.  doi:10.1103/PhysRevLett.119.030402.
  * Wang et al. [2017] Jianwei Wang, Stefano Paesani, Raffaele Santagati, Sebastian Knauer, Antonio A. Gentile, Nathan Wiebe, Maurangelo Petruzzella, Jeremy L. O’Brien, John G. Rarity, Anthony Laing, and et al.  Experimental quantum hamiltonian learning.  _Nature Physics_ , 13(6), March 2017.  ISSN 1745-2481.  doi:10.1038/nphys4074.
  * Kwon et al. [2020] Hee Young Kwon, H. G. Yoon, C. Lee, G. Chen, K. Liu, A. K. Schmid, Y. Z. Wu, J. W. Choi, and C. Won.  Magnetic Hamiltonian parameter estimation using deep learning techniques.  _Science Advances_ , 6(39), September 2020.  ISSN 2375-2548.  doi:10.1126/sciadv.abb0872.
  * Wang et al. [2020a] Dingchen Wang, Songrui Wei, Anran Yuan, Fanghua Tian, Kaiyan Cao, Qizhong Zhao, Yin Zhang, Chao Zhou, Xiaoping Song, Dezhen Xue, and Sen Yang.  Machine Learning Magnetic Parameters from Spin Configurations.  _Advanced Science_ , 7(16), August 2020a.  ISSN 2198-3844, 2198-3844.  doi:10.1002/advs.202000566.
  * Huang et al. [2020] Hsin-Yuan Huang, Richard Kueng, and John Preskill.  Predicting many properties of a quantum system from very few measurements.  _Nat. Phys._ , 16:1050––1057, 2020.  URL https://doi.org/10.1038/s41567-020-0932-7.
  * Caro [2022] Matthias C Caro.  Learning quantum processes and hamiltonians via the pauli transfer matrix.  _arXiv preprint arXiv:2212.04471_ , 2022.
  * Higgins et al. [2007] Brendon L Higgins, Dominic W Berry, Stephen D Bartlett, Howard M Wiseman, and Geoff J Pryde.  Entanglement-free heisenberg-limited phase estimation.  _Nature_ , 450(7168):393–396, 2007.
  * Kimmel et al. [2015] Shelby Kimmel, Guang Hao Low, and Theodore J Yoder.  Robust calibration of a universal single-qubit gate set via robust phase estimation.  _Physical Review A_ , 92(6):062315, 2015.
  * Huang et al. [2023] Hsin-Yuan Huang, Yu Tong, Di Fang, and Yuan Su.  Learning many-body hamiltonians with heisenberg-limited scaling.  _Physical Review Letters_ , 130(20):200403, 2023\.
  * Dutkiewicz et al. [2023] Alicja Dutkiewicz, Thomas E O’Brien, and Thomas Schuster.  The advantage of quantum control in many-body hamiltonian learning.  _arXiv preprint arXiv:2304.07172_ , 2023.
  * Krantz et al. [2019] Philip Krantz, Morten Kjaergaard, Fei Yan, Terry P Orlando, Simon Gustavsson, and William D Oliver.  A quantum engineer’s guide to superconducting qubits.  _Applied physics reviews_ , 6(2), 2019.
  * Blais et al. [2021] Alexandre Blais, Arne L Grimsmo, Steven M Girvin, and Andreas Wallraff.  Circuit quantum electrodynamics.  _Reviews of Modern Physics_ , 93(2):025005, 2021\.
  * Wang et al. [2020b] Jianwei Wang, Fabio Sciarrino, Anthony Laing, and Mark G Thompson.  Integrated photonic quantum technologies.  _Nature Photonics_ , 14(5):273–284, 2020b.
  * Aspelmeyer et al. [2014] Markus Aspelmeyer, Tobias J Kippenberg, and Florian Marquardt.  Cavity optomechanics.  _Reviews of Modern Physics_ , 86(4):1391, 2014\.
  * Metcalfe [2014] Michael Metcalfe.  Applications of cavity optomechanics.  _Applied Physics Reviews_ , 1(3), 2014.
  * Blais et al. [2004] Alexandre Blais, Ren-Shou Huang, Andreas Wallraff, Steven M Girvin, and R Jun Schoelkopf.  Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation.  _Physical Review A_ , 69(6):062320, 2004.
  * Clerk et al. [2010] Aashish A Clerk, Michel H Devoret, Steven M Girvin, Florian Marquardt, and Robert J Schoelkopf.  Introduction to quantum noise, measurement, and amplification.  _Reviews of Modern Physics_ , 82(2):1155, 2010\.
  * Chamberland et al. [2022] Christopher Chamberland, Kyungjoo Noh, Patricio Arrangoiz-Arriola, Earl T Campbell, Connor T Hann, Joseph Iverson, Harald Putterman, Thomas C Bohdanowicz, Steven T Flammia, Andrew Keller, et al.  Building a fault-tolerant quantum computer using concatenated cat codes.  _PRX Quantum_ , 3(1):010329, 2022.
  * Adhikari [2014] Rana X Adhikari.  Gravitational radiation detection with laser interferometry.  _Reviews of Modern Physics_ , 86(1):121, 2014\.
  * Qin et al. [2022] Wei Qin, Adam Miranowicz, and Franco Nori.  Beating the 3 db limit for intracavity squeezing and its application to nondemolition qubit readout.  _Physical Review Letters_ , 129(12):123602, 2022\.
  * Dassonneville et al. [2021] R Dassonneville, R Assouly, T Peronnin, AA Clerk, A Bienfait, and B Huard.  Dissipative stabilization of squeezing beyond 3 db in a microwave mode.  _PRX Quantum_ , 2(2):020323, 2021.
  * Kronwald et al. [2013] Andreas Kronwald, Florian Marquardt, and Aashish A Clerk.  Arbitrarily large steady-state bosonic squeezing via dissipation.  _Physical Review A_ , 88(6):063833, 2013.
  * Wiseman and Milburn [1994] HM Wiseman and GJ Milburn.  Squeezing via feedback.  _Physical Review A_ , 49(2):1350, 1994.
  * Nguyen et al. [2022] Nhung H Nguyen, Minh C Tran, Yingyue Zhu, Alaina M Green, C Huerta Alderete, Zohreh Davoudi, and Norbert M Linke.  Digital quantum simulation of the schwinger model and symmetry protection with trapped ions.  _PRX Quantum_ , 3(2):020324, 2022.
  * Halimeh et al. [2021] Jad C Halimeh, Haifeng Lang, Julius Mildenberger, Zhang Jiang, and Philipp Hauke.  Gauge-symmetry protection using single-body terms.  _PRX Quantum_ , 2(4):040311, 2021.
  * Hartmann et al. [2008] Michael J Hartmann, Fernando GSL Brandao, and Martin B Plenio.  Quantum many-body phenomena in coupled cavity arrays.  _Laser & Photonics Reviews_, 2(6):527–556, 2008\.
  * Porras and Cirac [2004] Diego Porras and J Ignacio Cirac.  Bose-einstein condensation and strong-correlation behavior of phonons in ion traps.  _Physical review letters_ , 93(26):263602, 2004\.
  * Shen et al. [2014] Chao Shen, Zhen Zhang, and L-M Duan.  Scalable implementation of boson sampling with trapped ions.  _Physical review letters_ , 112(5):050504, 2014\.
  * Watrous [2018] John Watrous.  _The theory of quantum information_.  Cambridge university press, 2018.
  * Ni et al. [2023] Hongkang Ni, Haoya Li, and Lexing Ying.  On low-depth algorithms for quantum phase estimation.  _arXiv preprint arXiv:2302.02454_ , 2023.
  * Jaynes and Cummings [1963] Edwin T Jaynes and Frederick W Cummings.  Comparison of quantum and semiclassical radiation theories with application to the beam maser.  _Proceedings of the IEEE_ , 51(1):89–109, 1963\.
  * Stark et al. [2017] Alexander Stark, Nati Aharon, Thomas Unden, Daniel Louzon, Alexander Huck, Alex Retzker, Ulrik L Andersen, and Fedor Jelezko.  Narrow-bandwidth sensing of high-frequency fields with continuous dynamical decoupling.  _Nature communications_ , 8(1):1–6, 2017.
  * Gordon et al. [2008] Goren Gordon, Gershon Kurizki, and Daniel A Lidar.  Optimal dynamical decoherence control of a qubit.  _Physical review letters_ , 101(1):010403, 2008\.
  * Fanchini et al. [2007] Felipe Fernandes Fanchini, JEM Hornos, and R d J Napolitano.  Continuously decoupling single-qubit operations from a perturbing thermal bath of scalar bosons.  _Physical Review A_ , 75(2):022329, 2007.
  * Bermudez et al. [2012] A Bermudez, Piet O Schmidt, Martin B Plenio, and A Retzker.  Robust trapped-ion quantum logic gates by continuous dynamical decoupling.  _Physical Review A_ , 85(4):040302, 2012.
  * Xu et al. [2012] Xiangkun Xu, Zixiang Wang, Changkui Duan, Pu Huang, Pengfei Wang, Ya Wang, Nanyang Xu, Xi Kong, Fazhan Shi, Xing Rong, et al.  Coherence-protected quantum gate by continuous dynamical decoupling in diamond.  _Physical review letters_ , 109(7):070502, 2012\.
  * Sekatski et al. [2016] Pavel Sekatski, Michalis Skotiniotis, and Wolfgang Dür.  Dynamical decoupling leads to improved scaling in noisy quantum metrology.  _New Journal of Physics_ , 18(7):073034, 2016\.
  * Schmitt et al. [2017] Simon Schmitt, Tuvia Gefen, Felix M Stürner, Thomas Unden, Gerhard Wolff, Christoph Müller, Jochen Scheuer, Boris Naydenov, Matthew Markham, Sebastien Pezzagna, et al.  Submillihertz magnetic spectroscopy performed with a nanoscale quantum sensor.  _Science_ , 356(6340):832–837, 2017.
  * Arenz et al. [2017] Christian Arenz, Daniel Burgarth, and Robin Hillier.  Dynamical decoupling and homogenization of continuous variable systems.  _Journal of Physics A: Mathematical and Theoretical_ , 50(13):135303, 2017.
  * Heinze and König [2019] Margret Heinze and Robert König.  Universal uhrig dynamical decoupling for bosonic systems.  _Physical review letters_ , 123(1):010501, 2019\.
  * Dür et al. [2014] Wolfgang Dür, Michalis Skotiniotis, Florian Froewis, and Barbara Kraus.  Improved quantum metrology using quantum error correction.  _Physical Review Letters_ , 112(8):080801, 2014\.
  * Arrad et al. [2014] Gilad Arrad, Yuval Vinkler, Dorit Aharonov, and Alex Retzker.  Increasing sensing resolution with error correction.  _Physical review letters_ , 112(15):150801, 2014\.
  * Kessler et al. [2014] Eric M Kessler, Igor Lovchinsky, Alexander O Sushkov, and Mikhail D Lukin.  Quantum error correction for metrology.  _Physical review letters_ , 112(15):150802, 2014\.
  * Herrera-Martí et al. [2015] David A Herrera-Martí, Tuvia Gefen, Dorit Aharonov, Nadav Katz, and Alex Retzker.  Quantum error-correction-enhanced magnetometer overcoming the limit imposed by relaxation.  _Physical review letters_ , 115(20):200501, 2015\.
  * Demkowicz-Dobrzański et al. [2017] Rafał Demkowicz-Dobrzański, Jan Czajkowski, and Pavel Sekatski.  Adaptive quantum metrology under general markovian noise.  _Physical Review X_ , 7(4):041009, 2017.
  * Zhou et al. [2017] Sisi Zhou, Mengzhen Zhang, John Preskill, and Liang Jiang.  Achieving the heisenberg limit in quantum metrology using quantum error correction.  _Nature Communications_ , 9, 2017.
  * Elben et al. [2022] Andreas Elben, Steven T Flammia, Hsin-Yuan Huang, Richard Kueng, John Preskill, Benoît Vermersch, and Peter Zoller.  The randomized measurement toolbox.  _arXiv preprint arXiv:2203.11374_ , 2022.
  * Gu et al. [2022b] Tianren Gu, Xiao Yuan, and Bujiao Wu.  Efficient measurement schemes for bosonic systems.  _arXiv preprint arXiv:2210.13585_ , 2022b.
  * Becker et al. [2022] Simon Becker, Nilanjana Datta, Ludovico Lami, and Cambyse Rouzé.  Classical shadow tomography for continuous variables quantum systems.  _arXiv preprint arXiv:2211.07578_ , 2022.
  * Gandhari et al. [2022] Srilekha Gandhari, Victor V Albert, Thomas Gerrits, Jacob M Taylor, and Michael J Gullans.  Continuous-variable shadow tomography.  _arXiv preprint arXiv:2211.05149_ , 2022.

##  Appendix A Robust frequency estimation

Our main tool to achieve the Heisenberg limit is an algorithm to estimate the
frequency from a complex-valued signal with Heisenberg-limited scaling. This
algorithm resembles the robust phase estimation algorithm in [38] but is
different in that we can deal with cases where the frequency we want is
encoded in the expectation value rather than the probability.

More precisely, we assume access to a signal $Z_{\delta}(t)$ that is close to
$e^{-i\omega t}$, where $|\omega|<W$, by a constant amount of error in both
the phase and the amplitude with probability $1-\delta$, where $\delta$ can be
tuned. It is also reasonable to assume that for smaller $\delta$, generating
the corresponding $Z_{\delta}(t)$ will be more costly, i.e., requiring longer
evolution time. Our algorithm then uses $Z_{\delta}(t)$ for different values
of $\delta$ and $t$ to refine the estimation of $\omega$ iteratively. In each
iteration, we use the result from the previous iteration to get an estimate
$\theta_{j}$ satisfying

$\omega/\widetilde{W}\in\left(\theta_{j}-\frac{\pi}{3\cdot
2^{j}},\theta_{j}+\frac{\pi}{3\cdot 2^{j}}\right),\leavevmode\nobreak\ \mod
2\pi,$ (24)

where $\widetilde{W}=3W/\pi$ is a normalization factor. A detailed algorithm
description can be found in Algorithm 1, adapted from [60].

Algorithm 1 Robust frequency estimation

Input: $\epsilon$: target accuracy, $\delta$: upper bound of the failure
probability, $W$: upper bound for $\omega$.

Let $\widetilde{W}=\frac{3W}{\pi}$ and
$J=\max\\{1,\lceil\log_{2}(4\pi\widetilde{W}/(3\epsilon))\rceil$}.

Set $\theta_{-1}=0$.

for $j=0,1,\ldots,J-1$ do

Set $t_{j}=2^{j}/\widetilde{W}$.

Calculate
$\delta_{j}=\frac{27\epsilon^{2}}{\pi^{2}\widetilde{W}^{2}}\frac{2^{3j-6}}{(2^{J}-1)}$.

Define a candidate set $S_{j}=\left\\{\frac{2k\pi-\arg
Z_{\delta_{j}}(t_{j})}{2^{j}}\right\\}_{k=0,\ldots,2^{j}-1}$.

$\theta_{j}=\operatorname{argmin}_{\theta\in
S_{j}}\left|\theta-\theta_{j-1}\right|_{2\pi}$.

end for

Add a proper integer multiple of $2\pi$ to $\theta_{J-1}$ such that
$|\theta_{J-1}|\leq\pi$.

Output: $\widetilde{W}\theta_{J-1}$ as an approximation to $\omega$.

In the following theorem, we analyze the performance of the above algorithm
for a fixed set of values for $\delta$ in each iteration. We will then
optimize these values to achieve the $\mathcal{O}(\epsilon^{-1})$ scaling in
Corollary 2.

######  Theorem 1.

Suppose that $|\omega|<W$ is known in advance and that we have access to a
signal $Z_{\delta}(t)$ such that

  * •

$|Z_{\delta}(t)|=1$,

  * •

$|Z_{\delta}(t)-e^{-i(\omega t+f(t))}|\leq\eta$ with probability at least
$1-\delta$, where $\sup_{t}|f(t)|\leq C_{f}<{\pi/3}$,

  * •

$2\arcsin\frac{\eta}{2}+C_{f}\leq\frac{\pi}{3}$,

  * •

generating $Z_{\delta}(t)$ requires evolution time
$C_{Z}t(\log(\delta^{-1})+1)$,

then we can produce an estimate $\hat{\omega}$ such that

$\mathbb{E}[|\hat{\omega}-\omega|^{2}]\leq\sum_{j=0}^{J-1}E_{j}^{2}\delta_{j}+\epsilon^{2}/4,$
(25)

with total evolution time at most

$\frac{\pi
C_{Z}}{\widetilde{W}}\sum_{j=0}^{J-1}2^{j}(\log(\delta_{j}^{-1})+1),$ (26)

where

$E_{0}=2\pi,\quad E_{j}=4\pi\widetilde{W}/(3\cdot 2^{j})\ \forall j\geq
1,\quad
J=\lceil\log_{2}(4\pi\widetilde{W}/(3\epsilon))\rceil,\quad\widetilde{W}=3W/\pi,$
(27)

and each $\delta_{j}\in(0,1]$ is arbitrarily chosen.

###### Proof.

Denote $\omega/\widetilde{W}=\frac{\omega\pi}{3W}$ by $\tilde{\omega}$, then
$|\tilde{\omega}|<\pi/3$. We proceed by choosing a sequence of $t_{j}$,
$j=0,1,\ldots,J-1$ and refining the estimation of $\tilde{\omega}$
progressively. In each iteration we generate a signal $Z_{\delta_{j}}(t_{j})$
for arbitrarily chosen $\delta_{j}$.

First, let $t_{0}=1/\widetilde{W}$, then with probability at least
$1-\delta_{0}$,

$|Z_{\delta_{0}}(t_{0})-e^{-i(\tilde{\omega}+f(t_{0}))}|\leq\eta,$

which yields

$|\tilde{\omega}-(-\arg Z_{\delta_{0}}(t_{0}))|\leq
2\arcsin\frac{\eta}{2}+C_{f}<\pi/3\mod 2\pi.$

Thus $\tilde{\omega}\in(-\arg Z_{\delta_{0}}(t_{0})-\pi/3+2k\pi,-\arg
Z_{\delta_{0}}(t_{0})+\pi/3+2k\pi)$ for some integer $k$. Let $\theta_{-1}=0$
and $S_{0}=\\{-\arg Z_{\delta_{0}}(t_{0})\\}$, then by choosing
$\theta_{0}=\operatorname{argmin}_{\theta\in
S_{0}}|\theta-\theta_{-1}|_{2\pi}=-\arg Z_{\delta_{0}}(t_{0})$, we obtain

$\tilde{\omega}\in(\theta_{0}-\frac{\pi}{3},\theta_{0}+\frac{\pi}{3}),\leavevmode\nobreak\
\mod 2\pi.$

Here $\left|\theta\right|_{2\pi}$, is defined to be the minimum distance to
$0$ modulo $2\pi$, i.e., $\left|\theta\right|_{2\pi}=\pi-|(\theta\mod
2\pi)-\pi|$.

At step $j$, we set $t_{j}=2t_{j-1}$, $S_{j}=\left\\{\frac{2k\pi-\arg
Z_{\delta_{j}}(t_{j})}{2^{j}}\right\\}_{k=0,\ldots,2^{j}-1}$, and
$\theta_{j}=\operatorname{argmin}_{\theta\in
S_{j}}\left|\theta-\theta_{j-1}\right|_{2\pi}$. Now we are ready to prove that
if

$|Z_{\delta_{j^{\prime}}}(t_{j^{\prime}})-e^{-i(\tilde{\omega}2^{j^{\prime}}+f(t_{j^{\prime}}))}|\leq\eta,$
(28)

for all $0\leq j^{\prime}\leq j$, then

$\tilde{\omega}\in(\theta_{j}-\frac{\pi}{3\cdot
2^{j}},\theta_{j}+\frac{\pi}{3\cdot 2^{j}}),\leavevmode\nobreak\ \mod 2\pi,$
(29)

for all $j$ by induction. The case $j=0$ is already proved. Suppose (29) holds
for $j-1$. Because of 28, we have

$|\tilde{\omega}2^{j}-(-\arg Z_{\delta_{j}}(t_{j}))|\leq
2\arcsin\frac{\eta}{2}+C_{f}<\pi/3\mod 2\pi.$

Thus

$\tilde{\omega}\in I_{k}:=\left(\frac{2k\pi-\arg
Z_{\delta_{j}}(t_{j})-\pi/3}{2^{j}},\frac{2k\pi-\arg
Z_{\delta_{j}}(t_{j})+\pi/3}{2^{j}}\right)\mod 2\pi,$

for some $k=0,1,\ldots 2^{j}-1$. Notice that
$\operatorname{dist}(I_{k},I_{k^{\prime}})\geq\frac{\pi}{2^{j-1}}-\frac{\pi}{3\cdot
2^{j-1}}=\frac{\pi}{3\cdot 2^{j-2}}$ for any $k\not=k^{\prime}$, and that the
length of the previous estimation $(\theta_{j-1}-\frac{\pi}{3\cdot
2^{j-1}},\theta_{j-1}+\frac{\pi}{3\cdot 2^{j-1}})$ is exactly
$\frac{\pi}{3\cdot 2^{j-2}}$, we can ensure that only one $k^{*}$ satisfies
$I_{k}\cap(\theta_{j-1}-\frac{\pi}{3\cdot
2^{j-1}},\theta_{j-1}+\frac{\pi}{3\cdot 2^{j-1}})\not=\varnothing,\mod 2\pi$.
Moreover, the corresponding $k^{*}$ satisfies $\frac{2k^{*}\pi-\arg
Z_{\delta_{j}}(t_{j})}{2^{j}}=\operatorname{argmin}_{\theta\in
S_{j}}|\theta-\theta_{j-1}|_{2\pi}$, since

$\left|\frac{2k^{*}\pi-\arg
Z_{\delta_{j}}(t_{j})}{2^{j}}-\theta_{j-1}\right|_{2\pi}\leq\left|\frac{2k^{*}\pi-\arg
Z_{\delta_{j}}(t_{j})}{2^{j}}-\tilde{\omega}\right|_{2\pi}+\left|\tilde{\omega}-\theta_{j-1}\right|_{2\pi}<\frac{\pi}{3\cdot
2^{j}}+\frac{\pi}{3\cdot 2^{j-1}}=\frac{\pi}{2^{j}},$

and

$\left|\frac{2k\pi-\arg
Z_{\delta_{j}}(t_{j})}{2^{j}}-\theta_{j-1}\right|_{2\pi}\geq\frac{\pi}{2^{j-1}}-\left|\frac{2k^{*}\pi-\arg
Z_{\delta_{j}}(t_{j})}{2^{j}}-\theta_{j-1}\right|_{2\pi}>\frac{\pi}{2^{j-1}}-\frac{\pi}{2^{j}}=\frac{\pi}{2^{j}}$

for any $k\not=k^{*}$. Now we have proved (29).

In the end, notice that (29) has an ambiguity of modulus $2\pi$, we add a
proper integer multiple of $2\pi$ to $\theta_{J-1}$ such that
$|\theta_{J-1}|\leq\pi$. We then choose this adjusted $\theta_{J-1}$ as our
estimate for $\tilde{\omega}$, and our estimate for $\omega$ is
$\widetilde{W}\theta_{J-1}=:\hat{\omega}$. From the above analysis we can see
that if (28) holds for $0\leq j^{\prime}\leq j-1$, which means that all the
iterations from $0$ to $j-1$ are successful, then by (29), $\tilde{\omega}$ is
contained in $(\theta_{j-1}-\pi/(3\cdot 2^{j-1}),\theta_{j-1}+\pi/(3\cdot
2^{j-1}))+2k\pi$ for some integer $k$, and our estimate $\theta_{J-1}$ is
contained in $(\theta_{j-1}-\pi/(3\cdot 2^{j-1}),\theta_{j-1}+\pi/(3\cdot
2^{j-1}))+2k^{\prime}\pi$ for some integer $k^{\prime}$. Since
$|\tilde{\omega}|<\pi/3$, we have $((\theta_{j-1}-\pi/(3\cdot
2^{j-1}),\theta_{j-1}+\pi/(3\cdot 2^{j-1}))+2k\pi)\subset(-\pi,\pi)$, and then
$((\theta_{j-1}-\pi/(3\cdot 2^{j-1}),\theta_{j-1}+\pi/(3\cdot
2^{j-1}))+2k^{\prime}\pi)\cap[-\pi,\pi]=\varnothing$ if $k^{\prime}\not=k$.
Hence we must have $k=k^{\prime}$ since $|\theta_{J-1}|\leq\pi$. Therefore the
error in the normalized $\tilde{\omega}$ is at most
$E_{j}/\widetilde{W}=4\pi/(3\cdot 2^{j})$, for $j=1,2,\cdots,J-1$. If the very
first iteration fails, the error is at most $E_{0}/\widetilde{W}=2\pi$. If all
the iterations are successful, then by (29) and the argument above,
$|\theta_{J-1}-\tilde{\omega}|\leq\pi/(3\cdot
2^{J-1})\leq\epsilon/(2\widetilde{W})$. From these observations, we will
compute the expected error.

We define the random variable $j_{\mathrm{fail}}$ to be the first iteration
that fails, i.e.,

$|Z_{\delta_{j^{\prime}}}(t_{j^{\prime}})-e^{-i(\tilde{\omega}2^{j^{\prime}}+f(t_{j^{\prime}}))}|\leq\eta,\
\forall
j^{\prime}<j_{\mathrm{fail}},\quad|Z_{\delta_{j_{\mathrm{fail}}}}(t_{j_{\mathrm{fail}}})-e^{-i(\tilde{\omega}2^{j_{\mathrm{fail}}}+f(t_{j_{\mathrm{fail}}}))}|>\eta.$
(30)

If such a $j_{\mathrm{fail}}$ cannot be found, i.e., all iterations are
successful, then we let $j_{\mathrm{fail}}=J$. From the above analysis,
conditional on $j_{\mathrm{fail}}=j<J$, the error will be at most
$E_{j}/\widetilde{W}$. In other words
$\mathbb{E}[|\tilde{\omega}-\theta_{J-1}|^{2}|j_{\mathrm{fail}}=j]\leq
E_{j}^{2}/\widetilde{W}^{2}$. If $j=J$, then the error is at most
$\epsilon/(2\widetilde{W})$. Also, we have

$\Pr[j_{\mathrm{fail}}=j]=(1-\delta_{0})(1-\delta_{1})\cdots(1-\delta_{j-1})\delta_{j}\leq\delta_{j}.$
(31)

Therefore the expected square error is

$\displaystyle\mathbb{E}[|\omega-\hat{\omega}|^{2}]$
$\displaystyle=\widetilde{W}^{2}\mathbb{E}[|\tilde{\omega}-\theta_{J-1}|^{2}]$
(32)
$\displaystyle=\widetilde{W}^{2}\sum_{j=0}^{J}\mathbb{E}[|\tilde{\omega}-\theta_{J-1}|^{2}|j_{\mathrm{fail}}=j]\Pr[j_{\mathrm{fail}}=j]$
$\displaystyle\leq\sum_{j=0}^{J-1}E_{j}^{2}\delta_{j}+\epsilon^{2}/4.$

This proves (25). Generation of each $Z_{\delta_{j}}(t_{j})$ requires an
evolution time of $C_{Z}t_{j}(\log(\delta_{j}^{-1})+1)$, and hence we have
total evolution time (26) by adding them up. ∎

In the theorem above, we have left a great deal of flexibility in choosing
$\delta_{j}$. Below, we will try to answer that if we want the MSE to satisfy
$\mathbb{E}[|\hat{\omega}-\omega|^{2}]\leq\epsilon^{2}$, how we should choose
the $\delta_{j}$ to minimize the total evolution time required. We first state
our result:

######  Corollary 2.

Suppose that $|\omega|<W$ is known in advance and that we have access to a
signal $Z_{\delta}(t)$ such that

  * •

$|Z_{\delta}(t)|=1$,

  * •

$|Z_{\delta}(t)-e^{-i(\omega t+f(t))}|\leq\eta$ with probability at least
$1-\delta$, where $\sup_{t}|f(t)|\leq C_{f}<{\pi/3}$,

  * •

$2\arcsin\frac{\eta}{2}+C_{f}\leq\frac{\pi}{3}$,

  * •

generating $Z_{\delta}(t)$ requires evolution time
$C_{Z}t(\log(\delta^{-1})+1)$,

then we can produce an estimate $\hat{\omega}$ such that
$\mathbb{E}[|\hat{\omega}-\omega|^{2}]\leq\epsilon^{2}$, with total evolution
time at most $\mathcal{O}(C_{Z}\epsilon^{-1})$.

###### Proof.

By (25) and (26), we essentially need to solve the following optimization
problem to get the optimal $\\{\delta_{j}\\}$:

$\displaystyle\underset{\\{\delta_{j}\\}}{\text{minimize}}\
\sum_{j=0}^{J-1}2^{j}\log(\delta_{j}^{-1})$ (33) $\displaystyle\text{subject
to}\ \sum_{j=0}^{J-1}E_{j}^{2}\delta_{j}\leq\frac{3}{4}\epsilon^{2}.$

This optimization problem can be easily solved using the concavity of the
logarithmic function, and the optimal $\delta_{j}$ is

$\delta_{j}=\frac{3\epsilon^{2}}{4E_{j}^{2}}\frac{2^{j}}{(2^{J}-1)}.$ (34)

Using this choice of $\delta_{j}$, we can then compute the total evolution
time required through 26.

$\frac{\pi
C_{Z}}{\widetilde{W}}\sum_{j=0}^{J-1}2^{j}(\log(\delta_{j}^{-1})+1)=\underbrace{\frac{\pi
C_{Z}}{\widetilde{W}}\sum_{j=0}^{J-1}2^{j}\log(\delta_{j}^{-1})}_{\mathrm{(I)}}+\underbrace{\frac{\pi
C_{Z}}{\widetilde{W}}(2^{J}-1)}_{\mathrm{(II)}}.$ (35)

For term $\mathrm{(II)}$, we have

$\frac{\pi C_{Z}}{\widetilde{W}}(2^{J}-1)<\pi C_{Z}\frac{8\pi}{3\epsilon}$
(36)

by our choice of $J$ given in 27. For $\mathrm{(I)}$, using our expression for
$\delta_{j}$ and the expression for $E_{j}$ in (27), we have

$\displaystyle\mathrm{(I)}$ $\displaystyle=\frac{\pi
C_{Z}}{\widetilde{W}}\sum_{j=0}^{J-1}2^{j}\log\left(\frac{4E_{j}^{2}}{3\epsilon^{2}}\frac{2^{J}-1}{2^{j}}\right)$
(37) $\displaystyle=\frac{\pi
C_{Z}}{\widetilde{W}}\log\left(\frac{16\pi^{2}\widetilde{W}^{2}}{3\epsilon^{2}}(2^{J}-1)\right)+\frac{\pi
C_{Z}}{\widetilde{W}}\sum_{j=1}^{J-1}2^{j}\log\left(\frac{64\pi^{2}\widetilde{W}^{2}}{27\epsilon^{2}}\frac{2^{J}-1}{2^{j}}\frac{1}{4^{j}}\right)$
$\displaystyle<\frac{\pi
C_{Z}}{\widetilde{W}}\log\left(\frac{64\pi^{3}\widetilde{W}^{3}}{9\epsilon^{3}}\right)+\frac{\pi
C_{Z}}{\widetilde{W}}\sum_{j=1}^{J-1}2^{j}\log\left(\frac{4}{3}8^{J-j}\right)$
$\displaystyle=\frac{\pi
C_{Z}}{\widetilde{W}}\log\left(\frac{64\pi^{3}\widetilde{W}^{3}}{9\epsilon^{3}}\right)+\frac{\pi
C_{Z}}{\widetilde{W}}\log\left(\frac{4}{3}\right)(2^{J}-2)+\frac{\pi
C_{Z}}{\widetilde{W}}\log(8)(2^{J+2}-2J-2)$
$\displaystyle\leq\mathcal{O}(\epsilon)+\mathcal{O}(C_{Z}\epsilon^{-1}).$

In the last line, we have used the fact that $\epsilon\leq W$ (as otherwise we
can simply estimate $\omega$ by $0$) to bound the first term on the second-to-
last line. Combining (35), (36), and (37), we can see that the total evolution
time of the entire procedure is $\mathcal{O}(C_{Z}\epsilon^{-1})$. ∎

##  Appendix B Learning an anharmonic oscillator

The basic building block of our algorithm is a method to learn a single
anharmonic oscillator of the form

$H_{\mathrm{AHO}}=\omega b^{\dagger}b+\frac{\xi}{2}n(n-1),$ (38)

where $n=b^{\dagger}b$. We will then outline the experiments we run to learn
the coefficients $\omega$ and $\xi$ from this Hamiltonian. We first start with
a coherent state

$\ket{\alpha}=e^{-|\alpha|^{2}/2}\sum_{k=0}^{\infty}\frac{\alpha^{k}}{\sqrt{k!}}\ket{k}.$
(39)

We then let the system evolve under the Hamiltonian $H_{\mathrm{AHO}}$ for
time $t$, and obtain the quantum state

$e^{-iH_{\mathrm{AHO}}t}\ket{\alpha}=e^{-|\alpha|^{2}/2}\sum_{k=0}^{\infty}\frac{\alpha^{k}}{\sqrt{k!}}e^{-i\omega
kt-i\frac{\xi}{2}k(k-1)t}\ket{k}.$ (40)

In the end, we perform POVM measurement in the eigenbasis of either
$X=(b+b^{\dagger})/\sqrt{2}$ or $P=i(b^{\dagger}-b)/\sqrt{2}$, and by taking
average we obtain $\braket{X}_{\alpha,t}$ and $\braket{P}_{\alpha,t}$, where
$\braket{\cdot}_{\alpha,t}$ means taking expectation with respect to the state
$e^{-iH_{\mathrm{AHO}}t}\ket{\alpha}$. With these, we can then obtain the
expectation value of $b$ through

$\braket{b}_{\alpha,t}=\frac{1}{\sqrt{2}}\left(\braket{X}_{\alpha,t}+i\braket{P}_{\alpha,t}\right).$
(41)

The expectation values $\braket{b}_{\alpha,t}$ for a certain set of $\alpha$
and $t$ will enable us to estimate $\omega$ and $\xi$, and we will demonstrate
this below. First we can compute $be^{-iH_{\mathrm{AHO}}t}\ket{\alpha}$ to be

$\displaystyle be^{-iH_{\mathrm{AHO}}t}\ket{\alpha}$
$\displaystyle=e^{-|\alpha|^{2}/2}\sum_{k=1}^{\infty}\frac{\alpha^{k}e^{-i\omega
kt}e^{-i\frac{\xi}{2}k(k-1)t}}{\sqrt{(k-1)!}}\ket{k-1}$ (42)
$\displaystyle=e^{-|\alpha|^{2}/2}\sum_{k=0}^{\infty}\frac{\alpha^{k+1}e^{-i\omega(k+1)t}e^{-i\frac{\xi}{2}k(k+1)t}}{\sqrt{k!}}\ket{k}.$

This yields a closed-form expression for $\braket{b}_{\alpha,t}$:

$\displaystyle\braket{b}_{\alpha,t}$
$\displaystyle=e^{-|\alpha|^{2}}\sum_{k}\frac{\alpha|\alpha|^{2k}e^{-i\omega
t}e^{-i\xi kt}}{k!}$ (43) $\displaystyle=\alpha e^{-|\alpha|^{2}}e^{-i\omega
t}e^{|\alpha|^{2}e^{-i\xi t}}.$

Now we are ready to estimate $\omega$ and $\xi$ with the help of Corollary 2.
To estimate $\omega$, we define

$Z(t)=\frac{\braket{b}_{\alpha,t}}{\left|\braket{b}_{\alpha,t}\right|}=e^{-i(\omega
t+|\alpha|^{2}\sin(\xi t))},$

then $Z(t)=e^{-i(\omega t+f(t))}$, where $f(t)=|\alpha|^{2}\sin(\xi t)$.
Therefore, $\sup_{t}|f(t)|\leq|\alpha|^{2}$. The exact value of $Z(t)$ is,
however, inaccessible in practice, and we need to find an approximation
$Z_{\delta}(t)$ such that $|Z_{\delta}(t)-Z(t)|\leq\eta$ with probability at
least $\delta$ if we want to utilize Corollary 2. In the following, we
decompose the approximation error into three parts and analyze them
separately.

Truncation error. The first part of the approximation error comes from the
truncation of the observables. In our protocol, we truncate the observables up
to a threshold $M$, which means rather than estimating the
$\braket{X}_{\alpha,t}$ and $\braket{P}_{\alpha,t}$ we estimate
$\braket{X\mathbf{1}_{\\{\\}}{X}{\leq M}}{}{\alpha,t}$ and
$\braket{P\mathbf{1}_{\\{\\}}{P}{\leq M}}{}{\alpha,t}$. Here
$\mathbf{1}_{\\{|X|\leq M\\}}$ and $\mathbf{1}_{\\{|P|\leq M\\}}$ are defined
to be

$\mathbf{1}_{\\{|X|\leq
M\\}}=\int_{-M}^{M}\ket{x}\bra{x}\mathrm{d}x,\quad\mathbf{1}_{\\{|P|\leq
M\\}}=\int_{-M}^{M}\ket{p}\bra{p}\mathrm{d}p.$

This is necessary for the robustness of our protocol. With the original
unbounded observables $X$ and $P$, any small error in the quantum state can
potentially be infinitely magnified in the expectation value. The use of
bounded observables will ensure that this does not happen.

In the following, we will ensure that the error introduced by this truncation
is acceptable for our protocol. From Chebyshev’s inequality, one has

$\mathbb{P}(|X|\geq
M)\leq\frac{\braket{X^{2}}_{\alpha,t}}{M^{2}}\leq\frac{2\braket{b^{\dagger}b}_{\alpha,t}+1}{M^{2}}=\frac{2|\alpha|^{2}+1}{M^{2}},$
(44)

where we have used the fact that
$\braket{b^{\dagger}b}_{\alpha,t}=|\alpha|^{2}$. Then, by Cauchy-Schwarz
inequality,

$\displaystyle\left|\braket{X}_{\alpha,t}-\braket{X\mathbf{1}_{\\{\\}}{X}{\leq
M}}{}{\alpha,t}\right|$
$\displaystyle=|\braket{X\mathbf{1}_{\\{\\}}{X}{>M}}{}{\alpha,t}|\leq\sqrt{\braket{X^{2}}_{\alpha,t}}\sqrt{\braket{\mathbf{1}_{\\{\\}}{X}{>M}^{2}}{}{\alpha,t}}$
(45)
$\displaystyle\leq\sqrt{2\braket{b^{\dagger}b}_{\alpha,t}+1}\sqrt{\mathbb{P}(|X|\geq
M)}\leq\frac{2|\alpha|^{2}+1}{{M}}.$

Similarly, one has

$\mathbb{P}(|P|\geq
M)\leq\frac{\braket{P^{2}}_{\alpha,t}}{M^{2}}\leq\frac{2\braket{b^{\dagger}b}_{\alpha,t}+1}{M^{2}}=\frac{2|\alpha|^{2}+1}{M^{2}},$
(46)

and

$\displaystyle\left|\braket{P}_{\alpha,t}-\braket{P\mathbf{1}_{\\{\\}}{P}{\leq
M}}{}{\alpha,t}\right|\leq\frac{2|\alpha|^{2}+1}{{M}}.$ (47)

Combining the error bounds for $X$ and $P$ truncations, we will have an error
bound for the truncated $b$ operator. Let

$Z_{M}(t)=\frac{1}{\sqrt{2}}\left(\braket{X\mathbf{1}_{\\{\\}}{X}{\leq
M}}{}{\alpha,t}+i\braket{P\mathbf{1}_{\\{\\}}{P}{\leq M}}{}{\alpha,t}\right),$
(48)

then

$\displaystyle\left|Z_{M}(t)-\braket{b}_{\alpha,t}\right|$
$\displaystyle=\sqrt{\left|\operatorname{Re}Z_{M}(t)-\operatorname{Re}\braket{b}_{\alpha,t}\right|^{2}+\left|\operatorname{Im}Z_{M}(t)-\operatorname{Im}\braket{b}_{\alpha,t}\right|^{2}}$
(49)
$\displaystyle=\sqrt{\frac{1}{2}\left(\left|\braket{X}_{\alpha,t}-\braket{X\mathbf{1}_{\\{\\}}{X}{\leq
M}}{}{\alpha,t}\right|^{2}+\left|\braket{P}_{\alpha,t}-\braket{P\mathbf{1}_{\\{\\}}{P}{\leq
M}}{}{\alpha,t}\right|^{2}\right)}$
$\displaystyle\leq\frac{2|\alpha|^{2}+1}{{M}}.$

Simulation error. In practice, the final state we obtained is different from
the ideal state $e^{-iH_{\mathrm{AHO}}t}\ket{\alpha}$. This is because in the
multi-mode situation, $H_{\mathrm{AHO}}$ is the Hamiltonian of the effective
dynamics, which differs from the actual dynamics by a small error. In this
sense, we only _simulate_ the effective dynamics, thus calling this error _the
simulation error_. We denote the expectation with respect to the real final
state obtained by $\braket{\cdot}_{\alpha,t,r}$, where $r$ stands for the
parameters used in the simulation. More precisely, as will be explained in
Section C.1, and in particular (76), $r$ is the number of random unitaries
that we insert during the time evolution. In Section E, we will show that for
any given $\eta_{0}>0$, there exists a choice of $r$ such that

$\left|\braket{O}_{\alpha,t,r}-\braket{O}_{\alpha,t}\right|\leq\left\|O\right\|\eta_{0},$
(50)

for any bounded observable $O$. In particular, for any given $\eta_{0}>0$,
there is a choice of $r$ such that

$\left|\braket{X\mathbf{1}_{\\{\\}}{X}{\leq
M}}{}{\alpha,t,r}-\braket{X\mathbf{1}_{\\{\\}}{X}{\leq
M}}{}{\alpha,t}\right|\leq M\eta_{0},\leavevmode\nobreak\
\left|\braket{P\mathbf{1}_{\\{\\}}{P}{\leq
M}}{}{\alpha,t,r}-\braket{P\mathbf{1}_{\\{\\}}{P}{\leq
M}}{}{\alpha,t}\right|\leq M\eta_{0}.$ (51)

Define

$Z_{M,r}(t)=\frac{1}{\sqrt{2}}\left(\braket{X\mathbf{1}_{\\{\\}}{X}{\leq
M}}{}{\alpha,t,r}+i\braket{P\mathbf{1}_{\\{\\}}{P}{\leq
M}}{}{\alpha,t,r}\right),$ (52)

then

$\displaystyle\left|Z_{M,r}(t)-Z_{M}(t)\right|$ (53)
$\displaystyle=\sqrt{\frac{1}{2}\left(\left|\braket{X\mathbf{1}_{\\{\\}}{X}{\leq
M}}{}{\alpha,t,r}-\braket{X\mathbf{1}_{\\{\\}}{X}{\leq
M}}{}{\alpha,t}\right|^{2}+\left|\braket{P\mathbf{1}_{\\{\\}}{P}{\leq
M}}{}{\alpha,t,r}-\braket{P\mathbf{1}_{\\{\\}}{P}{\leq
M}}{}{\alpha,t}\right|^{2}\right)}$ $\displaystyle\leq M\eta_{0}.$

Statistical error. In practice, homodyne measurement generates samples
corresponding to the quadrature operator $X$. By discarding the samples with
norm larger than $M$, we obtain samples
$\hat{x}_{1},\hat{x}_{2},\cdots,\hat{x}_{L}$ corresponding to
$X\mathbf{1}_{\\{|X|\leq M\\}}$. We then approximate
$\braket{X\mathbf{1}_{\\{\\}}{X}{\leq M}}{}{\alpha,t,r}$ through the average
$\bar{x}=(\hat{x}_{1}+\hat{x}_{2}+\cdots+\hat{x}_{L})/L$. Similarly, we can
generate $\hat{p}_{1},\hat{p}_{2},\cdots,\hat{p}_{L}$ corresponding to
$P\mathbf{1}_{\\{|P|\leq M\\}}$, and use
$\bar{p}=(\hat{p}_{1}+\hat{p}_{2}+\cdots+\hat{p}_{L})/L$ to approximate
$\braket{P\mathbf{1}_{\\{\\}}{P}{\leq M}}{}{\alpha,t,r}$. It is clear that
$\hat{x}$ and $\hat{p}$ are unbiased estimates for
$\braket{X\mathbf{1}_{\\{\\}}{X}{\leq M}}{}{\alpha,t,r}$ and
$\braket{P\mathbf{1}_{\\{\\}}{P}{\leq M}}{}{\alpha,t,r}$. Define

$\bar{Z}=\frac{1}{\sqrt{2}}(\bar{x}+i\bar{p}),$ (54)

then

$\displaystyle\left|\bar{Z}-Z_{M,r}(t)\right|=\sqrt{\frac{1}{2}\left(\left|\bar{x}-\braket{X\mathbf{1}_{\\{\\}}{X}{\leq
M}}{}{\alpha,t,r}\right|^{2}+\left|\bar{p}-\braket{P\mathbf{1}_{\\{\\}}{P}{\leq
M}}{}{\alpha,t,r}\right|^{2}\right)}$ (55)
$\displaystyle\leq\max\left\\{\left|\bar{x}-\braket{X\mathbf{1}_{\\{\\}}{X}{\leq
M}}{}{\alpha,t,r}\right|,\left|\bar{p}-\braket{P\mathbf{1}_{\\{\\}}{P}{\leq
M}}{}{\alpha,t,r}\right|\right\\}.$

Thus by the union bound and Hoeffding’s inequality, we have

$\displaystyle\mathbb{P}\left(|\bar{Z}-Z_{M,r}(t)|\geq{\eta_{1}}\right)$
$\displaystyle\leq\mathbb{P}\left(\left|\bar{x}-\braket{X\mathbf{1}_{\\{\\}}{X}{\leq
M}}{}{\alpha,t,r}\right|\geq\eta_{1}\right)+\mathbb{P}\left(\left|\bar{p}-\braket{P\mathbf{1}_{\\{\\}}{P}{\leq
M}}{}{\alpha,t,r}\right|\geq\eta_{1}\right)$ (56) $\displaystyle\leq
2e^{-\frac{L\eta_{1}^{2}}{2M^{2}}}+2e^{-\frac{L\eta_{1}^{2}}{2M^{2}}}=4e^{-\frac{L\eta_{1}^{2}}{2M^{2}}}.$

Putting the three types of error together, we have

$\displaystyle\left|\bar{Z}-\braket{b}_{\alpha,t}\right|$
$\displaystyle\leq\left|\bar{Z}-Z_{M,r}(t)\right|+\left|Z_{M,r}(t)-Z_{M}(t)\right|+\left|Z_{M}(t)-\braket{b}_{\alpha,t}\right|$
(57) $\displaystyle\leq\eta_{1}+M\eta_{0}+\frac{2|\alpha|^{2}+1}{{M}},$

with probability at least $1-4e^{-\frac{L\eta_{1}^{2}}{2M^{2}}}$. Define

$Z_{\delta}(t)=\frac{\bar{Z}}{\left|\bar{Z}\right|},$ (58)

then

$\displaystyle\left|Z_{\delta}(t)-Z(t)\right|$
$\displaystyle=\left|\frac{\bar{Z}}{\left|\bar{Z}\right|}-\frac{\braket{b}_{\alpha,t}}{\left|\braket{b}_{\alpha,t}\right|}\right|\leq\frac{2\left|\bar{Z}-\braket{b}_{\alpha,t}\right|}{\left|\braket{b}_{\alpha,t}\right|}$
(59)
$\displaystyle=\frac{2\left|\bar{Z}-\braket{b}_{\alpha,t}\right|}{|\alpha|e^{|\alpha|^{2}(\cos(\xi
t)-1)}}\leq
2|\alpha|^{-1}e^{2|\alpha|^{2}}\left|\bar{Z}-\braket{b}_{\alpha,t}\right|.$

Hence,

$\left|Z_{\delta}(t)-Z(t)\right|\leq
2|\alpha|^{-1}e^{2|\alpha|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha|^{2}+1}{{M}}\right),$
(60)

with probability at least $1-4e^{-\frac{L\eta_{1}^{2}}{2M^{2}}}$. In order for
the condition $2\arcsin\frac{\eta}{2}+C_{f}\leq\frac{\pi}{3}$ in Theorem 1 to
hold, we need

$2\arcsin\left(|\alpha|^{-1}e^{2|\alpha|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha|^{2}+1}{{M}}\right)\right)+|\alpha|^{2}\leq\frac{\pi}{3}.$
(61)

In order for $1-4e^{-\frac{L\eta_{1}^{2}}{2M^{2}}}\geq 1-\delta$ to hold, we
need

$L\geq\frac{2M^{2}}{\eta_{1}^{2}}\log\frac{4}{\delta}.$ (62)

In conclusion, we have constructed a signal $Z_{\delta}(t)$ to estimate the
parameter $\omega$ that satisfies the conditions required by Corollary 2.

######  Lemma 3.

Define $Z_{\delta}(t)=\bar{Z}/\left|\bar{Z}\right|$, where
$\bar{Z}=\frac{1}{\sqrt{2}}(\bar{x}+i\bar{p})$, and $(\bar{x},\bar{p})$ are
the average values computed from $L$ measurement results each for
$\braket{X\mathbf{1}_{\\{\\}}{X}{\leq M}}{}{\alpha,t,r}$ and
$\braket{P\mathbf{1}_{\\{\\}}{P}{\leq M}}{}{\alpha,t,r}$, respectively. Here
$\braket{\cdot}_{\alpha,t,r}$ denotes the expectation with respect to the real
final state obtained by a simulation using $r$ randomly inserted unitaries,
which is an approximation of the state $e^{-iH_{\mathrm{AHO}}t}\ket{\alpha}$
satisfying
$\left|\braket{O}_{\alpha,t,r}-\braket{O}_{\alpha,t}\right|\leq\left\|O\right\|\eta_{0}$
for any bounded operator $O$. Then $Z_{\delta}(t)$ satisfies the conditions of
Corollary 2 for the estimation of $\omega$ if

$\displaystyle|\alpha|^{2}<\frac{\pi}{3},\leavevmode\nobreak\
M>\frac{e^{2|\alpha|^{2}}(2|\alpha|^{2}+1)}{|\alpha|\sin(\pi/6-|\alpha|^{2}/2)},$
(63)
$\displaystyle\eta_{0}<\frac{1}{M}\left(|\alpha|e^{-2|\alpha|^{2}}\sin\left(\frac{\pi}{6}-\frac{|\alpha|^{2}}{2}\right)-\frac{2|\alpha|^{2}+1}{{M}}\right),$
$\displaystyle\eta_{1}\leq|\alpha|e^{-2|\alpha|^{2}}\sin\left(\frac{\pi}{6}-\frac{|\alpha|^{2}}{2}\right)-\frac{2|\alpha|^{2}+1}{{M}}-M\eta_{0},$
$\displaystyle L\geq\frac{2M^{2}}{\eta_{1}^{2}}\log\frac{4}{\delta}.$

As a result, $\alpha$, $M$, $\eta_{0}$ and $\eta_{1}$ can be chosen as
$\mathcal{O}(1)$ constants and the total runtime needed in producing
$Z_{\delta}(t)$ is $\mathcal{O}(t(\log(1/\delta)+1))$.

######  Remark 4.

When choosing the parameters in practice, one can follow the order in (63),
i.e., first decide the value of $\alpha$, then choose $M$, $\eta_{0}$,
$\eta_{1}$ and $L$ accordingly.

Next, we will build the signal $Z_{\delta}(t)$ for the estimation of $\xi$
with the help of the result above. In particular, when (57) and (61) hold, one
can deduce that

$|\bar{Z}-\braket{b}_{\alpha,t}|\leq|\alpha|e^{-2|\alpha|^{2}}\sin\left(\frac{\pi}{6}-\frac{|\alpha|^{2}}{2}\right)\leq\frac{1}{2}|\alpha|e^{-2|\alpha|^{2}}\leq\frac{1}{2}|\braket{b}_{\alpha,t}|.$
(64)

We first observe that $\cos(\xi t)$ can be obtained by

$\cos(\xi
t)=\frac{1}{|\alpha|^{2}}\log\frac{\left|\braket{b}_{\alpha,t}\right|}{\alpha}+1.$
(65)

Therefore, when when (57) and (61) hold, the error in the estimation of
$\cos(\xi t)$ caused by using $\bar{Z}$ instead of $\braket{b}_{\alpha,t}$ is

$\displaystyle\left|\left(\frac{1}{|\alpha|^{2}}\log\frac{\left|\braket{b}_{\alpha,t}\right|}{\alpha}+1\right)-\left(\frac{1}{|\alpha|^{2}}\log\frac{\left|\bar{Z}\right|}{\alpha}+1\right)\right|$
(66) $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{1}{|\alpha|^{2}}\left|\log\frac{\left|\bar{Z}\right|}{\left|\braket{b}_{\alpha,t}\right|}\right|=\frac{1}{|\alpha|^{2}}\left|\log\left(1+\frac{\left|\bar{Z}\right|-\left|\braket{b}_{\alpha,t}\right|}{\left|\braket{b}_{\alpha,t}\right|}\right)\right|$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ \frac{2\log
2}{|\alpha|^{2}}\left|\frac{\left|\bar{Z}\right|-\left|\braket{b}_{\alpha,t}\right|}{\left|\braket{b}_{\alpha,t}\right|}\right|$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ 2\log
2|\alpha|^{-3}e^{2|\alpha|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha|^{2}+1}{{M}}\right),$

where we have used the concavity of $\log x$ and (64) in the third line. For
estimating $\sin(\xi t)$, we use two different values for $\alpha$. The ratio
between $\braket{b}_{\alpha_{1},t}$ and $\braket{b}_{\alpha_{2},t}$ is

$\frac{\braket{b}_{\alpha_{1},t}}{\braket{b}_{\alpha_{2},t}}=\frac{\alpha_{1}}{\alpha_{2}}e^{(|\alpha_{2}|^{2}-|\alpha_{1}|^{2})(1-e^{-i\xi
t})}.$ (67)

Let
$Z_{\alpha_{1},\alpha_{2}}=\frac{\braket{b}_{\alpha_{1},t}/\braket{b}_{\alpha_{2},t}}{\left|\braket{b}_{\alpha_{1},t}/\braket{b}_{\alpha_{2},t}\right|}$
and $\beta=\left|\alpha_{2}\right|^{2}-\left|\alpha_{1}\right|^{2}$. Assume
that $\left|\beta\right|<\frac{\pi}{2}$, then

$\displaystyle\sin(\xi t)$
$\displaystyle=\frac{1}{\beta}\arcsin\left(\operatorname{Im}Z_{\alpha_{1},\alpha_{2}}\right),$
(68)

Now we analyze the error in the estimate of $\sin(\xi t)$ caused by
approximation. We assume that (61) holds for both $\alpha_{1}$ and
$\alpha_{2}$, and we condition on the event that

$\displaystyle\left|\bar{Z}_{1}-\braket{b}_{\alpha_{1},t}\right|\leq\eta_{1}+M\eta_{0}+\frac{2|\alpha_{1}|^{2}+1}{{M}},$
(69)
$\displaystyle\left|\bar{Z}_{2}-\braket{b}_{\alpha_{2},t}\right|\leq\eta_{1}+M\eta_{0}+\frac{2|\alpha_{2}|^{2}+1}{{M}}.$

Then

$\displaystyle\left|\operatorname{Im}\frac{\bar{Z}_{1}/\bar{Z}_{2}}{\left|\bar{Z}_{1}/\bar{Z}_{2}\right|}-\operatorname{Im}\frac{\braket{b}_{\alpha_{1},t}/\braket{b}_{\alpha_{2},t}}{\left|\braket{b}_{\alpha_{1},t}/\braket{b}_{\alpha_{2},t}\right|}\right|$
(70) $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\left|\frac{\bar{Z}_{1}/\bar{Z}_{2}}{\left|\bar{Z}_{1}/\bar{Z}_{2}\right|}-\frac{\braket{b}_{\alpha_{1},t}/\braket{b}_{\alpha_{2},t}}{\left|\braket{b}_{\alpha_{1},t}/\braket{b}_{\alpha_{2},t}\right|}\right|$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\frac{2\left|\bar{Z}_{1}/\bar{Z}_{2}-\braket{b}_{\alpha_{1},t}/\braket{b}_{\alpha_{2},t}\right|}{\left|\braket{b}_{\alpha_{1},t}/\braket{b}_{\alpha_{2},t}\right|}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
2\frac{\left|\braket{b}_{\alpha_{2},t}\right|}{\left|\braket{b}_{\alpha_{1},t}\right|}\left[\frac{\left|\bar{Z}_{1}-\braket{b}_{\alpha_{1},t}\right|}{\left|\bar{Z}_{2}\right|}+\frac{\left|\braket{b}_{\alpha_{1},t}\right|}{\left|\braket{b}_{\alpha_{2},t}\right|}\frac{\left|\bar{Z}_{2}-\braket{b}_{\alpha_{2},t}\right|}{\left|\bar{Z}_{2}\right|}\right]$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
4\left[\frac{\left|\bar{Z}_{1}-\braket{b}_{\alpha_{1},t}\right|}{\left|\braket{b}_{\alpha_{1},t}\right|}+\frac{\left|\bar{Z}_{2}-\braket{b}_{\alpha_{2},t}\right|}{\left|\braket{b}_{\alpha_{2},t}\right|}\right]$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
4|\alpha_{1}|^{-1}e^{2|\alpha_{1}|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha_{1}|^{2}+1}{{M}}\right)+4|\alpha_{2}|^{-1}e^{2|\alpha_{2}|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha_{2}|^{2}+1}{{M}}\right).$

Here we have used the fact that
$\left|\bar{Z}\right|\geq\frac{1}{2}\left|\braket{b}_{\alpha_{,}t}\right|$ in
the fifth line, which can be deduced from (61). Now, if we further assume that
$\left|\beta\right|\leq\frac{\pi}{3}$, then since the function $\arcsin$ is
$2$-Lipschitz on $[-\sin\frac{\pi}{3},\sin\frac{\pi}{3}]$, we have

$\displaystyle\left|\frac{1}{\beta}\arcsin\left(\operatorname{Im}\frac{\bar{Z}_{1}/\bar{Z}_{2}}{\left|\bar{Z}_{1}/\bar{Z}_{2}\right|}\right)-\frac{1}{\beta}\arcsin(\operatorname{Im}Z_{\alpha_{1},\alpha_{2}})\right|$
(71) $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\frac{2}{\beta}\left|\operatorname{Im}\frac{\bar{Z}_{1}/\bar{Z}_{2}}{\left|\bar{Z}_{1}/\bar{Z}_{2}\right|}-\operatorname{Im}Z_{\alpha_{1},\alpha_{2}}\right|$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\frac{8}{\beta}\left(|\alpha_{1}|^{-1}e^{2|\alpha_{1}|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha_{1}|^{2}+1}{{M}}\right)+|\alpha_{2}|^{-1}e^{2|\alpha_{2}|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha_{2}|^{2}+1}{{M}}\right)\right)$

Combining (66) and (71), we have

$\displaystyle\left|e^{i\xi
t}-\frac{\hat{c}+i\hat{s}}{\left|\hat{c}+i\hat{s}\right|}\right|\leq 4\log
2|\alpha|^{-3}e^{2|\alpha|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha|^{2}+1}{{M}}\right)$
(72)
$\displaystyle+\frac{16}{\beta}\left(|\alpha_{1}|^{-1}e^{2|\alpha_{1}|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha_{1}|^{2}+1}{{M}}\right)+|\alpha_{2}|^{-1}e^{2|\alpha_{2}|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha_{2}|^{2}+1}{{M}}\right)\right),$

where
$\hat{c}=\frac{1}{|\alpha|^{2}}\log\frac{\left|\bar{Z}\right|}{\alpha}+1$ and
$\hat{s}=\frac{1}{\beta}\arcsin\left(\operatorname{Im}\frac{\bar{Z}_{1}/\bar{Z}_{2}}{\left|\bar{Z}_{1}/\bar{Z}_{2}\right|}\right)$,
and the condition in Corollary 2 reads

$\displaystyle 1\geq 4\log
2|\alpha|^{-3}e^{2|\alpha|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha|^{2}+1}{{M}}\right)$
(73)
$\displaystyle+\frac{16}{\beta}\left(|\alpha_{1}|^{-1}e^{2|\alpha_{1}|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha_{1}|^{2}+1}{{M}}\right)+|\alpha_{2}|^{-1}e^{2|\alpha_{2}|^{2}}\left(\eta_{1}+M\eta_{0}+\frac{2|\alpha_{2}|^{2}+1}{{M}}\right)\right).$

In particular, we can take $\alpha=\alpha_{1}$ and obtain the following
result.

######  Lemma 5.

Define $Z_{\delta}(t)=\frac{\hat{c}+i\hat{s}}{\left|\hat{c}+i\hat{s}\right|}$,
where
$\hat{c}=\frac{1}{|\alpha_{1}|^{2}}\log\frac{\left|\bar{Z}_{1}\right|}{\alpha_{1}}+1$,
$\hat{s}=\frac{1}{\beta}\arcsin\left(\operatorname{Im}\frac{\bar{Z}_{1}/\bar{Z}_{2}}{\left|\bar{Z}_{1}/\bar{Z}_{2}\right|}\right)$,
and $(\bar{Z}_{1},\bar{Z}_{2})$ are defined in the same way as in Lemma 3 for
$\alpha_{1}$ and $\alpha_{2}$, respectively. Then $Z_{\delta}(t)$ satisfies
the conditions of Corollary 2 for the estimation of $\xi$ if

$\displaystyle|\alpha_{1}|^{2}<\frac{\pi}{3},\leavevmode\nobreak\
|\alpha_{2}|^{2}<\frac{\pi}{3},\leavevmode\nobreak\
\beta:=\left||\alpha_{1}|^{2}-|\alpha_{2}|^{2}\right|<\frac{\pi}{2},$ (74)
$\displaystyle M>\left(4\log
2|\alpha_{1}|^{-3}+\frac{16}{\beta}|\alpha_{1}|^{-1}\right)e^{2|\alpha_{1}|^{2}}\left({2|\alpha_{1}|^{2}+1}\right)+\frac{16}{\beta}|\alpha_{2}|^{-1}e^{2|\alpha_{2}|^{2}}\left({2|\alpha_{2}|^{2}+1}\right),$
$\displaystyle\eta_{0}<\frac{M-\left(4\log
2|\alpha_{1}|^{-3}+\frac{16}{\beta}|\alpha_{1}|^{-1}\right)e^{2|\alpha_{1}|^{2}}\left({2|\alpha_{1}|^{2}+1}\right)-\frac{16}{\beta}|\alpha_{2}|^{-1}e^{2|\alpha_{2}|^{2}}\left({2|\alpha_{2}|^{2}+1}\right)}{M^{2}\left(\left(4\log
2|\alpha_{1}|^{-3}+\frac{16}{\beta}|\alpha_{1}|^{-1}\right)e^{2|\alpha_{1}|^{2}}+\frac{16}{\beta}|\alpha_{2}|^{-1}e^{2|\alpha_{2}|^{2}}\right)},$
$\displaystyle\eta_{1}\leq\frac{M-\left(4\log
2|\alpha_{1}|^{-3}+\frac{16}{\beta}|\alpha_{1}|^{-1}\right)e^{2|\alpha_{1}|^{2}}\left({2|\alpha_{1}|^{2}+1+M^{2}\eta_{0}}\right)-\frac{16}{\beta}|\alpha_{2}|^{-1}e^{2|\alpha_{2}|^{2}}\left({2|\alpha_{2}|^{2}+1+M^{2}\eta_{0}}\right)}{M\left(\left(4\log
2|\alpha_{1}|^{-3}+\frac{16}{\beta}|\alpha_{1}|^{-1}\right)e^{2|\alpha_{1}|^{2}}+\frac{16}{\beta}|\alpha_{2}|^{-1}e^{2|\alpha_{2}|^{2}}\right)},$
$\displaystyle L\geq\frac{2M^{2}}{\eta_{1}^{2}}\log\frac{8}{\delta}.$

As a result, $\alpha_{1}$, $\alpha_{2}$, $M$, $\eta_{0}$ and $\eta_{1}$ can be
chosen as $\mathcal{O}(1)$ constants and the total runtime needed in producing
$Z_{\delta}(t)$ is $\mathcal{O}(t(\log(1/\delta)+1))$.

##  Appendix C Learning two coupled anharmonic oscillators

In this section, we consider a system consisting of two coupled anharmonic
oscillators, and the Hamiltonian is of the following form:

$H=\omega_{1}b_{1}^{\dagger}b_{1}+\omega_{2}b_{2}^{\dagger}b_{2}+h_{12}b_{1}^{\dagger}b_{2}+h_{21}b_{2}^{\dagger}b_{1}+\frac{\xi_{1}}{2}n_{1}(n_{1}-1)+\frac{\xi_{2}}{2}n_{2}(n_{2}-1)$
(75)

The goal is to learn all the coefficients $\omega_{1}$, $\omega_{2}$,
$\xi_{1}$, $\xi_{2}$, and $h_{12}$ ($h_{21}=h^{*}_{12}$).

###  C.1 Single-mode coefficients

We first focus on learning the single-mode coefficients $\omega_{1}$,
$\omega_{2}$, $\xi_{1}$, and $\xi_{2}$. To do this, we will insert random
unitaries during time evolution to decouple the bosonic modes from each other.
In other words, the time evolution operator undergoes the following
transformation

$e^{-iHt}\mapsto\prod_{j=1}^{r}U_{j}^{\dagger}e^{-iH\tau}U_{j}=\prod_{j=1}^{r}e^{-iU_{j}^{\dagger}HU_{j}\tau},$
(76)

where the $U_{j}$, $j=1,2,\cdots,r$, are the random linear optics unitaries
that we insert, $r=t/\tau$, and the product goes from right to left. Each
$U_{j}$ is independently drawn from a distribution that we denote by
$\mathcal{D}$. In the limit of $\tau\to 0$, the dynamics can be described by
an effective Hamiltonian

$H_{\mathrm{effective}}=\mathbb{E}_{U\sim\mathcal{D}}U^{\dagger}HU.$ (77)

This can be seen by considering the Taylor expansion of the time-evolved state
in a small time step:

$\displaystyle\mathbb{E}_{U\sim\mathcal{D}}[e^{-iU^{\dagger}HU\tau}\rho
e^{iU^{\dagger}HU\tau}]$
$\displaystyle=\rho-i\tau\mathbb{E}_{U\sim\mathcal{D}}[[U^{\dagger}HU,\rho]]+\mathcal{O}(\tau^{2})$
(78) $\displaystyle=e^{-i\mathbb{E}_{U\sim\mathcal{D}}[U^{\dagger}HU]\tau}\rho
e^{i\mathbb{E}_{U\sim\mathcal{D}}[U^{\dagger}HU]\tau}+\mathcal{O}(\tau^{2}).$

The above is not a rigorous proof because the $\mathcal{O}(\tau^{2})$ residue
is an unbounded operator. We will provide a rigorous bound of how far the
actual dynamics deviate from the limiting effective dynamics with finite
$\tau>0$ in Section E.

To learn all the single mode coefficients, we let the unitary $U$ drawn from
the distribution $\mathcal{D}$ be

$U=e^{-i\theta b_{1}^{\dagger}b_{1}},\quad\theta\sim\mathcal{U}([0,2\pi]).$
(79)

Here $\mathcal{U}([0,2\pi])$ is the uniform distribution over $[0,2\pi]$. We
can then compute the effective Hamiltonian

$H_{\mathrm{effective}}=\frac{1}{2\pi}\int_{0}^{2\pi}e^{i\theta
b_{1}^{\dagger}b_{1}}He^{-i\theta
b_{1}^{\dagger}b_{1}}\mathrm{d}\theta=\omega_{1}b_{1}^{\dagger}b_{1}+\omega_{2}b_{2}^{\dagger}b_{2}+\frac{\xi_{1}}{2}n_{1}(n_{1}-1)+\frac{\xi_{2}}{2}n_{2}(n_{2}-1).$
(80)

In other words, the coupling term
$h_{12}b_{1}^{\dagger}b_{2}+h_{21}b_{2}^{\dagger}b_{1}$ is cancelled in the
process, due to the equality

$\frac{1}{2\pi}\int_{0}^{2\pi}e^{i\theta b_{1}^{\dagger}b_{1}}b_{1}e^{-i\theta
b_{1}^{\dagger}b_{1}}\mathrm{d}\theta=\frac{1}{2\pi}\int_{0}^{2\pi}e^{i\theta}b_{1}\mathrm{d}\theta=0.$
(81)

We can interpret this procedure as enforcing a particle number conservation on
the first bosonic mode.

The effective Hamiltonian has the desirable feature that the two bosonic modes
are no longer coupled together. Therefore we can apply the learning algorithm
described in Section B to learn the parameters of the two modes separately.

###  C.2 The coupling coefficient

Next, we consider learning the coupling coefficient $h_{12}$. We observe that
the coupling term can be transformed into a local one under a single-particle
basis transformation. This is done through the following two operators

$U_{x}(\theta)=e^{i\theta(b_{1}^{\dagger}b_{2}+b_{2}^{\dagger}b_{1})},\quad
U_{y}(\theta)=e^{\theta(b_{1}^{\dagger}b_{2}-b_{2}^{\dagger}b_{1})},$ (82)

which correspond to Pauli-$X$ and $Y$ rotations. They transform the
annihilation operators in the following way

$\begin{pmatrix}U_{x}(\theta)b_{1}U_{x}^{\dagger}(\theta)\\\
U_{x}(\theta)b_{2}U_{x}^{\dagger}(\theta)\end{pmatrix}=\begin{pmatrix}\cos(\theta)&i\sin(\theta)\\\
i\sin(\theta)&\cos(\theta)\end{pmatrix}\begin{pmatrix}b_{1}\\\
b_{2}\end{pmatrix},\
\begin{pmatrix}U_{y}(\theta)b_{1}U_{y}^{\dagger}(\theta)\\\
U_{y}(\theta)b_{2}U_{y}^{\dagger}(\theta)\end{pmatrix}=\begin{pmatrix}\cos(\theta)&\sin(\theta)\\\
-\sin(\theta)&\cos(\theta)\end{pmatrix}\begin{pmatrix}b_{1}\\\
b_{2}\end{pmatrix}.$ (83)

We first perform the Pauli-$Y$ rotation and define

$\tilde{b}_{1}=U_{y}(\pi/4)b_{1}U_{y}^{\dagger}(\pi/4),\quad\tilde{b}_{2}=U_{y}(\pi/4)b_{2}U_{y}^{\dagger}(\pi/4).$
(84)

Through (83) we have

$\begin{pmatrix}b_{1}\\\
b_{2}\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&-1\\\
1&1\end{pmatrix}\begin{pmatrix}\tilde{b}_{1}\\\ \tilde{b}_{2}\end{pmatrix}$
(85)

We will then rewrite the Hamiltonian (75) in terms of $\tilde{b}_{1}$ and
$\tilde{b}_{2}$. The quadratic part of $H$ can be written as

$\tilde{\omega}_{1}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}+\tilde{\omega}_{2}\tilde{b}_{2}^{\dagger}\tilde{b}_{2}+\tilde{h}_{12}\tilde{b}_{1}^{\dagger}\tilde{b}_{2}+\tilde{h}_{21}\tilde{b}_{2}^{\dagger}\tilde{b}_{1},$
(86)

where

$\begin{pmatrix}\tilde{\omega}_{1}&\tilde{h}_{12}\\\
\tilde{h}_{21}&\tilde{\omega}_{2}\end{pmatrix}=\frac{1}{2}\begin{pmatrix}1&1\\\
-1&1\end{pmatrix}\begin{pmatrix}\omega_{1}&h_{12}\\\
h_{21}&\omega_{2}\end{pmatrix}\begin{pmatrix}1&-1\\\ 1&1\end{pmatrix}.$ (87)

In particular, we have

$\tilde{\omega}_{1}=\frac{\omega_{1}+\omega_{2}}{2}+\operatorname{Re}h_{12}.$
(88)

For the quartic part, we have

$\displaystyle\frac{\xi_{1}}{2}n_{1}(n_{1}-1)$
$\displaystyle=\frac{\xi_{1}}{2}b_{1}^{\dagger}b_{1}^{\dagger}b_{1}b_{1}=\sum_{ijkl=1}^{2}\xi^{(1)}_{ijkl}\tilde{b}^{\dagger}_{i}\tilde{b}^{\dagger}_{j}\tilde{b}_{k}\tilde{b}_{l},$
(89) $\displaystyle\frac{\xi_{2}}{2}n_{2}(n_{2}-1)$
$\displaystyle=\frac{\xi_{2}}{2}b_{2}^{\dagger}b_{2}^{\dagger}b_{2}b_{2}=\sum_{ijkl=1}^{2}\xi^{(2)}_{ijkl}\tilde{b}^{\dagger}_{i}\tilde{b}^{\dagger}_{j}\tilde{b}_{k}\tilde{b}_{l}.$

In particular

$\xi^{(1)}_{1111}=\frac{\xi_{1}}{4},\quad\xi^{(2)}_{1111}=\frac{\xi_{2}}{4}.$
(90)

Combining (86) and (89), the Hamiltonian $H$ can be written in terms of
$\tilde{b}_{1}$ and $\tilde{b}_{2}$ as

$H=\tilde{\omega}_{1}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}+\tilde{\omega}_{2}\tilde{b}_{2}^{\dagger}\tilde{b}_{2}+\tilde{h}_{12}\tilde{b}_{1}^{\dagger}\tilde{b}_{2}+\tilde{h}_{21}\tilde{b}_{2}^{\dagger}\tilde{b}_{1}+\sum_{ijkl=1}^{2}(\xi^{(1)}_{ijkl}+\xi^{(2)}_{ijkl})\tilde{b}^{\dagger}_{i}\tilde{b}^{\dagger}_{j}\tilde{b}_{k}\tilde{b}_{l}.$
(91)

The above expression is much more complicated than the original expression in
(75), but we will use random unitaries to produce a much simpler effective
Hamiltonian. This time, the random unitary we use will be

$U=e^{-i\theta\tilde{b}_{1}^{\dagger}\tilde{b}_{1}},\quad\theta\sim\mathcal{U}([0,2\pi]).$
(92)

With the same derivation as in (78), we can obtain the effective Hamiltonian
as $\mathbb{E}[U^{\dagger}HU]$. Note that in conjugating with
$e^{i\theta\tilde{b}_{1}^{\dagger}\tilde{b}_{1}}$, each $\tilde{b}_{1}$ in the
Hamiltonian acquires a phase $e^{i\theta}$, and each $\tilde{b}_{1}^{\dagger}$
acquires a phase $e^{-i\theta}$. If in a Hamiltonian term $\tilde{b}_{1}$ and
$\tilde{b}_{1}^{\dagger}$ do not appear the same number of times, then the
term will acquire a phase $e^{ic\theta}$ with $c\in\\{-2,-1,1,2\\}$, and
integrating over $\theta$ will cancel out this term. For example

$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}e^{i\theta\tilde{b}_{1}^{\dagger}\tilde{b}_{1}}\tilde{b}_{1}^{\dagger}\tilde{b}_{2}e^{-i\theta\tilde{b}_{1}^{\dagger}\tilde{b}_{1}}\mathrm{d}\theta=\frac{1}{2\pi}\int_{0}^{2\pi}e^{i\theta}\tilde{b}_{1}^{\dagger}\tilde{b}_{2}\mathrm{d}\theta=0,$
(93)
$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}e^{i\theta\tilde{b}_{1}^{\dagger}\tilde{b}_{1}}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}^{\dagger}\tilde{b}_{2}\tilde{b}_{2}e^{-i\theta\tilde{b}_{1}^{\dagger}\tilde{b}_{1}}\mathrm{d}\theta=\frac{1}{2\pi}\int_{0}^{2\pi}e^{2i\theta}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}^{\dagger}\tilde{b}_{2}\tilde{b}_{2}\mathrm{d}\theta=0.$

In other words, only the terms that conserve the particle number on the first
bosonic mode are preserved in the effective Hamiltonian. We can then write the
effective Hamiltonian as

$H_{\mathrm{effective}}=\tilde{\omega}_{1}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}+\tilde{\omega}_{2}\tilde{b}_{2}^{\dagger}\tilde{b}_{2}+(\xi^{(1)}_{1111}+\xi^{(2)}_{1111})\tilde{n}_{1}(\tilde{n}_{1}-1)+(A\tilde{n}_{1}+B\tilde{n}_{2}+C)\tilde{n}_{2},$
(94)

where $\tilde{n}_{1}=\tilde{b}_{1}^{\dagger}\tilde{b}_{1}$, and
$\tilde{n}_{2}=\tilde{b}_{2}^{\dagger}\tilde{b}_{2}$.

Recall that our goal is to learn the coupling coefficient $h_{12}$, whose real
part can be derived from $\tilde{\omega}_{1}$, $\omega_{1}$, and $\omega_{2}$
through (88), and $\omega_{1}$, and $\omega_{2}$ can be learned using the
procedure outlined in Section C.1. We, therefore, only need to estimate
$\tilde{\omega}_{1}$ from the effective Hamiltonian.

To do this, we start with a product state $\ket{\alpha}\ket{0}$ on the two
bosonic modes. Then we apply $U_{y}(\pi/4)$ to this state to get the initial
state of our time evolution

$\ket{\Phi(0)}=U_{y}(\pi/4)\ket{\alpha}\ket{0}.$ (95)

This state is the tensor product of the coherent states of $\tilde{b}_{1}$ and
$\tilde{b}_{2}$ because one can verify that, using (84),

$\displaystyle\tilde{b}_{1}\ket{\Phi(0)}=\tilde{b}_{1}U_{y}(\pi/4)\ket{\alpha}\ket{0}$
$\displaystyle=U_{y}(\pi/4)b_{1}\ket{\alpha}\ket{0}=\alpha
U_{y}(\pi/4)\ket{\alpha}\ket{0}$ (96)
$\displaystyle\tilde{b}_{2}\ket{\Phi(0)}=\tilde{b}_{2}U_{y}(\pi/4)\ket{\alpha}\ket{0}$
$\displaystyle=U_{y}(\pi/4)b_{2}\ket{\alpha}\ket{0}=0.$

Because of the above equation, we can see that there is no particle in the
bosonic mode $\tilde{b}_{2}$ in this state $\ket{\Phi(0)}$. As the effective
Hamiltonian in (94) conserves the particle number on both bosonic modes, the
particle number on the mode $\tilde{b}_{2}$ will stay $0$. Consequently, any
term that involves $\tilde{n}_{2}$ will not affect the dynamics. Therefore we
can safely discard these terms and get a new effective Hamiltonian

$H_{\mathrm{effective}}^{\prime}=\tilde{\omega}_{1}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}+(\xi^{(1)}_{1111}+\xi^{(2)}_{1111})\tilde{n}_{1}(\tilde{n}_{1}-1).$
(97)

Note that this Hamiltonian only acts non-trivially on the bosonic mode
$\tilde{b}_{1}$. Therefore we can use the single-mode protocol in Section B to
learn the coefficient $\tilde{\omega}_{1}$. As guaranteed in (96), we start
from the $\alpha$-coherent state for $\tilde{b}_{1}$. In the time evolution,
the expectation value $\braket{\tilde{b}_{1}}$ contains the information to
determine $\tilde{\omega}_{1}$. The expectation value $\braket{\tilde{b}_{1}}$
can be extracted through homodyne measurement with two quadrature operators.
Note that we need to convert this homodyne measurement into homodyne
measurement for $b_{1}$ or $b_{2}$. This can be easily done because
$\tilde{b}_{1}=U_{y}(\pi/4)b_{1}U_{y}^{\dagger}(\pi/4)$. We can therefore
apply the unitary $U_{y}^{\dagger}(\pi/4)$ at the end of the time evolution
and then perform homodyne measurement for $(b_{1}+b_{1}^{\dagger})/\sqrt{2}$
and $i(b_{1}-b_{1}^{\dagger})/\sqrt{2}$, which combined yields the expectation
value $\braket{\tilde{b}_{1}}$.

Let us now briefly summarize the whole procedure. We start from a state
$\ket{\alpha}\ket{0}$, apply $U_{y}(\pi/4)$, let the system evolve for time
$t=r\tau$, while applying random
$e^{-i\theta\tilde{b}_{1}^{\dagger}\tilde{b}_{1}}$ with interval $\tau$, and
in the end apply $U_{y}^{\dagger}(\pi/4)=U_{y}(-\pi/4)$, after which we
perform homodyne measurement for $b_{1}$. The quantum state right before the
measurement is applied is

$U_{y}(-\pi/4)\prod_{j=1}^{r}(e^{i\theta_{j}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}}e^{-iH\tau}e^{-i\theta_{j}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}})U_{y}(\pi/4)\ket{\alpha}\ket{0},$
(98)

for randomly sampled $\theta_{j}$, $j=1,2,\cdots,r$. Note that
$e^{-i\theta_{j}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}}=e^{-i(\theta_{j}/2)(n_{1}+n_{2})}U_{x}(-\theta_{j}/2)$,
and $H$ commute with $n_{1}+n_{2}$ because the particle number is conserved.
We therefore have

$e^{i\theta_{j}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}}e^{-iH\tau}e^{-i\theta_{j}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}}=U_{x}(\theta/2)e^{-iH\tau}U_{x}(-\theta/2).$
(99)

Consequently we can replace all
$e^{-i\theta_{j}\tilde{b}_{1}^{\dagger}\tilde{b}_{1}}$ with
$U_{x}(-\theta_{j}/2)$. The quantum state we get in the end is, therefore

$U_{y}\left(-\frac{\pi}{4}\right)\prod_{j=1}^{r}\left(U_{x}\left(\frac{\theta_{j}}{2}\right)e^{-iH\tau}U_{x}\left(-\frac{\theta_{j}}{2}\right)\right)U_{y}\left(\frac{\pi}{4}\right)\ket{\alpha}\ket{0}.$
(100)

Note that the adjacent $U_{x}(-\theta_{j}/2)$ and $U_{x}(\theta_{j-1}/2)$ can
be merged into $U_{x}(-(\theta_{j}-\theta_{j-1})/2)$, so that we only need to
apply one $X$ rotation in each time step instead of two.

In the above procedure, we estimate $\tilde{\omega}_{1}$, which through (88)
we can estimate $\operatorname{Re}h_{12}$. For $\operatorname{Im}h_{12}$, we
can instead define

$\tilde{b}_{1}=U_{x}(\pi/4)b_{1}U_{x}^{\dagger}(\pi/4),\quad\tilde{b}_{2}=U_{x}(\pi/4)b_{2}U_{x}^{\dagger}(\pi/4),$
(101)

and then (14) will become

$\tilde{\omega}_{1}=\frac{\omega_{1}+\omega_{2}}{2}+\operatorname{Im}h_{12}.$
(102)

We can then change the whole procedure accordingly to estimate
$\operatorname{Im}h_{12}$, and the corresponding state before the measurement
is

$U_{x}\left(-\frac{\pi}{4}\right)\prod_{j=1}^{r}\left(U_{y}\left(-\frac{\theta_{j}}{2}\right)e^{-iH\tau}U_{y}\left(\frac{\theta_{j}}{2}\right)\right)U_{x}\left(\frac{\pi}{4}\right)\ket{\alpha}\ket{0}.$
(103)

##  Appendix D Using a divide-and-conquer approach to learn an $N$-mode system

In this section, we consider the general case, where the Hamiltonian is of the
form:

$H=\sum_{\braket{i,j}}h_{ij}b_{i}^{\dagger}b_{j}+\sum_{i}\omega_{i}b_{i}^{\dagger}b_{i}+\frac{\xi_{i}}{2}\sum_{i}n_{i}(n_{i}-1).$
(104)

We will use a divide-and-conquer approach to learn the coefficients in this
Hamiltonian. Specifically, we will insert random unitaries during time
evolution to decouple the system into clusters containing one or two modes
that do not interact with each other and learn the coefficients in each
cluster in parallel.

We assume that the bosonic modes are arranged on a graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is the set
containing all vertices, each of which corresponds to a bosonic mode, and
$\mathcal{E}$ contains all edges. $\sum_{\braket{i,j}}$ means summation over
all vertices linked by an edge.

We consider decoupling the system with the help of a graph
$\mathcal{L}=(\mathcal{E},\mathcal{E}_{\mathcal{L}})$ that is the link graph
of $\mathcal{G}$. The set $\mathcal{E}$ is the set of all edges in
$\mathcal{G}$, and $\mathcal{E}_{\mathcal{L}}$ is the set of edges of
$\mathcal{L}$, which we will now define. For any two edges
$e,e^{\prime}\in\mathcal{E}$, we have $(e,e^{\prime})\in\mathcal{E}$ if and
only if they share a vertex in $\mathcal{V}$.

Next, we color the graph $\mathcal{L}$ with the following rule: any two
vertices in $\mathcal{E}$ must be colored differently if they are at most
distance $2$ from each other. The number of colors needed for this coloring is
at most $\chi=\operatorname{deg}(\mathcal{L})^{2}+1$, and such a coloring can
be easily found by a greedy algorithm: we can simply color a vertex by any
color that its neighbors or next-neighbors have not used, and such a color is
always available because there are at most $\chi-1$ neighbors and next-
neighbors. For a graph $\mathcal{G}$ with degree $D$,
$\operatorname{deg}(\mathcal{L})\leq 2(D-1)$, and therefore $\chi\leq
4(D-1)^{2}+1$. This coloring yields a decomposition of the edges

$\mathcal{E}=\bigsqcup_{c=1}^{\chi}\mathcal{E}_{c},$ (105)

where $\mathcal{E}_{c}$ is the set of edges with color $c$.

For each color $c=1,2,\cdots,\chi$, we then learn all the coefficients
associated with this color. We denote by $\mathcal{V}_{c}$ all the vertices
(bosonic modes) that are contained in an edge in $\mathcal{E}_{c}$. During
time evolution, we apply random unitaries of the form

$U=\prod_{i\in\mathcal{V}\setminus\mathcal{V}_{c}}e^{-i\theta_{i}b_{i}^{\dagger}b_{i}},\quad\theta_{i}\sim\mathcal{U}([0,2\pi]).$
(106)

Here $\theta_{i}$, $i\in\mathcal{V}\setminus\mathcal{V}_{c}$, are independent
random variables. Following the derivation in (78), we can see that the
effective Hamiltonian is

$H_{\mathrm{effective}}=\prod_{i\in\mathcal{V}\setminus\mathcal{V}_{c}}\left(\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\theta_{i}\right)e^{-i\sum_{i\in\mathcal{V}\setminus\mathcal{V}_{c}}\theta_{i}n_{i}}He^{i\sum_{i\in\mathcal{V}\setminus\mathcal{V}_{c}}\theta_{i}n_{i}}.$
(107)

We can then examine the effect of this transformation on each term. For a term
$b_{k}^{\dagger}b_{l}$, $k\neq l$, if $k$ is in
$\mathcal{V}\setminus\mathcal{V}_{c}$ but $l$ is not, then

$\prod_{i\in\mathcal{V}\setminus\mathcal{V}_{c}}\left(\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\theta_{i}\right)e^{-i\sum_{i\in\mathcal{V}\setminus\mathcal{V}_{c}}\theta_{i}n_{i}}b_{k}^{\dagger}b_{l}e^{i\sum_{i\in\mathcal{V}\setminus\mathcal{V}_{c}}\theta_{i}n_{i}}=\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\omega_{k}e^{i\omega_{k}}b_{k}^{\dagger}b_{l}=0$
(108)

The same is true if $l$ is in $\mathcal{V}\setminus\mathcal{V}_{c}$ but $k$ is
not. When both $k,l\in\mathcal{V}\setminus\mathcal{V}_{c}$, then

$\displaystyle\prod_{i\in\mathcal{V}\setminus\mathcal{V}_{c}}\left(\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\theta_{i}\right)e^{-i\sum_{i\in\mathcal{V}\setminus\mathcal{V}_{c}}\theta_{i}n_{i}}b_{k}^{\dagger}b_{l}e^{i\sum_{i\in\mathcal{V}\setminus\mathcal{V}_{c}}\theta_{i}n_{i}}$
(109)
$\displaystyle=\frac{1}{(2\pi)^{2}}\int_{0}^{2\pi}\mathrm{d}\omega_{k}\int_{0}^{2\pi}\mathrm{d}\omega_{l}e^{i(\omega_{k}-\omega_{l})}b_{k}^{\dagger}b_{l}=0.$

In other words, for any coupling term $b^{\dagger}_{k}b_{l}$, the above
procedure will cancel it out if either $k$ or $l$ is in
$\mathcal{V}\setminus\mathcal{V}_{c}$. All other terms are preserved because
they commute with $n_{i}$ for any $i\in\mathcal{V}\setminus\mathcal{V}_{c}$.
The only possible $b^{\dagger}_{k}b_{l}$ terms left are those with
$k,l\in\mathcal{V}_{c}$. This also means that $(k,l)\in\mathcal{E}_{c}$
because of the following argument: first by definition of $\mathcal{V}_{c}$
there must exists $k^{\prime}$ and $l^{\prime}$ such that
$(k,k^{\prime})\in\mathcal{E}_{c}$ and $(l,l^{\prime})\in\mathcal{E}_{c}$. We
must have $(k,l)\in\mathcal{E}$, as otherwise, this coupling term would not
exist at all. This means that unless $k^{\prime}=l^{\prime}$, the two edges
$(k,k^{\prime})$ and $(l,l^{\prime})$ as vertices in $\mathcal{L}$ are next-
neighbors, which is not allowed in our coloring. Therefore
$k^{\prime}=l^{\prime}$ and we have $(k,l)\in\mathcal{E}_{c}$. Consequently,
the effective Hamiltonian is

$H_{\mathrm{effective}}=\sum_{(i,j)\in\mathcal{E}_{c}}h_{ij}b_{i}^{\dagger}b_{j}+\sum_{i}\omega_{i}b_{i}^{\dagger}b_{i}+\frac{\xi_{i}}{2}\sum_{i}n_{i}(n_{i}-1).$
(110)

Next, we will show that the above Hamiltonian is decoupled into clusters of
sizes at most $2$. We will do this by showing that any bosonic mode $i$
interacts with at most one other bosonic mode in the above Hamiltonian. This
can be proved by contradiction: if $i$ interacts with both $j$ and $k$ in the
above Hamiltonian, then $(i,j)\in\mathcal{E}_{c}$ and
$(i,k)\in\mathcal{E}_{c}$, which makes $(i,j)$ and $(i,k)$ neighbors as
vertices in $\mathcal{L}$, and this is forbidden in our coloring.

With the decoupled Hamiltonian in (110), we can then learn the coefficients in
each one- or two-mode cluster independently and in parallel using the
algorithms described in Sections B and C. Looping over all colors
$c\in\\{1,2,\cdots,\chi\\}$, we will obtain all the coefficients in the
Hamiltonian.

##  Appendix E Deviation from the effective dynamics

In this section, we consider the error introduced by simulating the effective
dynamics with the insertion of random unitaries, as mentioned in Section C.1.
Suppose $\mathcal{D}$ is a distribution over the set of unitaries, and the
initial state of the system is represented by the density matrix $\rho(0)$.
The actual final state obtained after the insertion of $r$ random unitaries is

$\mathbb{E}_{U_{j}\sim\mathcal{D}}\left(\prod_{1\leq j\leq
r}^{\leftarrow}U_{j}^{\dagger}e^{-i\tau
H}U_{j}\right)\rho(0)\left(\prod_{1\leq j\leq
r}^{\rightarrow}U_{j}^{\dagger}e^{i\tau H}U_{j}\right),$ (111)

where each $U_{j}$ is inserted after time $\tau=\frac{t}{r}$. On the other
hand, the desired final state, which facilitates the subsequent steps of the
learning process, is

$e^{-itH_{\mathrm{effective}}}\rho(0)e^{itH_{\mathrm{effective}}},$ (112)

where $H_{\mathrm{effective}}$ is the effective Hamiltonian:

$H_{\mathrm{effective}}=\mathbb{E}_{U\sim\mathcal{D}}U^{\dagger}HU.$

In this section, we provide an analysis of the difference between the two
dynamics for a certain class of Hamiltonians and thereby complete the analysis
of approximation errors investigated in Section B. For the sake of the
Hamiltonians studied in this paper, we consider the Hamiltonians of the
following form:

$H=\sum_{\braket{i,j}}h_{ij}b_{i}^{\dagger}b_{j}+\sum_{i}\omega_{i}n_{i}+\frac{1}{2}\sum_{\braket{jklm}}\xi_{jklm}b_{j}^{\dagger}b_{k}^{\dagger}b_{l}b_{m},$
(113)

where in the last term we denote by $\braket{jklm}$ the index quadruples such
that $\\{j,k,l,m\\}$ form a connected subgraph in the underlying graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ of bosonic modes. We begin with a
lemma describing the action of these Hamiltonians on the product of coherent
states.

######  Lemma 6.

Let

$H=\sum_{\braket{i,j}}h_{ij}b_{i}^{\dagger}b_{j}+\sum_{i}\omega_{i}n_{i}+\frac{1}{2}\sum_{\braket{jklm}}\xi_{jklm}b_{j}^{\dagger}b_{k}^{\dagger}b_{l}b_{m},$
(114)

and

$\phi=\bigotimes_{i\in\mathcal{V}}\left(e^{-|\alpha_{i}|^{2}/2}\sum_{k=0}^{\infty}\frac{\alpha_{i}^{k}e^{-i\zeta_{i,k}}}{\sqrt{k!}}\ket{k}_{i}\right),$
(115)

where $\alpha_{i}$ is a complex number of magnitude $O(1)$, and
$\zeta_{i,k}\in\mathbb{R}$ can be any real number. Then

$\|H\phi\|=O(N\max\\{|\xi_{jklm}|,|\omega_{i}|,|h_{i,j}|\\}),$ (116)

and

$\|H^{2}\phi\|=O(N^{2}(\max\\{|\xi_{jklm}|,|\omega_{i}|,|h_{i,j}|\\})^{2}),$
(117)

where $N=|\mathcal{V}|+|\mathcal{E}|$.

###### Proof.

It suffices to prove the result for
$H/\max\\{|\xi_{jklm}|,|\omega_{i}|,|h_{i,j}|\\}$. Therefore we assume
$\max\\{|\xi_{jklm}|,|\omega_{i}|,|h_{i,j}|\\}=1$ without loss of generality.
Notice that $H$ is the sum of $O(N)$ terms, and each term takes the form
$b^{\dagger}_{p}b_{q}$ or $b^{\dagger}_{p}b_{q}b^{\dagger}_{r}b_{s}$, where
$p$, $q$, $r$, $s$ may be repeated. We will prove that each term acting on
$\ket{\phi}$ yields a state whose norm is $O(1)$. We first demonstrate this
for $\left\|b^{\dagger}_{p}b_{q}b^{\dagger}_{r}b_{s}\phi\right\|$. Simple
calculation shows

$\displaystyle b^{\dagger}_{p}b_{q}b^{\dagger}_{r}b_{s}\phi$ (118)
$\displaystyle=$
$\displaystyle\bigotimes_{i\notin\\{p,q,r,s\\}}\left(e^{-|\alpha_{i}|^{2}/2}\sum_{k=0}^{\infty}\frac{\alpha_{i}^{k}e^{-i\zeta_{i,k}}}{\sqrt{k!}}\ket{k}_{i}\right)\otimes\bigotimes_{j\in\\{p,q,r,s\\}}\left(e^{-|\alpha_{j}|^{2}/2}\sum_{k=0}^{\infty}\frac{\alpha_{j}^{k}e^{-i\zeta_{j,k}}}{\sqrt{k!}}\sqrt{P_{j}(k)}\ket{k+\sigma_{j}}_{j}\right),$

where $P_{j}$’s are polynomials with $\sum_{j\in\\{p,q,r,s\\}}\deg P_{j}=4$,
and $\sigma_{j}$ is an integer determined by the numbers of $b^{\dagger}_{j}$
and $b_{j}$ in $b^{\dagger}_{p}b_{q}b^{\dagger}_{r}b_{s}$. For example, if
$p=q=r=1$, $s=2$, then $P_{1}(k)=(k+1)^{3}$, $\sigma_{1}=1$, $P_{2}(k)=k$,
$\sigma_{2}=-1$. Straight calculations can show that

$\displaystyle\left\|e^{-|\alpha_{j}|^{2}/2}\sum_{k=0}^{\infty}\frac{\alpha_{j}^{k}e^{-i\zeta_{j,k}}}{\sqrt{k!}}\sqrt{P_{j}(k)}\ket{k+\sigma_{j}}_{j}\right\|^{2}$
(119) $\displaystyle={}$ $\displaystyle
e^{-|\alpha_{j}|^{2}}\sum_{k=0}^{\infty}\frac{|\alpha_{j}|^{2k}}{k!}P_{j}(k)=Q_{j}(|\alpha_{j}|^{2})=O(1).$

where $Q_{j}$ is a polynomial that can be determined by $P_{j}$, but we do not
care about its explicit form. Therefore we have shown that

$\left\|b^{\dagger}_{p}b_{q}b^{\dagger}_{r}b_{s}\phi\right\|=\sqrt{\prod_{j\in\\{p,q,r,s\\}}Q_{j}(|\alpha_{j}|^{2})}=O(1).$
(120)

Similarly, we can show that $\|b^{\dagger}_{p}b_{q}\phi\|=O(1)$. Therefore
(116) is established.

Next, we will prove (117). We can fully expand $H^{2}$ into $O(N^{2})$ terms,
each of which has the form
$b^{\dagger}_{p}b_{q}b^{\dagger}_{p^{\prime}}b_{q^{\prime}}$,
$b^{\dagger}_{p}b_{q}b^{\dagger}_{r}b_{s}b^{\dagger}_{p^{\prime}}b_{q^{\prime}}$,
$b^{\dagger}_{p}b_{q}b^{\dagger}_{p^{\prime}}b_{q^{\prime}}b^{\dagger}_{r^{\prime}}b_{s^{\prime}}$,
or
$b^{\dagger}_{p}b_{q}b^{\dagger}_{r}b_{s}b^{\dagger}_{p^{\prime}}b_{q^{\prime}}b^{\dagger}_{r^{\prime}}b_{s^{\prime}}$.
Again, we may go through a similar process as above and conclude that each
term acting on $\phi$ yields a state of magnitude $O(1)$. ∎

Assume that $\ket{\phi_{0}}=\bigotimes_{i}\ket{\alpha_{i}}$ is a product of
coherent states, and $\ket{\phi_{t}}$ is the state obtained by evolving under
the effective dynamics ${H_{\mathrm{effective}}^{(1)}}$ for time $t$, i.e.,
$\ket{\phi_{t}}=e^{-it{H_{\mathrm{effective}}^{(1)}}}\ket{\phi_{0}}$, then
$\ket{\phi_{t}}$ is a state of the form described in (115) for the
distribution $\mathcal{D}$ used in previous sections. Using density matrices,
the effective dynamics with the Hamiltonian ${H_{\mathrm{effective}}^{(1)}}$
starts from the state $\rho(0):=\ket{\phi_{0}}\bra{\phi_{0}}$ and end up in
the state $\rho(t):=\ket{\phi_{t}}\bra{\phi_{t}}$ at time $t$, while the
actual final state obtained is given by (111). To bound its distance from the
desired state $\rho(t)$, we define the following density operators:

$\rho^{(\ell)}(t)=\mathbb{E}\left(\prod_{1\leq
j\leq\ell}^{\leftarrow}U_{j}^{\dagger}e^{-i\tau
H}U_{j}\right)\rho(t-\ell\tau)\left(\prod_{1\leq
j\leq\ell}^{\rightarrow}U_{j}^{\dagger}e^{i\tau H}U_{j}\right).$ (121)

Then $\rho^{(0)}(t)=\rho(t)$ and $\rho^{(r)}(t)$ is the density operator in
(111). Now consider the distance between $\rho^{(L-1)}(t)$ and
$\rho^{(L)}(t)$. Define

$Q^{(L)}=\prod_{1\leq j\leq L-1}^{\rightarrow}U_{j}^{\dagger}e^{i\tau
H}U_{j},$ (122)

then by the independence of $U_{j}$, we have

$\displaystyle\|\rho^{(L)}(t)-\rho^{(L-1)}(t)\|_{*}$ (123)
$\displaystyle=\left\|\mathbb{E}_{Q^{(L)}}\left[Q^{(L)}\left(\mathbb{E}_{U}\left(U^{\dagger}e^{-i\tau
H}U\rho(t-L\tau)U^{\dagger}e^{i\tau
H}U-e^{-i\tau{H_{\mathrm{effective}}^{(1)}}}\rho(t-L\tau)e^{-i\tau{H_{\mathrm{effective}}^{(1)}}}\right)\right)(Q^{(L)})^{\dagger}\right]\right\|_{*}$
$\displaystyle\leq\mathbb{E}_{Q^{(L)}}\left\|Q^{(L)}\left(\mathbb{E}_{U}\left(U^{\dagger}e^{-i\tau
H}U\rho(t-L\tau)U^{\dagger}e^{i\tau
H}U-e^{-i\tau{H_{\mathrm{effective}}^{(1)}}}\rho(t-L\tau)e^{-i\tau{H_{\mathrm{effective}}^{(1)}}}\right)\right)(Q^{(L)})^{\dagger}\right\|_{*}$
$\displaystyle=\mathbb{E}_{Q^{(L)}}\left\|\mathbb{E}_{U}\left(U^{\dagger}e^{-i\tau
H}U\rho(t-L\tau)U^{\dagger}e^{i\tau
H}U-e^{-i\tau{H_{\mathrm{effective}}^{(1)}}}\rho(t-L\tau)e^{-i\tau{H_{\mathrm{effective}}^{(1)}}}\right)\right\|_{*}$
$\displaystyle=\left\|\mathbb{E}_{U}\left(U^{\dagger}e^{-i\tau
H}U\rho(t-L\tau)U^{\dagger}e^{i\tau
H}U-e^{-i\tau{H_{\mathrm{effective}}^{(1)}}}\rho(t-L\tau)e^{-i\tau{H_{\mathrm{effective}}^{(1)}}}\right)\right\|_{*},$

where$\|\cdot\|_{*}$ denotes the trace norm (nuclear norm). The fourth line
follows from the property of trace norm and the fact that $Q^{(L)}$ is
unitary. From the Taylor expansion, one can obtain

$\displaystyle\mathbb{E}\left(U^{\dagger}e^{-i\tau
H}U\rho(t-L\tau)U^{\dagger}e^{i\tau H}U\right)-\rho(t-L\tau)$ (124)
$\displaystyle=\mathbb{E}\left(e^{-i\tau U^{\dagger}HU}\rho(t-L\tau)e^{i\tau
U^{\dagger}HU}\right)-\rho(t-L\tau)$
$\displaystyle=\mathbb{E}\left(-i\tau[U^{\dagger}HU,\rho(t-L\tau)]-\int_{0}^{\tau}e^{-isU^{\dagger}HU}[U^{\dagger}HU,[U^{\dagger}HU,\rho(t-L\tau)]]e^{isU^{\dagger}HU}(\tau-s)\right)$
$\displaystyle=-i\tau[\mathbb{E}(U^{\dagger}HU),\rho(t-L\tau)]-\mathbb{E}\left(\int_{0}^{\tau}e^{-isU^{\dagger}HU}[U^{\dagger}HU,[U^{\dagger}HU,\rho(t-L\tau)]]e^{isU^{\dagger}HU}(\tau-s)\right)$
$\displaystyle=-i\tau[{H_{\mathrm{effective}}^{(1)}},\rho(t-L\tau)]-\mathbb{E}\left(\int_{0}^{\tau}e^{-isU^{\dagger}HU}[U^{\dagger}HU,[U^{\dagger}HU,\rho(t-L\tau)]]e^{isU^{\dagger}HU}(\tau-s)\right).$

Similarly, one has

$\displaystyle\mathbb{E}\left(e^{-i\tau{H_{\mathrm{effective}}^{(1)}}}\rho(t-L\tau)e^{i\tau{H_{\mathrm{effective}}^{(1)}}}\right)-\rho(t-L\tau)$
(125)
$\displaystyle=-i\tau[{H_{\mathrm{effective}}^{(1)}},\rho(t-L\tau)]-\int_{0}^{\tau}e^{-is{H_{\mathrm{effective}}^{(1)}}}[{H_{\mathrm{effective}}^{(1)}},[{H_{\mathrm{effective}}^{(1)}},\rho(t-L\tau)]]e^{is{H_{\mathrm{effective}}^{(1)}}}(\tau-s).$

Combining (124) and (125), one obtains

$\displaystyle\left\|\mathbb{E}\left(U^{\dagger}e^{-i\tau
H}U\rho(t-L\tau)U^{\dagger}e^{i\tau
H}U-e^{-i\tau{H_{\mathrm{effective}}^{(1)}}}\rho(t-L\tau)e^{-i\tau{H_{\mathrm{effective}}^{(1)}}}\right)\right\|_{*}$
(126)
$\displaystyle\leq\left\|\mathbb{E}\left(\int_{0}^{\tau}e^{-isU^{\dagger}HU}[U^{\dagger}HU,[U^{\dagger}HU,\rho(t-L\tau)]]e^{isU^{\dagger}HU}(\tau-s)\right)\right\|_{*}$
$\displaystyle+\left\|\int_{0}^{\tau}e^{-is{H_{\mathrm{effective}}^{(1)}}}[{H_{\mathrm{effective}}^{(1)}},[{H_{\mathrm{effective}}^{(1)}},\rho(t-L\tau)]]e^{is{H_{\mathrm{effective}}^{(1)}}}(\tau-s)\right\|_{*}$
$\displaystyle\leq\frac{1}{2}\tau^{2}\bigg{(}\sup_{U}\left\|[U^{\dagger}HU,[U^{\dagger}HU,\rho(t-L\tau)]]\right\|_{*}+\left\|[{H_{\mathrm{effective}}^{(1)}},[{H_{\mathrm{effective}}^{(1)}},\rho(t-L\tau)]]\right\|_{*}\bigg{)}$

One only needs to bound
$\left\|[U^{\dagger}HU,[U^{\dagger}HU,\rho(t-L\tau)]]\right\|_{*}$ and
$\left\|[{H_{\mathrm{effective}}^{(1)}},[{H_{\mathrm{effective}}^{(1)}},\rho(t-L\tau)]]\right\|_{*}$.
By a direct calculation, one sees that

$\displaystyle\left\|[{H_{\mathrm{effective}}^{(1)}},[{H_{\mathrm{effective}}^{(1)}},\rho(t-L\tau)]]\right\|_{*}$
(127)
$\displaystyle\leq\|{H_{\mathrm{effective}}^{(1)}}^{2}\rho(t-L\tau)\|_{*}+2\|{H_{\mathrm{effective}}^{(1)}}\rho(t-L\tau){H_{\mathrm{effective}}^{(1)}}\|_{*}+\|\rho(t-L\tau){H_{\mathrm{effective}}^{(1)}}^{2}\|_{*}$
$\displaystyle=2\|\phi_{t-L\tau}\|\|{H_{\mathrm{effective}}^{(1)}}^{2}\phi_{t-L\tau}\|+2\|{H_{\mathrm{effective}}^{(1)}}\phi_{t-L\tau}\|^{2}$
$\displaystyle\leq CN^{2}\max\\{|\xi_{jklm}|,|\omega_{i}|,|h_{i,j}|\\}^{2},$

where $C=\mathcal{O}(1)$ is a constant, and we have used the property of the
trace norm for rank-$1$ matrices. In the last step, we are using 6 with
$H={H_{\mathrm{effective}}^{(1)}}$ and $\phi=\phi_{t-L\tau}$. Similarly, one
can obtain

$\displaystyle\left\|[U^{\dagger}HU,[U^{\dagger}HU,\rho(t-L\tau)]]\right\|_{*}$
(128)
$\displaystyle=2\|\phi_{t-L\tau}\|\|U^{\dagger}H^{2}U\phi_{t-L\tau}\|+2\|U^{\dagger}HU\phi_{t-L\tau}\|^{2}$
$\displaystyle=2\|H^{2}U\phi_{t-L\tau}\|+2\|HU\phi_{t-L\tau}\|^{2}$
$\displaystyle\leq CN^{2}\max\\{|\xi_{jklm}|,|\omega_{i}|,|h_{i,j}|\\}^{2}.$

In the last step, we are using 6 with $H=H$ and $\phi=U\phi_{t-L\tau}$. As a
result, we have proved the following:

######  Theorem 7.

For a Hamiltonian of the form described in (113) and a product of coherent
states $\ket{\phi_{0}}=\otimes_{i}\ket{\alpha_{i}}$ such that $\alpha_{i}$ are
$\mathcal{O}(1)$ constants, we have

$\displaystyle\Bigg{\|}\mathbb{E}_{U_{j}\sim\mathcal{D}}\left(\prod_{1\leq
j\leq r}^{\leftarrow}U_{j}^{\dagger}e^{-i\tau
H}U_{j}\right)\rho(0)\left(\prod_{1\leq j\leq
r}^{\rightarrow}U_{j}^{\dagger}e^{i\tau
H}U_{j}\right)-e^{-itH_{\mathrm{effective}}}\rho(0)e^{itH_{\mathrm{effective}}}\Bigg{\|}_{*}$
(129) $\displaystyle\leq
CN^{2}\frac{t^{2}}{r}\max\\{|\xi_{jklm}|,|\omega_{i}|,|h_{i,j}|\\}^{2},$

where $\rho(0)=\ket{\phi_{0}}\bra{\phi_{0}}$,
$H_{\mathrm{effective}}=\mathbb{E}_{U\sim\mathcal{D}}U^{\dagger}HU$, $C$ is a
$\mathcal{O}(1)$ constant, $N=|\mathcal{V}|+|\mathcal{E}|$ and
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ is the underlying graph of bosonic
modes.

###### Proof.

The left-hand side of (129) can be expressed by
$\|\rho^{(r)}(t)-\rho^{(0)}(t)\|_{*}$, where $\rho^{(r)}(t)$ and
$\rho^{(0)}(t)$ are defined in (121). Thus

$\displaystyle\|\rho^{(r)}(t)-\rho^{(0)}(t)\|_{*}\leq\sum_{L=1}^{r}\|\rho^{(L)}(t)-\rho^{(L-1)}(t)\|_{*}$
(130)
$\displaystyle\leq\sum_{L=1}^{r}CN^{2}\tau^{2}\max\\{|\xi_{jklm}|,|\omega_{i}|,|h_{i,j}|\\}^{2}=CN^{2}\frac{t^{2}}{r}\max\\{|\xi_{jklm}|,|\omega_{i}|,|h_{i,j}|\\}^{2},$

where we have used (126), (127) and (128) in the second inequality. ∎