collected by the agent during training, and we also change the loss function
from the Huber loss (Eq. 4.59) to the mean squared error so as to allow for
divergence in gradient. In our experiments, we find that whether the DQN
algorithm diverges or not is generally task-dependent, and it has a larger
probability to diverge if the task is more difficult. The result for Atari
2600 game Space Invaders is shown in Fig. 4.9. It can be seen that while the
DQN algorithm diverges, the C-DQN algorithm learns stably and its learning
speed is only slightly reduced. This confirms that the C-DQN algorithm is
convergent regardless of the properties of the training data.

Figure 4.9: Training performance and training loss on Atari 2600 games Space
Invaders when half of the data are randomly discarded.

Figure 4.10: Training performance and training loss on Atari 2600 games Space
Invaders when the replay memory adopts a random replacement strategy (left)
and when the size of the replay memory is reduced by a factor of 10 and adopts
different strategies (middle and right).

The same situation arises if when the replay memory (i.e. the dataset) is
full, one does not use the first-in-first-out (FIFO) strategy to replace old
data by new data, but randomly chooses old transition data to be replaced by
new data. In this case, the dataset can also contain incomplete trajectories
of data. In Fig. 4.10, we show that the DQN algorithm can actually diverge in
this simple setting. In existing literature on DQN algorithms, the replacement
strategy used with the dataset is often ignored, while here we find that it
can be an important detail that affects the final results of reinforcement
learning. In practice, the FIFO strategy is almost always used; nevertheless,
the FIFO strategy makes the data in the dataset less diverse, and it also
increases the risk of non-convergence due to the oscillation of the co-
evolvement of the learned policy and the dataset. As a consequence, a large
size of the replay memory is often necessary in reinforcement learning. In
Fig. 4.10, it can be seen that when we reduce the size of the replay memory by
a factor of $10$, the C-DQN algorithm can utilize the random replacement
strategy to achieve a higher performance, while the DQN algorithm cannot. Note
that the DQN algorithm does not show divergence in this experiment. We
conjecture that it is because the replay memory is small, it contains more
recent data and more complete trajectories of experience, which alleviates
divergence. As the C-DQN algorithm can be trained stably with data that comes
from an arbitrary distribution, the result opens up a new possibility of
reinforcement learning of only storing and learning important data, which is
not possible with the conventional DQN algorithm.

##### Difficult Games in Atari 2600

Here we consider difficult games in the Atari 2600 benchmark, which require
the use of large discount factors $\gamma$. Although the DQN algorithm becomes
unstable and often diverges when $\gamma$ becomes increasingly close to $1$,
the convergence property of the C-DQN algorithm does not depend on $\gamma$
and it can work with any $\gamma$ in principle. Nevertheless, we notice that a
large $\gamma$ does not necessarily lead to better performance, because a
large $\gamma$ requires the agent to learn to predict rewards that are distant
in future time steps, which is often not necessary and is irrelevant for
learning the task. Therefore, a large $\gamma$ can reduce the learning
efficiency. We also notice that if we have $\gamma\geq 0.9999$, the order of
magnitude of the term $(1-\gamma)Q_{\theta}$ becomes close to the inherent
noise in the gradient descent optimization algorithm due to the finite
learning rate, and the learning process can stagnate. Therefore, we do not
expect $\gamma$ to be larger than $0.9999$ and we only require $\gamma$ to
satisfy $0.99\leq\gamma\leq 0.9998$. Because the appropriate discount factor
$\gamma$ can be different for different problems, we also use a heuristic
algorithm to evaluate the frequency of reward signals in each problem so as to
determine $\gamma$ for each problem separately. Details concerning this
strategy are presented in Ref. [51] in the appendix. We also normalize the Q
functions before training using the evaluated mean and scale of the reward
signals, and we follow the techniques in Ref. [99] to transform the Q
functions approximately by the square root, so that the Q functions always
have appropriate orders of magnitude.

With the C-DQN algorithm and large values of $\gamma$, several difficult
problems which could not be solved by simple variants of the DQN algorithm can
now be solved, as shown in Fig. 4.11. Especially for Atari 2600 games Skiing,
Private Eye and Venture, the agent significantly benefits from large values of
$\gamma$ and achieves a higher best performance during training, despite the
fact that the games Private Eye and Venture are only partially observable
problems and therefore not fully learnable, which results in unstable
performance.

Figure 4.11: Training performance for several difficult games in the Atari
2600 benchmark. Each line shows the performance in a single repetition of the
experiment and the shaded regions show the standard deviation. The discount
factors $\gamma$ are shown in the titles. The DQN algorithm fails to learn
these tasks and shows significant instability, and the DQN loss increases up
to around $10^{5}\sim 10^{8}$ in these experiments, while the C-DQN loss stays
below $1$.

In the following, we compare the test performance of the C-DQN algorithm with
other works. To evaluate the test performance, we pick the best-performing
agent during training and test its performance using 400 trials, using the
$\epsilon$-greedy policy with $\epsilon=0.01$ and no-op starts666No-op starts
mean that whenever an episode begins, the no-operation action is executed
randomly for $1$ to $30$ frames, so that the agent does not always starts at
exactly the same state. [62]. The average of the test performances of the 3
different repetitions of our experiments are shown in Table 4.1 with the
standard error, compared with existing works and the human performance.

Table 4.1: Test performance on difficult Atari 2600 games. The results for the DQN algorithm are obtained by us using the same experimental settings as the C-DQN algorithm. Human results and results for Agent57 are due to Ref. [88], and results for Rainbow DQN are due to Ref. [69]. Note that the human results only correspond to reasonably adequate performance and they are not the highest possible performance of human. Task | C-DQN | DQN | Human | Rainbow DQN | Agent57 (SOTA)
---|---|---|---|---|---
Skiing | -3697 $\pm$ 157 | -29751 $\pm$ 224 | -4337 | -12958 | -4203 $\pm$ 608
Tennis | 10.9 $\pm$ 6.3 | -2.6 $\pm$ 1.4 | -8.3 | 0.0 | 23.8 $\pm$ 0.1
Private Eye | 14730 $\pm$ 37 | 7948 $\pm$ 749 | 69571 | 4234 | 79716 $\pm$ 29545
Venture | 893 $\pm$ 51 | 386 $\pm$ 85 | 1188 | 5.5 | 2624 $\pm$ 442

As we have followed the standard procedure of training the DQN agent on the
Atari 2600 benchmark as in Ref. [62] and Ref. [69], the performance obtained
by our C-DQN agent allows for a fair comparison with the results of the
Rainbow DQN algorithm in Ref. [69].777Precisely speaking, a fair comparison
with the Rainbow DQN algorithm cannot be made on the game Skiing. This is
because the reward clipping strategy adopted by the Rainbow DQN algorithm does
not permit learning the game Skiing. However, this does not affect our
conclusion. In Table 4.1, we see that the Rainbow DQN algorithm fails to make
progress in learning for these four difficult Atari 2600 games, and the C-DQN
algorithm can achieve performance higher than the Rainbow DQN algorithm and it
has non-trivial learning behaviour. The results of the Agent57 algorithm are
only for reference [88], representing the currently known highest overall
performance on the Atari 2600 benchmark. The results of the Agent57 algorithm
do not allow for a fair comparison with the C-DQN algorithm, because it
requires considerably more computation, sophisticated methods, and larger and
more complicated neural networks. Notably, we find that our result on the game
Skiing is exceptional, where the C-DQN algorithm achieves the state-of-the-art
performance despite its simplicity, which is discussed in following.

##### The Atari Game Skiing

77footnotetext: https://github.com/mgbellemare/Arcade-Learning-Environment

In Table 4.1, an exceptional result is that the C-DQN algorithm achieves a
performance that is higher than the Agent57 algorithm on the game Skiing, in
fact, utilizing an amount of computation budget less than $0.1\%$ of that of
the Agent57 algorithm. We find that this is the highest performance so far and
therefore achieves the state-of-the-art (SOTA) for this task. To elucidate the
reason, we describe the game in the following.

Figure 4.12: A screenshot of game Skiing in the Atari 2600 benchmark. The
program is provided under the GNU General Public License v2.0 in Ref.
[50].99footnotemark: 9

Figure 4.13: Training performance of the C-DQN algorithm on Atari 2600 game
Skiing when the learning rate is reduced by half, following the same
experimental procedure as in Fig. 4.11. The standard deviation is shown as the
shaded region.

A screenshot of Atari 2600 game Skiing is shown in Fig. 4.12. This game is
basically a racing game, where the player needs to go downhill as fast as
possible, and the time elapsed before reaching the goal is regarded as the
minus reward. For each time step, the agent receives a small amount of minus
reward which represents the time elapsed, until it reaches the goal and the
game terminates. Additionally, the player is required to pass through the
gates along the way, which are shown by the two small blue flags in Fig. 4.12.
If the player fails to pass through a gate, at the moment when the player
reaches the goal, a 5-second penalty is added to the elapsed time. The number
of gates that are yet to pass are shown on the top of the screen in the game.
With the standard setting in Ref. [62], the number of time steps for an
episode in this game is $\sim 1300$ for the random policy, is $\sim 4500$ if
the player slows down, and is $\sim 500$ if the policy is near-optimal.

Because the penalty for not passing through a gate is given to the agent only
at the end of the game, it is necessary for the agent to learn the relation
between the penalty at the end of the game and the events that occur early in
the game. Therefore, the time horizon for planning must be long enough, and
the discount factor $\gamma$ should be at least $1-\frac{1}{500}$ to allow for
learning. However, learning may stagnate if $\gamma$ is not even larger than
$1-\frac{1}{500}$, because if $\gamma$ is small, the agent would prefer
spending longer time before reaching the goal, so that the penalty at the end
of the game is delayed and the Q values for the states in the early game are
increased, which will further increase the number of time steps for an episode
and make learning difficult. Therefore, we have tuned our hyperparameter
setting so that we have $\gamma\approx 1-\frac{1}{5000}$ on this game. The
C-DQN agent learns with this value of $\gamma$ successfully and produces a new
record on this task. The large fluctuations on the learning curves shown in
Fig. 4.11 are mainly due to noise coming from the large learning rate, which
have been confirmed in Fig. 4.13 by repeating the experiment using a smaller
learning rate. However, we find that with the small learning rate, the policy
can easily get trapped in local optima and the test performance is actually
worse, and therefore we still use the large learning rate in our experiments.
It is worth noting that in fact, we cannot compare this result fairly with
that of the Agent57 algorithm, because our hyperparameters have been tuned so
that the discount factor $\gamma$ is suitable for this task, while the Agent57
algorithm adopts a bandit algorithm to dynamically determine $\gamma$, which
is more general.

###  4.5 Conclusions and Outlook

In this chapter, we have discussed the inefficiency issues of the RG
algorithm, and proposed a convergent DQN (C-DQN) algorithm to address the
long-standing problem of non-convergence of Q-learning. We have discussed the
properties of the C-DQN algorithm and demonstrated the effectiveness of the
C-DQN algorithm on the standard Atari 2600 benchmark for deep reinforcement
learning. With the stability of the C-DQN algorithm, we can tune the discount
factor $\gamma$ freely without sacrificing stability, and we can consider the
possibility of only learning important pieces of data to improve efficiency.
The C-DQN algorithm has a better stability and convergence property than the
DQN algorithm, and it can be applied to difficult tasks for which the DQN
algorithm fails due to instability. The idea of the C-DQN algorithm can be
combined with other reinforcement learning strategies involving target
networks and potentially improve their stability.

Many outstanding issues exist concerning the C-DQN algorithm. The C-DQN loss
is non-smooth, and it is still not clear how the non-smoothness affects the
optimization and the learning process. When the task is stochastic, the MSBE
loss $L_{\textit{MSBE}}$ used in the C-DQN algorithm does not converge exactly
to the optimal Q function, and therefore, it would be desirable if the C-DQN
algorithm can be improved so that stochastic tasks can be learned correctly
without bias. It would be interesting to investigate how the C-DQN loss
interplays with several other DQN extensions such as distributional DQN and
soft Q-learning [91, 92] and how the gradient descent optimization algorithm
affects the learning dynamics of the C-DQN algorithm. It is also an
interesting problem how the target network in the C-DQN algorithm can be
updated smoothly as in the DDPG algorithm [100].

##  Chapter 5 Control of Continuous Quantum Systems with Convergent Deep
Q-Learning

###  5.1 Introduction

In this chapter, we apply our convergent deep Q network (C-DQN) algorithm [51]
to quantum measurement-feedback control problems in continuous space, namely,
measurement-feedback cooling of a quantum-mechanical quartic oscillator and
measurement-feedback cooling of a trapped quantum-mechanical rigid body in
numerical simulation, and we compare the results obtained using C-DQN with
those obtained using the conventional DQN algorithm to demonstrate the
advantages of C-DQN. In Section 5.2, we present our model of the controlled
quantum quartic oscillator, and show that the C-DQN algorithm performs
significantly more stably compared with the conventional DQN algorithm, which
suffers much less randomness in the final results. In Section 5.3, we first
introduce the background of the quantum-mechanical rigid body and review the
derivation of the Hamiltonian following Ref. [101] in Section 5.3.1, and then,
we present our model of the trapped and controlled quantum-mechanical rigid
body and derive its Hamiltonian and the time-evolution equation in Section
5.3.2. We then apply the C-DQN and the DQN algorithms to the cooling problem
of this system, and compare them with the standard linear-quadratic-Gaussian
control strategy that involves approximation in Section 5.3.3.

###  5.2 Cooling of a Quantum Quartic Oscillator

In this section, we consider the problem of measurement-feedback cooling of a
one-dimensional quantum quartic oscillator, which is a minimal model of a
nonlinear quantum system in continuous space. We consider a situation in which
the system is subject to continuous position measurement, and we consider a
controllable external force for use in feedback control to reduce the energy
of the system, thereby cooling the system. We show that while the conventional
DQN algorithm exhibits instability and a large variance in learning the task,
the C-DQN algorithm is much more stable. We have mostly followed our previous
work [102] concerning the settings of the quantum system.

####  5.2.1 Significance of the Quartic System

In contrast to a quantum quartic oscillator, a quantum harmonic oscillator is
simple and its optimal control strategy can be analytically obtained [22]. The
Hamiltonian of a harmonic oscillator is given by

$\hat{H}=\frac{\hat{p}^{2}}{2m}+\frac{k}{2}\hat{x}^{2},$ (5.1)

where $\hat{p}$ is the momentum operator and $\hat{x}$ is the position
operator. The Hamiltonian is therefore quadratic with respect to the operators
$\hat{x}$ and $\hat{p}$, and the time-evolution equations of the expectation
values $\langle\hat{x}\rangle$ and $\langle\hat{x}\rangle$ are linear and
given by

$d\langle\hat{x}\rangle=\frac{1}{m}\langle\hat{p}\rangle\,dt,\qquad
d\langle\hat{p}\rangle=k\langle\hat{x}\rangle\,dt.$ (5.2)

When the position of the system is continuously measured as discussed in
Section 2.2.3, the time evolution of the position $\langle\hat{x}\rangle$ and
the momentum $\langle\hat{p}\rangle$ is subject to a Gaussian noise, and the
standard linear-quadratic-Gaussian (LQG) control is applicable. The LQG
control regards $\langle\hat{x}\rangle$ and $\langle\hat{p}\rangle$ as the
system variables and it effectively minimizes the term

$\mathbb{E}\left[\int\left(\frac{1}{2m}\langle\hat{p}\rangle^{2}+\frac{k}{2}\langle\hat{x}\rangle^{2}\right)dt\right],$
(5.3)

which corresponds to the minimization of the expectation of the total energy

$\mathbb{E}\left[\int\langle\hat{H}\rangle\,dt\right]=\mathbb{E}\left[\int\left(\frac{1}{2m}\langle\hat{p}\rangle^{2}+\frac{k}{2}\langle\hat{x}\rangle^{2}+\frac{1}{2m}(\langle\hat{p}^{2}\rangle-\langle\hat{p}\rangle^{2})+\frac{k}{2}(\langle\hat{x}^{2}\rangle-\langle\hat{x}\rangle^{2})\right)\,dt\right],$
(5.4)

since the variances $\langle\hat{p}^{2}\rangle-\langle\hat{p}\rangle^{2}$ and
$\langle\hat{x}^{2}\rangle-\langle\hat{x}\rangle^{2}$ in the above equation
are known to converge to steady values under continuous measurement [22], and
the state becomes a Gaussian state. Therefore, the quantum harmonic oscillator
under continuous measurement is effectively classical, in the sense that the
position and momentum variables $\langle\hat{x}\rangle$ and
$\langle\hat{p}\rangle$ are sufficient to describe the state in the time
evolution, and the control strategy can be conveniently derived in the same
way as for classical systems. The control force of the LQG control is linear
with respect to the variables $\langle\hat{x}\rangle$ and
$\langle\hat{p}\rangle$.

Nevertheless, for a quartic oscillator of which the Hamiltonian is given by

$\hat{H}=\frac{\hat{p}^{2}}{2m}+\lambda\hat{x}^{4},$ (5.5)

the time evolution of $\langle\hat{x}\rangle$ and $\langle\hat{p}\rangle$ is
given by

$d\langle\hat{x}\rangle=\frac{\langle\hat{p}\rangle}{m},\qquad
d\langle\hat{p}\rangle=-4\lambda\langle\hat{x}^{3}\rangle,$ (5.6)

which involves the cubic term $\langle\hat{x}^{3}\rangle$. Therefore, the
skewness
$\left\langle\left(\hat{x}-\langle\hat{x}\rangle\right)^{3}\right\rangle$ is
also relevant in the dynamics. The time evolution of the skewness is given by
[102]

