Thermally induced entanglement of atomic oscillators

Laser cooled ions trapped in a linear Paul trap are long-standing ideal candidates for realizing quantum simulation, especially of many-body systems. The properties that contribute to this also provide the opportunity to demonstrate unexpected quantum phenomena in few-body systems. A pair of ions interacting in such traps exchange vibrational quanta through the Coulomb interaction. This linear interaction can be anharmonically modulated by an elementary coupling to the internal two-level structure of one of the ions. Driven by thermal energy in the passively coupled oscillators, which are themselves coupled to the internal ground states of the ions, the nonlinear interaction autonomously and unconditionally generates entanglement between the mechanical modes of the ions. We examine this counter-intuitive thermally induced entanglement for several experimentally feasible model systems, and propose parameter regimes where state of the art trapped ion systems can produce such phenomena. In addition, we demonstrate a multiqubit enhancement of such thermally induced entanglement.


I. INTRODUCTION
Quantum superposition phenomena are typically obtained by precisely and coherently driving nonlinear systems well-isolated from the environment at low temperatures. Quantum technologies taking advantage of such phenomena will be more practical if the quantum components instead operate autonomously [1]. To do so, they require resources such as quantum superposition and entanglement generated without an external coherent drive, and emerging purely from thermal energy. While counter-intuitive, as thermal energy is typically destructive to such resources, thermally induced entanglement can be generated quite straightforwardly. A prototypical example is a pair of two-level systems, one in the ground state and the other in a thermal state, which, through unitary interactions can exhibit greater entanglement for greater temperatures. Theoretical extensions of such cases to the interaction between a two-level system and an oscillator appeared some time ago [2], as also the case of entanglement between a pair of qubits in pure states induced by coupling to a thermal field [3]. Further studies show the manner in which separable pure states of pairs of oscillators or qubits generate enhanced entanglement after being exposed either to a thermal environment [4] or when interfaced directly with a single-mode thermal field [5][6][7]. Similar results hold for multiqubit collections of atoms [8] as well. However, deterministic thermally enhanced entanglement arising from heating passively coupled oscillators with a Hilbert space of infinite dimension, coupled to cold atoms, remains hitherto unexplored. It is a nontrivial complementary task, as an * pradip.laha@upol.cz † slodicka@optics.upol.cz ‡ darren.moore@upol.cz § filip@optics.upol.cz excited two-level system only absorbs less than a single energy quantum from coupled thermal oscillators.
With passive interactions, thermal oscillators do not entangle spontaneously. Therefore, oscillator-oscillator entanglement driven by thermal energy can be extremely useful. For example, conditional thermally induced entanglement can lead to Bell-inequality violations [9]. After many years of experimental development such experiments have become available using highly nonlinear controlled-swap interactions [10,11]. A simpler approach can be used to observe the thermally induced mutual coherence of atoms [12]. However, a deterministic version of thermally induced entanglement or mutual coherence for oscillators is still hard to imagine. An early step towards this intriguing possibility was the thermally induced unconditional generation of nonclassical oscillator states from the elementary Jaynes-Cummings (JC) interaction of the oscillator with a ground-state two-level system [13,14]. The nonclassicality emerges from a complex anharmonic modulation of the oscillator energy distribution by a simple, coherently interacting, two-level system initially in the ground state. In this method, the smallest nonlinear unit interacts with an oscillator to produce such phenomena. Using the thermal energy already present in the oscillator contrasts with the external environment used in reservoir engineering experiments [15]. Autonomous and unconditional generation of entanglement or coherence are the important resources to consider in scaling such behaviour to more complex systems. Here we focus on inducing entanglement in a pair of thermal oscillators using such methods.
An array of laser driven ultracold ions confined in a linear Paul trap (see Fig. 1) is an excellent and experimentally viable platform available for the implementation of a wide range of quantum protocols involving quantum simulation [16,17], for example of Ising spin systems [18][19][20][21], the Bose-Hubbard (BH) model [22][23][24] and the Jaynes-Cummings-Hubbard (JCH) model . A schematic of an ion trap setup that can realise the proposed deterministic generation of entanglement by increasing the temperature difference between two oscillator modes. Two ions in a linear trap align along the axial (ωz) direction and interact via the Coulomb force (κBS), sharing excitations among the radial vibrational modes. Alternatively, the internal energy of a single ion can be addressed in order to stimulate a Jaynes-Cummings interaction (κJC) between a pair of radial modes. A single qubit is prepared in the ground state, while the relevant radial modes (see main text) are prepared in initially uncorrelated thermal states. By selecting one of the interaction arrangements one of our models is implemented. From the full system dynamics of each model we discard the qubit contribution to calculate the entanglement of the joint two-field subsystem in terms of the logarithmic negativity. Note that system 1M (see text) is implemented with just a single ion with both radial modes coupled to the internal state of the ion, while system 1S requires the radial modes of two ions. 31]. Recently, a single ion has been used to investigate quantum phase transitions in mechanical motion [32]. These systems possess several striking advantages such as the presence of long coherence times, the ability to efficiently manipulate the internal and motional (centre of mass) degrees of freedom, and individual access to each ion allowing implementation of quantum gates between specific ions [33,34]. Spontaneous emission of the internal states and damping of the motion of the ions are effectively minimised in state of the art experiments by tuning the appropriate metastable transitions of the ions stored in ultrahigh precision spectroscopy and maintaining low pressures for long times such that collisions with background atoms and other incoherent couplings are negligible [35]. The high level of control afforded by such systems combined with their anticipated scalability in simulating many-body effects provide us with a framework in which the autonomous quantum behaviour that may become a necessity of larger-scale quantum technologies can be investigated. While we focus entirely on trapped ion systems in our analysis, we also note the potential application of our scheme to superconducting circuits [36] which share many properties we make use of below.

