Multi-phonon relaxation and generation of quantum states in a nonlinear mechanical oscillator

The dissipative quantum dynamics of an anharmonic oscillator is investigated theoretically in the context of carbon-based nano-mechanical systems. In the short-time limit, it is known that macroscopic superposition states appear for such oscillators. In the long-time limit, single and multi-phonon dissipation lead to decoherence of the non-classical states. However, at zero temperature, as a result of two-phonon losses the quantum oscillator eventually evolves into a non-classical steady state. The relaxation of this state due to thermal excitations and one-phonon losses is numerically and analytically studied. The possibility of verifying the occurrence of the non-classical state is investigated and signatures of the quantum features arising in a ring-down setup are presented. The feasibility of the verification scheme is discussed in the context of quantum nano-mechanical systems.


Introduction
Nanoelectro-mechanical (NEM) resonators such as cantilevers, membranes and beams are important systems for the study of quantum phenomena in macroscopic, mechanical man-made objects [1]. For about a decade, the goal has been to cool the NEM resonators to the ground state and to verify the achievement of this state. This has recently been accomplished [2][3][4], and the next step is the generation, detection and coherent manipulation of more intricate quantum states, e.g., superpositions of macroscopic states (Schrödinger cat states) [5] and Fock states [6].
While a harmonic quantum oscillator displays a behavior analogous to its classical counterpart, the same does not hold for nonlinear oscillators. The presence of a conservative (Kerr) nonlinearity facilitates the creation of non-classical states. Although top-down fabricated NEM-resonators typically possess a finite Kerr constant µ, it is usually too small to alone suffice for non-classical state generation. In contrast, carbon nanotubes (CNTs) and graphene membranes can exhibit strong enough nonlinearities for non-classical states to emerge during free evolution [5,7]. Using CNTs or graphene also alleviates the need for active cooling due to their low masses and high frequencies [1,5]. For this reason, carbon based NEM-resonators are highly promising for the study of nonlinear mechanical systems in the quantum regime.
If a system possesses conservative nonlinearities, it is natural to also take into account higher orders in the system-environment coupling. This results in nonlinear damping (NLD) [8][9][10][11]. In carbon based NEM-resonators it has recently been shown that the dissipative mechanisms can indeed be nonlinear in the oscillator coordinate [12,13]. Hence, to fully understand the quantum dynamics of carbon based NEM resonators it is essential to comprehend the interplay between conservative nonlinearities and NLD. In general, NLD will create interesting quantum features rather than remove them [14][15][16][17][18]. For example, as is well known in quantum optics [14][15][16]19], for a harmonic oscillator, the steady state ρ obtained from NLD due to two-photon losses at zero temperature (T = 0) is a mixed state with Tr[ρ 2 ] < 1. Moreover, because of the parity conserving nature of the NLD, this steady state density matrix has coherences in terms of nonzero off-diagonal elements ρ 01 in the number basis. Recently, it was proposed to use this effect in order to generate and preserve Schrödinger cat states in a double well potential realized by a SQUID ring [18]. In an optomechanical setup non-thermal and squeezed states were shown to be obtainable by two-phonon cooling [20].
In this paper we study a situation more relevant to the coherent quantum dynamics of carbon based NEM-resonators, namely, a nonlinear oscillator with NLD. In this field there is little to be found in the literature [10]. Our main goals are to characterize the resulting steady states along with the corresponding features of the relaxation dynamics. For the nonlinear oscillator with NLD we find that the parity conserving nature of NLD helps to preserve phase coherence also in the presence of a finite conservative nonlinearity. However, we also find that for sufficiently large ratio between the Kerr constant µ and the NLD strength γ 2 , the coherence of the steady state is suppressed.
Further, at finite temperatures, with NLD, the final state is not a standard thermal state with a reduced density matrix ρ ∝ e −H/k B T unless linear damping (LD) is present. The corresponding Wigner distribution has a dip in the center, whose depth depends on the initial conditions. In the presence of both NLD and LD, the relaxation is governed by a crossover between two different relaxation rates. In the long-time limit, the decay rate Γ is independent of the initial oscillator amplitude for NLD. Instead it is a function only of temperature T . To facilitate experimental determination of the governing damping mechanism, we also present the corresponding signatures in the time evolution of the position coordinate variance.
This paper is organized as follows: In section 2, we present the underlying theoretical frameworks used for analytical (rotating wave approximation [RWA]) and numerical (Wangsness-Bloch-Redfield [21][22][23]) calculations. The latter are used to verify the validity of the RWA-calculations. Then, in section 3, we present results pertaining to the steady states for various scenarios involving NLD and LD. Results on the relaxation towards the steady states are found in section 4, followed by a section with summary and conclusions.