$d\left\langle\left(\hat{x}-\langle\hat{x}\rangle\right)^{3}\right\rangle=\frac{3\left\langle(\hat{x}-\langle\hat{x}\rangle)(\hat{p}-\langle\hat{p}\rangle)(\hat{x}-\langle\hat{x}\rangle)\right\rangle}{m}dt$
(5.7)

and the time evolution of
$\left\langle(\hat{x}-\langle\hat{x}\rangle)(\hat{p}-\langle\hat{p}\rangle)(\hat{x}-\langle\hat{x}\rangle)\right\rangle$
is given by

$\begin{split}d\left\langle(\hat{x}-\langle\hat{x}\rangle)(\hat{p}-\langle\hat{p}\rangle)(\hat{x}-\langle\hat{x}\rangle)\right\rangle&=\frac{2\left\langle(\hat{p}-\langle\hat{p}\rangle)(\hat{x}-\langle\hat{x}\rangle)(\hat{p}-\langle\hat{p}\rangle)\right\rangle}{m}dt\\\
&\quad-4\lambda\left\langle(\hat{x}^{3}-\langle\hat{x}^{3}\rangle)(\hat{x}-\langle\hat{x}\rangle)^{2}\right\rangle
dt.\end{split}$ (5.8)

For a Gaussian state, of which the skewness, the odd central moments and the
excess kurtosis are all zero, the term
$\left\langle(\hat{x}^{3}-\langle\hat{x}^{3}\rangle)(\hat{x}-\langle\hat{x}\rangle)^{2}\right\rangle$
is found to be [102]

$\begin{split}\left\langle(\hat{x}^{3}-\langle\hat{x}^{3}\rangle)(\hat{x}-\langle\hat{x}\rangle)^{2}\right\rangle=9\left(\langle\hat{x}^{2}\rangle-\langle\hat{x}\rangle^{2}\right)\langle\hat{x}\rangle,\end{split}$
(5.9)

which implies that a state that is initially Gaussian in a quartic potential
gradually develops nonzero skewness. The dynamics of $\langle\hat{x}\rangle$
and that of $\langle\hat{p}\rangle$ are affected by the profile of the wave
function, and therefore, the state cannot be merely characterized by the
expectation values $\langle\hat{x}\rangle$ and $\langle\hat{p}\rangle$ and the
variances, and the dynamics exhibits genuinely quantum-mechanical effects, as
shown in Figs. 5.1 and 5.2.

Figure 5.1: Snapshots of the time evolution of a state in a quartic potential.
The left panels shows the initial state which is Gaussian, and after several
oscillations in the quartic potential, the state becomes non-Gaussian and it
is shown in the right panel. The blue and the orange curves show the real and
the imaginary parts of the wave functions, and the red curves show the
probability densities. The grey curves show the potential, and the scale of
the potential is arbitrary.  Figure 5.2: The time evolution of the expectation
value of the position $\langle x\rangle$ of the state evolving in the quartic
potential as in Fig. 5.1. It can be seen that the expectation value of the
position ceases to oscillate, while the energy of the system remains high.

It has been known that a one-dimensional quartic oscillator corresponds to the
one-dimensional $\phi^{4}$ theory [103], and the system is difficult to
analyse. The appropriate control strategy of this system is unknown, and
therefore, we consider cooling of the quartic oscillator as our first
nontrivial example of continuous quantum control, and we compare the result
obtained using the conventional DQN algorithm with the result obtained using
the C-DQN algorithm.

####  5.2.2 Model

Following Section 2.2.3, the stochastic time-evolution equation of the state
subject to continuous position measurement is given by

$\displaystyle
d|\psi\rangle=\left[\left(-\frac{i}{\hbar}\hat{H}-\frac{\gamma}{4}(\hat{x}-\langle\hat{x}\rangle)^{2}\right)dt+\sqrt{\dfrac{\gamma}{2}}(\hat{x}-\langle\hat{x}\rangle)dW\right]|\psi\rangle,$
(5.10)
$\displaystyle\hat{H}=\frac{\hat{p}^{2}}{2m}+\lambda\hat{x}^{4}-F_{\text{con}}\hat{x},$
(5.11)

where $\gamma$ is the measurement strength and $dW$ is a Wiener increment,
which is a Gaussian random variable satisfying $\mathbb{E}[dW]=0$ and
$\mathbb{E}[dW^{2}]=dt$, and $F_{\text{con}}$ is the external control force.
Due to the Gaussian property of the position measurement, the system tends to
behave like a classical particle and to have a Gaussian profile when the
measurement strength is large. Therefore, we choose a sufficiently small
measurement strength so that quantum effects can be significant. The system
parameters that we use in our numerical simulation are given in Table 5.1.

|  $m$ ($m_{c}$) |  $\lambda$ $\left(\frac{m^{2}_{c}\omega^{3}_{c}}{\hbar}\right)$ |  $\gamma$ ($\frac{m_{c}\omega_{c}^{2}}{\hbar}$) |  $F_{\text{max}}$ ($\sqrt{\hbar m_{c}\omega_{c}^{3}}$) |  $x_{\text{max}}$ ($\sqrt{\frac{\hbar}{m_{c}\omega_{c}}}$)
---|---|---|---|---|---
quartic oscillator | $\dfrac{1}{\pi}$ | $\dfrac{\pi}{25}$ | $\dfrac{\pi}{100}$ | $3\pi$ | 8.5
Table 5.1: System parameters of the quartic oscillator used in our numerical
experiments, in terms of a reference angular momentum $\omega_{c}$ and a
reference mass $m_{c}$. $F_{\text{max}}$ is the maximum of the control force
$F_{\text{con}}$ that we allow. $x_{\text{max}}$ is the boundary of the 1D
space that we simulate, or the maximal distance away from the center of the
potential, in our numerical simulation.

We assume that the measurement efficiency is unity. As the measurement can
purify an arbitrary mixed state and the state is continuously measured, we
assume that the state, or the wave function, is already known by the external
observer and is therefore available. Every $\frac{1}{18\omega_{c}}$ time, the
controller determines the control force $F_{\text{con}}$ on the basis of the
information about the instantaneous wave function $|\psi\rangle$, and the
force $F_{\text{con}}$ is kept constant during a time of
$\frac{1}{18\omega_{c}}$. The control loss, or the minus reward in the setting
of reinforcement learning, is the energy of the state, given by
$\left\langle\frac{\hat{p}^{2}}{2m}+\lambda\hat{x}^{4}\right\rangle$.

To meet the requirement of deep Q-learning, the space of actions, namely, the
space of possible choices of the control force $F_{\text{con}}$, must be
discrete. We therefore discretize the continuous interval
$\left[-F_{\text{max}},+F_{\text{max}}\right]$ into 21 points, and the set of
control actions is given by

$\\{F_{\text{con}}\ |\ F_{\text{con}}=n\times 0.3\pi\sqrt{\hbar
m_{c}\omega_{c}^{3}},\quad-10\leq n\leq 10,\quad n\in\mathbb{Z}\\}.$ (5.12)

In the setting of reinforcement learning, the agent, or the controller,
determines its action $F_{\text{con}}$ based on the current state
$|\psi\rangle$. After an evolution time of $\frac{1}{18\omega_{c}}$ of the
state, the agent experiences its next time step and it observes the state
again to make decision on the control action. The reward is given by the minus
energy of the state at each time step.

The deep neural network $Q_{\theta}$ used in reinforcement learning as
discussed in Chapter 3 takes information of the state $|\psi\rangle$ as its
input, and it outputs the Q values for each $F_{\text{con}}$ in the action set
in Eq. (5.12). The neural network has 4 layers in total, with hidden units
512, 512, and 256. The neural network is trained using the Adam optimizer [64]
with a minibatch size 512, learning each experience data for 8 times on
average, and the update period of the target network is set to be 300 gradient
descent steps. The discount factor $\gamma$ in Q-learning discussed in Chapter
3 is set to be $0.99$.

To efficiently give the necessary information of the state to the neural
network, we use distribution moments of the Wigner quasiprobability
distribution of the state as the input of the neural network. Specifically, we
include $\langle\hat{x}\rangle$, $\langle\hat{p}\rangle$,
$\left\langle\left(\hat{x}-\langle\hat{x}\rangle\right)^{2}\right\rangle$,
$\left\langle\left(\hat{p}-\langle\hat{p}\rangle\right)^{2}\right\rangle$,
$\text{Re}\left[\left\langle\left(\hat{x}-\langle\hat{x}\rangle\right)\left(\hat{p}-\langle\hat{p}\rangle\right)\right\rangle\right]$,
$\left\langle\left(\hat{x}-\langle\hat{x}\rangle\right)^{3}\right\rangle$, and
$\left\langle\left(\hat{x}-\langle\hat{x}\rangle\right)\left(\hat{p}-\langle\hat{p}\rangle\right)\left(\hat{x}-\langle\hat{x}\rangle\right)\right\rangle$,
etc. in the input of the neural network, totally up to all fifth central
distribution moments with respect to $\hat{x}$ and $\hat{p}$. The motivation
for using distribution moments as the input is that the LQG control only needs
$\langle\hat{x}\rangle$ and $\langle\hat{p}\rangle$ to find the optimal
control force for a harmonic oscillator, and naturally by including higher
distribution moments, the controller may become aware of additional nontrivial
details of the state and better control strategies may be found. These
distribution moments are also physically relevant, as they are physical
observables.

The quantum state is approximately simulated in discrete space, and the
spacing between adjacent discrete sites is set to be
$0.1\sqrt{\frac{\hbar}{m_{c}\omega_{c}}}$, and the time step used in the
numerical integration of the time-evolution equation is set to be
$\frac{1}{1440\omega_{c}}$. To numerically evaluate the term
$\dfrac{\partial}{\partial x}$ by finite difference methods, we use the
formula

$\begin{split}\frac{\partial}{\partial
x}f(x_{0})=&\frac{3f(x_{0}-4h)-32f(x_{0}-3h)+168f(x_{0}-2h)-672f(x_{0}-h)}{840h}\\\
&+\frac{672f(x_{0}+h)-168(x_{0}+2h)+32(x_{0}+3h)-3(x_{0}+4h)}{840h}+O(h^{8})\end{split}$
(5.13)

and to evaluate the term $\dfrac{\partial^{2}}{\partial x^{2}}$, we use

$\begin{split}\frac{\partial^{2}}{\partial
x^{2}}f(x_{0})=&\frac{-9f(x_{0}-4h)+128f(x_{0}-3h)-1008f(x_{0}-2h)+8064f(x_{0}-h)-14350f(x_{0})}{5040h}\\\
&+\frac{8064f(x_{0}+h)-1008(x_{0}+2h)+128(x_{0}+3h)-9(x_{0}+4h)}{5040h}+O(h^{8}),\end{split}$
(5.14)

so that we can obtain accurate results efficiently in our numerical
simulation. We use the implicit 1.5 order strong scheme [104] for numerical
integration of the stochastic differential equation (5.10), and additionally,
in order to prevent numerical divergence due to the high energy part of the
quartic potential, we include high-order terms of the time evolution of the
Hamiltonian in our numerical integration, including additional terms
$\frac{(-\frac{i}{\hbar}\,dt\,\hat{H})^{n}}{n!}$ up to $n=6$. The programming
codes of our numerical experiments have been publicized for
reference,111https://github.com/Z-T-
WANG/PhDThesis/tree/main/quartic%20oscillator where all details can be found.
For discussion concerning numerical integration of stochastic differential
equations, see appendix in Ref. [102].

####  5.2.3 Evaluation and Results

After establishing the model, we train the reinforcement learning agent to
reduce the energy of the state using the control force $F_{\text{con}}$
through trial and error. In the following, we describe how we have set the
task for the agent and how we evaluate the performance of the agent, and the
results are presented, comparing the DQN algorithm with the C-DQN algorithm.

##### Initialization

To make sure that the initial state at the beginning of the control is a
typical non-Gaussian state in the quartic potential, at the initialization, we
first initialize the state as a Gaussian wave packet at the center of the
potential with a random momentum ranging from $-0.3\pi\sqrt{\hbar
m_{c}\omega_{c}}$ to $+0.3\pi\sqrt{\hbar m_{c}\omega_{c}}$, and then, we let
the state evolve during a time which changes randomly from
$\dfrac{15}{\omega_{c}}$ to $\dfrac{20}{\omega_{c}}$ with $F_{\text{con}}=0$,
and lastly we use the resulting state as the initial state for the
reinforcement learning agent to start to control. This procedure ensures that
the initial state from the perspective of the controller is sufficiently non-
Gaussian, and the initial energy of the state approximately lies in
$5\hbar\omega_{c}$ and $7\hbar\omega_{c}$. We re-initialize the state if its
energy exceeds $7\hbar\omega_{c}$.

##### Setting of Episodes

Following convention in reinforcement learning, we set an episode of the
control to be $\dfrac{100}{\omega_{c}}$ time in simulation. In other words, we
maximally allow the state to evolve for a time interval of
$\dfrac{100}{\omega_{c}}$ under control, and after that evolution, we stop the
simulation and re-initialize the state to start over again. To evaluate the
performance of the controller, we calculate the average energy of the
controlled state starting from time $\frac{30}{\omega_{c}}$ to
$\frac{100}{\omega_{c}}$, so that transient behaviour of the control when the
state is initialized is ignored.

In order to make sure that the state does not move out of the boundary of the
space that is simulated, we stop the simulation and end the episode whenever
the probability distribution around the boundary of the space is not
negligible. Also, regarding the energy, we end the episode whenever the energy
of the state exceeds $12\hbar\omega_{c}$. In other words, $12\hbar\omega_{c}$
is the maximal energy that we allow during the control, beyond which we stop
the control and regard the control as having failed in this episode. In these
two cases of failure, the evaluated performance, i.e. the average energy, is
considered to be $12\hbar\omega_{c}$. If an episode ends with failure, we set
the Q value that is to be learned at the end of the episode to be the final
energy of the state divided by $\frac{1}{1-\gamma}$, which is a large value,
so that the reinforcement learning agent should learn to avoid getting close
to these failure cases.

At the beginning of training, the control of the agent fails quickly and
frequently, and as the learning makes progress, the state is kept at a low
energy for a longer period of time, and finally the state can be stabilized
and cooled. In training, we first train the agent until it can stabilize the
state for a time interval of $\frac{100}{\omega_{c}}$, i.e., until it
completes a whole episode without failure, and then, we train the agent for
11000 more episodes, during which we gradually reduce the learning rate and
the $\epsilon$ hyperparameter of the $\epsilon$-greedy policy in Q-learning.

##### Results

Figure 5.3: Learning curves of the quartic cooling problem for the DQN
algorithm and the C-DQN algorithm, with the same experimental settings, for 5
different random seeds. The abscissa shows the simulated time of the evolution
of the quartic oscillator, which represents the number of learned data and the
training time. The ordinate shows the average energy, which shows the
performance of the controller, and a smaller value represents a better
performance. Gaussian smoothing with the standard deviation of 40 is applied
to smooth the performance data. The failure of stabilization and control
corresponds to an energy of 12.

We repeat the experiment for 5 times using 5 different random seeds, and the
learning curves of the DQN algorithm and the C-DQN algorithm in the training
process are shown in Fig. 5.3, 5.4 and 5.5, for each of the random seeds
separately, and summarized in Fig. 5.6.

In Fig. 5.3, we see that although both the DQN and the C-DQN algorithm can
learn to reduce the energy of the system, the performance of C-DQN is
consistently better than DQN during training. In addition, whereas the
performance of C-DQN improves steadily and consistently during training, the
performance of DQN does not improve consistently, and the performance
fluctuates. Especially, the performance can get worse after the agent has
learned how to cool the system well. To see more details, we show the
performance data with less smoothing in Fig. 5.5. In Fig. 5.5, it is clear
that the performance of C-DQN is very stable, while the performance of DQN
strongly fluctuates in training when it does not perform well. Specifically,
if we check the failure rate of the controller in training, i.e. the
probability of an episode to end with control failure, we see that the
deterioration of the performance of DQN is strongly correlated with control
failure, as shown in Fig. 5.5, except for the case of the random seed 5.
Therefore, we see that the DQN algorithm does not perform stably throughout
the training process, and it tends to get unstable and occasionally fails even
at a later stage of training, and such instability is hard to remove and does
not completely disappear after an extended period of training.

Figure 5.4: Learning curves of the quartic cooling problem for the DQN
algorithm and the C-DQN algorithm, for 5 different random seeds, which shows
the same data as in Fig. 5.3 but using Gaussian smoothing with the standard
deviation of 4 to show the fluctuations in performance. Figure 5.5: Failure
rate of the controller learned by the DQN algorithm and the C-DQN algorithm in
the process of training, for 5 different random seeds, evaluated using a
Gaussian window with the standard deviation of 20.

Meanwhile, by comparing the results of the repetitions of the experiment with
different random seeds, we see that the performance of DQN has a large
variance with respect to the random seeds, while that of C-DQN almost has no
variance, as shown in Fig. 5.6. The learning curves of DQN can be
qualitatively different for different random seeds, and especially, the
results for the random seed 5 are qualitatively different compared with the
others, as can be seen in Fig. 5.3, 5.4 and 5.5. The agent does not cool the
system to an energy that is as low as the others, but it performs stably
throughout training and does not encounter many failures, which implies that
the agent has learns differently from the others. These results show that the
instability of DQN increases the randomness in its results, and therefore,
reduces the reproducibility and reliability of the results, especially if the
task is relatively difficult. The significant randomness in results
considerably increases the difficulty of improving and fine-tuning the AI
algorithm, and requires experimenters to repeat the experiments many times in
order to confirm a result, and the result may change qualitatively due to
minor details which are supposed to only slightly perturb the training [45].
Also, in practical scenarios, as one usually wants a controller that is
trained to have a performance as high as possible, one usually needs to run
the experiments many times in order to pick the best trained controller, and
therefore, the obtained performance eventually would highly depend on the
number of repetitions of the experiment.

