Quantum simulation of the dynamical Casimir effect with trapped ions

To perform the mapping of the Hamiltonian $\hat{H}$ to an ion chain that has discrete positions, we first need to express it in discretized variables. This can be achieved by dividing the real axis ${\mathbb{R}}$ in suitable intervals $[{z}_{j},{z}_{j+1}]$, with ${z}_{j}\lt {z}_{j+1},j\in {\mathbb{Z}}$. For simplicity, we describe here the case of equidistant ion spacings, where the intervals are of equal length $d={z}_{j}-{z}_{j+1}$. It is straightforward to generalize the subsequent discussion to intervals of non-equal length. This permits one to take a non-equidistant distribution of the ions and a resulting variation of nearest-neighbor coupling strengths into account.

To arrive at a discretized Hamiltonian, we introduce the coarse-grained operators

$$\begin{eqnarray}&&{\hat{A}}_{j}=\displaystyle \frac{1}{\sqrt{d}}{\displaystyle \int }_{{z}_{j}}^{{z}_{j+1}}\hat{A}(z){\rm{d}}z,\end{eqnarray} \tag{ 23a }$$

$$\begin{eqnarray}&&{\hat{{\rm{\Pi }}}}_{j}=\displaystyle \frac{1}{\sqrt{d}}{\displaystyle \int }_{{z}_{j}}^{{z}_{j+1}}\hat{{\rm{\Pi }}}(z){\rm{d}}z,\end{eqnarray} \tag{ 23b }$$

which fulfill the commutation relations

$$\begin{eqnarray}&&[{\hat{A}}_{i},{\hat{A}}_{j}]=0,\end{eqnarray} \tag{ 24a }$$

$$\begin{eqnarray}&&[{\hat{{\rm{\Pi }}}}_{i},{\hat{{\rm{\Pi }}}}_{j}]=0,\end{eqnarray} \tag{ 24b }$$

$$\begin{eqnarray}&&[{\hat{A}}_{i},{\hat{{\rm{\Pi }}}}_{j}]={\rm{i}}{\hslash }{\delta }_{i,j}.\end{eqnarray} \tag{ 24c }$$

By using these operators, we obtain a discretized version of the Hamiltonian $\hat{H}$

$$\begin{eqnarray}&&{\hat{H}}^{{\rm{d}}}={\hat{H}}_{0}^{{\rm{d}}}+{\hat{H}}_{1}^{{\rm{d}}}(t),\end{eqnarray} \tag{ 25 }$$

with

$$\begin{eqnarray}&&{\hat{H}}_{0}^{{\rm{d}}}=\displaystyle \frac{1}{2}\displaystyle \sum _{i\in {\mathbb{Z}}}[{\hat{{\rm{\Pi }}}}_{i}^{2}+{d}^{-2}{({\hat{A}}_{i+1}-{\hat{A}}_{i})}^{2}]\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}&&{\hat{H}}_{1}^{{\rm{d}}}(t)=\displaystyle \frac{1}{2}\displaystyle \sum _{i\in {\mathbb{Z}}}{c}_{1}^{i}(t){{\hat{A}}_{i}}^{2},\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}&&{c}_{1}^{i}(t)=\displaystyle \frac{1}{d}{\displaystyle \int }_{{z}_{i}}^{{z}_{i+1}}{c}_{1}(t,z){\rm{d}}z.\end{eqnarray} \tag{ 28 }$$

In the above, we assumed that the modes of most importance are slowly varying on the length scale induced by the interval lengths d, which is well satisfied for the low-energetic modes. In the limit $d\to 0$, we recover the dynamics induced by the original Hamiltonian $\hat{H}$.