Model
To model the dynamics of the anharmonic oscillator subject to linear and nonlinear damping, we consider a single mechanical mode with frequency Ω and Kerr constant µ. The mode is coupled linearly and nonlinearly to a reservoir of harmonic oscillators with mode frequencies ω k . The mode mass and the reduced Planck constant are set to be m = 1 and = 1, respectively. Correspondingly, the total Hamiltonian is where the mode position and momentum operators are q = q 0 (a † + a)/ √ 2, p = i(a † − a)/ √ 2q 0 , and similarly the bath operators are given by Here q 0 = 1/Ω and x 0k = 1/m k ω k are the zero point amplitudes of the oscillator and the bath, and satisfy the canonical commutation relations. The linear and the nonlinear coupling constants of the k'th reservoir mode to the resonator mode are denoted by η 1k and η 2k , respectively. These couplings lead to amplitude-independent (linear) and amplitude-dependent (non-linear) [8][9][10][11] damping of the oscillator motion. The former coupling involves the exchange of single vibration quanta (vibrons) with the reservoir, while in the latter case two vibrons are transferred simultaneously.
We consider the following scenario: The oscillator mode is assumed to be initially in the ground state |0 . Then, by applying a voltage pulse for a finite time, the oscillator is displaced. The corresponding state is, by definition [24], a coherent state, Ψ(t = 0) = D(α) |0 ≡ |α . Here, D(α) = exp(αa † − α * a) is the displacement operator and the amplitude α depends on the strength and duration of the voltage pulse. During dissipation-free evolution under H S only, a coherent state will evolve into a Yurke-Stoler cat state (|α +i|−α )/ √ 2. This state will periodically re-appear with a frequency set by the Kerr constant µ. Typically the state's emergence is slow compared to the system's eigenfrequency Ω [25]. When taking into account the interaction with the reservoir, the created cat states will be influenced by the dissipation and the oscillator will eventually relax into a steady state. The state of the mechanical mode during the relaxation is described by the reduced density matrix (RDM) ρ = tr B ρ T , which is obtained from the state of the total system ρ T by tracing over bath degrees of freedom.
The time-evolution of the reduced density matrix is described by a (generalized) quantum master equation (QME) [26]. Here we follow the standard approach using the Born-Markov approximation to obtain the QME from the von Neumann equation for the total density matrix [26]. In this approximation the QME in the interaction picture with respect to H S is given by ∂ ∂t The operators A m are given by the expansion of the interaction Hamiltonian in the interaction picture where Assuming the reservoir to be initially in thermal equilibrium, the reservoir correlation functions are given by where n B (ω) = [e ω/k B T − 1] −1 is the Bose distribution function, ρ B is the thermal state and κ ml (ω) = κ lm (ω) = 2π k (η mk η lk q m+l 0 )δ(ω − ω k ) is the spectral density of the environment coupling.
In general, the frequency dependence of κ ml is determined by the specific geometry of the system. As long as the spectral density is smooth and slowly varying around ω = Ω and ω = 2Ω the detailed dependence on ω is not important [27] (see also 2.2). To be specific, we consider an Ohmic spectral density [26], i.e., we take The values of Γ 11 and Γ 22 are used to specify the strength of LD and NLD, respectively. The off-diagonal value Γ 12 can, in general, be chosen independently of Γ 11 and Γ 22 . If the oscillator couples to two independent baths, the off-diagonal spectral density vanishes. In the RWA, which is discussed in 2.2, the dissipation terms associated with Γ 12 can be neglected. Here we set Γ 12 = √ Γ 11 Γ 22 , which vanishes if either linear or nonlinear damping is absent.