Compared with the DQN algorithm, using exactly the same experimental setting,
the C-DQN algorithm does not encounter any difficulty encountered by the DQN
algorithm, and the C-DQN algorithm steadily approaches a satisfactory
performance without instability and with little randomness, as shown in Fig.
5.6. This clearly demonstrates the stability of C-DQN and the consistency of
its results, which translates into the reliability of the results of C-DQN in
solving scientific problems.

Figure 5.6: Learning curves of the quartic cooling problem for the DQN
algorithm (left) and the C-DQN algorithm (right). Each curve represent a
repetition of the experiment with a different random seed, summarizing the
results in Fig. 5.3.

As the satisfactory stability and consistency of C-DQN have been confirmed as
above for the nontrivial quantum control problem, we achieve our purpose of
developing a stable and reliable deep reinforcement learning algorithm for
physical control problems. In the following, we apply the C-DQN algorithm to a
more complicated system which is of realistic relevance.

###  5.3 Cooling of a Quantum-Mechanical Rigid Body

With rapid development in synthesis and manipulation of nanoparticles,
nanorotors have recently been experimentally realized and studied [47]. When
the nanoparticles are sufficiently isolated from their environment and
sufficiently cooled, quantum-mechanical effects of rotations are expected to
be observed, and these quantum-mechanical rotors are expected to find their
applications in sensing devices and fundamental tests of physical principles,
because they allow an unprecedented accuracy of sensing using their rotational
degrees of freedom. However, in experiments, the rotational degrees of freedom
have yet to be successfully cooled to approach a quantum regime so far, and
only the center-of-mass degree of freedom has been successfully cooled to
approach the ground state [105, 106]. The dynamics of rotation is essentially
nonlinear, and it is much harder to manipulate compared to the dynamics in
position space.

In the following, we consider the problem of measurement-feedback cooling of a
quantum-mechanical rigid body, which is both of realistic significance and of
theoretical interest. We derive the Hamiltonian of the system, describe our
model, and present the results of the control learned by the DQN and the C-DQN
algorithms, and compare the results with those obtained using the standard LQG
control which employs a linear approximation.

####  5.3.1 Derivation of the Rotational Hamiltonian

In this section, we derive the Hamiltonian of a free quantum rigid body in
terms of the angles which define the orientation of the rigid body, following
Ref. [101].

##### Euler Angles

A standard representation of the orientation of a rigid body is given by the
Euler angles $(\alpha,\beta,\gamma)$. There are several different conventions
on the Euler angles, and here we follow Ref. [101] and use the z-x-z
convention for the Euler rotation (see Eq. (5.15) below). In this case, the
rotation $(\alpha,\beta,\gamma)$ represents first rotating the rigid body
around the $z$ axis of the laboratory frame through angle $\gamma$, then
rotating around the $x$ axis of the laboratory frame through angle $\beta$,
and finally rotating around the $z$ axis of the laboratory frame through angle
$\alpha$. In terms of angular momentum operators
$\hat{L}_{x},\hat{L}_{y},\hat{L}_{z}$, this Euler rotation is given by

$e^{-\frac{i}{\hbar}\alpha\hat{L}_{z}}e^{-\frac{i}{\hbar}\beta\hat{L}_{x}}e^{-\frac{i}{\hbar}\gamma\hat{L}_{z}}.$
(5.15)

If we assume that before the rotation, the rigid body is at the default
position and the local $x$, $y$ and $z$ axes of the rigid body align with the
$x$, $y$ and $z$ axes of the laboratory frame, then, the rotation operator
$e^{-\frac{i}{\hbar}\alpha\hat{L}_{z}}e^{-\frac{i}{\hbar}\beta\hat{L}_{x}}e^{-\frac{i}{\hbar}\gamma\hat{L}_{z}}$
also represents the orientation of the rigid body in the laboratory frame
after the rotation. The angle $\beta$ is called the polar angle or the zenith
angle, the angle $\alpha$ is called the azimuthal angle, and the angle
$\gamma$ corresponds to self-rotation of the object around its local $z$ axis.

Although the representation of Euler angles is physically straightforward, it
suffers from the well-known problem called gimbal lock, or, coordinate
singularity, when $\beta$ is equal to $0$ and $\pi$. The cases of $\beta$
being equal to $0$ and $\pi$ correspond to the north pole and the south pole
of a spherical coordinate. When we have $\beta=0$, we have

$e^{-\frac{i}{\hbar}\alpha\hat{L}_{z}}e^{-\frac{i}{\hbar}\beta\hat{L}_{x}}e^{-\frac{i}{\hbar}\gamma\hat{L}_{z}}=e^{-\frac{i}{\hbar}(\alpha+\gamma)\hat{L}_{z}},$
(5.16)

and therefore $\alpha$ and $\gamma$ become equivalent and essentially one
degree of freedom disappears. A similar situation occurs for $\beta=\pi$, in
which case we have

$e^{-\frac{i}{\hbar}\alpha\hat{L}_{z}}e^{-\frac{i}{\hbar}\beta\hat{L}_{x}}e^{-\frac{i}{\hbar}\gamma\hat{L}_{z}}=e^{-\frac{i}{\hbar}(\alpha-\gamma)\hat{L}_{z}}e^{-\frac{i}{\hbar}\beta\hat{L}_{x}},$
(5.17)

and the remaining degree of freedom becomes $\alpha-\gamma$. This issue leads
to singularity when one attempts to find the wave function of a quantum state
in the space of Euler angles and attempts to differentiate with respect to
$\alpha$, $\beta$ and $\gamma$. Also, if one considers a smooth movement
across the point $\beta=0$, it is clear that the angles $\alpha$ and $\gamma$
flip and change discontinuously, which shows that a continuous physical
movement does not correspond to a continuous change in the coordinate
representation using Euler angles.

##### Quaternions

To overcome the problem of coordinate singularity, we use quaternions in the
following. The quaternion variables are related to the Euler angles by

$\begin{split}\xi=\cos\dfrac{\alpha-\gamma}{2}\sin\frac{\beta}{2},\\\
\eta=\sin\frac{\alpha-\gamma}{2}\sin\frac{\beta}{2},\\\
\zeta=\sin\frac{\alpha+\gamma}{2}\cos\frac{\beta}{2},\\\
\chi=\cos\frac{\alpha+\gamma}{2}\cos\frac{\beta}{2}.\end{split}$ (5.18)

It can be checked that at $\beta=0$, we have $\sin\frac{\beta}{2}=0$ and
$\xi=\eta=0$, and the remaining degree of freedom becomes $\alpha+\gamma$; at
$\beta=\pi$, we have $\cos\frac{\beta}{2}=0$ and $\zeta=\chi=0$, and the
remaining degree of freedom becomes $\alpha-\gamma$. All coordinate
singularities are removed by the quaternion representation. Nevertheless,
there is a redundant degree of freedom, and we have
$\xi^{2}+\eta^{2}+\zeta^{2}+\chi^{2}\equiv 1$, and the rotation of
$(\xi,\eta,\zeta,\chi)$ and that of $(-\xi,-\eta,-\zeta,-\chi)$ represent the
same physical rotation.

As an extension of complex numbers, in quaternion computation, the quaternions
are conventionally written as

$\chi+\xi\,i+\eta\,j+\zeta\,k,$ (5.19)

where $\chi$ represents a real number, and the calculation rules are given by

$i^{2}=j^{2}=k^{2}=-1,$ (5.20)

and

$i\times j=-j\times i=k,\quad j\times k=-k\times j=i,\quad k\times i=-i\times
k=j.$ (5.21)

The inverse of a quaternion is given by

$\frac{1}{\chi+\xi\,i+\eta\,j+\zeta\,k}=\frac{\chi-\xi\,i-\eta\,j-\zeta\,k}{\xi^{2}+\eta^{2}+\zeta^{2}+\chi^{2}}.$
(5.22)

The quaternions are intimately related to rotations in three-dimensional
space. For a vector $\vec{r}=(r_{x},r_{y},r_{z})$, after rotating around a
unit vector $\vec{a}=(a_{x},a_{y},a_{z})$ by angle $\theta$, the resulting
vector can be represented by

$\begin{split}&\left(\cos\frac{\theta}{2}+\sin\frac{\theta}{2}a_{x}i+\sin\frac{\theta}{2}a_{y}j+\sin\frac{\theta}{2}a_{z}k\right)\times\left(r_{x}i+r_{y}j+r_{z}k\right)\\\
&\times\left(\cos\frac{\theta}{2}-\sin\frac{\theta}{2}a_{x}i-\sin\frac{\theta}{2}a_{y}j-\sin\frac{\theta}{2}a_{z}k\right)\\\
={}&\left[\left(\cos\theta+a_{x}^{2}(1-\cos\theta)\right)r_{x}+\left(a_{x}a_{y}(1-\cos\theta)-a_{z}\sin\theta\right)r_{y}+\left(a_{z}a_{x}(1-\cos\theta)+a_{y}\sin\theta\right)r_{z}\right]i\\\
&+\left[\left(a_{x}a_{y}(1-\cos\theta)+a_{z}\sin\theta\right)r_{x}+\left(\cos\theta+a_{y}^{2}(1-\cos\theta)\right)r_{y}+\left(a_{y}a_{z}(1-\cos\theta)-a_{x}\sin\theta\right)r_{z}\right]j\\\
&+\left[\left(a_{z}a_{x}(1-\cos\theta)-a_{y}\sin\theta\right)r_{x}+\left(a_{y}a_{z}(1-\cos\theta)+a_{x}\sin\theta\right)r_{y}+\left(\cos\theta+a_{z}^{2}(1-\cos\theta)\right)r_{z}\right]k,\end{split}$
(5.23)

where the factors before $i$, $j$ and $k$ are the $x$, $y$ and $z$ coordinates
of the vector, respectively. Therefore, quaternions can be used to compute a
combination of rotations. The result of rotation 2,
$\chi_{2}+\xi_{2}\,i+\eta_{2}\,j+\zeta_{2}\,k$, after rotation 1,
$\chi_{1}+\xi_{1}\,i+\eta_{1}\,j+\zeta_{1}\,k$, can be computed as

$(\chi_{2}+\xi_{2}\,i+\eta_{2}\,j+\zeta_{2}\,k)\times(\chi_{1}+\xi_{1}\,i+\eta_{1}\,j+\zeta_{1}\,k),$
(5.24)

since the vector $\vec{r}$ after the rotation 1 and rotation 2 is given by

$\begin{split}&(\chi_{2}+\xi_{2}\,i+\eta_{2}\,j+\zeta_{2}\,k)\times(\chi_{1}+\xi_{1}\,i+\eta_{1}\,j+\zeta_{1}\,k)\times\left(r_{x}i+r_{y}j+r_{z}k\right)\\\
&\times(\chi_{1}-\xi_{1}\,i-\eta_{1}\,j-\zeta_{1}\,k)\times(\chi_{2}-\xi_{2}\,i-\eta_{2}\,j-\zeta_{2}\,k).\end{split}$
(5.25)

##### Angular Momentum Operators

Consider an infinitesimal rotation around the $z$ axis of the laboratory frame
by $\epsilon$, denoted by $\hat{R}_{z;\epsilon}$. In the quaternion
representation, the rotation is given by

$\hat{R}_{z;\epsilon}=1+\sin\frac{\epsilon}{2}\,k+O(\epsilon^{2}).$ (5.26)

Therefore, given an initial orientation
$\hat{R}_{0}=\chi_{0}+\xi_{0}i+\eta_{0}j+\zeta_{0}k$, after rotation
$\hat{R}_{z;\epsilon}$, the orientation becomes

$\begin{split}\hat{R}_{z;\epsilon}\hat{R}_{0}=(\chi_{0}-\frac{\epsilon}{2}\zeta_{0})+(\xi_{0}-\frac{\epsilon}{2}\eta_{0})i+(\eta_{0}+\frac{\epsilon}{2}\xi_{0})j+(\zeta_{0}+\frac{\epsilon}{2}\chi_{0})k.\end{split}$
(5.27)

That is, the rotation $\hat{R}_{z;\epsilon}$, or in terms of the angular
momentum operator, $e^{-\frac{i}{\hbar}\epsilon\hat{L}_{z}}$, shifts the
quaternion variables $(\xi,\eta,\zeta,\chi)$ by an amount of
$\frac{\epsilon}{2}$. Therefore, for functions defined in the space of the
quaternion variables $(\xi,\eta,\zeta,\chi)$, the angular momentum operator
$\hat{L}_{z}$ that generates the rotations around the $z$ axis is given by

$\hat{L}_{z}=-\frac{i\hbar}{2}\left(-\eta\frac{\partial}{\partial\xi}+\xi\frac{\partial}{\partial\eta}+\chi\frac{\partial}{\partial\zeta}-\zeta\frac{\partial}{\partial\chi}\right).$
(5.28)

Similarly, the other angular momentum operators are given by

$\hat{L}_{x}=-\frac{i\hbar}{2}\left(+\chi\frac{\partial}{\partial\xi}-\zeta\frac{\partial}{\partial\eta}+\eta\frac{\partial}{\partial\zeta}-\xi\frac{\partial}{\partial\chi}\right),$
(5.29)
$\hat{L}_{y}=-\frac{i\hbar}{2}\left(+\zeta\frac{\partial}{\partial\xi}+\chi\frac{\partial}{\partial\eta}-\xi\frac{\partial}{\partial\zeta}-\eta\frac{\partial}{\partial\chi}\right).$
(5.30)

To confirm the commutation relations among the angular momentum operators,
using the bilinearity of the commutator, we can easily obtain

$\begin{split}[\hat{L}_{x},\hat{L}_{y}]&=\hat{L}_{x}\hat{L}_{y}-\hat{L}_{y}\hat{L}_{x}\\\
&=-\frac{\hbar^{2}}{2}\left(-\chi\frac{\partial}{\partial\zeta}+\zeta\frac{\partial}{\partial\chi}+\eta\frac{\partial}{\partial\xi}-\xi\frac{\partial}{\partial\eta}\right)\\\
&=i\hbar\hat{L}_{z},\end{split}$ (5.31)

where $[\cdot,\cdot]$ is the commutator, defined as $[A,B]:=AB-BA$. Similarly
we have

$=i\hbar\hat{L}_{x},\qquad[\hat{L}_{z},\hat{L}_{x}]=i\hbar\hat{L}_{y},$ (5.32)

which are the standard commutation relations of the angular momentum
operators.

Defining the total angular momentum operator $\hat{L}^{2}$ by

$\hat{L}^{2}:=\hat{L}_{x}^{2}+\hat{L}_{y}^{2}+\hat{L}_{z}^{2},$ (5.33)

we have

$[\hat{L}_{x},\hat{L}^{2}]=[\hat{L}_{y},\hat{L}^{2}]=[\hat{L}_{z},\hat{L}^{2}]=0.$
(5.34)

Besides the operators $\hat{L}_{x}$, $\hat{L}_{y}$ and $\hat{L}_{z}$, there
are several other important operators. A rigid body has body-fixed principal
axes which rotate together with the rigid body, and these local axes of the
rigid body are associated with the principal moments of inertia
$I_{x},I_{y},I_{z}$. To find the rotational energy of the rigid body, we need
the angular momenta projected onto these principal axes, and we consider these
principal axes to be the $x$, $y$ and $z$ axes of the local frame of the rigid
body.

The angular momentum operators in the directions of the body-fixed local axes
are associated with rotations around those local axes. Given an orientation
$\hat{R}_{0}$, the orientation $\hat{R}_{0}$ can be interpreted as the
orientation after rotating by $\hat{R}_{0}$ from the default position
$(\xi=0,\eta=0,\zeta=0,\chi=1)$, i.e., $\alpha=\beta=\gamma=0$, at which point
the body-fixed local axes align with the laboratory axes. Therefore, rotating
around the local $z$ axis of a rigid body by $\epsilon$ results in the
orientation $\hat{R}_{0}\hat{R}_{z;\epsilon}$. This is because the action of
rotating around the local $z$ axis of a rigid body is denoted by
$\hat{R}_{0}\hat{R}_{z;\epsilon}\hat{R}_{0}^{-1}$, which represents the
rotation around the rotated $z$ axis, rotated by $\hat{R}_{0}$, and we have
$\hat{R}_{0}\hat{R}_{z;\epsilon}\hat{R}_{0}^{-1}\hat{R}_{0}=\hat{R}_{0}\hat{R}_{z;\epsilon}$.
The result means that it is the same for the body to first rotate by
$\hat{R}_{0}$ and then rotate around its local $z$ axis, and to first rotate
around the laboratory $z$ axis at the default position and then rotate by
$\hat{R}_{0}$. The orientation after the rotation is given by

$\begin{split}\hat{R}_{0}\hat{R}_{z;\epsilon}&=(\xi_{0},\eta_{0},\zeta_{0},\chi_{0})\cdot(0,0,\frac{\epsilon}{2},1)\\\
&=(\xi_{0}+\frac{\epsilon}{2}\eta_{0},\ \eta_{0}-\frac{\epsilon}{2}\xi_{0},\
\zeta_{0}+\frac{\epsilon}{2}\chi_{0},\
\chi_{0}-\frac{\epsilon}{2}\zeta_{0}).\end{split}$ (5.35)

Therefore, the corresponding angular momentum operator $\hat{Q}_{z}$ which
generates the rotation around the local $z$ axis is given by

$\hat{Q}_{z}=-\frac{i\hbar}{2}\left(\eta\frac{\partial}{\partial\xi}-\xi\frac{\partial}{\partial\eta}+\chi\frac{\partial}{\partial\zeta}-\zeta\frac{\partial}{\partial\chi}\right),$
(5.36)