To implement the Hamiltonian ${\hat{H}}^{{\rm{d}}}$, we map it to the radial motion of a linear ion chain. Hereby, the position and momentum of each ion represent, respectively, the fields $\hat{A}(z)$ and $\hat{{\rm{\Pi }}}(z)$ averaged over one of the intervals $[{z}_{i},{z}_{i+1}]$. We consider a linear chain of N ions, confined in a suitable trapping potential ${V}_{\mathrm{trap}}$ and with equilibrium positions ${{\bf{R}}}_{1},{{\bf{R}}}_{2},...,{{\bf{R}}}_{N}$. If we assume that the deviations from the equilibrium positions are small, we can apply a second-order tailor expansion around the equilibrium positions. In this harmonic approximation, the motional degrees of freedom along different symmetry directions are uncoupled. In the following, we focus on the radial motion along the x-direction, which is described by the Hamiltonian

$$\begin{eqnarray}&&{\hat{H}}^{\mathrm{ions}}=\displaystyle \frac{1}{2m}\displaystyle \sum _{i=1}^{N}\;{\hat{P}}_{i}^{2}-\displaystyle \frac{{Z}^{2}{e}^{2}}{8\pi {\epsilon }_{0}}\displaystyle \sum _{i\gt j}\;\displaystyle \frac{{({\hat{X}}_{i}-{\hat{X}}_{j})}^{2}}{| | {{\bf{R}}}_{i}-{{\bf{R}}}_{j}| {| }^{3}}+\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}\;\displaystyle \frac{{\partial }^{2}{V}_{\mathrm{trap}}(t,{{\bf{R}}}_{i})}{\partial {x}^{2}}{\hat{X}}_{i}^{2}.\end{eqnarray} \tag{ 29 }$$

By applying the canonical transformation

$$\begin{eqnarray}&&{\hat{X}}_{i}\to {(-1)}^{i}{\hat{X}}_{i},\end{eqnarray} \tag{ 30a }$$

$$\begin{eqnarray}&&{\hat{P}}_{i}\to {(-1)}^{i}{\hat{P}}_{i},\end{eqnarray} \tag{ 30b }$$

we obtain

$$\begin{eqnarray}&&{\hat{H}}^{\mathrm{ions}}=\displaystyle \frac{1}{2m}\displaystyle \sum _{i=1}^{N}\;{\hat{P}}_{i}^{2}-\displaystyle \frac{1}{2}\displaystyle \sum _{i\gt j}\;{(-1)}^{i-j}{k}_{i,j}{({\hat{X}}_{i}-{\hat{X}}_{j})}^{2}+\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}\;{\chi }_{i}(t){\hat{X}}_{i}^{2},\end{eqnarray} \tag{ 31 }$$

with

$$\begin{eqnarray}&&{k}_{i,j}=\displaystyle \frac{{Z}^{2}{e}^{2}}{4\pi {\epsilon }_{0}}| | {{\bf{R}}}_{i}-{{\bf{R}}}_{j}| {| }^{-3},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\chi }_{i}{(t)=\displaystyle \frac{{\partial }^{2}{V}_{\mathrm{trap}}(t,{R}_{i})}{\partial {x}^{2}}-\displaystyle \sum _{j\ne i}(1-(-1)}^{i-j}){k}_{i,j}.\end{eqnarray} \tag{ 33 }$$

The mapping of the dynamics of the radiation field in a variable-length cavity to the dynamics of the ion chain is achieved by the formal similarity of equations (25) and (31). If we restrict the phonon Hamiltonian to nearest-neighbor interactions and consider an equidistant distribution of the ions, ${\hat{H}}^{\mathrm{ions}}$ indeed reproduces the Hamilton ${\hat{H}}^{{\rm{d}}}$, if we establish the following relations

$$\begin{eqnarray}&&{\hat{{\rm{\Pi }}}}_{i}={\hat{P}}_{i}/\sqrt{m},\end{eqnarray} \tag{ 34a }$$

$$\begin{eqnarray}&&{\hat{A}}_{i}={\hat{X}}_{i}\sqrt{m},\end{eqnarray} \tag{ 34b }$$

$$\begin{eqnarray}&&{d}^{-2}={k}_{i+1,i}/m,\end{eqnarray} \tag{ 34c }$$

