††thanks: Now at Oxford Quantum Circuits, Thames Valley Science Park, Reading
RG2 9LH<EMAIL_ADDRESS>
# Characterisation of spatial charge sensitivity in a multi-mode
superconducting qubit

J. Wills  G. Campanaro  S. Cao  S. D. Fasciati  P. J. Leek  B. Vlastakis
Clarendon Laboratory, Department of Physics, University of Oxford, Oxford OX1
3PU, United Kingdom

###### Abstract

Understanding and suppressing sources of decoherence is a leading challenge in
building practical quantum computers. In superconducting qubits, low frequency
charge noise is a well-known decoherence mechanism that is effectively
suppressed in the transmon qubit. Devices with multiple charge-sensitive modes
can exhibit more complex behaviours, which can be exploited to study charge
fluctuations in superconducting qubits. Here we characterise charge-
sensitivity in a superconducting qubit with two transmon-like modes, each of
which is sensitive to multiple charge-parity configurations and charge-offset
biases. Using Ramsey interferometry, we observe sensitivity to four charge-
parity configurations and track two independent charge-offset drifts over hour
timescales. We provide a predictive theory for charge sensitivity in such
multi-mode qubits which agrees with our results. Finally, we demonstrate the
utility of a multi-mode qubit as a charge detector by spatially tracking
local-charge drift.

##  I INTRODUCTION

Realising practical quantum computation requires the low-error operation of
many fully controlled quantum bits with individual state readout and
initialisation DiVincenzo (2000). Superconducting circuits are well
established as a potential platform for quantum computation Devoret and
Schoelkopf (2013); Arute  _et al._ (2019). One superconducting qubit variant,
the transmon Koch  _et al._ (2007), has found enduring success due to its
resilience to decoherence and simplicity in design. However the development of
alternatives is an active area of research Manucharyan  _et al._ (2009); Kou
_et al._ (2017); Groszkowski  _et al._ (2018); Gyenis  _et al._ (2021). One
particular issue for the transmon is its fixed and sizeable coupling to other
qubits and circuit modes due to its large dipole moment, which can cause
unwanted crosstalk and non-negligible computational errors Sundaresan  _et
al._ (2020). This issue can be addressed by using tuneable couplers Roy  _et
al._ (2017); Srinivasan  _et al._ (2011), level-structure engineering Richer
_et al._ (2017); Noguchi  _et al._ (2020), or auxiliary circuit modes Smith
_et al._ (2020).

Here we investigate the multi-mode transmon qubit. Such qubits allow for
strongly coupled, and far detuned modes, that can assist in preventing cross-
talk Finck  _et al._ (2021) and limit contributions to decay due to the
Purcell effect Gambetta  _et al._ (2011). The energy level structure can be
engineered to allow for generous resonance conditions and fast entangling
operations Roy  _et al._ (2018). In addition, this can create mechanisms for
photon-dephasing suppression Richer  _et al._ (2017); Zhang  _et al._ (2017)
or tuneable coupling Gambetta  _et al._ (2011), making multi-mode qubits of
potential use in superconducting quantum processors.

Since the modes of these multi-mode systems are transmon-like, they are
susceptible to the same decoherence mechanisms, and can exhibit more complex
sensitivities to noise Pashkin  _et al._ (2003). One such decoherence
mechanism that affects superconducting circuits is low frequency, $1/f$-like
charge noise Gustafsson  _et al._ (2013). This has been well characterised and
studied in single-mode transmons, which are typically designed to operate in a
charge-insensitive regime Schreier  _et al._ (2008). Alternatively, devices
can be deliberately designed to be sensitive to charge noise, and operate as
detectors of it, in order to better understand and characterise this source of
decoherence Christensen  _et al._ (2019); Wilen  _et al._ (2021); Tennant  _et
al._ (2021). The origins of charge noise in superconducting circuits are
theorised to be surface drifts, patch potentials, voltage fluctuations in
control electronics, or the effects of the absorption of cosmic radiation
within substrates Christensen  _et al._ (2019); Wilen  _et al._ (2021);
Martinis (2021); Tennant  _et al._ (2021).

In this paper, we investigate the effect of charge noise in a two-mode
transmon design, depicted in Fig. 1. We present a description of the mode
structure and device operation, as well as a predictive theory for charge
sensitivity in such devices. We observe sensitivity to the four charge-parity
configurations that arise as a result of the two degrees of freedom in a
charge-sensitive two-mode transmon qubit. We also demonstrate effective
suppression of charge sensitivity in an alternative design two-mode device
with high $E_{J}/E_{C}$ ratio, in accordance with our predictive theory.
Finally, using Ramsey interferometry we track charge-offset fluctuations
showing a proof-of-concept detector for spatial drifts in charge noise over
$100~{}\rm{\mu m}$ length scales.

Understanding decoherence mechanisms, such as charge noise, is vital to the
application of superconducting qubits in a low-error quantum processor. In
turn, a charge-sensitive multi-mode transmon could prove a useful tool for
identifying the origins of such charge fluctuations in high-coherence quantum
devices.

Figure 1: Device Description. (a) A simplified schematic of the device shows a
coaxial pad geometry with a split inner conductor resulting in three
superconducting islands. Two near-identical Josephson junctions with phase
$\varphi_{1}$, $\varphi_{2}$ are used resulting in two non-degenerate modes
($\Delta,\Sigma$). A lumped element LC resonator (blue) on the opposing side
of substrate allows for dispersive readout of modes. (b) An equivalent circuit
of the two-mode coaxial transmon. The two inner superconducting islands are
connected by a coupling capacitance ($C_{m}$), and connected to the outer
superconducting island by Josephson junction ($E_{J}$) and capacitance ($2C$).
(c) Energy level diagram of the two-mode coaxial transmon illustrating
anharmonicities ($\eta_{\Sigma,\Delta}$), and state-dependent shift
($\chi_{\Sigma\Delta}$). States labelled as $\ket{nm}$ where $n(m)$
corresponds to the number of excitations in the $\Sigma$-mode ($\Delta$-mode).
The V-shaped qutrit subspace is highlighted in grey. (d) (e) Simulation of
charge bias induced by a point charge (1e) in the summation ($n_{g\Sigma}$,
left), and difference ($n_{g\Delta}$, right) charge configurations on the two
inner islands.

##  II A TWO-MODE COAXIAL TRANSMON

Figure 2: Charge dispersion and parity configurations in a multimode system.
(a) The
$V(\varphi_{\Sigma},\varphi_{\Delta})=-2E_{J}\cos{\frac{\varphi_{\Sigma}}{2}}\cos{\frac{\varphi_{\Delta}}{2}}$
periodic potential energy landscape of the device, where local minima are
$2\pi$ periodic along the diagonals. White dashed arrows indicate the four
possible tunnelling routes to nearest neighbouring lattice sites. The
probability of tunnelling to a neighbouring well in each of the four
directions is equal. (b) Dispersion of eigenenergies of the four lowest energy
levels of the qubit Hamiltonian as a function of gate charge offset $n_{g1}$,
with constant gate charge offset $n_{g2}$, showing the four possible parity
configurations (OO, EO, OE, EE). (c) Comparison of numerical (solid) and tight
binding model (dashed) calculations of charge dispersion for the lowest four
energy levels of the qubit Hamiltonian, as a function of the ratio
$E_{J}/E_{C}$. The right vertical scale gives the charge dispersion in MHz for
a lowest transition frequency of 5 GHz. Vertical dotted (dashed) line shows
the regime of the charge sensitive Device A, $E_{J}/E_{C}$ = 22 (charge
suppressed Device B, $E_{J}/E_{C}$ = 70).

A transmon qubit is a simple superconducting circuit consisting of a
capacitance in parallel with a Josephson junction, and is insensitive to
charge fluctuations across the capacitor in the regime of large Josephson-to-
charging energy ratio, $E_{J}/E_{C}$ Koch  _et al._ (2007). Here we work with
a circuit with two transmon-like modes, built from three superconducting
islands and two Josephson junctions, depicted in Fig. 1(b).

The two-mode circuit has the Hamiltonian:

$\displaystyle\hat{H}=$ $\displaystyle
4E_{C}(\hat{n}_{1}-n_{g1})^{2}+4E_{C}(\hat{n}_{2}-n_{g2})^{2}$ (1)
$\displaystyle+4E_{p}(\hat{n}_{1}-n_{g1})(\hat{n}_{2}-n_{g2})$ $\displaystyle-
E_{J}\cos{\hat{\varphi}_{1}}-E_{J}\cos{\hat{\varphi}_{2}}$

where the charging energy $E_{C}=e^{2}C^{*}/2(C^{*2}-C_{m}^{2})$, with
$C^{*}=C(\frac{C+2C_{m}}{C+C_{m}})$ the total shunting capacitance of each
island. $E_{J}$ is the Josephson energy of each junction, and
$E_{p}=e^{2}C_{m}/({C^{*}}^{2}-C_{m}^{2})$ is the coupling energy between the
two modes Pashkin  _et al._ (2003). The Hamiltonian of this system is
identical to that of two resonantly-coupled transmons, the eigenmodes of which
are sum and difference modes with
$\hat{\varphi}_{\Sigma}=\hat{\varphi}_{1}+\hat{\varphi}_{2}$ and
$\hat{\varphi}_{\Delta}=\hat{\varphi}_{1}-\hat{\varphi}_{2}$ respectively. In
this basis, the Hamiltonian of the circuit becomes:

$\displaystyle\hat{H}=$ $\displaystyle
8E_{C_{\Sigma}}(\hat{n}_{\Sigma}-n_{g\Sigma})^{2}+8E_{C_{\Delta}}(\hat{n}_{\Delta}-n_{g\Delta})^{2}$
(2)
$\displaystyle-2E_{J}\cos{\frac{\hat{\varphi}_{\Sigma}}{2}}\cos{\frac{\hat{\varphi}_{\Delta}}{2}}$

where $E_{C_{\Sigma}}=E_{C}+E_{p}/2=e^{2}/2(C^{*}-C_{m})$ and
$E_{C_{\Delta}}=E_{C}-E_{p}/2=e^{2}/2(C^{*}+C_{m})$ are the charging energies
of the sum and difference modes respectively.

In this representation, it is clearer that the first two excited states of
this system correspond to an excitation in either the in-phase
($\Sigma$-mode), or out-of-phase oscillations ($\Delta$-mode) of the two inner
superconducting islands, shown in Fig. 1 (a). The antisymmetric mode is lower
in frequency, has an electric dipole moment and couples well to electric
fields polarised in the plane of the device. The symmetric mode is higher in
frequency, has an electric quadrupole moment and so couples well to radial
fields and coaxial control ports. This difference in polarisation symmetry can
be used for Purcell protection Gambetta  _et al._ (2011).

Starting from the Hamiltonian of Eqn. 2, we can perform a second quantisation
by expanding the cosine potential terms to fourth order, and assuming a weakly
non-linear oscillator-like behaviour. Keeping only counter-rotating terms
produces:

$\displaystyle\hat{H}/\hbar$
$\displaystyle=\omega_{\Sigma}\hat{a}^{\dagger}_{\Sigma}\hat{a}_{\Sigma}+\omega_{\Delta}\hat{a}^{\dagger}_{\Delta}\hat{a}_{\Delta}$
(3)
$\displaystyle-\eta_{\Sigma}/2*\hat{a}^{\dagger}_{\Sigma}{}^{2}\hat{a}_{\Sigma}{}^{2}-\eta_{\Delta}/2*\hat{a}^{\dagger}_{\Delta}{}^{2}\hat{a}_{\Delta}{}^{2}$
$\displaystyle-\chi_{\Sigma\Delta}*\hat{a}^{\dagger}_{\Sigma}\hat{a}_{\Sigma}\hat{a}^{\dagger}_{\Delta}\hat{a}_{\Delta}$

where $\hat{a}_{\Sigma(\Delta)}^{(\dagger)}$ is the annihilation (creation)
operator, $\omega_{\Sigma(\Delta)}$ is the frequency and
$\eta_{\Sigma(\Delta)}$ is the anharmonicity of the $\Sigma$-mode
($\Delta$-mode). $\chi_{\Sigma\Delta}$ is the non-linear coupling between the
modes. This Hamiltonian is derived in Appendix A.