and the other angular momentum operators are given by

$\hat{Q}_{x}=-\frac{i\hbar}{2}\left(\chi\frac{\partial}{\partial\xi}+\zeta\frac{\partial}{\partial\eta}-\eta\frac{\partial}{\partial\zeta}-\xi\frac{\partial}{\partial\chi}\right),$
(5.37)
$\hat{Q}_{y}=-\frac{i\hbar}{2}\left(-\zeta\frac{\partial}{\partial\xi}+\chi\frac{\partial}{\partial\eta}+\xi\frac{\partial}{\partial\zeta}-\eta\frac{\partial}{\partial\chi}\right).$
(5.38)

The commutation relations among these operators are given by

$\begin{split}[\hat{Q}_{x},\hat{Q}_{y}]&=\hat{Q}_{x}\hat{Q}_{y}-\hat{Q}_{y}\hat{Q}_{x}\\\
&=-\frac{\hbar^{2}}{2}\left(+\chi\frac{\partial}{\partial\zeta}-\zeta\frac{\partial}{\partial\chi}+\eta\frac{\partial}{\partial\xi}-\xi\frac{\partial}{\partial\eta}\right)\\\
&=-i\hbar\hat{Q}_{z},\end{split}$ (5.39)

and

$\displaystyle=-i\hbar\hat{Q}_{x},\qquad[\hat{Q}_{z},\hat{Q}_{x}]=-i\hbar\hat{Q}_{y},$
(5.40)
$\displaystyle\hat{Q}^{2}:=\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}+\hat{Q}_{z}^{2}\equiv\hat{L}^{2},$
(5.41)
$\displaystyle[\hat{Q}_{x},\hat{Q}^{2}]=[\hat{Q}_{y},\hat{Q}^{2}]=[\hat{Q}_{z},\hat{Q}^{2}]=0.$
(5.42)

We have the relation $\hat{Q}^{2}\equiv\hat{L}^{2}$, indicating that the total
angular momentum is the same for both the local frame and the laboratory
frame.

Comparing Eqs. (5.31) and (5.32) with Eqs. (5.39) and (5.40), one notices that
the angular momentum operators $\hat{Q}_{x}$, $\hat{Q}_{y}$ and $\hat{Q}_{z}$
have a different sign in the commutation rules compared with the operators
$\hat{L}_{x}$, $\hat{L}_{y}$ and $\hat{L}_{z}$. This is because the rotations
generated by $\hat{Q}_{x}$, $\hat{Q}_{y}$ and $\hat{Q}_{z}$ are combined on
the right of $\hat{R}_{0}$, while the rotations generated by $\hat{L}_{x}$,
$\hat{L}_{y}$ and $\hat{L}_{z}$ are combined on the left. Other important
results are

$\begin{split}[\hat{L}_{z},\hat{Q}_{z}]&=-\frac{\hbar^{2}}{4}\left(\left[-\eta\frac{\partial}{\partial\xi}+\xi\frac{\partial}{\partial\eta},\,+\eta\frac{\partial}{\partial\xi}-\xi\frac{\partial}{\partial\eta}\right]+\left[\chi\frac{\partial}{\partial\zeta}-\zeta\frac{\partial}{\partial\chi},\,\chi\frac{\partial}{\partial\zeta}-\zeta\frac{\partial}{\partial\chi}\right]\right)\\\
&=0,\end{split}$ (5.43)

and

$[\hat{Q}_{x},\hat{L}_{y}]=0,$ (5.44)

that is,

$[\hat{Q}_{x/y/z},\hat{L}_{x/y/z}]=0,$ (5.45)

and any of $\hat{Q}_{x}$, $\hat{Q}_{y}$ and $\hat{Q}_{z}$ commutes with all of
$\hat{L}_{x}$, $\hat{L}_{y}$ and $\hat{L}_{z}$. Therefore, the state of a
quantum-mechanical rigid body can be characterized by three quantum numbers
$(\hat{L}^{2},\hat{L}_{z},\hat{Q}_{z})$.

###### Physical Interpretation

Unlike a shapeless particle, e.g., an elementary particle, which can be
specified by $(\hat{L}^{2},\hat{L}_{z})$, the state of a rigid body needs
three quantum numbers $(\hat{L}^{2},\hat{L}_{z},\hat{Q}_{z})$ to be specified
because besides the angular momentum of the rotational motion, there are
remaining degrees of freedom concerning how the rotation aligns with the
principal axes, i.e., the shape, of the rigid body. For example, a rod can
rotate around its axis of symmetry and align with the laboratory $z$ axis,
having both large $\hat{L}_{z}$ and $\hat{Q}_{z}$ angular momenta; however,
the axis of symmetry of the rod may rotate on the $x$-$y$ plane of the
laboratory frame as well, so that the rod can have a large $\hat{L}_{z}$
angular momentum while having a zero $\hat{Q}_{z}$ angular momentum.
Therefore, it is necessary to specify both the angular momentum in the
laboratory frame and how the angular momentum aligns with the shape, i.e., the
local axes, of the rigid body. This explains the reason why $\hat{L}_{z}$ and
$\hat{Q}_{z}$ commute but they are not independent: they result from the same
angular momentum projected onto the laboratory $z$ axis and onto the rigid
body $z$ axis, and therefore both of them are constrained by the total angular
momentum $\hat{L}^{2}$ and have the same dimension in their matrix
representations.

##### Symmetry and Conservation of Angular Momenta

The rotational energy of a classical rigid body is given by

$\frac{1}{2}\left(I_{x}\omega_{x}^{2}+I_{y}\omega_{y}^{2}+I_{z}\omega_{z}^{2}\right),$
(5.46)

where $\omega_{x}$, $\omega_{y}$ and $\omega_{z}$ are the angular velocities
projected onto the local $x$, $y$ and $z$ axes of the rigid body, and
therefore the Hamiltonian of a quantum-mechanical rigid body is given by

$\hat{H}=\frac{\hat{Q}_{x}^{2}}{2I_{x}}+\frac{\hat{Q}_{y}^{2}}{2I_{y}}+\frac{\hat{Q}_{z}^{2}}{2I_{z}}.$
(5.47)

As $\hat{H}$ commutes with $\hat{L}_{x}$, $\hat{L}_{y}$ and $\hat{L}_{z}$, the
angular momenta in the laboratory frame and the total angular momentum
$\hat{L}^{2}=\hat{Q}^{2}$ are conserved.

The conservation of $\hat{Q}_{x}$, $\hat{Q}_{y}$ and $\hat{Q}_{z}$ depends on
the symmetry of the rigid body, i.e., the relations among $I_{x}$, $I_{y}$ and
$I_{z}$. For a spherically symmetric rigid body, we have
$I:=I_{x}=I_{y}=I_{z}$ and therefore

$\hat{H}=\dfrac{\hat{Q}^{2}}{2I}=\dfrac{\hat{L}^{2}}{2I}$ (5.48)

and all angular momenta are conserved.

For an axially symmetric rigid body, we have $I_{\perp}:=I_{x}=I_{y}$ and
$I_{\parallel}:=I_{z}$ and

$\begin{split}\hat{H}&=\frac{\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}}{2I_{\perp}}+\frac{\hat{Q}_{z}^{2}}{2I_{\parallel}}\\\
&=\frac{1}{2}\left(\frac{1}{I_{\parallel}}-\frac{1}{I_{\perp}}\right)\hat{Q}_{z}^{2}+\frac{\hat{Q}^{2}}{2I_{\perp}},\end{split}$
(5.49)

and therefore we have

$[\hat{Q}_{z},\hat{H}]=\left[\hat{Q}_{z},\,\frac{1}{2}\left(\frac{1}{I_{\parallel}}-\frac{1}{I_{\perp}}\right)\hat{Q}_{z}^{2}\right]+\left[\hat{Q}_{z},\,\frac{\hat{Q}^{2}}{2I_{\perp}}\right]=0,$
(5.50)

showing that the angular momentum $\hat{Q}_{z}$ is conserved. In this case,
the angular momentum operators $\hat{L}^{2}$, $\hat{L}_{z}$ and $\hat{Q}_{z}$
fully diagonalize the Hamiltonian. For a generic asymmetric rigid body,
$\hat{H}$ does not commute with $\hat{Q}_{x}$, $\hat{Q}_{y}$ or $\hat{Q}_{z}$
and none of the angular momenta $\hat{Q}_{x}$, $\hat{Q}_{y}$ and $\hat{Q}_{z}$
is conserved. These results are consistent with the case for a classical rigid
body.

##### Expressions of the Angular Momentum Operators and the Hamiltonian in
terms of Angles

To find the expressions in terms of straightforwardly physically relevant
variables, we need to take the results from the coordinates of
$(\xi,\eta,\zeta,\chi)$ to the coordinates of angles and distances. First, to
compensate for the redundant degree of freedom in the quaternion
parametrization, we introduce an additional degree of freedom $r$ and define

$\begin{split}\xi&=r\sin\frac{\beta}{2}\cos\dfrac{\alpha-\gamma}{2}\\\
\eta&=r\sin\frac{\beta}{2}\sin\frac{\alpha-\gamma}{2}\\\
\zeta&=r\cos\frac{\beta}{2}\sin\frac{\alpha+\gamma}{2}\\\
\chi&=r\cos\frac{\beta}{2}\cos\frac{\alpha+\gamma}{2},\end{split}$ (5.51)

so that the total number of degrees of freedom becomes 4, and we
have$\xi^{2}+\eta^{2}+\zeta^{2}+\chi^{2}=r^{2}$. Next, for simplicity, we
define

$\mu:=\frac{\alpha-\gamma}{2},\quad\nu:=\frac{\alpha+\gamma}{2},$ (5.52)

and we have

$\alpha\equiv\nu+\mu,\quad\gamma\equiv\nu-\mu,$ (5.53)

and the relations between coordinates $(\xi,\eta,\zeta,\chi)$ and
$(r,\beta,\mu,\nu)$ are given by

$\begin{split}\xi&=r\sin\frac{\beta}{2}\cos\mu,\\\
\eta&=r\sin\frac{\beta}{2}\sin\mu,\\\ \zeta&=r\cos\frac{\beta}{2}\sin\nu,\\\
\chi&=r\cos\frac{\beta}{2}\cos\nu,\end{split}$ (5.54)

and

$\begin{split}r&=\sqrt{\xi^{2}+\eta^{2}+\zeta^{2}+\chi^{2}},\\\
\beta&=2\arcsin\sqrt{\frac{\xi^{2}+\eta^{2}}{\xi^{2}+\eta^{2}+\zeta^{2}+\chi^{2}}},\\\
\mu&=\arctan\frac{\eta}{\xi},\\\ \nu&=\arctan\frac{\zeta}{\chi}.\end{split}$
(5.55)

To transform from the coordinates from $(\xi,\eta,\zeta,\chi)$ to
$(r,\beta,\mu,\nu)$, we need to replace differentials like
$\dfrac{\partial}{\partial\xi}$ with $\dfrac{\partial
r}{\partial\xi}\dfrac{\partial}{\partial
r}+\dfrac{\partial\beta}{\partial\xi}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\mu}{\partial\xi}\dfrac{\partial}{\partial\mu}+\dfrac{\partial\nu}{\partial\xi}\dfrac{\partial}{\partial\nu}$
by utilizing the following relations:

$\begin{split}\dfrac{\partial
r}{\partial\xi}=\frac{\xi}{r},\qquad\dfrac{\partial\beta}{\partial\xi}&=\frac{2}{\tan\frac{\beta}{2}}\frac{\xi}{r^{2}},\qquad\
\
\dfrac{\partial\mu}{\partial\xi}=-\frac{\sin\mu}{r\sin\frac{\beta}{2}},\quad\dfrac{\partial\nu}{\partial\xi}=0,\\\
\dfrac{\partial
r}{\partial\eta}=\frac{\eta}{r},\qquad\dfrac{\partial\beta}{\partial\eta}&=\frac{2}{\tan\frac{\beta}{2}}\frac{\eta}{r^{2}},\qquad\
\
\dfrac{\partial\mu}{\partial\eta}=\frac{\cos\mu}{r\sin\frac{\beta}{2}},\quad\
\ \dfrac{\partial\nu}{\partial\eta}=0,\\\ \dfrac{\partial
r}{\partial\zeta}=\frac{\zeta}{r},\qquad\dfrac{\partial\beta}{\partial\zeta}&=-2\tan\frac{\beta}{2}\,\frac{\zeta}{r^{2}},\quad\dfrac{\partial\mu}{\partial\zeta}=0,\qquad\qquad\dfrac{\partial\nu}{\partial\zeta}=\frac{\cos\nu}{r\cos\frac{\beta}{2}},\\\
\dfrac{\partial
r}{\partial\chi}=\frac{\chi}{r},\qquad\dfrac{\partial\beta}{\partial\chi}&=-2\tan\frac{\beta}{2}\,\frac{\chi}{r^{2}},\quad\dfrac{\partial\mu}{\partial\chi}=0,\qquad\qquad\dfrac{\partial\nu}{\partial\chi}=-\frac{\sin\nu}{r\cos\frac{\beta}{2}},\end{split}$
(5.56)

where for simplicity, $\xi$, $\eta$, $\zeta$, and $\chi$ are considered as
functions of $(r,\beta,\mu,\nu)$ on the right-hand sides of the equations.

Then, we can expand the expression of the angular momentum operator
$\hat{Q}_{z}$

$\begin{split}\hat{Q}_{z}=&-\frac{i\hbar}{2}\left(+\eta\frac{\partial}{\partial\xi}-\xi\frac{\partial}{\partial\eta}+\chi\frac{\partial}{\partial\zeta}-\zeta\frac{\partial}{\partial\chi}\right),\\\
=&-\frac{i\hbar}{2}\left[\eta\left(\dfrac{\partial
r}{\partial\xi}\dfrac{\partial}{\partial
r}+\dfrac{\partial\beta}{\partial\xi}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\mu}{\partial\xi}\dfrac{\partial}{\partial\mu}+\dfrac{\partial\nu}{\partial\xi}\dfrac{\partial}{\partial\nu}\right)\right.\\\
&-\xi\left(\dfrac{\partial r}{\partial\eta}\dfrac{\partial}{\partial
r}+\dfrac{\partial\beta}{\partial\eta}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\mu}{\partial\eta}\dfrac{\partial}{\partial\mu}+\dfrac{\partial\nu}{\partial\eta}\dfrac{\partial}{\partial\nu}\right)\\\
&+\chi\left(\dfrac{\partial r}{\partial\zeta}\dfrac{\partial}{\partial
r}+\dfrac{\partial\beta}{\partial\zeta}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\mu}{\partial\zeta}\dfrac{\partial}{\partial\mu}+\dfrac{\partial\nu}{\partial\zeta}\dfrac{\partial}{\partial\nu}\right)\\\
&\left.-\zeta\left(\dfrac{\partial r}{\partial\chi}\dfrac{\partial}{\partial
r}+\dfrac{\partial\beta}{\partial\chi}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\mu}{\partial\chi}\dfrac{\partial}{\partial\mu}+\dfrac{\partial\nu}{\partial\chi}\dfrac{\partial}{\partial\nu}\right)\right].\end{split}$
(5.57)

The terms involving $\dfrac{\partial}{\partial r}$ in the expression of
$\hat{Q}_{z}$ are found to be the following which cancel out

$\eta\xi\frac{\partial}{\partial r}-\xi\eta\frac{\partial}{\partial
r}+\chi\zeta\frac{\partial}{\partial r}-\zeta\chi\frac{\partial}{\partial
r}=0,$ (5.58)

and therefore, we can show that the degree of freedom $r$ is irrelevant to
rotations. The same holds true for $\hat{Q}_{x}$, $\hat{Q}_{y}$,
$\hat{L}_{x}$, $\hat{L}_{y}$ and $\hat{L}_{z}$, all of which are independent
of $r$. After calculation, the angular momentum operators $\hat{Q}_{x}$,
$\hat{Q}_{y}$ and $\hat{Q}_{z}$ are given by

$\begin{split}\hat{Q}_{z}=&-\frac{i\hbar}{2}\left[\eta\left(\dfrac{\partial\beta}{\partial\xi}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\mu}{\partial\xi}\dfrac{\partial}{\partial\mu}\right)-\xi\left(\dfrac{\partial\beta}{\partial\eta}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\mu}{\partial\eta}\dfrac{\partial}{\partial\mu}\right)\right.\\\
&\left.+\chi\left(\dfrac{\partial\beta}{\partial\zeta}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\nu}{\partial\zeta}\dfrac{\partial}{\partial\nu}\right)-\zeta\left(\dfrac{\partial\beta}{\partial\chi}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\nu}{\partial\chi}\dfrac{\partial}{\partial\nu}\right)\right]\\\
=&-\frac{i\hbar}{2}\left(-\frac{\partial}{\partial\mu}+\frac{\partial}{\partial\nu}\right),\end{split}$
(5.59)
$\begin{split}\hat{Q}_{x}=&-\frac{i\hbar}{2}\left[\chi\left(\dfrac{\partial\beta}{\partial\xi}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\mu}{\partial\xi}\dfrac{\partial}{\partial\mu}\right)+\zeta\left(\dfrac{\partial\beta}{\partial\eta}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\mu}{\partial\eta}\dfrac{\partial}{\partial\mu}\right)\right.\\\
&\left.-\eta\left(\dfrac{\partial\beta}{\partial\zeta}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\nu}{\partial\zeta}\dfrac{\partial}{\partial\nu}\right)-\xi\left(\dfrac{\partial\beta}{\partial\chi}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\nu}{\partial\chi}\dfrac{\partial}{\partial\nu}\right)\right]\\\
=&-\frac{i\hbar}{2}\left(2\cos(\nu-\mu)\dfrac{\partial}{\partial\beta}+\frac{\sin(\nu-\mu)}{\tan\frac{\beta}{2}}\frac{\partial}{\partial\mu}+\tan\frac{\beta}{2}\,\sin(\nu-\mu)\frac{\partial}{\partial\nu}\right)\end{split}$
(5.60)