II. THERMALLY INDUCED ENTANGLEMENT
A collection of ions can be confined in three dimensions by a trap that is harmonic at the level of the motion of the ions. That is, the ions are confined by harmonic potentials of frequencies ω α , α = x, y, z. If ω z ω x,y then the ions form a chain by aligning themselves along the axial (z) direction, with equilibrium positions (0, 0, z 0 i ). The vibrational motion of the ions occurs around equilibrium points situated along this chain.
In the following derivation we focus on a qubit interacting with a pair of oscillator modes. In the case of a beamsplitter interaction two ions are required, however we emphasise that a single ion can be leveraged to implement an alternate configuration involving only Jaynes-Cummings interactions (see schematics of Fig. 2). A pair of confined ions undergo vibrational motion in the radial directions x and y. The Coulomb interaction can induce the hopping of vibrational quanta, phonons, between the corresponding radial modes. This can straightforwardly be seen through a Hamiltonian analysis of the system. The generic Hamiltonian including both the trapping potential and the Coulomb interaction between the ions is given by where α i and p α,i are respectively the α components of the displacement from the equilibrium position and the momentum of the ith ion, and m and e are respectively the mass and charge of an ion, while the r i = (x i , y i , z i ) are the positions of the ions relative to the equilibrium position. The distance term in the Coulomb potential can be expanded around the displacements from equilibrium to include up to harmonic terms. This approximation leads to the decoupling of the individual modes which is a good approximation for axial potential frequencies much smaller than radial ones and for small thermal energies of ions. Now, decomposing the dynamics into one of the transverse directions, say x, we have the Hamiltonian [37,38] (see the Supplementary Material for a detailed derivation) The Coulomb effect therefore induces frequency shifts in the oscillators, ω 2 With this the Hamiltonian can be written in second quantization by promoting x i and p i to dimensionless quadrature opera- . Neglecting the zero-point energy of the oscillators, we may write In a frame rotating with respect to the free motion, the interaction Hamiltonian is composed of resonant and fast-rotating terms. In the rotating wave approximation we retain only the resonant terms which constitute a beamsplitter-like interaction allowing excitation transfer between the modes, Simultaneously, if one considers the internal energy levels of an ion, addressed by an external electromagnetic wave, then a discrete-variable nonlinearity can be introduced [39]. Addressing a pair of energy levels leads to an interaction Hamiltonian in the rotating wave approximation of the form where δ is the detuning of the driving from the frequency of the two-level system and Ω is the coupling strength.
Assuming the Lamb-Dicke approximation corresponding to η x i 2 << 1 and setting the laser to the resonance with the first red motional sideband, δ = −ω x , gives where κ JCi = Ω η. A coherent JC coupling to the second radial mode can simultaneously be implemented with a bichromatic laser excitation that addresses two red sidebands simultaneously. While the oscillator interaction is passively linear, one may take advantage of the nonlinear coupling inherent in the JC model in order to use the temperature difference of the oscillator pair to induce entanglement. Various configurations can be imagined and we analyse and contrast two basic instances below.
The entanglement production is driven by the difference between the initial temperatures of the subsystems. Thermal states of the oscillator and qubit, in thermal equilibrium with nonzero temperatures T O i and T q respectively, are characterised by the Boltzmann distributions where |g and |e are the energy eigenstates of the qubit. The mean excitationn i of the ith oscillator is related to The probability of finding the qubit in the ex- . In trapped ion systems, the characteristic decoherence rates are small compared to the short interaction times required to generate entanglement [35], and hence we assume unitary dynamics (see, however, the Supplementary Material).

