# Automated machine learning for secure key rate in discrete-modulated
continuous-variable quantum key distribution

Zhi-Ping Liu  National Laboratory of Solid State Microstructures, School of
Physics and Collaborative Innovation Center of Advanced Microstructures,
Nanjing University, Nanjing 210093, China  Min-Gang Zhou  National Laboratory
of Solid State Microstructures, School of Physics and Collaborative Innovation
Center of Advanced Microstructures, Nanjing University, Nanjing 210093, China
Wen-Bo Liu  National Laboratory of Solid State Microstructures, School of
Physics and Collaborative Innovation Center of Advanced Microstructures,
Nanjing University, Nanjing 210093, China  Wen-Bo Liu  National Laboratory of
Solid State Microstructures, School of Physics and Collaborative Innovation
Center of Advanced Microstructures, Nanjing University, Nanjing 210093, China
Jie Gu  National Laboratory of Solid State Microstructures, School of Physics
and Collaborative Innovation Center of Advanced Microstructures, Nanjing
University, Nanjing 210093, China  Hua-Lei Yin<EMAIL_ADDRESS>National
Laboratory of Solid State Microstructures, School of Physics and Collaborative
Innovation Center of Advanced Microstructures, Nanjing University, Nanjing
210093, China  Zeng-Bing Chen<EMAIL_ADDRESS>National Laboratory of Solid
State Microstructures, School of Physics and Collaborative Innovation Center
of Advanced Microstructures, Nanjing University, Nanjing 210093, China

###### Abstract

Continuous-variable quantum key distribution (CV QKD) with discrete modulation
has attracted increasing attention due to its experimental simplicity, lower-
cost implementation and compatibility with classical optical communication.
Correspondingly, some novel numerical methods have been proposed to analyze
the security of these protocols against collective attacks, which promotes key
rates over one hundred kilometers of fiber distance. However, numerical
methods are limited by their calculation time and resource consumption, for
which they cannot play more roles on mobile platforms in quantum networks. To
improve this issue, a neural network model predicting key rates in nearly real
time has been proposed previously. Here, we go further and show a neural
network model combined with Bayesian optimization. This model automatically
designs the best architecture of neural network computing key rates in real
time. We demonstrate our model with two variants of CV QKD protocols with
quaternary modulation. The results show high reliability with secure
probability as high as $99.15\%-99.59\%$, considerable tightness and high
efficiency with speedup of approximately $10^{7}$ in both cases. This
inspiring model enables the real-time computation of unstructured quantum key
distribution protocols’ key rate more automatically and efficiently, which has
met the growing needs of implementing QKD protocols on moving platforms.

##  I Introduction

In recent decades, machine learning (ML) has gained impressive breakthroughs
that deeply impact both industry and academia, including autonomous driving
Grigorescu et al. (2020); Levinson et al. (2011), natural language processing
Deng et al. (2013); Young et al. (2018), protein structure prediction Jumper
et al. (2021) and even proving mathematical conjectures Davies et al. (2021).
ML aims to recognize patterns in data, especially multidimensional data, and
generalize them to new instances, which contributes to automating tasks and
reveals hidden patterns beyond humans intuition. This modern information-
processing technology also benefits solving intractable quantum tasks, since
quantum tasks are usually counterintuitive and involve high dimensions.
Several significant advances have been made by applying ML to quantum physics,
from classifying quantum states Gao et al. (2018); Yang et al. (2019); Ahmed
et al. (2021), quantum control Bukov et al. (2018); Lumino et al. (2018); Niu
et al. (2019) to quantum metrology Hentschel and Sanders (2011).

Quantum key distribution (QKD) enables unconditional security between two
legitimate users (Alice and Bob) against any eavesdropper called Eve Bennett
and Brassard (1984); Ekert (1991), which is guaranteed by quantum mechanics
laws Shor and Preskill (2000). According to different detection methods, QKD
is currently divided into two categories: discrete-variable (DV) QKD Lo et al.
(2014); Gisin et al. (2002) and continuous-variable (CV) QKD Grosshans and
Grangier (2002); Lance et al. (2005); Huang et al. (2015); Yin et al. (2019).
Between these two categories, CV QKD has unique edges on a higher secret key
rate and excellent compatibility with standard communication components
Fossier et al. (2009); Huang et al. (2016); Jouguet et al. (2012), which
enables CV QKD to be competitive at a metropolitan distance Pirandola et al.
(2020). To enhance the practicality of CV QKD, several works introduce machine
learning-based methodologies to the CV QKD area, such as developing a novel CV
QKD scheme Jin et al. (2021); Liao et al. (2020), parameter prediction Liu et
al. (2018) and detecting quantum attacks Mao et al. (2020).

CV QKD protocols with discrete modulation have attracted increasing attention
for decades. Its appealing advantages include easier experimental
implementation and higher error-correction efficiencies which promote CV QKD
over longer distances Xu et al. (2020); Leverrier and Grangier (2009, 2011);
Zhao et al. (2009). These properties bring potential advantages in large-scale
deployment in quantum-secured networks Simon (2017). However, the security
analysis of discrete-modulated CV QKD protocols is more complicated owing to
the lack of symmetry Coles et al. (2016). Recently, some novel numerical
approaches Lin et al. (2019); Winick et al. (2018) have been proposed to
analyze the security of discrete-modulation protocols against collective
attacks, where key rate calculation involves minimizing a convex function over
all eavesdropping attacks that are consistent with the experimental data.
These numerical approaches achieve much higher key rates over significantly
longer distances compared with previous security analyses. Based on these
numerical approaches, a neural network model was presented to quickly predict
the secure key rate of discrete-modulated CV QKD with high reliability (secure
probability as high as $99.2\%$). This neural network model learns the mapping
between input parameters and key rates from datasets generated by numerical
methods, which supports the computation of secure key rates in real time Zhou
et al. (2021). However, the mapping complexity between input parameters and
key rates depends on the solving complexity of discrete-modulated protocols’
key rates through numerical approaches Hu et al. (2021). Selecting
architectures and hyperparameters plays a critical role in the performance of
a neural network. Therefore, to learn different mappings from different
protocols, the architectures of neural networks and the corresponding
hyperparameters should be adjusted carefully by humans, which comes at a great
price Yu and Zhu (2020).

Here, we propose a more flexible and automatic neural network model combined
with Bayesian optimization Shahriari et al. (2015), which maintains extremely
high reliability and efficiency and reduces complicated manual adjustment. Our
method is universal for a variety of unstructured QKD protocols that lack
analytical tools and rely on numerical methods. We apply our model to two
variants of discrete-modulated CV QKD protocols and acquire high secure key
rates with considerable tightness in both cases. We then compare the time
consumption of our model with the numerical method proposed in Ref. Lin et al.
(2019), which shows a great speedup of approximately $10^{7}$.

This paper is organized as follows. In section II, we introduce the numerical
method for CV QKD with discrete modulation proposed in Ref. Lin et al. (2019),
and we rely on it to collect a dataset to train and test the model. In section
III, we introduce more details about the Bayesian optimization used in this
paper. In section IV, we demonstrate all the main results of this paper.
Section V provides a discussion and concludes this paper.

##  II Numerical method for CV QKD with discrete modulation

In this work, we apply the model in two discrete-modulated CV QKD protocols
with different detection techniques to demonstrate the generalizability of our
model. One is the quadrature phase-shift-keying (QPSK) heterodyne detection
protocol Lin et al. (2019), and the other is an improved QPSK homodyne
detection protocol Liu et al. (2021). To collect a dataset for training neural
networks, we generate secure key rates of both protocols by applying the same
numerical method Lin et al. (2019); Hu et al. (2021). In the following, we
briefly introduce how computing key rates can be transformed into a relevant
convex objective function for numerical optimization. A more detailed
description can be found in Ref.Lin et al. (2019).

Here, we consider a CV QKD protocol with quaternary modulation that involves
two parties: a sender Alice and a receiver Bob. During each time in an
iteration of N rounds, Alice randomly prepares one of the four coherent states
$\left|\alpha_{k}\right\rangle=\left||\alpha|e^{i(2k\pi/4+\pi/4)}\right\rangle$,
where $k\in\\{0,1,2,3\\}$, and sends it to Bob via an untrusted quantum
channel. Then, Bob uses either homodyne or heterodyne detection to estimate k.
The secret key rate under collective attacks in the asymptotic limit is given
by the following expression according to the Devetak-Winter formula Devetak
and Winter (2005)

$R^{\infty}=p_{\text{pass }}[\min_{\rho\in\mathbf{S}}H(\mathbf{Z}\mid
E)-\delta_{\mathrm{EC}}]$ (1)

where $H\left(\mathbf{Z}\mid E\right)$ is conditional von Neumann entropy,
which describes the uncertainty of the string $\mathbf{Z}$ in Eve’s view.
Eve’s maximal knowledge of Bob’s string $\mathbf{Z}$ requires the minimum
uncertainty of $\mathbf{Z}$ under a certain density matrix $\rho$. Therefore,
we need to find the optimum $\rho^{*}$ in feasible domain $\mathbf{S}$ to
minimize $H\left(\mathbf{Z}\mid E\right)$, $p_{\text{pass }}$ is the sifting
probability, and $\delta_{\mathrm{EC}}$ is the actual amount of information
leakage per signal in the error-correction step. To turn this problem into a
convex optimization problem, the above expression can be reformulated as

$R^{\infty}=\min_{\rho_{AB}\in\mathrm{S}}D\left(\mathcal{G}\left(\rho_{AB}\right)\|\mathcal{Z}\left[\mathcal{G}\left(\rho_{AB}\right)\right]\right)-p_{\text{pass
}}\delta_{\mathrm{EC}}$ (2)

in which
$D(\rho\|\sigma)=\operatorname{Tr}\left(\rho\log_{2}\rho\right)-\operatorname{Tr}\left(\rho\log_{2}\sigma\right)$.
As shown in Ref.Winick et al. (2018), $\mathcal{G}$ is a completely positive
and trace nonincreasing map that describes the postprocessing of different
quadratures. $\mathcal{Z}$ is a pinching quantum channel that reads out the
key information.

Since the term $p_{\text{pass }}\delta_{\mathrm{EC}}$ in formula
$\left(2\right)$ is easy to compute, we can only consider the following
relevant optimization problem:

$\displaystyle\operatorname{minimize}D\left(\mathcal{G}\left(\rho_{AB}\right)\|\mathcal{Z}\left[\mathcal{G}\left(\rho_{AB}\right)\right]\right)$
(3) subject to
$\displaystyle\operatorname{Tr}\left[\rho_{AB}\left(|k\rangle\left\langle\left.k\right|_{A}\otimes\hat{q}\right)\right]=p_{k}\langle\hat{q}\rangle_{k}\right.$
(4)
$\displaystyle\operatorname{Tr}\left[\rho_{AB}\left(|k\rangle\left\langle\left.k\right|_{A}\otimes\hat{p}\right)\right]=p_{k}\langle\hat{p}\rangle_{k}\right.$
(5)
$\displaystyle\operatorname{Tr}\left[\rho_{AB}\left(|k\rangle\left\langle\left.k\right|_{A}\otimes\hat{n}\right)\right]=p_{k}\langle\hat{n}\rangle_{k}\right.$
(6) $\displaystyle\operatorname{Tr}[\rho_{AB}(|k\rangle\left\langle
k\right|_{A}\otimes\hat{d})]=p_{k}\langle\hat{d}\rangle_{k}$ (7)
$\displaystyle\operatorname{Tr}\left[\rho_{AB}\right]=1$ (8)
$\displaystyle\rho_{AB}\geq 0$ (9)
$\displaystyle\operatorname{Tr}_{B}\left[\rho_{AB}\right]=\sum_{i,j=0}^{3}\sqrt{p_{i}p_{j}}\left\langle\alpha_{j}\mid\alpha_{i}\right\rangle|i\rangle\left\langle\left.j\right|_{A}\right.$
(10)

where $k\in\left\\{0,1,2,3\right\\}$, $\langle\hat{q}\rangle_{k}$,
$\langle\hat{p}\rangle_{k}$, $\langle\hat{n}\rangle_{k}$ and
$\langle\hat{d}\rangle_{k}$ denote the expectation values of corresponding
operators when Bob measures states labeled by $k$. These expectation values
can be obtained through homodyne or heterodyne measurements. These first four
constraints come from experimental outcomes. The next two constraints are
natural requirements since $\rho_{AB}$ is a density matrix. The last
constraint on the partial trace of system B comes from the fact that the
quantum channel cannot influence system A of Alice. We can handle the above
density matrix and operators in finite dimensions $N_{c}$ after imposing the
photon-number cutoff assumption on this optimization problem Ghorai et al.
(2019); Lin et al. (2019). Then, this problem can be solved numerically.
Eventually, we solve this minimization problem by the numerical method
proposed in Ref.Winick et al. (2018). The specific implementation of this
numerical method in our work can be found in Ref.Liu et al. (2021). This
method involves two steps:

1\. Find a solution that is close to optimal, which gives an upper bound on
the key rate.

2\. Convert this upper bound to a lower bound on the key rate by considering
its dual problem.

\begin{overpic}[width=256.0748pt]{f1.pdf} \end{overpic} Figure 1: Illustration
of the Bayesian optimization procedure. Bayesian optimization estimates the
true objective function with a probability model called a surrogate. The blue
real curve represents the true objective function. Red cross points are
sampling points for the true objective function $f(x)$. The black dotted curve
and purple shadow represent the mean and confidence intervals estimated with
the surrogate of the objective function.

##  III Bayesian optimization

In this section, we present a brief introduction to Bayesian optimization.
Bayesian optimization is a powerful strategy for global optimization of
objective functions that are expensive to evaluate Shahriari et al. (2015);
Bergstra et al. (2011). This method is gaining great popularity in
hyperparameter optimization. In particular, hyperparameter optimization in
machine learning can be represented as follows:

$\displaystyle x^{\star}=\arg\min_{x\in\mathcal{X}}f(x),$ (11)

where $f(x):\mathcal{X}\rightarrow\mathbb{R}$ is an objective function to
minimize, $x^{\star}$ is a hyperparameter vector yielding the lowest value of
$f$, and the dimension of domain $\mathcal{X}$ depends on the total type of
concerned hyperparameters. In practice, the evaluation of the objective
function is extremely costly, which leads to selecting proper hyperparameters
by hand becoming intractable. Beyond the manual tuning method, grid search and
random search Bergstra and Bengio (2012) are two common methods that perform
slightly better. However, these methods still waste a large amount of time
evaluating poor hyperparameters across the entire search space, which is
relatively inefficient. In contrast, Bayesian optimization estimates the true
objective function with a probability model. Then, it utilizes Bayes’ theorem
to update this model based on previous results and chooses the next promising
hyperparameters. In practice, this method can find better hyperparameters in
less time. Figure 1 illustrates the Bayesian optimization procedure.

Algorithm 1 Sequential Model-Based Optimization

1:$\mathcal{H}_{0}\leftarrow\emptyset$

2:for $n=1,2,\ldots,N$ do

3: $x_{n+1}=\arg\max_{x}\alpha_{n}(x,S_{n})$

4: evaluate $y_{n+1}=f(x_{n+1})$

5: update $\mathcal{H}_{n+1}=\\{\mathcal{H}_{n},(x_{n+1},y_{n+1})\\}$

6: fit a new model $S_{n}$ to $\mathcal{H}_{n+1}$

7:end for

8:Return $\mathcal{H}_{N}$

Sequential model-based optimization (SMBO) algorithms are formalizations of
Bayesian optimization. Bergstra et al. (2011) These algorithms have two key
ingredients:

1\. A probabilistic surrogate model $S$. SMBO approximates the objective
function $f$ with a probabilistic model called a surrogate, which is cheaper
to evaluate. This surrogate contains a prior distribution capturing beliefs
about the behavior of the objective function and is then updated sequentially
after each new trial.

2\. An acquisition function $\alpha:\mathcal{X}\rightarrow\mathbb{R}$. The
acquisition function is the criterion by which the next vector of
hyperparameters is chosen from the surrogate function.

For an SMBO algorithm at iteration n, the next location $x_{n+1}$ is selected
by optimizing $\alpha_{n}$ and to evaluate the true $f$ to obtain a result
$y_{n+1}=f(x_{n+1})$. The new tuple $(x_{n+1},y_{n+1})$ is appended to the
historical set $\mathcal{H}$. Then, the surrogate model $S$ is updated
incorporating the new results, which means that the prior is updated to
produce a more informative posterior distribution over the space of objective
functions. The pseudocode of this framework is summarized in Algorithm 1.

Figure 2: Schematic diagram of our neural network model combined with Bayesian
optimization. The dataset training neural network to predict key rate is
generated by some numerical approach. Here, the hyperparameters related to the
neural network’s architecture are not determined by humans but updated by
Bayesian optimization. Bayesian optimization primarily establishes a
probability model for the distribution of neural network hyperparameters and
performance. Then, according to the evaluation metric produced in each trial,
such as validation loss, Bayesian optimization updates the probabilistic
surrogate model and suggests the next choice of hyperparameters. After several
trials, we can automatically obtain the best-performing neural network. Table
1: Hyperparameter search space of the neural network under the TPE algorithm
for the QPSK heterodyne detection protocol Lin et al. (2019). The neural
network model we use here is a fully connected forward network. By fixing the
number of neurons in the input layer and output layer, we search this neural
network’s architecture for hidden layers. For each of the three hidden layers,
there are the number of neurons, activation function in this layer and the
ratio of dropout layer following it waiting to determine, where the dropout
technique Srivastava et al. (2014) is used to prevent overfitting. The batch
size of the training process and two essential hyperparameters $\gamma$ and
$\varepsilon$ are searched as well. The brace $\\{\\}$ refers to a finite set
that contains all possible discrete values. Bracket $()$ represents a
continuous range.

| Number of neurons | Activation function | Dropout | Batch size | $\gamma$ | $\varepsilon$
---|---|---|---|---|---|---
Input layer | 29(fixed) | - | - | $\left\\{32,64,128,256\right\\}$ | $\left(0.05,0.2\right)$ | $\left(0.8,0.95\right)$
Hidden layer 1 | $\left\\{512,1024\right\\}$ | {tanh,ReLU,sigmoid} | $\left(0,0.3\right)$ |  |  |
Hidden layer 2 | $\left\\{128,256,512\right\\}$ | {tanh,ReLU,sigmoid} | $\left(0,0.3\right)$ |  |  |
Hidden layer 3 | $\left\\{128,256,512\right\\}$ | {tanh,ReLU,sigmoid} | $\left(0,0.3\right)$ |  |  |
Output layer | 1(fixed) | Linear(fixed) | - |  |  |

Table 2: Hyperparameter search space of the neural network under the TPE
algorithm for the QPSK homodyne detection protocol Liu et al. (2021).
Different from the QPSK heterodyne detection protocol Lin et al. (2019), here,
we search the number of hidden layers in $3$ or $4$.

| Number of neurons | Activation function | Dropout | Batch size | $\gamma$ | $\varepsilon$
---|---|---|---|---|---|---
Input layer | 29(fixed) | - | - | $\left[32,64,128,256\right]$ | $\left(0.05,0.2\right)$ | $\left(0.8,0.99\right)$
Hidden layer 1 | $\left\\{128,256,512\right\\}$ | [tanh,ReLU,sigmoid] | $\left(0,0.3\right)$ |  |  |
Hidden layer 2 | $\left\\{128,256,512\right\\}$ | [tanh,ReLU,sigmoid] | $\left(0,0.3\right)$ |  |  |
Hidden layer 3 | $\left\\{128,256\right\\}$ | [tanh,ReLU,sigmoid] | $\left(0,0.3\right)$ |  |  |
Hidden layer 4(optional) | $\left\\{64,128\right\\}$ | [tanh,ReLU,sigmoid] | $\left(0,0.3\right)$ |  |  |
Output layer | 1(fixed) | Linear(fixed) | - |  |  |

The most common choice of acquisition function is expected improvement (EI):

$\mathrm{E}\mathrm{I}_{y^{*}}(x):=\int_{-\infty}^{\infty}\max\left(y^{*}-y,0\right)p_{S}\left(y\mid
x\right)dy$ (12)

Here $y^{*}$ is a threshold value of the objective function $f$, and
$p_{S}\left(y\mid x\right)$ represents the surrogate probability model. If
this expectation is positive, then the vector of hyprparameters $x$ is
expected to produce a better result than $y^{*}$. There are several different
strategies for constructing the surrogate model: a Gaussian process approach
Williams and Rasmussen (2006), random forests Breiman (2001) and a tree-
structured Parzen estimator(TPE) Bergstra et al. (2011). In this work, the TPE
approach is adopted, which supports continuous, categorical and conditional
parameters, as well as priors for each hyperparameter over which values are
expected to perform best Hutter et al. (2015). In contrast, the Gaussian
process approach and random forests only support one or two types of the above
parameters, which are not capable of our following task covering continuous,
categorical and conditional parameters. Instead of directly modeling
$p\left(y\mid x\right)$, this method models $p\left(x\mid y\right)$ using two
such densities over the configuration space $\mathcal{X}$:

$p\left(x\mid y\right)=\begin{cases}\ell(x)&\text{ if }y<y^{*}\\\ g(x)&\text{
if }y\geq y^{*}\end{cases}$ (13)

