---
abstract: 'We apply adaptive feedback for the refrigeration of a mechanical resonator, i.e. with the aim of simultaneously cooling the classical thermal motion of more than one vibrational degree of freedom. The feedback is obtained from a neural network trained via a reinforcement learning strategy to choose the correct sequence of actions from a finite set in order to reduce the total energy of all modes of vibration. The actions are realized either as optical modulations of spring constants or as radiation pressure induced momentum kicks. As a proof of principle we numerically show simultaneous cooling of four independent modes with an overall strong reduction of the total system temperature.'
author:
- Christian Sommer
- Muhammad Asjad
- Claudiu Genes
bibliography:
- 'bibfileReinforcement.bib'
title: Prospects of reinforcement learning for simultaneous damping of many mechanical modes
---

The radiation pressure effect of light onto the motion of mechanical resonators has been extensively employed to bring such macroscopic systems towards the quantum ground state [@Aspelmeyer2014cavity; @windey2019cavity; @deli2019cavity; @Rossi2017; @Clark2017; @Qiu2019; @Asenbaum2013; @Mancini1998; @Clemens2016; @Kiesel2013; @Millen2015]. In a standard approach, the aim is to isolate a single vibrational mode and bring it to its quantum ground state where the only relevant motion is given by the zero-point fluctuations. Cold-damping is one of the used techniques, where one detects motionally-induced phase changes in the cavity output and an electronic feedback loop is implemented to dynamically modify the cavity drive such as to produce an extra optical damping effect [@Genes2008Ground; @Steixner2005; @Bushev2006; @Rossi2018; @Cohadon1999; @Poggio2007; @Wilson2015; @Tebbenjohanns2018Cold]. Alternatively, in the good cavity limit where the photon loss rate is smaller than the mechanical frequency, the resolved sideband technique can be implemented by detuning the drive to the cooling sideband [@Gigan2006self; @Braginsky2007parametric; @Marquardt2007Quantum; @Wilson2007Theory; @Teufel2011Sideband]. As the effect stems from the inherent time delay between the action of the mechanical resonator onto the cavity field and the back-action of light, this can be seen as a sort of automatic cavity induced feedback. While both techniques are successfully applied to isolate and cool a single vibrational mode, they are not necessarily optimal to induce full refrigeration of the mechanical resonator, i.e. to simultaneously cool all the vibrational modes into which the thermal energy is distributed. The main impediment is that the detected output signal only gives information on generalized collective quadrature but not on all modes. This leads to efficient cooling of some collective mode (for example center of mass) while some collective modes become dark and remain in a high temperature state.\
Here, we propose a machine learning approach towards devising a strategy capable of providing refrigeration of the classical motion of a mechanical resonator based on the feedback obtained from the detection of a single optical mode. To this end we provide a proof-of-principle multi-mode numerical simulation using a neural network trained using a reinforcement learning algorithm to generate the feedback signal capable of simultaneously extracting thermal energy from four distinct modes of a single mechanical resonator.\
Machine learning techniques have been recently applied to various applications in quantum physics ranging from the identification of phases in many-body systems, predicting ground-state energies for electrostatic potentials, active learning approaches to propose and optimize experimental setup configurations and towards applications for quantum control and quantum-error correction [@Chen2014Fidelity; @Carrasquilla2017Machine; @Nieuwenburg2017Learning; @Dunjko2017Machine; @Carleo2017Solving; @Mills2017Deep; @Melnikov2018Active; @Bukov2018Reinforce; @Foesel2018Reinforce]. In particular, a few studies [@Chen2014Fidelity; @Bukov2018Reinforce; @Foesel2018Reinforce] successfully applied the technique of reinforcement learning with neural networks [@Russel2018Modern]. This approach originates from the idea, to let an intelligent agent that observes its environment choose an action, that is determined by a given policy trying to optimize a particular reward and minimize a punishment.\
Here, we employ such a technique for optically assisted cooling of the classical thermal state of a multi-mode mechanical resonator system [@Nielsen2017Multimode; @Piergentili2018; @PhysRevA.99.023851]. The learning technique allows one to acquire a nonlinear function that chooses a feedback action that will be applied on the dynamical system upon taking the full or partial measured state of the system as an input. The training of this function that is given by a dense neural network is obtain by trial and error and quantified by an increased reward that is obtained by successfully reducing the energy of the resonators.\

