# Continuous-Variable Quantum Key Distribution with a Real Local Oscillator
and without Auxiliary Signals

Sebastian Kleis, Max Rueckmann, Christian G. Schaeffer<EMAIL_ADDRESS>Helmut-Schmidt-Universität Hamburg

(August 27, 2024)

###### Abstract

Continuous-variable quantum key distribution (CV-QKD) is realized with
coherent detection and is therefore very suitable for a cost-efficient
implementation. The major challenge in CV-QKD is mitigation of laser phase
noise at a signal to noise ratio of much less than
$0\text{\,}\mathrm{dB}\text{/}$. So far, this has been achieved with a remote
local oscillator or with auxiliary signals. For the first time, we
experimentally demonstrate that CV-QKD can be performed with a real local
oscillator and without auxiliary signals which is achieved by applying Machine
Learning methods. It is shown that, with the most established discrete
modulation protocol, the experimental system works down to a quantum channel
signal to noise ratio of $-19.1\text{\,}\mathrm{dB}\text{/}$. The performance
of the experimental system allows CV-QKD at a key rate of
$9.2\text{\,}\mathrm{Mbit}\text{/}\mathrm{s}$ over a fiber distance of
$26\text{\,}\mathrm{km}\text{/}$. After remote local oscillator and auxiliary
signal aided CV-QKD, this could mark a starting point for a third generation
of CV-QKD systems that are even more attractive for a wide implementation
because they are almost identical to standard coherent systems.

##  I Introduction

With quantum computers threatening to break the security of today’s
cryptosystems, the field of quantum communications attracts increasing
attention. In this field, continuous variable quantum key distribution (CV-
QKD) is very attractive for a practical implementation because it is based on
coherent detection and promises to require only commercial off-the-shelf
components. However, one crucial difference to classical optical
communications is that the received power level is usually less than
$1\text{\,}\mathrm{ph}\text{/}\mathrm{bit}$ for a discrete modulation
[leverrier2011], translating to a signal to noise ratio (SNR) of less than
$0\text{\,}\mathrm{dB}\text{/}$. Additionally, residual phase noise
contributes to the excess noise [kleis2017], which is the critical performance
parameter for CV-QKD [grosshans2003]. The combination of ultra low SNR and
required accuracy makes carrier phase recovery the major challenge in
practical CV-QKD. The first approach to deal with this challenge was to avoid
phase noise by using a remote local oscillator (LO) [grosshans2003,
lodewyck2005, qi2007, jouguet2013]. The disadvantage is that the remote LO
compromises the security and limits the achievable distance [jouguet2013a].
Therefore, it is highly preferable that CV-QKD systems work with a real LO
generated by a separate laser source at the receiver site (Bob). Since 2015,
various real LO systems have been proposed and experimentally demonstrated
[kleis2015, huang2015, qi2015, soh2015, schrenk2016, marie2017, schrenk2017,
comandar2017, laudenbach2017, wang2018]. In all the proposed systems, carrier
phase estimation is based on auxiliary signals, also called pilots, that are
generated at the transmitter site (Alice). However, the pilots occupy
additional bandwidth and significantly increase the complexity of the system.
For example, the system in [laudenbach2017] uses only one polarization for the
quantum signal and the orthogonal one for a pilot tone. In [kleis2015], two
pilot tones are multiplexed with the quantum signal in the frequency domain
and occupy more bandwidth than the quantum signal itself. In both cases, the
spectral efficiency could at least be doubled without pilots. Another
important issue is that pilot tones are not included in the current security
proofs for CV-QKD [fossier2009, leverrier2011]. Despite all these issues, no
real LO system without pilot tones has been reported up to now. Regarding the
large number of pilot-based systems proposed, it appears as if pilot tones
were indispensable in CV-QKD systems with a real LO. However, no convincing
argument or experimental demonstration about the necessity of pilots has been
made yet. The present article addresses the question of whether or not
designing a CV-QKD system without pilot tones is possible with current fiber-
optic technology. Our approach is to design and experimentally investigate
such a CV-QKD system. In this, the quantum signal is discrete phase modulated
with an order of $M=2$ or $M=4$. The latter corresponds to the most
established CV-QKD protocol with a discrete modulation. A Bayesian particle
smoother, which is trained using Monte Carlo Markov chain (MCMC) methods, is
used for carrier phase estimation at the receiver. As particle smoothing
achieves optimum tracking of dynamic variables it is very suitable to
investigate the limits of carrier phase estimation in the ultra-low SNR
regime. In any CV-QKD protocol, it is absolutely necessary that Alice reveals
a fraction of her transmitted symbols via a public classical channel
[grosshans2003a]. We demonstrate that these already available symbols can also
be used to substantially improve the phase noise mitigation which enables
pilotless CV-QKD without any loss of efficiency.

##  II Results

Figure 1: Experimental setup and digital signal processing (DSP) routine for
the CV-QKD system without pilot tones. The quantum signal is discrete phase-
modulated with a modulation order $M=2$ or $M=4$ at a symbol rate of
$17\text{\,}\mathrm{GBd}\text{/}$. The received and digitized signal is
recorded for offline processing. For the DSP routine, the total signal is
split into blocks with index $n$. Each block contains
$L=$6.8\text{\times}{10}^{5}$$ symbols. Publicly revealing a randomly selected
part of Alice’s symbols $a_{k}$ is necessary in CV-QKD [lodewyck2005]. Carrier
phase estimation is performed by a particle smoother that takes into account
these available symbols.

The experimental setup and digital signal processing (DSP) routine is shown in
figure 1. The quantum signal is discrete phase-modulated in baseband with a
modulation order $M=2$ or $M=4$ at a symbol rate of
$17\text{\,}\mathrm{GBd}\text{/}$. Bob uses a balanced receiver to perform
heterodyne detection at an intermediate frequency of
$\delta\nu_{\mathrm{AB}}\approx$10\text{\,}\mathrm{GHz}\text{/}$$. In the DSP
routine, Bob first down-converts the signal to baseband and applies the
matched filter $h_{\mathrm{Q}}$. Then, the signal is down-sampled to 2 samples
per symbol (sps) and equalized by the FIR-filter $h_{\mathrm{Eq}}$, which is
tuned using the constant-modulus algorithm only once for a high SNR. Timing
synchronization is performed by the digital filter and square algorithm
[oerder1988] followed by an extended Kalman filter that tracks the argument of
the complex Fourier coefficient $X_{n}$, where $n$ is the block index of the
blockwise DSP procedure. After timing correction ($h_{\mathrm{shift}}$) and
down-sampling to one sample per symbol, carrier phase estimation is performed
by a particle smoother. In order to find the optimum parameters $\hat{\theta}$
of the state space model, an extended Kalman filter is trained using MCMC
methods [sarkka2013].

The achievable key rate and distance of any CV-QKD system is very sensitive to
the excess noise power $\xi_{b}$ normalized to shot noise units because it
could provide information to an eavesdropper. To obtain an accurate excess
noise estimate after the quantum communication, receiver noise calibration and
quantum signal power estimation are required in CV-QKD. Therefore, excess
noise estimation is also an important part of our experimental investigation.
One can write the mean power of the received symbols $b_{k}$ as

$\displaystyle P_{b}=P_{\mathrm{Q}}+P_{\mathrm{SN}}+P_{\mathrm{EN}}+\xi_{b},$
(1)

where $P_{\mathrm{Q}}$, $P_{\mathrm{SN}}$, $P_{\mathrm{EN}}$ and $\xi_{b}$ are
the quantum signal power, the shot noise power, the electrical receiver noise
power and the excess noise power in arbitrary units respectively. The excess
noise in shot noise units is calculated as

$\displaystyle\xi_{b}^{\prime}=2\frac{P_{b}-P_{\mathrm{Q}}-P_{\mathrm{SN}}-P_{\mathrm{EN}}}{P_{\mathrm{SN}}}.$
(2)

To calibrate $P_{\mathrm{SN}}$ and $P_{\mathrm{EN}}$, the quantum signal is
deactivated. The total noise power $P_{\mathrm{SN}}+P_{\mathrm{EN}}$ and the
electrical noise power $P_{\mathrm{EN}}$ are calibrated in separate
measurements with activated and deactivated LO respectively. This is done
after each quantum signal measurement. To estimate $P_{\mathrm{Q}}$ in our
experiments, we evaluate the correlation between $a_{k}$ and $b_{k}$.

Generally in CV-QKD, to estimate $P_{\mathrm{Q}}$ and other important
parameters such as the channel transmission and the mutual information shared
between Alice and Bob, Alice has to reveal a randomly selected subset of her
transmitted symbols $a_{k}$ after the quantum signal transmission
[grosshans2003a]. More precisely, she publicly reveals the vector $a_{k}\cdot
r_{k}$, where

$\displaystyle r_{k}=\begin{cases}1,&\text{with probability }p_{r}\\\\[0.0pt]
0,&\text{with probability }1-p_{r}.\end{cases}$ (3)