$$\begin{eqnarray}&&{c}_{1}^{i}(t)={\chi }_{i}(t)/m.\end{eqnarray} \tag{ 34d }$$

The translation table from simulated objects to simulating objects is summarized in table 1.

Table 1.  Summary of the connection between the simulated objects (radiation field in variable-length cavity) and the simulating objects (ion chain with time-dependent trapping potential).

Simulated objects	Simulating objects
Photons	Phonons
Field operators	Position and momentum operators
${\hat{A}}_{i}$ and ${\hat{{\rm{\Pi }}}}_{i}$	of the radial motion ${\hat{X}}_{i}$ and ${\hat{P}}_{i}$
Variable-length cavity	Spatial- and time-dependent trap-
modeled by ${c}_{1}^{i}(t)$	ping potential ${V}_{\mathrm{trap}}$, respectively ${\chi }_{i}$
Discretized	Hamiltonian ${\hat{H}}^{\mathrm{ions}}$ describing the
Hamiltonian ${\hat{H}}^{{\rm{d}}}$	radial motion of the ion chain

The additional interaction terms beyond nearest neighbors could, to a large extent, be reabsorbed in a different choice for the approximations leading to ${\hat{H}}^{{\rm{d}}}$. In fact, the discretization $\tfrac{1}{2}{\sum }_{j\in {\mathbb{Z}}}{d}^{-2}{({\hat{A}}_{j+1}-{\hat{A}}_{j})}^{2}$ of the term $\tfrac{1}{2}{\int }_{{\mathbb{R}}}{\left(\tfrac{\partial \hat{A}}{\partial z}(z)\right)}^{2}{\rm{d}}z$ is not unique. The chosen discretization is motivated by the following approximation of the first derivative

$$\begin{eqnarray}&&\displaystyle \frac{\partial \hat{A}}{\partial z}(z)\approx \displaystyle \frac{1}{d}[\hat{A}(z+d)-\hat{A}(z)].\end{eqnarray} \tag{ 35 }$$

There are, however, other possible approximations, and as such other possible versions of ${\hat{H}}^{{\rm{d}}}$ (which give proper results if d is sufficiently small). In principle, this liberty could be exploited to account for interaction terms beyond nearest neighbors, but we will see below in section 6 that this is not necessary to get good qualitative and even quantitative agreement using realistic parameters. The underlying reason is the fast decrease of dipolar interactions with distance. Since the integral over dipolar interactions in one dimension converges, any perturbative effect due to interactions beyond nearest neighbors will saturate quickly when increasing the system size.

The remaining ingredient to simulate the DCE is the proper choice of the coefficients ${\chi }_{i}$. By virtue of equations (28) and (34d), these are directly connected to the function ${c}_{1}(t,z)$. Since the outermost ions represent the mirrors while the central part of the ion chain is to simulate the vacuum within the cavity, the coefficients ${\chi }_{i}$ need to vary across the chain. This can be achieved by exploiting that the ${\chi }_{i}$ are not only determined by the confining potential ${V}_{\mathrm{trap}}$ along the x-direction, but also by the Coulomb repulsion between the neighboring ions, see equation (33). As we will show in the next section, by properly balancing both contributions, we can design suitable coefficients ${\chi }_{i}$ , even with a small number of electrode segments.

To induce the DCE, the coefficients ${\chi }_{i}$, moreover, need to vary in time during stage II, which can be realized by a time-dependent trapping potential ${V}_{\mathrm{trap}}$. Since a small spatial motion of the mirrors corresponds to a change of the ${\chi }_{i}$ only over a small number of ions, this requires some local addressability. A suitable modulation can be achieved by combining a time-independent electric potential ${V}_{{\rm{E}}}$ (including the RF potential), generated via the segmented Paul trap, with an additional time-dependent optical potential ${V}_{{\rm{O}}}$ derived from a laser that addresses only one or a few of the ions [31]. Using the time-dependent optical potential, we can vary the boundary of the cavity during stage II of the experiment. This completes all required ingredients for the simulation of the DCE. In the next section, we demonstrate that good agreement to the ideal model can be obtained already for about 20 ions and in present-day architectures.