A. Thermally induced entanglement via a single qubit
The simplest possible generalisation of the thermally induced entanglement protocols involves a pair of oscillator modes and a single qubit. The two elementary arrangements of these three systems involve the qubit coupled to one member of the oscillator mode pair which is coupled to the other member via the beamsplitter interaction, or a single qubit acting as mediator between a pair of uncoupled oscillator modes (Fig. 2 schematics). The interaction Hamiltonians of such systems can be expressed as . The qubit is initially in the ground state, while the coupling ratios are given by r1S = 0.18 (0.5). After a short latent period entanglement becomes visible, and increasing the temperature difference significantly increases the LN. (c) The logarithmic negativity for H1M with coupling ratio r1M = 1. The symmetry of H1M ensures that no different behaviour is observed by swapping the roles ofn1 andn2 whereas the asymmetry of H1S is visible by the reduced amount of entanglement produced when the thermal energy is concentrated in the oscillator distant from the qubit. This asymmetry also modifies the optimal coupling strength, requiring stronger coupling in order to maximise the induced entanglement when the thermal noise is initialised in oscillator 2.
Here S (M) stands for a qubit at the side (in the middle). H 1S , in taking advantage of the beamsplitter interaction requires the interaction of the radial modes of two ions, whereas H 1M requires a single qubit and is implemented using the radial modes of a single ion. To study the entanglement dynamics we measure the entanglement between the reduced state ρ of the oscillator pair by numerically estimating the logarithmic negativity (LN) [40,41], which is defined by where A is the partial transpose operation and ||ρ|| 1 denotes the trace norm of the operator ρ. Contrary to intuition, increasing the temperature difference of the oscillator modes enhances the entangle-ment produced by the unitary evolution of these systems. We consider these systems in turn below. In discussing system 1S it is useful to first define a single timescale between those of κ JC1 and κ BS12 by introducing r 1S = κBS 12 κJC 1 such that the dynamics occurs in units of τ = κ JC1 t. For system 1M we write r 1M = κJC 2 κJC 1 . The symmetry in H 1M means that r −1 1M simply represents a relabelling of the oscillators. It is possible, without loss of generality, to study only the effects of p e andn 2 with respect ton 1 with 0 < r 1M ≤ 1.
We first consider p e = 0, so that the atom is in its ground state and examine the cases where eithern 1 = 0 orn 2 = 0 with the remaining oscillator having a nonzero mean phonon number. Figs. 2 (a-c) illustrate the generation of LN over scaled time τ for different initial distributions of the thermal noise in the first and second oscillators. Firstly, for system 1S (top), we note that the general behaviour indicates that thermal noise in the first oscillator, the one directly coupled to the qubit, is most significant for generating entanglement. This is supported by the noise distribution in favour of oscillator 1 producing much more entanglement with a much weaker qubit-oscillator coupling rate r 1S , compared to the initial thermal noise in oscillator 2. For system 1M (bottom) the symmetry of the configuration means that the two variations are identical. The LN is noticeably less than that of system 1S. For both systems the entanglement detected by the LN emerges after a short delay, and then experiences periodic collapses and revivals. This delay indicates that short time approximations of the JC interaction will not produce nonzero LN of the oscillators and that this is instead a property of the full complex dynamics introduced by the qubit. Interestingly, the first rise in entanglement for system 1M takes approximately half the time of system 1S. Increasing the thermal energy increases the entanglement generated, but notably does not change the profiles of the collapses and revivals, which are monotonic with respect to the thermal noise, up to saturation; when the temperature becomes too high, the entanglement saturates.
In what follows we select the optimal evolution time as the first maximum in the entanglement dynamics to analyse the intermediate parameter regimes. In some cases ( Fig. 2 (b), for example), a later maximum can be higher than the first but we wish to minimise any potential detrimental effects from decoherence of the system. Figs. S1 (a-b) and Figs. S1 (e-f) illustrate the manner in which entanglement increases with an increase in the thermal noise. Regardless of the coupling strengths, increasing the thermal noise increases the entanglement, up to saturation and up to a critical value (black line) above which entanglement generation begins to be inhibited again. The second and fourth rows indicate that increasing the mixedness in the remaining components of both systems decreases the capacity to generate entanglement. In fact, ifn 1 =n 2 there is no entanglement generation even with the qubit in the ground state. Both the ground state of the qubit and the temperature difference are important for entanglement generation. As said however, thermal noise in oscillator 2 for system 1S is less effective in generating entanglement, as seen in panels (b) and (d), so that while the temperature difference drives the entanglement generation, the distribution of the thermal noise must also be such that it interacts directly with the qubit induced nonlinearity. This is supported by panel (d) which shows that adding thermal noise to the more effective oscillator 1 does not increase the entanglement, but rather is consistent with our interpretation that it is the temperature difference that drives the entanglement generation with one exception. Fig. S1 (d) shows that for very small values ofn 2 , increasingn 1 is useful. This supports the argument that it is both the temperature difference and the direct interfacing of the nonlinearity with the thermal energy that are important for entanglement generation. Larger values of n 2 obviate this effect. In contrast to the other results, panel (b) shows a crossing in the LN. Similar to increasingn 1 , strong coupling can compensate for the incorrect distribution of the thermal noise for lower values ofn 2 . The symmetry of 1M means that when the roles of oscillators 1 and 2 are swapped the results do not change. This is further supported by the fact that the optimal coupling ratio is the balanced r 1M = 1.