Relaxation dynamics in the rotating frame
Typically, for NEMS the time scales associated with Ω and µ are well separated with Ω µ (see also section 5), which justifies an analytic treatment in the RWA. Within this approximation the QME (2) in the Schrödinger representation becomes [28] where the Lindblad superoperators are given by The superoperators L 0 , L 1 and L 2 describe the coherent evolution due to H S , linear and the nonlinear interaction with the reservoir, respectively. The linear and nonlinear dissipation rates, γ ± and γ 2± , are Fourier transforms of the correlation function (5), and are given by We would like to note, that the RWA equations given above are only valid for a weak nonlinearity [28]. This is reflected by the dependence of the rates on the transition frequency of the harmonic part of the Hamiltonian H S and not on the frequencies of the full nonlinear H S . The validity of this approximation in the present case will be shown in the next sections by comparing to numerical calculations which include the exact transition frequencies (see also section 2.3). The solution to (7) is obtained by rewriting the matrix equation as a set of decoupled equations. To this end, ρ is projected onto the number basis ρ n,n+λ = n|ρ|n + λ , where λ = 0, 1, 2, . . .. It is convenient to define ρ n,n+λ (t) ≡ n! (n+λ)! ψ n (λ, t)e iλ(Ω+µ)t [14], which also removes oscillations due to the term linear in a † a in (7). For fixed λ the function ψ n (λ, t) gives the entries of the λ'th super diagonal of ρ. Hence, (7) becomes As indicated above, the equations for different λ do not couple. In particular, for λ = 0, we obtain a single equation for the populations ρ n,n (t) = ψ n (0, t), which is used in the following to analyze the relaxation to the steady state.
In general, when all coefficients in (10) are non-zero the time evolution of ψ n (λ, t) has to be found numerically. If either γ − or γ 2− is vanishing, i.e., there are only one or two vibron losses present, the time-dependence of ψ n (λ, t) can be calculated analytically (see [29][30][31] and [14]). In these cases it is also straight-forward to find the steady state solutions for all λ. An exact solution can also be obtained in terms of the positive P function where the nonlinear damping is equivalent to the real part of a complex Kerr constant [32]. At finite temperatures, when in addition to the losses, one or two vibron excitations (γ ± > 0 or γ 2± > 0) contribute, the stationary (time independent) probability distribution ρ n,n can be found using the detailed balance condition [33,34] (see Appendix B). If one and two-vibron processes compete (γ ± > 0 and γ 2± > 0), a solution for ρ n,n can be given in terms of a continued fraction [35] or (for γ 2+ = 0) in terms of confluent hypergeometric functions [36].
In the next sections we are in particular interested in the behavior of the position variance where • = tr S [•ρ]. Expressing q in terms of creation and annihilation operators in the rotating frame one finds where all oscillating contributions have been discarded. The expectation values of a † a and a † can be expressed in terms of ψ n (0, t) and ψ n (1, t), respectively, Knowing ψ n (0, t) and ψ n (1, t) therefore allows for an evaluation of the position variance.
In the following we focus on obtaining such solutions in the long-time limit.

Numerical Computation
In order to investigate the oscillator dynamics and to verify the results obtained by the RWA approximation, the QME (2) is solved numerically by Wangsness-Bloch-Redfield [21][22][23] calculations. To this end the operators A m are represented in the eigenbasis of the oscillator Hamiltonian H S and the Hilbert space is truncated after 30 states (40 states for α = 3), which yields converged results. The initial coherent state, ρ(0) = |α α|, is generated by applying the displacement operator D(α) to the ground state.