That means that for each individual quantum symbol, there is a probability
$p_{r}$ of being revealed after the quantum communication. As a side effect of
symbol revelation, the particle smoother can take advantage of the revealed
symbols to improve the accuracy of phase estimation. Depending on the value of
each element $r_{k}$, the particle smoother adapts its measurement model
accordingly. This means that phase estimation can only take place after
Alice’s symbol revelation. However, as the sequence $\tilde{b}_{k}$ is already
down-sampled to one sample per symbol, it is no additional effort to store the
sequence $\tilde{b}_{k}$ compared to $b_{k}$. And since the continuous
variable $b_{k}$ must be stored for error correction anyway, which takes place
after Alice’s symbols have been revealed [grosshans2003], this approach does
not affect the feasibility nor the efficiency of the CV-QKD system in any way.

The noise calibration not only enables the excess noise estimation but also
the estimation of the electrical to shot noise ratio
$\varepsilon_{\mathrm{el}}$ and the receiver efficiency $\eta$. These are
important characteristics as they not only have an impact on the secret key
rate but also on the relation between received optical power $N_{\mathrm{B}}$
in photons per symbol and the SNR of $b_{k}$ as

$\displaystyle
N_{\mathrm{B}}=\mathrm{SNR}_{b}(1+\varepsilon_{\mathrm{el}})/\eta.$ (4)

In the present setup, $\varepsilon_{\mathrm{el}}\approx$0.70$$ and
$\eta\approx$0.232$$ which results in
$N_{\mathrm{B}}\approx$7.33$\,\mathrm{SNR}_{b}$ corresponding to a penalty for
$\mathrm{SNR}_{b}$ of $8.65\text{\,}\mathrm{dB}\text{/}$ compared to an ideal
heterodyne receiver with $\varepsilon_{\mathrm{el}}=0$ and $\eta=1$.

To initialize the quantum communication, all symbols of the first signal block
are revealed by Alice. This helps the particle smoother to obtain an accurate
initial estimate $\hat{f}_{\mathrm{Q}}$ of the quantum channel frequency
$f_{\mathrm{Q}}$ even at ultra low SNR. After the first block, the down-
conversion frequency is updated with $\hat{f}_{\mathrm{Q}}$ and the quantum
communication begins.

A residual phase noise in $b_{k}$ induces excess noise that is proportional to
the received quantum signal power [kleis2019]. This means that with decreasing
$\mathrm{SNR}_{b}$ the impact of phase errors on the excess noise is reduced.
At ultra low $\mathrm{SNR}_{b}$, the excess noise can be relatively small,
even if the carrier phase estimation fails completely. However, in this case
the mutual information between Alice and Bob drops drastically which prevents
a successful key generation. Therefore, in order to verify a successful signal
demodulation at ultra low $\mathrm{SNR}_{b}$, we also evaluate the hard
decision mutual information $I_{\mathrm{AB}}$ and compare it to its
theoretical value.

Figure 2: Block diagram for the simulations. The measured phase was obtained
using the experimental setup where the SNR of the quantum signal was
$7.1\text{\,}\mathrm{dB}\text{/}$ and $p_{r}=1$. The transmission model
corresponds to an ideal AWGN channel. For phase estimation, the same particle
smoother as in the experiments is used. All other impairments that are present
in the experiments such as timing errors, chromatic dispersion, bandwidth
limitations and other linear distortions are excluded from the simulations.

In addition to the experiments, we carried out simulations to isolate the
impact of phase noise and to investigate the phase noise limited performance
of our system. A block diagram of the simulation model is shown in figure 2.
As only additive noise and laser phase noise are included in the simulation
model, it does indeed correspond to a phase noise limited system. The phase
noise sequence used in the simulations is taken from a measurement with the
experimental setup. Thus, the performance of the simulation model corresponds
to the best that could be achieved with the laser sources in use.

$-20$$-15$$-10$$-5$$0$$5$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$I_{\mathrm{AB}}$
[$\mathrm{bit}\text{/}\mathrm{sym}$]$p_{r}=$0$$$p_{r}=$0.05$$$p_{r}=$0.05$$,
fit$p_{r}=$1$$theory$-20$$-15$$-10$$-5$$0$$5$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$\mathrm{SNR}_{b}$
[$\mathrm{dB}\text{/}$]$\hat{\xi}_{b}^{\prime}$ [$\mathrm{SNU}\text{/}$] (a)
Experiments with $M=4$

$-35$$-30$$-25$$-20$$-15$$-10$$-5$$0$$5$$10$$-35$$-30$$-25$$-20$$-15$$-10$$-5$$0$$5$$10$$\mathrm{SNR}_{b}$
[$\mathrm{dB}\text{/}$] (b) Simulations with $M=4$

Figure 3: Resulting system performance for $M=4$ in terms of the mutual
information between Alice and Bob $I_{\mathrm{AB}}$ and the excess noise in
Bob’s received symbols $\xi_{b}^{\prime}$. Experimental results are shown on
the left hand side (a). Simulation results are shown on the right hand side
(b).

$-20$$-15$$-10$$-5$$0$$5$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$I_{\mathrm{AB}}$
[$\mathrm{bit}\text{/}\mathrm{sym}$]$p_{r}=$0$$$p_{r}=$0.05$$$p_{r}=$0.05$$,
fit$p_{r}=$1$$theory$-20$$-15$$-10$$-5$$0$$5$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$\mathrm{SNR}_{b}$
[$\mathrm{dB}\text{/}$]$\hat{\xi}_{b}^{\prime}$ [$\mathrm{SNU}\text{/}$] (a)
Experiments with $M=2$

$-30$$-25$$-20$$-15$$-10$$-5$$0$$5$$-30$$-25$$-20$$-15$$-10$$-5$$0$$5$$\mathrm{SNR}_{b}$
[$\mathrm{dB}\text{/}$] (b) Simulations with $M=2$

Figure 4: Resulting system performance for $M=2$ in terms of the mutual
information between Alice and Bob $I_{\mathrm{AB}}$ and the excess noise in
Bob’s received symbols $\xi_{b}^{\prime}$. Experimental results are shown on
the left hand side (a). Simulation results are shown on the right hand side
(b). The small and solid markers indicate negative values.

##  III Discussion

The experimental and simulation results for $M=4$ and $M=2$ are shown in
figures 3 and 4 respectively. A typical value of $p_{r}$ is $0.5$
[jouguet2013] but in principle it could also be lower. Bob’s phase estimation
should only use the transmitted symbols that are revealed anyway. Therefore,
we investigate our system for a low $p_{r}$ of $0.05$. Additionally, we
investigate the cases $p_{r}=0$ and $p_{r}=1$. While $p_{r}=1$ represents an
upper bound for the performance of phase estimation, $p_{r}=0$ shows the
achievable performance without symbol revelation. Also, the case $p_{r}=0$
could be interesting for other fields than CV-QKD where symbol revelation is
not possible.

With $p_{r}=0$, the receiver sensitvity in the experimental system is
$\mathrm{SNR}_{b,\mathrm{min}}=$-6.9\text{\,}\mathrm{dB}\text{/}$$. This is in
very good agreement with the simulation, where the signal demodulation
completely fails at an $\mathrm{SNR}_{b}$ of $-8\text{\,}\mathrm{dB}\text{/}$.
Thus, despite the narrow linewidth fiber lasers, the receiver sensitivity is
limited by phase noise. In terms of received optical power, at least
$1.5\text{\,}\mathrm{ph}\text{/}\mathrm{sym}$ would be required. This is too
high for the discrete phase modulation scheme with $M=4$, where the
transmitted optical power should usually be lower than
$0.5\text{\,}\mathrm{ph}\text{/}\mathrm{sym}$ [leverrier2011]. However, even a
small revelation probability of $p_{r}=$0.05$$ is sufficient to enable
successful signal demodulation down to at least an $\mathrm{SNR}_{b}$ of
$-19.1\text{\,}\mathrm{dB}\text{/}$. In the experiment, the SNR was not
decreased further because the receiver calibration was not accurate enough to
allow for a reliable evaluation of $\mathrm{SNR}_{b}$ and the excess noise in
that regime. The experimentally confirmed receiver sensitivity of
$\mathrm{SNR}_{b,\mathrm{min}}<$-19.1\text{\,}\mathrm{dB}\text{/}$$ is
sufficient for CV-QKD over $26\text{\,}\mathrm{km}\text{/}$ and more. In the
simulation results, it can be seen that the carrier phase estimation can be
successful down to an $\mathrm{SNR}_{b}$ of $-35\text{\,}\mathrm{dB}\text{/}$.
However, there is already an increased probability of failure above
$-35\text{\,}\mathrm{dB}\text{/}$ as can be seen in the simulations for $M=2$
shown in figure 4, where the minimum $\mathrm{SNR}_{b}$ was
$-33\text{\,}\mathrm{dB}\text{/}$. Therefore, we take this higher value as the
receiver sensitivity.