Both modes are very strongly coupled to each other through the junctions of
the device. Notably, the non-linear coupling ($\chi_{\Sigma\Delta}$) is larger
than their respective anharmonicities ($\eta_{\Sigma}$, $\eta_{\Delta}$). This
produces addressable transitions that make up an effective V-shaped qutrit
energy diagram Dumur  _et al._ (2015), as shown in Fig. 1 (c). The states are
labelled as $\ket{nm}$, where $n(m)$ corresponds to the number of excitations
in the $\Sigma$-mode ($\Delta$-mode). We have the potential to use the qutrit
system for all-microwave two-qubit gates Hazra  _et al._ (2020); Finck  _et
al._ (2021), using one transition as a computational bit and the other to
generate entanglement with other qutrits. This allows for far-detuned
computational transitions for minimised single-qubit-gate crosstalk, while
retaining gate speeds comparable to fixed coupled transmons. These features
make the two-mode transmon a potentially useful component for an extensible
quantum computing architecture.

We build this system in a coaxial geometry with out-of-plane coupling to a
lumped-element resonator for dispersive readout Rahamim  _et al._ (2017). The
designed circuit has a relative coupling $E_{p}/E_{C}$ $=2.5$, which is
approximately two orders of magnitude larger than typical qubit-qubit
couplings.

The Hamiltonian of Eqn. 2 shows a dependence on gate-charge offsets
$n_{g\Sigma}$ and $n_{g\Delta}$, corresponding to the sum and difference of
gate-charge offsets of the two inner islands, $n_{g1}$ and $n_{g2}$. Using
electrostatic simulations, we find the induced spatially dependant offset
charge due to a point charge on the surface of the substrate, shown in Fig. 1
(d) and (e). The sensitivity pattern differs for $n_{g\Sigma}$ and
$n_{g\Delta}$ due to the symmetric/anti-symmetric behaviour of each
configuration. This difference crucially allows us to spatially detect local
charge fluctuations.

We can understand the nature and behaviour of charge dispersion in Eqn. 2 in a
transmon regime ($E_{C},E_{p}\ll E_{J}$) by using a 2D tight-binding
approximation (see Appendix D). In Fig. 2 (a), the potential energy landscape
shows that for each energy minimum, there are four neighbouring lattice sites
dependant on $n_{g1}$ and $n_{g2}$, with an equal tunnel-barrier energy
between them. This will produce a dispersion relationship of the form
$E(n_{g1},n_{g2})\sim(\cos{n_{g1}}+\cos{n_{g2}})$.

This can be represented in terms of the sum and difference offset charges,
$n_{g\Sigma}$ and $n_{g\Delta}$, as:

$E_{mn}(n_{g\Sigma},n_{g\Delta})\approx\overline{E_{mn}}+\frac{\epsilon_{mn}}{4}\cos{\pi
n_{g\Sigma}}\cos{\pi n_{g\Delta}}$ (4)

where $\epsilon_{mn}$ is the maximum measured charge dispersion for the
$E_{mn}$ level, where $m$($n$) is the number of excitations in the
$\Sigma$-mode ($\Delta$-mode).

One source of decoherence is sudden changes in offset charge due to tunnelling
of quasiparticles across Josephson junctions from one superconducting island
to another Serniak  _et al._ (2018). This corresponds to jumps in either
$n_{g1}\rightarrow n_{g1}+0.5$, $n_{g2}\rightarrow n_{g2}+0.5$, otherwise
denoted as jumps in charge parity Odd (O) to Even (E) or Even to Odd. This
results in four different parity configurations, two for each mode in the
system. Fig. 2 (b) shows the energy-level dispersion with dependence on four
parity configurations (OO, EO, OE, EE), with a maximum dispersion of
$\epsilon_{mn}$.

In Fig. 2 (c) we show numerical calculation of charge dispersion
$\epsilon_{mn}$ as a function of $E_{J}$ and $E_{C}$ for $E_{p}=0.4E_{C}$,
obtained by calculating the eigenvalues of the Hamiltonian in Eqn. 1 in the
charge basis. We use a semi-analytical wavefunction approach Catelani  _et
al._ (2011) to derive an analytical form of this dependence:

$\displaystyle\epsilon_{mn}$ $\displaystyle\approx
A_{0}E_{J}\frac{2^{2(m+n)}}{m!n!}$ (5)
$\displaystyle\times\left(\frac{E_{J}}{E_{C}(1+E_{p}/2)}\right)^{m/2}\left(\frac{E_{J}}{E_{C}(1-E_{p}/2)}\right)^{n/2}$
$\displaystyle\times\exp{-\left(\sqrt{\frac{2E_{J}}{E_{C}(1+E_{p}/2)}}+\sqrt{\frac{2E_{J}}{E_{C}(1-E_{p}/2)}}\,\right)}$

This charge dispersion follows an exponential suppression, and in the limit
$E_{p}\to 0$, the exponent term tends towards the value for the standard
transmon. The factor $A_{0}$ is obtained empirically by fitting to the
numerical calculation of charge dispersion, shown in Fig. 2 (c).

Each mode is sensitive to both charge offsets, i.e. if there are fluctuations
in $n_{g\Sigma}$, this can be observed in dispersion of the $\Delta$-mode.
Importantly, this means that a single mode can detect fluctuations in both
offset charges simultaneously.

##  III METHODS

In this report, we present two devices. One device (A), is designed for charge
sensitivity, and the other (device B) is designed for charge-noise
suppression. We fabricate device A (B) through electron beam lithography (and
photolithography), patterning both sides of a $0.5~{}{\rm mm}$ sapphire
(silicon) substrate. Each device is mounted inside an aluminium sample holder
within a mu-metal magnetic shield, anchored to the 10 mK stage of a dilution
refrigerator, operating with a standard cQED experimental setup Spring  _et
al._ (2021). A comparison of the device designs is shown in Appendix B.

A state-dependent resonator frequency shift $\chi_{\Sigma(\Delta)r}$ exists
for each mode, allowing for simultaneous dispersive readout of the multi-mode
state of the device, as defined in Appendix A. Using qubit spectroscopy, we
find the transition frequencies for device A to be $\omega_{\Delta}/2\pi=5.51$
GHz, and $\omega_{\Sigma}/2\pi=6.71$ GHz, with anharmonicites of
$\eta_{\Delta}/2\pi=-0.38$ GHz and $\eta_{\Sigma}/2\pi=-0.34$ GHz. We identify
the $\ket{00}-\ket{11}$ transition at
$(\omega_{\Delta}+\omega_{\Sigma}-\chi_{\Sigma\Delta})/4\pi=5.86$ GHz, showing
an inter-modal state-dependent shift of $\chi_{\Sigma\Delta}/2\pi=0.50$ GHz.
These values are consistent with numerical solutions obtained with finite-
element simulation methods and energy-participation-ratio (EPR) analysis Minev
_et al._ (2021). A summary of measured parameters for both devices is shown in
Table 1, in Appendix B.

From these parameters, we use numerical methods to estimate values of
$E_{J}/h=11$ GHz, $E_{C}/h=0.5$ GHz, and $E_{p}/h=0.2$ GHz, as defined in
section II. Given these parameters, we estimate a charge dispersion of the
lowest two modes to be 4 MHz and 4.1 MHz, calculated numerically. This
predicted charge dispersion and $E_{J}/E_{C}$ regime is shown in Fig. 2 (c)
for both devices.

We next perform time-domain measurements of the energy relaxation time $T_{1}$
and spin-echo coherence time $T_{2E}$ of the two modes of the device. We
perform 100 repeated measurements over the course of 3 hours to find the
values reported in Table 1. The fact that $T_{1}^{\Delta}\gg T_{1}^{\Sigma}$
is likely to originate from the intentional difference in geometry between the
two modes, and hence their coupling to the coaxial output ports, as well as
the difference in detuning and coupling to the readout resonator Gambetta  _et
al._ (2011).

We use Ramsey interferometry in order to measure energy dispersion. A mode is
prepared in a superposition state using an $X_{\pi/2}$ pulse, allowed to idle
for time $\Delta t$, before a second $X_{\pi/2}$ pulse is applied, as shown in
the pulse sequence in Fig. 3 (a). We detune the frequency of the control
pulses from the average mode frequency by approximately 3.5 MHz to prevent
aliasing for large frequency excursions. At each $\Delta t$, we sample 2500
times, taking approximately 100 ms to acquire. Typical transmon devices
exhibit a quasiparticle tunnelling rate of $0.01~{}\rm{\mu s^{-1}}$ Serniak
_et al._ (2018), therefore we expect to average over all possible parity
configurations during the measurement. As a result, we expect to observe four
frequency components in a Ramsey oscillation, corresponding to the four parity
configurations, shown in Fig. 2 (b).

Figure 3: Ramsey interferometry based measurement of charge bias
configuration. Experimental measurement of charge dispersion and four parity
configurations in the coaxial multi-mode qubit, with data points in black and
fits in solid lines. (a) Pulse sequence for the Ramsey interferometry
measurement. (b) (c) Ramsey oscillation performed on the $\Delta$-mode. The
solid line in (b) shows a decaying oscillation with frequency components, and
decay constant, obtained from fitting the Fourier transform of the
oscillation, shown in (c). The dashed line shows an exponential decay envelope
with this same fitted decay constant. The Fourier transform of the measured
$\Delta$-mode Ramsey oscillations shows four frequency peaks, corresponding to
the four possible parity configurations. The data is fitted to four Lorentzian
peaks, symmetric about the average frequency. We define the separation of the
inner peaks to be $\Delta f_{1}$, and the separation of the outer peaks to be
$\Delta f_{2}$. From this separation of frequency peaks we can extract the
charge bias configuration. (d) (e) Ramsey oscillation and Fourier transform of
measured oscillation performed on the $\Sigma$-mode.

In Fig. 3 we show example Ramsey oscillations measured for both the $\Delta$
\- mode and $\Sigma$ \- mode. Fig. 3 (c) and (e) show the fast fourier
transform (FFT) of the Ramsey oscillations, in which we observe four distinct
frequency components. We fit the FFT data to four identical Lorenztian peaks,
which are symmetric about the average frequency,
$\overline{E_{0\Sigma(\Delta)}}$. The separation of the two inner peaks
closest to the symmetry point, is labelled as $\Delta f_{1}$, and the
separation of the two outermost peaks is labelled as $\Delta f_{2}$, as shown
in Fig. 3. Using our tight-binding model, we find that:

$\displaystyle\Delta f_{1}=\frac{\epsilon_{mn}}{h}\sin{\pi
n_{g\Sigma}}\sin{\pi n_{g\Delta}},$ (6) $\displaystyle\Delta
f_{2}=\frac{\epsilon_{mn}}{h}\cos{\pi n_{g\Sigma}}\cos{\pi n_{g\Delta}}$

This allows us to determine the charge configuration from the energy
dispersion. Repeating this measurement, we are able to track frequency
fluctuations due to correlated and anti-correlated charge noise dynamics over
extended periods of time.

Note that there are several technical shortcomings of our demonstration
experiment, which can be remedied in the future. Firstly, there is no gate
charge control in this current architecture, which prevents us from resolving
jumps or drifts larger than 0.5e for $n_{g1}$ and $n_{g2}$. This limits our
ability to determine a charge noise spectral density at this time Christensen
_et al._ (2019), but can be alleviated with the incorporation of local control
of static electric fields via gate electrodes. Secondly, the time taken to
acquire each Ramsey oscillation trace limits the ability to observe changes in
charge configuration faster than two minutes. This can be remedied using
higher fidelity readout with a parametric amplifier, or a more efficient
sampling of Ramsey delay times.

##  IV RESULTS

Figure 4: Time series measurement of charge dispersion. (a) Fourier transform
of repeated Ramsey oscillation experiments performed on the $\Sigma$-mode,
over a period of 2.5 hours as shown in Fig. 3. (b) Tracking of inner frequency
peak separation (solid line), $\Delta f_{1}$, and outer frequency peak
separation (dashed line), $\Delta f_{2}$, from time series measurement shown
in (a). Black vertical dashed lines highlight $T=0$ mins, $T=50$ mins and
$T=100$ mins. (c), (d), (e) Mapping of charge configuration to physical
location of charge on substrate at highlighted times shown in (b). Uncertainty
in position due to error in fitting of Fourier transform data, shown for
$1\sigma$ (solid line) and $2\sigma$ (dashed line), with the most likely
position indicated by the black dot. One quarter of the device is shown due to
the four-fold symmetry of the system.