Zero temperature
We begin with the RWA-analysis of the steady state ρ(∞) at zero temperature. To evaluate the applicability of the RWA, we compare the results to numerical calculations as described in section 2.3. If only single vibron losses are present in the system (γ + = 0, γ 2± = 0), one finds from the equation of motion (10) for ρ that any initial state ρ(0) will decay to the ground state, i.e., ρ(∞) = |0 0|. The associated position variance is ∆q 2 /q 2 0 = 1/2. On the other hand, if only nonlinear losses are present (γ 2+ = 0, γ ± = 0), a general initial state will relax to the steady state [14,15] where the matrix elements in (15) depend on the initial density matrix ρ(0). Since for two-vibron losses parity is conserved, the only nonzero matrix elements on the main diagonal of ρ(∞) are The steady state off-diagonal elements ρ 0,1 and ρ 1,0 are in general non-zero and can also be found from a conservation law. From (10) one gets where the coefficients C 2n are given in terms of Gamma functions (see Appendix A) This expression is a generalization of the result of Simaan and Loudon [14] to the case of an anharmonic oscillator (µ = 0). Note, that unless µ = 0, the off-diagonal element ρ 0,1 is time-dependent in the long time limit. Specifically, for a system which is initially in a coherent state |ψ(t = 0) = |α , i.e., with matrix elements the constants of motion are Here, 0 F 1 is the confluent hypergeometric limit function [37]. For µ = 0 the latter becomes a Bessel function of the first kind and the expression given in [14] is recovered.
To illustrate the differences between the numerically calculated states obtained by LD and NLD, their Wigner distributions [38] 15)). In general it can be shown that the purity of this state is Tr[ρ 2 ] = 1 − 2P even P odd + 2|ρ 0,1 | 2 < 1. Figures 1 (c) and (d) display typical snapshots of this steady state for two different values of the ratio between the nonlinear damping rate and the Kerr constant. By comparing the sizes of the stationary states obtained by LD and NLD, it is seen that the latter occupy a larger phase-space area. This suggests, that the position variance ∆q 2 is suitable for quantifying the differences between the states resulting from these damping mechanisms.
Using the definition (12) and the results above, one can evaluate the variance to be figure 2, the variance according to (21) is shown as function of the displacement amplitude α of initial coherent states for different values of γ 2− /µ (dashed lines). The numerical results (symbols) are in accordance with (21). For comparison, also the variance of the ground state is included (red line), indicating its independence on the initial displacement amplitude. For small initial displacements, (α 1), the behavior of the variance is governed by the coherences ρ 0,1 = ρ * 1,0 , since ρ 1,1 → 0. For large amplitudes, (α 1), and weak NLD, (γ 2− µ), the variance saturates at 1 (gray dotted line). This is due to the ground and the excited states being equally populated, ρ 0,0 = ρ 1,1 = 1/2, and the vanishing off-diagonal elements, since the hypergeometric function converges to zero. In the limit of strong NLD, (γ 2− µ), and large amplitudes, (α 1), the hypergeometric function converges to a non-zero value and the variance becomes ∆q 2 → 1 − 1/2π, which is larger than the variance of the ground state.