The simulations show the phase noise limited performance. At an
$\mathrm{SNR}_{b}$ of $10\text{\,}\mathrm{dB}\text{/}$, there is no difference
in the resulting excess noise between different probabilities of revelation.
Apparently, in that case the uncertainty about the carrier phase is not
increased by the discrete phase modulation. This is the SNR regime of
classical communications where symbol error rates are low and the modulation
could also be canceled effectively before carrier phase estimation. In the
quantum regime of $\mathrm{SNR}_{b}<$0\text{\,}\mathrm{dB}\text{/}$$, the
modulation clearly affects the carrier phase estimation and hence the excess
noise. With $p_{r}=1$, $\xi_{b}^{\prime}$ exhibits a constant slope of about
$3\text{\,}\mathrm{dB}\text{/}$ per decade. This is the result of a decreasing
quantum signal power decreasing the excess noise while at the same time
decreasing the SNR increases the carrier phase uncertainty which induces
additional excess noise. With $p_{r}=0$, in the regime of
$\mathrm{SNR}_{b}<$5\text{\,}\mathrm{dB}\text{/}$$, there is a net increase of
$\xi_{b}^{\prime}$ with decreasing $\mathrm{SNR}_{b}$ until the demodulation
fails completely. Interestingly, for $p_{r}=0.05$, there is not much value of
the revelation around
$\mathrm{SNR}_{b}\approx$0\text{\,}\mathrm{dB}\text{/}$$. In this regime, the
unrevealed symbols still provide significant information about the carrier
phase. But in the regime of
$\mathrm{SNR}_{b}<$-8\text{\,}\mathrm{dB}\text{/}$$, the slope of
$\xi_{b}^{\prime}$ approaches the one of the case $p_{r}=1$. In this regime,
the revealed symbols provide more information about the carrier phase than the
unrevealed ones and the difference between $p_{r}=1$ and $p_{r}=0.05$ can be
interpreted as a difference of the effective symbol rate. Due to this, there
is a relatively constant excess noise penalty of about a factor of $9$ in the
$\mathrm{SNR}_{b}$ regime between $-10\text{\,}\mathrm{dB}\text{/}$ and
$-26\text{\,}\mathrm{dB}\text{/}$. Below an $\mathrm{SNR}_{b}$ of
$-26\text{\,}\mathrm{dB}\text{/}$, the slope is flattened.

Looking at the experimental results for the excess noise, we observe a much
higher level than in the simulations. Also, there is a proportional relation
between $\hat{\xi}_{b}^{\prime}$ and $\mathrm{SNR}_{b}$ that slightly flattens
below an $\mathrm{SNR}_{b}$ of $-7\text{\,}\mathrm{dB}\text{/}$. This confirms
that the excess noise mainly originates from signal distortions. The signal to
distortion ratio can be quantified for
$\mathrm{SNR}_{b}=$7.1\text{\,}\mathrm{dB}\text{/}$$ as
$P_{\mathrm{Q}}/\xi_{b}\approx$31.82$$, which is relatively low and could be
improved by optimizing the equalization concept.

For $M=2$, the results are shown in figure 4. The main difference to $M=4$ is
that the receiver sensitivity for the case $p_{r}=0$ is much lower. In the
simulations it was
$\mathrm{SNR}_{b,\mathrm{min}}=$-17\text{\,}\mathrm{dB}\text{/}$$ and in the
experiment $\mathrm{SNR}_{b,\mathrm{min}}=$-13.2\text{\,}\mathrm{dB}\text{/}$$
compared to the $-6.9\text{\,}\mathrm{dB}\text{/}$ for $M=4$. The difference
between experiment and simulations seems quite large. However, there is still
a significant probability of successful demodulation below these values in the
experiments and simulations. This indicates that the probability of failure
increases less steep with decreasing $\mathrm{SNR}_{b}$ as in the case of
$M=4$. With $p_{r}=$0.05$$, the transition to relying only on the revealed
symbols for phase estimation is much wider for $M=2$ compared to $M=4$. Even
at $\mathrm{SNR}_{b}=$-30\text{\,}\mathrm{dB}\text{/}$$, the excess noise is
slightly lower with $M=2$ compared to $M=4$. However, there is no significant
difference between these two cases in terms of experimental excess noise which
is dominated by other signal distortions.

$0$$10$$20$$30$$40$$50$$60$$70$$80$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$M=4$$M=2$$\bigstar$fiber
length [$\mathrm{km}\text{/}$]key rate [$\mathrm{bit}\text{/}\mathrm{sym}$]no
excess noisesimulation$\bigstar$experiment, $M=4$ Figure 5: Achievable key
rates for $p_{r}=0.05$, taking into account the polynomial fits for the excess
noise in figure 3 as well as the experimental parameters $\hat{\eta}$,
$\hat{\varepsilon}_{\mathrm{el}}$ and $\mathrm{SNR}_{b,\mathrm{min}}$. A
reconciliation efficiency of $95\text{\,}\mathrm{\char 37\relax}\text{/}$ and
a fiber loss of $0.2\text{\,}\mathrm{dB}\text{/}\mathrm{km}$ was assumed. The
launch power was optimized based on the security analysis in [leverrier2011].

In order to address the question of whether or not the performance of the
investigated system with a real LO and without pilots is sufficient for CV-
QKD, we calculated asymptotic secret key rates. The calculations are based on
the well established security proofs [leverrier2010] and [leverrier2011].
These assume a linear Gaussian channel and are therefore not fully general.
Recently, a security proof that does not incorporate these assumptions has
been presented [ghorai2019]. However, the question of how to get rid of the
linear Gaussian assumption is still under discussion.

We optimized the transmitted optical power with respect to a maximum key rate.
For this, we used the polynomial fits for the excess noise in the case of
$p_{r}=$0.05$$ that are plotted in the figures 3 and 4. An additional
constraint of the optimization was that $\mathrm{SNR}_{b}$ is not allowed to
be lower than in the experiments and simulations respectively and also not
allowed to be lower than the receiver sensitivity. The measured receiver
efficiency $\hat{\eta}$, the electronic noise ratio
$\hat{\varepsilon}_{\mathrm{el}}$ as well as the fiber length are also
included in the secret key rate calculations. The reconciliation efficiency is
assumed to be $95\text{\,}\mathrm{\char 37\relax}\text{/}$.

The resulting secret key rates for $p_{r}=$0.05$$ are shown in figure 5. As a
reference, the case of zero excess noise is also plotted. The simulation
corresponds to the case of the excess noise being dominated by carrier phase
uncertainty. In this case, distances of up to $72\text{\,}\mathrm{km}\text{/}$
could be achieved with $M=4$ and $32\text{\,}\mathrm{km}\text{/}$ with $M=2$.
This shows that laser phase noise is not an insurmountable obstacle for
designing a CV-QKD system with a real LO and without pilots. Even for the
experimental system, where the excess noise is much higher due to other signal
distortions, CV-QKD could be performed successfully with $M=4$ at a key rate
of $5.7\text{\times}{10}^{-4}\text{\,}\mathrm{bit}\text{/}\mathrm{sym}$ over
the experimental distance of $26\text{\,}\mathrm{km}\text{/}$. With $M=2$, the
experimental performance is not sufficient for CV-QKD. The worse performance
of $M=2$ is due to stricter requirements regarding the transmitted power and
excess noise [leverrier2011]. The experimental key rate can be improved
towards the simulation by reducing signal distortions in the experiment. This
can be achieved by improving the equalization concept of the system. To
achieve the long distances that are possible as shown by the simulations, the
receiver calibration should be improved. We expect that this can be achieved
by reducing the time difference between calibration and quantum signal
transmission in order to mitigate fluctuations of the receiver noise and LO
power. If the system is optimized such that it is limited by phase noise, the
parameter $p_{r}$ can be increased which directly reduces the excess noise for
even higher key rates and longer distances. Also, increasing $p_{r}$ can
enable the use of lasers with stronger phase noise while keeping the
achievable distance.

To conclude, we investigated a fiber-based quantum communication system for
CV-QKD that employs a real LO and works without any pilot tones. As Bob’s
clock and LO are free running, he relies only on the modulated quantum signal
itself to perform carrier phase estimation and timing recovery. Except for the
fact that the signal is attenuated before it is transmitted, the physical
implementation is identical to classical coherent systems. For the first time
it could be demonstrated experimentally, that such a system can be feasible
for CV-QKD. An important factor to achieve this is the particle smoother that
is used for carrier phase estimation, which is optimized using Monte Carlo
Markov chain methods. Based on the experimental results, the achievable key
rate over a distance of $26\text{\,}\mathrm{km}\text{/}$ is
$5.7\text{\times}{10}^{-4}\text{\,}\mathrm{bit}\text{/}\mathrm{sym}$
corresponding to $9.2\text{\,}\mathrm{Mbit}\text{/}\mathrm{s}$. The simulation
results indicate that this performance can even be largely improved. Possible
ways to achieve this is by reducing other signal distortions than phase noise
and by improving the receiver calibration. The achieved results are an
important milestone for CV-QKD. Also, the achieved receiver sensitivities
without symbol revelation ($p_{r}=0$) of $-6.9\text{\,}\mathrm{dB}\text{/}$
for $M=4$ and $-13.2\text{\,}\mathrm{dB}\text{/}$ for $M=2$ have, to our
knowledge, never been reached before. Therefore, the novel techniques are also
beneficial in other coherent systems where the SNR can be extremely low, such
as optical satellite communications.