and

$\begin{split}\hat{Q}_{y}=&-\frac{i\hbar}{2}\left[-\zeta\left(\dfrac{\partial\beta}{\partial\xi}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\mu}{\partial\xi}\dfrac{\partial}{\partial\mu}\right)+\chi\left(\dfrac{\partial\beta}{\partial\eta}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\mu}{\partial\eta}\dfrac{\partial}{\partial\mu}\right)\right.\\\
&\left.+\xi\left(\dfrac{\partial\beta}{\partial\zeta}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\nu}{\partial\zeta}\dfrac{\partial}{\partial\nu}\right)-\eta\left(\dfrac{\partial\beta}{\partial\chi}\dfrac{\partial}{\partial\beta}+\dfrac{\partial\nu}{\partial\chi}\dfrac{\partial}{\partial\nu}\right)\right]\\\
=&-\frac{i\hbar}{2}\left(-2\sin(\nu-\mu)\dfrac{\partial}{\partial\beta}+\frac{\cos(\nu-\mu)}{\tan\frac{\beta}{2}}\frac{\partial}{\partial\mu}+\tan\frac{\beta}{2}\,\cos(\nu-\mu)\frac{\partial}{\partial\nu}\right).\end{split}$
(5.61)

The rotational Hamiltonian is then given by

$\hat{H}=\frac{\hat{Q}_{x}^{2}}{2I_{x}}+\frac{\hat{Q}_{y}^{2}}{2I_{y}}+\frac{\hat{Q}_{z}^{2}}{2I_{z}}.$
(5.62)

####  5.3.2 Model

In this section, we describe the model that we consider in our numerical
experiment of cooling of a quantum rigid body, and derive its time-evolution
equation.

We consider a axially symmetric trapped nanorod as in Refs. [107, 108, 109],
with $I_{x}=I_{y}$, and the axis of symmetry of the rigid body is
approximately aligned with the $z$ axis of the laboratory by a trapping
potential. The trap we consider is the optical dipole trap, or the optical
tweezer [110, 107], which uses optical fields to trap the nanoparticle, and
the potential is given by

$V\propto-\left(\vec{l}\cdot\vec{E}\right)^{2},$ (5.63)

where $\vec{l}$ is the vector of the axis of the rod, and $\vec{E}$ is the
electric field.

##### Hamiltonian in terms of Locally Flat Coordinates near $\beta=0$

As the rod is approximately aligned with the $z$ axis of the laboratory frame,
the angle $\beta$ is small, and therefore, we consider a set of coordinates
that approximately describe the local plane around $\beta=0$, i.e., around the
north pole of spherical coordinates. The coordinates that we consider should
straightforwardly describe the position of the head of the rod, as well as the
rotation of the rod around its axis of symmetry. Our choice of the coordinates
is given by the following:

$\begin{split}x:={}&\beta\sin(\mu+\nu)=\beta\sin\alpha,\\\
y:={}&-\beta\cos(\mu+\nu)=-\beta\cos\alpha,\\\ \theta:={}&2\nu.\end{split}$
(5.64)

We transform from the coordinates $(\mu,\nu,\beta)$ or $(\alpha,\beta,\gamma)$
to the coordinates $(x,y,\theta)$. Consider the local plane at the point of
$\beta=0$ of spherical coordinates, as $\beta$ is small, with the azimuthal
angle $\alpha$, the angles $\beta$ and $\alpha$ represent the radius and the
azimuth of the two-dimensional polar coordinates for the local plane around
$\beta=0$, whereas $x$ and $y$ correspond to the usual Euclidean coordinates.

The inverse relations are given by

$\begin{split}\beta&=\sqrt{x^{2}+y^{2}},\\\
\cos\alpha&=\cos(\mu+\nu)=-\frac{y}{\sqrt{x^{2}+y^{2}}},\\\
\sin\alpha&=\sin(\mu+\nu)=\frac{x}{\sqrt{x^{2}+y^{2}}},\\\
\nu&=\frac{\theta}{2},\\\
\gamma&=\nu-\mu=\theta+\arctan\dfrac{x}{y}.\end{split}$ (5.65)

To evaluate the partial differential operators, we use the following
relations:

$\begin{split}\dfrac{\partial
x}{\partial\beta}&=\sin(\mu+\nu),\quad\quad\dfrac{\partial
y}{\partial\beta}=-\cos(\mu+\nu),\ \
\dfrac{\partial\theta}{\partial\beta}=0,\\\ \dfrac{\partial
x}{\partial\mu}&=\beta\cos(\mu+\nu),\quad\dfrac{\partial
y}{\partial\mu}=\beta\sin(\mu+\nu),\quad\dfrac{\partial\theta}{\partial\mu}=0,\\\
\dfrac{\partial x}{\partial\nu}&=\beta\cos(\mu+\nu),\quad\dfrac{\partial
y}{\partial\nu}=\beta\sin(\mu+\nu),\quad\dfrac{\partial\theta}{\partial\nu}=2,\end{split}$
(5.66)

where $\mu$, $\nu$ and $\beta$ on the right-hand sides of the equations are
considered as functions of $x$, $y$ and $\theta$.

Then, we can calculate

$\begin{split}\hat{Q}_{x}=&-{i\hbar}\left[\cos(\nu-\mu)\frac{\partial}{\partial\beta}+\frac{\sin(\nu-\mu)}{2\tan\frac{\beta}{2}}\frac{\partial}{\partial\mu}+\frac{\tan\frac{\beta}{2}\,\sin(\nu-\mu)}{2}\frac{\partial}{\partial\nu}\right]\\\
=&-{i\hbar}\left[\cos(\nu-\mu)\left(\sin(\mu+\nu)\frac{\partial}{\partial
x}-\cos(\mu+\nu)\frac{\partial}{\partial y}\right)\right.\\\
&+\frac{\sin(\nu-\mu)}{2\tan\frac{\beta}{2}}\left(\beta\cos(\mu+\nu)\frac{\partial}{\partial
x}+\beta\sin(\mu+\nu)\frac{\partial}{\partial y}\right)\\\
&\left.+\frac{\tan\frac{\beta}{2}\,\sin(\nu-\mu)}{2}\left(\beta\cos(\mu+\nu)\frac{\partial}{\partial
x}+\beta\sin(\mu+\nu)\frac{\partial}{\partial
y}+2\frac{\partial}{\partial\theta}\right)\right]\\\
=&-{i\hbar}\left[\cos\gamma\left(\frac{x}{\sqrt{x^{2}+y^{2}}}\frac{\partial}{\partial
x}+\frac{y}{\sqrt{x^{2}+y^{2}}}\frac{\partial}{\partial
y}\right)+\sin\gamma\left(-\frac{y}{\sqrt{x^{2}+y^{2}}}\frac{\partial}{\partial
x}+\frac{x}{\sqrt{x^{2}+y^{2}}}\frac{\partial}{\partial y}\right)\right.\\\
&\left.+\left(\frac{1}{2\tan\frac{\beta}{2}}-\frac{1}{\beta}+\frac{\tan\frac{\beta}{2}}{2}\right)\sin\gamma\left(-y\frac{\partial}{\partial
x}+x\frac{\partial}{\partial
y}\right)+\tan\frac{\beta}{2}\sin\gamma\frac{\partial}{\partial\theta}\right].\end{split}$
(5.67)

For simplicity, we define

$g(\beta):=\frac{1}{2\tan\frac{\beta}{2}}-\frac{1}{\beta}+\frac{\tan\frac{\beta}{2}}{2},$
(5.68)

and use

$\beta\equiv\sqrt{x^{2}+y^{2}},$ (5.69)

and we have

$\begin{split}\hat{Q}_{x}=&-{i\hbar}\left[\cos\gamma\left(\frac{x}{\sqrt{x^{2}+y^{2}}}\frac{\partial}{\partial
x}+\frac{y}{\sqrt{x^{2}+y^{2}}}\frac{\partial}{\partial
y}\right)+\sin\gamma\left(-\frac{y}{\sqrt{x^{2}+y^{2}}}\frac{\partial}{\partial
x}+\frac{x}{\sqrt{x^{2}+y^{2}}}\frac{\partial}{\partial y}\right)\right.\\\
&\left.+g(\beta)\sin\gamma\left(-y\frac{\partial}{\partial
x}+x\frac{\partial}{\partial
y}\right)+\tan\frac{\beta}{2}\sin\gamma\frac{\partial}{\partial\theta}\right].\\\
=&-{i\hbar}\left[\left(\frac{x\cos\gamma-y\sin\gamma}{\sqrt{x^{2}+y^{2}}}-g(\beta)y\sin\gamma\right)\frac{\partial}{\partial
x}+\left(\frac{x\sin\gamma+y\cos\gamma}{\sqrt{x^{2}+y^{2}}}+g(\beta)x\sin\gamma\right)\frac{\partial}{\partial
y}\right.\\\
&\left.+\tan\frac{\beta}{2}\sin\gamma\frac{\partial}{\partial\theta}\right]\\\
={}&{i\hbar}\left\\{\left[\sin\theta+y(x\cos\theta+y\sin\theta)\frac{g(\beta)}{\beta}\right]\frac{\partial}{\partial
x}-\left[\cos\theta+x(x\cos\theta+y\sin\theta)\frac{g(\beta)}{\beta}\right]\frac{\partial}{\partial
y}\right.\\\
&\left.-(x\cos\theta+y\sin\theta)\frac{\tan\frac{\beta}{2}}{\beta}\frac{\partial}{\partial\theta}\right\\}.\end{split}$
(5.70)

The angular momentum operator $\hat{Q}_{y}$ is given by

$\begin{split}\hat{Q}_{y}={}&{i\hbar}\left\\{\left[\cos\theta+y(y\cos\theta-x\sin\theta)\frac{g(\beta)}{\beta}\right]\frac{\partial}{\partial
x}+\left[\sin\theta-x(y\cos\theta-x\sin\theta)\frac{g(\beta)}{\beta}\right]\frac{\partial}{\partial
y}\right.\\\
&\left.-(y\cos\theta-x\sin\theta)\frac{\tan\frac{\beta}{2}}{\beta}\frac{\partial}{\partial\theta}\right\\},\end{split}$
(5.71)

whereas $\hat{Q}_{z}$ is simply given by

$\begin{split}\hat{Q}_{z}=&-{i\hbar}\frac{\partial}{\partial\theta}.\end{split}$
(5.72)

It can be confirmed that the commutation relations

$[\hat{Q}_{x},\hat{Q}_{y}]=-i\hbar\hat{Q}_{z},\qquad[\hat{Q}_{y},\hat{Q}_{z}]=-i\hbar\hat{Q}_{x},\qquad[\hat{Q}_{z},\hat{Q}_{x}]=-i\hbar\hat{Q}_{y}$
(5.73)

indeed hold.

To calculate the Hamiltonian, we need to evaluate the squared operators
$\hat{Q}_{x}^{2}$, $\hat{Q}_{y}^{2}$ and $\hat{Q}_{z}^{2}$ which involve
multiple differentiations. To proceed, for simplicity, we define

$h:=x\cos\theta+y\sin\theta,\qquad k:=y\cos\theta-x\sin\theta,$ (5.74)

and

$g_{1}:=\frac{1}{x}\left(\frac{\partial}{\partial
x}\left(\frac{g}{\beta}\right)\right)\equiv\frac{1}{y}\left(\frac{\partial}{\partial
y}\left(\frac{g}{\beta}\right)\right),\qquad
g_{2}:=\frac{1}{x}\left(\frac{\partial}{\partial
x}\left(\frac{\tan\frac{\beta}{2}}{\beta}\right)\right)\equiv\frac{1}{y}\left(\frac{\partial}{\partial
y}\left(\frac{\tan\frac{\beta}{2}}{\beta}\right)\right),$ (5.75)

and the expressions for $\hat{Q}_{x}^{2}$ and $\hat{Q}_{y}^{2}$ are explicitly
given by

$\begin{split}\hat{Q}_{x}^{2}=&-\hbar^{2}\left\\{\left[\sin\theta+yh\frac{g}{\beta}\right]^{2}\partial_{x}^{2}+\left[\cos\theta+xh\frac{g}{\beta}\right]^{2}\partial_{y}^{2}+h^{2}\left(\frac{\tan\frac{\beta}{2}}{\beta}\right)^{2}\partial_{\theta}^{2}\right.\\\
&\qquad-2\left[\sin\theta\cos\theta+\left((y^{2}+x^{2})\sin\theta\cos\theta+xy\right)\frac{g}{\beta}+xyh^{2}\left(\frac{g}{\beta}\right)^{2}\right]\partial_{x}\partial_{y}\\\
&\qquad-\left(\sin\theta+yh\frac{g}{\beta}\right)h\frac{2\tan\frac{\beta}{2}}{\beta}\partial_{x}\partial_{\theta}+\left(\cos\theta+xh\frac{g}{\beta}\right)h\frac{2\tan\frac{\beta}{2}}{\beta}\partial_{y}\partial_{\theta}\\\
&\qquad-h\left[\cos\theta\left(\frac{g}{\beta}+\frac{\tan\frac{\beta}{2}}{\beta}\right)+\left(\left(x^{2}-y^{2}\right)\cos\theta+2xy\sin\theta\right)\left(\frac{g}{\beta}\right)^{2}+yk\left(g_{1}+\frac{\tan\frac{\beta}{2}}{\beta}\frac{g}{\beta}\right)\right]\partial_{x}\\\
&\qquad-h\left[\sin\theta\left(\frac{g}{\beta}+\frac{\tan\frac{\beta}{2}}{\beta}\right)-\left(\left(x^{2}-y^{2}\right)\sin\theta-2xy\cos\theta\right)\left(\frac{g}{\beta}\right)^{2}-xk\left(g_{1}+\frac{\tan\frac{\beta}{2}}{\beta}\frac{g}{\beta}\right)\right]\partial_{y}\\\
&\left.\qquad-
hk\left[\frac{g\tan\frac{\beta}{2}}{\beta^{2}}-g_{2}-\left(\frac{\tan\frac{\beta}{2}}{\beta}\right)^{2}\right]\partial_{\theta}\right\\},\end{split}$
(5.76)
$\begin{split}\hat{Q}_{y}^{2}=&-\hbar^{2}\left\\{\left[\cos\theta+yk\frac{g}{\beta}\right]^{2}\partial_{x}^{2}+\left[\sin\theta-
xk\frac{g}{\beta}\right]^{2}\partial_{y}^{2}+k^{2}\left(\frac{\tan\frac{\beta}{2}}{\beta}\right)^{2}\partial_{\theta}^{2}\right.\\\
&\qquad+2\left[\sin\theta\cos\theta+\left((y^{2}+x^{2})\sin\theta\cos\theta-
xy\right)\frac{g}{\beta}-xyk^{2}\left(\frac{g}{\beta}\right)^{2}\right]\partial_{x}\partial_{y}\\\
&\qquad-\left(\cos\theta+yk\frac{g}{\beta}\right)k\frac{2\tan\frac{\beta}{2}}{\beta}\partial_{x}\partial_{\theta}-\left(\sin\theta-
xk\frac{g}{\beta}\right)k\frac{2\tan\frac{\beta}{2}}{\beta}\partial_{y}\partial_{\theta}\\\
&\qquad+k\left[\sin\theta\left(\frac{g}{\beta}+\frac{\tan\frac{\beta}{2}}{\beta}\right)+\left(\left(x^{2}-y^{2}\right)\sin\theta-2xy\cos\theta\right)\left(\frac{g}{\beta}\right)^{2}+yh\left(g_{1}+\frac{\tan\frac{\beta}{2}}{\beta}\frac{g}{\beta}\right)\right]\partial_{x}\\\
&\qquad-k\left[\cos\theta\left(\frac{g}{\beta}+\frac{\tan\frac{\beta}{2}}{\beta}\right)-\left(\left(x^{2}-y^{2}\right)\cos\theta+2xy\sin\theta\right)\left(\frac{g}{\beta}\right)^{2}+xh\left(g_{1}+\frac{\tan\frac{\beta}{2}}{\beta}\frac{g}{\beta}\right)\right]\partial_{y}\\\
&\left.\qquad+hk\left[\frac{g\tan\frac{\beta}{2}}{\beta^{2}}-g_{2}-\left(\frac{\tan\frac{\beta}{2}}{\beta}\right)^{2}\right]\partial_{\theta}\right\\}.\end{split}$
(5.77)

We consider an axially symmetric rigid body, so that we have
$I_{\perp}:=I_{x}=I_{y}$ and $I_{\parallel}:=I_{z}$, and the Hamiltonian is
given by

$\hat{H}=\frac{\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}}{2I_{\perp}}+\frac{\hat{Q}_{z}^{2}}{2I_{\parallel}}+V.$
(5.78)

Here the term $\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}$ is calculated to be

$\begin{split}\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}=&-\hbar^{2}\left\\{\left[1+2y^{2}\frac{g}{\beta}+y^{2}g^{2}\right]^{2}\partial_{x}^{2}+\left[1+2x^{2}\frac{g}{\beta}+x^{2}g^{2}\right]^{2}\partial_{y}^{2}+\left(\tan\frac{\beta}{2}\right)^{2}\partial_{\theta}^{2}\right.\\\
&\qquad-2xy\left[\frac{2g}{\beta}+g^{2}\right]\partial_{x}\partial_{y}-y\left(1+g\beta\right)\frac{2\tan\frac{\beta}{2}}{\beta}\partial_{x}\partial_{\theta}+x\left(1+g\beta\right)\frac{2\tan\frac{\beta}{2}}{\beta}\partial_{y}\partial_{\theta}\\\
&\left.\qquad-x\left(\frac{g}{\beta}+\frac{\tan\frac{\beta}{2}}{\beta}+g^{2}\right)\partial_{x}-y\left(\frac{g}{\beta}+\frac{\tan\frac{\beta}{2}}{\beta}+g^{2}\right)\partial_{y}\right\\},\end{split}$
(5.79)