Finite temperature
Also at finite temperature ρ(∞) obtained by NLD will differ from the state obtained by LD. For LD only, ρ(∞) is fully characterized by the values of the diagonal matrix 3.5 elements of the density matrix, where Here, k B is the Boltzmann constant and T is the temperature of the bath. This result is found from (10) by taking λ = 0 and setting Figure 2. Zero temperature position variance of the steady state resulting from NLD as a function of the displacement amplitude α of the initial coherent state. The variance is averaged over one period 2π/µ. Dashed lines show the behavior according to (21) for different values of the ratio γ 2− /µ. Symbols denote the corresponding numerical results. The lower dotted line (red) indicates the variance of the ground state, resulting from LD. The upper dotted line (gray) indicates the upper variance limit, resulting from weak NLD for α → ∞.
the time-derivative to zero. One obtains the detailed balance condition, which leads to (22). As we have pointed out in section 2.2, the RWA equations (7) are valid for weak nonlinearities since the bath is only probed at the frequency Ω. Consequently, at finite temperatures the stationary probability distribution calculated within the RWA will not depend on the Kerr constant. However, for sufficiently small temperatures, i.e., µ a † a Ω + µ, the approximation is valid as we will also show in the following. Since for nonlinear damping the two-vibron exchange conserves parity, one finds that detailed balance holds for each parity sector individually in this case. It can be shown (Appendix B) that Note that the stationary probability distribution at finite temperature also depends on the initial state. This implies that expectation values will in general also depend on P even and P odd . For example, the average number of vibrons in the oscillator mode is given by which is derived in Appendix B. Typical features of a stationary state obtained from nonlinear damping at finite temperature are displayed in figure 3. Figure 3 (a) shows the simulated Wigner distribution (see (20)), with parameters γ 2− = µ, α = 1, k B T /Ω = 1. A sector is cut out for better visualisation of its features. As discussed in section 3.2, the thermal probability distribution of a nonlinearly damped state is split into an odd and an even parity sector. Each of the sectors has a Boltzmann-like distribution (see (23)). The shape of the Wigner distribution reflects this parity dependence by having a dip centered at the origin. It can be shown that W (0, 0) = (P even − P odd )/2π. Therefore, the dip in the Wigner distribution only depends on the initial state and not on temperature. Hence, in the case of an initial coherent state |α , the dip will only depend on the initial oscillator displacement. For large displacement amplitudes the difference between the sectors becomes smaller, and W (0, 0) → 0. If the initial state had an odd parity, the dip of the Wigner distribution would attain a negative value. Figure 3 (b) shows the Boltzmann-like distributions of the odd and the even sectors of the stationary state discussed above. The narrow bars represent the distribution of the thermal state obtained by linear damping (see (22)). The insets show the Wigner distributions of the even and odd sectors.
As the Wigner distribution of the finite temperature stationary state displays quantum signatures, so does its variance. Figure 4 shows the stationary state variance as function of the displacement amplitude α of the initial coherent state for several  temperatures (k B T /Ω ∈ [0, 1/4, 1/2, 3/4, 1] and γ 2− = µ). Dashed lines correspond to the analytical expression given by (33). Numerical calculations (symbols) are shown for the same parameters. The variance is seen to always be larger than the zero temperature variance, displayed in figure 2. This is due to an additional contribution from the Boltzmann distribution and the vanishing off-diagonal elements of the density matrix (see section 4.3). The variance is still dependent on the initial displacement and saturates at ∆q 2 /q 2 0 = 1 + 2n B > 1 for large α, which is in contrast to k B T /Ω = 0 for which ∆q 2 /q 2 0 ≤ 1.

Zero temperature, NLD
We now investigate the relaxation dynamics of initial coherent states which evolve into the steady states discussed in the previous section. Figure 5 (a) shows the numerically calculated variance ∆q 2 as function of rescaled time τ = µt/π for several values of the ratio γ 2− /µ. The variance is sampled at multiples of the period Ω/(2π), and therefore fast oscillations are not seen. For weak NLD (γ 2− /µ 10), the short time dynamics (τ 1) is governed by the Kerr constant µ. This is reflected in the initial increase of variance. The maximum variance, at which the Yurke-Stoler cat state is first expected to appear [5,25], is attained at τ ≈ 1/2. The variance of this state is given by ∆q 2 cat = q 2 0 (4α 2 + 1)/2, which is always larger than the variance of a coherent state for a nonzero α. As seen in the figure 5 (a), the shape of the variance peaks changes in time as the cat states decohere by NLD. Figure 5 (b) shows the Wigner distribution sampled at τ = 1/2, displaying a cat state under influence of NLD. Some of its original features are still visible, like the remnants of the negative interference domains. After τ = 1, the system gradually relaxes to the state (15). In figure 5 (a) this is reflected by steady oscillations around a mean value of ∆q 2 which is, as discussed in section 3.1, larger than 1/2. The oscillations in ∆q 2 are caused by the complex valued, oscillating coherences of ρ 0,1 . The relaxation dynamics described above becomes faster with increasing ratio γ 2− /µ. Figure 5 (c) shows the Wigner distribution of the state sampled at τ = 11/2, where the resemblance to the steady state in figure 1 (c) is clearly seen.