##  IV Methods

###  IV.1 Experimental Details

The transmitter laser and LO are continuous wave DFB fiber lasers of type NKT
Koheras E15. Only one polarization of the dual polarization modulator is used,
the orthogonal one is biased to the zero transmission point. The baseband
quantum signal $m_{i}$ contains Alice’s symbols $a_{k}$ that are pulse-shaped
by the root raised cosine filter $h_{\mathrm{Q}}$ with a roll-off factor of
$0.1$ and a bandwidth of $17\text{\,}\mathrm{GHz}\text{/}$.

For each experimental scenario, defined by the transmitted optical power and
modulation order, Bob’s receiver is calibrated after quantum signal
transmission. For this, the quantum signal is deactivated to record the total
receiver noise. After that, the LO is also deactivated to record only the
electrical receiver noise. The noise signals undergo the same DSP routine as
the quantum signal but with deactivated timing recovery and phase estimation
because these methods do not alter the evaluated mean power of the received
noise sequence. The signals are recorded with a time difference of about
$3\text{\,}\mathrm{s}\text{/}$. For each quantum and noise signal, a total
number of 100 blocks are evaluated corresponding to
$6.55\text{\times}{10}^{6}$ symbols $b_{k}$. This does not include the first
signal block that is used to initialize $\hat{f}_{\mathrm{Q}}$. The
transmitted optical power in the initialization block is the same as for the
subsequent quantum communication.

###  IV.2 Simulation Details

In the simulations, the number of evaluated quantum symbols is the same as in
the experiment. The true phase $\varphi_{k}$ is directly accessible. Excess
noise due to inaccuracies in the estimated phase $\hat{\varphi_{k}}$ can be
calculated as [kleis2019]

$\displaystyle\xi_{b}^{\prime}=2\frac{P_{\mathrm{Q}}}{P_{\mathrm{SN}}}\langle\sin^{2}(\varphi_{k}-\hat{\varphi_{k}})\rangle,$
(5) $\displaystyle P_{\mathrm{Q}}=\langle|a_{k}|^{2}\rangle.$ (6)

In the simulations, only the total noise power $P_{\mathrm{TN}}$ is specified.
To calculate $\xi_{b}^{\prime}$, one must assume a specific electrical to shot
noise ratio $\varepsilon_{\mathrm{el}}$ to calculate

$\displaystyle\xi_{b}^{\prime}=2\,\mathrm{SNR}_{b}\,(1+\varepsilon_{\mathrm{el}})\,\langle\sin^{2}(\varphi_{k}-\hat{\varphi_{k}})\rangle.$
(7)

Here, $\varepsilon_{\mathrm{el}}$ is set to $0.7$ which corresponds to the
value that was measured in the experiments.

###  IV.3 Particle Smoother

The particle smoother is a Bayesian method that requires a measurement model
and a dynamic model formulated as probability densities. A detailed general
description of particle smoothing can be found in [sarkka2013]. The
transmitted symbols are $M$-ary phase modulated, meaning that

$\displaystyle
a_{k}\in\left\\{a^{(1)},\ldots,a^{(M)}\left|a^{(i)}=\mathrm{e}^{\,\mathrm{j}2\pi\frac{i}{M}}\right.\right\\}.$
(8)

In an AWGN channel, Bob’s measurement signal at the input of the particle
smoother can be written as

$\displaystyle\tilde{b}_{k}=a_{k}\mathrm{e}^{\,\mathrm{j}\varphi_{k}}+n_{k},$
(9)

where $n_{k}$ is a complex white Gaussian noise sequence with mean power
$P_{\mathrm{TN}}$ and $\varphi_{k}$ is the dynamic laser phase. The
normalization of Bob’s signal such that it matches $|a_{k}|=1$ is performed
using the noise power that is known from receiver calibration. Thus, for the
case $r_{k}=0$, the measurement probability density can be written as

$\displaystyle p\bigl{(}\tilde{b}_{k}|\varphi_{k}\bigr{)}$
$\displaystyle=\frac{1}{M}\sum\limits_{i=1}^{M}p\bigl{(}\tilde{b}_{k}|a^{(i)},\varphi_{k}\bigr{)}$
(10)
$\displaystyle=\frac{1}{M}\sum\limits_{i=1}^{M}\mathcal{N}\bigl{(}a^{(i)}\mathrm{e}^{\,\mathrm{j}\varphi_{k}},P_{\mathrm{TN}}/2\cdot
I_{2}\bigr{)}.$ (11)

Here, $\mathcal{N}(\mu,\Sigma)$ is a Gaussian probability density with mean
$\mu$ and covariance matrix $\Sigma$. $I_{2}$ denotes the 2-dimensional
identity matrix. If $r_{k}=1$, the transmitted symbol $a_{k}$ is known to Bob
and the measurement probability density reduces to
$p\bigl{(}\tilde{b}_{k}|\varphi_{k}\bigr{)}=\mathcal{N}\bigl{(}a_{k}\mathrm{e}^{\,\mathrm{j}\varphi_{k}},P_{\mathrm{TN}}/2\cdot
I_{2}\bigr{)}$. There are no unknown parameters in the measurement model.

The state space model $p(\mathbf{x}_{k}|\mathbf{x}_{k-1})$ describes the
dynamics of the variable $\mathbf{x}_{k}$. We define it as

$\displaystyle p(\mathbf{x}_{k}|\mathbf{x}_{k-1})$
$\displaystyle=p(\varphi_{k}|\varphi_{k-1},\Omega_{k-1})\cdot
p(\Omega_{k}|\Omega_{k-1})$ (12)
$\displaystyle=\mathcal{N}(\varphi_{k-1}+\Omega_{k-1},\sigma_{\varphi}^{2})\cdot\mathcal{N}(\Omega_{k-1},\sigma_{\Omega}^{2}).$
(13)

The variable $\Omega_{k}$ is a normalized frequency that models a drift of the
differential laser frequency as a random walk with variance
$\sigma^{2}_{\Omega}$ [piels2015]. Additionally, $\varphi_{k}$ is affected by
a random walk with variance $\sigma_{\varphi}^{2}$. Based on the stated
measurement and dynamic model, the particle smoother is implemented as a
bootstrap filter with $N=$200$$ particles in combination with a backward-
simulation particle smoother with $10$ trajectories. The resampling condition
of the bootstrap filter is $N_{\mathrm{eff}}<N/5$, where $N_{\mathrm{eff}}$ is
the effective number of particles.

###  IV.4 Bayesian Parameter Optimization

The particle smoother is capable of performing optimum phase estimation, given
that the state space model is an accurate representation of the measurement
and laser dynamics. Thus, optimizing the set of unknown parameters
$\theta=[\sigma^{2}_{\Omega},\sigma_{\varphi}^{2}]$ is essential. For this, we
used Monte Carlo Markov chain methods which are based on minimizing the energy
function $\Phi_{K}(\theta)$. The energy function has the property
$\mathrm{e}^{-\Phi_{K}(\theta)}\propto p(\theta|\tilde{b}_{1:K})$. Therefore,
minimizing it leads to the parameters $\hat{\theta}$ that provide the most
accurate description of the true phase dynamics. We performed the minimization
based on a received signal $\tilde{b}_{k}$ with
$\mathrm{SNR}_{b}=$11.5\text{\,}\mathrm{dB}\text{/}$$ and $p_{r}=1$ using an
extended Kalman filter to calculate $\Phi_{K}(\theta)$ for
$K=$8.1\text{\times}{10}^{6}$$. The extended Kalman filter is based on the
same state space model as the particle smoother. The minimum was found using
the simplex search method [lagarias1998]. The resulting optimized parameters
are
$\hat{\sigma}_{\Omega}^{2}=$1.66\text{\times}{10}^{-16}\text{\,}{\mathrm{rad}}^{2}\text{/}$$
and
$\hat{\sigma}_{\varphi}^{2}=$6.36\text{\times}{10}^{-9}\text{\,}{\mathrm{rad}}^{2}\text{/}$$.

###  IV.5 Timing Recovery