which is the main result of our derivation. The expression does not involve
explicit functions of $\theta$, which is consistent with the fact that the
Hamiltonian commutes with $\hat{Q}_{z}$.

The last term $\hat{Q}_{z}^{2}$ in the Hamiltonian is given by

$\hat{Q}_{z}^{2}=-\hbar^{2}\partial_{\theta}^{2}.$ (5.80)

##### Properties of the Hamiltonian

To see the properties of the Hamiltonian, we calculate the leading terms of
the Taylor series of $g$ and $\tan\frac{\beta}{2}$. They are given by

$g(\beta)=\frac{\beta}{6}+\frac{7}{360}\beta^{3}+\frac{31}{15120}\beta^{5}+\frac{127}{604800}\beta^{7}+\frac{73}{3421440}\beta^{9}+O(\beta^{11}),$
(5.81)

and

$\tan\frac{\beta}{2}=\frac{\beta}{2}+\frac{\beta^{3}}{24}+\frac{\beta^{5}}{240}+\frac{17}{40320}\beta^{7}+\frac{31}{725760}\beta^{9}+O(\beta^{11}),$
(5.82)

showing that the high-order terms decay rapidly, as the angle $\beta$ is the
polar angle of spherical coordinates, which necessarily satisfies
$\beta\in[0,\pi]$.

Up to the second order of $\beta$, the term $\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}$
is given by

$\begin{split}\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}=&-\hbar^{2}\left[\left(1+\frac{y^{2}}{3}\right)\frac{\partial^{2}}{\partial
x^{2}}+\left(1+\frac{x^{2}}{3}\right)\frac{\partial^{2}}{\partial
y^{2}}+\frac{x^{2}+y^{2}}{4}\frac{\partial^{2}}{\partial\theta^{2}}-\frac{2xy}{3}\frac{\partial^{2}}{\partial
x\,\partial y}\right.\\\ &\left.-y\frac{\partial^{2}}{\partial
x\,\partial\theta}+x\frac{\partial^{2}}{\partial
y\,\partial\theta}-\frac{2x}{3}\frac{\partial}{\partial
x}-\frac{2y}{3}\frac{\partial}{\partial y}\right]+O(\beta^{3}),\end{split}$
(5.83)

where we have used $x\sim O(\beta)$ and $y\sim O(\beta)$ because of
$\beta=\sqrt{x^{2}+y^{2}}$. When $\beta$ is small, the state of the rigid body
is localized around $\beta=0$ in a small region, and, therefore, $x$ and $y$
become small, whereas $\partial_{x}$ and $\partial_{y}$ can become large. If
we keep the terms involving $\partial_{\theta}$ and ignore the other small
terms, the above equation becomes

$\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}\approx-\hbar^{2}\left(\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial
y^{2}}+\frac{x^{2}+y^{2}}{4}\frac{\partial^{2}}{\partial\theta^{2}}-y\frac{\partial^{2}}{\partial
x\,\partial\theta}+x\frac{\partial^{2}}{\partial y\,\partial\theta}\right),$
(5.84)

which is equivalent to the Hamiltonian of a charged particle in a magnetic
field. As the Hamiltonian commutes with
$\hat{Q}_{z}\equiv-{i\hbar}{\partial_{\theta}}$, we assume that the state we
investigate is an eigenstate of $\hat{Q}_{z}$ so that we can regard
$\hat{Q}_{z}$ as a constant. Then, we map Eq. (5.84) to the Hamiltonian of a
charged particle in a magnetic field. Since the Hamiltonian of a charged
particle in a magnetic field is given by

$H=\frac{1}{2m}\left(\vec{p}-\frac{Q}{c}\vec{A}\right)^{2},$ (5.85)

the term $\frac{Q}{c}\vec{A}$ in our case should be

$\frac{Q}{c}\vec{A}=\left(-\frac{y}{2}\hat{Q}_{z},\,+\frac{x}{2}\hat{Q}_{z},\,0\right),$
(5.86)

and we have

$\nabla\times\frac{Q}{c}\vec{A}=\left(0,\,0,\,\hat{Q}_{z}\right),$ (5.87)

which corresponds to a magnetic field that drives the particle to rotate
counter-clockwise, given a positive value of $\hat{Q}_{z}$. Therefore, up to
the lowest order with respect to a small $\beta$, the system of the quantum-
mechanical rigid body behaves in the same way as a charged particle in a
magnetic field, and the strength of the magnetic field is determined by the
angular momentum $\hat{Q}_{z}$.

##### Metric, Hermiticity and Momentum Operators

The space of Euler angles is non-Euclidean, and when we integrate a function
$f$ over the space, we should include the factor resulting from the metric,
and the integration is given by

$\int_{0}^{2\pi}\int_{0}^{2\pi}\int_{0}^{\pi}f\sin\beta\,d\beta\,d\alpha\,d\gamma.$
(5.88)

In terms of the coordinates $(\beta,\mu,\nu)$, the integration is given by

$\int_{0}^{\pi}\int_{\mu}^{\mu+2\pi}\int_{0}^{\pi}f\cdot
2\sin\beta\,d\beta\,d\nu\,d\mu.$ (5.89)

The determinant of the Jacobian matrix given by Eq. (5.66) is equal to
$2\beta$, and therefore, in terms of the coordinates $(x,y,\theta)$, the
integration is given by

$\int_{0}^{2\pi}\left(\iint_{\sqrt{x^{2}+y^{2}}\leq\pi}f\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx\,dy\right)\,d\theta,$
(5.90)

and we see that the metric term
$\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}$ is approximately equal to
$1$ when both $x$ and $y$ are small, showing that the space is indeed
approximately flat near $\beta=0$.

Due to the nontrivial metric term that one needs to take into account when one
performs the integration, the inner product $\langle\phi|\psi\rangle$ should
also be calculated using

$\int_{0}^{2\pi}\iint_{\sqrt{x^{2}+y^{2}}\leq\pi}\phi^{*}\psi\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx\,dy\,d\theta,$
(5.91)

and, therefore, usual operators like $-{i}{\hbar}\frac{\partial}{\partial x}$
and $-{i}{\hbar}\frac{\partial}{\partial y}$ are no longer hermitian, as one
generally has

$\begin{split}\int_{x_{0}}^{x_{1}}\phi^{*}\left(\frac{\partial}{\partial
x}\psi\right)\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx\neq-\int_{x_{0}}^{x_{1}}\left(\frac{\partial}{\partial
x}\phi^{*}\right)\psi\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx\end{split}$
(5.92)

even if $\phi$ and $\psi$ vanish at $x_{0}$ and $x_{1}$, subject to the
vanishing boundary condition. Instead, integration by parts yields

$\int_{x_{0}}^{x_{1}}\phi^{*}\left(\frac{\partial}{\partial
x}\psi\right)\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx=-\int_{x_{0}}^{x_{1}}\left(\frac{\partial}{\partial
x}\left(\phi^{*}\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\right)\right)\psi\,dx,$
(5.93)

and an additional term $\phi^{*}\left(\frac{\partial}{\partial
x}\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\right)$ appears.

With the nontrivial metric, the hermitian momentum operators conjugate to the
position operators $x$ and $y$ can be given by

$\hat{p}_{x}:=-i\hbar\left(\frac{\partial}{\partial
x}+\frac{1}{2}\left(\frac{\partial}{\partial
x}\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\right)\right),$ (5.94)
$\hat{p}_{y}:=-i\hbar\left(\frac{\partial}{\partial
y}+\frac{1}{2}\left(\frac{\partial}{\partial
y}\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\right)\right).$ (5.95)

To confirm the hermiticity, using the vanishing boundary condition and
integration by parts, we calculate

$\begin{split}&\int_{x_{0}}^{x_{1}}\phi^{*}\left(\hat{p}_{x}\psi\right)\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx\\\
=&\int_{x_{0}}^{x_{1}}\phi^{*}\left(-i\hbar\left(\frac{\partial}{\partial
x}\psi+\frac{1}{2}\left(\frac{\partial}{\partial
x}\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\right)\psi\right)\right)\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx\\\
=&\int_{x_{0}}^{x_{1}}i\hbar\left(\frac{\partial}{\partial
x}\left(\phi^{*}\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\right)\right)\psi\,dx-\frac{i\hbar}{2}\int_{x_{0}}^{x_{1}}\phi^{*}\psi\left(\frac{\partial}{\partial
x}\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\right)\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx\\\
=&\int_{x_{0}}^{x_{1}}i\hbar\left(\frac{\partial}{\partial
x}\phi^{*}\right)\psi\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx+\frac{i\hbar}{2}\int_{x_{0}}^{x_{1}}\phi^{*}\psi\left(\frac{\partial}{\partial
x}\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\right)\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx\\\
=&\int_{x_{0}}^{x_{1}}i\hbar\left(\frac{\partial}{\partial
x}\phi^{*}+\frac{1}{2}\left(\frac{\partial}{\partial
x}\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\right)\phi^{*}\right)\psi\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx\\\
=&\int_{x_{0}}^{x_{1}}\left(\hat{p}_{x}\phi\right)^{*}\psi\cdot\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\,dx,\end{split}$
(5.96)

which confirms the hermiticity.

The momentum operators $\hat{p}_{x}$ and $\hat{p}_{y}$ satisfy the usual
commutation relations

$\left[\hat{p}_{x},\hat{p}_{y}\right]=0,\quad\left[x,\hat{p}_{x}\right]=i\hbar,\quad\left[y,\hat{p}_{y}\right]=i\hbar,\quad\left[x,\hat{p}_{y}\right]=\left[y,\hat{p}_{x}\right]=0.$
(5.97)

We use these operators to compute the distribution moments of the wave
function in our numerical algorithm in the following sections.

##### Control Fields

In the following, we proceed to describe the model of the trapped and
controlled rigid body that we consider.

We consider a nanorod trapped by an electromagnetic field through dipole-
induced dipole interaction using a laser, and the potential is given by Eq.
(5.63), i.e.

$V\propto-\left(\vec{l}\cdot\vec{E}\right)^{2},$ (5.98)

where $\vec{l}$ is the vector of the axis of the rod, and $\vec{E}$ is the
electric field. In order to control the rotation of the rigid body,
specifically, to control the movement of the rigid body in $x$ and $y$
coordinates discussed above, we consider two additional lasers which can shift
the center of the trapping potential in the space spanned by $x$ and $y$. The
trapping laser has a fixed intensity $E_{z}$ with the polarization direction
$z$, and the two control lasers have tunable intensities $E_{x}$ and $E_{y}$
with the polarization directions $x$ and $y$. The two control lasers are
arranged perpendicularly, and the three laser beams are phased-locked and
arranged on the same plane, the $x$-$y$ plane, and they intersect at the
position of the rigid body, creating a field $\vec{E}=(E_{x},E_{y},E_{y})$ at
the rigid body. Then, we use $E_{x}$ and $E_{y}$ as the relevant control
variables to control the time evolution of the system and cool the system.

The position of the head of the nanorod relative to the center of the rod is
given by

$\left(\frac{xl}{2}\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}},\,\frac{yl}{2}\frac{\sin\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}},\,\frac{l}{2}\cos\sqrt{x^{2}+y^{2}}\right),$
(5.99)

where $l$ is the length of the rigid body. Therefore, the potential is given
by

$V=-\left(x\frac{\sin\beta}{\beta}E_{x}+y\frac{\sin\beta}{\beta}E_{y}+\cos\beta
E_{z}\right)^{2},\qquad\beta\equiv\sqrt{x^{2}+y^{2}},$ (5.100)

where the other coefficients have been absorbed into the coefficients $E_{x}$,
$E_{y}$ and $E_{z}$. With $E_{x}=E_{y}=0$, up to the lowest order in $\beta$,
we have

$V=-E_{z}^{2}+E_{z}^{2}\beta^{2}=\frac{k}{2}(x^{2}+y^{2})-E_{z}^{2}+O(\beta^{4}),\qquad
k:=2E_{z}^{2},$ (5.101)

where we have used $x\sim O(\beta)$ and $y\sim O(\beta)$, and $k$ is the
strength of the trapping potential. Therefore, in the neighbourhood of
$\beta=0$, the trapping potential is the standard harmonic potential up to a
shift of the constant $-E_{z}^{2}$.

With nonzero $E_{x}$ and $E_{y}$, we have

$V=E_{z}^{2}\left[\left(x-\frac{E_{x}}{E_{z}}\right)^{2}+\left(y-\frac{E_{y}}{E_{z}}\right)^{2}\right]-E_{z}^{2}-E_{x}^{2}-E_{y}^{2}+O(\beta^{3}),$
(5.102)

where we have assumed $\dfrac{E_{x}}{E_{z}}\sim O(\beta)$ and
$\dfrac{E_{y}}{E_{z}}\sim O(\beta)$. The result shows that the center of the
potential is shifted by $\dfrac{E_{x}}{E_{z}}$ and $\dfrac{E_{y}}{E_{z}}$ in
the $x$ and $y$ directions.

##### Measurements

To our knowledge, so far there is not any widely accepted model of quantum-
mechanical measurement for a rigid body system, and in our model, we simply
assume that the measurement is Gaussian with respect to the $x$ and $y$
coordinates, which represent the position of the head of the quantum-
mechanical rod, and we assume that the measurement efficiency is unity.

The measurement can be regarded as Gaussian if the following assumptions hold
true: (1) the polar angle $\beta$ is small so that the space spanned by $x$
and $y$ is isotropic; (2) the measurement is week and repetitive so that it
can be regarded as approximately continuous; (3) the state of the rigid body
is sufficiently localized and the variance of measurement outcomes is not
large, so that one can take the average of measurement outcomes without
ambiguity in the space of angles.

One simple example of the measurement that measures the position of the head
of the nanorod is given in Ref. [107], where a probe light shines at the rigid
body in the $z$ direction and the direction of the scattered light is probed.
Because the rigid body is axially symmetric, the scattered light does not
provide information on the angle $\theta$, and it only provides information
about $x$ and $y$, which represent the position of the head of the trapped
rod. As both $x$ and $y$ are measured, the measurement backaction perturbs the
state in both the $x$ and $y$ directions simultaneously, and we regard the
perturbations in the $x$ and $y$ directions as independent.

The model of measurement completes the description of our quantum-mechanical
model of the controlled rigid body. The stochastic time-evolution equation of
the state is given by

$\displaystyle
d|\psi\rangle=\left[\left(-\frac{i}{\hbar}\hat{H}-\frac{\gamma}{4}(\hat{x}-\langle\hat{x}\rangle)^{2}-\frac{\gamma}{4}(\hat{y}-\langle\hat{y}\rangle)^{2}\right)dt+\sqrt{\dfrac{\gamma}{2}}(\hat{x}-\langle\hat{x}\rangle)dW_{1}+\sqrt{\dfrac{\gamma}{2}}(\hat{y}-\langle\hat{y}\rangle)dW_{2}\right]|\psi\rangle,$
(5.103)
$\displaystyle\hat{H}=\frac{\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}}{2I_{\perp}}+\frac{\hat{Q}_{z}^{2}}{2I_{\parallel}}+V,$
(5.104)

where $\gamma$ is the measurement strength, $I_{\perp}$ and $I_{\parallel}$
are the moments of inertia, $dW_{1}$ and $dW_{2}$ are independent Wiener
increments, i.e. random variables satisfying $dW_{1}\sim\mathcal{N}(0,dt)$ and
$dW_{2}\sim\mathcal{N}(0,dt)$, following the convention discussed in Section
2.2.3. The terms $\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}$, $\hat{Q}_{z}^{2}$ and $V$
are given by Eqs. (5.79), (5.80) and (5.100), respectively.

##### Settings of the Numerical Algorithm

We use Eq. (5.103) in our numerical simulation of the controlled quantum-
mechanical nanorod, and we do not use any approximation with respect to the
polar angle $\beta$, or, with respect to $x$ and $y$. We assume that the state
is an eigenstate of $\hat{Q}_{z}$ and we set the value of
$\hat{Q}_{z}\equiv-{i\hbar}{\partial_{\theta}}$ to be a constant, and we
simulate the time evolution of the wave function in the two-dimensional space
spanned by $x$ and $y$. As in Section 5.2, we use finite difference methods,
discretizing the continuous space and time into discrete sites and time steps,
so as to simulate the time evolution of the state.

As discussed in the section of the properties of the Hamiltonian, when $\beta$
is very small, the system is approximately linear, because if we define

$\hat{p}_{x}:=-i\hbar\frac{\partial}{\partial{x}},\qquad\hat{p}_{y}:=-i\hbar\frac{\partial}{\partial{y}},$
(5.105)

we can rewrite Eq. (5.84) into

$\begin{split}\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}&\approx-\hbar^{2}\left(\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial
y^{2}}+\frac{x^{2}+y^{2}}{4}\frac{\partial^{2}}{\partial\theta^{2}}-y\frac{\partial^{2}}{\partial
x\,\partial\theta}+x\frac{\partial^{2}}{\partial y\,\partial\theta}\right)\\\
&=\hat{p}_{x}^{2}+\hat{p}_{y}^{2}+\frac{x^{2}+y^{2}}{4}\hat{Q}_{z}^{2}-y\hat{p}_{x}\hat{Q}_{z}+x\hat{p}_{y}\hat{Q}_{z},\end{split}$
(5.106)