Through repeated Ramsey interferometry of device A, monitored over the course
of 10 hours, we find a maximum dispersion of $\epsilon_{10}/h=4$ MHz
($\epsilon_{01}/h=8$ MHz) for the $\Sigma$-mode ($\Delta$-mode). This is
consistent with the order of magnitude of our predicted values of 4 MHz (4.1
MHz). The larger difference between $\epsilon_{10}$ and $\epsilon_{01}$ is due
to asymmetries in the Josephson energies of the junctions, caused by
fabrication imperfections.

In device B, we perform Ramsey experiments on the more charge sensitive
$\ket{01}-\ket{11}$ transition, and find no frequency beating up to a
resolution of 10 kHz, demonstrating that this device design iteration has a
suppressed sensitivity to charge noise, by at least a factor of 400 compared
to device A. This measurement is shown in Appendix E.

However, this suppression of charge sensitivity in device B comes at a cost of
lower mode anharmonicities, and state-dependent shifts, as shown in Table I.
This reduces the maximum speed at which single-qubit gates, and entangling
operations between modes can be completed.

##  V SPATIALLY RESOLVED CHARGE DETECTION

We now use device A to demonstrate a proof-of-concept detector of localised
charge fluctuations. By repeating the Ramsey interferometry we can track the
frequency dispersion $\Delta f_{1}$ and $\Delta f_{2}$, and using our tight-
binding model, we are able to infer a time-series measurement of charge
configuration of $n_{g\Sigma}$ and $n_{g\Delta}$. We find the expected spatial
charge sensitivity of these sum and difference modes using an electrostatic
simulation, shown in Fig. 1 (d) and (e). As the induced gate-charge offsets
$n_{g\Sigma}$ and $n_{g\Delta}$ exhibit different spatial sensitivities, we
can deduce the position of a potential surface-charge on the device, using a
biangulation method. For a given value of $n_{g\Sigma}$, we use the simulation
data shown in Fig. 1 (d) to identify an area where a charge ($1e$) would
induce a gate-charge offset of that value. We repeat this for a corresponding
value of $n_{g\Delta}$ with the simulation data shown in Fig. 1 (e). The
overlap of these two individually obtained areas allows us to identify the
location of a charge ($1e$) on the surface of the substrate. However, the use
of two modes produces an ambiguity in the quadrant in which any surface charge
is located. Future devices could use three or more islands, and incorporate a
symmetry breaking geometry, in order to more accurately triangulate surface
charge position.

We perform a demonstrative tracking experiment using repeated Ramsey
interferometry measurement over a 150 minute period, as shown in Fig. 4. We
observe both slow frequency drifts corresponding to fluctuating charge
configuration, as well as a singular large frequency jump, indicative of non-
equilibrium charge dynamics Martinis (2021). We convert the measured frequency
dispersion to give a spatial estimation of charge location at time T = 0, T =
50, and T = 100 minutes, indicated by the dashed vertical lines on Fig. 4 (b).

From $T=0$ to $T=50$ minutes, we observe a slow drift in surface charge moving
outwards away from the inner islands of the device. The uncertainty in
position from $T=50$ onwards is high, as the device is not sensitive to
spatial fluctuations far from the two inner islands, shown by the lower
gradient in the simulation of induced charge in Fig. 1. In this period we also
observe a much clearer signal in the raw FFT data of Fig. 4 (a), consistent
with the charge configuration remaining stable throughout the individual
Ramsey oscillation measurement.

From $T=100$ onwards, we observe a shift in the surface charge distribution
towards the two inner islands of the device. In this region, we have a much
higher spatial resolution, shown by the larger gradient in simulation of
induced charge from Fig. 1. As such, any small movement in surface charge
distribution will cause a larger shift in charge configuration of the device.
We observe this as a noisier measurement of frequencies shown in Fig. 4 (b).
This is also shown in the raw FFT data where the peaks are less clear in this
region, consistent with the charge configuration changing during the course of
the measurement. With an improved measurement rate, this detector could be
used for observing spatial charge fluctuations with a resolution of less than
$100$ $\mu m$.

##  VI CONCLUSION

In this work, we have investigated charge sensitivity in a two-mode transmon.
Our results show observations of multiple charge parity configurations, and
show agreement with a predictive theory for charge sensitivity using a tight
binding approximation. In addition, we show this sensitivity can be suppressed
for high coherence quantum computing applications.

We use a two-mode device in the charge-sensitive regime (device A) to
demonstrate proof-of-principle spatially-resolved charge detection. Combining
a similar device with high fidelity state readout in future would enable
detection of parity jumps between three-superconducting islands, and determine
how quasi-particle tunnelling events affect energy dissipation in multiple
modes. If local control of static electric field were incorporated in a future
device via gate electrodes, it would become possible to unambiguously
translate frequency fluctuations to precise charge configurations, for precise
detection of local differential charge-noise vs. global charge-noise. We hence
propose such multi-mode charge-sensitive qubits as potentially powerful tools
for ongoing investigations into sources of decoherence in superconducting
quantum devices.

###### Acknowledgements.

This work has received funding from the United Kingdom Engineering and
Physical Sciences Research Council under Grants No. EP/J013501/1,
EP/M013243/1, EP/N015118/1, EP/T001062/1 and from Oxford Quantum Circuits
Limited. S.D.F. acknowledges support from the Swiss Study Foundation and the
Bakala Foundation. B.V. acknowledges support from an EU Marie Sklodowska-Curie
fellowship.

##  VII APPENDIX

###  VII.1 Circuit Quantisation

The two-mode coaxial transmon has the Hamiltonian:

$\hat{H}=8E_{C_{\Sigma}}\hat{n}_{\Sigma}^{2}+8E_{C_{\Delta}}\hat{n}_{\Delta}^{2}-2E_{J}\cos{\frac{\hat{\varphi_{\Sigma}}}{2}}\cos{\frac{\hat{\varphi_{\Delta}}}{2}}$
(A1)

where $E_{C_{\Sigma}}=E_{c}+E_{p}/2=e^{2}/2(C^{*}-C_{m})$ and
$E_{C_{\Delta}}=E_{c}-E_{p}/2=e^{2}/2(C^{*}+C_{m})$ are the charging energy of
the sum and difference modes respectively, as described in the main text. This
shows independent charge variables $\hat{n}_{\Sigma},~{}\hat{n}_{\Delta}$ with
a coupled potential energy term
$U(\varphi_{\Sigma},\varphi_{\Delta})=-2E_{J}\cos{\frac{\varphi_{\Sigma}}{2}}\cos{\frac{\varphi_{\Delta}}{2}}$.
We can predict the behaviour of this system in the transmon regime $E_{p}\ll
E_{c}\ll E_{J}$, by expanding the cosine potential terms. To fourth order in
$\varphi_{\Sigma},~{}\varphi_{\Delta}$, we find:

$\displaystyle\hat{H}=8E_{C_{\Sigma}}\hat{n}_{\Sigma}^{2}+8E_{C_{\Delta}}\hat{n}_{\Delta}^{2}$
(A2)
$\displaystyle+2E_{J}(1/8\hat{\varphi}_{\Sigma}^{2}+1/8\hat{\varphi}_{\Delta}^{2}$
$\displaystyle-1/{24*4^{4}}\hat{\varphi}_{\Sigma}^{4}-1/{24*4^{4}}\hat{\varphi}_{\Delta}^{4}$
$\displaystyle-1/{16}\hat{\varphi}_{\Sigma}^{2}\hat{\varphi}_{\Delta}^{2})$

Assuming a weakly nonlinear oscillator-like behaviour, we take
$\hat{\varphi}_{\Sigma(\Delta)}=\sqrt{\frac{2E_{J}}{E_{C_{\Sigma}(\Delta)}}}(\hat{a}_{\Sigma(\Delta)}^{\dagger}+\hat{a}_{\Sigma(\Delta)})$,
and
$\hat{n}_{\Sigma(\Delta)}=\frac{i}{2}\sqrt{\frac{E_{C_{\Sigma}(\Delta)}}{2E_{J}}}(\hat{a}_{\Sigma(\Delta)}^{\dagger}-\hat{a}_{\Sigma(\Delta)})$,
where $\hat{a}_{\Sigma(\Delta)}^{(\dagger)}$ is the annihilation (creation)
operator for the $\Sigma$-mode ($\Delta$-mode). Keeping only counter-rotating
terms produces:

$\displaystyle\hat{H}/\hbar$
$\displaystyle=\omega_{\Sigma}\hat{a}^{\dagger}_{\Sigma}\hat{a}_{\Sigma}+\omega_{\Delta}\hat{a}^{\dagger}_{\Delta}\hat{a}_{\Delta}$
(A3)
$\displaystyle-\eta_{\Sigma}/2*\hat{a}^{\dagger}_{\Sigma}{}^{2}\hat{a}_{\Sigma}{}^{2}-\eta_{\Delta}/2*\hat{a}^{\dagger}_{\Delta}{}^{2}\hat{a}_{\Delta}{}^{2}$
$\displaystyle-4\sqrt{\eta_{\Sigma}\eta_{\Delta}}*\hat{a}^{\dagger}_{\Sigma}\hat{a}_{\Sigma}\hat{a}^{\dagger}_{\Delta}\hat{a}_{\Delta}$

where $\eta_{\Sigma(\Delta)}$ is the anharmonicity of the $\Sigma$-mode
($\Delta$-mode). Notice that the state dependent shift term
$4\sqrt{\eta_{\Sigma}\eta_{\Delta}}\equiv\chi_{\Sigma\Delta}$ can be greater
than the individual mode anharmonicities $\eta_{\Sigma(\Delta)}$.

The full Hamiltonian of the system shown in Fig. 1 (a) is given by
$\hat{H}_{full}=\hat{H}+\hat{H}_{r}$, where $\hat{H}_{r}$ describes the
coupling between each mode and the readout resonator, given by:

$\displaystyle\hat{H}_{r}/\hbar$
$\displaystyle=\omega_{r}\hat{a}^{\dagger}_{r}\hat{a}_{r}+g_{\Sigma(\Delta)r}(\hat{a}_{r}^{\dagger}+\hat{a}_{r})(\hat{a}_{\Sigma(\Delta)}^{\dagger}+\hat{a}_{\Sigma(\Delta)})$
(A4)

where $\hat{a}^{(\dagger)}_{r}$ is the annihilation (creation) operator for
the resonator mode, $\omega_{r}$ is the resonator mode frequency, and
$g_{\Sigma(\Delta)r}$ is the coupling between the resonator mode and the
$\Sigma$-mode ($\Delta$-mode). This weak coupling between each transmon mode
and the resonator mode leads to a mode-state-dependent resonator frequency
shift $\chi_{\Sigma(\Delta)r}$ as defined in ref. Rahamim  _et al._ (2017).
This allows us to perform simultaneous dispersive readout of the multi-mode
state of the device.

###  VII.2 Device Parameters

In Table 1 we show the measured parameters of the devices used in these
experiments. Fig. 5 shows the geometry of both the charge sensitive (A) and
insensitive (B) two-mode coaxial transmons.

Table 1: Device Parameters | Device A | Device B
---|---|---
LC Resonator |  |
Frequency $f_{r}$ [GHz] | 9.72 | 9.32
Linewidth $\kappa_{r}/2\pi$ [MHz] | 2.8 | 1.9
Dispersive Shift $2\chi_{\Delta r}/2\pi$ [MHz] | 3.9 | 2.1
Dispersive Shift $2\chi_{\Sigma r}/2\pi$ [MHz] | 4.9 | 1.7
$\Delta$ \- Mode |  |
Transition Frequency $\omega_{\Delta}/2\pi$ [GHz] | 5.51 | 4.58
Anharmonicity $\eta_{\Delta}/2\pi$ [MHz] | 380 | 104
$T_{1}$ [$\mu s$] | 40.7 | 50.9
$T_{2}$ Echo [$\mu s$] | 33.8 | 24.8
$\Sigma$ \- Mode |  |
Transition Frequency $\omega_{\Sigma}/2\pi$ [GHz] | 6.71 | 5.74
Anharmonicity $\eta_{\Sigma}/2\pi$ [MHz] | 340 | 144
$T_{1}$ [$\mu s$] | 13.5 | 26.5
$T_{2}$ Echo [$\mu s$] | 13.1 | 19.4
Cross Kerr Shift $\chi_{\Sigma\Delta}/2\pi$ [MHz] | 500 | 269
Figure 5: Comparison of device designs. Device A (left) has an inner island
radius of $125~{}\rm\mu m$, and inner island gap of $120~{}\rm\mu m$. Device B
(right) has an inner island radius of $220~{}\rm\mu m$, and inner island gap
of $125~{}\rm\mu m$. Both devices have the same outer island geometry, with
inner radius of $389.5~{}\rm\mu m$, and outer radius of $489.5~{}\rm\mu m$.