For signal processing, the total signal is split into consecutive blocks of
length $L=$6.8\text{\times}{10}^{5}$$ symbols with index $n$. The digital
filter and square algorithm for timing recovery is implemented as described in
[oerder1988] and calculates the complex Fourier coefficient $X_{n}$ using one
complete block. From the argument of $X_{n}$, the current timing offset can be
obtained. The SNR of $X_{n}$ scales approximately as
$\mathrm{SNR}_{X}\propto\mathrm{SNR}_{b}^{2}L$ in the regime of
$\mathrm{SNR}_{b}\ll 1$. Thus, the accuracy of timing recovery can be improved
by increasing $L$. However, the timing experiences a relatively constant drift
of $0.032$ symbol periods per block due to a clock offset between Alice and
Bob of $4.7\text{\times}{10}^{-8}\text{\,}\mathrm{Hz}\text{/}\mathrm{Hz}$
which means that the block length cannot be arbitrarily increased. At the
given block length, pulse shape and at an $\mathrm{SNR}_{b}$ of
$-20\text{\,}\mathrm{dB}\text{/}$, the resulting $\mathrm{SNR}_{X}$ is only
$-0.75\text{\,}\mathrm{dB}\text{/}$ which is too low to obtain sufficiently
low excess noise. Therefore, our approach is to track the argument $X_{n}$
using an extended Kalman filter. As this is very similar to the problem of
carrier phase tracking, the same state space model is used.

###  IV.6 Estimation of Mutual Information and its Theoretical Value

For mutual information estimation, Bob performs a hard decision on his
symbols, meaning that each $b_{k}$ is mapped to the symbol $a^{(i)}$ of the
discrete alphabet with the smallest euclidean distance. Bob’s symbols after
this mapping are denoted as $\hat{a}_{k}$. Based on the measured probabilities
$p(\hat{a}_{k}=a^{(i)},a_{k}=a^{(j)})$ we calculate the mutual information,
where all transmitted and received symbols are taken into account.

Assuming an ideal AWGN channel, the received symbols can be written as
$b_{k}=\sqrt{\mathrm{SNR}_{b}}\cdot a_{k}+n_{k}$, where
$a_{k}=\operatorname{e}^{\mathrm{j}\alpha_{k}}$ is Alice’s transmitted symbol
and $n_{k}$ is a white Gaussian noise sequence with a mean power of 1. The
probability density of the phase $\beta_{k}=\arg(b_{k})$ is then

$\displaystyle
p(\beta_{k})=\frac{1}{2\pi}\operatorname{e}^{-\mathrm{SNR}_{b}}+\frac{\sqrt{\mathrm{SNR}_{b}}}{2\sqrt{\pi}}\cos(\alpha_{k}-\beta_{k})\operatorname{e}^{-\mathrm{SNR}\sin^{2}(\alpha_{k}-\beta_{k})}\operatorname{erfc}\left(-\sqrt{\mathrm{SNR}_{b}}\cos(\alpha_{k}-\beta_{k})\right).$
(14)

By numerical integration of $p(\beta_{k})$ over each decision region and for
each transmitted symbol, we obtain the theoretical probabilities
$p(\hat{a}_{k}=a^{(i)},a_{k}=a^{(j)})$ that are required to calculate the
theoretical value of $I_{\mathrm{AB}}$.

###  IV.7 Estimation of Quantum Signal Power

For excess noise estimation, the quantum signal power $P_{\mathrm{Q}}$ must be
estimated. In our experiments, this is done by calculating

$\displaystyle
P_{\mathrm{Q}}=\left|\frac{1}{K}\sum\limits_{k=1}^{K}a_{k}b^{*}_{k}\right|^{2},$
(15)

where $K$ is the total number of evaluated symbols. For the estimation of
$P_{\mathrm{Q}}$, all transmitted and received symbols are used to obtain an
estimate that is accurate regardless of $p_{r}$.


# Molecular Structure Optimization based on Electrons-Nuclei Quantum Dynamics
Computation

Hirotoshi Hirai<EMAIL_ADDRESS>Takahiro Horiba  Soichi Shirai
Toyota Central R&D Labs., Inc., 41-1, Yokomichi, Nagakute, Aichi 480-1192,
Japan  Keita Kanno  Keita Arimitsu  Yuya O. Nakagawa  Sho Koh  QunaSys Inc.,
Aqua Hakusan Building 9F, 1-13-7 Hakusan, Bunkyo, Tokyo 113-0001, Japan

###### Abstract

A new concept of the molecular structure optimization method based on quantum
dynamics computations is presented. Nuclei are treated as quantum mechanical
particles, as are electrons, and the many-body wave function of the system is
optimized by the imaginary time evolution method. A demonstration with a
2-dimensional H${}_{2}^{+}$ molecule shows that the optimized nuclear
positions can be specified with a small number of observations. This method is
considered to be suitable for quantum computers, the development of which will
realize its application as a powerful method.

Introduction. Accurate quantum chemistry computations remain a challenge on
classical computers, especially for molecules with industrially relevant
sizes, despite significant efforts and developments of recent years. The
computational cost of exact methods for quantum chemistry on classical
computers grows exponentially with the molecular size [1, 2], whereas that can
be suppressed in polynomial scaling on quantum computers [3]. For this reason,
quantum chemistry computations have been considered to be a promising
application of quantum computers. By manipulating quantum states of matter and
taking advantage of their unique features, such as superposition and
entanglement, quantum computers promise to efficiently deliver accurate
results for many important problems in quantum chemistry, such as the
electronic structure of molecules [4, 5].

In quantum chemistry computations, the Born-Oppenheimer (BO) approximation [6]
has typically been used to speed up the computation of molecular wave
functions and other properties for large molecules. The BO approximation is
the assumption that the wave functions of atomic nuclei and electrons in
molecular systems can be treated separately, based on the fact that the nuclei
are much heavier than the electrons. In many cases of quantum chemistry
computations by classical computers, nuclei have been treated as point charges
(classical particles). Many quantum algorithms for quantum chemistry
computations on quantum computers are also typically based on the BO
approximation: for example, the quantum phase estimation (QPE) method [7, 8]
for a fault-tolerant quantum computer (FTQC) and the variational quantum
eigensolver (VQE) method [9, 10] for a noisy intermediate-scale quantum (NISQ)
[11] device are often employed to solve systems under the BO approximation. On
the other hand, it is also possible to think of methods beyond the BO
approximation with quantum computers, and such methods have been proposed in
literature [4, 12, 13, 14]. Kassal et al. [13], for instance, reported that a
completely non-adiabatic grid-based method on FTQCs, where nuclei are treated
as quantum particles and correlated with electrons, is not only more accurate
but also faster and more efficient than the methods based on the BO
approximation.

Geometry optimization of molecules is an important process to obtain the
equilibrium molecular structures in quantum chemistry computations [15]. Since
the physical and chemical properties of molecules are dependent on their
specific geometrical structures, elucidation of the optimized structures
enables the prediction of properties and identification of chemical products.
In the conventional molecular structure optimization methods, which are based
on the BO approximation and treat nuclei as classical particles, the
electronic states should be computed to obtain the forces acting on nuclei and
update the molecular structure for each optimization step. These conventional
methods can be accelerated through the use of quantum computers to perform the
electronic state calculations.

However, a large number of measurements (circuit executions) are required
because it is necessary to repeat the electronic state computations on the
quantum computer at each iteration in the geometry optimization process.
Furthermore, as gradient-driven descent routines are often used to optimize
the molecular structure, the system tends to be relaxed to the closest local
minimum from an initial structure for the geometry optimization. It is thus a
difficult task to find the global minimum structure in the systems which have
multiple local minima (i.e. many isomers), such as alloy cluster systems [15].

In this study, we propose a method for optimizing molecular structures based
on quantum dynamics computations with working on an FTQC in mind. In our
method, the many-body wave functions of nuclei and electrons are directly
treated as wavepackets and optimized by the imaginary time evolution method.
The many-body wave function on the lattice can be mapped to a state on a
quantum computer with the number of qubits growing only logarithmically with
respect to the number of lattice points [12, 5]. Although the imaginary time
evolution is a non-unitary operation, probabilistic methods for FTQC are
proposed [16, 17, 18], and the method of linear combination of unitaries could
be also be useful [19, 20]. A variational method for NISQ devices is also
proposed [21]. We leave a concrete implementation of our proposal on quantum
computers as a future work. Since nuclei are treated as quantum particles, the
nuclear quantum effects are naturally included in this method. If we start
with an initial many-body wave function that covers all molecular structures,
it is expected theoretically that the global minimum of the molecular
structure can be obtained. The optimized many-body wave functions have large
stochastic amplitudes at the most stable structure; therefore, the nuclei
positions can be specified, i.e. the optimization of the molecular structure
is possible, with a small number of observations. We demonstrate the concept
of the molecular structure optimization method based on quantum dynamics
computations with a classical computer and the size of the molecules that can
be handled.

Quantum imaginary time evolution. In non-relativistic quantum dynamics
computations, the time evolution of quantum systems can be described by the
time-dependent Schrödinger equation:

$i\hbar\frac{\partial}{\partial t}\psi(t)=H\psi(t),$ (1)

where $H$ is the Hamiltonian that describes the system. The solution of the
time-dependent Schrödinger equation is formally written by

