Quantum Simulation of Non-perturbative Cavity QED with Trapped Ions

We discuss the simulation of non-perturbative cavity-QED effects using systems of trapped ions. Specifically, we address the implementation of extended Dicke models with both collective dipole-field and direct dipole-dipole interactions, which represent a minimal set of models for describing light-matter interactions in the ultrastrong and deep-strong coupling regime. We show that this approach can be used in state-of-the-art trapped ion setups to investigate excitation spectra or the transition between sub- and superradiant ground states, which are currently not accessible in any other physical system. Our analysis also reveals the intrinsic difficulty of accessing this non-perturbative regime with larger numbers of dipoles, which makes the simulation of many-dipole cavity QED a particularly challenging test case for future quantum simulation platforms.


I. INTRODUCTION
Quantum electrodynamics (QED) is our fundamental theory for describing the dynamics of charges coupled to the quantized electromagnetic field. In contrast to quantum chromodynamics for the strong force, QED is a weakly-interacting theory, which is characterized by the small value of the electromagnetic finestructure constant, α fs 1/137. Apart from its implications in particle physics, this property affects as well many processes relevant in our daily life, for example, the way in which light interacts with atoms, molecules and solid matter. Specifically, the smallness of α fs implies that the coupling strength g between a single elementary dipole and a single photon of frequency ω c is constrained to g/ω c √ 2πα fs 1. [1][2][3] As a consequence, the coupling between matter and photons can typically be treated as a small perturbation on top of the absolute energy scales and does not considerably alter the overall structure of ground-and excited states.
In recent years there has been a growing interest in the physics of light-matter interactions beyond this conventional coupling regime. [4,5] In many experiments with dense ensembles of electrons, excitons or molecules it is now possible to reach ultrastrong-coupling (USC) conditions, [6] where the collective coupling, G = √ N g, between N 1 dipoles and a single cavity mode reaches a considerable fraction of the bare photon frequency. Moreover, in the field of circuit QED, [7][8][9] analogue models for light-matter interactions can be implemented by coupling artificial atoms, i.e., superconducting two-level systems, with microwave photons. In this case the bound mentioned above can be overcome by using high-impedance resonators or galvanic coupling schemes, such that g/ω c 1 can be realized even with a single qubit. [10,11] In this, often called deep-strong-coupling (DSC), regime, [12] the interaction between atoms and * tuomas.jaako@tuwien.ac.at photons has a non-perturbative effect on the energy-level structure of the combined system.
Despite a considerable amount of work on this subject, non-perturbative effects in cavity and circuit QED are still little understood. In the past, a lot of theoretical studies in this field have been devoted to the quantum Rabi model and variations thereof. [4] This model, however, only describes the coupling of a single dipole to a cavity mode and does not capture cavity-mediated interaction effects. In turn, collective interactions and phase transitions are usually discussed in the opposite limit of a large number of dipoles, N 1. In this case the relevant coupling parameter per atom, g/ω c , is typically assumed to be small, such that most effects can be understood in terms of conventional electrodynamics. [3] Thus, the most intriguing regime, where both non-perturbative and many-body effects play a role, remains largely unexplored, which is also related to the fact that the combined conditions g/ω c 1 and N > 1 have not been demonstrated in any of the mentioned experimental platforms so far. This motivates the search for alternative quantum simulation schemes, where the USC physics can be explored, independently of any technological constraints and under fully controlled conditions. [13][14][15][16][17][18][19] In this work we investigate the use of systems of trapped ions as a quantum simulator for multi-dipole cavity QED systems in the USC regime. This platform is naturally suited for this purpose since experimental techniques for implementing Jaynes-Cummings-, Rabiand Dicke-type couplings between the internal atomic states and motional modes (which represent the photons in the effective model) are already well-established. [20][21][22][23] However, for N > 1 such models provide very restricted or inconsistent descriptions of quantum electrodynamics beyond the weak-coupling regime and must be complemented by additional interaction terms. [3] This includes all-to-all and short-range dipole-dipole interactions, which account for the "P 2 -term" [24][25][26] and electrostatic forces, respectively. Both contributions play a dominant role for cavity QED systems in the regime of very large coupling. Here we describe how such extended Dicke models can be simulated in a chain of trapped ions Note that compared to Ref. [3], here an independent scaling of J0 and g is assumed.
with a large tunability of the effective model parameters.
Specifically, we demonstrate that characteristic USC effects in the excitation spectra and in the ground state of this system can be probed with N ≈ 2 − 10 ions in stateof-the-art experimental setups. Therefore, small-scale trapped-ion quantum simulators can already be used to explore the few-dipole USC regime of cavity QED, which is not accessible in any other physical platform today.
II. ULTRASTRONG-COUPLING CAVITY QED Figure 1(a) depicts a generic cavity QED setup, where N two-level dipoles with transition frequency ω 0 are coupled to a single electromagnetic mode of frequency ω c . Under the assumption that the electric field of the dynamical mode is sufficiently homogeneous, such a scenario is described by the Hamiltonian ( = 1) [3] H cQED =ω c a † a + ω 0 S z + g(a † + a)S x Here the σ α i , where α = x, y, z, are the Pauli-operators for the i-th dipole, S α = 1/2 i σ α i are the corresponding collective spin operators and a (a † ) is the annihilation (creation) operator for the electromagnetic mode. The first line in Eq. (1) is the usual Dicke model, [27] where g denotes the coupling strength of a single dipole. This model describes well the collective interaction between N dipoles and a common field mode in the weak coupling regime, G = g √ N ω c . At larger coupling strengths, the two additional spinspin coupling terms in the second line of in Eq. (1) must be taken into account. The first contribution ∼ S 2 x is the so-called depolarization or P 2 -term. Its origin is related to the fact that H cQED is derived in the dipole gauge, [3,26,28] where the canonical momentum variable of the electromagnetic mode is the displacement field, D = ε 0 E + P , [24] with E being the electric field and P ∼ S x the polarization density. Therefore, when expanding the electric field energy, ∼ E 2 ∼ (D − P ) 2 , we obtain both the dipole-field interaction together with the accompanying S 2 x -term, which therefore should be interpreted as part of the field energy. In contrast, the last term in Eq. (1) accounts for the actual dipoledipole interactions, which exist independently of the dynamical mode. In free space we would simply obtain J ij ∼ 1/|r i − r j | 3 , where r i are the positions of the dipoles. However, in the presence of metallic boundaries, screening effects, etc., the actual dependence may be substantially modified, [3] and can also be engineered to be infinite-ranged in circuit QED systems. [29,30] Equation (1) shows that even at a minimal level, models of cavity QED involve collective interactions between spins and a bosonic mode as well as direct spin-spin interactions with different spatial dependencies. These terms are not completely independent of each other and in particular the strength of the P 2 -term must match the dipole-field coupling to ensure consistency with basic electrodynamics. The static dipole-dipole interactions depend on the system and geometry under consideration and will in general introduce short-range interactions, which compete with the collective dipole-field coupling. Depending on the ratio of g/ω c and the sign and strength of the couplings J ij , different normal (paraelectric), superradiant (ferroelectric) and subradiant (antiferroelectric) phases can occur. This behaviour is illustrated in Figure 1(b) for the simplified case J ij = J 0 /N , where due to symmetry the model can still be solved numerically for small and moderate numbers of dipoles. [3] For short-range interactions, and depending on the geometry, different other types of phases may exists, but due to its computational complexity, little is still known about the ground and excited states of H cQED in such general scenarios.

