Single-molecule topological actuator for nano-mechanical modes

We investigate theoretically the dynamics of two quasi-degenerate orthogonal mechanical modes of a suspended nanowire coupled to the two-level system of a single-fluorescent molecule by Stark effect. We show that by driving the molecular two-level system with a laser field one can engineer the effective mechanical spectrum leading to an exceptional degeneracy point where the two mechanical modes coalesce. It allows the topological actuation of the modes by adiabatically varying the laser frequency and intensity. This is supported by quantum Monte-Carlo simulations to account for strong quantum fluctuations induced by spontaneous emission. We finally propose a frequency-modulation detection scheme based on single-molecule spectroscopy.

We investigate theoretically the dynamics of two quasi-degenerate orthogonal mechanical modes of a suspended nanowire coupled to the two-level system of a single-fluorescent molecule by Stark effect. We show that by driving the molecular two-level system with a laser field one can engineer the effective mechanical spectrum leading to an exceptional degeneracy point where the two mechanical modes coalesce. It allows the topological actuation of the modes by adiabatically varying the laser frequency and intensity. This is supported by quantum Monte-Carlo simulations to account for strong quantum fluctuations induced by spontaneous emission. We finally propose a frequency-modulation detection scheme based on single-molecule spectroscopy.

√
Hz force sensitivity [20], are then highly desirable to develop more sensitive vectorial probes. One interesting possibility is to exploit the non-conservative force induced by the detection system. Indeed, open systems can exhibit intriguing degeneracy points in the analytic continuation of their spectra, associated with the coalescence of the eigenstates and known as exceptional points (EPs). Recently, EPs have allowed efficient topological energy transfers between two harmonic modes of a membrane placed in the middle of an optical cavity [21].
Two-level systems (TLSs) constitute minimal quantum systems to detect and manipulate mechanical motions [22,23]. For example, their coupling to flexural modes allowed the localisation of emitters randomly distributed in micropillars [24,25]. It was also shown that single-molecule TLSs are sensitive local probes to measure through Stark effect the small displacements of a charged nanotube [26,27]. There are two fundamental differences compared to cavity optomechanics. (i) A strongly pumped optical cavity has a linear behaviour, whereas the quantum nature of a TLS is intrinsically nonlinear. (ii) In highly-populated optical cavities the poissonian fluctuations are negligible, whereas the TLS experience strong quantum fluctuations due to spontaneous emission. Whether it is possible to observe EPs in the electromechanical spectrum of a cantilever coupled to a TLS naturally appears as a fundamental question.
In this paper we show that adiabatically varying along a closed path the frequency and the intensity of a laser driving a TLS coupled to two flexural modes of a cantilever it is possible to induce a change in the state of the oscillator that depends on the topology of the path. We show that this behaviour is due to the presence of an EP in the mean-field description of the electromechanical spectrum. Using quantum Monte-Carlo simulation we prove that this property holds in presence of strong TLS fluctuations. The presence of the EP allows to generate elliptic mechanical eigenmodes, whose axis angles can be controlled. Finally, we propose a detection scheme to probe the topological switch between the two quasidegenerate flexural modes in single-molecule spectroscopy.
Electro-mechanical system. -We consider the system defined by the Hamiltonian: (1) The operators σ x and σ z are Pauli matrices describing a TLS of splitting ω TLS ( = 1). The TLS is coupled to a laser of frequency ω L with Rabi frequency Ω L . It also couples with strength g i two quasi-degenerate orthogonal flexural modes of a suspended nanowire of frequencies ω i and described by destruction (creation) operators b i (b † i ) [26]. In Hamiltonian (4) we omit the coupling of the TLS to the electromagnetic environment leading to incoherent spontaneous emission with decay rate Γ. Also the mechanical modes are subject to dissipation with damping rates γ i ( ω i ). This generic model can describe different physical systems [24,25,28]. For the figures and estimates below, we refer to a suspended carbon nanotube coupled by Stark effect to the static dipole moment of a single-molecule in a solid-state matrix, as described in Fig. 1. Typical parameter values are Γ/(2π) ≈ 10 MHz, ω i /(2π) ≈ 1-10 MHz, γ i /ω i ≈ 10 −3 − 10 −5 , and g i /(2π) ≤ 10 3 MHz [29]. Non-Hermitian dynamics. -In usual cavity optomechanics, the coupling between the oscillator and the cavity can be linearised to solve the problem exactly. This is not possible with a TLS, due to its intrinsic non-linear nature. We thus proceed with an approximation based on the time-scale separation γ i Γ [1,30]. The expected value of the displacement 2mω i denotes the zeropoint fluctuations displacement, satisfies a Langevin equation: The last term describes a coupling between the mechanical modes mediated by the retarded response of the TLS, namely is the fluctuations of the population difference. With . . . 0 we indicate averages taken for g i = 0. We evaluate them by solving the Born-Markov master equations for the driven TLS in the presence of dissipation, and using the quantum regression theorem [29]. In Fourier space this leads: where the polynomial P(z) = 3 n=0 a n z n is defined by the real coefficients a 0 = Γ(Γ 2 + (2δ) 2 + 2Ω 2 L )/4, a 1 = (5Γ 2 + (2δ) 2 + 4Ω 2 L )/4, a 2 = 2Γ, a 3 = 1, and δ = ω L − ω TLS is the laser detuning. The fluctuating forces δF i describe the thermal fluctuations and the non-equilibrium stochastic fluctuations due to spontaneous emission. We begin by neglecting these fluctuations and the average ofσ z that only shifts the equilibrium position of the oscillator. We then Fourier transform Eq. (2) and linearise it close to the poles ±ω i , assuming sufficiently small coupling so that the frequency renormalization remains smaller than the bare frequency. Since x i (t) is real we consider only the positive pole. Fourier-transforming back to the time domain, we obtain an equation of motion for the displacement vector x = 2Re[X], where x = (x 1 , x 2 ) and X(t) describes the positive-frequency oscillations of the displacement [29]. It formally obeys a Schrödinger equation iẊ = HX with non-Hermitian effective Hamiltonian where S R is evaluated at ω 1 ≈ ω 2 , since ∆ω = ω 1 − ω 2 Γ. EP in the electromechanical spectrum. -To study the effective non-Hermitian dynamics, we first perform the unitary transformation H = UHU † . Here U = e i π 4 τ y e i π 4 τ z , and τ i are Pauli matrices acting in the mechanical-mode space (Bloch sphere). It leads to where h ± = ω 1 − ω 2 − i(γ 1 − γ 2 )/2 − (g 2 1 − g 2 2 ± 2ig 1 g 2 ) S R , and h 0 = TrH. The electromechanical spectrum is characterized by two eigenvalues satisfying 2λ For non-vanishing coupling (g 1 g 2 S R 0), quantities h + and h − cannot vanish simultaneously. It is then clear that the degeneracy points S EP± correspond to the branch points of a complex square root and are EPs [31][32][33]. Therefore, the resonance frequencies Ω ± = Re[λ ± ] and the damping rates Γ ± = − Im[λ ± ] support a Riemann-surface representation in the vicinity of the EP, as illustrated in Fig. 2 (a). If δ and Ω L are varied such that S R performs a loop around S EP± in the complex plane, the eigenvalues λ + and λ − switch from one to the other due to the Riemann-surface topology. This is shown in Figs. 2 (b) and (c). For g 1 g 2 the condition for the EP appearance in the positive imaginary plane reads S EP i∆ω/(2g 2 i ). Since S EP is nearly pure imaginary, to vary S R around this value, one must be able to change the sign of the real part of S R and to increase its modulus up to |S EP |. From Eq. (3), these two conditions give ω i > Γ/2 and Γ∆ω 2g 2 i , using the maximum of S R is obtained when δ ∼ Ω L ∼ Γ. This implies a typical coupling of g i /(2π) ≈ 0.3 MHz. In the case of the dipolar TLS of a dibenzoterrylene molecule in 2,3-dibromonaphthanlene crystal, it is four orders of magnitude below the discharge coupling between the nanotube and the substrate [29]. For instance, it can be realized by positioning the tip of a carbon nanotube 100 nm away from the molecule with a 100 µV bias.
Evolution of electromechanical modes. -EPs also lead to intriguing properties for the eigenstates. We write the orthogonal left-and right-eigenstates of the non-hermitian matrix H as Y ± (θ) = (±Z 1 4 , where Z = −(2S EP e −iθ /ρ + 1) and ρe iθ = S R − S EP . When S R encircles n times the EP clockwise, the multivalued eigenstates change as X ± (θ+2nπ) = −i n X ± (θ) if n is even, and as X ± (θ+2nπ) = −i n X ∓ (θ) otherwise [29]. It implies that, similarly to the eigenvalues, the eigenstates also switch after one cycle. Such a switch cannot happen in Hermitian systems, where the real eigenvalues cannot exchange without passing through a degeneracy. Besides, since either h + or h − vanishes at the EP, Eq. (5) readily shows that H supports a Jordan matrix representation, meaning the two eigenstates coalesce [34]. Now let us see how the eigenstate multivaluation affects the flexural modes of the cantilever in real space. The eigenstates are solutions of the effective Shrödinger equation: , the real-space dynamics satisfies R e S R Im SR ± ± R e S R Im SR Here Thus, each electromechanical mode exhibits a damped elliptical motion. The ellipse semi-axes α and β are ϕ-rotated with respect to the orthogonal polarizations of the uncoupled mechanical modes.
Topological and chiral actuation. -Importantly, equation (7) provides an explicit relation between the rotation angle ϕ of the elliptical modes and the phase θ of the TLS retarded response, as specified by the colormap in Fig.2 We find that, if S R encircles n times S EP , the semi-axes lengths remain unchanged: α(θ + 2nπ) = α(θ) and β(θ + 2nπ) = β(θ). However, the semi-axes rotate as ϕ(θ + 2nπ) = ϕ(θ) − nπ/2. This four-fold invariance arises from the complex fourth-root multivaluation of the eigenstates X ± . It implies in particular that the bare orthogonal modes undergo a π/2-rotation when S R loops once around S EP . This rotation comes with the switch of the resonance frequencies and damping rates. Thus, a bare flexural mode initially activated as where ∆φ is the phase accumulated along the loop. This operation only depends on whether S R encircles S EP , regardless of the precise loop geometry, and is therefore topological. This intriguing feature is exemplified in Fig.2 (d). Interestingly, the eccentricity of the elliptical motion vanishes when approaching the EP where α = β: the two electromechanical modes merge into a single damped mode of circular polarization, as shown in Fig.2 (e). This results from the eigenstate coalescence of the Jordan matrix representation of H at the EP.
The topological actuation from one mode to another cru-cially relies on the adiabatic transport of the eigenstates X ± by varying the laser parameters. Nevertheless, gain and loss can break the adiabatic theorem in open systems [35][36][37]. Here we investigate this issue by solving numerically the equation iẊ = HX [29]. For an adiabatic eigenstate, the effect of the dissipation along a loop of period T in parameter space is γ ± T = Γ ± (t)dt. If state +(-) initially corresponds to state 1(2), and assuming γ − < γ + , then only can state -experience adiabatic transport from state 2 to state 1, because state + is more damped than state -by a factor e (γ − −γ + )T along the loop. Reversing the loop orientation, i.e. reversing the arrows in Fig. 2 (c), implies γ − > γ + . Then, for the same reason, only can state + experience adiabatic transport from state 1 to state 2. Thus, both flexural modes can experience the topological swap, but for two opposite orientations of the loop. This chiral behavior around the EP shows that the topological mechanical actuation mediated by the TLS is reversible. Back-action and fluctuations. -The present picture is very appealing but, in contrast with optical cavities coupled to mechanical oscillators that can be exactly solvable, it relies on the approximate Eq. (2). It neglects the change of the detuning induced by the displacement δ[x(t)] = δ − i g i x i (t)/x zpf i , but even more important, the TLS dynamics generates strong force fluctuations. Defining n = (1 − σ z )/2 as the occupation operator of the exited state, the variance squared of the force reads (∆F i ) 2 = (g i /x zpf i ) 2 n (1 − n ). This should be compared with the average force squared (g i /x zpf i ) 2 n . Since n < 1/2, one obtains that (∆F i ) 2 / F i 2 > 1, and actually diverges for n → 0. In contrast, for the cavity one finds that the same ratio reads 1/ n c 1, where n c is the number of circulating photons in the cavity. Besides, Eq. (2) neglects the one may reduce the thermal fluctuations by cooling the mechanical oscillator through the TLS [26,29].
To test the transparent mean-field description based on Eq. (2), we perform simulations with the Monte-Carlo quantum-jump approach [38][39][40]. This implies solving on short time scales ∆t 1/Γ the differential equations where we introduced the TLS wave-function |ψ(t) = a g (t)|g + a e (t)|e , the thermal white-noise force ∆F th with variance 2mγ i k B T i , and σ z (t) = |a g (t)| 2 −|a e (t)| 2 . At each time step one randomly allows either a transition to the ground state with probability Γ∆t|a e | 2 , or normalize the wave-function and proceeds with the time evolution.
We begin by verifying that one can cool the nanowire down to temperatures of the order of Γ (0.5mK) for ω i = Γ, ∆ω = 10 −3 Γ, γ i = 10 −5 Γ, g i = 0.1Γ, δ = −Γ, and Ω L = Γ [29]. We then simulate the evolution of the oscillator with initial amplitude 10x zpf for mode 1 or 2, by adiabatically varying δ from 0 to -8Γ over T = 10 4 /Γ for fixed values of Ω L . The mean dynamics of the oscillator, averaged over 10 2 simulated trajectories, shows quasi-periodic elliptical oscillations in real space. We thus identify the semi axes in each period and show in Fig. (3) its rotation angle ϕ at a given value of δ and Ω L . In panel (a) we present the prediction of the mean-field description, clearly showing the presence of the EP at δ −0.85Γ. Panel (b) shows the rotation angle for the same parameters of panel (a), but calculated with the Monte-Carlo simulations. Though fluctuations are present, it is quite clear that the axis of the oscillator performs on average a rotation around this point, and that starting with an oscillation in the horizontal direction (ϕ = 0) one ends with a perpendicular oscillation (ϕ = π/2), showing that the adiabatic picture remains valid. Of course the rotation implies a transfer of energy from mode 2 to mode 1, but continuously following the evolution of the axis allows to prove without doubts that this transfer is purely due to the adiabatic evolution predicted by the mean-field approach, and not to dissipation or stochastic effects. In Fig. (3) panel (c) we show the same evolution from mode 1, for which no energy transfer is realized, as expected from the mean-field adiabatic transport.
Oscillator's frequency detection. -After the adiabatic loop, the energy of one mode transfered to the other one, implying that the resonator oscillation frequency changed. Here we show that detecting such a small change can be performed by laser frequency modulation. When the laser frequency is modulated at the frequency ω λ , it induces an oscillation that adds up to the contribution of the mechanical oscillator in the Hamiltonian in the modulated rotating frame: δ(t) = δ − (β λ /ω λ ) cos(ω λ t) − g i x i (t) [29]. We assume only mode i is exited at the end of the loop. It is then clear that when ω λ = ω i the laser frequency modulation and the mechanical oscillations interfere. The interference is possible if |ω λ − ω i | < γ i . It leads to the sharp peaks in the luminescence spectrum as a function of ω λ in Fig. (3) (d). Importantly, the width of these peaks is γ i ∆ω, so that the two eigenmodes can be resolved clearly in single-molecule experiments.
Conclusion. -We have shown that coupling the flexural modes of a nanowire to a laser-driven TLS it is possible to induce a topological (and chiral) swap between the two quasidegenerate orthogonal modes. The phenomenon relies on the presence of an EP in the electro-mechanical spectrum as predicted by a mean-field approach. Remarkably, the solution of the problem by quantum Monte-Carlo simulations shows that the prediction is robust to the strong non-linearities and fluctuations of the TLS. The prediction of the existence of such an EP in a non-linear system does not only has a theoretical interest by itself, but also opens the way to exploit topological properties in the manipulations of nano-mechanical tips with promises for high-resolution applications in vectorialforce microscopy. For the choice of typical parameters, we refer to the electromechanical system in Fig.1 in the main text. It has recently been shown that single molecules, such as dibenzoterrylene in 2,3-dibromonaphthanlene crystal, can gain a permanent electric dipole δµ up to 2D [41]. The electromagnetic environment coupled to the molecule leads to the typical decay rate Γ/(2π) 8-10 MHz for experiments performed in the temperature range 10 mK-1 K. As suspended nanowires, a singly-clamped carbon nanotube (NT) exhibit two fundamental flexural modes that are orthogonal to one another. Their nearly degenerate frequencies ω i are typically in the range 1-10 MHz and the NT mass is m 10 −20 kg [42]. Besides, we focus on the under-damped regime and consider typical quality factor Q i = 10 3 -10 5 .

Retarded response of the TLS
To determine the effective coupling between the flexural modes we evaluate S R through Born-Markov methods. After tracing out the electromagnetic environment, the reduced density matrix of the TLS ρ = (ρ 11 , ρ 12 , ρ 21 , ρ 22 ) obeys a Liouville-von Neumann equationρ = Lρ associated to the super-operator where δ = ω L − ω TLS defines the laser detuning. This corresponds to the optical Bloch equations that are known to provide a realistic description of electric dipoles in single-molecule experiments [44,45]. In the stationary regime Lρ 0 =0, the quantum regression theorem leads to the autocorrelation function where w 0 = (1, 0, 0, 1) denotes the kernel left-hand eigenvector of L, and M z = diag(1, 1, −1, −1) [27]. The power spectral density S (ω) = dt e iωt S (t) characterises the absorption (ω>0) and emission (ω<0) of the dissipative driven molecule. The frequency asymmetry of the quantum noise is controlled by the laser detuning δ. It relates to the imaginary part of the response function through 2 Im S R (ω) = S (ω) − S (−ω). From Kramers-Kronig relations, we finally obtain where the polynomial P(z) = 3 n=0 a n z n has real coefficients a 0 = Γ(Γ 2 + (2δ) 2 + 2Ω 2 L )/4, a 1 = (5Γ 2 + (2δ) 2 + 4Ω 2 L )/4, a 2 = 2Γ, and a 3 = 1. The behaviour of S R ≡ S R (ω 0 ) is shown in Fig. S4 as a function of the laser detuning and amplitude. It vanishes at resonance and for strong detuning, so that the flexural modes are uncoupled and have orthogonal linear polarizations. This explains the oriented loops outlined by S R in Figs. S5 (a) and (e), as well as in Fig. 2(b) in the main text, when red detuning the laser between resonance and −8Γ. Note that, contrary to optical cavities, the saturation of the TLS for Ω L /Γ 1 leads to S R → 0, thus decoupling the mechanical modes.

Non-Hermitian dynamics
The expectation value of the displacement operator obeys Eq. (2) in the main text. We begin with neglecting the fluctuating forces δF i and σ z 0 that only shifts the equilibrium position of the oscillator. In Fourier space, this leads to Since we consider the under-damped regime (γ i ω i ), the bare mechanical poles verify: ω i± ±ω i − iγ i /2. Then we can linearize Eq. (S11) in the vicinity of the mechanical frequencies ω i± , and show that the positive-and negative-frequency parts x i± of the displacement x i = x i+ + x i− satisfy the following set of conjugated equations: where S R ≡ S R (ω 0 ) and ω 0 = (ω 1 + ω 2 )/2 is the mean mechanical frequency. Since where we have used 2mω i x (S14) The retarded response of the TLS mediates an electromechanical coupling between the orthogonal flexural modes as long as g 1 g 2 S R 0.

Topological swap of the eigenstates
First, we introduce the unitary operator It represents two subsequent π/2-rotations of axes y and z on the Bloch sphere of the mechanical modes. Then we perform the transformation H = UHU † = i h i τ i /2 and find (S16) We can rewrite it in a more compact matrix form where h ± = h 1 ± ih 2 . This leads to Eq. (5) in the main text. Here we are interested in the eigenstates in the presence of the exceptional points S EP± in the electromechanical spectrum, where (S18) In the main text, we focus on the exceptional point S EP− ≡ S EP that lies in the positive imaginary plane. It is then convenient to set S EP as new origin of the complex plane. Introducing the variable z = S R − S EP , we find where we have defined G ± ≡ −(g 2 1 − g 2 2 ± i2g 1 g 2 ) and ∆S ≡ S EP − S EP+ . The electromechanical spectrum is characterized by the two eigenvalues Thus, it is explicit that the exceptional point S EP (z = 0) is nothing but the branch point of the complex square root √ z in the spectrum, where the effective Hamiltonian reduces to a Jordan matrix form.
The bi-orthogonal left and right instantaneous eigenstates of the effective Hamiltonian can be chosen as where we assume polar coordinates (z = ρe iθ ) so that Let us assume that S R can be varied smoothly along a ρ-radius loop that encloses only S EP (ρ < |∆S |). Then the eigenstates, which meet the condition of parallel transport Y ± · ∇ θ X ± = 0, change as follows: Thus, we expect the two eigenstates to swap after one loop, to pick up a geometrical Berry phase π after two loops, and only after four loops come the eigenstates back onto the initial states without any geometrical phase. The eigenstate swap clearly results from the multivaluation of the eigenstates associated to the EP. The four-fold multivaluation comes from the complex fourth root in the eigenstate definitions (S21). Note that if the loop encircles the two exceptional points (|∆S | < ρ), it is then possible to make the eigenstates single-valued and X ± (2πn) = +1 X ± (0) for any integer n. The eigenstates neither pick any geometrical phase nor swap in this case.

Chiral swap of the eigenstates
The swap introduced above is an intrinsic property of the instantaneous eigenstates X ± . In order to observe the effects of their multivaluation around an EP, one can study their adiabatic transport. In Hermitian systems, the adiabatic theorem ensures that one can neglect the non-adiabatic transitions over some typical time scale T 1/|λ + − λ − |. In open systems, however, this is no longer true for all the eigenstates [21,35,36]. For a two-state system, in particular, only is the least dissipative state expected to be transported adiabatically and swap [37]. Here we investigate this issue by numerically solving Eq. (S13).
We consider a fixed laser amplitude Ω L = 0.3Γ, and we ramp the detuning linearly between δ = 0 and δ = −8Γ within T = 10000/Γ. This way, the response function S R performs a loop around S EP (cf. Figs. S5 (a) and (e), and Fig.2 (b) in the main text), resulting in the swap of the eigenvalues and eigenstates. We describe the exact dynamics by |X ± (t)| 2 = |Y ± (t)U(t, 0)X(0)| 2 , where U denotes the time-evolution operator associated to H. We choose the initial states as satisfying either X(0) = X + = X 1 or X(0) = X − = X 2 , where ω 1 = 1.000Γ and ω 2 = 0.999Γ. We then compare their exact dynamics to the adiabatic evolution described by |Ad ± (t) The results are shown in Fig. S5 for two orientations of the S R loop around S EP . Varying the detuning from 0 to −8Γ leads to a clockwise loop, as shown in Fig S5(a). Figures S5(b)-(d) show that only the least dissipative eigenstate along loop, namely X − (in blue), can be transported adiabatically: From Fig. S5(e)-(h) we can check that the least dissipative state becomes the most dissipative one when reversing the orientation of the loop, thus preventing this state from following its adiabatic evolution: This asymmetry in the adiabatic breakdown with respect to the orientation of the loop around the EP reveals the chiral nature of the eigenstate swap.

Polarization rotation of the electromechanical modes
The instantaneous eigenstates obey the Shrödinger-like equation iẊ(t) ± = HX ± (t) and explicitly read X ± (t) = x ± (0) e −iλ ± t X ± . In real space, this leads to the eigenmodes where when ignoring the π/4 phase shift. We then introduce the counter-clockwise rotation matrix of angle ϕ so that This characterizes Eq. (7) in the main text. The flexural modes, which are linearly polarized along transverse directions in the absence of coupling (S R = 0), become elliptically polarized with orthogonal semi axes when coupled to the molecular TLS. If S R is varied clockwise on a loop around S EP , we can show from Eq. (S22) that Z 1/4 (θ = −2π) = iZ 1/4 , which leads to Thus, using Eq. (S30), we obtain and similar relations for β. This shows that the semi-axes α and β do not change after the loop, whereas their elliptical polarization has been rotated through the angle π/2 (cf. Fig. 2 (d) in the main text). Finally, note that in the limit of the loop radius vanishes (R → 0), Eq. (S27) leads to a + a − and φ + φ − . Thus the excentricity vanishes (α β) and the two eigenmodes coalesce into a single circularly polarized flexural mode (cf. Fig. 2 (e) in the main text).

Back-action cooling
The power spectrum of the TLS is obtained from Eq. (S9), as S (ω) = +∞ −∞ dte iωt S (t). It is convenient to introduce the symmetric-and asymmetric-in-frequency parts of the power spectrum: The optical damping then satisfies where ω 0 = (ω 1 + ω 2 )/2 is the mean mechanical frequency. Mechanical cooling via the TLS requires γ TLS γ 0 , so that For the typical values of the parameters we consider, ω 0 γ 0 is at least two order of magnitude smaller than critical coupling g c . Regarding the temperature of the TLS, it verifies Then one can show that the effective temperature satisfies  Figure S6 (a) represents the effective temperature as a function of δ and Ω L . It shows that the mechanical system can be cooled even below Γ( 0.5 mK). We can very this possibility by means of quantum Monte Carlo simulations. Each mechanical mode is assumed to be initially in thermal equilibrium with its bath at temperature T n . Thus, we begin with randomly generating the position and momentum of each mode according to Boltzmann distribution. We then simulate the dynamics of the two modes coupled to the TLS using the quantum Monte Carlo approach described in the main text. The results for δ = −Γ and Ω L = Γ are shown in Fig. S6 (b) after averaging over 1000 dynamics. The figure presents the evolution of dimensionless energies E n . They relate to the mechanical energies as follows: In the beginning of the dynamics we assume thermal equilibrium with E n = 4ω n Γ k B T n Γ = 80, which corresponds to T n 10 mK. Figure S6 (c) shows that, at the end of the mean dynamics, the system has been cooled below Γ ( 0.5 mK), down to about T 50 µK.

Laser frequency modulation
The eigenstate swap performed when S R loops around the EP is also accompanied by the exchange of the resonant frequencies. The eigenstate swap can then be probed by measuring the mechanical frequency after the loop. Such frequencies can be obtained from the luminescence excitation spectrum in single-molecule experiments [26]. It has been shown to exhibit side-band peaks of width Γ centred on multiples of the flexural frequency [26]. Here, however, we need to resolve two quasi-degenerate frequencies, and the side-band shift cannot be resolved since |ω 1 − ω 2 | Γ. In order to resolve it accurately, we propose to modulate the laser frequency, whose interferences with the mechanical oscillations can be revealed in the luminescence spectrum. At the semi-classical level, the Hamiltonian of the TLS driven by a frequency modulated laser can be written as where F(t) = t 0 dτ [ω L + f (τ)] and f (t) = a λ cos(ω λ t + φ λ ). This describes the frequency modulation with respect to the laser frequency ω L . Let ψ(t) obey the time-dependent Shrödinger equation Using a rotating wave approximation and introducing the gauge transformation φ(t) = U † (t)ψ(t) based on the unitary operator we obtain where the effective Hamiltonian isH Here the time-dependent detuning is The luminescence excitation spectrum scales linearly with the excited-state population of the TLS. The equations of motion for the reduced density-matrix components σ 12 = σ * 21 and σ 22 = 1 − σ 11 are obtained fromH TLS within the Born-Markov approximation. They correspond to the optical Bloch equations: They can be rewritten as where X = σ 12 = σ * 21 , Y = σ 22 = 1 − σ 11 , τ = Γt, and = Ω L /Γ. Here we look for perturbative solutions of the type X = n n X n and Y = n n Y n in the limit 1.

Order n=0
The Bloch equations lead to Y 0 = Y 0 (τ 0 ) e −(τ−τ 0 ) , In the limit τ 0 → −∞, they reduce to The Bloch equations lead to The Bloch equations lead to and the solutions generically read In the limit τ 0 → −∞, they reduce to .

(S53)
At the end of the loop around the EP, either the flexural modes have flipped or they have not. In both cases, the flexural dynamics is of the type x(t) = Ae −γ i t cos(ω i t + φ i ). From now on we consider φ i = 0, which fixes an arbitrary origin of time. The luminescence excitation spectrum scales linearly with where β λ = a λ /ω λ , and we introduced the dimensionless parameter β = 2gA/ω i . In the limit β 1, the luminescence excitation spectrum can be approximated by where the sum involves positive and negative values of the integers m and n. Besides, J k denotes the k-th order Bessel function of the first kind, and A 0 = e −i(m−n)ω λ t (δ + nω λ − i Γ 2 )((m − n)ω λ + iΓ) .
The term A 0 does not depend on β and describes the luminescence excitation spectrum in the absence of electromechanical coupling. We focus on ∆σ 22 = σ 22 g 0 − σ 22 g=0 . We can then accumulate the luminescence over time, which results in