![image](fig1a.pdf){width="2.0\columnwidth"}

The physical systems considered are depicted in Fig. \[fig1\]. The mechanical resonator is subject to environmental noise described by a standard Brownian motion stochastic force leading to thermalization at some equilibrium temperature $T$. The feedback action is implemented via the radiation pressure force, i.e. photon kicks either from one or two sides. The induced damping is straightforward in the two-sided kicking case \[illustrated in Fig. \[fig1\]a\]: the read-out of motion is followed by appropriate kicking action from the side towards which the resonator is moving. However, one-sided kicking \[illustrated in Fig. \[fig1\]b\] already suffices as the oscillator can be displaced by a constant force and damped around the modified equilibrium position. The typical weak free space photon-phonon interaction can be drastically increased by the filtering of the action field through a high-finesse optical cavity as shown in Fig. \[fig1\]c. Such a situation is characterized by a linear coupling of the photon number to the membrane’s displacement and has been extensively studied in single mode cooling via cavity time delayed effects [@Genes2008Ground] or by implementation of cold damping techniques [@Genes2008Ground] especially in the bad cavity regime. The membrane-in-the-middle [@Thompson_2008; @Jayich_2008; @Asjad2014Robust] scenario in Fig. \[fig1\]d,e corresponds to a quadratic coupling in displacement leading to the possibility of optically modulating the mechanical oscillation frequency [@Asjad2014Robust]. We describe in Fig. \[fig1\]e a possible approach for feedback cooling via cavity field detection and neural network assisted feedback.\
We will focus in the following on the bad cavity case where the cavity back-action is negligible and where standard cold damping techniques are used for cooling. In such a case, the situations described in Fig. \[fig1\]b and Fig. \[fig1\]c are physically equivalent with the difference that in Fig. \[fig1\]c the action of a single photon is multiplied by a large number roughly proportional to the finesse of the cavity. We will therefore distinguish physically distinct situations such as parametric cooling and linear cooling. First we analyze the performance of a neural network suggested set of actions onto the cooling of a single mode via parametric modulation of the oscillation frequency: we describe the shape of the action and numerically show the efficient reduction of energy from the initial thermal distribution. We then apply the technique to the linear cooling of four distinct modes of the resonator and find a more complex set of actions required for efficient simultaneous cooling of all four modes (with limitations arising due to the numerical complexity of the simulations).

**Model and equations** — We consider a membrane resonator with a few modes of oscillations of frequencies $\omega_j$ (where $j=1,...N$). We start with a quantum formulation of the system’s dynamics aimed at future treatments of cooling in the presence of quantum noise. However the current formulation aims only at the reduction of classical thermal noise and is therefore obtained as a classical averaging of the quantum equations of motion. The Hamiltonian for the collection of modes is written as $H_m =\sum_{j=1} \hbar\omega_j/2\left(p_{j}^{2} + q_{j}^{2} \right),$ in terms of dimensionless position and momentum quadratures $q_{j}$ and $p_{j}$ for each independent membrane oscillation mode. The effect of the thermal reservoir can be easily included in a set of equations of motion supplemented with the proper input stochastic noise terms:

\[FreeEq.1\] $$\begin{aligned}
\dot{q}_j &=\omega_{j} p_j \\
\dot{p}_j &= -\omega_{j}q_j - \gamma_{j} p_j + \xi_j + F_{j}(t).\end{aligned}$$