III. EFFECTIVE CAVITY QED MODELS WITH TRAPPED IONS
For the implementation of H cQED as an effective, but fully controllable, model we consider a system of N trapped ions in a linear Paul trap as shown in Figure  2(a). At low enough temperatures the ions will arrange themselves in a one-dimensional (1D) chain and by writing the position of the i-th ion as r i = r 0 i + u i , we can linearize the motional dynamics around the equilibrium positions r 0 i . The residual displacements u i can be quantized and written in terms of a set of bosonic annihilation (creation) operators b α,n (b † α,n ) as Here M is the mass of the ions and ν α,n and ξ α n denote the frequency and mode function of the n-th vibrational eigenmode in direction α, respectively. A typical phonon spectrum is shown in Figure 2(b) for the case ν z < ν x < ν y , where the ν α are the trapping frequencies along the three principal axes.
The ions are driven by two pairs of laser beams, which are slightly tuned to the red (r) and the blue (b) of the transition frequency ω eg between the long-lived electronic states |g and |e . As indicated in Figure 2(a), one pair of lasers is directed along the x-axis and the other pair along the y-axis and we assume that the beams are sufficiently broad such that they can be treated as plane waves. The Hamiltonian of the whole ion chain is then given by where ω α,l and Ω α,l are the frequency and the Rabifrequency of laser l ∈ {r, b} with wavevector k α,l = k α e α . Since in the considered configuration the motion along the z-axis remains unaffected, we can restrict α ∈ {x, y} and assume e ik α,l ·r 0 i 1.