$\psi(t)=e^{-i\frac{H}{\hbar}t}\psi(0).$ (2)

That is, if an appropriate initial wave function $\psi(0)$ can be prepared,
then the time evolution of the system can be obtained by application of the
time evolution operator,

$U(t)=e^{-i\frac{H}{\hbar}t}.$ (3)

If $\psi(0)$ is the eigenstate of $H$, i.e. $\psi(0)=\psi_{i}$ and
$H\psi_{i}=E_{i}\psi_{i}$, then Eq. (3) becomes

$U(t)\psi(0)=e^{-i\frac{E_{i}}{\hbar}t}\psi_{i}.$ (4)

If not (in general cases), then the exponential operator must be applied to
the wave function to compute the time evolution of the system. A number of
methods have been developed to apply the exponential operator to the wave
function [22, 23]. We here explain the method based on the second-order
Suzuki-Trotter decomposition [5, 4, 24].

Denoting the kinetic energy term of $H$ by $T$ and the potential energy term
by $V$, the second-order Suzuki-Trotter decomposition is

$U(t)=e^{-i\frac{T}{\hbar}\frac{t}{2}}e^{-i\frac{V}{\hbar}t}e^{-i\frac{T}{\hbar}\frac{t}{2}}+O(t^{3}).$
(5)

The error arises because $T$ and $V$ are non-commutative. To reduce the error,
we set $t=N_{\mathrm{step}}dt$ and express the time evolution operator as the
product of $N_{\mathrm{step}}$ operators in time increments of $dt$,

$U(t)=(U(dt))^{N_{\mathrm{step}}}=(e^{-i\frac{H}{\hbar}dt})^{N_{\mathrm{step}}}.$
(6)

Because $T$ is diagonal in wavenumber space, when we apply
$e^{-i\frac{T}{\hbar}\frac{dt}{2}}$ in the second-order Suzuki-Trotter
decomposition, the wave function is expressed in wavenumber space by
performing the fast Fourier transformation (FFT) for the wave funciton in real
space, when running on classical computers. Similarly, $V$ is diagonal in real
space, so the operator $e^{-i\frac{V}{\hbar}dt}$ is applied to the wave
function in real space after application of the inverse FFT to the wave
function in wavenumber space, $e^{-i\frac{T}{\hbar}\frac{dt}{2}}\psi$. Using
these procedures, the application of the time evolution operator $U(dt)$ can
be computed without expanding the exponential function operator. By repeating
these calculations $N_{\mathrm{step}}$ times, the wave function at the desired
time $t$ can be obtained.

Now let us convert time into an imaginary number, $it\rightarrow\tau$. The
solution of the time-dependent Schrödinger equation then becomes

$\psi(\tau)=e^{-\frac{H}{\hbar}\tau}\psi(0).$ (7)

When we formally expand $\psi(0)$ by the eigenfunctions of $H$,
$\\{\phi_{i}\\}_{i}$, we have

$\psi(\tau)=\sum_{i}c_{i}e^{-\frac{E_{i}}{\hbar}\tau}\phi_{i},$ (8)

where $c_{i}\in\mathbb{C}$ is a coefficient for the expandsion
$\psi(0)=\sum_{i}c_{i}\phi_{i}$, and $E_{i}$ is the eigenvalue of $\phi_{i}$.
If we increase $\tau$, the factor $e^{-\frac{E_{i}}{\hbar}\tau}$ attenuates
faster for large $E_{i}$, and the ground state $\phi_{0}$, which has the
smallest eigenvalue $E_{0}$, remains until the end as long as the coefficient
(overlap between the initial state $\psi(0)$ and the ground state) $c_{0}$ is
non-zero. This is how we obtain the ground state by imaginary time evolution.
Note that the norm of the wave function changes during the imaginary time
evolution, so it is necessary to renormalize the wave function after the
evolution.

Geometry optimization with fully quantum wave function for nuclei and
electrons. In this study, in contrast to the conventional method based on the
BO approximation, we propose a geometry optimization method where the atomic
nuclei are regarded as quantum mechanical particles as well as the electrons
and treated by wave functions. We consider the many-body wave function
$\psi(R,r,\tau=0)$, including the degrees of freedom of nuclei and electrons,
where $R$ represents the coordinates of $N_{\mathrm{nuc}}$ nuclei,
$R={R_{1},R_{2},...,R_{N_{\mathrm{nuc}}}}$, and $r$ represents the coordinates
of $n$ electrons, $r={r_{1},r_{2},...,r_{n}}$. Starting from any initial wave
function (ultimately, it may be a constant over the entire simulation region
in real space, but should have non-zero overlap with the desired ground
state), the most stable state (ground state) of the nuclei and the electrons
can be obtained by the imaginary time evolution of the initial wave function
under the Hamiltonian

$\begin{split}H&=-\sum_{i}^{n}\frac{\hbar^{2}}{2m_{i}}\nabla_{r_{i}}^{2}-\sum_{i}^{N_{\mathrm{nuc}}}\frac{\hbar^{2}}{2M_{i}}\nabla_{R_{i}}^{2}+\sum_{i<j}\frac{q_{e}^{2}}{|r_{i}-r_{j}|}\\\
&+\sum_{i<j}\frac{Q_{i}Q_{j}}{|R_{i}-R_{j}|}-\sum_{i}^{N_{\mathrm{nuc}}}\sum_{j}^{n}\frac{q_{e}Q_{i}}{|R_{i}-r_{j}|},\end{split}$
(9)

where $Q_{i}$ is a electric charge of $i$th nuclear and $q_{e}(>0)$ is a
charge of electron. The optimized wave function for the nuclei, which can be
obtained by integrating out the electronic degrees of freedom, does not spread
like electrons and has sharp peaks (cusps) at the most stable positions;
therefore, the position of each nuclei should be identifiable with several
repetitions of the projective measurement for the positions of the nuclei. The
stable structure can be obtained with a high probability by performing the
measurement to the state after imaginary time evolution. This calculation
includes quantum effects such as the zero-point oscillation effect and the
non-adiabatic effects that are ignored when the atomic nuclei are treated as
point charges [25]. If we take an initial state (for the nuclei part of the
wave function) as the superposition state for all possible molecular
structures, the most stable structure among all isomers of the molecule would
be obtained. Because the search for the most stable isomer generally becomes
difficult as the number of atoms increases [15], it is an appealing feature of
this method that treats the atomic nuclei as quantum mechanical particles.

Such an approach could not be taken to date due to the difficulty of
describing many-body wave functions with a classical computer. In contrast, a
many-body wave function can be described with a decent number of qubits on a
quantum computer. This is discussed in more detail later.

Numerical demonstration of proof-of-concept for our method. We present
numerical demonstrations of our method by taking a two-dimensional (2D)
H${}_{2}^{+}$ molecule as an example. This is a three-particle system with two
protons and an electron that is naively described as a six-dimensional system,
which will be difficult to calculate with a classical computer. For example,
if we take a grid of $2^{5}=32$ per dimension in that space, the wave function
has $2^{30}=1,073,741,824$ elements, and the Hamiltonian has
$2^{60}=1,152,921,504,606,846,976$ elements in a naive estimate. On the other
hand, in a quantum computer, a wave function can be expressed with 30 qubits
because $N$ qubits can represent $2^{N}$ variables, i.e. $2^{N}$ wave function
elements on grids. In addition, the number of qubits scales linearly with the
number of particles.

In numerical demonstrations, we perform a dimensional reduction to the system
so that the calculation can be performed on a classical computer. First, the
center of mass of the system is fixed at the origin $(0,0)$. As a result, the
degrees of freedom of translational motion of the molecule can be ignored. Now
that the purpose is to optimize the molecular structure, we only need to focus
on the internal degrees of freedom. Classically, the rotational motion around
the center of mass has also been frozen, but quantum mechanically, the
rotational motion and the vibration mode are coupled, so that the rotational
motion cannot be separated. Even if the rotational motion is left, observing
the optimized wave function gives only one molecular structure in a specific
orientation of the rotational degrees of freedom. When we obtain multiple
results for the measurements of the positions of the nuclei, the obtained
structures are described in terms of internal coordinates which allows
rotational-free description of molecules. The mean of each nucleus position is
then calculated, which gives the molecular structure. Next, the relative
coordinates $\vec{R}$ between the two nuclei are introduced to further reduce
the dimension. The relationship between the reduced coordinates and the
original coordinates is given as

$\vec{R}=\vec{r}_{2}-\vec{r}_{1},\>\vec{r}_{c}=\vec{r}_{e}-\frac{m_{1}\vec{r_{1}}+m_{2}\vec{r_{2}}}{M_{n}},$
(10)
$\vec{R}_{\mathrm{cm}}=\frac{m_{1}\vec{r_{1}}+m_{2}\vec{r_{2}}+m_{e}\vec{r_{e}}}{M}\>(\rightarrow\text{set
to }(0,0)),$ (11)