Zero temperature, LD and NLD
As shown in the previous section, NLD at T = 0 leads to the formation of the nonclassical state (15). In the additional presence of LD the generation of this state will be affected and eventually the ground state will be reached. In the following we discuss the resulting decay of the variance, which can be used to extract information about the nature of the damping mechanism and the state initially created by NLD. We introduce the ratio γ − /γ 2− as a measure of the relative dissipation strength, and set γ 2− = µ.
It is seen that only ρ 0,0 is non-zero in the stationary limit, i.e., ρ 1,1 and ψ n decay for sufficiently large times. In the limit of strong NLD, γ − γ 2− , one finds that the populations ρ n,n for n > 1 rapidly decay due to two-vibron losses, but P 1 is only decreased via one-vibron losses. Hence,ρ 1,1 ≈ −γ − ρ 1,1 and ρ 1,1 (t) ≈ e −γ − t P odd . Here it is assumed that the initial decay on the time-scale 1/γ 2− is not influenced by the LD and ρ 0,0 (0) ≈ P even and ρ 1,1 (0) ≈ P odd . Similarly, one finds that ψ 0 (1, t) ≈ exp[(iµ − γ − /2)t]ψ even (0). Therefore, the position variance will decay to the minimum uncertainty value 1/2 according to Figure 6 (a) shows the time evolution of the numerically calculated ∆q 2 for different values of the ratio γ − /γ 2− and µ = γ 2− . For γ − /γ 2− 1/10, the evolution of the variance for short times τ 1 resembles the result shown in figure 5. The shorttime dynamics is again influenced by the Kerr constant, while the long-time behavior, (τ 1) is showing the relaxation of the state created by NLD. In this stage the variance is exponentially decaying with the decay being faster for larger values of γ − /γ 2− . In figure 6 (b) the time evolution of (1 − ρ 0,0 ) is shown, confirming the decay to the ground state. The parameters are the same as in 6 (a). For equally strong LD and NLD, (γ − = γ 2− ), an exponential decay is seen for the entire time interval. For LD weaker than NLD, (γ − < γ 2− ) and τ 1/2, the decay is initially not exponential, whereas for τ 1 the ground state population displays exponential increase. Figure 6 (a) suggests the concept of sequential relaxation for weak LD. The first stage of the sequence (0 < τ 1), is mainly dominated by NLD and the Kerr constant, which leads to the formation of the state in (15). The second stage (τ 1), is mainly influenced by the relaxation to the ground state because of LD. In other words, the state in (15) is initially created by NLD, and then decays due to LD. Accordingly, the average variance in the second relaxation stage is expected to be given by (27), which is displayed in figure 6 (a) (dashed lines). It can be seen that the idea of sequential relaxation is quantitatively correct for γ − γ 2− . Furthermore, (27) implies that the decay rate, which is equal to γ − , is independent of the initial amplitude α. This is shown in figure 6 (c) where the time evolution of ∆q 2 (τ ) is plotted for α = 1, 3/2 and 3 (γ 2− = µ, γ − /γ 2− = 1/10).