A. Interaction engineering
In the Lamb-Dicke regime, where the residual motion of the ions is small compared to the laser wavelength, we can expand the exponentials in Eq. (3) up to first order in the displacements u i . [20] Under this approximation and by changing to a frame rotating with ω eg , the relevant laser-induced coupling between internal and external de-grees of freedom reduces to where As long as the couplings η α n Ω α,l are small compared to the spacing between the vibrational modes, the detunings δ α,l can be chosen to resonantly enhance the interaction with a specific phonon mode, while the coupling to other modes as well as direct transitions between the internal states are strongly suppressed. This general scheme is frequently used in trapped ion systems to engineer different types of spin-phonon and spin-spin interactions and references to some of the relevant previous works in this field will be given in the following discussion. However, the cavity QED Hamiltonian (1) involves several different types of interactions with a finetuned relation between the coupling parameters. Therefore, here our goal is to show how these general techniques can be combined to engineer H cQED with a large degree of control over all the parameters.

Collective dipole-field coupling
We first select one of the vibrational modes to represent the photonic mode in the cavity QED model. Here we choose the transverse center-of-mass (COM) mode along the x-axis, which has a homogeneous mode profile and a frequency ν x,COM = ν x . For the two laser beams along the x-direction we assume equal amplitudes, Ω x,l = Ω x , and detunings δ In the interaction picture with respect to the phonon modes and keeping only near-resonant terms we then obtain and after eliminating the time-dependence via the unitary transformation U (t) = e −i(ωca † a+ω0Sz)t we obtain the effective Hamiltonian This and closely related schemes have already been discussed in many previous works for implementing effective Rabi-and Dicke models, [13][14][15][16][17][21][22][23] with the crucial benefit that the ratio between g, ω c and ω 0 is fully controlled by laser or microwave detunings, rather than by the bare physical parameters.

The P 2 -term
The remaining terms in H cQED contain direct spin-spin interactions ∼ σ x i σ x j , with both constant and spatially varying prefactors. To implement the collective P 2 -term we use the two lasers along the y-direction to address the COM mode b y,COM with ν y,COM = ν y . By assuming equal Rabi frequencies Ω y, = Ω y and the detunings δ y,r/b = ∓µ+ω 0 we obtain the interaction-picture Hamiltonian where g y n = Ω y η y n . To realize pure spin-spin interactions, we consider the regime where g y n µ − ν y,n , such that the phonon modes are only virtually populated. The dynamics of the spins are then controlled by the spin Hamiltonian [31][32][33][34][35][36] When the beat-note is chosen to be close to the COM motional frequency, µ = ν y + ∆ y , with |∆ y | being small compared to the mode spacing, the dipole couplings are approximately constant, This value can then be tuned to match g 2 /ω c in order to reproduce the correct P 2 -term. Note that due to the presence of the local field ω 0 in Eq. (8) terms proportional to b † y,n b y,n σ z i will be generated, which, however, are suppressed by a factor of ω 0 /(µ − ν y,n ) and can be neglected in the parameter regime of interest. [37]