where $\vec{r}_{1}$ and $\vec{r}_{2}$ ($m_{1}$ and $m_{2}$) are the
coordinates (masses) of the two protons, $\vec{r}_{e}$ ($m_{e}$) are the
coordinates (mass) of the electron, $\vec{r}_{c}$ are the coordinates of the
electron as seen from the center of mass of the nuclei, and
$\vec{R}_{\mathrm{cm}}$ is the center of mass of the system and fixed at the
origin. We here also define the total mass of the system $M=m_{1}+m_{2}+m_{e}$
and that of the protons $M_{n}=m_{1}+m_{2}$. The reduced coordinates have four
degrees of freedom, $\vec{R}=(R_{x},R_{y})$ and $\vec{r}_{c}=(r_{cx},r_{cy})$.
The reduced masses for $\vec{R}$ and $\vec{r}_{c}$ are
$\mu_{R}^{-1}=m_{1}^{-1}+m_{2}^{-1}$ and
$\mu_{r_{c}}^{-1}=m_{e}^{-1}+M_{n}^{-1}$, respectively. By using the reduced
coordinates, the Hamiltonian of the system becomes

$\begin{split}H&=-\frac{1}{2\mu_{R}}\nabla_{R}^{2}-\frac{1}{2\mu_{r_{c}}}\nabla_{r_{c}}^{2}\\\
&+\frac{1}{|\vec{R}|}-\frac{1}{|\vec{r}_{c}+\frac{m2}{M_{n}}\vec{R}|}-\frac{1}{|\vec{r}_{c}-\frac{m1}{M_{n}}\vec{R}|}.\end{split}$
(12)

Here the atomic unit is introduced for simplicity ($\hbar=1$, $Q_{p}=1$,
$q_{e}=1$). We take an evenly spaced isotropic grid of $2^{5}$ elements per
dimension in each two-dimensional spaces for $\vec{R},\vec{r}_{c}$. See
Supplemental Material for more details on computational setup.

We present results of numerical demonstrations result for the method using the
2D H${}_{2}^{+}$ molecule. The conditional probability density distribution of
the relative distance $\vec{R}$ between the protons when the electron is near
the origin (0.0375, 0.0375) (note that the grid used in this study does not
contain (0, 0)), $|\psi(R_{x},R_{y},r_{cx}=0.0375,r_{cy}=0.0375)|^{2}$ is
shown in FIG. 1.

Figure 1: Conditional probability density distribution of the relative
coordinates of protons $\vec{R}$ for the H${}_{2}^{+}$ molecule when an
electron is near the origin point (0.0375, 0.0375).

We find that the probability density distribution of the relative coordinates
$\vec{R}$ between the protons has a donut-shaped structure. This means that
the H${}_{2}^{+}$ molecule can be oriented in any direction because it has
rotational degrees of freedom. This probability density distribution must be
donut-shaped and uniform for all directions; however, the probability is large
and small in specific directions in FIG. 1. The reason for this is considered
to be that the roughness of the grid space used in this study cannot represent
an appropriate distance for the equilibrium internuclear distance in specific
directions and the probability density distributions for such directions would
be small, and the probability density distributions for the specific
directions where the grid is well-matched with the equilibrium internuclear
distance would be large. We also calculate the expectation value (mean) of the
proton-proton distance $|\vec{R}|$, defined as

$\int|\vec{R}||\psi(\vec{R},\vec{r}_{c})|^{2}d\vec{R}d\vec{r}_{c},$ (13)

for our simulation data in FIG. 1. The result is 0.41 bohr, which is slightly
larger than the value of 0.37 bohr where the protons were treated as point
charges (see Supplemental Material). This larger expectation value is
considered to be due to the spread of the nuclear wave function and its
binding potential, which is very tight for smaller $|\vec{R}|$ but loose for
larger $|\vec{R}|$.

Next, we estimate how many measurements for the wave function would be
required to optimize the molecular structure within an acceptable error.
$|\psi(\vec{R},\vec{r_{c}})|^{2}$ gives the simultaneous probability of
finding relative proton-to-proton coordinates in $\vec{R}$ and electrons in
$\vec{r_{c}}$. To optimize the molecular structure, it is enough to measure
the positions of the nuclei coordinates ($\vec{R}$ in this case). We simulate
the measurements (observations) of $\vec{R}$ for the optimized wave function
by randomly sampling $\vec{R}$ with the probability
$p(\vec{R})=\int|\psi(\vec{R},\vec{r_{c}})|^{2}d\vec{r}_{c}$ using pseudo-
random numbers generated by the Mersenne Twister method [26].

Figure 2:  Mean of the sampled bond distance $|\vec{R}|$ drawn from the
optimized wave function. We plot the number of observations $N_{\mathrm{obs}}$
versus the mean up to $N_{\mathrm{obs}}$ observations,
$M(N_{\mathrm{obs}})=\frac{1}{N_{\mathrm{obs}}}\sum_{i=1}^{N_{\mathrm{obs}}}|\vec{R}_{i}|$,
where $\vec{R}_{i}$ is an $i$-th sample of the observed $\vec{R}$. The red
(blue) dotted line represents the exact mean (standard deviation) of
$|\vec{R}|$.

Figure 2 shows how the mean of $|\vec{R}|$ for $N_{\mathrm{obs}}$ observations
converges as we increase $N_{\mathrm{obs}}$. The exact mean and the standard
deviation (corresponding for the zero-point vibration of the ground state)
calculated from $p(\vec{R})$ are also shown in dotted lines. We see that the
mean of the sampled $|\vec{R}|$ gets close to the exact mean with sufficient
accuracy after as small as 200 observations. The mean $|\vec{R}|$ obtained
from 200 observations was 0.4082 bohr while exact one is 0.41 bohr as
explained earlier. The standard deviation for 200 observations was 0.0611 bohr
and the mode was 0.4142 bohr. The reason why the converged bond length
$\vec{R}$ was obtained with such a small number of observations is that the
wave function has a peak (cusp) at the most stable nuclei positions and that
the measurements will pick it up with high probability. It should be noted
that as the size of the molecule increases, the number of vibrational modes
also increases. As a result, the fluctuation width of the zero-point vibration
also increases, which may affect the required number of measurements to
determine the molecular structure with a certain precision. As discussed in
the Supplemental Material, the standard deviation of the zero-point
fluctuation scales as $O(\sqrt{N})$, where $N$ is the number of nuclei.
Thereby the required number of measurements for a certain precision will scale
as $O(N)$.

Finally, we illustrate how the presented result above changes when we consider
heavier nuclei. It is expected that as a nucleus becomes heavier, the less
wave-like properties it has and the closer it is to a classical particle. It
is predicted that the peaks at the nuclear positions of the wave function will
be even sharper and converged structures of molecules will be obtained with a
smaller number of observations. To see this, we consider the deuterium
molecular ion D${}_{2}^{+}$, in which the protons in the hydrogen molecular
ion are replaced with deuterons (the mass is twice that of proton), and the
tritium molecular ion T${}_{2}^{+}$, in which tritons (the mass is three times
that of proton) are substituted. We perform the same calculations as for
H${}_{2}^{+}$ for these isotopes. Figure 3 show the conditional probability
density distributions of the relative coordinates $\vec{R}$ for deuterons and
tritons, respectively, when the electron is near the origin (0.0375, 0.0375).

(a) D${}_{2}^{+}$ molecule

(b) T${}_{2}^{+}$ molecule

Figure 3: Conditional probability density distributions of the relative
coordinates of nuclei for (a) D${}_{2}^{+}$ and (b) T${}_{2}^{+}$ molecules
when the electron is near the origin point (0.0375, 0.0375).

The absolute values of the peaks are higher than those for H${}_{2}^{+}$, and
the spreads of the wave functions are narrowed. This tendency was more
significant for the heavier T${}_{2}^{+}$. In addition, the peak bias
increases for specific directions with the heavier isotopes D${}_{2}^{+}$ and
T${}_{2}^{+}$ because as the wave functions become narrower, it becomes
difficult for the grid points to exist at the positions where the peaks in the
continuous limit (infinite number of meshes of the gird) because the grid size
is not fine enough. Subsequently, as in the case of H${}_{2}^{+}$, we simulate
how many observations were required to optimize the molecular structure of
D${}_{2}^{+}$ and T${}_{2}^{+}$ within an acceptable error. Figure 4 plots the
same quantity as Fig. 2, the mean of the sampled $|\vec{R}|$ from the
optimized wave function, for D${}_{2}^{+}$ and T${}_{2}^{+}$ molecules. The
exact mean and the standard deviation are also shown in dotted lines.

(a) D${}_{2}^{+}$ molecule

(b) T${}_{2}^{+}$ molecule