###  VII.3 Wavefunctions

Fig. 6 shows the wavefunctions of the six lowest energy eigenstates of the
two-mode coaxial transmon Hamiltonian in Eqn. A1. This Hamiltonian was
constructed and solved in the charge basis and then transformed to the phase
basis to obtain the wavefunctions plotted, using the package scqubits
Groszkowski and Koch (2021).

Figure 6: Wavefunctions of the six lowest energy eigenstates of the two-mode
coaxial transmon. Diagonal dashed lines indicate the summation
($(\varphi_{1}+\varphi_{2})$, orange) and difference
($(\varphi_{1}-\varphi_{2})$, green) phases. States labelled as $\ket{nm}$
where $n(m)$ corresponds to the number of excitations in the $\Sigma$-mode
($\Delta$-mode).

###  VII.4 Tight-Binding Model

The Hamiltonian of Eqn. A1 can be written as $\hat{H}=\hat{T}+\hat{U}$, where
the kinetic term
$\hat{T}=8E_{C_{\Sigma}}\hat{n}_{\Sigma}^{2}+8E_{C_{\Delta}}\hat{n}_{\Delta}^{2}$,
and the potential term can be rewritten as:

$\displaystyle U(\varphi_{\Sigma},\varphi_{\Delta})$
$\displaystyle=-2E_{J}\cos{\frac{\varphi_{\Sigma}}{2}}\cos{\frac{\varphi_{\Delta}}{2}}$
(A5)
$\displaystyle=-2E_{J}\cos{\frac{\varphi_{\Sigma}}{2}}-2E_{J}\cos{\frac{\varphi_{\Delta}}{2}}$
$\displaystyle-4E_{J}\sin^{2}{\frac{\varphi_{\Sigma}}{4}}\sin^{2}{\frac{\varphi_{\Delta}}{4}}$
$\displaystyle=U_{0}(\varphi_{\Sigma},\varphi_{\Delta})+U_{1}(\varphi_{\Sigma},\varphi_{\Delta})$

where
$U_{0}=-2E_{J}\cos{\frac{\varphi_{\Sigma}}{2}}-2E_{J}\cos{\frac{\varphi_{\Delta}}{2}}$,
is the potential term for two uncoupled transmons, and
$U_{1}=-4E_{J}\sin^{2}{\frac{\varphi_{\Sigma}}{4}}\sin^{2}{\frac{\varphi_{\Delta}}{4}}$,
is a perturbation coupling the two transmon modes together. The potential
$U_{0}$ is $4\pi$ periodic in $\varphi_{\Sigma}$ and $\varphi_{\Delta}$, and
the additional potential term $U_{1}$ introduces lattice sites at
$(\varphi_{\Sigma},\varphi_{\Delta})=(\pm 2\pi,\pm 2\pi)$.

Given this lattice structure, we use Bloch’s theorem to pick an approximate
solution to the Schrödinger equation $\hat{H}\ket{\psi}=E\ket{\psi}$ as:

$\displaystyle\psi(\vec{\varphi})$
$\displaystyle=\frac{1}{\sqrt{2}}(\psi^{(1)}(\vec{\varphi})+\psi^{(2)}(\vec{\varphi}))$
(A6)
$\displaystyle=\frac{1}{\sqrt{2}}e^{i\vec{k}.\vec{\varphi}}(u_{k}(\vec{\varphi})+e^{i\vec{k}.\vec{a}}u_{k}(\vec{\varphi}-\vec{a}))$

where $\vec{k}$ is the wave vector of the wavefunction $\psi$, $\vec{\varphi}$
is the position vector given by $(\varphi_{\Sigma},\varphi_{\Delta})$,
$\vec{a}$ denotes the lattice site introduced by the additional potential
$U_{1}$, and $u_{k}$ is a $4\pi$ periodic function. Due to the symmetry of the
system, we only need to consider the lattice sites located at
$\vec{\varphi}=(0,0)$ and $\vec{\varphi}=(2\pi,2\pi)$.

We calculate the energies of the system with this approximate wavefunction
using:

$\displaystyle
E(k)=\frac{\bra{\psi_{k}}\hat{H}\ket{\psi_{k}}}{\bra{\psi_{k}}\ket{\psi_{k}}}$
(A7)

leading to the dispersion relation:

$\displaystyle E(n_{g\Sigma},n_{g\Delta})\approx E_{0}+\cos{\pi
n_{g\Sigma}}\cos{\pi n_{g\Delta}}(\gamma-\alpha\beta)$ (A8)

where the wave vector $\vec{k}=(n_{g\Sigma}/2,n_{g\Delta}/2)$, and the
quantities $\alpha$, $\gamma$ and $\beta$ are the tight binding coefficients
given by:

$\displaystyle\alpha=4\int_{0}^{2\pi}u_{k}(\vec{\varphi})u_{k}(\vec{\varphi}-2\pi)\,d^{2}\vec{\varphi}$
(A9)

$\displaystyle\beta=-16E_{J}\int_{0}^{2\pi}|u_{k}(\vec{\varphi})|^{2}\sin^{2}{\frac{\varphi_{\Sigma}}{4}}\sin^{2}{\frac{\varphi_{\Delta}}{4}}\,d^{2}\vec{\varphi}$
(A10)

$\displaystyle\gamma=-16E_{J}\int_{0}^{2\pi}u_{k}(\vec{\varphi})u_{k}(\vec{\varphi}-2\pi)\sin^{2}{\frac{\varphi_{\Sigma}}{4}}\sin^{2}{\frac{\varphi_{\Delta}}{4}}\,d^{2}\vec{\varphi}$
(A11)

The dominant term here is $\gamma$, which describes the bond energy between
wavefunctions at adjacent lattice sites, also known as the two center
integral. The $\alpha$ and $\beta$ terms describe the overlap integral between
wavefunctions on adjacent lattice sites, and the energy shift due to the
potential on neighbouring lattice sites respectively. These two terms are
small compared to $\gamma$ and so can be neglected, leading to the approximate
form of the maximum charge dispersion $\epsilon_{mn}/4\approx\gamma$.

This form of the charge dispersion $\epsilon_{mn}$ can either be calculated
numerically, using the wavefunctions shown in Appendix C, or analytically
using a semi analytical wavefunction approach to obtain the functional form of
$u_{k}$ Catelani  _et al._ (2011). Using this approach, we arrive at the
analytical form for the maximum charge dispersion shown in Eqn. 5.

###  VII.5 Ramsey Experiments On The Charge Insensitive Two-Mode Transmon

In order to demonstrate a suppression of offset charge sensitivity in device
B, we perform Ramsey interferometry measurements on the $\ket{01}-\ket{11}$
transition. From our predictive model, we expect a maximum charge dispersion
of $\epsilon_{11}/h\approx 10$ kHz.

In order to measure this dispersion, we first prepare the system in the
$\ket{01}$ state using a $X_{\pi}$ pulse, at a frequency $\omega_{\Delta}$. We
then place the system in a superposition using a $X_{\pi/2}$ pulse, at a
frequency $\omega_{\Sigma}-\chi$. The system is then left to idle for time
$\Delta t$, before a second $X_{\pi/2}$ pulse is applied, as shown in the
pulse scheme, shown in Fig. 7 (a). As with the previous Ramsey interferometry
measurements, the experiment is sampled many times, and so we expect to
average over all possible parity configuration.

In Fig. 7, we show the measured Ramsey oscillation, and FFT, for this
transition. In contrast to the measurements performed on the charge sensitive
device A, shown in Fig. 3, the oscillation is comprised of a single frequency
component. This is demonstrated in the FFT of the decaying oscillation, shown
in Fig. 7 (c), where we are unable to distinguish a further frequency
component, to a resolution of 10 kHz. This demonstrates this device design
iteration has a suppressed sensitivity to charge noise.

Figure 7: Ramsey interferometry on the $\ket{01}-\ket{11}$ transition of the
charge insensitive two-mode transmon (device B). (a) Pulse sequence for the
Ramsey interferometry measurement. (b)(c) Ramsey oscillation and Fourier
transform of the measurement, with data points in black, and fits in solid
lines. The data in (b) is fitted to a decaying oscillation, and the dashed
lines show an exponential decay envelope with the same decay constant. The
data in (c) is fitted to a single Lorentzian peak.