Dipole-dipole interactions
To implement additional short-range dipole-dipole interactions, we generalize the scheme from above to laser beams with multiple modulation sidebands with slightly different frequencies. This can be accounted for by substituting in Eq. (4) Ω y e iδ y,l t → m Ω y,m e iδ y,l,m t .
The resulting Hamiltonian is that of Eq. (8) with a sum over the different modulation sidebands. As long as the detunings between different modulation frequencies remain large, i.e., η y n Ω m |δ m − δ l |, the resulting cross terms from lasers with different beat-note frequencies δ m are rapidly oscillating and can be neglected. The dynamics of the spins is then determined by an effective Hamiltonian as in Eq. (9), where the generalized interaction matrix [31][32][33][34][35][36] can be engineered in a flexible manner by combining multiple near-resonant and/or far-detuned lasers. For example, by adding a laser which is far detuned from all modes, the coupling matrix will acquire an additional component which decays approximately as ∼ |i − j| −3 in the limit of very large detuning. [31,33] In contrast, when addressing one of the modes in the middle of the phonon band one obtains a coupling with an alternating sign, which typically leads to frustration and related phenomena. [34] B. Non-uniform couplings The expression for H cQED as given in Eq. (1) is based on the usual assumption that the field profile is homogeneous over the extent of the ensemble of dipoles. However, the model can easily be generalized to situations where the coupling strength g i is different for each dipole, by making also the corresponding substitution for the P 2 term, The interaction engineering schemes discussed above allow the implementation of such non-uniform models by considering non-uniform mode profiles for the driving lasers, Ω x/y → Ω x/y (r 0 i ). Since in our scheme we use the nearly uniform COM modes for engineering both the dipole-field coupling and the P 2 -term, the correct relation between the two terms is automatically guaranteed when the same mode profiles for the lasers along the x and the y direction are assumed. Therefore, with the use of spatial light modulators or other experimental techniques, arbitrary mode profiles g i can be engineered while still retaining physically consistent models. However, for concreteness we will focus on the homogeneous case, g i = g, in the remainder of the discussion.

C. Accessible parameter regimes
In summary, by addressing different phonon branches of the ion chain, both collective spin-photon and spin-spin interactions can be engineered independently.
where the D ij are given by Eq. (12). While in theory this approach provides full control over all relevant model parameters, the hierarchy of frequency scales and the singlemode addressability assumed in the derivation of the effective interactions still impose practical limitations on the accessible parameter regimes.