The trap we consider consists of only six DC electrodes. This turns out to be sufficient to form a suitable electric potential ${V}_{{\rm{E}}}$, which we compute by using the framework presented in [37] and by applying the gapless plane approximation. For the sake of simplicity, we assume that the extension of the RF-electrodes along the z-axis as well as the length of the DC-electrodes along the x-axis ${l}_{{\rm{e}}}$ are infinite. Inspired by [32], we assume a possible setup with $w=80\;\mu {\rm{m}},{h}_{0}=80\;\mu {\rm{m}},{d}_{{\rm{g}}}=230\;\mu {\rm{m}}$, and we consider singly-charged ions. We set the voltages of the DC electrodes to the values

$$\begin{eqnarray}&&{\phi }_{1}={\phi }_{2}={\phi }_{5}={\phi }_{6}=-5.61\;{\rm{V}},\\ &&{\phi }_{3}={\phi }_{4}=1.75\;{\rm{V}}\end{eqnarray} \tag{ 36 }$$

and use the RF electrodes to induce a confining potential that corresponds to a trapping frequency

$$\begin{eqnarray}&&{\omega }_{\mathrm{RF}}=\sqrt{7.40\bar{k}/m}.\end{eqnarray} \tag{ 37 }$$

Here, $\bar{k}=\tfrac{{e}^{2}}{4\pi {\epsilon }_{0}{\overline{{\rm{\Delta }}R}}^{3}}$ is the average nearest-neighbor coupling strength, with

$$\begin{eqnarray}&&\overline{{\rm{\Delta }}R}=| | {{\bf{R}}}_{20}-{{\bf{R}}}_{1}| | /19=4.00\;\mu {\rm{m}}\end{eqnarray} \tag{ 38 }$$

being the average nearest-neighbor distance. For calcium ions (${}^{40}{\mathrm{Ca}}^{+}$), for example, we obtain

$$\begin{eqnarray}&&\sqrt{\bar{k}/m}=2\pi \times 1.17\;\mathrm{MHz}\end{eqnarray} \tag{ 39 }$$

and ${\omega }_{\mathrm{RF}}=2\pi \times 3.18\;\mathrm{MHz}$. All these values lie in the range of existing experimental setups.

The calculation of the coefficients ${\chi }_{i}$ for these parameters yields the result depicted in figure 2. Here, we took the non-equidistant distribution of the ions for this trapping potential as well as all possible interactions (beyond nearest neighbors) into account.

Zoom In Zoom Out Reset image size

The trapping potential is adjusted such that the ${\chi }_{i}$ deviate significantly from zero only in the outer regions of the chain. Since the electric field experiences an exponential damping within the mirrors, a rather small number of ions proves sufficient to model the space inside the mirrors, in this example ions 1–4 and 17–20. The field inside the cavity is represented by the inner part of the ion chain, i.e., the ions 5–16.

We subject this ion chain to the three-stage protocol defined in section 2. The time dependence in stage II that we consider corresponds to a periodically oscillating left mirror, i.e.