which is quadratic with respect to $x$, $y$, $\hat{p}_{x}$ and $\hat{p}_{y}$,
and the time evolutions of $x$, $y$, $\hat{p}_{x}$ and $\hat{p}_{y}$ are thus
linear. Especially, as the high-order terms in Eq. (5.79) are relatively small
as discussed above, we need to ensure that $x$ and $y$ are moderately large so
that nonlinear behaviour of the system can be observed. Therefore, we use the
parameter regime $x,y\sim 0.5$ in our numerical experiments, and we describe
our system parameters in the following. Note that the coordinates $x$ and $y$
are dimensionless because they represent angles.

We simulate the space spanned by $x$ and $y$ in the region of $-1.29\leq x\leq
1.29$ and $-1.29\leq y\leq 1.29$, which approximately corresponds to a polar
angle $\beta$ of $74$ degrees. The spacing between adjacent discrete sites in
the simulated discretized space is set to be $0.03$, and therefore, we
simulate on the $87\times 87$ grid.

Instead of directly setting the system parameters $k$, $I_{\perp}$ and
$\gamma$, we set them indirectly through the parameters in Table 5.2.

|  $\omega$ ($\omega_{c}$) | $\sigma_{g}$ | $\dfrac{\gamma}{k}$ |  $\hat{Q}_{z}$ ($\hbar$) | $\dfrac{\max\left(\sqrt{E_{x}^{2}+E_{y}^{2}}\right)}{E_{z}}$
---|---|---|---|---|---
rigid body | $\pi$ | $0.1$ | $1,2,3$ | $5,80$ | 1
Table 5.2: System parameters of the quantum-mechanical rigid body that we use
in our numerical experiments, in terms of a reference angular momentum
$\omega_{c}$. $\omega$ is the angular frequency of the harmonic potential
given by Eq. (5.101) ignoring the high-order terms and $\sigma_{g}$ is the
standard deviation of the probability distribution of the corresponding ground
state, which is a Gaussian state, of the harmonic potential. We use $\omega$
and $\sigma_{g}$ to find the parameters $k$, $E_{z}$ and $I_{\perp}$ of the
system of the rigid body.

Using the angular frequency $\omega$ of the potential $V$ at the lowest-order
approximation given by Eq. (5.101), and the standard deviation $\sigma_{g}$ of
the probability distribution of the corresponding ground-state wave function,
as given in Table 5.2, we can find the values of $I_{\perp}$ and $k$, given by

$I_{\perp}=\frac{\hbar}{2\omega\sigma_{g}^{2}}$ (5.107)

and

$k=\omega^{2}I_{\perp},$ (5.108)

and we have

$E_{z}=\sqrt{\frac{k}{2}},$ (5.109)

following Eq. (5.101).

Because $k$ and the measurement strength $\gamma$ are of the same physical
dimension, instead of directly setting the value of $\gamma$, we tune the
ratio $\frac{\gamma}{k}$, and we try three different values $1,2$ and $3$ in
our experiments. The angular momentum $\hat{Q}_{z}$ is set to be $5$ or $80$.
When we have $\hat{Q}_{z}=5$, assuming
$I_{\parallel}\approx\dfrac{I_{\perp}}{10}$, the energy of the rotation around
the local $z$ axis, given by $\dfrac{\hat{Q}_{z}}{2I_{\parallel}}$, is
comparable to the rest of the energy of the system, given by
$\dfrac{\hat{Q}_{x}^{2}+\hat{Q}_{y}^{2}}{2I_{\perp}}+V$; when we have
$\hat{Q}_{z}=80$, the rotational energy around the local $z$ axis is hundreds
of times larger than the rest of the energy of the system. We therefore try
the two different cases in our experiments.

The control forces we allow are given by the following discrete set

$\left\\{(E_{x},E_{y})\ |\ \sqrt{E_{x}^{2}+E_{y}^{2}}\leq
E_{z},\quad\frac{E_{x}}{E_{z}}=0.2n_{x},\quad\frac{E_{y}}{E_{z}}=0.2n_{y},\quad
n_{x},n_{y}\in\mathbb{Z}\right\\},$ (5.110)

including 81 different choices of control forces in total.

Whenever the controller determines a choice of control forces, the control
forces are kept constant during a time of $\dfrac{1}{5\omega_{c}}$. After an
evolution time of $\dfrac{1}{5\omega_{c}}$ of the state, the energy
$\langle\hat{H}\rangle$ with $E_{x}=E_{y}=0$ is calculated as the minus reward
for the reinforcement learning agent, and the controller determines a new
choice of control forces on the basis of the information about the
instantaneous wave function.

Most of the other settings of the numerical algorithm are the same as those in
Section 5.2. The discount factor $\gamma$ in Q-learning is tuned to be 0.96,
and we change the numbers of the hidden units in the deep neural network to be
512, 1024, 512.

The input of the neural network is set to be the distribution moments of the
wave function as in our experiments of the case of a quartic oscillator. In
the case of a rigid body, the distribution moments are computed with respect
to the 4 operators $x$, $y$, $\hat{p}_{x}$ and $\hat{p}_{y}$, where
$\hat{p}_{x}$ and $\hat{p}_{y}$ are given by Eqs. (5.94) and (5.95), and we
compute up to the fifth central moments, totally constituting 125 real
numbers, and use them as the input of the neural network. Additionally, we
rescale $x$ and $y$ by a factor of $\sqrt{\frac{2}{k}}$ and rescale
$\hat{p}_{x}$ and $\hat{p}_{y}$ by a factor of $\sqrt{2I_{\perp}}$ to ensure
that the results are of the order of $O(1)$.

As in the case of the simulation of a quartic oscillator, we use the finite
difference methods given by Eqs. (5.13) and (5.14) to evaluate the partial
differential operators. To numerically integrate the stochastic time-evolution
equation given by Eq. (5.103), we use a time step of
$0.00125\times\frac{1}{\omega_{c}}$, and we use the explicit 1.5 order strong
scheme [104] with several modifications. First, we additionally include terms
obtained using the 4th-order Runge-Kutta method to include the high-order
terms of the deterministic time evolution of the state; next, we do not
include the 1.5 order terms of the stochastic evolution that are contributed
by the cross terms of the stochastic noise, $\iiint
dW_{1/2}\,dW_{1/2}\,dW_{1/2}$, because they are very hard calculate, and they
have been confirmed to contribute little to the time evolution, only making a
difference smaller than $0.1\%$ to the state; lastly, we ignore the parts of
the wave function that are too far away from the center of the probability
distribution, setting the amplitudes to be zero for $|x-\langle
x\rangle|>\frac{\pi}{3}$ and $|y-\langle y\rangle|>\frac{\pi}{3}$. The cross
term $\iint dW_{1}\,dW_{2}$ which is necessary for the numerical integration
of the differential equation is evaluated approximately based on series
expansion using the Legendre series, following Ref. [111] (pp. 761), up to 500
terms. Our numerical codes have been publicized and all details can be found
therein.222https://github.com/Z-T-WANG/PhDThesis/tree/main/rigid%20body

####  5.3.3 Evaluation and Results

The settings of episodes and the evaluation of performances are the same as in
the case of our experiments of a quartic oscillator, described in Section
5.2.3, except for several minor differences.

The initial state we use in our numerical experiments is set to be a Gaussian
wave packet centered at a position randomly picked in the region $-0.4\leq
x\leq 0.4$ and $-0.4\leq y\leq 0.4$, and the momentum of the wave packet is
set such that the velocity is zero, i.e. $\dfrac{d\langle x\rangle}{dt}=0$ and
$\dfrac{d\langle y\rangle}{dt}=0$, on the basis of the approximation given by
Eq. (5.84). The highest energy we allow in our simulation is $30\hbar\omega$,
or, $30\pi\hbar\omega_{c}$, beyond which we stop the simulation and consider
the control as having failed. We also stop the simulation when the wave
function gets close to the boundary of the simulated space, when the
probability of being around the boundary is larger than $1.5\times 10^{-3}$.
We use the 4 discrete sites counted from the boundary inwards to calculate
this probability. Regarding training, the agent is trained for 9000 episodes
after it successfully stabilizes the state for a time evolution of
$\dfrac{100}{\omega_{c}}$.

##### The LQG Control

To compare the results obtained using the DQN algorithm and the C-DQN
algorithm with standard control strategies, we consider the linear-quadratic-
Gaussian (LQG) control on the basis of the approximation given by Eqs. (5.84)
and (5.102). The control force is chosen such that after an evolution time of
$\dfrac{1}{5\omega_{c}}$ of the state, the state is expected to satisfy

$\frac{v_{x}}{\langle
x\rangle}=-\sqrt{\frac{k}{I_{\perp}}},\qquad\frac{v_{y}}{\langle
y\rangle}=-\sqrt{\frac{k}{I_{\perp}}},$ (5.111)

where

$v_{x}:=\dfrac{d\langle x\rangle}{dt},\qquad v_{y}:=\dfrac{d\langle
y\rangle}{dt},$ (5.112)

and the LQG controller effectively treats the state as a classical particle,
only taking the variables $\langle x\rangle$, $\langle y\rangle$, $v_{x}$ and
$v_{y}$ into consideration. This controller pushes the particle onto the
physical trajectory that corresponds to the Lagrangian

$L=T-V=\frac{1}{2I_{\perp}}\left(v_{x}^{2}+v_{y}^{2}\right)-\left(-\frac{k}{2}\left(x^{2}+y^{2}\right)\right)=\frac{1}{2I_{\perp}}\left(v_{x}^{2}+v_{y}^{2}\right)+\frac{k}{2}\left(x^{2}+y^{2}\right),$
(5.113)

which is exactly equal to the energy of the system. Therefore, this control
minimizes the energy of the system integrated with respect to time according
to the Hamilton principle.

Although the LQG control is provably optimal [15] when the dynamics is linear
and Gaussian noise is present, it can fail to stabilize the rigid body due to
nonlinearity of the system. Specifically, because the interaction potential is
given by

$V\propto-\left(\vec{l}\cdot\vec{E}\right)^{2},$ (5.114)

when the direction of the control field and the direction of the controlled
rod are perpendicular, the control force becomes $0$, and when the angle
between the two directions exceeds $\frac{\pi}{2}$, the interaction becomes
repulsive. However, the LQG control replies on the linear approximation and
considers that larger values of $E_{x}$ and $E_{y}$ should always drive the
particle in the $x$ and $y$ directions more strongly, which is not necessarily
true, as the particle can move away from the center and the angle between
$\vec{l}$ and $\vec{E}$ can become large. Therefore, in order to prevent the
systematic failure of the LQG control in controlling the rigid body, we put
constraints on $E_{x}$ and $E_{y}$ such that we have $\left|\langle
x\rangle-\frac{E_{x}}{E_{z}}\right|\leq\frac{\pi}{4}$ and $\left|\langle
y\rangle-\frac{E_{y}}{E_{z}}\right|\leq\frac{\pi}{4}$, and we do not put these
constraints on the AI-based controllers. With these constraints, the LQG
control can stabilize the rigid body for a long time in our numerical
experiments. For a fair comparison with the AI-based controllers, we also
discretize the control forces of the LQG control, mapping the control forces
of the LQG control to the nearest choice in the set given by Eq. (5.110).

##### Data Augmentation

Because the time-evolution equation given by Eq. (5.103) has central symmetry,
given a trajectory of the time evolution of a quantum state, we can flip the
$x$ and $y$ directions and obtain an equally valid trajectory, i.e., employing
the transformation $(x,y)\rightarrow(-x,-y)$ and
$(E_{x},E_{y})\rightarrow(-E_{x},-E_{y})$. Therefore, we use both the data of
the directly simulated state and the data of the symmetric counterpart of the
simulated state produced by the flip of the $x$ and $y$ directions for the AI
to learn. This method is called data augmentation in machine learning
literature because it increases the number of data for the AI to learn, and in
our case the number of data effectively becomes twice. In our numerical
experiments, we randomly flip the $x$ and $y$ directions of the data. The
learning curves with and without the data augmentation technique are shown in
Fig. (5.7), showing that the data augmentation is indeed beneficial for
learning.

Figure 5.7: Learning curves of the cooling problem of a rigid body for the
C-DQN algorithm, with and without the data augmentation technique. The
ordinate and the abscissa are the same in the two figures. The left panel
shows the results including the episodes that end with control failure,
applying Gaussian smoothing with the standard deviation of 40, and the right
panel shows the results excluding the episodes that end with control failure,
applying Gaussian smoothing with the standard deviation of 20. The system
parameters are given by $\dfrac{\gamma}{k}=3$ and $\hat{Q}_{z}=5$.

##### Experimental Results

As discussed above, we follow the same procedure as in Section 5.2 with
several modifications to train the AI using the DQN and the C-DQN algorithms
on the task of cooling a quantum-mechanical rigid body. The learning curves
for different values of system parameters $\gamma$ and $\hat{Q}_{z}$ are shown
in Figs. 5.8, 5.9, 5.10, 5.11 and 5.12, where the performances of the LQG
control are shown for comparison. In these figures, we see that both the DQN
and the C-DQN algorithms can learn the cooling problem and their performances
are comparable, which are also comparable with that of the LQG control,
showing that the AI algorithms can indeed successfully cool the quantum-
mechanical rigid body.

Figure 5.8: Learning curves of the cooling problem of a rigid body for the DQN
and the C-DQN algorithms, with the system parameters $\dfrac{\gamma}{k}=3$ and
$\hat{Q}_{z}=80$, and the horizontal dashed line shows the performance of the
LQG control. The left panel shows the results including the episodes that end
with control failure, applying Gaussian smoothing with the standard deviation
of 40, and the right panel shows the results excluding the episodes that end
with control failure, applying Gaussian smoothing with the standard deviation
of 10. The figure in the right panel is enlarged. Figure 5.9: Learning curves
of the cooling problem of a rigid body for the DQN and the C-DQN algorithms as
in Fig. 5.8, with the system parameters $\dfrac{\gamma}{k}=2$ and
$\hat{Q}_{z}=80$. Figure 5.10: Learning curves of the cooling problem of a
rigid body for the DQN and the C-DQN algorithms as in Fig. 5.8, with the
system parameters $\dfrac{\gamma}{k}=1$ and $\hat{Q}_{z}=80$. Figure 5.11:
Learning curves of the cooling problem of a rigid body for the DQN and the
C-DQN algorithms as in Fig. 5.8, with the system parameters
$\dfrac{\gamma}{k}=3$ and $\hat{Q}_{z}=5$. Figure 5.12: Learning curves of the
cooling problem of a rigid body for the DQN and the C-DQN algorithms as in
Fig. 5.8, with the system parameters $\dfrac{\gamma}{k}=1$ and
$\hat{Q}_{z}=5$.

As shown in the above figures, compared with the DQN algorithm, while the
C-DQN algorithm learns faster at the beginning of training when the behaviour
of the controller is unstable on this task, it learns more slowly at a later
stage. Given an equal amount of time of training, the performances of the
C-DQN algorithm can be marginally lower than that of the DQN algorithm, and as
we gradually reduce the learning rate throughout training, the DQN algorithm
quickly reaches a stable performance while the learning of the C-DQN algorithm
is slowed down and does not exactly reach the performance of the DQN algorithm
at the end of training. To see the results for a prolonged training period, we
double the training time of the C-DQN algorithm for the experiments for the
system parameters $\frac{\gamma}{k}=1,\hat{Q}_{z}=80$ and
$\frac{\gamma}{k}=3,\hat{Q}_{z}=5$, and the results are shown in Figs. 5.13
and 5.14. The results show that the performance of the C-DQN algorithm indeed
approaches that of the DQN algorithm given a longer time of training, and
confirms the effectiveness of the learning of the C-DQN algorithm.

The lower performance of the C-DQN algorithm compared with the DQN algorithm
can be attributed to the behaviour of the C-DQN algorithm on the stochasticity
of the task. As discussed in Section 4.3.2 and experimentally shown in Section
4.4.1, if the time evolution of the system is stochastic, the C-DQN algorithm
does not necessarily converge to the optimal solution of the Bellman equation,
and it may converge somewhere between the solution of the residual gradient
algorithm and the DQN algorithm, and as a result, the performance may be
slightly worse than that of the optimal control, which is a disadvantage of
the C-DQN algorithm. Nevertheless, as discussed in Section 4.4.1 and shown in
Figs. 5.13 and 5.14 above, after a prolonged period of training, the
performance of C-DQN can improve to be approximately equal to that of DQN,
and, therefore, the disadvantage of C-DQN is not significant and can easily be
remedied on this task.

Figure 5.13: Learning curves of the cooling problem of a rigid body for the
DQN and the C-DQN algorithms as in Fig. 5.10, with a doubled time of training
for the C-DQN algorithm. Figure 5.14: Learning curves of the cooling problem
of a rigid body for the DQN and the C-DQN algorithms as in Fig. 5.11, with a
doubled time of training for the C-DQN algorithm.

Concerning properties of the system of the controlled rigid body, comparing
Figs. 5.8, 5.9 and 5.10, we see that with the reduced measurement strength
$\gamma$, the fluctuations in performance are reduced and the system becomes
more stable, and the controllers achieve lower energies of cooling and the
reinforcement learning algorithms learn faster. Comparing Figs. 5.8 and 5.11,
when the angular momentum $\hat{Q}_{z}$ is reduced, with a large measurement
strength, the stochasticity of the system increases, and the rate of the
failure of control becomes significant for the case of $\frac{\gamma}{k}=3$
and $\hat{Q}_{z}=5$. This is because a large angular momentum $\hat{Q}_{z}$
behaves effectively like a magnetic field which localizes the state, and with
a localized state, the measurement backaction of the position measurement is
weaker, which leads to less stochastic perturbation and better stability of
the system. Therefore, the control is easier for the case of
$\frac{\gamma}{k}=3$ and $\hat{Q}_{z}=80$ compared with the control for the
case of $\frac{\gamma}{k}=3$ and $\hat{Q}_{z}=5$. When the measurement
strength is relatively weak, as shown in Figs. 5.10 and 5.12, the change of
$\hat{Q}_{z}$ may not make a significant difference.