Ultrastrong coupling regime
As a specific example we consider a chain of N = 10 trapped 40 Ca + ions with a phonon spectrum as shown in Figure 2(b). In this case the relevant Lamb-Dicke parameter is η x COM = 0.043 and for Ω x = 2π × 15.4 kHz, ∆ x,b = 2π × 0.41 kHz, and ∆ x,r = 0 we obtain ω c = ω 0 = 2π × 0.21 kHz and a coupling parameter of g/ω c = 1. For the implementation of the P 2 -term we follow the scheme in Sec. III A 2 and use two lasers with Rabi-frequency Ω y,1 = 2π × 139 kHz to drive the COM mode with detuning ∆ y,1 = 2π × 14 kHz. For η y COM = 0.041 this results in a collective S 2 x -coupling of strength g 2 /ω c = D = 2π × 0.21 kHz. Since in a real trap the mode function is not completely homogeneous and the laser will also weakly couple to all other y-modes, the exact evaluation of the coupling matrix D ij in Eq. (12) will result in small spatial variations, D ij ∼ |i − j| −0.16 . The resulting residual dipole-dipole interactions, J ij , are plotted in Figure 3(a). On average, each dipole feels a residual fieldJ i = j =i J ij , with a variance (∆J) 2 = i (J i ) 2 /N across the chain. For the current set of parameters ∆J/ω 0 ≈ 0.19, and to a good approximation the dipoles can be considered non-interacting. Note that imperfections in D ij scale with ∼ g 2 and become negligible for weaker couplings.
To add additional short-range dipole-dipole interactions in a controlled manner, we now consider a second pair of lasers along the y-axis with strength Ω y,2 = 2π × 1.0 MHz and detuning ∆ y,2 = 2π × 1.7 MHz. In this far-detuned limit the coupling to all the phonon modes is roughly the same and the resulting dipole-dipole couplings scale approximately as J ij ≈ J 0 /|i − j| α , where α 1.98 and J 0 = 2π × 80 Hz, see Figure 3(d). In this case the total mean fieldJ ≈ 2π ×0.20 kHz is comparable to ω 0 . The exact coupling matrix J ij , including the residual imperfections from the P 2 -term, is shown in Figure  3(b). For this simple driving scheme, the value of α 2 is limited by the ratio between the phonon bandwidth and the detuning. For a 1D chain the interactions can be considered as mid-range. [38] However, since the second laser must be detuned far to the blue, ∆ y,2 > 0, the effective interactions are necessarily repulsive. To implement an equivalent model with attractive interactions, we can simply invert the sign of all the other terms in H cQED , [39] which can be done by replacing Ω x → −Ω x and changing the sign of the detunings, ∆ x,r/b and ∆ y,1 . As a result we obtain the model −H cQED with J 0 < 0. For the purpose of quantum simulation, the overall minus sign is unimportant. The plot in Figure 3(c) shows the resulting coupling matrix for ∆ y,1 = −2π × 11 kHz and Ω y,1 = 2π × 112 kHz, which, apart from the sign, leads essentially to the same effective parameters as above.
In summary, this example shows that trapped ions can be used to engineer few-body cavity QED models with coupling parameters g/ω c ∼ O(1) and absolute frequency scales of a few hundreds of Hz. This is still fast compared to simulation times of tens of milliseconds available in state-of-the-art trapped-ion experiments. [37,40,41]

Non-perturbative regime
In the previous example the collective coupling G = √ N g already exceeds the cavity frequency by a factor of three. In the recent literature, [4,5] this regime is very generally called the DSC regime, without distinguishing between the collective and the single-dipole coupling constant. However, as indicated in the phase diagram in Figure 1(b) and discussed in more detail in Ref. [3], significant non-perturbative changes in the physical properties of the cavity QED system are only expected beyond a value of g/ω c ≈ 2 − 3 of the single-dipole coupling parameter, approximately independent of N . For simulating this regime, two main difficulties arise. First of all, by assuming a fixed value of D = g 2 /ω c ≈ 2π × 200 Hz as above, the frequencies ω 0 ≈ ω c ≈ 2π × 20 Hz must be reduced by a factor of about ten to reach these high values of the coupling parameter. Second, the reduced value of ω 0 also means that any residual deviations of the actual coupling matrix, ∆J/ω 0 ∼ 2, have now a much stronger impact on the bare model. These problems will in general become worse for larger N , where the conditions for single-mode resolution become more stringent and lead to a competition between the time-scales of the simulation and the quality of the model, i.e. the level of control over the interaction matrix D ij . However, this is not a fundamental limitation and in particular for small and moderate numbers of ions, there are still many interesting effects that can be explored under those constraints.

IV. EXAMPLES
The parameters estimated above show that systems of trapped ions can be used to simulate otherwise unaccessible parameter regimes in cavity QED. In this section we discuss two basic examples, which also illustrate different measurement techniques that one can apply to extract interesting information about this system.