This algorithm chooses $y^{*}$ to be some quantile $\gamma$ of the observed y
values, which means $p\left(y<y^{*}\right)=\gamma$. So the
$EI_{y^{*}}(x)=\frac{\gamma
y^{*}\ell(x)-\ell(x)\int_{-\infty}^{y^{*}}p(y)dy}{\gamma\ell(x)+(1-\gamma)g(x)}\propto\left(\gamma+\frac{g(x)}{\ell(x)}(1-\gamma)\right)^{-1}$.
The tree-structured form of $\ell$ and $g$ makes it easy to draw many
candidates according to $g(x)/\ell(x)$. On each iteration, the algorithm
returns the candidate $x$ with the greatest EI. We implement this algorithm
for the hyperparameter optimization of the neural networks predicting CV QKD
key rates, by using a Python library called Hyperopt Bergstra et al. (2013).

##  IV Method

Artificial neural networks can approximate arbitrary bounded continuous
mapping on a given domain, according to universal approximation theorem Hornik
et al. (1989). Therefore, we expect that the neural network can learn the
mapping between input variables defined in the constraints of Eq.(3) and
output key rates, which avoids solving the time-consuming optimization problem
and computes key rates with low latency. We demonstrated this possibility of
using a neural network to predict the key rates of discrete-modulated CV QKD
in previous work Zhou et al. (2021). In that work, we built a four-layer fully
connected forward neural network holding a loss function designed specifically
to predict the key rates of discrete-modulated CV QKD with homodyne detection.
The objective loss function is the key ingredient to keep the output key rates
reliable and tight. We retain it in this work but utilize the TPE algorithm to
search other parts of the neural network to improve the network’s overall
performance. The specific formula of the loss function is as follows:

$\displaystyle\mathcal{L}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\gamma\left(e_{i}^{*2}+\max\left(e_{i}^{*},-\log_{10}(\varepsilon)\right)\right)$
(14) $\displaystyle-(1-\gamma)\left(\min\left(e_{i}^{*},0\right)\right)$

Table 3: Resulting structure of neural networks of QPSK heterodyne detection
protocol Lin et al. (2019)

| Number of neurons | Activation function | Dropout | Batch size | $\gamma$ | $\varepsilon$
---|---|---|---|---|---|---
Input layer | 29 | - | - | $64$ | $0.0539$ | $0.8727$
Hidden layer 1 | $1024$ | sigmoid | $0.0769$ |  |  |
Hidden layer 2 | $256$ | tanh | $0.0362$ |  |  |
Hidden layer 3 | $256$ | sigmoid | $0.0481$ |  |  |
Output layer | 1 | Linear | - |  |  |

Table 4: Resulting structure of neural networks of QPSK homodyne detection
protocol Liu et al. (2021)

| Number of neurons | Activation function | Dropout | Batch size | $\gamma$ | $\varepsilon$
---|---|---|---|---|---|---
Input layer | 29 | - | - | $64$ | $0.1227$ | $0.8784$
Hidden layer 1 | $512$ | tanh | $0.2210$ |  |  |
Hidden layer 2 | $128$ | tanh | $0.2361$ |  |  |
Hidden layer 3 | $256$ | sigmoid | $0.0657$ |  |  |
Hidden layer 4 | $128$ | tanh | $0.0036$ |  |  |
Output layer | 1 | Linear | - |  |  |

Figure 3: Relative deviation distribution of key rates predicted securely in
respective test set for two protocols: (a)Protocol 1(with heterogyne
detection) Lin et al. (2019). In these secure results, the relative deviations
falling in $\left[-20\%,0\right]$ account for $92.49\%$ and falling in
$\left[-40\%,0\right]$ account for $96.48\%$. (b)Protocol 2 with impoved
homodyne detection Liu et al. (2021). The relative deviations falling in
$\left[-20\%,0\right]$ account for $90.26\%$ and falling in
$\left[-40\%,0\right]$ account for $98.39\%$. Figure 4: Comparison between NN-
predicted results and numerical results in key rates for two protocols. (a)The
predicted key rates vs. transmission distance for the QPSK heterodyne
detection protocol Lin et al. (2019) for different values of the excess noise,
from top to bottom, $\xi=0.005,0.010,0.015,0.020,0.025,0.030,0.035,0.040$. The
signal state amplitude of the QPSK heterodyne detection protocol is optimized
in the range $\left[0.62,0.72\right]$ with a step of $0.01$. (b)The predicted
key rates vs. transmission distance for the QPSK homodyne detection protocol
Liu et al. (2021) for different values of excess noise, from top to bottom,
$\xi=0.015,0.020,0.025,0.030,0.035,0.040$. The amplitude of the QPSK homodyne
detection protocol is optimized in the range $\left[0.62,1.03\right]$ with a
step of $0.01$. (c) and (d)Relative deviation between NN-predicted results and
numerical results, respectively, for the QPSK heterodyne detection protocol
and QPSK homodyne detection protocol, in which the corresponding NN-predicted
results are shown in (a) and (b). For all points, we set reconciliation
efficiency $\beta=0.95$, postselection $\Delta=0$, transmittance
$\eta=10^{-\frac{aL}{10}}$ in the distance L with $a=0.2$ dB/km, photon-number
cutoff $N_{c}=12$ and the maximal $N_{i}=300$ iteration of the first step in
the numerical method.  Figure 5: Running time comparison between neural-
network method and numerical method. (a)Time comparison of protocol 1 for
excess noise $\xi=0.005,0.015,0.025$ is shown as circles, triangles and
diamonds respectively. (b)Time comparison of protocol 2 for excess noise
$\xi=0.015,0.025,0.035$ is shown as circles, triangles and diamonds
respectively. Each point refers to a ratio of time consumed by numerical
method and neural-network method. The values of $\beta$, $\Delta$, $\eta$,
$N_{c}$ and $N_{i}$ remain the same as those mentioned in Fig. 4.

For training inputs $\left\\{\vec{x}_{i}\right\\}$ and corresponding labels
$\left\\{{y}_{i}\right\\}$, here $n$ is the size of
$\left\\{\vec{x}_{i}\right\\}$, $e_{i}^{*}=y_{i}^{*p}-y_{i}^{*}$ is the
residual error between preprocessed label $y_{i}^{*}$ and the corresponding
output of the neural network $y_{i}^{*p}$, where
$y_{i}^{*}=-\log_{10}\left(y_{i}\right)$. There are two significant
hyperparameters $\gamma$ and $\varepsilon$ contained in this loss function,
the choices of which are crucial to a model’s performance, as we presented in
Ref.Zhou et al. (2021). The meaning of hyperparameter $\gamma$ is to force the
predicted key rate as information-theoretically secure as possible, and
$\varepsilon$ is to force the predicted key rate as close to numerical results
as possible. Here, apart from the input layer and output layer, we do not fix
the structure of the neural network but utilize the TPE algorithm to search it
efficiently in a set configuration space. An illustration of our model is
shown in Fig. 2.

##  V Result

After this training under TPE searching is complete, we obtain the resulting
structures of neural networks in both cases, which are shown in Table 3-4.
Then, we use the selected and trained network to predict key rates on the test
set for both protocols. The predicted key rates that show security achieve as
high as $99.15\%$ for the QPSK heterodyne detection protocol and $99.59\%$ for
the QPSK homodyne detection protocol, which suggests that our method combining
a neural network with Bayesian optimization is highly reliable. For those key
rates predicted securely, namely, predicted results are lower than the true
values, we plot their relative deviation distributions for both protocols in
Fig. 3. Figure 3 suggests that our method has good tightness.

Before training the neural network under the TPE method, we generate datasets
for two different protocols by the aforementioned numerical approach. To
obtain datasets with diversity, for the QPSK heterodyne detection protocol, we
generate $36$ sets of data from excess noise $\xi=0.0045-0.0405$. Each dataset
contains $80$ random samplings for $\xi$ from an interval of length $0.001$,
for example $\left[0.0045,0.0055\right]$. Under each random sampling, we
generate data every $5$km with the transmission distance $L$ up to $200$km. At
each distance, we generate data from amplitude $\alpha=0.62-0.72$ in a step of
$0.01$. The total datasets contain $809,600$ input instances
$\left\\{\vec{x}_{i}\right\\}$ and corresponding labels
$\left\\{{y}_{i}\right\\}$. For the QPSK homodyne detection protocol, excess
noise is sampled randomly from $\xi=0.014-0.042$, where the length of sampling
interval is $0.002$, for example $\left[0.014,0.016\right]$, and amplitude
$\alpha$ is sampled from $\left[0.60,1.05\right]$. The size of total datasets
is $368,116$. For both protocols, each
$\vec{x}_{i}\in\left\\{\vec{x}_{i}\right\\}$ represents a vector of $29$
variables, there are $16$ variables that are the right parts of the first four
constraints of Eq. 3, $12$ variables are nondiagonal elements of the right
side matrix of the last constraint of Eq. 3, and the remaining variable is
excess noise $\xi$. Label $y_{i}\in\left\\{\vec{y}_{i}\right\\}$ represents
the corresponding key rate. Therefore, we fix the neurons of the network’s
input layer in $29$ and output layer in $1$, and the search space of other
hyperparameters can be found in Table 1-2.

Before feeding data into neural networks, we split data into a training set
and a test set and implement data preprocessing as in Ref.Zhou et al. (2021).
For the QPSK heterodyne detection protocol, the training set contains
$769,120$ data instances, and the test set contains $40,480$ data instances.
For the QPSK homodyne detection protocol, the training set contains $327,636$
data instances, and the test set contains $17,244$ data instances. For both
cases, there is $10\%$ of the training data split as the validation set. We
generate the dataset on the blade cluster system of the High Performance
Computing Center of Nanjing University. We consume over $250,000$ core hours,
and the node we use contains 4 Intel Xeon Gold 6248 CPUs, which involves
immense computational power. Under the TPE algorithm with max iteration $10$,
the Adam algorithm Kingma and Ba (2014) is used to train neural networks for
$200$ epochs, and the initial learning rate is set to $0.001$. It takes
roughly $53$ hours for the QPSK heterodyne detection protocol and $23$ hours
for the QPSK homodyne detection protocol on an Nvidia A100 GPU.

Here, we also compare the predicted results with numerical results in key
rates versus transmission distance for two protocols. The comparison is shown
in Fig. 4. For this plot, we implement the same numerical approach to compute
the best key rates of two protocols for different excess noises by optimizing
the amplitude $\alpha$ of signal states in the range $\left[0.62,0.72\right]$
and $\left[0.62,1.03\right]$ with a step of $0.01$. The choice of the excess
noise range is consistent with the sampling interval of previous training
data. The photon-number cutoff $N_{c}$ is $12$, and the maximal iteration
number of the first step in the numerical approach $N_{i}=300$. We record the
corresponding $29$ variables producing the best key rates as neural networks’
inputs to predict key rates. As shown in Figs. 4(c) and (d), the predicted
results are all secure and remain tight with relative deviations between
$10\%$ and $20\%$ when the transmission distance is below $150$ km for both
protocols.