B. Experimental feasibility
Experimental progress in the last decade towards implementing the BH or JCH model for trapped ions has resulted in the observation of phonon-hopping between radial modes of two trapped ions [42,43] or between the pair of radial modes in a single ion [44]. The significance of such models is typically held to lie in their many-body qubit properties such as Mött insulator and superconducting phases [31,37]. Such nonlinear interactions along with the ability to directly address specific ions also allow a collection of trapped ions to act as a quantum simulator. In this article, we have rather focused on applying these developments to a small-scale few-body quantum simulator of nonlinearly induced effects in linear oscillators. A similar effort can be seen in a recent experiment on quantum phase transitions with a single ion [32]. These experiments are well-suited to investigate the autonomously generated entanglement we have studied here. In particular, system 1S has already seen experimental realization in [42] and [43], demonstrating the hopping of excitations and the Hong-Ou-Mandel effect respectively. Two 40 Ca + ions, confined axially by a DC electric field and radially by radio frequency electric fields, share radial mode phonons mediated by the Coulomb interaction. The axial frequency is engineered to be much lower than the radial frequencies, so that the ions align. An extra DC field is applied to lift the degeneracy of the two radial modes, which can then be addressed distinctly. The internal state of one of the ions can be addressed by external lasers, initiating the JC interaction and a pulse on the blue Rabi sideband can prepare the qubit in the ground state.
In system 1S, the largest entanglement generation occurs for r 1S = 0.18. In Ref. [42], the ratio of such a coupling was r 1S 0.32, which is well within the range of our results. However, since this is above the optimal value forn 1 , it might be more effective to place the thermal energy in the second oscillator, with a nonzeron 2 . This corresponds more closely to the gold line in panel (b) of Fig. S1. Alternatively, the ratio can be further decreased to the optimal values by increasing the interion distance |z 0 1 − z 0 2 |. This is mainly controlled by the axial frequency ω z . If this can be reduced while main- Asymmetrically increasing the thermal noise of one oscillator increases the entanglement up to a saturation point. Additionally, increasing the coupling r1S or r1M up to a critical value increases the entanglement (black lines). Notably, H1M has a critical value at the symmetric value r1M = 1. The effects of increasing the mixedness of the remaining subsystems for the critical values of r1S (c-d) and r1M (g) result in a decrease in entanglement generation. H1S generates relatively less entanglement (b) and requires larger coupling to achieve it. Additionally, it shows crossings in LN that break the pattern of monotonic increase (with regard to r1S) seen otherwise. Panel (d) shows the importance of the thermal noise directly interfacing with the nonlinearity since, in contrast to the typical entanglement dynamics, adding thermal noise to oscillator 1 (decreasing the temperature difference) increases the entanglement for smalln2.
taining the conditions required for the approximations of the model then the more effective regimes of entanglement generation may be achievable. To realize the initial temperature difference in the oscillator pair, Doppler and sideband cooling can be applied to cool the vibrational modes towards the ground states, reaching average phonon numbers as low asn = 0.08. For thermally induced entanglement, it is the temperature difference that matters, and hence cooling of one of the two vibrational modes is sufficient to institute the required initial state.
It is possible to cool only one of the modes since the ions are individually addressable. However, the cooling time is typically a substantial proportion of the evolution time.
For this reason, we recommend that the initial states of the ions are prepared with a large inter-ion distance, then brought together, so that the excitation transfer dynamics only begins once the ions are prepared. In our results we have optimised the evolution time which usually requires only a few Rabi oscillations to produce nontrivial correlations. Therefore the entanglement is generated on a timescale in which dissipative effects are typically negligible. Notwithstanding this condition we give some attention to important decoherence effects in the Discussion and the Supplementary Material. Alternatively, system 1M is implemented in Ref. [44] in order to study NOON states. In this example, a single 171 Yb + ion is confined in a similar trap and an internal hyperfine transition is addressed to couple it to the radial motion. In this experiment the ratio takes the value r 1M 1, replicating our optimal coupling conditions. Similar remarks regarding the coherence times and state preparation apply here.

C. Multi-qubit thermal enhancement of entanglement
The JC model is readily generalised to the Tavis-Cummings model [45,46], in which a set of identical two-level systems interact with a single oscillator mode, and without a rotating wave approximation generalises to the Dicke model [47,48]. Such systems exhibit phase transitions that amount to collective enhancement of coherent quantum phenomena such as lasing and superradiance [49]. Motivated by such examples we briefly extend our analysis to demonstrate the potential enhancement of thermally induced entanglement via the coherent action of multiple qubits. Since we are interested only in demonstrating the existence of such an enhancement and the number of possible configurations of the system increases combinatorially with the number of qubits, we restrict ourselves to an analysis of two qubits interacting with an oscillator pair. Physically such systems can be conceived of as addressing a second ion with a driving electromagnetic field, inducing the JC interaction between this ion and the vibrational mode involved in the Coulomb-induced excitation exchange. Such a setup can be envisioned as involving a second narrow-waist laser directed at the second ion, or a broad-waist laser that addresses both ions, simultaneously exciting the internal motion. We again assume a ground state of the qubit, and analyse several possible configurations of the ionoscillator interactions.
To distinguish qubit labels from oscillator labels we will use the lower-case Latin alphabet to label the qubits. In analysing the various configurations of qubits and oscillators we introduce a few more ratio parameters by defining r BS 2S = κBS 12 κJC a1 for system 2S, and and r a2 2M = κJC a2 κJC a1 for system 2M (Fig. S2). Again, we study the dynamics in units of the scaled time τ = κ JCa1 t. In Fig. S2 we show a schematic of Hamiltonians 2S and 2M, defined by H 2S is symmetric in the qubit labels and reduces to H 1S when either κ JCa1 = 0 or κ JC b1 = 0, while H 2M , due to the qubit coupling symmetry, reduces to that of H 1M if either κ JCa1 = κ JCa2 = 0 or κ JC b1 = κ JC b2 = 0. In Fig. S2 the black curves indicate the cases in which these systems reduce to their relevant single qubit cases. Therefore, any parameter regime which surpasses the black curve indicates an enhancement over the singlequbit case. For system 2S, panels (a) and (b) indicate that enhancement occurs for the symmetric coupling of the qubits to the oscillator (red curve). In the previous section we noted the importance of coupling symmetry, alongside nonlinearity and temperature difference, for inducing thermal entanglement. This is strongly reflected here. However, we also see a striking difference from system 1S in that the significance of the temperature difference is inverted, so that the thermally excited distant oscillator is more significant for entanglement generation. Furthermore, crossing of the entanglement is visible and shows that there must be a sufficiently high temperature difference in order to gain enhancement. Panel (c) also demonstrates enhancement of entanglement for system 2M. As expected, the symmetry of the system means that the behaviour with respect ton 1 andn 2 is identical. The same requirement of symmetry occurs in order to achieve an enhanced generation of entanglement, and the same crossing occurs but with the requirement that the temperature difference is higher.
Finally we remark that a third configuration is given by H 1S−1S , where each qubit separately interacts with each oscillator. We write Similar to H 2S , H 1S−1S also reduces to H 1S when either κ JCa1 = 0 or κ JC b1 = 0. Our analysis reveals that H 1S−1S fails to produce any enhancement of entanglement than what is already present in H 1S in most parameter regimes.