A. Few-dipole excitation spectrum
As a first example we consider the measurement of the excitation spectrum of a few-dipole cavity QED system in the parameter regime g/ω c 1. The USC regime of cavity QED can be identified in the cavity excitation spectrum, as the region where the splitting between the two polariton modes, ∆ω, starts to deviate from the initial linear scaling ∆ω G. In typical experiments in the optical and THz regime the condition G ∼ ω c is only accessible with a very large number of dipoles, where g/ω c 1 and only linearized collective excitations can be probed. [4,5] Corrections to this linear spectrum are expected to become observable for N 10, [26] but reaching this regime presents a notable experimental challenge. Superconducting circuits can more easily enter the USC regime, but the condition g ω c has so far only been achieved with single flux qubits, due to the complexity of controlling and measuring multiple such devices. Finally, in all natural cavity QED systems, the dipole-dipole interactions are usually fixed or difficult to control. All these limitations are absent in our trapped-ion quantum simulator.
To measure the few-dipole excitation spectrum, the ions are initialized in state |g and the photon mode is cooled to its ground state. Then, all the coupling terms are gradually increased from zero to their final value such that the system is adiabatically prepared in the ground state |G of H cQED . Finally, a weak perturbation of the form H p (t) ∼ Ae iωt +A † e −iωt is applied for a time T p . In the limit T p → ∞, the amount of excitations created by such a perturbation will be proportional to the excitation spectrum where the average is take over the actual state ρ 0 ≈ |G G| after the adiabatic preparation. In practice, a measurement of the expectation value A † A before and after the applied perturbation will already provide an accurate estimate of S(ω) for finite T p . For A ≡ a this procedure provides a measurement of the cavity spectrum. In the following we focus instead on the case A ≡ σ − 3 , where S(ω) also contains information about the so-called dark polariton states, which are excitations of the dipoles that are decoupled from the cavity mode. Figure 4 shows the numerically simulated result of such an experiment for the case of N = 6 ions but otherwise similar parameters as discussed in Sec. III C 1. Since from these simulations we found that a fully adi-abatic preparation of the ground state requires a too long time of hundreds of milliseconds, we use a nonadiabatic bang-bang scheme, similar to what has been used previously. [23,42,43] With this procedure detailed in App. A the ground state can be prepared with a fidelity of F = G| ρ 0 |G 0.8 in a time T prep 7 ms. Starting from this state, we use Eq. (15) to evaluate the excitation spectrum for different frequencies ω, where we simply assume a common phenomenological decay rate of Γ = 2π × 4 Hz for all excited states, corresponding to an experimental runtime of T p ≈ 40 ms.
In Figure 4(a) the resulting spectrum is first plotted for non-interacting dipoles, where up to the residual imperfections described in Figure 3, J ij ≈ 0. For small G one observes the expected Rabi-splitting ∆ω G between the two bright polariton states while other excitations of the dipoles are decoupled from the cavity and remain almost unaffected. Note that since we consider the response of a single dipole, the signal is sensitive to all excitation modes, but the overlap with the collective polariton modes is reduced by a factor 1/N . Thus, the ability to see both collective and single-particle effects is a specifically interesting feature of the considered few-dipole regime. At large couplings the influence of the P 2 -term is no longer negligible and the spectrum starts to deviate from the predictions of the usual Dicke model. [6,25,[44][45][46] In particular, the frequency of the lower polariton mode stabilizes at a non-zero value for all couplings. [6,25] A somewhat unexpected observation is the downward shift of the dark polariton modes, which is not predicted by a purely linear theory. It arises from the fact that the ground state energy E G (g) increases with increasing g. When one of the dipoles is now promoted to a decoupled mode, less energy is needed. Finally, a finite splitting between the dark modes indicates residual dipole-dipole interactions J ij due to a nonuniform matrix D ij . As a next step we switch to a system with strong attractive dipole-dipole interactions, J ij ≈ J 0 /|i − j| α and J 0 < 0. In this case the dipoles can undergo a transition into a ferroelectric state at a critical coupling J c 0 . For an infinite system and α = 2 a value of J c 0 /ω 0 −0.4 is predicted. [47] In the presence of the cavity mode this value is expected to decrease as due to the dressing of the dipoles by virtual photons. [3] Of course, for the considered small number of dipoles, N = 6, there is only a smooth crossover between the normal and the ferroelectric phase. However, the two phases can still be distinguished by looking at the probability distribution p(m x ) = G| P mx |G , where P mx = s P s,mx and P s,mx is the projector on states with S x |ψ = m x |ψ and total spin s. The phase boundary can then be defined as the line, where this function changes from a single to a bi-modal distribution. In Figure 4(b) this boundary is sketched for the current model parameters and for different values of the light-matter coupling g. We see that within the accessible parameter range the transition line can be crossed in two different ways: Either in the conventional sense, by increasing |J 0 |, or by keeping |J 0 | < |J c 0 | fixed, but varying the coupling to the cavity. The corresponding spectra are shown in Figure  4(d) and (c). We see that in both cases the frequency of the lowest excited mode goes to zero around the expected transition point, where one must keep in mind that for a finite system the energy gap at J c 0 (g) remains finite. The difference between the two tuning schemes is clearly visible in the excited states. In the first case, avoided crossings around ω c indicate strong hybridization with the cavity mode, while the same features are absent when only J 0 is varied.
We emphasize that the strong reduction of the lower polariton frequency observed in Figure 4(c), which indicates the transition from the normal to the ferroelectric phase, is a true non-perturbative effect. It arises from to the renormalization of the dipole frequency ω 0 , while the strength of static dipole-dipole interactions remains fixed. This is in contrast to the celebrated superradiant phase transition in the Dicke model. [27] From Eq. (1) we see that when interpreted in the context of cavity QED, the Dicke model corresponds to the case of a ferroelectric ensemble of dipoles with J ij = −g 2 /ω c . As a consequence, when increasing g in this model, the system undergoes a regular ferroelectric phase transition, where the coupling to the cavity mode only introduces minor modifications. [3] B. Subradiant ground state While already in the regime g/ω c 1 first nonperturbative corrections are observable, the ground state is still determined primarily by the competition between ω 0 and J ij and the influence of the cavity mode is mainly seen in the excited states. This changes drastically in the regime g/ω c > 2, where apart from the renormalization of the transition frequency, the cavity also induces effective anti-ferroelectric interactions, [29] For small |J 0 |, these effective interactions compete with the short-range couplings and favor so-called subradiant ground states with completely anti-aligned dipoles that are decoupled from the cavity mode.
The simplest experimental setting where this effect can be explored is the case of N = 2 ions. In this case the matrices D ij = D and J ij = J 0 have only one relevant entry, which allows us to relax some of the detuning constraints and consider values of D = 2π × 2 kHz and corresponding values of ω c,0 ≈ 2π × 100 − 500 Hz to access coupling parameters up to g/ω c = 4. For these parameters we and ω0 = ωc. The ground state is prepared by adiabatically turning on the dipole-field coupling g and dipole-dipole interactions J0 on a timescale of Tprep 10 ms, starting from the ground state of the noninteracting system. plot in Figure 5(a) and (b) the expectation value of the photon number a † a in the simulated ground state, ρ 0 , for J 0 = 0 and J 0 /ω c = −3.5. For small g we see in both plots the expected increase of the photon number due to a hybridization between the dipoles and the photons. For the ferroelectric case, the photon number then increases rapidly after g/ω c ≈ 1, which is the characteristic signature of a superradiant phase. In contrast, for non-interacting dipoles this trend turns around after g 2ω c and the cavity returns back to its ground state for very large couplings. While for the considered preparation time T prep ≈ 10 ms the simulated photon number still differs from that of the true ground state, the characteristic maximum, which is the key signature for entering a subradiant ground state, [3,29] is clearly visible. In a trapped ion quantum simulator this effect can be verified independently by performing a full tomography of the internal state. For the maximal coupling we find an overlap of Tr{ρ 0 |T T |} ≈ 0.99, where |T = (|ee − |gg )/ √ 2 is the maximally entangled state that minimizes H AF .

C. Validity of the effective model and decoherence
All the results presented in this sections are based on numerical simulations of the effective model, taking all imperfections of the coupling matrix J ij into account, but neglecting the weak admixture of other phonon modes. To ensure validity of the effective model, we have chosen parameters such that This means that the occupation of the COM mode, which is used to implement the P 2 -term should be at most a few percent, b † b < 0.1. We have explicitly verified this estimate by performing numerical simulations where the dynamics of the COM mode is included. In these simulations we do not see any significant changes in the dynamics of the system when compared to the effective one-mode model.
In real experiments the system will also be affected by decoherence of the internal states and heating of the phonon modes. In all the examples discussed above the time for preparing the ground state, T prep , is chosen to be at most 10 ms. In state-of-the-art ion traps the heating rates can be as low as 1-10 quanta per second, [48,49] meaning that the number of added phonons during the simulations is less than ten percent. In addition, for g/ω c > 1 the spin states start to decouple from the oscillator mode and, thus, the internal state is even less affected by heating. A remaining source of error is the dephasing of the internal states by magnetic field fluctuations or laser phase noise. For the adiabatic ground state preparation scheme, the ions are initially in an eigenstate of σ z and, therefore, dephasing is only relevant in the final part of the protocol. A master equation simulation, including also a heating rate of 10 quanta per second for the phonon mode, of the ground state preparation protocol of Figure 5 for g/ω c = 4 shows that for a dephasing time as low as T 2 = 10 ms, the fidelity of the final state only changes by ∆F = 0.15 and by only ∆F = 0.02 for realistic dephasing times of T 2 = 100 ms. [50] These findings are consistent with other quantum simulation experiments, where simulation times of > 50 ms have been demonstrated. [40,41] V. DISCUSSION AND CONCLUSIONS In summary, we have presented a comprehensive analysis on the suitability of ion traps for simulating the extended Dicke model (1), which captures the essential nonperturbative effects of cavity QED systems in the USC and DSC regime. Compared to real cavity or circuit QED systems, such simulators provide a flexible way to tune independently the coupling between the dipoles to a dynamical cavity mode and direct electrostatic interactions.
The two examples discussed in more detail in Sec. IV illustrate different possibilities for exploring characteristic signatures of non-perturbative light-matter interactions in the ground and excited states of few-body cavity QED systems.
The analysis of this work has been restricted to a small number of ions, where adiabatic and non-adiabatic ground-state preparation schemes can still be benchmarked by a comparison with exact numerics. Similar to many closely related proposals, [13][14][15][16][17][18][19] in this case the simulation aspect primarily lies in the ability to study coupling regimes that are fundamentally not accessible in systems of atomic or molecular dipoles. However, in principle, the same ideas can be generalized to several tens of ions or multiple cavity modes in order to explore non-perturbative effects far beyond the reach of classical simulation capabilities. From our numerical studies we find that the main practical difficulty in doing so arises from the collective P 2 -term, ∼ g 2 /ω c S 2 x , which becomes the dominant contribution in the regime g/ω c 1. This term is implemented by selectively addressing the center-of-mass mode, which becomes more and more difficult as the number of ions increases and significantly prolongs the experimental timescales. This feature makes the extended Dicke model a particular challenge for trapped-ion systems and other quantum-simulation platforms.
However, we envision that with improved motional heating and spin coherence times in future ion traps, simulation timescales of several seconds will become possible. It also has been shown [51,52] that with full single-site addressability, stroboscopic techniques and numerical optimization the design of the coupling matrix D ij can be considerably improved to reduce residual imperfections while retaining a high coupling strength. Therefore, with further experimental and theoretical work along these lines, also the simulation of non-perturbative effects in cavity QED systems with tens of dipoles and multiple modes is achievable, where currently neither analytic predictions nor numerical simulations are available.