Finite temperature, NLD
While at T = 0 two-vibron losses facilitate the creation of the non-classical state (15), thermal two-vibron excitations will contribute to its decoherence. From (10) one sees that a finite rate γ 2+ couples the states (n = 0, 1) to other excited states. These thermal excitations will eventually destroy the coherence of the non-classical states and the off-diagonal matrix elements will vanish. To quantify this decay we consider low temperatures, such that γ 2+ γ 2− . To obtain the relaxation rate to lowest order in γ 2+ it is sufficient to study excitations from n = 0 to n = 2. Consequently, one gets from (10) the coupled equationṡ supplemented by initial conditions for ψ 0 (1, 0) and ψ 2 (1, 0). For small γ 2+ we can assume that these initial values correspond to the values of the state in (15), which is the steady state solution for γ 2+ = 0. This amounts to setting ψ 0 (1, t = 0) = ψ even (0), ψ 2 (1, t = 0) = 0 .
The system (28) can be solved for ψ 0 (1, t), (30) describes the relaxation of a non-classical state due to weak thermal excitations. The relaxation rate, which is determined by χ, depends on the strength of the Kerr constant µ and temperature via γ 2± . It is interesting to consider the limiting cases of strong and weak NLD. For strong NLD, γ 2− /µ → ∞, one finds χ → γ 2− + (11γ 2+ )/2 and ψ 0 (1, t) becomes Similarly, for weak NLD, γ 2− /µ → 0, χ → −iµ + 4γ 2+ + γ 2− and the time evolution is Obviously, the decoherence rate increases with increasing strength of the Kerr constant. This can be understood from the fact that the accumulated phase shift in the excited state is proportional to µ. Therefore, for finite µ there is an additional dephasing, which contributes to the decay of ψ 0 (1, t). Altogether, the position variance becomes Figure 7 (a) shows the time dependence of the numerically calculated variance for the same parameters as in figure 3 (d) with the initial displacement amplitude α = 1. The T = 0 variance is included as a reference (blue) and is almost identical to the curve for k B T /Ω = 1/4 (green). When comparing the graphs in figure 7 (a) with figure 5 (a) and 6 (a), one finds a qualitatively similar behavior. It is seen that also for finite temperatures the relaxation proceeds sequentially. The short-time dynamics show an increase in variance, corresponding to the formation of a Schrödinger cat state. For longer times the variance is increasing with time, and the rate of approaching the stationary value is faster for larger temperatures. In order to quantify the finite temperature relaxation rate, the time evolution of | ∆q 2 (τ ) − ∆q 2 (∞) | T >0 is shown in Figure 7 (b), where the stationary state variance ∆q 2 (∞) T >0 is subtracted from ∆q 2 (τ ) T >0 in order to remove the finite asymptotic value. The numerically obtained variance (symbols) is averaged over one period t = 2π/µ, which removes the oscillations due to the finite Kerr constant µ. The temperatures are the same as in figure 7 (a). The figure confirms the idea of different relaxation stages, with the decay in the second stage being exponential and the rate being temperature dependent.
For comparison, figure 7 (c) shows the time evolution of |ρ 0,1 (τ )| 2 for k B T /Ω = 1/2 and α ∈ [1, 3/2, 3]. It can be seen that in the long-time limit the decay rates do not depend on the initial oscillator displacement α. As discussed above, the off-diagonal  (35). (c) Time dependence of |ρ 0,1 (τ )| 2 at k B T /Ω = 1/2 for α = 1 (•), α = 3/2 ( ) and α = 3 ( ). The long-time limit relaxation rate is according to (34). matrix elements ρ 0,1 decay due to thermal excitations. For low temperatures, their evolution is given by (30). When the Kerr constant and the NLD are equally strong, (µ = γ 2− ), one finds χ → γ 2− + (19γ 2+ )/4 − i(µ − 3γ 2+ /4) and the time evolution is given by The resulting behavior is shown in figure 7 (c) and it can be seen that (34) correctly describes the decay of ρ 0,1 . From figures 7 (a) and (b) it can be concluded that the variance increase is mainly due to the decay of |ρ 0,1 | 2 . Consequently, we make the which will be valid in the second stage of the relaxation. From (34) we find Γ(T ) ≈ 5γ 2+ (2Ω). This Ansatz is displayed in figure 7 (b) (dashed lines) and a good quantitative agreement with the numerical results is found for all temperatures. Using (35) one can therefore infer the value of |ρ 0,1 | 2 at the end of the first relaxation stage, which is dominated by the NLD. This allows for a reconstruction of the properties of the nonclassical state (15), which is created at the end of the first relaxation stage.

Finite temperature, LD and NLD
Finally, in this section we consider the time evolution of the position variance at finite temperature under influence of both LD and NLD. Figure 8 (a) shows the timedependence of the numerically obtained ∆q 2 for γ 2− = µ, γ − /γ 2− = 1/10, α = 1 and different temperatures. The behavior for T = 0, which corresponds to setting γ − = 0, is included for reference (blue). As in the previous section we subtract the stationary state variance from ∆q 2 to study the relaxation towards the stationary state. This is shown in figure 8 (b). For long times an exponential decay is seen for temperatures k B T /Ω < 1. To quantify the decay we use the Ansatz in (35) and consider where σ and Γ are estimated by a least-square fit to the numerical data. To this end the variance is averaged over one period t = 2π/µ which removes the oscillations due to the finite Kerr constant µ. The fitting for T ≤ 3/4 is done for the time interval 3 < τ < 11 and for T = 1 in the time interval of 2 < τ < 4. The resulting behavior of (36) is displayed as dashed lines in figure 8 (b). In figures 8 (c) and (d) the extracted parameters σ and Γ are shown as functions of initial displacement amplitude α and temperature, respectively. Due to the sequential relaxation, the parameter σ approximately corresponds to the variance of the nonclassical state (15), which is created in the initial stage. The analytic expression for T = 0, given by equation (21), is shown in figure 8 (c) for comparison. The extracted low temperature results are in good agreement with the analytic curve. For higher temperatures deviations from the analytic behavior can be seen. However, the qualitative dependence on the initial amplitude is still similar to the one given by (21).
In figure 8 (d) the temperature dependence of the extracted decay rates Γ is compared to the rates of one-vibron loss γ − and two-vibron excitations γ 2+ , which determine the decay rate in the situations investigated in sections 3.2 and 4.2. It can be seen, that also in the present case the decay rate is independent of the initial coherent state displacement amplitude. For low temperatures, k B T /Ω < 0.3, one-vibron loss is the dominating damping mechanism. For higher temperatures the interplay of onevibron and two-vibron processes leads to a more complicated behavior. At first the rate follows the behavior given by the combined rates γ − + 5γ 2+ /2 but starts to deviate at k B T /Ω ≈ 0.3. With increasing temperature it stays between γ − + γ 2+ and γ − + 5γ 2+ /2. The results presented in this section show that a ring-down measurement of the variance can be used to identify the quantum features of an initially created non-classical state, provided the sequential relaxation regime can be accessed. Moreover, the analysis of the decay rates can provide information about the dominating dissipation mechanisms in the system under investigation.