When the stochasticity is large and the energy is relatively high, as shown in
Fig. 5.11 and 5.14, the rate of the failure of control is non-negligible, and
the overall performance of the DQN algorithm and that of the C-DQN algorithm
are better than the performance of the LQG control. If we exclude the cases of
control failure, the energy of the control of the C-DQN algorithm is higher
than that of DQN which is then higher than that of the LQG control. This
result implies that although the energy of the system cooled by the LQG
controller is often relatively low, the LQG control fails more frequently than
the DQN controller and the C-DQN controller, and the DQN controller fails more
frequently than the C-DQN controller. This result is consistent with our
argument regarding the nonlinearity of the system. When the energy of the
state is high and the wave function moves away from the center of the trap,
the nonlinearity becomes significant and the LQG control tends to have an
inferior performance. The result also shows that C-DQN tends to avoid failure
and has a more stable behaviour compared with DQN.

###  5.4 Conclusion and Future Perspectives

In this section, we summarize and discuss the results, and we present the
conclusions and future perspectives.

##### Summary and Discussion

In the experiments of cooling of a quantum quartic oscillator, we confirm that
both of the DQN and the C-DQN algorithms can learn to sufficiently cool the
nonlinear quantum system, in which case conventional control strategies do not
work satisfactorily as the state is non-Gaussian in general [52]. However, the
DQN algorithm exhibits significant instability and does not perform stably
throughout the course of training, and it has a large variance in its final
performance, which is likely due to the complexity of the nonlinear quantum
system. On the other hand, the C-DQN algorithm always performs stably for this
complicated system, and as a consequence, it achieves a better performance at
the end of training and has a negligible variance in its final performance.
The results demonstrate the stability and reliability of the C-DQN algorithm
for physical control problems.

In the experiments of cooling of a quantum-mechanical rigid body, we confirm
that both of the DQN and the C-DQN algorithms can learn to stabilize and cool
the system of the trapped quantum-mechanical nanorod, which is a non-Euclidean
multi-dimensional system, and that when the linear approximation holds, the
constrained conventional LQG controller also cools the system efficiently.
Within the regime of the system parameters we have considered, both the DQN
and C-DQN algorithms perform comparably with the LQG controller. The final
performances of the C-DQN algorithm are marginally lower than those of the DQN
algorithm given the same computational budget, because the C-DQN algorithm
generally learns more slowly and possibly converges to a solution that is
suboptimal when the system is stochastic. Nevertheless, our experimental
results show that the performances of the C-DQN algorithm can match those of
the DQN algorithm if the AI is trained for a longer period of time. Compared
with the results of the cooling of a quartic oscillator, we see that the DQN
algorithm may be preferred if it does not encounter any instability, in which
case the problem is often relatively easy and the AI learns quickly and
converges to an optimal solution. If the problem to be solved is difficult and
complicated, the DQN algorithm may suffer from significant instabilities, in
which case the C-DQN algorithm clearly performs better. The C-DQN algorithm
tends to have more stable and reproducible behaviour in the learning process,
which is beneficial for research and fine-tuning in general.

##### Future Perspectives

Concerning future perspectives, as the C-DQN algorithm may not converge to the
optimal solution, it is worthwhile to investigate the case where the C-DQN and
the DQN algorithms are used in combination. For example, one may train the AI
using the C-DQN algorithm at an early stage to avoid instability, and change
to the DQN algorithm to improve the final performance at a later stage. It is
also worthwhile to consider possible modifications of the C-DQN algorithm so
that it can satisfactorily deal with stochasticity.

Regarding the study of a quantum-mechanical rigid body, the model of
measurement we considered is only approximate, and models of more convenient
and realistic measurement protocols are desired. For example, one may consider
the measurement of the position of the head of the nanorod in the three-
dimensional space, by attaching particles or charges that are easily
measurable at the heads of the nanorod. Although different models of
measurement should result in different behaviour of the system, we believe the
our results concerning the DQN and the C-DQN algorithms are universal and hold
true for any measurement model. As experimental techniques continue
developing, the controls that we have investigated may be applied to real
experiments to cool the state of a rigid body to a quantum regime.

For simplicity, the rigid body we considered in our numerical experiments is
axially symmetric. It is also worthwhile to investigate the more complicated
case where the rigid body is asymmetric in general. The nonlinearity of an
asymmetric rigid body is more significant, and as $\hat{Q}_{z}$ does not
commute with the Hamiltonian, one would be able to cool all rotational degrees
of freedom using control fields in the $x$ and $y$ directions only. However,
experimentally, it can be more difficult to measure the orientation of an
asymmetric rigid body, and the computational cost is also considerably higher.

##  Chapter 6 Conclusions

In this thesis, we have reviewed the formulations of continuous quantum
measurement and deep reinforcement learning in Chapters 2 and 3. We have
discussed the non-convergence issue of conventional Q-learning strategies and
the inefficiency issues of existing convergent approaches in the reinforcement
learning literature, and developed a new convergent deep Q-learning algorithm
in Chapter 4, which we call the convergent deep Q network (C-DQN) algorithm,
as an alternative to the conventional deep Q network (DQN) algorithm. The
C-DQN algorithm is provably convergent, scalable and efficient, and we have
demonstrated its effectiveness on standard benchmarks in the reinforcement
learning literature, namely, the Atari 2600 benchmark [50]. Finally, in
Chapter 5, we have applied the C-DQN algorithm to the measurement-feedback
cooling problems of a quantum-mechanical quartic oscillator and a trapped
quantum-mechanical rigid body. We presented the physical models and analysed
the properties of the systems, and showed that although both of the DQN and
the C-DQN algorithms can learn to cool the systems, the C-DQN algorithm learns
stably and has better performances if the DQN algorithm suffers from
instability when the task is difficult; however, the C-DQN algorithm learns
relatively more slowly when the task is sufficiently simple such that the DQN
algorithm can work stably and quickly. Because the performances of the DQN
algorithm can have large variances and lack consistency from trial to trial if
the underlying task is difficult, the C-DQN algorithm can be a better choice
for researches on complicated physical control problems.

Our contribution is twofold: we have investigated the non-convergence issue of
the standard reinforcement learning algorithm, Q-learning, and developed a new
convergent algorithm and examined the properties of our algorithm; we have
established the quantum-mechanical model of the trapped and controlled rigid
body, and demonstrated the effectiveness of our control strategies for the
measurement-feedback cooling problem of this system.

Regarding future directions, we may consider the combination of the DQN and
the C-DQN algorithms so that we can obtain both the stability of the C-DQN
algorithm and the high performance of the final result of the DQN algorithm.
It is also desired if the C-DQN algorithm can be improved to deal with
stochasticity satisfactorily and to converge to an optimal solution of the
Bellman equation in the presence of stochasticity. Concerning the study of a
quantum rigid body, the control strategies we have investigated may be applied
to real experiments, using application-specific integrated circuits which
embed the control strategies to control the lasers to reduce the energy of the
trapped rigid bodies. The control strategies we have considered can help cool
the system of a trapped rigid body so that a quantum regime may be realized,
which has applications in sensing devices and fundamental physical research
[47]. It is also possible to extend our research to the case of a more
complicated asymmetric rigid body, which has highly nonlinear dynamics.

This thesis contributes to the field of the interdisciplinary study of quantum
control and machine learning, and we hope that our work helps the development
of the use of machine learning technologies for physical problems and the
development of better control strategies in the microscopic quantum world.

## Bibliography

  * [1] Daoyi Dong and Ian R Petersen.  Quantum control theory and applications: a survey.  IET Control Theory & Applications, 4(12):2651–2671, 2010.
  * [2] Warren S. Warren, Herschel Rabitz, and Mohammed Dahleh.  Coherent control of quantum dynamics: The dream is alive.  Science, 259(5101):1581–1589, 1993.
  * [3] Steven Chu.  Cold atoms and quantum control.  Nature, 416(6877):206, 2002.
  * [4] Moshe Shapiro and Paul Brumer.  Quantum control of molecular processes.  John Wiley & Sons, 2012.
  * [5] Marcos Dantus and Vadim V Lozovoy.  Experimental coherent laser control of physicochemical processes.  Chemical reviews, 104(4):1813–1860, 2004.
  * [6] Vlasta Bonacić-Kouteckỳ and Roland Mitrić.  Theoretical exploration of ultrafast dynamics in atomic clusters: Analysis and control.  Chemical reviews, 105(1):11–66, 2005.
  * [7] Peter Lodahl, Sahand Mahmoodian, and Søren Stobbe.  Interfacing single photons and single quantum dots with photonic nanostructures.  Rev. Mod. Phys., 87:347–400, May 2015.
  * [8] Ronald Hanson, Leo P Kouwenhoven, Jason R Petta, Seigo Tarucha, and Lieven MK Vandersypen.  Spins in few-electron quantum dots.  Reviews of modern physics, 79(4):1217, 2007.
  * [9] Iulia M Georgescu, Sahel Ashhab, and Franco Nori.  Quantum simulation.  Reviews of Modern Physics, 86(1):153, 2014.
  * [10] H. Häffner, C.F. Roos, and R. Blatt.  Quantum computing with trapped ions.  Physics Reports, 469(4):155 – 203, 2008.
  * [11] Lucio Robledo, Lilian Childress, Hannes Bernien, Bas Hensen, Paul FA Alkemade, and Ronald Hanson.  High-fidelity projective read-out of a solid-state spin quantum register.  Nature, 477(7366):574, 2011.
  * [12] Andreas Wallraff, David I Schuster, Alexandre Blais, Luigi Frunzio, R-S Huang, Johannes Majer, Sameer Kumar, Steven M Girvin, and Robert J Schoelkopf.  Strong coupling of a single photon to a superconducting qubit using circuit quantum electrodynamics.  Nature, 431(7005):162, 2004.
  * [13] Shuntaro Takeda and Akira Furusawa.  Perspective: Toward large-scale fault-tolerant universal photonic quantum computing.  arXiv preprint arXiv:1904.07390, 2019.
  * [14] Markus Aspelmeyer, Tobias J. Kippenberg, and Florian Marquardt.  Cavity optomechanics.  Rev. Mod. Phys., 86:1391–1452, Dec 2014.
  * [15] Brian D. O. Anderson and John B. Moore.  Optimal Control: Linear Quadratic Methods.  Prentice-Hall, Inc., Upper Saddle River, NJ, USA, 1990.
  * [16] J Werschnik and E K U Gross.  Quantum optimal control theory.  Journal of Physics B: Atomic, Molecular and Optical Physics, 40(18):R175–R211, sep 2007.
  * [17] Patrick Doria, Tommaso Calarco, and Simone Montangero.  Optimal control technique for many-body quantum dynamics.  Phys. Rev. Lett., 106:190501, May 2011.
  * [18] Navin Khaneja, Timo Reiss, Cindie Kehlet, Thomas Schulte-Herbrüggen, and Steffen J. Glaser.  Optimal control of coupled spin dynamics: design of nmr pulse sequences by gradient ascent algorithms.  Journal of Magnetic Resonance, 172(2):296 – 305, 2005.
  * [19] Jens Jakob WH Sørensen, Mads Kock Pedersen, Michael Munch, Pinja Haikka, Jesper Halkjær Jensen, Tilo Planke, Morten Ginnerup Andreasen, Miroslav Gajdacz, Klaus Mølmer, Andreas Lieberoth, et al.  Exploring the quantum speed limit with computer games.  Nature, 532(7598):210–213, 2016.
  * [20] Ehsan Zahedinejad, Sophie Schirmer, and Barry C. Sanders.  Evolutionary algorithms for hard quantum control.  Phys. Rev. A, 90:032310, Sep 2014.
  * [21] Navin Khaneja, Roger Brockett, and Steffen J. Glaser.  Time optimal control in spin systems.  Phys. Rev. A, 63:032308, Feb 2001.
  * [22] A. C. Doherty and K. Jacobs.  Feedback control of quantum systems using continuous state estimation.  Phys. Rev. A, 60:2700–2711, Oct 1999.
  * [23] Leigh Martin, Felix Motzoi, Hanhan Li, Mohan Sarovar, and K. Birgitta Whaley.  Deterministic generation of remote entanglement with active quantum feedback.  Phys. Rev. A, 92:062321, Dec 2015.
  * [24] F. Motzoi, J. M. Gambetta, P. Rebentrost, and F. K. Wilhelm.  Simple pulses for elimination of leakage in weakly nonlinear qubits.  Phys. Rev. Lett., 103:110501, Sep 2009.
  * [25] David Guéry-Odelin, Andreas Ruschhaupt, Anthony Kiely, Erik Torrontegui, Sofia Martínez-Garaot, and Juan Gonzalo Muga.  Shortcuts to adiabaticity: concepts, methods, and applications.  Reviews of Modern Physics, 91(4):045001, 2019.
  * [26] S. Meiboom and D. Gill.  Modified spin-echo method for measuring nuclear relaxation times.  Review of Scientific Instruments, 29(8):688–691, 1958.
  * [27] Lorenza Viola, Emanuel Knill, and Seth Lloyd.  Dynamical decoupling of open quantum systems.  Phys. Rev. Lett., 82:2417–2421, Mar 1999.
  * [28] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E. Hinton.  Imagenet classification with deep convolutional neural networks.  In Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 1, NIPS’12, pages 1097–1105, USA, 2012\. Curran Associates Inc.
  * [29] Florian Schroff, Dmitry Kalenichenko, and James Philbin.  Facenet: A unified embedding for face recognition and clustering.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 815–823, 2015.
  * [30] Ewan Dunbar, Xuan Nga Cao, Juan Benjumea, Julien Karadayi, Mathieu Bernard, Laurent Besacier, Xavier Anguera, and Emmanuel Dupoux.  The Zero Resource Speech Challenge 2017.  arXiv e-prints, page arXiv:1712.04313, Dec 2017.
  * [31] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [32] Christopher Berner, Greg Brockman, Brooke Chan, Vicki Cheung, Przemysław Dębiak, Christy Dennison, David Farhi, Quirin Fischer, Shariq Hashme, Chris Hesse, et al.  Dota 2 with large scale deep reinforcement learning.  arXiv preprint arXiv:1912.06680, 2019.
  * [33] Oriol Vinyals, Timo Ewalds, Sergey Bartunov, Petko Georgiev, Alexander Sasha Vezhnevets, Michelle Yeo, Alireza Makhzani, Heinrich Küttler, John Agapiou, Julian Schrittwieser, John Quan, Stephen Gaffney, Stig Petersen, Karen Simonyan, Tom Schaul, Hado van Hasselt, David Silver, Timothy Lillicrap, Kevin Calderone, Paul Keet, Anthony Brunasso, David Lawrence, Anders Ekermo, Jacob Repp, and Rodney Tsing.  StarCraft II: A New Challenge for Reinforcement Learning.  arXiv e-prints, page arXiv:1708.04782, Aug 2017.
  * [34] David Silver, Thomas Hubert, Julian Schrittwieser, Ioannis Antonoglou, Matthew Lai, Arthur Guez, Marc Lanctot, Laurent Sifre, Dharshan Kumaran, Thore Graepel, Timothy Lillicrap, Karen Simonyan, and Demis Hassabis.  A general reinforcement learning algorithm that masters chess, shogi, and go through self-play.  Science, 362(6419):1140–1144, 2018.
  * [35] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10684–10695, 2022.
  * [36] Yuki Fujimoto, Kenji Fukushima, and Koichi Murase.  Methodology study of machine learning for the neutron star equation of state.  Phys. Rev. D, 98:023019, Jul 2018.
  * [37] Tomohiro Mano and Tomi Ohtsuki.  Phase diagrams of three-dimensional anderson and quantum percolation models using deep three-dimensional convolutional neural network.  Journal of the Physical Society of Japan, 86(11):113704, 2017.
  * [38] E. Barberio, B. Le, E. Richter-Was, Z. Was, J. Zaremba, and D. Zanzi.  Deep learning approach to the higgs boson $cp$ measurement in $H\rightarrow\tau\tau$ decay and associated systematics.  Phys. Rev. D, 96:073002, Oct 2017.
  * [39] Talitha Weiss and Oriol Romero-Isart.  Quantum motional state tomography with non-quadratic potentials and neural networks.  arXiv preprint arXiv:1906.08133, 2019.
  * [40] Yue Liu, Tianlu Zhao, Wangwei Ju, and Siqi Shi.  Materials discovery and design using machine learning.  Journal of Materiomics, 3(3):159 – 177, 2017.  High-throughput Experimental and Modeling Research toward Advanced Batteries.
  * [41] Murphy Yuezhen Niu, Sergio Boixo, Vadim N. Smelyanskiy, and Hartmut Neven.  Universal quantum control through deep reinforcement learning.  npj Quantum Information, 5(1):33, 2019.
  * [42] Marin Bukov, Alexandre G. R. Day, Dries Sels, Phillip Weinberg, Anatoli Polkovnikov, and Pankaj Mehta.  Reinforcement learning in different phases of quantum control.  Phys. Rev. X, 8:031086, Sep 2018.
  * [43] Dalgaard Mogens, Motzoi Felix, Jens Jakob Sørensen, and Sherson Jacob.  Global optimization of quantum dynamics with alphazero deep exploration.  NPJ Quantum Information, 6(1), 2020.
  * [44] Thomas Fösel, Petru Tighineanu, Talitha Weiss, and Florian Marquardt.  Reinforcement learning with neural networks for quantum feedback.  Phys. Rev. X, 8:031084, Sep 2018.