Figure 4:  Mean of the sampled bond distance $|\vec{R}|$ drawn from the
optimized wave function for (a) D${}_{2}^{+}$ and (b) T${}_{2}^{+}$ molecules.
We plot the number of observations $N_{\mathrm{obs}}$ versus the mean up to
$N_{\mathrm{obs}}$ observations,
$M(N_{\mathrm{obs}})=\frac{1}{N_{\mathrm{obs}}}\sum_{i=1}^{N_{\mathrm{obs}}}|\vec{R}_{i}|$,
where $\vec{R}_{i}$ is an $i$-th sample of the observed $\vec{R}$. The red
(blue) dotted line represents the exact mean (standard deviation) of
$|\vec{R}|$.

Besides, the means, standard deviations, and modes obtained from 200
observations are summarized in TABLE. 1 along with the result for
H${}_{2}^{+}$ for comparison.

Table 1: Mean, standard deviation and mode of the 200 times observation for $|\vec{R}|$ in bohr. System | Mean | Standard Deviation | Mode
---|---|---|---
H${}_{2}^{+}$ | 0.4082 | 0.0611 | 0.4142
D${}_{2}^{+}$ | 0.3968 | 0.0494 | 0.3712
T${}_{2}^{+}$ | 0.3778 | 0.0355 | 0.3712

These results in FIG. 4 and TABLE. 1 inidicate that the variation of the
observed $\vec{R}$ becomes smaller and the convergence becomes faster as the
mass increases. This is because of the smaller fluctuation widths of the zero-
point vibrations for heavier nuclei. Since all the nuclei after the He atom in
the periodic table are heavier than triton, it is expected that molecular
structures obtained by measurements will converge even faster in term of the
number of the measurements. We briefly comment on the bond length $|\vec{R}|$
before ending the explanation of numerical demonstrations. The mean of the
bond length approaches the value for point charges (0.37 bohr) as the mass
increases. The tendency for the bond length to become shorter as the mass
increases has also been confirmed for the 3D hydrogen molecule and its isotope
molecules [27]. The modes of the sampled $|\vec{R}|$ are the same for
D${}_{2}^{+}$ and T${}_{2}^{+}$ because the values between 0.3712 and 0.4142
cannot be taken due to rough grid used in our calculations. This result
indicates that the step size of the grid must be sufficiently smaller than the
desired resolution of the molecular structures. For calculations on a
classical computer, making the grid finer leads to a large increase in
computational cost, but on a quantum computer it will scales in polynomial
time [13].

Summary and outlook. In summary, the concept of a molecular structure
optimization method using quantum dynamics computation has been presented. The
nuclei have been treated as quantum mechanical particles as were the
electrons, and the many-body wave function of the system was optimized using
the imaginary time evolution method. The optimized wave function has a large
probability amplitude at the most stable structure of the nuclei, which allows
us to determine the optimized nuclear positions with a small number of
observations. It is expected that our method has a favorable feature that the
most stable isomer structure can be obtained for complex systems without being
trapped in local minima by virtue of the imaginary time evolution. This method
may be particularly promising for searching for the most stable structure for
systems with many isomers, such as metal alloy clusters. Ultimately, it may
also be useful to solve protein folding problems. Another aspect of our method
is that it includes nuclear quantum effects (zero-point oscillation, non-
adiabatic correction, etc.) that are typically ignored by conventional
methods. The coulomb interactions between multiple particles are explicitly
incorporated and there is no approximation in nucleus-nucleus, electron-
nucleus and electron-electron interactions. Although executing our method for
industrially-interesting large molecules with classical computers is difficult
because of huge computational costs in treating the fully quantum wave
function for nuclei and electrons, quantum computers (possibly FTQC) can be an
appealing candidate to run our method. Our proposal can give a new insight for
quantum chemistry computations on quantum computers.

## References

  * [1] P.J. Knowles and N.C. Handy.  A new determinant-based full configuration interaction method.  Chemical Physics Letters, 111(4):315–321, 1984.
  * [2] Jeppe Olsen, Björn O. Roos, Poul Jo/rgensen, and Hans Jo/rgen Aa. Jensen.  Determinant based configuration interaction algorithms for complete and restricted configuration interaction spaces.  The Journal of Chemical Physics, 89(4):2185–2192, 1988.
  * [3] Michael A. Nielsen and Isaac L. Chuang.  Quantum Computation and Quantum Information: 10th Anniversary Edition.  Cambridge University Press, New York, NY, USA, 10th edition, 2011.
  * [4] Christof Zalka.  Simulating quantum systems on a quantum computer.  Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 454(1969):313–322, 1998.
  * [5] Christof Zalka.  Efficient simulation of quantum systems by quantum computers.  Fortschritte der Physik: Progress of Physics, 46(6-8):877–879, 1998\.
  * [6] Jean-Michel Combes, Pierre Duclos, and Ruedi Seiler.  The born-oppenheimer approximation.  In Rigorous atomic and molecular physics, pages 185–213. Springer, 1981.
  * [7] A Yu Kitaev.  Quantum measurements and the abelian stabilizer problem.  arXiv preprint quant-ph/9511026, 1995.
  * [8] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca.  Quantum algorithms revisited.  Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 454(1969):339–354, 1998.
  * [9] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O’Brien.  A variational eigenvalue solver on a photonic quantum processor.  Nature Communications, 5(1):4213, 2014.
  * [10] Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M Chow, and Jay M Gambetta.  Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.  Nature, 549(7671):242–246, 2017.
  * [11] John Preskill.  Quantum computing in the nisq era and beyond.  Quantum, 2:79, 2018.
  * [12] Stieve Wiesner.  Simulations of many-body quantum systems by a quantum computer.  arXiv:quant-ph/9603028, 1998.
  * [13] Ivan Kassal, Stephen P Jordan, Peter J Love, Masoud Mohseni, and Alán Aspuru-Guzik.  Polynomial-time quantum algorithm for the simulation of chemical dynamics.  Proceedings of the National Academy of Sciences, 105(48):18681–18686, 2008.
  * [14] Pauline J. Ollitrault, Guglielmo Mazzola, and Ivano Tavernelli.  Non-adiabatic molecular quantum dynamics with quantum computers.  Physical Review Letters, 125(26):260511, December 2020.  arXiv: 2006.09405.
  * [15] Frank Eckert, Peter Pulay, and Hans-Joachim Werner.  Ab initio geometry optimization for large molecules.  Journal of computational chemistry, 18(12):1473–1483, 1997.
  * [16] Colin P Williams.  Probabilistic nonunitary quantum computing.  In Quantum Information and Computation II, volume 5436, pages 297–306. International Society for Optics and Photonics, 2004.
  * [17] Hiroaki Terashima and Masahito Ueda.  Nonunitary quantum circuit.  International Journal of Quantum Information, 3(04):633–647, 2005\.
  * [18] Tong Liu, Jin-Guo Liu, and Heng Fan.  Probabilistic nonunitary gate in imaginary time evolution.  Quantum Information Processing, 20(6):1–21, 2021.
  * [19] Andrew M Childs and Nathan Wiebe.  Hamiltonian simulation using linear combinations of unitary operations.  arXiv preprint arXiv:1202.5822, 2012.
  * [20] Dominic W Berry, Andrew M Childs, Richard Cleve, Robin Kothari, and Rolando D Somma.  Simulating hamiltonian dynamics with a truncated taylor series.  Physical review letters, 114(9):090502, 2015.
  * [21] Sam McArdle, Tyson Jones, Suguru Endo, Ying Li, Simon C Benjamin, and Xiao Yuan.  Variational ansatz-based quantum simulation of imaginary time evolution.  npj Quantum Information, 5(1):1–6, 2019.
  * [22] Claude Leforestier, RH Bisseling, Charly Cerjan, MD Feit, Rich Friesner, A Guldberg, A Hammerich, G Jolicard, W Karrlein, H-D Meyer, et al.  A comparison of different propagation schemes for the time dependent schrödinger equation.  Journal of Computational Physics, 94(1):59–80, 1991.
  * [23] W Van Dijk and FM Toyama.  Accurate numerical solutions of the time-dependent schrödinger equation.  Physical Review E, 75(3):036707, 2007.
  * [24] Thomas Barthel and Yikang Zhang.  Optimized lie–trotter–suzuki decompositions for two and three non-commuting terms.  Annals of Physics, 418:168165, 2020.
  * [25] Hiromi Nakai.  Nuclear orbital plus molecular orbital theory: Simultaneous determination of nuclear and electronic wave functions without born–oppenheimer approximation.  International Journal of Quantum Chemistry, 107(14):2849–2869, 2007\.
  * [26] Makoto Matsumoto and Takuji Nishimura.  Mersenne twister: A 623-dimensionally equidistributed uniform pseudo-random number generator.  ACM Trans. Model. Comput. Simul., 8(1):3–30, January 1998.
  * [27] Steven L Mielke, Bruce C Garrett, and Kirk A Peterson.  A hierarchical family of global analytic born–oppenheimer potential energy surfaces for the h+ h 2 reaction ranging in quality from double-zeta to the complete basis set limit.  The Journal of chemical physics, 116(10):4142–4161, 2002.