Summary and conclusions
In this work we have investigated a nonlinear oscillator in the quantum regime. The system is subject to nonlinear damping (two-vibron exchange), which resembles the twophoton decay known from quantum optics. The initially displaced oscillator relaxes to a steady state, which is determined by the relative strength of the linear and the nonlinear damping mechanisms and the temperature of the thermal bath. For sufficiently strong NLD we find that the relaxation is sequential: In the first stage, the initial density matrix decays to a non-classical state (15), followed by the second stage, where this state further decays to the final steady state. The non-classical state depends on the initial state, while the decay rates depend on temperature and the strength of LD and NLD.
At zero temperature, the sequential relaxation allows for a reconstruction of the properties of the non-classical state even in the presence of LD, i.e., when the stationary state is always the ground state. This can be achieved by performing a ring-down measurement of the position variance. For finite temperatures, NLD and thermal excitations lead to a different steady state -a stationary thermal state with a parity dependent Boltzmann-like distribution (23). The parity dependence is a signature of the NLD at finite temperatures. In the sequential relaxation regime the properties of the non-classical state can again be extracted from variance measurements.
In the following we briefly discuss the possibility of observing the relaxation of the non-classical state (15) in experiments with carbon-based nanomechanical resonators. A recent experiment [12] indicates significant NLD in micrometer sized, carbon-based resonators (nanotubes and graphene) with resonance frequencies of several hundreds of MHz. These systems are nonlinear and can be described by the classical Duffing equation of motion [39] where q is the oscillator position, m and Ω are the oscillator mass and frequency, α 0 is the nonlinearity strength, γ is the LD rate and η is the NLD rate. By comparing (37) with the corresponding quantum mechanical equation, the relation between the parameters (µ, γ − , γ 2− ) and the parameters of the Duffing equation can be established, The temperature for reaching the quantum regime without cooling is of the order of T 10 mK for resonators with frequencies in the 100 MHz range. The values obtained in [12] correspond to µ/Ω 10 −3 , meaning that the RWA-analysis works well. The magnitude of the nonlinear damping rate is found to yield γ 2− ≈ µ, which is promising for the formation of a non-classical state. The condition to find the sequential relaxation regime amounts to where q 0 is the zero point amplitude. Consequently, for systems with a large zero point amplitude it is easier to reach the sequential relaxation regime. In this regard, carbon-based resonators are very promising due to their low mass and tunable frequency.
In present experiments the condition γ − < γ 2− is not yet fulfilled, but the ongoing experimental and technological progress in the field will help to improve the parameters in order to meet the condition (39).
In conclusion, our study shows that exploiting NLD present in (carbon-based) NEMS resonators is a promising route to generation and investigation of non-classical effects in NEMS. Studies of the sequential relaxation in multi-partite systems, such as coupled oscillators, would be of interest with respect to the generation of quantum entanglement.

Acknowledgments
We thank J. M. Kinaret and L. Y. Gorelik for helpful discussion. The research leading to this article has received funding from the Swedish Research Council.