## References

  * DiVincenzo (2000) D. P. DiVincenzo, Fortschritte der Physik 48, 771 (2000).
  * Devoret and Schoelkopf (2013) M. H. Devoret and R. J. Schoelkopf, Science 339, 1169 LP (2013).
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, and et al., Nature 574, 505 (2019).
  * Koch  _et al._ (2007) J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A 76, 042319 (2007).
  * Manucharyan  _et al._ (2009) V. E. Manucharyan, J. Koch, L. I. Glazman, and M. H. Devoret, Science 326, 113 (2009).
  * Kou  _et al._ (2017) A. Kou, W. C. Smith, U. Vool, R. T. Brierley, H. Meier, L. Frunzio, S. M. Girvin, L. I. Glazman, and M. H. Devoret, Phys. Rev. X 7, 031037 (2017).
  * Groszkowski  _et al._ (2018) P. Groszkowski, A. D. Paolo, A. L. Grimsmo, A. Blais, D. I. Schuster, A. A. Houck, and J. Koch, New Journal of Physics 20, 043053 (2018).
  * Gyenis  _et al._ (2021) A. Gyenis, P. S. Mundada, A. Di Paolo, T. M. Hazard, X. You, D. I. Schuster, J. Koch, A. Blais, and A. A. Houck, PRX Quantum 2, 010339 (2021).
  * Sundaresan  _et al._ (2020) N. Sundaresan, I. Lauer, E. Pritchett, E. Magesan, P. Jurcevic, and J. M. Gambetta, PRX Quantum 1, 020318 (2020).
  * Roy  _et al._ (2017) T. Roy, S. Kundu, M. Chand, S. Hazra, N. Nehra, R. Cosmic, A. Ranadive, M. P. Patankar, K. Damle, and R. Vijay, Phys. Rev. Applied 7, 054025 (2017).
  * Srinivasan  _et al._ (2011) S. J. Srinivasan, A. J. Hoffman, J. M. Gambetta, and A. A. Houck, Phys. Rev. Lett. 106, 083601 (2011).
  * Richer  _et al._ (2017) S. Richer, N. Maleeva, S. T. Skacel, I. M. Pop, and D. DiVincenzo, Phys. Rev. B 96, 174520 (2017).
  * Noguchi  _et al._ (2020) A. Noguchi, A. Osada, S. Masuda, S. Kono, K. Heya, S. P. Wolski, H. Takahashi, T. Sugiyama, D. Lachance-Quirion, and Y. Nakamura, Phys. Rev. A 102, 062408 (2020).
  * Smith  _et al._ (2020) W. C. Smith, A. Kou, X. Xiao, U. Vool, and M. H. Devoret, npj Quantum Information 6, 8 (2020).
  * Finck  _et al._ (2021) A. D. K. Finck, S. Carnevale, D. Klaus, C. Scerbo, J. Blair, T. G. McConkey, C. Kurter, A. Carniol, G. Keefe, M. Kumph, and O. E. Dial, (2021), arXiv:2105.11495 [quant-ph] .
  * Gambetta  _et al._ (2011) J. M. Gambetta, A. A. Houck, and A. Blais, Phys. Rev. Lett. 106, 030502 (2011).
  * Roy  _et al._ (2018) T. Roy, M. Chand, A. Bhattacharjee, S. Hazra, S. Kundu, K. Damle, and R. Vijay, Phys. Rev. A 98, 052318 (2018).
  * Zhang  _et al._ (2017) G. Zhang, Y. Liu, J. J. Raftery, and A. A. Houck, npj Quantum Information 3, 1 (2017).
  * Pashkin  _et al._ (2003) Y. A. Pashkin, T. Yamamoto, O. Astafiev, Y. Nakamura, D. V. Averin, and J. S. Tsai, Nature 421, 823 (2003).
  * Gustafsson  _et al._ (2013) M. V. Gustafsson, A. Pourkabirian, G. Johansson, J. Clarke, and P. Delsing, Phys. Rev. B 88, 245410 (2013).
  * Schreier  _et al._ (2008) J. A. Schreier, A. A. Houck, J. Koch, D. I. Schuster, B. R. Johnson, J. M. Chow, J. M. Gambetta, J. Majer, L. Frunzio, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. B 77, 180502 (2008).
  * Christensen  _et al._ (2019) B. G. Christensen, C. D. Wilen, A. Opremcak, J. Nelson, F. Schlenker, C. H. Zimonick, L. Faoro, L. B. Ioffe, Y. J. Rosen, J. L. DuBois, B. L. T. Plourde, and R. McDermott, Phys. Rev. B 100, 140503 (2019).
  * Wilen  _et al._ (2021) C. D. Wilen, S. Abdullah, N. A. Kurinsky, C. Stanford, L. Cardani, G. D’Imperio, C. Tomei, L. Faoro, L. B. Ioffe, C. H. Liu, and et al., Nature 594, 369 (2021).
  * Tennant  _et al._ (2021) D. M. Tennant, L. A. Martinez, C. D. Wilen, R. McDermott, J. L. DuBois, and Y. J. Rosen, (2021), arXiv:2106.08406 [quant-ph] .
  * Martinis (2021) J. M. Martinis, (2021), arXiv:2012.06137 [quant-ph] .
  * Dumur  _et al._ (2015) E. Dumur, B. Küng, A. K. Feofanov, T. Weissl, N. Roch, C. Naud, W. Guichard, and O. Buisson, Phys. Rev. B 92, 020515 (2015).
  * Hazra  _et al._ (2020) S. Hazra, K. V. Salunkhe, A. Bhattacharjee, G. Bothara, S. Kundu, T. Roy, M. P. Patankar, and R. Vijay, Applied Physics Letters 116, 152601 (2020).
  * Rahamim  _et al._ (2017) J. Rahamim, T. Behrle, M. J. Peterer, A. Patterson, P. A. Spring, T. Tsunoda, R. Manenti, G. Tancredi, and P. J. Leek, Applied Physics Letters 110, 222602 (2017).
  * Serniak  _et al._ (2018) K. Serniak, M. Hays, G. de Lange, S. Diamond, S. Shankar, L. D. Burkhart, L. Frunzio, M. Houzet, and M. H. Devoret, Phys. Rev. Lett. 121, 157701 (2018).
  * Catelani  _et al._ (2011) G. Catelani, R. J. Schoelkopf, M. H. Devoret, and L. I. Glazman, Phys. Rev. B 84, 064517 (2011).
  * Spring  _et al._ (2021) P. A. Spring, S. Cao, T. Tsunoda, G. Campanaro, S. D. Fasciati, J. Wills, V. Chidambaram, B. Shteynas, M. Bakr, P. Gow, L. Carpenter, J. Gates, B. Vlastakis, and P. J. Leek, (2021), arXiv:2107.11140 [quant-ph] .
  * Minev  _et al._ (2021) Z. K. Minev, Z. Leghtas, S. O. Mundhada, L. Christakis, I. M. Pop, and M. H. Devoret, (2021), arXiv:2010.00620 [quant-ph] .
  * Groszkowski and Koch (2021) P. Groszkowski and J. Koch, (2021), arXiv:2107.08552 [quant-ph] .


# Optimum phase estimation with two control qubits

Peyman Najafi  Donders Institute for Neuroscience, Radboud University,
Nijmegen, the Netherlands  School of Mathematical and Physical Sciences,
Macquarie University, Sydney, NSW 2109, Australia  Pedro C. S. Costa  School
of Mathematical and Physical Sciences, Macquarie University, Sydney, NSW 2109,
Australia  Dominic W. Berry  Electronic mail<EMAIL_ADDRESS>School
of Mathematical and Physical Sciences, Macquarie University, Sydney, NSW 2109,
Australia

###### Abstract

Phase estimation is used in many quantum algorithms, particularly in order to
estimate energy eigenvalues for quantum systems. When using a single qubit as
the probe (used to control the unitary we wish to estimate the eigenvalue of),
it is not possible to measure the phase with a minimum mean-square error. In
standard methods, there would be a logarithmic (in error) number of control
qubits needed in order to achieve this minimum error. Here show how to perform
this measurement using only two control qubits, thereby reducing the qubit
requirements of the quantum algorithm. Our method corresponds to preparing the
optimal control state one qubit at a time, while it is simultaneously consumed
by the measurement procedure.

##  I Introduction

Quantum phase estimation was originally applied in quantum algorithms for the
task of period finding, as in Shor’s algorithm Shor (1994). Later, quantum
phase estimation was applied to the task of estimating eigenvalues for
Hamiltonians in quantum chemistry Aspuru-Guzik  _et al._ (2005). The
appropriate way to perform quantum phase estimation is different between these
applications, due to the costing of the operations. In particular, for
estimating eigenvalues, the cost of Hamiltonian simulation is (at least)
proportional to the time of evolution, so the phase estimation procedure
should attempt to minimise the total evolution time. At the same time the
mean-square error in the estimate should be minimised.

As part of the phase estimation, the inverse quantum Fourier transform is
used. This operation can be decomposed into a ‘semiclassical’ form Griffiths
and Niu (1996), where one performs measurements on the control qubits in
sequence, with rotations controlled according to the results of previous
measurements. In the form of phase estimation as in Shor’s algorithm, the
control qubits would be an equal superposition state, which is just a tensor
product of $|+\rangle$ states on the individual qubits. In that scenario, only
one control qubit need be used at a time, because it can be prepared in the
$|+\rangle$ state, used as a control, then rotated and measured before the
next qubit is used.

This procedure with the control qubits in $|+\rangle$ states gives a
probability distribution for the error as a sinc function, which has a
significant probability for large errors. That is still suitable for Shor’s
algorithm, because it is possible to take large powers of the operators with
relatively small cost, which suppresses the phase measurement error. On the
other hand, for quantum chemistry where there is a cost of Hamiltonian
simulation proportional to time, the large error of the sinc is a problem.
Then it is more appropriate to use qubits in an entangled state Babbush  _et
al._ (2018), which was originally derived in an optical context in 1996 Luis
and Peřina (1996).

In 2000 we analysed the problem of how to perform measurements on these states
in a Mach-Zehnder interferometer Berry and Wiseman (2000). The same year, Jon
Dowling introduced NOON states in the context of lithography Boto  _et al._
(2000), and then in 2002 showed how NOON states may be used in interferometry
for phase measurement Kok  _et al._ (2002); Lee  _et al._ (2002). A drawback
to using NOON states is that they are highly sensitive to loss. In 2010 one of
us (DWB) visited Jon Dowling’s group to work on the problem of how to generate
states that are more resistant to loss and effectively perform measurements
with them. This resulted in the publication (separately from Jon) Dinani and
Berry (2014), followed by our first joint publication Dinani  _et al._ (2016).
We continued collaborating with Jon for many years on phase measurement Huang
_et al._ (2017), as well as state preparation Motes  _et al._ (2016), and
Boson-sampling inspired cryptography Huang  _et al._ (2021).

In separate work, we showed how to combine results from multiple NOON states
in order to provide highly accurate phase measurements suitable for quantum
algorithms Higgins  _et al._ (2007). Phase measurement via NOON states is
analogous to taking a $|+\rangle$ state and performing a controlled $U^{N}$ on
a target system in quantum computing. The photons in the arms of the
interferometer are analogous to the control qubit in quantum computing, with
the phase shift from $U^{N}$ instead arising from an optical phase shift
between the arms of the interferometer. The NOON state gives very high
frequency variation of the probability distribution for the phase, rather than
a probability distribution with a single peak. In 2007 we showed how to
combine the results from NOON states with different values of $N$ in order to
provide a phase measurement analogous to the procedure giving a sinc
distribution in quantum algorithms Higgins  _et al._ (2007). (It was
experimentally demonstrated with multiple passes through the phase shift
rather than NOON states.)

A further advance in Higgins  _et al._ (2007) was to show how to use an
adaptive procedure, still with individual $|+\rangle$ states, in order to give
the ‘Heisenberg limited’ phase estimate. That is, rather than the mean-square
error scaling as it does for the sinc, it scales as it does for the optimal
(entangled) control state. This procedure still only uses a single control
qubit at a time, so is suitable for using in quantum algorithms where the
number of qubits available is strongly limited; this is why it was used, for
example, in Kivlichan  _et al._ (2020). On the other hand, although it gives
the optimal scaling, the constant factor is not optimal, and improved
performance is provided by using the optimal entangled state.

In this paper we show how to achieve the best of both worlds. That is, we show
how to provide the optimal phase estimate (with the correct constant factor),
while only increasing the number of control qubits by one. It is therefore
suitable for quantum algorithms with a small number of qubits, while enabling
the minimum complexity for a given required accuracy.

In Section II we discuss the optimal state for phase estimation and how its
usage can be combined with the semiclassical quantum Fourier transform. Then
in Section III, we introduce a orthogonal basis of states for subsets of
qubits, and prove a recursive form. Finally, in Section IV we show how the
recursive form can be translated into a sequence of two-qubit unitaries to
create the optimal state.

##  II Phase measurement using optimal quantum states

###  II.1 The optimal states

The optimal states for phase estimation from Luis and Peřina (1996) are of the
form

$|\psi_{\rm
opt}\rangle=\sqrt{\frac{2}{N+2}}\sum_{n=0}^{N}\sin\left(\frac{\pi(n+1)}{N+2}\right)|n\rangle,$
(1)

where $N$ is the total photon number in two modes, and $n$ is the photon
number in one of the modes, as for example in a Mach-Zehnder interferometer.
It is also possible to consider the single-mode case where $N$ is a maximum
photon number and $|n\rangle$ a Fock state.

In either case a physical phase shift of $\phi$ results in a state of the form

$|\psi_{\rm
opt}\rangle=\sqrt{\frac{2}{N+2}}\sum_{n=0}^{N}e^{in\phi}\sin\left(\frac{\pi(n+1)}{N+2}\right)|n\rangle.$
(2)

The ideal ‘canonical’ phase measurement is then a positive operator-valued
measure (POVM) with elements Sanders  _et al._ (1997)

$\frac{N+1}{2\pi}|\check{\phi}\rangle\langle\check{\phi}|\,d\phi,$ (3)

where

$|\check{\phi}\rangle=\frac{1}{\sqrt{N+1}}\sum_{n=0}^{N}e^{in\check{\phi}}|n\rangle.$
(4)

Here we are using $\check{\phi}$ for the result of the measurement, as
distinct from the actual phase $\phi$. Such a canonical measurement typically
cannot be implemented using standard linear optical elements, though it can be
approximated with adaptive measurements Berry and Wiseman (2000).

It is easily seen that the error distribution after the measurement is then

$\frac{1}{\pi(N+2)}\left(\frac{\cos((\check{\phi}-\phi)(1+N/2))\sin(\pi/(2+N))}{\cos(\pi/(2+N))-\cos(\check{\phi}-\phi)}\right)^{2}.$
(5)

In contrast, if one were to use the state with an equal distribution over
basis states, then the error probability distribution would be close to a sinc

$\frac{1}{2\pi(N+1)}\frac{\sin^{2}((N+1)(\check{\phi}-\phi)/2)}{\sin^{2}((\check{\phi}-\phi)/2)}.$
(6)

The error distributions for these two states are shown in Figure 1. The
central peak for the equal superposition state is a little narrower, but it
has large tails in the distribution, whereas the probabilities of large errors
for the optimal state are strongly suppressed.