III. DISCUSSION
We have demonstrated the autonomous generation of entanglement in a pair of linearly interacting oscillators, enhanced by the temperature difference of the oscillators and driven by the anharmonic modulation of their motion via discrete variable (DV) systems. This effect is of a different nature than previous examples of thermally driven entanglement generation between DV systems, driven either by discrete or continuous variable (CV) systems. The entangling of the oscillators is a nontrivial outworking of the nonlinear system as evidenced by the emergence of large-dimensional entanglement only after a significant interaction time. While the temperature difference drives the entanglement generation, this eventually saturates. A two-qubit maximally entangled state has LN = 1 and a pair of continuous variable systems are capable of storing much more than that [50], however these require additional resources beyond thermal energy. With the oscillators initially in their respec-  . Either both qubits are directly coupled to one of the oscillator pair, or both qubits together mediate the interaction between the oscillator pair. Each curve is calculated for an optimal time as described in the text. r BS 2S and r a2 2M are fixed at unity while r b 2S and r bi 2M are varied. For simplicity, we maintain the symmetry of r b1 2M = r b2 2M . (a-b) The results show the effect of increasing the thermal energy in one or other of the oscillators on the entanglement generation for H2S. (c) The symmetry of system 2M means that the role ofn1 andn2 can be interchanged for specific ratio parameters. The black curves denote the degenerate cases when the systems 2S and 2M recapitulate systems 1S and 1M.
tive ground states, a population inverted qubit (p e = 1) indeed maximally generates 1 ebit in the oscillators, suggesting that the entanglement of the system is saturated by the capacity of the qubit. With a thermal qubit (0 ≤ p e ≤ 1 2 ) and ground state oscillators however this is limited to LN ≈ 0.27, which is much lower. If the thermal resources are associated with the oscillators however, we find in at least one case (see Fig. S1 (a)) that LN = 0.5 which is almost twice that available from a thermally excited qubit. That said, thermal resources can coherently push beyond a thermal equilibrium state [51]. A thermal oscillator interacting with a ground state qubit via the JC coupling can induce partial population inversion with p e ≈ 0.63. A qubit prepared in such a way, beyond thermal equilibrium, can subsequently be coupled to a ground state oscillator pair and induce entanglement with LN ≈ 0.44, larger than LN ≈ 0.27 achievable for incoherent thermal equilibrium. Whereas LN ≈ 0.27 can be overcome by both 1S and 1M, LN = 0.44 is overcome only by 1S (see the Supplementary Material for a detailed analysis). This demonstrates the superiority of a scheme based on the thermal states of the oscillators and a CV beamsplitter class of dynamics.
Since thermal resources are exhausted, what remains is that the nonlinear resource of the qubit is over capacity. This is supported by the enhancement effects of multiqubit systems, which only take effect once the thermal noise is sufficiently high. However, other nonlinear interactions with oscillators may have even stronger potential for thermally induced nonlinearities. We note in general that symmetry in interactions, where available, supports the generation of entanglement. This contrasts with the asymmetry required in the thermal energy of the oscillators. Our demonstration of thermally induced entanglement, spread over a large dimension in the twooscillator Hilbert space has remained unnoticed for a long time, despite being a natural step from previous analyses. A reason might be that the natural method in platforms such as trapped ions is to dynamically transfer the entanglement to the qubits where readout methods are more developed [52]. We investigated the possibility of using a general qubit-oscillator interaction, to transfer entanglement from the oscillator pair to a qubit pair. This appears to fail, consistent with older results [53], and suggests that the large-dimensional entanglement induced by the qubits is curiously incompatible with qubit readout.
The coherence of the qubit interaction with the oscillator is absolutely necessary for entanglement generation, which we have verified in the Supplementary Material by analysing decoherence effects, particularly dephasing of the qubit. Large qubit dephasing is detrimental for entanglement generation. For small dephasing, however, the behaviour is qualitatively the same but with the LN reduced. In general the decoherence rates for trapped ion systems are below the timescales required to generate entanglement, so that while quantum features are typically sensitive to such processes they are expected to have negligible effects in our scheme.
Using thermal states to induce entanglement indicates that an equilibrium probability distribution over the Fock basis has the possibility to induce quantum effects in linear oscillators when coupled to an ancillary ground-state nonlinear system. While inducing entanglement by thermal states is the main result of this article, we have also investigated coherent states and phase randomized coherent states. Both of them will not produce any entanglement by a linear passive coupling in the oscillator, but they are less noisy than thermal states. Taking coherent states as examples, we find that the generation of entanglement is greatly enhanced when using such resources as demonstrated in the Supplementary Material. Similarly, phase randomised coherent states (PRCS), which implement a classical Poisson distribution over the Fock states rather than a Boltzmann distribution, also produce larger entanglement than thermal states although still less than coherent states. Both the PRCS and the thermal states are mixtures over the Fock basis and cannot produce entanglement with only linear coupling. For the same mean phonon number, Poissonian statistics perform significantly better than Boltzmann statistics due to less noise in the phonon number. As mentioned, if autonomous quantum systems are to be achieved both entanglement and the coherent states must be available unconditionally and autonomously. If quantum coherence can enhance entanglement generation as we have seen with coherent states, then thermally induced quantum coherence, if possible, will be a particularly powerful tool for autonomous quantum systems. However, that is a separate problem for the oscillator equivalent to the entire problem of shot-noise-limited lasing without an external drive.
The enhancement of the entanglement generation by complex DV systems bears some relation to existing generalisations of the JC interaction such as the Tavis-Cummings model or the Dicke model. Increasing the number of two-level systems is one method to address increasing the capacity of the nonlinearity. Another is to increase the dimension of the DV system by addressing more than two energy levels in the motion induced by the dipole interaction. The scaling of the entanglement dynamics with multilevel subsystems may unveil the details of how hybrid quantum systems eventually build deterministic nonlinearities for the oscillator. Adding more ions also adds more sets of passively coupled oscillators, allowing for the possibility to autonomously generate multipartite entanglement beyond Gaussian approximations [54]. This possibility is a natural further extension of our investigation and bears strongly on the potential existence of new autonomous quantum effects in many-body systems. Alternatively, the procedure developed here can be repeated multiple times with a newly prepared ground state qubit in order to accumulate entanglement. Such a technique has already been successfully used to accumulate thermally induced nonclassical features of a single oscillator [55].
When an experiment becomes available to demonstrate these phenomena, it will stimulate further development of thermally induced entanglement of oscillators on a larger scale. In support of that end, a small quantum simulator of such effects will unlock a deeper understanding of these phenomena and the techniques required to certify their feasibility. Any form of entanglement is, by definition, conclusive proof that two systems interacted in the past in a genuinely quantum mechanical way, outside the class of local operations and classical communication (LOCC) [56]. The growth in entanglement from separable states then provides evidence that a more extensive exchange of quanta occurred than that which occurs between thermal linear oscillators. That is, a genuine quantum mechanical interaction must have taken place, such as coupling to two-level systems. This hysteresis is essential for developing platforms for quantum technology and can be simulated and investigated using a small-scale hybrid quantum simulator. The results of such an investigation will aid new platforms involving hybrid quantum physics [57][58][59][60][61][62][63][64] to navigate through the discoveries of quantum interactions between oscillators.