To show the efficiency of our method, we compare the running time between the
neural network method and the numerical method on a high-performance personal
computer with a 3.3 GHz AMD Ryzen 9 4900H 16 GB of RAM, as shown in Fig. 5.
The results suggest that the neural network method is generally 6-8 orders of
magnitude of the numerical methods. For example, when $\xi=0.025$, the
numerical method consumes approximately $850$ seconds to calculate the key
rate at $50$ km for the QPSK heterodyne detection protocol. When $\xi=0.035$,
the numerical method consumes approximately $1260$ seconds at $25$ km to
calculate the key rate for the QPSK homodyne detection protocol. However, we
can use a trained neural network to obtain results in approximately $0.0001$
seconds, which is almost real time.

##  VI Discussion and conclusion

To summarize, we develop a neural network model combined with Bayesian
optimization to directly extract key rates with high reliability, considerable
tightness and great efficiency. Beyond designing the neural network
architecture by human and troublesome manual tuning of hyperparameters, we
utilize a special Bayesian optimization method called the TPE algorithm to
automatically search the structure and hyperparameters that are the best fit
for a given dataset. We exemplify our method on two promising discrete-
modulated CV QKD protocols varied by different detection techniques across a
large range of excess noise and transmission distances. For both protocols,
the neural networks selected by the TPE algorithm predict the information-
theoretically secure key rates with great high probability(up to $99.15\%$ for
the QPSK heterodyne detection protocol and $99.59\%$ for the QPSK homodyne
detection protocol), and the results present considerable tightness.

We show that our method can achieve approximately $10^{7}$ faster than the
numerical method, which completely satisfies the requirement of the QKD system
in practice. In contrast, the numerical method takes several minutes to
calculate a point of key rate, which is intolerable since many free-space
sessions, such as satellite-ground or handheld QKD might have a window of only
minutes. While collecting enough data based on the numerical method to train
the model consumes a large amount of computing power, we can consider these
large computations offline. Once we obtain the trained neural network, it can
be deployed on a certain device to infer key rates online in milliseconds by
giving new inputs from the experiment. Ref.Wang and Lo (2019) demonstrated
that a neural network method for parameter optimization of QKD can be deployed
on various mobile low-power systems, which brings advantages of more power
efficiency and low latency. We can also forecast that our neural network
method combined with Bayesian optimization will play an essential role in
free-space QKD scenarios such as handheld Mélen et al. (2017), drone-based
Hill et al. (2017) or satellite-ground QKD Liao et al. (2017). Several works
have focused on machine learning for optimal parameters in QKD Wang and Lo
(2019); Liu et al. (2018); Lu et al. (2019); Ding et al. (2020). However, our
work predicts secure key rates directly by automatically designed neural
networks, which goes further than our previous work Zhou et al. (2021).

Based on our model, there are several directions worthy of investigation for
future work. Up to now, we have only covered computing the asymptotic key
rates. However, finite-size effects are practical issues considered in
discrete-modulated CV-QKD Leverrier et al. (2010). Note that a recent work has
analyzed the security and performance of discrete-modulated CV-QKD under a
finite-size scenario Almeida et al. (2021), which inspires us to improve our
model. To address these issues, we also consider applying our model to other
protocols in future work. Moreover, the issue of post-processing (notably the
error correction part) still limits the overall time acceleration for a
discrete-modulated continuous-variable QKD system. Note that the error
correction involving binary or quaternary error-correcting codes is less
complex compared with the situation of Gaussian modulation. Therefore, we also
consider developing an effective error-correction protocol for CV QKD with
discrete modulation using machine learning techniques in the future.

## Acknowledgments

We gratefully acknowledge the support from the Natural Science Foundation of
Jiangsu Province (No. BK20211145), the Fundamental Research Funds for the
Central Universities (No. 020414380182), the Key Research and Development
Program of Nanjing Jiangbei New Aera (No. ZDYD20210101), the Program for
Innovative Talents and Entrepreneurs in Jiangsu(No. JSSCRC2021484), and the
Key-Area Research and Development Program of Guangdong Province (No.
2020B0303040001). The authors would like to thank the High Performance
Computing Center of Nanjing University for the numerical calculations.