Figure 1:  The probability distribution for the error in phase measurements
with $N=10$ and the optimal state (1) (blue) and the equal superposition state
(orange). The left shows the linear scale and the right shows a log plot.

The optimal state (1) is optimal for minimising a slightly different measure
of error than usual. The Holevo variance for a distribution can be taken as
Holevo (1984)

$|\langle e^{i\check{\phi}}\rangle|^{-2}-1.$ (7)

This measure has the advantages that it is naturally modulo $2\pi$, as is
suitable for phase, and approaches infinity for a flat distribution (with no
phase information). Moreover it approaches the usual variance for suitably
narrowly peaked distributions. To eliminate biased estimate, one can
alternatively use the measure

$\langle\cos(\check{\phi}-\phi)\rangle^{-2}-1.$ (8)

This measure is analogous to the mean-square error. One could also take the
measure, as in Luis and Peřina (1996),

$2[1-\langle\cos(\check{\phi}-\phi)\rangle],$ (9)

and the optimisation problem is equivalent. The optimal state (1) gives a
minimam Holevo variance of

$\tan^{2}\left(\frac{\pi}{N+2}\right).$ (10)

It is also possible to consider minimisation of the mean-square error, but
there is not a simple analytic solution Berry  _et al._ (2012).

###  II.2 Phase measurement with the inverse Fourier transform

In the case of phase measurements in quantum computing, $\phi$ would instead
be obtained from a unitary operator $U$ with eigenvalue $e^{i\phi}$. If the
target system is in the corresponding eigenstate of $U$, denoted
$|\phi\rangle$, then if state $|n\rangle$ is used to control application of
$U^{n}$, then the $\phi$-dependent state from Eq. (2) is again obtained. In
practice, the integer $n$ is represented in binary in ancilla qubits. Then the
most-significant bit, $n_{1}$, is used to control $U$, the next most
significant bit, $n_{2}$, is used to control $U^{2}$, and so forth. In
general,

$n=\sum_{j=1}^{m}n_{j}2^{j-1},$ (11)

$|n_{j}\rangle$ is used to control $U^{2^{j-1}}$. This procedure is depicted
in Figure 2.

@R=2em @C=1em &—n_1⟩ ⋯ 4
—n_2⟩ ⋯ 3
⋮
—n_m⟩ 1 ⋯
—ϕ⟩ U^2^m-1 ⋯ U^2 U^1

Figure 2:  The circuit for a controlled power of $U$, where $|n\rangle$ on the
control qubits gives $U^{n}$ on the target register. With the target prepared
in an eigenstate $|\phi\rangle$ of $U$, the phase shift $e^{in\phi}$ is
obtained.

Here we have taken $m$ to be the number of bits. In practice, it is convenient
to take $N-1$ to be a power of 2, so $N=2^{m}-1$. In order to estimate the
phase, one wishes to perform the canonical measurement on the ancilla qubits.
To explain this, it is convenient to consider the POVM with $N+1$ states
$|\check{\phi}_{j}\rangle$ with $\check{\phi}_{j}=2\pi j/(N+1)$ for $j=0$ to
$N$. Then the states $|\check{\phi}\rangle$ are mutually orthogonal. Such a
projective measurement can then be obtained if one can perform the unitary
operation

$\sum_{j=0}^{N}|j\rangle\langle\check{\phi}_{j}|.$ (12)

That is, it maps the state $|\check{\phi}_{j}\rangle$ to a computational basis
state $|j\rangle$, so a measurement in the computational basis gives the
result for the phase. This operation is the inverse of the usual quantum
Fourier transform, which would map from $|j\rangle$ to
$|\check{\phi}_{j}\rangle$.

If one aims to obtain the original POVM, one can randomly (with equal
probability) select $\delta\phi\in[0,2\pi/(N+1)]$, and choose the states with
$\check{\phi}_{j}=2\pi j/(N+1)+\delta\phi$. Then perform a measurement in the
basis $|\check{\phi}_{j}\rangle$ with this randomly chosen offset. The
complete measurement, including the random choice of $\delta\phi$, is then
equivalent to the POVM with the set of outcomes over a continuous range of
$\check{\phi}$. This approach can be used in order to give a measurement that
is covariant (has an error distribution independent of the system phase
$\phi$). In practice it is not usually needed, so we will not consider it
further in this paper.

In order to obtain the estimate for the phase, one should therefore perform
the inverse quantum Fourier transform on the control qubits. The inverse
quantum Fourier transform can be performed in a semiclassical way, by
performing measurements on successive qubits followed by controlled rotations
Griffiths and Niu (1996). The usual terminology is the ‘semiclassical Fourier
transform’, though this is the inverse transform. An example with three qubits
is given in Figure 3. The bottom (least significant qubit) is measured first.
The result is used to control a phase rotation on the middle qubit. Then the
middle qubit is measured, and the results of both measurements are used to
control phase rotations on the top qubit. The net result is the same as
performing the inverse quantum Fourier transform and measuring in the
computational basis.

@R=2em @C=1em & R(π4) R(π2) H
R(π2) [-1] H [-1]
H [-1]

Figure 3:  An example of the semiclassical Fourier transform on three qubits,
where the bottom qubit corresponds to the least significant bit.

A further advantage of this procedure is that the fact that the controlled $U$
operations are also performed in sequence means that the sequences can be
matched. That is, we have the combined procedure as shown in Figure 4. In the
case where control registers are prepared in an equal superposition state,
then they are unentangled. This means that preparation of each successive
qubit can be delayed until it is needed, as shown in Figure 5.

@R=2em @C=1em & 3 R H
2 R(π2) H [-1]
1 H [-1] [-1]
—ϕ⟩ U^4 U^2 U

Figure 4:  The combined procedure with controlled $U^{2^{j}}$ operations to
give phase kickback, together with the semiclassical Fourier transform. The
final controlled phase rotation is just shown as $R$, because the angle of
rotation is controlled by the combined results of the first two measurements.

@R=2em @C=1em & —0⟩ H 3 R H
—0⟩ H 2 R(π2) H [-1]
—0⟩ H 1 H [-1] [-1]
—ϕ⟩ U^4 U^2 U

Figure 5:  The combined procedure with preparation of control ancillas delayed
until they are needed.

What this means is that only one control qubit need be used at once. The
preparation of the next control qubit can be delayed until after measurement
of the previous one, and that qubit can be reset and reused. That is useful in
quantum algorithms with a limited number of qubits available, and is also
useful in quantum phase estimation. In that case, one can replace the control
qubits with NOON states with photon numbers that are powers of 2. Then these
NOON states can be measured in sequence to give a canonical measurement of
phase, even though a canonical measurement of phase would not be possible on a
single two-mode state. In Higgins  _et al._ (2007) we demonstrated this, using
multiple passes through a phase shift rather than NOON states.

The drawback now is that, even though it is possible to perform the canonical
measurement, a suboptimal state is being used. We would like to be able to
perform measurements achieving that minimum Holevo phase variance. In Higgins
_et al._ (2007) we showed that, by using multiple NOON states of each number
it is possible to obtain the desired scaling with total photon number, even
though there is a different constant factor so the true minimum error is not
achieved.

###  II.3 Performing phase measurement with two control qubits

Up until this point this section has been revision of prior work. What is new
here is that we show how to prepare the optimal state for phase measurement in
a sequential way, so the number of qubits that need be used at once is
minimised. We will show how the optimal state can be prepared using a sequence
of two-qubit operations, as in Figure 6.

@R=2em @C=1em —0⟩ & ⋯ 1U_1
—0⟩ ⋯ U_1
⋮
—0⟩ 1U_m-2 ⋯
—0⟩ 1U_m-1 U_m-2 ⋯
—0⟩ U_m-1 ⋯

Figure 6:  Preparation of the optimal state with a sequence of two-qubit
operations.

When the optimal state is prepared in this way, its preparation may be delayed
until the qubits are needed, as shown in Figure 7. This is illustrated with
three control qubits, where introduction of the third qubit can be delayed
until the first qubit is measured. In general, with more control qubits,
introduction of each additional qubit can be delayed until after measurement
of the qubit two places down, so only two control qubits need be used at once.

@R=2em @C=1em & —0⟩ 1U_1 3 R H
—0⟩ 1U_2 U_1 2 R(π2) H [-1]
—0⟩ U_2 1 H [-1] [-1]
—ϕ⟩ U^4 U^2 U

Figure 7:  The combined procedure with preparation of control ancillas in the
optimal state delayed until they are needed.

The reason why it is possible to prepare the optimal state in this way is that
it is a superposition of two unentangled states. The sine is a combination of
a positive and negative complex exponential as

$|\psi_{\rm
opt}\rangle=\frac{1}{2i}\sqrt{\frac{2}{N+2}}\sum_{n=0}^{N}\left[\exp\left(i\frac{\pi(n+1)}{N+2}\right)-\exp\left(-i\frac{\pi(n+1)}{N+2}\right)\right]|n\rangle.$
(13)

When $N=2^{m}-1$, we can write this as

$\displaystyle|\psi_{\rm opt}\rangle$
$\displaystyle=\frac{e^{i\pi/M}}{2i}\sqrt{\frac{2}{M}}\sum_{n_{1},\cdots,n_{m}=0}^{1}e^{i\pi\sum_{j=1}^{m}n_{j}2^{j-1}/M}|n_{1}\cdots
n_{m}\rangle-\frac{e^{-i\pi/M}}{2i}\sqrt{\frac{2}{M}}\sum_{n_{1},\cdots,n_{m}=0}^{1}e^{-i\pi\sum_{j=1}^{m}n_{j}2^{j-1}/M}|n_{1}\cdots
n_{m}\rangle$
$\displaystyle=\frac{e^{i\pi/M}}{2i}\sqrt{\frac{2}{M}}\sum_{n_{1},\cdots,n_{m}=0}^{1}\bigotimes_{j=0}^{m}e^{i\pi
n_{j}2^{j-1}/N}|n_{j}\rangle-\frac{e^{-i\pi/N}}{2i}\sqrt{\frac{2}{M}}\sum_{n_{1},\cdots,n_{m}=0}^{1}\bigotimes_{j=0}^{m}e^{-i\pi
n_{j}2^{j-1}/M}|n_{j}\rangle$
$\displaystyle=\frac{e^{i\pi/M}}{2i}\sqrt{\frac{2^{m+1}}{M}}\bigotimes_{j=0}^{m}\left(\sum_{n_{j}=0}^{1}\frac{e^{i\pi
n_{j}2^{j-1}/M}}{\sqrt{2}}|n_{j}\rangle\right)-\frac{e^{-i\pi/M}}{2i}\sqrt{\frac{2^{m+1}}{M}}\bigotimes_{j=0}^{m}\left(\sum_{n_{j}=0}^{1}\frac{e^{-i\pi
n_{j}2^{j-1}/M}}{\sqrt{2}}|n_{j}\rangle\right)$
$\displaystyle=\sqrt{\frac{2^{m-1}}{M}}\bigotimes_{j=0}^{m}\left(\sum_{n_{j}=0}^{1}\frac{e^{-i\pi(-1)^{n_{j}}2^{j-2}/M}}{\sqrt{2}}|n_{j}\rangle\right)+\sqrt{\frac{2^{m-1}}{M}}\bigotimes_{j=0}^{m}\left(\sum_{n_{j}=0}^{1}\frac{e^{i\pi(-1)^{n_{j}}2^{j-2}/M}}{\sqrt{2}}|n_{j}\rangle\right),$
(14)

where $M=N+2=2^{m}+1$. In the last line we have used

$\displaystyle\sum_{n_{j}=0}^{1}\frac{e^{i\pi n_{j}2^{j-1}/M}}{\sqrt{2}}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(|0\rangle+e^{i\pi
n_{j}2^{j-1}/M}\right)$ $\displaystyle=e^{i\pi
n_{j}2^{j-2}/M}\frac{1}{\sqrt{2}}\left(e^{-i\pi
n_{j}2^{j-2}/M}|0\rangle+e^{i\pi n_{j}2^{j-2}/M}|1\rangle\right)$ (15)

and then