I. DERIVATION OF THE RADIAL HAMILTONIAN
In a typical rf or Paul trap, the trapping potential is taken to be harmonic in all three spatial directions (α = x, y, z) with respective secular frequencies ω x , ω y and ω z . Therefore, for a single ion with mass m and charge e, the effective trapping potential is V tr = α 1 2 mω 2 α α 2 . For a pair of trapped ions with equal mass and charge, the effective potential energy V has an extra Coulombic contribution V C , i.e., V = V tr + V C , where where r i = (x i , y i , z i ) with i = 1, 2. Since each ion is identical and the trap is assumed to be homogeneous, the secular frequencies are also identical for each ion. Since there are two ions, we can simply choose the z-axis to be along the line connecting the ions and solve for the z-axis equilibrium points by setting ∂V ∂zi = 0. In general, if the trapping potential is arranged such that ω x , ω y ω z , then the ions will distribute along the z-axis, with equilibrium positions V C can be rewritten in terms of the displacements from equilibrium and the motion momentarily collected into relative motions X = x 1 − x 2 , Y = y 1 − y 2 and Z = z 1 − z 2 . The equilibrium position of Z is inherited from the local coordinates as Z 0 = z 0 1 − z 0 2 . Then V C is expressed as where the equilibrium positions of the X and Y relative motions are zero by construction. Note that the harmonic terms from V tr simply decompose into oscillations of the center-of-mass and relative motion coordinates and do not interfere with the following harmonic expansion of V C . This expression can be expanded around equilibrium up to harmonic terms, where the zeroth order term is a constant irrelevant to the dynamics and the first order term is zero due to the equilibrium condition. The second order term is Evaluating this produces the expression Thus the motion along the three axes is decoupled and the Hamiltonian for the motion along x direction is which reproduces Eq. (2) of the main text.