## References

  * Grigorescu et al. (2020) S. Grigorescu, B. Trasnea, T. Cocias, and G. Macesanu, Journal of Field Robotics 37, 362 (2020).
  * Levinson et al. (2011) J. Levinson, J. Askeland, J. Becker, J. Dolson, D. Held, S. Kammel, J. Z. Kolter, D. Langer, O. Pink, V. Pratt, et al., in 2011 IEEE intelligent vehicles symposium (IV) (IEEE, 2011), pp. 163–168.
  * Deng et al. (2013) L. Deng, G. Hinton, and B. Kingsbury, in 2013 IEEE international conference on acoustics, speech and signal processing (IEEE, 2013), pp. 8599–8603.
  * Young et al. (2018) T. Young, D. Hazarika, S. Poria, and E. Cambria, ieee Computational intelligenCe magazine 13, 55 (2018).
  * Jumper et al. (2021) J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko, et al., Nature 596, 583 (2021).
  * Davies et al. (2021) A. Davies, P. Veličković, L. Buesing, S. Blackwell, D. Zheng, N. Tomašev, R. Tanburn, P. Battaglia, C. Blundell, A. Juhász, et al., Nature 600, 70 (2021).
  * Gao et al. (2018) J. Gao, L.-F. Qiao, Z.-Q. Jiao, Y.-C. Ma, C.-Q. Hu, R.-J. Ren, A.-L. Yang, H. Tang, M.-H. Yung, and X.-M. Jin, Phys. Rev. Lett. 120, 240501 (2018).
  * Yang et al. (2019) M. Yang, C.-l. Ren, Y.-c. Ma, Y. Xiao, X.-J. Ye, L.-L. Song, J.-S. Xu, M.-H. Yung, C.-F. Li, and G.-C. Guo, Phys. Rev. Lett. 123, 190401 (2019).
  * Ahmed et al. (2021) S. Ahmed, C. S. Muñoz, F. Nori, and A. F. Kockum, Phys. Rev. Research 3, 033278 (2021).
  * Bukov et al. (2018) M. Bukov, A. G. Day, D. Sels, P. Weinberg, A. Polkovnikov, and P. Mehta, Phys. Rev. X 8, 031086 (2018).
  * Lumino et al. (2018) A. Lumino, E. Polino, A. S. Rab, G. Milani, N. Spagnolo, N. Wiebe, and F. Sciarrino, Phys. Rev. Appl. 10, 044033 (2018).
  * Niu et al. (2019) M. Y. Niu, S. Boixo, V. N. Smelyanskiy, and H. Neven, npj Quantum Inf. 5, 33 (2019).
  * Hentschel and Sanders (2011) A. Hentschel and B. C. Sanders, Phys. Rev. Lett. 107, 233601 (2011).
  * Bennett and Brassard (1984) C. H. Bennett and G. Brassard, in Conf. on Computers, Systems and Signal Processing (Bangalore, India (1984), vol. 175.
  * Ekert (1991) A. K. Ekert, Phys. Rev. Lett. 67, 661 (1991).
  * Shor and Preskill (2000) P. W. Shor and J. Preskill, Phys. Rev. Lett. 85, 441 (2000).
  * Lo et al. (2014) H.-K. Lo, M. Curty, and K. Tamaki, Nat. Photonics 8, 595 (2014).
  * Gisin et al. (2002) N. Gisin, G. Ribordy, W. Tittel, and H. Zbinden, Rev. Mod. Phys. 74, 145 (2002).
  * Grosshans and Grangier (2002) F. Grosshans and P. Grangier, Phys. Rev. Lett. 88, 057902 (2002).
  * Lance et al. (2005) A. M. Lance, T. Symul, V. Sharma, C. Weedbrook, T. C. Ralph, and P. K. Lam, Phys. Rev. Lett. 95, 180503 (2005).
  * Huang et al. (2015) D. Huang, P. Huang, D. Lin, C. Wang, and G. Zeng, Opt. Lett. 40, 3695 (2015).
  * Yin et al. (2019) H.-L. Yin, W. Zhu, and Y. Fu, Sci. Rep. 9, 49 (2019).
  * Fossier et al. (2009) S. Fossier, E. Diamanti, T. Debuisschert, A. Villing, R. Tualle-Brouri, and P. Grangier, New J. Phys. 11, 045023 (2009).
  * Huang et al. (2016) D. Huang, P. Huang, H. Li, T. Wang, Y. Zhou, and G. Zeng, Opt. Lett. 41, 3511 (2016).
  * Jouguet et al. (2012) P. Jouguet, S. Kunz-Jacques, T. Debuisschert, S. Fossier, E. Diamanti, R. Alléaume, R. Tualle-Brouri, P. Grangier, A. Leverrier, P. Pache, et al., Opt. Express 20, 14030 (2012).
  * Pirandola et al. (2020) S. Pirandola, U. L. Andersen, L. Banchi, M. Berta, D. Bunandar, R. Colbeck, D. Englund, T. Gehring, C. Lupo, C. Ottaviani, et al., Advances in Optics and Photonics 12, 1012 (2020).
  * Jin et al. (2021) D. Jin, Y. Guo, Y. Wang, Y. Li, and D. Huang, Phys. Rev. A 104, 012616 (2021).
  * Liao et al. (2020) Q. Liao, G. Xiao, H. Zhong, and Y. Guo, New J. Phys. 22, 083086 (2020).
  * Liu et al. (2018) W. Liu, P. Huang, J. Peng, J. Fan, and G. Zeng, Phys. Rev. A 97, 022316 (2018).
  * Mao et al. (2020) Y. Mao, W. Huang, H. Zhong, Y. Wang, H. Qin, Y. Guo, and D. Huang, New J. Phys. 22, 083073 (2020).
  * Xu et al. (2020) F. Xu, X. Ma, Q. Zhang, H.-K. Lo, and J.-W. Pan, Rev. Mod. Phys. 92, 025002 (2020).
  * Leverrier and Grangier (2009) A. Leverrier and P. Grangier, Phys. Rev. Lett. 102, 180504 (2009).
  * Leverrier and Grangier (2011) A. Leverrier and P. Grangier, Phys. Rev. A 83, 042312 (2011).
  * Zhao et al. (2009) Y.-B. Zhao, M. Heid, J. Rigas, and N. Lütkenhaus, Phys. Rev. A 79, 012307 (2009).
  * Simon (2017) C. Simon, Nat. Photonics 11, 678 (2017).
  * Coles et al. (2016) P. J. Coles, E. M. Metodiev, and N. Lütkenhaus, Nat. Commun. 7, 1 (2016).
  * Lin et al. (2019) J. Lin, T. Upadhyaya, and N. Lütkenhaus, Phys. Rev. X 9, 041064 (2019).
  * Winick et al. (2018) A. Winick, N. Lütkenhaus, and P. J. Coles, Quantum 2, 77 (2018).
  * Zhou et al. (2021) M.-G. Zhou, Z.-P. Liu, W.-B. Liu, C.-L. Li, J.-L. Bai, Y.-R. Xue, Y. Fu, H.-L. Yin, and Z.-B. Chen, arXiv preprint arXiv:2108.02578 (2021).
  * Hu et al. (2021) H. Hu, J. Im, J. Lin, N. Lütkenhaus, and H. Wolkowicz, arXiv preprint arXiv:2104.03847 (2021).
  * Yu and Zhu (2020) T. Yu and H. Zhu, arXiv preprint arXiv:2003.05689 (2020).
  * Shahriari et al. (2015) B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. De Freitas, Proceedings of the IEEE 104, 148 (2015).
  * Liu et al. (2021) W.-B. Liu, C.-L. Li, Y.-M. Xie, C.-X. Weng, J. Gu, X.-Y. Cao, Y.-S. Lu, B.-H. Li, H.-L. Yin, and Z.-B. Chen, PRX Quantum 2, 040334 (2021).
  * Devetak and Winter (2005) I. Devetak and A. Winter, Proceedings of the Royal Society A: Mathematical, Physical and engineering sciences 461, 207 (2005).
  * Ghorai et al. (2019) S. Ghorai, P. Grangier, E. Diamanti, and A. Leverrier, Phys. Rev. X 9, 021059 (2019).
  * Bergstra et al. (2011) J. Bergstra, R. Bardenet, Y. Bengio, and B. Kégl, Advances in neural information processing systems 24 (2011).
  * Bergstra and Bengio (2012) J. Bergstra and Y. Bengio, Journal of machine learning research 13 (2012).
  * Williams and Rasmussen (2006) C. K. Williams and C. E. Rasmussen, Gaussian processes for machine learning, vol. 2 (MIT press Cambridge, MA, 2006).
  * Breiman (2001) L. Breiman, Machine learning 45, 5 (2001).
  * Hutter et al. (2015) F. Hutter, J. Lücke, and L. Schmidt-Thieme, KI-Künstliche Intelligenz 29, 329 (2015).
  * Bergstra et al. (2013) J. Bergstra, D. Yamins, and D. D. Cox, in Proceedings of the 12th Python in science conference (Citeseer, 2013), vol. 13, p. 20.
  * Srivastava et al. (2014) N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, The journal of machine learning research 15, 1929 (2014).
  * Hornik et al. (1989) K. Hornik, M. Stinchcombe, and H. White, Neural networks 2, 359 (1989).
  * Kingma and Ba (2014) D. P. Kingma and J. Ba, arXiv preprint arXiv:1412.6980 (2014).
  * Wang and Lo (2019) W. Wang and H.-K. Lo, Phys. Rev. A 100, 062334 (2019).
  * Mélen et al. (2017) G. Mélen, P. Freiwang, J. Luhn, T. Vogl, M. Rau, C. Sonnleitner, W. Rosenfeld, and H. Weinfurter, in Quantum Information and Measurement (Optical Society of America, 2017), pp. QT6A–57.
  * Hill et al. (2017) A. D. Hill, J. Chapman, K. Herndon, C. Chopp, D. J. Gauthier, and P. Kwiat, Urbana 51, 61801 (2017).
  * Liao et al. (2017) S.-K. Liao, W.-Q. Cai, W.-Y. Liu, L. Zhang, Y. Li, J.-G. Ren, J. Yin, Q. Shen, Y. Cao, Z.-P. Li, et al., Nature 549, 43 (2017).
  * Lu et al. (2019) F.-Y. Lu, Z.-Q. Yin, C. Wang, C.-H. Cui, J. Teng, S. Wang, W. Chen, W. Huang, B.-J. Xu, G.-C. Guo, et al., JOSA B 36, B92 (2019).
  * Ding et al. (2020) H.-J. Ding, J.-Y. Liu, C.-M. Zhang, and Q. Wang, Quantum Inf. Process. 19, 1 (2020).
  * Leverrier et al. (2010) A. Leverrier, F. Grosshans, and P. Grangier, Phys. Rev. A 81, 062343 (2010).
  * Almeida et al. (2021) M. Almeida, D. Pereira, N. J. Muga, M. Facão, A. N. Pinto, and N. A. Silva, Opt. Express 29, 38669 (2021).


11institutetext: Tianshu Luo, Yimu Guo 22institutetext: Institute of Solid
Mechanics, Department of Applied Mechanics, Zhejiang University,
Hangzhou, Zhejiang, 310027, P.R.China 22email<EMAIL_ADDRESS>33institutetext: Yimu Guo 44institutetext: Institute of Solid Mechanics,
Department of Applied Mechanics, Zhejiang University,
Hangzhou, Zhejiang, 310027, P.R.China 44email<EMAIL_ADDRESS>
# A Sort of Relation Between a Dissipative Mechanical System and Conservative
Ones

Tianshu Luo  Yimu Guo

(Received: date / Accepted: date)

###### Abstract

In this paper we proposed a proposition: for any nonconservative classical
mechanical system and any initial condition, there exists a conservative one;
the two systems share one and only one common phase curve; the Hamiltonian of
the conservative system is the sum of the total energy of the nonconservative
system on the aforementioned phase curve and a constant depending on the
initial condition. Hence, this approach entails substituting an infinite
number of conservative systems for a dissipative mechanical system
corresponding to varied initial conditions. One key way we use to demonstrate
these viewpoints is that by the Newton-Laplace principle the nonconservative
force can be reasonably assumed to be equal to a function of a component of
generalized coordinates $q_{i}$ along a phase curve, such that a
nonconservative mechanical system can be reformulated as countless
conservative systems. Utilizing the proposition, one can apply the method of
Hamiltonian mechanics or Lagrangian mechanics to dissipative mechanical
system. The advantage of this approach is that there is no need to change the
definition of canonical momentum and the motion is identical to that of the
original system.

###### Keywords:

Hamiltonian, dissipation, non-conservative system, damping, symplectic
algorithm

††journal: China Science Bulletein

##  1 Introduction

In general, Hamiltonian mechanics and Lagrangian mechanics are applied to
conservative classical mechanical system or conservative quantum-mechanical
system. In this paper we attempt to find a sort of relationship between a
dissipative classical mechanical system between nonconservative classical
mechanical ones, then we might apply some methods derived from symplectic
geometry to dissipative classical mechanical system.

Some researchers attempt to represent a dissipative system as Hamiltonian
formalism or Lagrangian formalism. For instances, about half a century ago,
CalirolaPCalirola1941 ,Kana1948PThPh…3..440K  adopted the Hamiltonian

$H_{ck}(q,p)=\frac{1}{2}\left(e^{-2\eta t}p^{2}+e^{2\eta
t}\omega^{2}q^{2}\right),$ (1)

which leads exactly to the classical equation of motion of a damped harmonic
oscillator,

$\ddot{x}+2\eta\dot{x}+\omega^{2}x^{2}=0,\ \ \eta>0$ (2)

In this Hamiltonian-description, the canonical momentum is defined as

$p_{ck}=e^{2\eta t}p$

In 1940s Morse and Feshbachbook3  gave an example of an artificial Hamiltonian
for a damped oscillator based on a “mirror-image” trick, incorporating a
second oscillator with negative friction. The resulting Hamiltonian is
unphysical: it is unbounded from below and under time reversal the oscillator
is transformed into its “mirror-image”. By this arbitrary trick dissipative
systems can be handled as though they were conservative. BatemanPhysRev.38.815
proposed a similar approach. For the system (2), we have

$\displaystyle\ddot{x}+2\eta\dot{x}+\omega^{2}x^{2}=0\ \ (original)$ (3)
$\displaystyle\ddot{y}-2\eta\dot{x}+\omega^{2}x^{2}=0\ \ (mirror-image).$ (4)

Correspondingly, there is Bateman(-Morse-Feshbach) Lagrangian:

$L_{B}(x,\dot{x},y,\dot{y})=\dot{x}\dot{y}+\eta(x\dot{y}-\dot{x}y)-m\omega^{2}xy$
(5)

Rajeevart7  considered that a large class of dissipative systems can be
brought to a canonical form by introducing complex coordinates in phase space
and a complex-valued Hamiltonian. Rajeevart7  indicated that Eq.(2) can be
brought to diagonal form by a linear transformation:

$z=A\left[-\mathrm{i}(p+\eta
x)+\omega_{1}x\right],\frac{\mathrm{d}z}{\mathrm{d}t}=\left[-\gamma+\mathrm{i}\omega_{1}\right]z,$
(6)

where

$\omega_{1}=\sqrt{\omega_{2}-\gamma^{2}},$ (7)

and the constant $A=1/\sqrt{2\omega_{1}}$. Then art7  defined the complex-
valued function as Hamiltonian

$\mathcal{H}=(\omega_{1}+\mathrm{i}\eta)zz^{*},$ (8)

which satisfied

$\frac{\mathrm{d}z}{\mathrm{d}t}=\left\\{\mathcal{H},z\right\\},\ \
\frac{\mathrm{d}z^{*}}{\mathrm{d}t}=\left\\{\mathcal{H},z^{*}\right\\}$

By reviewing the works of PCalirola1941 1948PThPh…3..440K PhysRev.38.815 book3
art7 , we can find that they attempt to transform a dissipative system into a
conservative system entirely and these approaches might be suitable for
Hamiltonian representation of one-dimensional damped oscillators (weak non-
Lagrangian systems) and quantization. Because by observing Eq.(1), Eq.(4),
Eq.(5) and the transformation (6), one can find that the damping coefficient
is independent of other particles, and art7  had wrote: ’ These complex
coordinates are the natural variable(normal modes) of the system. ’

In area of quantum mechanics, PhysRevA.81.022112 ; Kochan2010219  attempts to
quantize dissipative forces in terms of the two form $\Omega$ ( an analog of
$\mathrm{d}p\wedge\mathrm{d}q-\mathrm{d}H\wedge\mathrm{d}t$), avoiding to
obtain Hamiltonian or Lagrangian formulation of non-Lagrangian system.

Marsden Marsden2007  and other researchers applied the equations as below to
the problem of stability of dissipative systems

$\displaystyle\dot{p}_{i}$ $\displaystyle=$ $\displaystyle-\frac{\partial
H}{\partial q_{i}}+\bm{F}\left(\frac{\partial{r}}{\partial q_{i}}\right)$
$\displaystyle\dot{q}_{i}$ $\displaystyle=$ $\displaystyle\frac{\partial
H}{\partial p_{i}},$ (9)

where the position vector $r$ depends on the canonical variable $\\{q,p\\}$,
i.e. $r(q,p)$, $H$ denotes Hamiltonian, and
$\bm{F}(\partial{r}/\partial{q_{i}})$ denotes a generalized force in the
direction $i$, $i=1,\dots,n$. Marsden considered that Eqs.(9) was composed of
a conservative part and a non-conservative part. Eq.(9) apparently is a
representation of dissipative mechanical systems in the phase space. Although
one can utilize the approaches discussed in some papersPCalirola1941
1948PThPh…3..440K PhysRev.38.815 book3 art7  to convert Eq.(9) into a
conservative system, one must first change the definition of the canonical
momentum of the system. If one uses numerical algorithms to solve the
Hamiltonian system, the numerical solution will lose the physical
characteristics of the original system, because the phase flow of the original
system is different from that of the new system. We need a Hamiltonian system
that shares common phase flow or solution with the original system. But this
demand cannot be satisfied, because it conflicts with Louisville’s theorem.
Therefore, we would have to attempt to find other relationship between
dissipative systems and conservative ones.

Based on Eq.(9), in this paper we will attempt to demonstrate that a
dissipative mechanical system shares a single common phase curve with a
conservative system. In the light of this property, we will propose an
approach to substitute a group of conservative systems for a dissipative
mechanical system. In the following section, we will illustrate the
relationship between a dissipative mechanical system and a conservative one.

##  2 Relationship between a Dissipative Mechanical System and a Conservative
One

###  2.1 A Proposition

Under general circumstances, the force $\bm{F}$ is a damping force that
depends on the variable set
$q_{1},\cdots,q_{n},\dot{q}_{1},\cdots,\dot{q}_{n}$. $F_{i}$ denotes the
components of the generalized force $\bm{F}$.

$F_{i}(q_{1},\cdots,q_{n},\dot{q}_{1},\cdots,\dot{q}_{n})=\bm{F}\left(\frac{\partial{r}}{\partial
q_{i}}\right).$ (10)

Thus we can reformulate Eq.(9) as follows:

$\displaystyle\dot{p}_{i}$ $\displaystyle=$ $\displaystyle-\frac{\partial
H}{\partial q_{i}}+F_{i}(q_{1},\cdots,q_{n},\dot{q}_{1},\cdots,\dot{q}_{n})$
$\displaystyle\dot{q}_{i}$ $\displaystyle=$ $\displaystyle\frac{\partial
H}{\partial p_{i}}.$ (11)

Suppose the Hamiltonian quantity of a conservative system without damping is
$\hat{H}$. Thus we may write a Hamilton’s equation of the conservative system
:

$\displaystyle\dot{p}_{i}$ $\displaystyle=$
$\displaystyle-\frac{\partial{\hat{H}}}{\partial q_{i}}$
$\displaystyle\dot{q}_{i}$ $\displaystyle=$
$\displaystyle\frac{\partial\hat{H}}{\partial p_{i}}.$ (12)

We do not intend to change the definition of momentum in classical mechanics,
but we do require that a special solution of Eq.(12) is the same as that of
Eq.(11). We may therefore assume a phase curve $\gamma$ of Eq.(11) coincides
with that of Eq.(12). The phase curve $\gamma$ corresponds to an initial
condition $q_{i0},p_{i0}$. Consequently by comparing Eq.(11) and Eq.(12), we
have

$\displaystyle\left.\frac{\partial{\hat{H}}}{\partial{q_{i}}}\right|_{\gamma}$
$\displaystyle=$ $\displaystyle\left.\frac{\partial H}{\partial
q_{i}}\right|_{\gamma}-\left.F_{i}(q_{1},\cdots,q_{n},\dot{q}_{1},\cdots,\dot{q}_{n})\right|_{\gamma}$
$\displaystyle\left.\frac{\partial{\hat{H}}}{\partial{p_{i}}}\right|_{\gamma}$
$\displaystyle=$ $\displaystyle\left.\frac{\partial H}{\partial
p_{i}}\right|_{\gamma},$ (13)

where
$\left.\frac{\partial{\hat{H}}}{\partial{q_{i}}}\right|_{\gamma},\left.\frac{\partial
H}{\partial
q_{i}}\right|_{\gamma},\left.\frac{\partial{\hat{H}}}{\partial{p_{i}}}\right|_{\gamma}and\left.\frac{\partial
H}{\partial p_{i}}\right|_{\gamma}$ denote the values of these partial
derivatives on the phase curve $\gamma$ and
$\left.F_{i}(q_{1},\cdots,q_{n},\dot{q}_{1},\cdots,\dot{q}_{n})\right|_{\gamma}$
denotes the value of the force $F_{i}$ on the phase curve $\gamma$. In
classical mechanics the Hamiltonian $H$ of a conservative mechanical system is
mechanical energy and can be written as:

$H=\int_{\gamma}\left(\frac{\partial{H}}{\partial{q_{i}}}\right)\mathrm{d}q_{i}+\int_{\gamma}\left(\frac{\partial
H}{\partial p_{i}}\right)\mathrm{d}p_{i}+const_{1},$ (14)

where $const_{1}$ is a constant that depends on the initial condition
described above. If $q_{i}=0,p_{i}=0$, then $const_{1}=0$. The mechanical
energy $H$ of the system (11) can be evaluated via Eq. (14) too. The Einstein
summation convention has been used this section. Thus an attempt has been made
to find $\left.\hat{H}\right|_{\gamma}$ through line integral along the phase
curve $\gamma$ of the dissipative system

$\displaystyle\int_{\gamma}\frac{\partial{\hat{H}}}{\partial{q_{i}}}\mathrm{d}q_{i}$
$\displaystyle=$ $\displaystyle\int_{\gamma}\left[\frac{\partial H}{\partial
q_{i}}-F_{i}(q_{1},\cdots,q_{n},\dot{q}_{1},\cdots,\dot{q}_{n})\right]\
\mathrm{d}q_{i}$ $\displaystyle\int_{\gamma}\frac{\partial\hat{H}}{\partial
p_{i}}\mathrm{d}p_{i}$ $\displaystyle=$
$\displaystyle\int_{\gamma}\frac{\partial H}{\partial p_{i}}\mathrm{d}p_{i}.$
(15)

Analogous to Eq.(14), we have

$\left.\hat{H}\right|_{\gamma}=\int_{\gamma}\frac{\partial{\hat{H}}}{\partial{q_{i}}}\mathrm{d}q_{i}+\int_{\gamma}\frac{\partial{\hat{H}}}{\partial
p_{i}}\mathrm{d}p_{i}+const_{2},$ (16)

where $const_{2}$ is a constant which depends on the initial condition.
Substituting Eq.(14)(15) into Eq.(16), we have

$\left.\hat{H}\right|_{\gamma}=H-\int_{\gamma}F_{i}(q_{1},\cdots,q_{n},\dot{q}_{1},\cdots,\dot{q}_{n})\mathrm{d}q_{i}+const.$
(17)

where $const=const_{2}-const_{1}$, and $H=\left.H\right|_{\gamma}$ because $H$
is mechanical energy of the nonconservative system(11). According to the
physical meaning of Hamiltonian, $const_{1}$, $const_{2}$ and $const$ are
added into Eq.(14)(16)(17) respectively such that the integral constant
vanishes in the Hamiltonian quantity. ArnoldArnold1997  had presented the
Newton-Laplace principle of determinacy as, ’This principle asserts that the
state of a mechanical system at any fixed moment of time uniquely determines
all of its (future and past) motion.’ In other words, in the phase space the
position variable and the velocity variable are determined only by the time
$t$. Therefore, we can assume that we have already a solution of Eq.(11)

$\displaystyle q_{i}$ $\displaystyle=$ $\displaystyle q_{i}(t)$
$\displaystyle\dot{q_{i}}$ $\displaystyle=$ $\displaystyle\dot{q_{i}}(t),$
(18)

where the solution satisfies the initial condition. We can divide the whole
time domain into a group of sufficiently small domains and in these domains
$q_{i}$ is monotone, and hence we can assume an inverse function $t=t(q_{i})$.
If $t=t(q_{i})$ is substituted into the nonconservative force
$\left.F_{i}\right|_{\gamma}$, we can assume that:

$\left.F_{i}(q_{1}(t(q_{i})),\cdots,q_{n}(t(q_{i})),\dot{q}_{1}(t(q_{i})),\cdots,\dot{q}_{n}(t(q_{i})))\right|_{\gamma}=\mathcal{F}_{i}(q_{i}),$
(19)

where $\mathcal{F}_{i}$ is a function of $q_{i}$ alone. In Eq.(19) the
function $F_{i}$ is restricted on the curve $\gamma$, such that a new function
$\mathcal{F}_{i}(q_{i})$ yields. Thus we have

$\displaystyle\int_{\gamma}F_{i}\mathrm{d}q_{i}$ $\displaystyle=$
$\displaystyle\int_{q_{i0}}^{q_{i}}\mathcal{F}_{i}(q_{i})\mathrm{d}q_{i}=W_{i}(q_{i})-W_{i}(q_{i0}).$
(20)

According to Eq.(20) the function $\mathcal{F}_{i}$ is path independent, and
therefore $\mathcal{F}_{i}$ can be regarded as a conservative force. For that
Eq.(19) represents an identity map from the nonconservative force $F_{i}$ on
the curve $\gamma$ to the conservative force $\mathcal{F}_{i}$ which is
distinct from $F_{i}$. It must be noted, that Eq.(19) is tenable only on the
phase curve $\gamma$. Consequently the function form of $\mathcal{F}_{i}$
depends on the aforementioned initial condition; from other initial conditions
$\mathcal{F}_{i}$ with different function forms will yield.

According to the physical meaning of Hamiltonian, $const$ is added to Eq.(17)
such that the integral constant vanishes in Hamiltonian quantity. Hence
$const=-W_{i}(q_{i0})$. Substituting Eq.(20) and $const=-W_{i}(q_{i0})$ into
Eq.(17), we have

$\left.\hat{H}\right|_{\gamma}=H-W_{i}(q_{i})$ (21)

where $-W_{i}(q_{i})$ denotes the potential of the conservative force
$\mathcal{F}_{i}$ and $W_{i}(q_{i})$ is equal to the sum of the work done by
the nonconservative force $F$ and $const$. In Eq.(21) $\hat{H}$ and $H$ are
both functions of $q_{i},p_{i}$ and $W_{i}(q_{i})$ a function of $q_{i}$.
Eq.(21) and Eq.(17) can be thought of as a map from the total energy of the
dissipative system(11) to the Hamiltonian of the conservative system(12).
Indeed, $\left.\hat{H}\right|_{\gamma}$ and the total energy differ in the
constant $const=-W_{i}(q_{i0})$. When the conservative system takes a
different initial condition, if one does not change the function form of
$\left.\hat{H}\right|_{\gamma}$, one can consider
$\left.\hat{H}\right|_{\gamma}$ as a Hamiltonian quantity $\hat{H}$,

$\hat{H}=\left.\hat{H}\right|_{\gamma}=H-W_{i}(q_{i})$ (22)

and the conservative system(12) can be thought of as an entirely new
conservative system.

Based on the above, the following proposition is made:

###### Proposition 1

For any nonconservative classical mechanical system and any initial condition,
there exists a conservative one; the two systems share one and only one common
phase curve; the value of the Hamiltonian of the conservative system is equal
to the sum of the total energy of the nonconservative system on the
aforementioned phase curve and a constant depending on the initial condition.

###### Proof

First we must prove the first part of the Proposition 1, i.e. that a
conservative system with Hamiltonian presented by Eq.(22) shares a common
phase curve with the nonconservative system represented by Eq.(11). In other
words the Hamiltonian quantity presented by Eq.(22) satisfies Eq.(13) under
the same initial condition. Substituting Eq.(22) into the left side of
Eq.(13), we have

$\displaystyle\frac{\partial{\hat{H}(q_{i},p_{i})}}{\partial{q_{i}}}$
$\displaystyle=$ $\displaystyle\frac{\partial
H(q_{i},p_{i})}{\partial{q_{i}}}-\frac{\partial{W_{j}(q_{j})}}{\partial{q_{i}}}$
$\displaystyle\frac{\partial{\hat{H}(q_{i},p_{i})}}{\partial{p_{i}}}$
$\displaystyle=$ $\displaystyle\frac{\partial
H(q_{i},p_{i})}{\partial{p_{i}}}-\frac{\partial{W_{j}(q_{j})}}{\partial{p_{i}}}.$
(23)

It must be noted that although $q_{i}$ and $p_{i}$ are considered as distinct
variables in Hamilton’s mechanics, we can consider $q_{i}$ and $\dot{q_{i}}$
as dependent variables in the process of constructing of $\hat{H}$. At the
trajectory $\gamma$ we have

$\displaystyle\frac{\partial{{W_{j}(q_{j})}}}{\partial{q_{i}}}$
$\displaystyle=$
$\displaystyle\frac{\partial{(\int_{q_{j0}}^{q_{j}}\mathcal{F}_{j}(q_{j})\mathrm{d}q_{j}+W_{i}(q_{i0}))}}{\partial{q_{i}}}=\mathcal{F}_{i}(q_{i})$
$\displaystyle\frac{\partial{{W_{j}(q_{j})}}}{\partial{p_{i}}}$
$\displaystyle=0,$ (24)

where $\mathcal{F}_{i}(q_{i})$ is equal to the damping force $F_{i}$ on the
phase curve $\gamma$. Hence under the initial condition $q_{0},p_{0}$, Eq.(13)
is satisfied. As a result, we can state that the phase curve of Eq.(12)
coincides with that of Eq.(11) under the initial condition; and $\hat{H}$
represented by Eq.(22) is the Hamiltonian of the conservative system
represented by Eq.(12).

Then we must prove the second part of Proposition 1: the uniqueness of the
common phase curve.

We assume that eq.(12) shares two common phase curves, $\gamma_{1}$ and
$\gamma_{2}$, with eq.(11). Let a point of $\gamma_{1}$ at the time $t$ be
$z_{1}$, a point of $\gamma_{2}$ at the time $t$ $z_{2}$, and $g^{t}$ the
Hamiltonian phase flow of eq.(12). Suppose a domain $\Omega$ at $t$ which
contains only points $z_{1}$ and $z_{2}$, and $\Omega$ is not only a subset of
the phase space of the nonconservative system(11) but also that of the phase
space of the conservative system(12). Hence there exists a phase flow
$\hat{g}^{t}$ composed of $\gamma_{1}$ and $\gamma_{2}$, and $\hat{g}^{t}$ is
the phase flow of eq.(11) restricted by $\Omega$. According to the following
Louisville’s theoremArnold1978 :

###### Theorem 2.1

The phase flow of Hamilton’s equations preserves volume: for any region $D$ we
have

$volume\ of\ g^{t}D=volume\ of\ D$

where $g^{t}$ is the one-parameter group of transformations of phase space

$g^{t}:(p(0),q(0))\longmapsto:(p(t),q(t))$

$g^{t}$ preserves the volume of $\Omega$. This implies that the phase flow of
eq.(11) $\hat{g}^{t}$ preserves the volume of $\Omega$ too. But the system
(11) is not conservative, which conflicts with Louisville’s theorem; hence
only a phase curve of eq.(12) coincides with that of eq.(11).

∎∎

In the next section three examples is given to demonstrate Proposition 1.

##  3 Examples

In this section, first two simple analytical examples are given, then a pro
forma example is given.

###  3.1 One-dimensional Analytical Example

Consider a special one-dimensional simple mechanical system:

$\ddot{x}+c\dot{x}=0,$ (25)

where $c$ is a constant. The exact solution of the equation above is

$x=A_{1}+A_{2}e^{-ct},$ (26)

where $A_{1},A_{2}$ are constants. From the equation above, we derived the
velocity:

$\dot{x}=-cA_{2}e^{-ct}.$ (27)

From the initial condition $x_{0},\dot{x}_{0}$, we find
$A_{1}=x_{0}+\dot{x}_{0}/c,A_{2}=-\dot{x}_{0}/c$. From Eq.(26)

$t=-\frac{1}{c}\ln\frac{x-A_{1}}{A_{2}}$ (28)

Substituting the equation above into Eq.(27), such we have

$\dot{x}=-c(x-A_{1})=-c(x-A_{1})$ (29)

The dissipative force $F$ in the dissipative system (25) is

$F=c\dot{x}.$ (30)

Substituting Eq.(29) into Eq.(30), such we have the conservative force
$\mathcal{F}$

$\mathcal{F}=-c^{2}(x-A_{1});$ (31)

Clearly the conservative force $\mathcal{F}$ depends on the initial condition
of the dissipative system (25), in other words an initial condition determine
a conservative force. Consequently a new conservative system yields

$\ddot{x}+\mathcal{F}=0\rightarrow\ddot{x}-c^{2}(x-A_{1})=0.$ (32)

The stiffness coefficient of the equation above must be negative. One can
readily verify that the particular solution (26) of the dissipative system can
satisfy the conservative one (32). This point agrees with Proposition (1).

The potential of the conservative system（32）is

$V=\int_{0}^{x}\left[-c^{2}(x-A_{1})\right]=-\frac{c^{2}}{2}x^{2}+c^{2}A_{1}x$

If $t\rightarrow\infty$，$x\rightarrow A_{1}$ and $\dot{x}\rightarrow 0$. This
implies that the kinetic energy of the corresponding conservative system would
tend to $0$ and the potential a constant $C^{2}A_{1}^{2}/2$ which is equal to
the energy loss of the original system. Both the mechanical energy of the
conservative system (32) at initial instance and $t\rightarrow\infty$ are
$c^{2}A_{1}^{2}/2$.

###  3.2 Two-dimensional Analytical Example

Let us consider a special two-dimensional mechanical system

$\displaystyle\ddot{x}+\dot{x}-\dot{y}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\ddot{y}-\dot{x}+\dot{y}$ $\displaystyle=$ $\displaystyle 0.$
(33)

The exact solution of the equation above with initial initial condition
$x_{0},y_{0},\dot{x}_{0},\dot{y}_{0}$ is

$\displaystyle x(t)$ $\displaystyle=$
$\displaystyle-\frac{\dot{y}_{0}-\dot{x}_{0}-4x_{0}}{4}+\frac{{e}^{-2t}(\dot{y}_{0}-\dot{x}_{0})}{4}+\frac{t(\dot{y}_{0})}{2}+\frac{t(\dot{x}_{0})}{2}$
$\displaystyle y(t)$ $\displaystyle=$
$\displaystyle\frac{\dot{y}_{0}-\dot{x}_{0}+4y_{0}}{4}-\frac{{e}^{-2t}(\dot{y}_{0}-\dot{x}_{0})}{4}+\frac{t(\dot{y}_{0})}{2}+\frac{t(\dot{x}_{0})}{2}$
(34)

For convenience to obtain $t=t(x),t=t(y)$，let $\dot{x}_{0}+\dot{y}_{0}=0$,
then simplify the particular solution above to

$\displaystyle
x(t)=-\frac{\dot{y}_{0}-\dot{x}_{0}-4x_{0}}{4}+\frac{{e}^{-2t}(\dot{y}_{0}-\dot{x}_{0})}{4}$
$\displaystyle
y(t)=\frac{\dot{y}_{0}-\dot{x}_{0}+4y_{0}}{4}-\frac{{e}^{-2t}(\dot{y}_{0}-\dot{x}_{0})}{4},$
(35)

From the equation above, we derived the velocity:

$\displaystyle\dot{x}$ $\displaystyle=$
$\displaystyle-\frac{{e}^{-2t}(\dot{y}_{0}-\dot{x}_{0})}{2},$ (36)
$\displaystyle\dot{y}$ $\displaystyle=$
$\displaystyle\frac{{e}^{-2t}(\dot{y}_{0}-\dot{x}_{0})}{2}$ (37)

Let the phase curve be denoted as $\gamma$. From Eq.(35), we obtain the
inverse functions

$\displaystyle t$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\ln\left[\frac{4}{\dot{y}_{0}-\dot{x}_{0}}(x-x_{0})+1\right]$
(38) $\displaystyle t$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\ln\left[-\frac{4}{\dot{y}_{0}-\dot{x}_{0}}(y-y_{0})+1\right]$
(39)

Substituting Eq.(38)(39) into Eq.(36), we have the map at $\gamma$ from $x,y$
to $\dot{x}$:

$\displaystyle\dot{x}(x)$ $\displaystyle=$
$\displaystyle-2x-\frac{\dot{y}_{0}-\dot{x}_{0}}{2}+2x_{0}$ (40)
$\displaystyle\dot{x}(y)$ $\displaystyle=$ $\displaystyle
2y-\frac{\dot{y}_{0}-\dot{x}_{0}}{2}-2y_{0}$ (41)

Substituting Eq.(38)(39) into Eq.(37), we have the map at $\gamma$ from $x,y$
to $\dot{y}$:

$\displaystyle\dot{y}(y)$ $\displaystyle=$
$\displaystyle-2y+\frac{\dot{y}_{0}-\dot{x}_{0}}{2}+2y_{0}$ (42)
$\displaystyle\dot{y}(x)$ $\displaystyle=$ $\displaystyle
2x+\frac{\dot{y}_{0}-\dot{x}_{0}}{2}-2x_{0}$ (43)

The components of nonconservative $\bm{F}$ in the system (33) are

$\displaystyle F_{1}$ $\displaystyle=$ $\displaystyle\dot{x}-\dot{y}$ (44)
$\displaystyle F_{2}$ $\displaystyle=$ $\displaystyle-\dot{x}+\dot{y}$ (45)

Substituting Eq.(40)(43) into $F_{1}$(44), then take the quantity as the first
component the conservative force $\mathcal{F}$:

$\mathcal{F}_{1}(x)=-4x-(\dot{y}_{0}-\dot{x}_{0})+4x_{0}.$ (46)

Substituting Eq.(41)(42) into $F_{2}$(45), then take the quantity as the
second component the conservative force $\mathcal{F}$:

$\mathcal{F}_{2}(y)=-4y+(\dot{y}_{0}-\dot{x}_{0})+4y_{0}$ (47)

Since $\partial\mathcal{F}_{1}/\partial y=\partial\mathcal{F}_{2}/\partial
x=0$, $\mathcal{F}$ must be conservative. Consequently we obtain a new
conservative system:

$\displaystyle\ddot{x}$ $\displaystyle=$ $\displaystyle-\mathcal{F}_{1}$
$\displaystyle=$ $\displaystyle 4x+(\dot{y}_{0}-\dot{x}_{0})-4x_{0}$
$\displaystyle\ddot{y}$ $\displaystyle=$ $\displaystyle-\mathcal{F}_{2}$ (48)
$\displaystyle=$ $\displaystyle 4y-(\dot{y}_{0}-\dot{x}_{0})-4y_{0}.$

We can readily prove that the particular solution (35) can satisfy Eq.(48)
too. In this case, this point agrees with Proposition 1 too.

###  3.3 A Formell Example in Vibration Mechanics

Take an $n$-dimensional oscillator with damping as an example, the governing
equation of which is as below:

$\ddot{\bm{q}}+\mathsfsl{C}\dot{\bm{q}}+\mathsfsl{K}\bm{q}=0,$ (49)

where $\bm{q}=\left[q_{1},\dots,q_{n}\right]^{T}$, superscript $T$ denotes a
matrix transpose,

$\mathsfsl{C}=\left[\begin{array}[]{ccc}C_{11}&\dots&C_{1n}\\\
\vdots&\ddots&\vdots\\\
C_{n1}&\dots&C_{nn}\end{array}\right],\mathsfsl{K}=\left[\begin{array}[]{ccc}K_{11}&\dots&K_{12}\\\
\vdots&\ddots&\vdots\\\ K_{21}&\dots&K_{22}\end{array}\right]$

, and $C_{ij}$ and $K_{ij}$ are constants.

It is complicated to solve Eq.(49). If Eq.(49) is higher dimensional, it is
almost impossible to solve Eq.(49) analytically. Therefore we assume that a
solution exists already.

$\bm{q}=\bm{q}(t)=\left[q_{1}(t),\dots,q_{n}(t)\right].$ (50)

Suppose a group of inverse functions

$t=t(q_{1}),\dots,t=t(q_{n}).$ (51)

As in Eq.(19) we can consider that the damping forces are equal to some
conservative force under an initial condition

$\begin{array}[]{ccc}c_{11}\dot{q}_{1}=\varrho_{11}(q_{1})&\dots&c_{1n}\dot{q}_{n}=\varrho_{1n}(q_{1})\\\
\vdots&\ddots&\vdots\\\
c_{n1}\dot{q}_{1}=\varrho_{21}(q_{n})&\dots&c_{nn}\dot{q}_{n}=\varrho_{nn}(q_{n}),\end{array}$
(52)

where $\varrho_{ij}(q_{i})$ is a function of $q_{i}$. For convenience, these
conservative forces can be defined as functions which are analogous to elastic
restoring forces:

$\begin{array}[]{ccc}\varrho_{11}(q_{1})=\kappa_{11}(q_{1})q_{1}&\dots&\varrho_{1n}(q_{1})=\kappa_{1n}(q_{1})q_{1}\\\
\vdots&\ddots&\vdots\\\
\varrho_{n1}(q_{1})=\kappa_{n1}(q_{n})q_{n}&\dots&\varrho_{nn}(q_{n})=\kappa_{nn}(q_{n})q_{n},\end{array}$
(53)

where $\kappa_{ij}(q_{i})$ is a function of $q_{i}$. An equivalent stiffness
matrix $\mathsfsl{\tilde{K}}$ is obtained, which is a diagonal matrix

$\mathsfsl{\tilde{K}}_{ii}=\sum_{l=1}^{n}\kappa_{il}(q_{l}).$ (54)

Consequently an $n$-dimensional conservative system is obtained

$\bm{\ddot{q}}+(\mathsfsl{K}+\mathsfsl{\tilde{K}})\bm{q}=0$ (55)

which shares a common phase curve with the $n$-dimensional damping system(49).
The Hamiltonian of Eqs.(55) is

$\hat{H}=\frac{1}{2}\bm{p}^{T}\bm{p}+\frac{1}{2}\bm{q}^{T}\mathsfsl{K}\bm{q}+\int_{\bm{0}}^{\bm{q}}(\tilde{\mathsfsl{K}}\bm{q})^{T}\mathrm{d}\bm{q},$
(56)

where $\bm{0}$ is a zero vector, $\bm{p}=\dot{\bm{q}}$. $\hat{H}$ in Eq.(56)
is the mechanical energy of the conservative system(55), because
$\int_{\bm{0}}^{\bm{q}}(\tilde{\mathsfsl{K}}\bm{q})^{T}\mathrm{d}\bm{q}$ is a
potential function such that $\hat{H}$ doest not depend on any path.

###  3.4 Discussion

Based on the above, we can outline the relationship between a dissipative
mechanical system and a group of conservative systems by means of Fig. 1. The
relationship can be stated from two perspectives:

Figure 1: A Dissipative Mechanical System and Conservative Systems

If one explains the relationship from a geometrical perspective, one can
obtain Proposition 1. In this paper the conservative systems (12) and (55) are
called the substituting systems. Although a substituting system shares a
common phase curve with the original system, under other initial conditions
the substituting system exhibits different phase curves. Therefore the phase
flow of the substituting system differs from that of the original system, it
follows that the substituting systems is not equal to the original system.
According to Louisville’s theorem (2.1), the phase flow of the original
dissipative system Eq.(11) certainly does not preserve its phase volume, but
the phase flow of the substituting conservative Eq.(12) does.

One also could explain the relationship from a mechanical perspective. It is
known that there are non-conservative forces in a nonconservative system. The
total energy of the nonconservative system consists of the work done by
nonconservative forces. Hence the function form of the total energy depends on
a phase curve i.e. under an initial condition. If one constrains the total
energy function to a phase curve $\gamma$, the total energy function can be
converted into a function of $q,p$. One take $\hat{H}$ consisting of this new
function and a constant as a Hamiltonian quantity, such that a Hamilton’s
system (i.e., a conservative system) is obtained. Under the initial condition
mentioned above, the solution curve of the conservative system is the same as
that of the original nonconservative system; under other initial conditions
the solution curve of the conservative is different from that of the original
nonconservative system. Since one defines the forces(19,52,53,54) in the new
system, the Hamiltonian quantity of the conservative can be thought of as the
mechanical energy of the new conservative system as Eq.(56).

One might doubt that the orbit of a dissipative dynamical system must be
asymptotic, can the asymptotic orbit coincide with one of a conservative
mechanical system. In some literatureSunyishui_book_e_2008 , a conservative
system defined a system with the behavior of the preservation of phase volume.
HasselblattHasselblatt_Katok2003  had explained the question: ’A key to
understanding this difference is given by a property that is not directly
observed by looking at individual orbits but by considering the evolution of
large sets of initial conditions simultaneously, the preservation of phase
volume.’ This point agrees with the second part of the proof of the
Proposition 1.

The Hamiltonians of the new conservative systems in general are not
analytically integrable, unless the original mechanical system is integrable.
The reason is that the work done by damping force depends on the phase curve.
If the system is integrable, then the phase curve can be explicitly written
out, the system has an analytical solution, and therefore the work done by
damping force can be explicitly integrated. Subsequently, the Hamiltonian
$\hat{H}$ can be explicitly expressed. Most systems do not have an analytical
solution. Despite this, the Hamilton quantity, coordinates and momentum must
satisfy Eq.(12) under a certain initial condition. Why had KleinKlein1928
written, ”Physicists can make use of these theories only very little, an
engineers nothing at all”? The answer: when one is seeking an analytical
solution to a classical mechanics problem by utilizing Hamiltonian formalism,
in fact one must inevitably convert the problem back to Newtonian formalism.
This means that an explicit form of Hamiltonian quantity is not necessary for
classical mechanics. What is important is the relationship between $q,p$ and
the Hamiltonian quantity embodied in the Hamilton’s Equation.

##  4 conclusions

We can conclude that a dissipative mechanical system has such properties: for
any nonconservative classical mechanical system and any initial condition,
there exists a conservative one, the two systems share one and only one common
phase curve; the Hamiltonian of the conservative system is the sum of the
total energy of the nonconservative system on the aforementioned phase curve
and a constant depending on the initial condition. We can further conclude,
that a dissipative problem can be reformulated as an infinite number of non-
dissipative problems, one corresponding to each phase curve of the dissipative
problem. One can avoid having to change the definition of the canonical
momentum in the Hamilton formalism, because under a certain initial condition
the motion of one of the group of conservative systems is the same as the
original dissipative system.

## References

  * (1) Arnold., V.I.: Mathematical Methods of classical Mechanics, second edition.  Springer-Verlag, Berlin (1978)
  * (2) Arnold., V.I.: Mathematical aspects of classical and celestial mechanics.  Springer-Verlag, Berlin (1997)
  * (3) Bateman, H.: On dissipative systems and related variational principles.  Phys. Rev. 38(4), 815–819 (1931).  DOI 10.1103/PhysRev.38.815
  * (4) Caldirola, P.: Forze non conservative della meccanica quantistica.  Nuovo Cim 18, 393–400 (1941)
  * (5) F.Klein: Entwickelung der Mathematik im 19 Jahrhundert.  Teubner (1928)
  * (6) Hasselblatt, B., Katok, A.: A FIRST COURSE IN DYNAMICS with a Panorama of Recent Developments.  AMBRIDGE UNIVERSITY PRESS (2003)
  * (7) Kanai, E.: On the Quantization of the Dissipative Systems.  Progress of Theoretical Physics 3, 440–442 (1948)
  * (8) Kochan, D.: Functional integral for non-lagrangian systems.  Phys. Rev. A 81(2), 022,112 (2010).  DOI 10.1103/PhysRevA.81.022112
  * (9) Kochan, D.: How to quantize forces (?): An academic essay on how the strings could enter classical mechanics.  Journal of Geometry and Physics 60(2), 219 – 229 (2010).  DOI DOI: 10.1016/j.geomphys.2009.09.014.  URL http://www.sciencedirect.com/science/article/B6TJ8-4XDCHN8-1/2/ed88d72d5e7e026557c477b9416a744e
  * (10) Krechetnikov, R., Marsden, J.E.: Dissipation-induced instabilities in finite dimensions.  Reviews of Modern Physics 79, 519–553 (2007).  DOI 10.1103/RevModPhys.79.519
  * (11) P.Morse, Feshbach, H.: Methods of Theoretical Physics.  McGraw-Hill, New York (1953)
  * (12) Rajeev, S.: A canonical formulation of dissipative mechanics using complex-valuedhamiltonians.  ANNALS of PHYSICS 322(3), 1541–1555 (2007)
  * (13) Sun, Y., Zhou, Y.: Introduction to Modern Celestial Mechanics.  Higher Education Press (2008)