$\displaystyle e^{i\pi/M}\prod_{j=0}^{m}e^{i\pi n_{j}2^{j-2}/M}$
$\displaystyle=e^{i\pi/M}e^{i\pi n_{j}\sum_{j=0}^{m}2^{j-2}/M}$
$\displaystyle=e^{i\pi/M}e^{i\pi n_{j}(2^{m}-1)/2M}$ $\displaystyle=e^{i\pi
n_{j}(2^{m}+1)/2M}$ $\displaystyle=e^{i\pi n_{j}/2}$ $\displaystyle=i.$ (16)

In order to write the optimal state in a more compact way we define the
following

$\displaystyle c$ $\displaystyle\coloneqq\sqrt{\frac{2^{m-1}}{M}},$
$\displaystyle|\phi_{j}^{+}\rangle$ $\displaystyle\coloneqq\frac{e^{i\pi
2^{j-2}/M}|0\rangle+e^{-i\pi 2^{j-2}/M}|1\rangle}{\sqrt{2}},$
$\displaystyle|\phi_{j}^{-}\rangle$ $\displaystyle\coloneqq\frac{e^{-i\pi
2^{j-2}/M}|0\rangle+e^{i\pi 2^{j-2}/M}|1\rangle}{\sqrt{2}}.$ (17)

Then the optimal state in Section II.3 can be written as

$|\psi_{\rm
opt}\rangle=c\bigotimes_{j=1}^{m}|\phi_{j}^{+}\rangle+c\bigotimes_{j=1}^{m}|\phi_{j}^{-}\rangle.$
(18)

That is, it is an equally weighted superposition of two states, which are each
unentangled between all qubits. What this means is that any bipartite split of
the state will have Schmidt number 2, so the entanglement across the bipartite
split can be represented on a single qubit on one side. We use that principle
in the state preparation. At any stage, after performing the two-qubit
operation between qubit $j$ and $j+1$, there will be the correct bipartite
entanglement in the split between qubits up to $j$ and qubits from $j+1$ to
$m$. However, at that stage qubits from $1$ to $j-1$ have not been initialised
yet, so the entanglement across the bipartite split (for qubits $1$ to $j$) is
represented just on qubit $j$.

##  III Recursive construction of the optimum state

In this section, we show how to create the optimal state Eq. 18 recursively.
We introduce the partial tensor product states

$\displaystyle|\phi_{[\ell]}^{+}\rangle$
$\displaystyle\coloneqq\bigotimes_{j=1}^{\ell}|\phi_{j}^{+}\rangle,\qquad|\phi_{[\ell]}^{-}\rangle\coloneqq\bigotimes_{j=1}^{\ell}|\phi_{j}^{-}\rangle.$
(19)

Because $|\psi_{\rm opt}\rangle$ is a linear combination of
$|\phi_{[m]}^{\pm}\rangle$, the state of qubits 1 to $\ell$ an be represented
as a linear combination of $|\phi_{[\ell]}^{\pm}\rangle$. In order to describe
the operations needed to prepare the state $|\psi_{\rm opt}\rangle$, we need
to describe the state of qubits 1 to $\ell$ in terms of orthogonal states,
which we will denote by $|\Phi_{[\ell]}^{\pm}\rangle$. These orthogonal (but
not normalised) states

$\displaystyle|\Phi_{[\ell]}^{+}\rangle$
$\displaystyle\coloneqq|\phi_{[\ell]}^{+}\rangle+|\phi_{[\ell]}^{-}\rangle,\qquad|\Phi_{[\ell]}^{-}\rangle\coloneqq|\phi_{[\ell]}^{+}\rangle-|\phi_{[\ell]}^{-}\rangle.$
(20)

It is possible to prove that these states are orthogonal as in the following
Lemma.

######  Lemma 1.

The states $|\Phi_{[\ell]}^{+}\rangle$ and $|\Phi_{[\ell]}^{-}\rangle$ defined
in Eq. 20, are orthogonal:

$\langle\Phi_{[\ell]}^{-}|\Phi_{[\ell]}^{+}\rangle=0.$ (21)

###### Proof.

From Section II.3 we have

$\displaystyle\langle\phi_{j}^{-}|\phi_{j}^{+}\rangle$ $\displaystyle=\cos(\pi
2^{j-1}/M),$ (22)

which is real. In turn, that implies
$\langle\phi_{[\ell]}^{-}|\phi_{[\ell]}^{+}\rangle$ is real, and equal to
$\langle\phi_{[\ell]}^{+}|\phi_{[\ell]}^{-}\rangle$. Moreover, because
$|\Phi_{j}^{\pm}\rangle$ are normalised, so are $|\Phi_{[\ell]}^{\pm}\rangle$.
Therefore we obtain

$\displaystyle\langle\Phi_{[\ell]}^{-}|\Phi_{[\ell]}^{+}\rangle$
$\displaystyle=\langle\phi_{[\ell]}^{+}|\phi_{[\ell]}^{+}\rangle-\langle\phi_{[\ell]}^{-}|\phi_{[\ell]}^{-}\rangle+\langle\phi_{[\ell]}^{+}|\phi_{[\ell]}^{-}\rangle-\langle\phi_{[\ell]}^{-}|\phi_{[\ell]}^{+}\rangle=0$
(23)

Here, the first two terms cancel because they are both 1 (due to
normalisation) and the second two terms cancel because they are equal (because
they are real). ∎

Next we wish to show that there is a simple recurrence relation for the states
$|\Phi_{[\ell]}^{\pm}\rangle$, in their normalised form

$|\widetilde{\Phi}_{[\ell]}^{\pm}\rangle\coloneqq\frac{|\Phi_{[\ell]}^{\pm}\rangle}{\sqrt{\langle\Phi_{[\ell]}^{\pm}|\Phi_{[\ell]}^{\pm}\rangle}}.$
(24)

We will use the recurrence relation to derive the sequence of two-qubit
operations to prepare the initial state. The result is as follows.

######  Lemma 2.

The states $|\widetilde{\Phi}_{[\ell]}^{+}\rangle$ and
$|\widetilde{\Phi}_{[\ell]}^{-}\rangle$ defined in Eq. 20 and Eq. 24, have
recurrence relation

$|\widetilde{\Phi}_{[\ell+1]}^{\pm}\rangle=\mu_{\ell}^{(0,\pm)}|\widetilde{\Phi}_{[\ell]}^{\pm}\rangle|+\rangle+\mu_{\ell}^{(1,\pm)}|\widetilde{\Phi}_{[\ell]}^{\mp}\rangle|-\rangle,$
(25)

where

$\displaystyle\mu_{\ell}^{(0,\pm)}$
$\displaystyle=\cos(2^{\ell-1}\pi/M)P^{0,\pm},$ (26)
$\displaystyle\mu_{\ell}^{(1,\pm)}$
$\displaystyle=i\sin(2^{\ell-1}\pi/M)P^{1,\pm},$ (27) $\displaystyle
P^{s,\pm}$
$\displaystyle=\sqrt{\frac{1\pm(-1)^{s}\prod_{j=1}^{\ell}\cos(2^{j-1}\pi/M)}{1\pm\prod_{j=1}^{\ell+1}\cos(2^{j-1}\pi/M)}}.$
(28)

###### Proof.

To prove this, we start by noting that

$\displaystyle\cos(2^{j-2}\pi/M)$
$\displaystyle=\langle+|\phi_{j}^{\pm}\rangle,$ $\displaystyle
i\sin(2^{j-2}\pi/M)$ $\displaystyle=\langle-|\phi_{j}^{+}\rangle,$
$\displaystyle-i\sin(2^{j-2}\pi/M)$
$\displaystyle=\langle-|\phi_{j}^{-}\rangle.$ (29)

Therefore we can see that

$\displaystyle\langle+|\Phi_{[\ell+1]}^{+}\rangle$
$\displaystyle=\cos(2^{\ell-1}\pi/M)|\Phi_{[\ell]}^{+}\rangle,$
$\displaystyle\langle-|\Phi_{[\ell+1]}^{+}\rangle$
$\displaystyle=i\sin(2^{\ell-1}\pi/M)|\Phi_{[\ell]}^{-}\rangle,$ (30)

which implies

$|\Phi_{[\ell+1]}^{+}\rangle=\cos(2^{\ell-1}\pi/M)|\Phi_{[\ell]}^{+}\rangle|+\rangle+i\sin(2^{\ell-1}\pi/M)|\Phi_{[\ell]}^{-}\rangle|-\rangle.$
(31)

Similarly, we find

$\displaystyle\langle+|\Phi_{[\ell+1]}^{-}\rangle$
$\displaystyle=\cos(2^{\ell-1}\pi/M)|\Phi_{[\ell]}^{-}\rangle,$
$\displaystyle\langle-|\Phi_{[\ell+1]}^{-}\rangle$
$\displaystyle=i\sin(2^{\ell-1}\pi/M)|\Phi_{[\ell]}^{+}\rangle,$ (32)

which implies

$|\Phi_{[\ell+1]}^{-}\rangle=\cos(2^{\ell-1}\pi/M)|\Phi_{[\ell]}^{-}\rangle|+\rangle+i\sin(2^{\ell-1}\pi/M)|\Phi_{[\ell]}^{+}\rangle|-\rangle.$
(33)

This gives us recurrence relations for $|\Phi_{[\ell]}^{\pm}\rangle$, which
can be written

$|\Phi_{[\ell+1]}^{\pm}\rangle=\cos(2^{\ell-1}\pi/M)|\Phi_{[\ell]}^{\pm}\rangle|+\rangle+i\sin(2^{\ell-1}\pi/M)|\Phi_{[\ell]}^{\mp}\rangle|-\rangle.$
(34)

Let us define the normalisation

${\mathcal{N}}_{[\ell]}^{\pm}=\sqrt{\langle\Phi_{[\ell]}^{\pm}|\Phi_{[\ell]}^{\pm}\rangle}.$
(35)

In terms of this, the recurrence relation for the normalised states is

$|\widetilde{\Phi}_{[\ell+1]}^{\pm}\rangle=\cos(2^{\ell-1}\pi/M)\frac{{\mathcal{N}}_{[\ell]}^{+}}{{\mathcal{N}}_{[\ell+1]}^{\pm}}|\Phi_{[\ell]}^{\pm}\rangle|+\rangle+i\sin(2^{\ell-1}\pi/M)\frac{{\mathcal{N}}_{[\ell]}^{-}}{{\mathcal{N}}_{[\ell+1]}^{\pm}}|\widetilde{\Phi}_{[\ell]}^{\mp}\rangle|-\rangle.$
(36)

The normalisation can be determined using

$\langle\phi_{[\ell]}^{-}|\phi_{[\ell]}^{+}\rangle=\prod_{j=1}^{\ell}\cos(\pi
2^{j-1}/M),$ (37)

which gives

$\displaystyle({\mathcal{N}}_{[\ell]}^{\pm})^{2}=\langle\Phi_{[\ell]}^{\pm}|\Phi_{[\ell]}^{\pm}\rangle=2\pm
2\prod_{j=1}^{\ell}\cos(\pi 2^{j-1}/M).$ (38)

That gives us the ratios of norms

$\frac{{\mathcal{N}}_{[\ell]}^{+}}{{\mathcal{N}}_{[\ell+1]}^{\pm}}=P^{0,\pm},\qquad\frac{{\mathcal{N}}_{[\ell]}^{-}}{{\mathcal{N}}_{[\ell+1]}^{\pm}}=P^{1,\pm}.$
(39)

Hence Eq. 37 is the form of the recurrence relation required. ∎

##  IV Preparing optimum state with two-qubit unitaries

In the previous section, we showed that it is possible to construct an
orthonormal basis for the state on qubits 1 to $\ell$ as
$|\widetilde{\Phi}_{[\ell]}^{\pm}\rangle$, which satisfy the recursive
relation

$\displaystyle|\widetilde{\Phi}_{[\ell+1]}^{+}\rangle$
$\displaystyle=\mu_{\ell}^{(0,+)}|\widetilde{\Phi}_{[\ell]}^{+}\rangle|+\rangle+\mu_{\ell}^{(1,+)}|\widetilde{\Phi}_{[\ell]}^{-}\rangle|-\rangle,$
$\displaystyle|\widetilde{\Phi}_{[\ell+1]}^{-}\rangle$
$\displaystyle=\mu_{\ell}^{(0,-)}|\widetilde{\Phi}_{[\ell]}^{-}\rangle|+\rangle+\mu_{\ell}^{(1,-)}|\widetilde{\Phi}_{[\ell]}^{+}\rangle|-\rangle\,,$
(40)