II. NOTE ON QUBIT-OSCILLATOR ENTANGLEMENT
In the main text, for various model systems we discussed at length the behavior of oscillator-oscillator entanglement, characterised by logarithmic negativity (LN). In this section, we give details of our investigation on other bipartite entanglement between two subsystems that constitute each model. For instance, the qubit-oscillator entanglement for system 1S (Fig. S1, left panel) roughly follows heuristic explanations of the entanglement generation. First, the oscillator with the noise becomes entangled with the qubit, then the nonclassical state of the oscillator induces entanglement between the oscillators through a linear interaction. Finally, the second oscillator becomes entangled with the atom which decreases the entanglement of the oscillators. There are residual correlations between the oscillators and the qubit, but they do not seem to be as significant as those created between the oscillators. System 1M is less easily explained in these terms (Fig. S1, right panel). However, since the qubit and one of the oscillators begin in the ground state, the first dynamics must be the excitation of the qubit by the noisy oscillator. The qubit can then drive nonclassicality in the oscillator pair, with it first becoming visible in the logarithmic negativity between the qubit and the less noisy oscillator. Subsequently the two oscillators become entangled, mediated by the qubit. In the multiqubit case the qubit-qubit and qubit-oscillator entanglement is even less significant compared to that generated between the oscillators (Fig. S2, left panel). , atom-atom (blue), atom 1-oscillator 1 (red), and atom 1-oscillator 2 (green), atom 2-oscillator 1 (yellow), and atom 2-oscillator 2 (pink) for both systems 2S (left) and 2M (right). The qubits in both the systems are initially kept the ground state, whilen1 = 2,n2 = 0. Red and yellow produce the same dynamics. These can be directly traced back to the symmetry of the systems at the critical value. Similarly, green and pink. For system 1S atom-atom entanglement is extremely negligible, mostly zero.

III. DECOHERENCE EFFECTS
In examining thermally induced entanglement among the vibrational modes of trapped ions we have assumed unitary dynamics. However small, there are still decoherence effects present in such systems. In this Section we give details of our numerical observations on the robustness of the entanglement preparation against such effects. In particular, we study the changes in the system dynamics due to the presence of qubit dephasing and other dissipative agents of open systems. Ignoring the contributions leading to a small renormalization of the system energy levels through the Lamb shift, the Lindblad master equation describing such environmental effects is given by Here ρ S (t) is the reduced density matrix of the system. The operators L k = √ λ k A k are the Lindblad jump operators, and the environment couples to the system through the operators A k with coupling rates λ k . We focus solely on H 1S here as the conclusions drawn are generic to the remaining systems.
We first regard the environment as a zero-temperature bath. The various sources of decoherence and their associated Lindblad jump operators are qubit dephasing A 1 = σ z , qubit decoherence A 2 = σ − , and mechanical decoherence A 3 = a 1 and A 4 = a 2 . Each have their associated decoherence rates, with λ 3 = λ 4 . The top panels of Fig. S3 display the inhibition of entanglement generation due to qubit dephasing for two different values ofn 1 and for various dephasing rates γ 1 . The effect of qubit relaxation alone on the dynamics is, however, less than that of dephasing as can be seen by comparing Figs. S3(b) and (c). Fig. S3(d) shows the cumulative effect of qubit dephasing and relaxation. This suggests that the coherence of the qubit is an important factor in the entanglement dynamics. On the other hand Fig. S3(e) corresponds to the dissipative effect arising due to both the oscillators interacting with the environment. Finally, in Fig. S3(f) we see the effect of all these environmental couplings. Now we consider a non-zero temperature thermal environment. For this relatively simple analysis we set the average number of phonons,n th , to be the same for each dissipation channel. For relaxation and mechanical dissipation effects the decoherence rates are modified to the pairs λ 2 (1 +n th ) and √ λ 2nth , and λ i (1 +n th ) and √ λ inth for i = 3, 4. The effect of increasingn th on the dynamics is displayed in Fig. S4 for various parameter regimes.

IV. COMPARISON: ENTANGLEMENT INDUCED BY COHERENT STATES
An alternative to driving entanglement generation in passively coupled oscillators using classical correlations (thermal states) is to consider quantum coherence, i.e. quantum correlations internal to a vibrational mode. As mentioned, two of the resources for quantum technologies are entanglement and coherence, and hence it is relevant for autonomous generation of these resources to examine how they intersect. As an example, coherent states which are distributed differently in the Fock basis, show quantum coherence, and are pure states. Here we briefly consider the entanglement generation when the oscillator pair consists of a single coherent state and the ground state.
The entanglement induced by a coherent state (CS) |α (expressed in the Fock basis as n c n |n where c n = e − 1 2 |α| 2 α n / √ n! and α ∈ C) of the oscillator is shown in Fig S5. Again, the dynamics is unitary and the qubit is initially in the ground state. It is clear that the entanglement thus generated is significantly higher than that induced thermally for both H 1S and H 1M . The critical value of the coupling r 1S is lower, although if r 1S is increased far away from the critical value the entanglement shows a decreasing trend with |α| 2 , after an initial increase. Once again, for H 1M the symmetry in the system is borne out in the dynamics and maximum LN is achieved for r 1M = 1. However, the entanglement generation is not monotonous with |α| 2 .

V. COMPARISON: ENTANGLEMENT INDUCED BY POISSONIAN NOISE
To distinguish between the quantum coherence of the coherent states and merely the change in distribution from Boltzmann to Poisson we also examine the entanglement induced via phase randomised coherent states. When the phase of a coherent state is unknown it becomes a mixed state and the photon number follows a Poisson distribution with a mean of |α| 2 . Mathematically phase randomised coherent states are defined by Therefore the number states are incoherent with each other.
The entanglement induced by such a state of the oscillator is shown in Fig S6. Again, the dynamics is unitary and the qubit is initially in the ground state. It is clear that the entanglement generated in this case is significantly higher than that induced thermally for both H 1S and H 1M but less than that of the coherent states. Again, the entanglement generation is not monotonous with |α| 2 . The qubit is prepared in the ground state. The entanglement is significantly higher compared to the case of an initially thermal state of the oscillator and less than that of standard coherent state. Left and right columns respectively correspond to oscillators 1 and 2 prepared as phase randomised coherent states.