The parameter $\gamma_{j}$ describes the damping of the $j$’s resonator mode. Its associated zero-averaged Gaussian stochastic noise term leading to thermalization with the environment can be fully described by the two-time correlation function: $$\begin{aligned}
\label{FreeEq.2}
\langle \xi_j(t) \xi_{j'}(t')\rangle &=& \frac{\gamma_{m_j}}{\omega_j}\int_{0}^{\Omega} \frac{d\omega}{2\pi}e^{-i\omega(t-t')} S_{\text{th}}(\omega)\delta_{jj'},\end{aligned}$$ where $\Omega$ is the frequency cutoff of the reservoir and the thermal noise spectrum is given by $S_{\text{th}}(\omega)=\omega[\coth\left(\hbar \omega/2k_B T\right) + 1]$. For sufficiently high temperatures $k_B T \gg \hbar \omega_j$, the correlation function becomes a standard white noise input with delta correlations both in frequency and time. Specifically, one can approximate $\langle \xi_j(t) \xi_{j'}(t')\rangle \approx (2\bar{n}_{j}+1)\gamma_{m_j}\delta(t-t')\delta_{jj'}$, where the occupancy of each vibrational mode is given by $\bar{n}_{j} = (\exp(\hbar \omega_j/k_B T)-1)^{-1} \approx k_B T/\hbar \omega_j$. For numerical simulations we generate a stochastic input noise as a delta-correlated Wiener increment with variance proportional to the integration time-step (see Appendix) and follow an approach described in Ref. [@Higham2001An]. For consistency we check (in the Appendix) that the thermal bath indeed correctly describes the expected thermalization of an initially cold oscillator towards the equilibrium temperature $T$ at a rate given by $\gamma$.\
The momentum kicks selected by the network are encompassed in the action of the force terms $F_{j}(t)$. This can be realized for example by the radiation pressure effect of the laser beam, modulated by a device like an AOM (acousto-optic modulator). Here, forces acting on different resonators given by $F_j$ and $F_{j'}$ for $j \neq j'$ differ only by a constant multiplication factor as they are all obtained from the same quantity (the output field).\
To amplify the effect of the action force onto the mechanical resonator one can utilize optical cavities. A cavity also allows control over the coupling by placing the membrane either in a node (quadratic coupling) or anti-node (linear coupling) of the cavity mode. The Hamiltonian is now modified by the addition of the free cavity mode $\hbar \omega_c a^{\dagger} a$, laser driving resonant to the cavity mode $i\hbar \mathcal{E}(t)\left(a^{\dagger} -a \right)$ (in a frame rotating at $\omega_c$) and optomechanical interaction of linear $\sum_j \hbar g^{(1)}_j a^{\dagger} a q_j$ or quadratic form $\sum_j \hbar g^{(2)}_j a^{\dagger} a q^2_j$. The amplification effect of the light field amplitude can be seen from the relation $\mathcal{E}(t) = \sqrt{2\mathcal{P}_{0}(t)\kappa_L /\hbar \omega_{c}}$ connecting the driving amplitude to the input laser power $\mathcal{P}_{0}(t)$ through the left mirror with losses at rate $\kappa_L$. For high-finesse cavities photons perform many round trips before leaking out through the mirrors resulting in a large momentum transfer onto the mirror: this can be seen by taking the limit of small $\kappa_L$ resulting in a large value of $\langle a(t) \rangle$ for a given $\mathcal{P}_{0}(t)$. Notice that we considered a double-sided cavity with left $\kappa_L$ and right $\kappa_R$ decay rates adding to the total loss rate $\kappa = \kappa_L + \kappa_R$. The coefficients $g^{(1)}_j$ and $g^{(2)}_j$ are the linear and quadratic per photon optomechanical coupling rates corresponding to the two situations depicted in Fig. \[fig1\]c and Fig. \[fig1\]d. While the cavity field amplitude inherently depends on the displacement of the mechanical mode, we will assume the unresolved sideband regime where this dependence is weak. Moreover, we are interested in the classical problem i.e. in simulating the proper set of actions that results in the shrinking of an initial large thermal distribution for the total energy of the oscillator. To this end we only consider the trivial dynamics of the cavity field classical amplitude $\alpha(t) = \langle a(t)\rangle$ which follows the driving field as $\dot{\alpha}(t) = -\kappa \alpha + \mathcal{E}(t)$. We can then reduce the dynamics of the system to

$$\begin{aligned}
\label{OscEq.2a}
\dot{q}_j &= \omega_j p_j \\
\label{OscEq.3a}
\dot{p}_j &= -\omega_jq_j - \gamma_{j} p_j + \xi_j - g^{(1)}_j|\alpha(t)|^2\\
\label{OscEq.4a}
\dot{\alpha} &= -\kappa \alpha + \mathcal{E}(t),\end{aligned}$$

which resemble Eqs. \[FreeEq.1\], where we can identify the action forces $F_j(t)=-g^{(1)}_j|\alpha(t)|^2$ (the cavity field $\alpha(t)$ playing the role of the action delivering the cooling momentum kicks to the mechanical oscillators). As noted before, as the actions are obtained from the same cavity field intensity, they only differ by the multiplicative $g^{(1)}_j$ factor. Notice also that this configuration strongly resembles a cold damping approach [@Genes2008Ground; @Steixner2005; @Bushev2006].\
In contrast, for a quadratic coupling Hamiltonian, the changes in the momentum are of a very different nature $$\begin{aligned}
\label{eq4}
\dot{p}_j &=& -\left[\omega_j + 2g^{(2)}_j|\alpha(t)|^2\right]q_j - \gamma_j p_j + \xi_j,\end{aligned}$$ as the cavity periodically modulates the oscillation frequencies of each mode.\
To provide the neural network feedback onto the motional dynamics, we use the inferred $q$ and $\dot{q}$ at a given time $t-\Delta t$ as input values for the neural network \[see Fig. \[fig1\]\]. The trained network then selects the appropriate action by choosing the value of $\mathcal{E}(t)$ (from a finite number of possible values) to be acted upon the system. The size of the time-step $\Delta t$ is chosen such that $1/\omega_{j} \gg \Delta t$ for all $j$ to minimize the error in the numerical integration. For a given drive amplitude, the set of actions on the different modes will be different according to the values of the optomechanical couplings (as they are proportional to $g^{(1)}_j|\alpha(t)|^2$ or $2g^{(2)}_j|\alpha(t)|^2$).We then use the Runge-Kutta fourth-order method (RK4) for the numerical integration of the dynamical system where we iteratively sum for each time step. Additionally, at each time step we inject the measured parameters of the dynamical system as input data into the nonlinear function formed by the neural network to predict on the action and thereby the momentum kick or cavity field strength suitable for the next time step, which is acquired from the output nodes (neurons) of the network.\
**Reinforcement learning**—The neural network provides a nonlinear function, that for some given input data, which harbor information about the oscillator states at a given time step $t$, predicts the correct action for the next time step $t+\Delta t$ that helps to reduce the overall energy of the dynamical system at later times. For the training of the neural network we employ the technique of reinforcement learning [@Sutton1998Reinforcement; @Russel2018Modern] and in particular a policy gradient approach [@Williams1992Simple]. Here, the network acts as an agent that by observing parameters of the environment (resonator) improves its probabilistic policy that chooses the right actions $a_t$ at a given time $t$ to increase an overall reward $R$ that is connected to the reduction of the energy of the resonator modes. The actions are chosen from a finite set (of values of different amplitudes) and realized as momentum kicks or translated into frequency shifts. As an input to the network we feed information about the state of the environment given by $s_t = (q(t) ,\dot{q}(t))$. The network outputs the probabilities $\pi_{\theta}(a_t|s_t)$ for the actions $a_t$ that could be applied to the dynamical system. Here, the parameter $\theta$ encompasses all the weights and biases of the network. We take the action with the highest probability and apply it in the next iteration of Eq. \[FreeEq.1\] up to Eq. \[OscEq.3a\]. The probabilities $\pi_{\theta}(a_t|s_t)$ can be optimized with respect to an increased reward return $R_t$ by employing an update rule for the weights and biases of the neural network, following $\theta \rightarrow \theta + \Delta \theta$ and $$\begin{aligned}
\label{ReinEq.1}
\nonumber
\Delta \theta_{j} &=& \eta \partial_{\theta_{j}} \mathbb{E}[R] = \eta\mathbb{E}\left[\sum_{t}(R_t-b)\partial_{\theta_{j}}(\ln \pi_{\theta}(a_t|s_t))\right],\\\end{aligned}$$ where $\mathbb{E}$ is the expectation value over all state and action sequences (full trajectories), which here is approximated by averaging over a large enough set of oscillator trajectories (training batch) and their corresponding action sequences which we have obtained from the iterative summation of the dynamical equations (RK4) and from the predictions of the neural network at each time step for various randomly chosen initial conditions. The learning rate is given by the parameter $\eta$ and $b$ is a baseline to suppress fluctuations of the reward gradient [@Williams1992Simple; @Weaver2001Optimal]. Here, the baseline is approximated by $b \approx b_n = (1/n)\sum_{i=1}^{n}R^{(n)}$, where $R^{(n)}$ is the average total Reward from the $n$’s learning epoch. Here, the training epoch is defined as the number of updates $\theta \rightarrow \theta + \Delta \theta$.\
The neural network which is represented by the array $\theta$ encompassing all weights and biases, consists of an input and output layer and two hidden layers whereby the number of input neurons depends on the number measured variables of the system while the number of output neurons depends on the number of possible output actions, respectively (see Appendix). The two hidden layers consists of up to $60$ to $100$ neurons each. The network is densely connected and we chose “relu” (rectified linear unit) as a nonlinear function acting on each neuron in the two hidden layers. The probabilities for each action given out by the output layer are obtained by using the “softmax” nonlinear function for the output neurons. From these probabilities the action is chosen by taking the neuron index with the highest probability value in the output.\

![*Single mode parametric cooling* (a) Time dynamics followed on two trajectories with initial conditions drawn from a Boltzmann distribution corresponding to an initial thermal state with average occupancy of $\bar{n} = 100$. The inset refers to the choice of cooling performed parametrically by modulation of the spring constant. (b) Corresponding phase space trajectories. (c) Action sequences chosen by the network for each trajectory (d) Average energy for a thermal ensemble of trajectories exposed to the actions chosen by the network and rescaled to give the value of the occupation number of the harmonic oscillator. (e) Zoom-in into the time dynamics of the actions and quadratures for two distinct trajectories. (f) Phase space comparison of initial (orange) and final (blue) distributions. (g) Corresponding histogram of energy distribution in the initial and final states. The parameters are $\gamma = 4\times 10^{-5}\, \omega$, $g^{(2)} = 1\times 10^{-8}\,\omega$, $|\alpha|^{2} \approx 0.5\times 10^{7}$ and $\Delta t = 0.05\,\omega^{-1}$.[]{data-label="fig2"}](fig2.pdf){width="0.97\columnwidth"}

**Single mode cooling** — In a first step we numerically simulate the time dynamics of a single oscillating mode of frequency $\omega$ initially in a thermal distribution imposed by its coupling to an environment at some temperature $T$. This corresponds to the following distribution of energies $$\begin{aligned}
\label{SinEq.1}
P(E) &=& Z^{-1}e^{-\beta E \hbar \omega},\end{aligned}$$ with a partition function $Z \approx (\beta \hbar \omega)^{-1}$ and the total normalized energy $E = (q(0)^2 + p(0)^2)/2$. We then randomly pick an initial energy value from the thermal distribution $$\begin{aligned}
\label{SinEq.2}
E &=& -(1/\beta \hbar \omega)\ln(1-s),\end{aligned}$$ by picking a random number $s$ between zero and one. From the equipartition theorem we deduce $q(0)$ and $p(0)$ and train the neural network by recursively injecting sequences of $(q(t) ,\dot{q}(t))|_{[0\dots T]}$, obtained by applying the terms of Eqs. \[OscEq.2a\],\[OscEq.4a\] and Eq. \[eq4\] recursively on the initial values, as training data into the network. A reward at each time step is only given when the action reduced the energy of the resonator at a given time step with respect to the previous time. The reward is defined by $$\begin{aligned}
\label{SinEq.3}
R_t = (E_{0} - E_{t})\theta(E_{t-1}-E_{t}),\end{aligned}$$ where $E_t$ is the total energy at time $t$. The reward gets larger when the energy separation between the current and initial energy increases therefore optimizing the effective cooling rate. The application of the reward to the single mode cooling is exemplified for the quadratic coupling configuration illustrated in Fig. \[fig1\]d. The cooling dynamics is exemplified both as amplitude decreases Fig. \[fig2\]a and in phase space Fig. \[fig2\]b on two trajectories corresponding to two different initial states randomly picked from a thermal initial distribution with average occupancy $\bar{n}=100$. The action sequences of the network in Fig. \[fig1\]c show a periodic structure matching the same frequency of the resonator more visible in the zoom-in plot provided in Fig. \[fig2\]e. There one can follow the time dynamics of the applied action and the effect onto both the position and momentum. The results of thousands of sample trajectories are presented as initial and final phase space distributions in Fig. \[fig2\]f and as a histogram of the energy distributions in Fig. \[fig2\]g. The injected thermal noise in all plots in Fig. \[fig2\] corresponds to a thermal occupation number of $\bar{n} \approx 100$ and a thermalization rate of $\gamma/\omega=4\times 10^{-5}$.\

![*Simultaneous cooling of 4 modes*. (a) Time dynamics of the oscillation amplitudes of four independent modes under the action of a collective force in the setup shown in the inset. (b) Corresponding time dynamics in phase space. (c) Neural network indicated sequence of actions leading to simultaneous cooling. (d) Magnification of the action and the first mode as a reference for the time scale. (e) Corresponding decrease of the total energy of all four modes towards very low temperatures. (f) Reduction of the initial thermal distribution (orange) towards a low temperature state (blue) for each independent oscillation mode after $2\times 10^4$ time steps. (g) Histogram of initial (orange) and final (blue) energy distribution for $4\times 10^3$ trajectories. The parameters for the simulation are given by $\omega_{2,3,4} = (0.8,\, 1.2,\, 0.6)\,\omega_1$, $\gamma_{1,2,3,4} = (4,\,3,\, 5,\, 2)\,\times 10^{-5} \omega_1$, $g_{1,2,3,4} = (0.3,\,0.2,\,0.4,\,0.3)$ and $\Delta t = 0.05\,\omega_1^{-1}$.[]{data-label="fig3"}](fig3b.pdf){width="0.95\columnwidth"}

**Simultaneous cooling of many modes** — We display the generality of this approach by applying the network to find a strategy to cool up to four modes simultaneously. As in the case of the single mode cooling, we apply a single force on the mirror: this poses a challenge as a good cooling strategy for a given mode might actually lead to the heating of the other modes. In general, owing to this challenge, a simultaneous cooling strategy has an increased complexity in the choice of the action sequences, which leads to overall slower cooling rates. The same principles for cooling a single resonator are applied to cooling four modes subjected to the same actions by the network as presented in Fig. \[fig3\]. Here, we use the setup configuration presented in Fig. \[fig1\]b. While the four modes have different frequencies $\omega_j$ and coupling strengths $g_j$ they are subjected to the same time sequence of actions delivered by intensity variations of an impinging laser beam in free space or via the field intensity $|\alpha(t)|^2$ in the cavity.

The input to the network is given by $s_t = (Q(t), \dot{Q}(t))$ with $Q(t) = \sum_{j} q_{j}(t)$ and $\dot{Q}(t) = \sum_{j} \omega_jp_{j}(t)$ as a collective position coordinate and its derivative. The quantities can be obtained for example from an interferometer that is sensitive to the fluctuations on the membrane (see Fig. \[fig1\]e). The derivation of the initial condition is described in the Appendix.\
Here, the agent needs to find a strategy that simultaneously cools the center of mass motion as well as all of the relative mode dynamics. As an example the trajectories of the four modes are presented in Fig. \[fig3\]a,b, where the cooling results from the corresponding actions presented in Fig. \[fig3\]c with magnified view shown in Fig. \[fig3\]d. In contrast to the action sequence imposed on the trajectory of a single resonator presented in Fig. \[fig2\]c, which basically shows a periodic signal matching the frequency of the oscillator, here we find a more complex signal with a quasi periodic pattern. The change of the average energy of the four resonators as a function of time is presented in Fig. \[fig3\]f. In Fig. \[fig3\]e the initial values obtained from a Boltzmann distribution and final phase space values of a thousand trajectories for all oscillators are presented. A histogram of the sum of their individual energies is given in Fig. \[fig3\]g. These results show that all four resonators can be simultaneously cooled down to lower temperatures, that differ by orders of magnitude from their initial values and thereby exemplify the strength of this adaptive approach.\
**Conclusions -** We have shown results of numerical simulations for the simultaneous cooling of a few degrees of freedom of a vibrating mechanical resonator. The feedback action has been realized via the reinforcement learning technique implemented on a neural network. Simultaneous cooling of a few modes indicates the possibility of full refrigeration of mechanical resonators via optical control. The technique could be easily extended to cool a number of oscillators or a number of particles trapped inside optical cavities or with tweezers. While the present treatment considers classical stochastic dynamics, a full quantum theory of neural network aided cooling will be tackled in the future that might also lead towards feedback production of squeezed or squashed states. In this regard, for both single and many oscillation modes, we plan to analyze the efficiency of neural network cooling in comparison with standard cold-damping optical cooling. It is expected that a Fourier analysis of the action function indicated by the network could hint towards feedback implementations that could surpass existing techniques, especially at the level of many degrees of freedom.\
**Acknowledgments -** We acknowledge financial support from the Max Planck Society. We acknowledge fruitful discussions with Michael Reitz. We are thankful for the comprehensive lecture notes on machine learning for physicists presented by Florian Marquardt at the University of Erlangen-Nuremberg [@Marquardt].

Appendices {#appendices .unnumbered}
==========

Initial conditions
------------------

We obtain the initial conditions for the oscillation modes from a Boltzmann distribution of energies $P(E_1, \dots, E_n) = Z^{-1}e^{-\beta(E_1+\dots+ E_n)}$ where $Z = 1/\beta^n$. By integration $$\begin{aligned}
\label{Ap:Eq.1}
\nonumber
\int^{\tilde{E}_1}_{0}\dots \int^{\tilde{E}_n}_{0} dE_{1}\dots dE_{n}P(E_1, \dots, E_n)= \prod _{j=1}^{n} \left[1 - e^{-\beta\tilde{E}_{j}} \right] =s,\end{aligned}$$ where $s$ is a value between $0$ and $1$. We set $b_j = \left[1 - e^{-\beta\tilde{E}_{j}} \right]$ and define that $b_j \in [0,1]$ for $0 \leq j \leq n$ which results in $$\begin{aligned}
\label{Ap:Eq2}
\tilde{E}_{j} &=& -\frac{1}{\beta}\ln(1-b_j).\end{aligned}$$ Since $p_{j}(0)^{2} + q_{j}(0)^{2} = 2\tilde{E}_{j}/\omega_j$ we obtain $p_j$ and $q_j$ from $$\begin{aligned}
\label{Ap:Eq3}
q_{j}(0) = \sqrt{\frac{2\tilde{E}_{j}}{\omega_j}}\cos(2\pi\phi_{j}) \qquad \text{and} \qquad p_{j}(0) = \sqrt{\frac{2\tilde{E}_{j}}{\omega_j}}\sin(2\pi\phi_{j}),\end{aligned}$$ where $\phi_{j}$ is a random number between $0$ and $1$ for all $j\in {1,\dots, n}$.

Thermalization dynamics
-----------------------

![*Thermalization dynamics of a single mode of vibration*. (a) Phase space coordinates from initialization with a Boltzmann distribution with $\bar{n} = 100$ (black points) and obtained from the evolution equation containing thermal noise. An initial distribution around the phase space origin (blue points) is evolved up to the steady state (red points). In (b) the corresponding energy distributions following the same color coding are presented. (c) Two example trajectories with starting energies close to zero are following the evolution guided by thermal noise. In (d) the average energy is presented as a function of time.[]{data-label="figA"}](Appfig1.pdf){width="0.6\columnwidth"}

Let us describe the numerical procedure for simulating the action of a thermal environment onto the state of the mechanical resonator. We consider the equations of motion $$\begin{aligned}
\label{Ap:Eq4}
\dot{q} = \omega p \qquad \text{and} \qquad \dot{p} = -\omega q - \gamma p + \xi,\end{aligned}$$ which we rewrite as equations of differential forms

\[Ap:Eq5\] $$\begin{aligned}
dq &= \omega dp \\
\label{Ap:Eq5a}
dp &= -\omega dq - \gamma dp + \sqrt{(2\bar{n}+1)\gamma}dW(t),\end{aligned}$$

where the noise $dW(t)$ is included as a Wiener process. Numerically, this can be realized by $\Delta W(t) \approx \sqrt{\Delta t}N(0,1)$, where $N(0,1)$ describes a normally distributed random variable of unit variance - consequently the Wiener increment is normally distributed with a variance equal to the numerical time increment $\Delta t$. As a numerical check we simulate the thermalization of an initially cold mechanical mode under the action of an environment with rate $\gamma=4 \times 10^{-5}\, \omega$ and at an effective occupancy $\bar{n}=100$. The analytical Boltzmann distribution for this occupancy is shown in Fig. \[figA\]a as black dots while the final state obtained from the numerical integration is represented by the red dots (blue dots in the middle are the initial state). In Fig. \[figA\]b the agreement between the numerical simulation and the Boltzmann distribution is illustrated as a histogram of energy states. The time evolution is shown in Fig. \[figA\]c as dynamics for the position and in Fig. \[figA\]d for the total energy showing it approaching $\bar{n}=100$ in the long time limit.

Network parameters
------------------

The simulations were run on a standard Laptop computer (CPU, Intel(R) Core(TM) i$7-5500$U CPU @ $2.40\,$GHz). We use the Keras package for Python and the Theano framework [@AlRfou2016Theano] to realize the neural network and the reinforcement learning procedure. The network and training parameters are presented in Tab. \[Ap:Tb1\].

                                                                Network parameters
  --------------------------------------- ----------------------------------------
  Training batch size                                                           80
  Learning rate $\eta$ (Adam optimizer)                0.00008 (single oscillator)
                                                         0.0006 (four oscillators)
  Neurons per layer                          $(2, 60, 60, 11)$ (single oscillator)
                                            $(2, 100, 100, 11)$ (four oscillators)
  Reward scale                                                                 1.0

  : *List of neural network and training parameters.*[]{data-label="Ap:Tb1"}