for constants $\mu_{\ell}^{(s,\pm)}$. Moreover, it is obvious from the
definitions that $|\widetilde{\Phi}_{[m]}^{+}\rangle=|\psi_{\rm opt}\rangle$.
Therefore, the optimum state may be constructed in a recursive way from
$|\widetilde{\Phi}_{[m-1]}^{\pm}\rangle$ on qubits 1 to $m-1$, which can be
constructed from states $|\widetilde{\Phi}_{[m-2]}^{\pm}\rangle$ on qubits 1
to $m-2$, and so forth.

In order to prepare the optimal state, we can apply a stepwise procedure where
the principle is to use a single qubit to flag which of
$|\Phi_{[\ell]}^{\pm}\rangle$ is to be prepared on the remaining qubits 1 to
$\ell$. We start from qubits $m$ and $m-1$ and work back to qubits 1 and 2. It
is convenient to describe the operations as acting on qubits initialised as
$|+\rangle$. Then we initially perform an operation on qubits $m$ and $m-1$
that maps

$U_{m-1}|+\rangle|+\rangle=\mu_{m-1}^{(0,+)}|+\rangle|+\rangle+\mu_{m-1}^{(1,+)}|-\rangle|-\rangle.$
(41)

The principle of this operation is that it corresponds to the recursion

$|\psi_{\rm
opt}\rangle=|\widetilde{\Phi}_{[m]}^{+}\rangle=\mu_{m-1}^{(1)}|\widetilde{\Phi}_{[m-1]}^{+}\rangle|+\rangle+\mu_{m-1}^{(2)}|\widetilde{\Phi}_{[m-1]}^{-}\rangle|-\rangle\,.$
(42)

At this stage we only have qubits, so the states
$|\widetilde{\Phi}_{[m-1]}^{\pm}\rangle$ on qubits 1 to $m-1$ need to be
represented by $|\pm\rangle$ on a single qubit.

It is trivial to see that unitary $U_{m-1}$ exists; it can explicitly be
performed by rotating qubit $m$ as

$|+\rangle\mapsto\mu_{m-1}^{(0,+)}|+\rangle+\mu_{m-1}^{(1,+)}|-\rangle.$ (43)

We can alternatively describe the operation as having the matrix form in the
$|\pm\rangle$ basis

$U_{m-1}\equiv\begin{pmatrix}\mu_{m-1}^{(0,+)}&*&*&*\\\ 0&*&*&*\\\ 0&*&*&*\\\
\mu_{m-1}^{(1,+)}&*&*&*\end{pmatrix},$ (44)

where $*$ indicates entries where the value is unimportant.

We then perform $U_{m-2}$ on qubits $m-1$ and $m-2$, down to $U_{1}$ on qubits
1 and 2. The unitary $U_{\ell}$ needs to map

$\displaystyle U_{\ell}|+\rangle|+\rangle$
$\displaystyle=\mu_{\ell}^{(0,+)}|+\rangle|+\rangle+\mu_{\ell}^{(1,-)}|-\rangle|-\rangle$
$\displaystyle U_{\ell}|+\rangle|-\rangle$
$\displaystyle=\nu_{\ell}^{(0,+)}|-\rangle|+\rangle+\nu_{\ell}^{(1,-)}|+\rangle|-\rangle.$
(45)

This corresponds to the recursion given in Section IV, and the states
$|\pm\rangle$ on qubit $\ell$ are being used to represent
$|\widetilde{\Phi}_{[\ell]}^{\pm}\rangle$ on qubits 1 to $\ell$. This
operation has the matrix entries

$U_{\ell}\equiv\begin{pmatrix}\mu_{\ell}^{(0,+)}&0&*&*\\\
0&\mu_{\ell}^{(1,-)}&*&*\\\ 0&\mu_{\ell}^{(0,-)}&*&*\\\
\mu_{\ell}^{(1,+)}&0&*&*\end{pmatrix}.$ (46)

This operation may be achieved in the following way. Define the single-qubit
rotations $V_{\ell}^{(0,+)}$ and $V_{\ell}^{(1)}$ to act as

$\displaystyle V_{\ell}^{(0,+)}|+\rangle$
$\displaystyle=\mu_{\ell}^{(0,+)}|+\rangle+\mu_{\ell}^{(1,+)}|-\rangle$
$\displaystyle V_{\ell}^{(0,-)}|+\rangle$
$\displaystyle=\mu_{\ell}^{(0,-)}|-\rangle+\mu_{\ell}^{(1,-)}|+\rangle\,.$
(47)

Then $U_{\ell}$ corresponds to the controlled operation

$U_{\ell}=V_{\ell}^{(0,+)}\otimes|+\rangle\bra{+}+V_{\ell}^{(0,-)}\otimes|-\rangle\bra{-}.$
(48)

This method could be used for $U_{m-1}$, though the method described above is
simpler.

After performing this sequence of unitaries, we then need to map $|\pm\rangle$
to $|\widetilde{\Phi}_{[1]}^{\pm}\rangle$ on qubit 1. This is a simple single-
qubit unitary operation, which can be combined with $U_{1}$ to give the
correct final state with a sequence of two-qubit unitary operations. Thus our
recursive expression for the states gives us a sequence of two-qubit unitaries
to create the optimal state.

To be more specific about what operation is needed,

$|\phi_{[1]}^{+}\rangle=\frac{e^{i\pi/2M}|0\rangle+e^{-i\pi/2M}|1\rangle}{\sqrt{2}},$
(49)

so that

$\displaystyle|\Phi_{[1]}^{+}\rangle$
$\displaystyle=\sqrt{2}\left(\cos(\pi/2M)|0\rangle+\cos(\pi/2M)|1\rangle\right),$
(50) $\displaystyle|\Phi_{[1]}^{-}\rangle$
$\displaystyle=i\sqrt{2}\left(\sin(\pi/2M)|0\rangle-\sin(\pi/2M)|1\rangle\right).$
(51)

That gives the normalised states

$\displaystyle|\widetilde{\Phi}_{[1]}^{+}\rangle$
$\displaystyle=|+\rangle,\qquad|\widetilde{\Phi}_{[1]}^{-}\rangle=i|-\rangle.$
(52)

Therefore the operation needed is an $i$ phase shift on $|-\rangle$.

##  V Conclusion

We have shown how to create the optimal state for phase estimation, in the
sense of minimising Holevo variance, using a sequence of two qubit operations.
When combining this sequential process with the semiclassical quantum Fourier
transform, we can entangle new qubits after measuring control qubits in such a
way that only two control qubits are needed at once. This means that the qubit
that is measured can be reset and used as the new qubit to be entangled,
minimising the need for ancilla qubits.

In quantum algorithms where phase estimation is needed with a small number of
logical qubits this is ideal. Previously the method used was either many
entangled qubits, increasing the size of the quantum computer needed, or a
single control qubit, which significantly increases the error. In our method
the number of control qubits is only increased by 1, while giving the minimal
error. Here the quantity being exactly minimised is the Holevo variance, which
is very close to the mean-square error (MSE) for sharply peaked distributions.
If one were interested in minimising MSE, then these states give the same
leading order term for MSE as the minimum MSE Berry  _et al._ (2012), so these
states are still suitable.

Our method of preparing the state, although it has been derived for the
specific case of the optimal state for minimising Holevo variance, could also
be applied to other states that are a superposition of two unentangled states.
The crucial feature is that the Schmidt number is 2 for any bipartite split
across the qubits. One could also consider states with larger Schmidt number,
and use a larger number of qubits as controls. That could potentially be used
for states that are optimal for minimising other measures of error. For
example, one could consider methods of approximating Kaiser windows or the
digital prolate spheroidal sequence, as is suitable for optimising confidence
intervals Berry  _et al._ (2022).

Another interesting question is whether this procedure could be demonstrated
with photons. A scheme with optimal phase states for $N=2$ using two photons
was demonstrated in Daryanoosh  _et al._ (2018). With the preparation scheme
we have outlined, it would potentially be possible to demonstrate these states
with larger $N$, though it would require entangling operations that might
require nonlinear optical elements.

###### Acknowledgements.

DWB worked on this project under a sponsored research agreement with Google
Quantum AI. DWB is also supported by Australian Research Council Discovery
Projects DP190102633, DP210101367, and DP220101602.

## Author Declarations

### Conflict of interest

The authors have no conflicts to disclose.

## Data availability

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

## References

  * Shor (1994) P. Shor, in  _Proceedings 35th Annual Symposium on Foundations of Computer Science_ (1994) pp. 124–134.
  * Aspuru-Guzik  _et al._ (2005) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704 (2005).
  * Griffiths and Niu (1996) R. B. Griffiths and C.-S. Niu, Physical Review Letters 76, 3228 (1996).
  * Babbush  _et al._ (2018) R. Babbush, C. Gidney, D. W. Berry, N. Wiebe, J. McClean, A. Paler, A. Fowler, and H. Neven, Physical Review X 8, 041015 (2018).
  * Luis and Peřina (1996) A. Luis and J. Peřina, Physical Review A 54, 4564 (1996).
  * Berry and Wiseman (2000) D. W. Berry and H. M. Wiseman, Physical Review Letters 85, 5098 (2000).
  * Boto  _et al._ (2000) A. N. Boto, P. Kok, D. S. Abrams, S. L. Braunstein, C. P. Williams, and J. P. Dowling, Physical Review Letters 85, 2733 (2000).
  * Kok  _et al._ (2002) P. Kok, H. Lee, and J. P. Dowling, Physical Review A 65, 052104 (2002).
  * Lee  _et al._ (2002) H. Lee, P. Kok, and J. P. Dowling, Journal of Modern Optics 49, 2325 (2002).
  * Dinani and Berry (2014) H. T. Dinani and D. W. Berry, Physical Review A 90, 023856 (2014).
  * Dinani  _et al._ (2016) H. T. Dinani, M. K. Gupta, J. P. Dowling, and D. W. Berry, Physical Review A 93, 063804 (2016).
  * Huang  _et al._ (2017) Z. Huang, K. R. Motes, P. M. Anisimov, J. P. Dowling, and D. W. Berry, Physical Review A 95, 053837 (2017).
  * Motes  _et al._ (2016) K. R. Motes, R. L. Mann, J. P. Olson, N. M. Studer, E. A. Bergeron, A. Gilchrist, J. P. Dowling, D. W. Berry, and P. P. Rohde, Physical Review A 94, 012344 (2016).
  * Huang  _et al._ (2021) Z. Huang, P. P. Rohde, D. W. Berry, P. Kok, J. P. Dowling, and C. Lupo, Quantum 5, 447 (2021).
  * Higgins  _et al._ (2007) B. L. Higgins, D. W. Berry, S. D. Bartlett, H. M. Wiseman, and G. J. Pryde, Nature 450, 393 (2007).
  * Kivlichan  _et al._ (2020) I. D. Kivlichan, C. Gidney, D. W. Berry, N. Wiebe, J. McClean, W. Sun, Z. Jiang, N. Rubin, A. Fowler, A. Aspuru-Guzik, H. Neven, and R. Babbush, Quantum 4, 296 (2020).
  * Sanders  _et al._ (1997) B. C. Sanders, G. J. Milburn, and Z. Zhang, Journal of Modern Optics 44, 1309 (1997).
  * Holevo (1984) A. S. Holevo, in  _Quantum Probability and Applications to the Quantum Theory of Irreversible Processes_ , edited by L. Accardi, A. Frigerio, and V. Gorini (Springer Berlin Heidelberg, Berlin, Heidelberg, 1984) pp. 153–172.
  * Berry  _et al._ (2012) D. W. Berry, M. J. W. Hall, M. Zwierz, and H. M. Wiseman, Physical Review A 86, 053813 (2012).
  * Berry  _et al._ (2022) D. W. Berry, Y. Su, C. Gyurik, R. King, J. Basso, A. D. T. Barba, A. Rajput, N. Wiebe, V. Dunjko, and R. Babbush, arXiv: 2209.13581 (2022).
  * Daryanoosh  _et al._ (2018) S. Daryanoosh, S. Slussarenko, D. W. Berry, H. M. Wiseman, and G. J. Pryde, Nature Communications 9, 4606 (2018).