VI. COMPARISON: THREE QUBIT CASE
In the main text we have considered a protocol wherein the passively coupled oscillators are prepared in thermal equilibrium, and this constitutes the resource (alongside coherent coupling with a two-level system) for inducing entanglement. Decreasing the purity of the qubit from the ground state negatively affects entanglement generation.
However, alternative strategies can be devised. We commented on them briefly in the discussion, and here we give a fuller treatment. One such alternative strategy is to consider a thermal equilibrium qubit and ground state oscillators. In this scenario there is a maximum of one excitation due to the nonzero probability of finding the qubit in the excited state. Thus, the full system may be well-approximated by a three-qubit system, where two of the qubits are prepared in their respective ground states and one in a thermal state. The models 1S and 1M in this approximation are equivalent and are described by the following Hamiltonian where numerical superscripts label the qubits. The label 1 always refers to the thermal qubit while 2 and 3 denote the proxies for the oscillator pair. For thermal qubits, max(p e ) = 1 2 and the maximum available entanglement generated among the remaining qubits is LN ≈ 0.27 (black curve in Fig. S7(a)). This contrasts with the nonequilibrium state p e = 1 which prepares maximally entangled pairs of qubits (LN = 1). Both system 1S and system 1M are sufficient to surpass the entanglement induced by a thermal qubit.
The possibility of approaching a maximally entangled state via p e > 1 2 , i.e. population inversion, motivates us to attempt to use thermal equilibrium states to generate such qubit states. By coupling a ground state qubit to a thermal oscillator, one may prepare a nonequilibrium qubit with p e ≈ 0.63, using only equilibrium resources (red curve in Fig. S7(b)). If this state is selected as the initial state in the three-qubit approximation, the maximum achievable entanglement is enhanced to LN ≈ 0.44. This is below the level of entanglement achievable with H 1S . With H 1M there is also the capacity to cross this threshold, but it occurs for the second entanglement peak, which we do not analyse.

VII. NOTE ON NUMERICAL METHODS
Our simulations take place in a truncated Hilbert space characterised by a Fock dimension n, which characterises the size of each subspace associated with a vibrational mode. To ensure stability and accuracy we adhere to an error threshold of 10 −6 , much higher than the precision typically available in experiment. Fig. S9 shows how the truncation dimension n scales with the error threshold. The error threshold parameter is referred to as the power scaling 10 − . As expected higher accuracy requires higher dimensions. In order to maintain this accuracy with increasingn 1 the dimension must rise quickly. The figure presented is generated for system I (with p e =n 2 = 0 and r 1S = 1), but the result is generic and holds qualitatively for all other systems.

VIII. IMPORTANCE OF THE GROUND STATE
In the main text we emphasised that the initial state of the qubit being the ground state is important for thermally induced entanglement. Here we expand on this by examining the state of the qubit throughout the entanglement generation process. The most clear-cut example occurs for the approximation in Section 5. The fidelity of the qubit with the ground state, alongside the evolution of the entanglement of the remaining pair of qubits is shown in Fig S8. For this unitary system, the behaviour is simple and oscillatory, showing collapses and revivals of the entanglement and of the excitation of the thermal qubit. Despite this simple behaviour, a full collapse and revival of the entanglement occurs while the thermal qubit stays close to the ground state, eluding simple explanations based on the transfer of a single excitation.
More complicated behaviour occurs for the main protocol of the article and is shown in Fig. S8. Here the correlation of the ground state of the qubit with the entangled states is much weaker. However, it is clear that the qubit undergoes some evolution before the entanglement emerges, reflecting the nontrivial dynamics generating the entanglement.

IX. EARLY ENTANGLEMENT DYNAMICS
In the main text we observe that the logarithmic negativity is zero for the early parts of the dynamics. LN is a sufficient condition of entanglement and therefore can be zero even when the state is entangled. For such largedimensional systems it is difficult to evaluate this region more clearly. However we can compare the entanglement dynamics of the LN with those of the concurrence (in the three qubits analogy, Eq. (S8)). The concurrence is shown in Fig. S10 against the LN. The concurrence has its major peaks synchronised with the LN, however for system 1M, if we were to follow our principle of taking the first peak of the entanglement dynamics for analysis of thermally induced effects we would miss the possibility of large entanglement that the LN captures.

X. ENTANGLEMENT POTENTIAL AND LN
The nonclassicality of a state ρ capable of being converted to entanglement can be characterised through its entanglement potential (EP), which characterises the state's capacity to become entangled via passive interactions with the ground state of an oscillator. This is quantified by the logarithmic negativity (LN), of the following state: where U BS = e π 4 (c † a−ca † ) and c is an auxiliary mode prepared in the ground state. The EP for a qubit coupled to a thermal oscillator is shown in red in Fig. S11. In contrast, the black curve shows the corresponding LN (see also Fig. 2 of the main text) when the thermal oscillator is simultaneously coupled to a ground state oscillator. The EP of the thermal oscillator increases before the LN of the full system, supporting our interpretation that first the qubit must drive the thermal oscillator to a nonclassical state before it can entangle with the second oscillator. However the LN eventually surpasses the EP. This suggests that the back action from the qubit continues to drive the nonclassicality of the thermal oscillator while it generates entanglement with the second oscillator.