$$\begin{eqnarray}&&r(t)={r}_{0},\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}l(t)={l}_{0}+\delta \left\{\begin{array}{ll}0 & t\in [{t}_{0},{t}_{1}),\\ {\mathrm{sin}}^{2}({\omega }_{{\rm{D}}}(t-{t}_{1})/2) & t\in [{t}_{1},{t}_{2}),\\ {\mathrm{sin}}^{2}({\omega }_{{\rm{D}}}({t}_{2}-{t}_{1})/2) & t\in [{t}_{2},\infty ),\end{array}\right.\end{eqnarray} \tag{ 41 }$$

where r₀ and l₀ denote the initial positions of the right and left mirrors, respectively. Further, $\delta /2$ is the amplitude of the variation and ${\omega }_{{\rm{D}}}$ the driving frequency. This choice of the time dependence of the boundaries leads to an efficient photon production [38, 39].

In order to simulate these mirror trajectories, we use a laser beam to change the radial confinement of ion i = 5 such that

$$\begin{eqnarray}{\omega }_{O}^{2}(t)=\alpha \displaystyle \frac{\bar{k}}{m}\left\{\begin{array}{ll}0 & t\in [{t}_{0},{t}_{1}),\\ {\mathrm{sin}}^{2}({\omega }_{{\rm{D}}}(t-{t}_{1})/2) & t\in [{t}_{1},{t}_{2}),\\ {\mathrm{sin}}^{2}({\omega }_{{\rm{D}}}({t}_{2}-{t}_{1})/2) & t\in [{t}_{2},\infty ),\end{array}\right.\end{eqnarray} \tag{ 42 }$$

where we choose $\alpha =0.6$. It is hereby not a strict requirement to address precisely a single ion. Addressing several neighboring ions due to a larger beam waist just corresponds to a larger variation δ of the position of the left mirror. Similarly, a different depth of the optical potential $\alpha \bar{k}/m$ also just amounts to a different δ.

In the following section, we will numerically compare the dynamics for the ideally conducting mirrors modeled according to Moore [3] with our ion-chain quantum simulation. For evaluating the photon number, we choose the canonical set of mode functions for the experimental stage III, defined by equations (9a)–(9c). We order the mode frequencies ${\omega }_{{\ell }}$ as ${\omega }_{1}\lt {\omega }_{2}\lt \cdots$ and denote with ${\hat{n}}_{1},{\hat{n}}_{2}...$ the corresponding photon number operators. For better comparison, we choose ${r}_{0},{l}_{0}$, and δ such that the frequency of the lowest instantaneous eigenmode in the cavity matches the frequency of the lowest instantaneous vibrational mode of the ion chain at time instances corresponding to the maximal and minimal cavity length, i.e., ${\omega }_{{\rm{D}}}(t-{t}_{1})=0$ and ${\omega }_{{\rm{D}}}(t-{t}_{1})=\pi$. This matching is obtained by

$$\begin{eqnarray}&&{r}_{0}-{l}_{0}=15.22d,\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&\delta =0.72d,\end{eqnarray} \tag{ 44 }$$

where $d={(\bar{k}/m)}^{-1/2}$ denotes the length of the discretization intervals $[{z}_{i},{z}_{i+1}]$ used in the previous section. The result for the 'cavity length', ${r}_{0}-{l}_{0}$, can be understood by recalling that each ion represents the averaged field in one of the intervals $[{z}_{i},{z}_{i+1}]$ and that the field inside the cavity is roughly represented by the 12 inner ions. The deviation between $15.22d$ and the length $12d$ expected from this simple consideration is caused by the non equidistant distribution of the ions and the discretization of the field.

We have now all the necessary parameters to numerically simulate the DCE as can be studied in a realistic ion-trap experiment. We assume that during the experimental stage I the ion chain resides in its vibrational ground state. Since this initial state corresponds to the vacuum of the radiation field, all photons measured in stage III are those that have been produced in stage II. In figure 3, we show the final photon number in the modes 1 and 2 for 20 periods of mirror oscillations during the stage II. Since Hamiltonian (31) is a quadratic bosonic theory, exact predictions for the ion quantum simulator can be calculated by solving the linear Heisenberg equations of motion for the annihilation and creation operators, while the results for Moore's model were evaluated numerically by using the method of images discussed in [3].

Zoom In Zoom Out Reset image size

As a function of the driving frequency ${\omega }_{{\rm{D}}}$, one finds peaks of high photon production centered around integer multiples of the frequency ${\omega }_{1}$. This finding is in accordance with well known analytical results [39]. The main contribution to these peaks stems from single-mode and two-mode squeezing, which is connected to the resonance condition ${\omega }_{{\rm{D}}}=\langle {\omega }_{{{\ell }}_{1}}{\rangle }_{{\rm{T}}}+\langle {\omega }_{{{\ell }}_{2}}{\rangle }_{{\rm{T}}}$, with $\langle {\omega }_{{\ell }}{\rangle }_{{\rm{T}}}$ denoting the time average of the instantaneous eigenfrequency ${\omega }_{{\ell }}(t)$ (averaged over one oscillation period of the mirror). The peaks depicted in figure 3 are slightly shifted from integer multiples of ${\omega }_{1}$, because ${\omega }_{{\ell }}\ne \langle {\omega }_{{\ell }}{\rangle }_{{\rm{T}}}$ and due to artefacts caused by the discretization of the model.

Additionally, the trapped-ion quantum simulator allows one to monitor the photon production over time. Figure 4 displays the corresponding results for the average photon number in mode 1, for a system driven with the frequency ${\omega }_{{\rm{D}}}=2\langle {\omega }_{1}{\rangle }_{{\rm{T}}}$ (with $\langle {\omega }_{1}{\rangle }_{{\rm{T}}}\approx 1.028{\omega }_{1}$). For this frequency, which lies is at the first peak in figure 3(a), the photon production is dominated by single-mode squeezing. When the resonance condition for single-mode squeezing, ${\omega }_{{\rm{D}}}=2\langle {\omega }_{{\ell }}{\rangle }_{{\rm{T}}}$, is met, the average photon number is approximately given by [5]

$$\begin{eqnarray}&&\langle {n}_{{\ell }}\rangle ={\mathrm{sinh}}^{2}\left[\langle {\omega }_{{\ell }}{\rangle }_{{\rm{T}}}\displaystyle \frac{\delta }{4({r}_{0}-{l}_{0})}({t}_{2}-{t}_{1})\right],\end{eqnarray} \tag{ 45 }$$

valid if $\delta /({r}_{0}-{l}_{0})$ is sufficiently small and the duration of the periodic driving is sufficiently short. At short times, this approximate expression indeed coincides with the results for the ion trap as well as for the ideal Moore's model. At larger times, the curves start to deviate but the qualitative agreement remains satisfactory. For reasons of comparison, we also evaluate the phonon production for an improved choice of the time dependence of the optical trapping potential ${\omega }_{O}^{2}(t)$, which minimizes discretization effects connected to the fact that just a single ion was used to represent the motion of the left mirror. Hereby, we chose ${\omega }_{O}^{2}(t)$ such, that the instantaneous eigenfrequency of the lowest vibrational mode matches the instantaneous eigenfrequency of the lowest optical mode in the idealized cavity setup for all times. The optimized trajectory lies much closer to the prediction of Moore's model than the simple sine wave, which shows that the deviations are mainly artefacts from using a discretized representation for moving the mirror.

Zoom In Zoom Out Reset image size

The additional small peaks of the blue curves in figures 3(a) and (b) are also artefacts connected to use of a single ion to represent the motion of the left mirror. These artefacts can be reduced by smoothing the mirror motion, i.e., by increasing the number of ions that experience a periodic modulation during stage II. Furthermore, by increasing the number of ions representing the radiation field inside the cavity, the distribution of the low-lying eigenfrequencies further approaches the equidistant distribution of the eigenfrequencies in an ideal cavity. In this way, it is possible to reduce the slight shifts of the peaks seen in figure 3 that are caused by discretization effects.

As the above results show, an ion chain with realistic parameters can indeed simulate the photon production in the DCE, where we find a good agreement to Moore's model [3] already for 20 ions. Small deviations do appear due to the limited number of ions. This is no fundamental limitation, however, and by increasing the number of ions, or by using additional electrodes, it will be possible to reduce these artefacts and further improve the simulation of the DCE. As mentioned previously, the protocol is also rather resilient towards a change in the (time-dependent) trapping potential. Insufficient control of the spatial dependence will simply model a slightly modified cavity.

Additionally, in a realistic experiment, one has to make sure that the observed phonons are not generated by a heating of the ion chain. In the experimental setup of [32], which has similar parameters to the ones discussed above, a spectral density of electric field noise of ${S}_{{\rm{E}}}(\omega )\leqslant 3.8\cdot {10}^{-13}\tfrac{{V}^{2}}{{m}^{2}}\;\mathrm{Hz}$ has been reported for $\omega =2\pi \times 1.38\;\mathrm{MHz}$. Under the assumption that heating is dominated by electric field noise and that ${S}_{{\rm{E}}}(\omega )\omega$ is approximately frequency independent, we obtain a heating rate of the lowest-lying mode (with ${\omega }_{1}\approx 0.21{(\bar{k}/m)}^{1/2}=2\pi \times 0.241\;\mathrm{MHz}$) of $1.31\mathrm{quanta}\;{\mathrm{ms}}^{-1}$. Hereby, we could neglect the cross coupling between the RF and noise fields because of the relatively small frequency ${\omega }_{1}$. For the higher modes ${\omega }_{2}$, ${\omega }_{3},...$ , the effect of the electric field noise is even smaller, since the heating rate it causes (ignoring cross coupling between RF and noise fields) scales as $1/{\omega }^{2}$. Another source of heating are scattered photons from the laser beam. The corresponding heating rates can be suppressed by using sufficiently intense and sufficiently detuned standing-wave laser fields. For a sufficiently far detuned trapping laser, the main source of decoherence is parametric heating caused by intensity fluctuations of the laser beam. Starting from the ground state of the vibrational ground state, parametric heating causes a transition to a two photon state with one photon being in mode ℓ₁ and the other being in mode ℓ₂ with the rate (see [40, 41] for a detailed theoretical description)

$$\begin{eqnarray}{R}_{\mathrm{2\leftarrow 0}}=\displaystyle \frac{\pi {\omega }_{{\rm{O}}}^{4}}{8{\omega }_{{{\ell }}_{1}}{\omega }_{{{\ell }}_{2}}}{|{\beta }_{{{\ell }}_{1}}|}^{2}{|{\beta }_{{{\ell }}_{2}}|}^{2}{S}_{\epsilon }({\omega }_{{{\ell }}_{1}}+{\omega }_{{{\ell }}_{2}})\left\{\begin{array}{cc}1 & {{\ell }}_{1}={{\ell }}_{2}\\ 2 & {{\ell }}_{1}\ne {{\ell }}_{2}\end{array}\right.\end{eqnarray} \tag{ 46 }$$

with

$$\begin{eqnarray}&&{S}_{\epsilon }(\omega )=\displaystyle \frac{2}{\pi {I}_{0}^{2}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}\tau \mathrm{cos}(\omega \tau )\langle (I(t)-{I}_{0})(I(t+\tau )-{I}_{0}){\rangle }_{{\rm{T}}},\end{eqnarray} \tag{ 47 }$$

being the one-sided noise power spectrum of the fractional intensity noise, I(t) being the time dependent intensity of the trapping laser, and ${I}_{0}=\langle I(t){\rangle }_{{\rm{T}}}$ being the time averaged laser intensity. The parameter ${\beta }_{{\ell }}$ (with $| {\beta }_{{\ell }}{| }^{2}\leqslant 1$) takes the overlap of the normalized vibrational mode ℓ with the sites of the ion chain affected by the laser drive into account. In our particular example, we obtain for mode 1 $| {\beta }_{1}| \approx 0.11$ (roughly scales with $| {\beta }_{{\ell }}{| }^{2}\propto 1/N$). Experimental parameters for a possible realistic implementation are discussed in [31] in case of ${}^{24}{\mathrm{Mg}}^{+}$ ions. It is estimated that with a retro-reflected beam from a far off detuned commercial laser source an optical trapping potential with an oscillation frequency of ${\omega }_{O}=2\pi \times 2.4\;\mathrm{MHz}$ could be realized. This oscillation frequency is well beyond our requirements (see equation (42)) even if one takes the higher mass of ${}^{40}{\mathrm{Ca}}^{+}$ into account. It is also estimated, that this configuration would lead to a heating rate of $0.075\;\mathrm{quanta}\;{\mathrm{ms}}^{-1}$ (caused by intensity and frequency noise of the laser), which is well below our estimated heating rate caused by the electric field noise of the ion trap. It should also be possible to achieve similar experimental parameters for other species of ions.

These heating rates have to be compared to the relevant experimental time scales. The data at the first peak in figure 3(a) corresponds to a duration of the experimental stage II of about $({t}_{2}-{t}_{1})=0.041\;\mathrm{ms}$. Thus, heating is expected to increase the average phonon number of the 1st mode by roughly 0.054 phonons, which is one order of magnitude smaller than the number of phonons generated by the simulated DCE. Even more, one could further reduce the effect of ion heating due to electric field noise by decreasing the average nearest-neighbor distance between the ions. This would increase the frequencies of the modes and hence decrease the time needed to run the experiment. Therefore, according to these numbers it will be possible to cleanly observe the first as well as higher peaks in experiment.

As discussed above, the main idea of our ion-chain quantum simulation of the DCE is to map the photons of the radiation field on the phonons of the radial ion motion. Thus, for probing the radiation generated by the DCE, we have to measure the generated phononic excitations in stage III of the experiment. This can be done with high temporal resolution and high accuracy on the single phonon level by using the methods available for ion chains [28], which is one of the main advantages of our ion-chain quantum simulation compared to other schemes [8–10].

One possibility for evaluating the number of phonons populating mode ℓ is to drive the corresponding red or blue detuned sideband for a short time period ${\rm{\Delta }}t$ by addressing a single ion with a laser beam. Hereby, one should choose an ion that takes part in the collective motion described by the mode ℓ. By repeating this experiment for several runs, the probability for exciting the ion, ${P}_{{\rm{e}}}({\rm{\Delta }}t)$, can be determined. For driving the blue detuned sideband one obtains

$$\begin{eqnarray}&&{P}_{{\rm{e}}}({\rm{\Delta }}t)=\displaystyle \sum _{n}\;{p}_{{\ell }}(n){\mathrm{sin}}^{2}(\sqrt{n+1}\;{{\rm{\Omega }}}_{\mathrm{0,1}}{\rm{\Delta }}t),\end{eqnarray} \tag{ 48 }$$

where ${p}_{{\ell }}(n)$ is the probability of finding n phonons in mode ℓ and ${{\rm{\Omega }}}_{\mathrm{0,1}}$ is a characteristic Rabi frequency determined by the intensity of the applied laser beam and the corresponding Lamb–Dicke parameter. For short time periods ${\rm{\Delta }}t$, the above expression simplifies to ${P}_{{\rm{e}}}({\rm{\Delta }}t)\approx {{\rm{\Omega }}}_{0,1}^{2}{\rm{\Delta }}{t}^{2}\langle 1+{\hat{n}}_{{\ell }}\rangle +O({\rm{\Delta }}{t}^{4})$, which allows us to evaluate the average phonon number [42, 43]. By measuring ${P}_{{\rm{e}}}({\rm{\Delta }}t)$ for several ${\rm{\Delta }}t$ over a longer time period, it is possible to determine the phonon number distribution ${p}_{{\ell }}(n)$ by calculating the Fourier transform of ${P}_{{\rm{e}}}({\rm{\Delta }}t)$ [44, 45]. In this way, we can probe not only the average photon number, but even the detailed photon statistics of the radiation generated by the DCE.

In principle, it is also possible to apply other methods developed for probing the quantum state of motion and accessing other observables [28]. The choice of the most suitable method depends on the experimental parameters of the specific setup.