Particle pair creation by inflation of quantum vacuum fluctuations in an ion trap

TL;DR

This paper demonstrates an experimental analog of cosmological particle creation using trapped ions, creating phonon pairs from quantum vacuum fluctuations to explore early universe phenomena and quantum information applications.

Contribution

It introduces a method to simulate particle pair creation during inflation with trapped ions, linking cosmology and quantum information science.

Findings

Successful creation of phonon pairs from vacuum fluctuations

Preparation of entangled motional states in ions

Potential for optimizing quantum simulation of cosmological processes

Abstract

The creation of matter and structure in our universe is currently described by an intricate interplay of quantum field theory and general relativity. Signatures of this process during an early inflationary period can be observed, while direct tests remain out of reach. Here, we study an experimental analog of the process based on trapped atomic ions. We create pairs of phonons by tearing apart quantum vacuum fluctuations. Thereby, we prepare ions in an entangled state of motion. Controlling timescales and the coupling to environments should permit optimizing efficiencies while keeping the effect robust via established tools in quantum information processing (QIP). This might also permit to cross-fertilize between concepts in cosmology and applications of QIP, such as, quantum metrology, experimental quantum simulations and quantum computing.

Equations

$${\ddot{\varphi }}_k+\left[c^2{k}^2+a^2(t)\frac{m^2{c}^4}{{\hslash }^2}\right]{\varphi }_k={\ddot{\varphi }}_k+{\Omega }_k^2(t){\varphi }_k=0,$$

$$\displaystyle ds^{2}=c^{2}d\tau^{2}-a^{2}(\tau)d\mbox{\boldmath$r$}^{2}$$

$$\displaystyle\ddot{\phi}_{\mbox{\boldmath$\scriptstyle k$}}+\left[\frac{m^{2}c^{4}}{\hbar^{2}}+\frac{c^{2}\mbox{\boldmath$k$}^{2}}{a^{2}}-\frac{3}{2}\left(\frac{\ddot{a}}{a}+\frac{1}{2}\frac{\dot{a}^{2}}{a^{2}}\right)\right]\phi_{\mbox{\boldmath$\scriptstyle k$}}=0\,.$$

$$\displaystyle ds^{2}=a^{2}(t)\left[c^{2}dt^{2}-d\mbox{\boldmath$r$}^{2}\right]\,.$$

$$\displaystyle\ddot{\phi}_{\mbox{\boldmath$\scriptstyle k$}}+\left[a^{2}\,\frac{m^{2}c^{4}}{\hbar^{2}}+c^{2}\mbox{\boldmath$k$}^{2}-\frac{\ddot{a}}{a}\right]\phi_{\mbox{\boldmath$\scriptstyle k$}}=0\,.$$

$$\displaystyle\ddot{A}_{\mbox{\boldmath$\scriptstyle k$},\sigma}+\left[a^{2}\,\frac{m^{2}c^{4}}{\hbar^{2}}+c^{2}\mbox{\boldmath$k$}^{2}\right]A_{\mbox{\boldmath$\scriptstyle k$},\sigma}=0\,.$$

$$\displaystyle\left[\partial_{t}^{2}-c^{2}\left(\mbox{\boldmath$\nabla$}-i\frac{q}{\hbar}\mbox{\boldmath$A$}\right)^{2}+\frac{m^{2}c^{4}}{\hbar^{2}}\right]\Phi=0\,,$$

$$\displaystyle\ddot{\phi}_{\mbox{\boldmath$\scriptstyle k$}}+\left[c^{2}\left(\mbox{\boldmath$k$}-\frac{q}{\hbar}\mbox{\boldmath$A$}\right)^{2}+\frac{m^{2}c^{4}}{\hbar^{2}}\right]\phi_{\mbox{\boldmath$\scriptstyle k$}}=0\,.$$

$$\gamma =\frac{\omega mc}{qE}\sim \frac{mc}{qA},$$

$$d{s}^2=2c{e}^{\kappa t}dtdU-e^{2\kappa t}d{U}^2,$$

$$d{s}^2=\frac{2cdTdU-d{U}^2}{(\kappa T{)}^2},$$

$$\left[{\partial }_T^2+2c{\partial }_T{\partial }_U+a^2(T)\frac{m^2{c}^4}{{\hslash }^2}\right]\Phi =0.$$

$$\left[{\partial }_T^2+c^2{K}^2+a^2(T)\frac{m^2{c}^4}{{\hslash }^2}\right]{\phi }_K(T)=0,$$

$$\delta {\ddot{q}}_i+{\omega }_{\mathrm{rad}}^2(t)\delta q_i+\sum _j{M}_{ij}\delta q_j=0,$$

$$\delta {\ddot{q}}_I+{\omega }_{\mathrm{rad}}^2(t)\delta q_I+{\lambda }_I\delta q_I=0,$$

$${\rho }_{\text{fin}}=\hat{S}(r)\hat{D}(\alpha ){\rho }_{\text{ini}}{\hat{D}}^{†}(\alpha ){\hat{S}}^{†}(r).$$

$$\hat{H}(t)$$

$$=\hslash \omega (t)\left({\hat{a}}^{†}\hat{a}+\frac{1}{2}\right)-\frac{i\hslash }{4}\frac{d\ln \left(\omega (t)\right)}{dt}\left[{\hat{a}}^2-{\left({\hat{a}}^{†}\right)}^2\right].$$

$${\hat{H}}_F(t)=x_0{F}_0(t)(\hat{a}+{\hat{a}}^{†}),$$

$$|\text{ini}⟩$$

$$=\sqrt{1-e^{-2\beta }}\sum _n{e}^{-\beta n}{|n⟩}_{\text{A}}{|n⟩}_{\text{B}},$$

$$|\text{fin}⟩$$

$$={|0⟩}_{\text{A}}{|0⟩}_{\text{B}}-\frac{\xi }{2}{|1⟩}_{\text{A}}{|1⟩}_{\text{B}}-\frac{\xi }{\sqrt{8}}\left({|0⟩}_{\text{A}}{|2⟩}_{\text{B}}+{|2⟩}_{\text{A}}{|0⟩}_{\text{B}}\right)+\mathcal{O}(r^2).$$

$$E_F(\chi )=\frac{(1/2+\chi {)}^2}{2\chi }\ln \left(\frac{(1/2+\chi {)}^2}{2\chi }\right)-\frac{(1/2-\chi {)}^2}{2\chi }\ln \left(\frac{(1/2-\chi {)}^2}{2\chi }\right),$$

$${\lambda }_1$$

$${\lambda }_2$$

Full text

Particle pair creation by inflation of quantum vacuum fluctuations in an ion trap

M. Wittemer1∗, F. Hakelberg1, P. Kiefer1, J.-P. Schröder1,

C. Fey2, R. Schützhold3,4,5, U. Warring1, T. Schaetz1

1Physikalisches Institut, Universität Freiburg, Hermann-Herder-Straße 3, 79104 Freiburg, Germany

2Zentrum für Optische Quantentechnologien, Universität Hamburg,

Fachbereich Physik, Luruper Chaussee 149, 22761 Hamburg, Germany

3Fakultät für Physik, Universität Duisburg-Essen, Lotharstrasse 1, 47057 Duisburg, Germany

4Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstrasse 400, 01328 Dresden, Germany

5Institut für Theoretische Physik, Technische Universität Dresden, 01062 Dresden, Germany

∗To whom correspondence should be addressed; E-mail: matthias.wittemer@physik.uni-freiburg.de

The creation of matter and structure in our universe is currently described by an intricate interplay of quantum field theory and general relativity. Signatures of this process during an early inflationary period can be observed, while direct tests remain out of reach. Here, we study an experimental analog of the process based on trapped atomic ions. We create pairs of phonons by tearing apart quantum vacuum fluctuations. Thereby, we prepare ions in an entangled state of motion. Controlling timescales and the coupling to environments should permit optimizing efficiencies while keeping the effect robust via established tools in quantum information processing (QIP). This might also permit to cross-fertilize between concepts in cosmology and applications of QIP, such as, quantum metrology, experimental quantum simulations and quantum computing.

The quantum vacuum is not empty, but filled with ubiquitous fluctuations, as implied by the Heisenberg uncertainty principle. These fluctuations remain invisible but have observable consequences such as the Lamb shift of atomic spectral lines  (?), van-der-Waals and Casimir forces  (?), or the “stimulation” of spontaneous emission  (?, ?). They can be pictured as the random creation and annihilation of virtual pairs of particles and anti-particles. Intriguingly, extreme conditions such as a rapid cosmic expansion can tear these virtual particles apart and, thereby, turn them into real particles  (?, ?). These particle pairs represent lasting excitations of the quantum field and their creation generates quantum entanglement at large distances. Similar mechanisms cause the Sauter-Schwinger effect  (?) and Hawking radiation  (?). In addition, pair-creation can be realized by parametric driving, e.g., via the dynamical Casimir effect (?). According to our standard model of cosmology, the tearing apart of fluctuations by an expanding space-time during the inflationary period of the early universe explains the creation of the seeds for structure formation  (?). Realizing such extreme conditions in the laboratory remains elusive, but some analog features have been observed in several experimental platforms  (?, ?, ?, ?, ?, ?, ?, ?, ?). However, preserving the fragile quantum dynamics described above requires fast control of the system on the level of single quanta, as well as close to ideal isolation from the environment. Trapped atomic ions are well suited to study fundamental quantum dynamics as they feature unique fidelities in preparation, control, and detection of quantum states  (?, ?).

Here, we present the creation of pairs of particles, more precisely phonons, by inflating quantum vacuum fluctuations of the motion of ions in their storage field. This is accompanied by the creation of spatial entanglement. We explain the process as an experimental analog to an inflationary period of an early universe  (?, ?, ?).

For simplicity of the analog, let us consider a scalar (e.g., inflaton or Higgs) field $\Phi$ in 1+1 dimensions (other fields behave analogously, see  (?)). In terms of the conformal time $t$, the expanding space-time during cosmic inflation can be described by the Friedmann-Lemaître-Robertson-Walker metric $ds^{2}=a^{2}(t)\left[c^{2}dt^{2}-dx^{2}\right]$, where the time-dependent scale parameter $a(t)$ governs the cosmic expansion. After a spatial Fourier transform $\Phi(t,x)\to\phi_{k}(t)$, the equation of motion for a specific mode $k$ reads

[TABLE]

where $c$ denotes the speed of light, $\hbar$ is the reduced Planck constant, while $k$ labels the related harmonic potential and identifies the momentum $\hbar k$ of one excitation with mass $m$. The squared frequency $\Omega_{k}^{2}(t)$ contains the internal (propagating) contribution $c^{2}k^{2}$ as well as the external dynamics $\propto a^{2}(t)$ in the second term. Initially ($\Omega_{k}\approx ck$), the internal dynamics dominates and the mode $\phi_{k}$ oscillates freely, staying close to its ground state. During the course of inflation, $a(t)$ grows and the cosmic expansion stretches the physical wavelength $\lambda(t)=2\pi a(t)/k$ of the mode. If the external variation $a^{2}(t)$ is sufficiently fast to dominate the dynamics, i.e., the temporal change $\dot{\Omega}_{k}$ is not small compared to $\Omega_{k}^{2}$, the quantum state departs from the ground state, cf. Figs. 1A and B. The mode transforms to an excited, squeezed state, which corresponds to the creation of pairs of particles. Thereby, the cosmic expansion tears apart the initial quantum vacuum fluctuations and turns them into pairs of particles with opposite momenta $\pm\hbar k$, generating quantum entanglement at large distances  (?), cf. Fig. 1C. Relevant parameters for the number of particle pairs created are the scale parameter $a(t)$ and the relative rate of frequency change $\dot{\Omega}_{k}/\Omega_{k}^{2}$. The above mechanism is universal and applies to any quantum harmonic oscillator subjected to large and rapid change of its frequency $\omega_{a}(t)$, which causes squeezing of its wavefunction $\psi_{a}$. In an ion trap, the motional (phonon) modes represent the $\phi_{k}$, while changes of the trapping potential mimic the cosmic expansion $a(t)$ and, thus, stimulate the creation of pairs of phonons  (?).

In our experiments we study the motional modes of two atomic magnesium ions confined in a linear Paul trap, details on the experimental setup are given in  (?, ?). The two ions align along the axial direction of the trap and in the following, we focus on their harmonic motion along the weaker of the radial directions. We distinguish an in-phase (center-of-mass) mode from an out-of-phase (rocking) mode with frequency $\omega_{a}$. The latter is intrinsically robust against homogeneous noise fields, and thus, well suited to implement the dynamics according to Eq. (1). Our apparatus allows for real-time control of a trapping voltage $U$ that tunes $\omega_{a}$, as depicted in Fig. 2A. We define an analog scale parameter, postponing the internal dynamics of Eq. (1), according to $a_{a}(t)=\omega_{a}(t)/\omega_{a}(0)$, see  (?). In our experiments we ramp $\omega_{a}/(2\pi)$, down and up, spanning $2.50(1)\text{\,}\mathrm{MHz}$ to $0.50(1)\text{\,}\mathrm{MHz}$, where each ramp corresponds to an inflation of space by about 1.6 $e$-foldings, see Fig. 2B. We ramp $U$ within $t_{\text{ramp}}=$1\text{,}\mathrm{\SIUnitSymbolMicro s}$$, yielding a non-adiabatic evolution of $\psi_{a}$, characterized by $\dot{\omega}_{a}/\omega_{a}^{2}\approx 5$, cf. Fig. 2C. Numerical simulations  (?) illustrate the generation of squeezing with an amplitude $r$ during a single ramp of $\omega_{a}$, see Fig. 2D, which corresponds to an instantaneous average particle number $\bar{n}_{\text{sq}}=\sinh^{2}\left(r\right)$. Due to the non-instantaneous switching of $\omega_{a}$, the Wentzel-Kramers-Brillouin (WKB) phase $\varphi_{a}(t)=\int_{0}^{t}\omega_{a}(t^{\prime})dt^{\prime}$ of $\psi_{a}$ evolves significantly during $t_{\text{ramp}}$ and leads to an oscillatory accumulation of $r$ and $\bar{n}_{\text{sq}}$.

In a first experimental realization (protocol depicted in Fig. 3A) we cool all modes close to their motional ground state, for the rocking mode with a measured, residual thermal $\bar{n}_{\text{th}}=0.03(6)$, see Fig. 3B. In order to increase the number of created particle pairs, we apply a sequence of two ramps to form one pulse: After ramping down $\omega_{a}$ and an appropriate duration $t_{\text{hold}}$, $\omega_{a}$ is ramped back up. Finally, we readout the motional state by mapping it onto the electronic state of one ion  (?, ?), and reconstruct the individual phonon number distribution $P_{n}$ (data points in Figs. 3B and C). Typically, experimental sequences are repeated $40000$ times for a given set of parameters. In comparison to the initial state, the data show increased populations for higher phonon number states $n$, indicating mode excitation due to the pulse. To evaluate our results, we decompose the excitation into contributions by $r$, a coherent displacement $\alpha$, and a thermal share $\bar{n}_{\text{th}}$. We fit the data with a parametrized phonon number distribution $P_{n}^{\text{par}}(r,\alpha,\bar{n}_{\text{th}})$ of a Gaussian state (bar chart in Fig. 3C), with $r=0.54(8)$ and $|\alpha|=0.88(6)$. Since all heating effects remain negligible within experimental sequences we keep $\bar{n}_{\text{th}}$ at the initial value 0.03. We attribute the coherent displacement to residual static differential stray fields in our setup  (?). Both excitation amplitudes, $r$ and $|\alpha|$, depend on the phase $\varphi_{a}$ accumulated during the pulse, i.e., they are determined by $t_{\text{ramp}}$, $t_{\text{hold}}$, and $\omega_{a}(t)$, see  (?).

In order to reverse the coherent displacement, we perform purifying echo sequences consisting of two pulses, separated by a duration of free evolution $t_{\text{free}}$, see Fig. 4A. To illustrate the echo effect of the second pulse we depict the corresponding phase space animation in Fig. 4B. The initial thermal state (i) is squeezed and coherently displaced by the first pulse (ii), similar to Fig. 3C. During $t_{\text{free}}$ the state oscillates, i.e., rotates on a circular path in phase space (iii). $\psi_{a}$ accumulates two distinct phases corresponding to excitations $r$ and $\alpha$ with frequencies $2\omega_{a}(t_{\text{free}})$ and $\omega_{a}(t_{\text{free}})$, respectively. Ideally, the second pulse is applied after a duration $\pi/\omega_{a}(t_{\text{free}})$ (or odd multiples), where the coherent displacement has evolved from $\alpha$ to $-\alpha$ (iv), and the final state after the second pulse (v) is characterized by $\alpha\approx 0$ and enhanced squeezing  (?). We experimentally realize the two-pulse sequence for variable $t_{\text{free}}$, perform the motional state analysis, and depict resulting squeezing and coherent excitations in Fig. 4C. As indicated by the data, for $t_{\text{free}}\approx$30.2\text{,}\mathrm{\SIUnitSymbolMicro s}, the squeezing is significantly increased, while the coherent excitation is reduced with respect to Fig. 3C. Accordingly, in Fig. 4D we depict the $P_{n}$ for $t_{\text{free}}=$30.2\text{\,}\mathrm{\SIUnitSymbolMicro s}. Here the squeezing is directly evidenced by increased populations $P_{n}$ for even states only, while odd states remain nearly unpopulated. This corresponds to the creation of pairs of phonons. For example, we detect single phonon pairs ($n=2$) in $P_{2}\approx$20\text{,}\mathrm{\char 37\relax} of our realizations ($\approx 8000$ in total), and the Wigner function indicates a suppression of the variance $\left(\Delta x\right)^{2}$ by $\approx$7.2\text{\,}\mathrm{d}\mathrm{B}.

Moreover, features of entanglement of the particle pairs can be witnessed in the spatial degrees of freedom of our two-ion crystal  (?). To first order, the wavefunction’s squeezing excitation can be visualized in a simplifying way with a term $\propto r\left(\ket{{}^{\bullet}\,{}_{\bullet}}+\ket{{}_{\bullet}\,{}^{\bullet}}\right)$, resembling a Schrödinger cat state. Here, $\ket{{}^{\bullet}\,{}_{\bullet}}$ and $\ket{{}_{\bullet}\,{}^{\bullet}}$ indicate the ions’ non-classical anti-correlation, that can be identified with particles of momenta $+\hbar k_{a}$ and $-\hbar k_{a}$, respectively. We note that the creation of multiple particle pairs per state ($n=4,6,8$) that we observe evidences the bosonic character of the underlying quantum statistics of phonons  (?). By employing the criterion in  (?) for Gaussian states we quantify the entanglement of formation $E_{F}\approx 0.41$, significantly increased with respect to the intrinsic vacuum state entanglement $E_{F}\approx 10^{-5}$, see  (?). In future studies, we may be able to effectively decouple both ions, while preserving their non-classically correlation, and individually map the entanglement onto the ions’ electronic degree of freedom  (?, ?).

In general, many fundamental predictions rely on the premise that the laws of quantum theory remain valid, not only, across very different length and energy scales, but also during drastic changes of the external conditions. Here, we tested this premise in our trapped ion analog. We can continue these studies for generalized analogies involving scalar or vector fields in different dimensions. Depending on the interpretation of the time coordinate (e.g., transformation from proper to conformal time), we may simulate different dynamics of the related scale parameter  (?) and investigate its consequences. For this, we can further tune durations and shapes of the ramps of our potential or add parametric driving, e.g., to simulate inflation followed by a cosmological (p)re-heating phase corresponding to an oscillatory particle creation. Our platform may allow studying the causal connections of squeezing, pair creation, and entanglement, e.g., in the context of Hawking radiation due to surface gravity, the crossing of cosmic horizons during inflation, and the Sauter-Schwinger effect relating to high-field lasers  (?). Considering the controlled coupling of our closed quantum system to environments and noise fields might permit simulating realistically extended analogs. Further, our method represents a novel tool to create squeezing that can assist in gaining sensitivity for quantum metrology applications, cf.  (?). In addition, squeezing has also been proposed to substantially enhance effective spin-spin interactions, required for experimental simulations of quantum spin models  (?, ?), and has recently been used to implement qubits for QIP in the motional states of trapped ions  (?). Finally, our results regarding squeezing and purification from coherent excitations should be considered when implementing (multi-ion) entangling gates on multiplex trap architectures, where rapid changes of potential landscapes are required for scaling towards a universal quantum computer  (?, ?).

Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft (DFG) [237456450 and 278162697 (SFB 1242)].

Supplementary Materials

Materials and Methods

Supplementary Text

Figures S1 to S4

References *(37-47)

1 Analogies

In this Section, we provide a brief (and incomplete) list of scenarios for particle pair creation by tearing apart quantum vacuum fluctuations and show how they can be mapped onto harmonic oscillators with time-dependent frequencies which undergo large and rapid changes. We focus here on emphasizing analogies that might help to understand causal connections, similarities, and differences of the effect of particle creation in different contexts. For a detailed description, please refer to  (?). At the end of this Section we provide a summary in a dictionary, translating the physics and its corresponding terminology into the language of trapped ions.

Proper Time

In 3+1 dimensions, the expanding space-time during cosmic inflation can be described by the Friedmann-Lemaître-Robertson-Walker metric

[TABLE]

in terms of the proper time $\tau$, where the time-dependent scale parameter $a(\tau)$ governs the cosmic expansion.

After a spatial Fourier transform \Phi(\tau,\mbox{\boldmathr})\to\Phi_{\mbox{\boldmath\scriptstyle k}}(\tau) and a field redefinition \phi_{\mbox{\boldmath\scriptstyle k}}(\tau)=\Phi_{\mbox{\boldmath\scriptstyle k}}(\tau)\sqrt{a^{3}(\tau)}, the equation of motion $(\Box+m^{2}c^{2}/\hbar^{2})\Phi=0$ of a scalar field simplifies to

[TABLE]

In this representation, one can directly read off the stretching of physical wavelengths \lambda=2\pi a/|\mbox{\boldmathk}| due to the cosmic expansion in the term c^{2}\mbox{\boldmathk}^{2}/a^{2}.

The Hubble rate $H$ measuring the rapidity of the cosmic expansion is given by $H=\dot{a}/a$. For example, the de Sitter space-time is described by a constant Hubble rate $H$ and corresponds to an exponential growth of the scale parameter $a(\tau)\propto\exp(H\tau)$. In this case, the last term in the square brackets in (S2) is simply given by $9H^{2}/4$. Its competition with the contribution c^{2}\mbox{\boldmathk}^{2}/a^{2} from the internal (propagating) dynamics marks the event of horizon crossing: For wavelengths \lambda=2\pi a/|\mbox{\boldmathk}| much smaller than the Hubble radius $c/H$ (which determines the size of the cosmic horizon), the modes oscillate nearly freely. However, at some point in time the wavelength has been stretched enough by the cosmic expansion such that it exceeds the Hubble radius. Afterwards, the modes cannot oscillate anymore due to loss of causality and thus they cannot adapt to the changing external conditions, i.e., adiabaticity breaks down.

Conformal Time

In terms of the conformal time coordinate $t$ used in the manuscript, the 3+1 dimensional Friedmann-Lemaître-Robertson-Walker metric reads

[TABLE]

Again, a field redefinition \phi(t,\mbox{\boldmathr})=a(t)\Phi(t,\mbox{\boldmathr}) and a spatial Fourier transform yield

[TABLE]

Both Equations (S2) and (S4) are analog to harmonic oscillators with time-dependent frequencies, but the explicit expressions vary due to the different time coordinates. Note that also the explicit functional form of the time-dependence of the scale parameter $a$ changes with such a coordinate transformation. For example, the de Sitter space-time corresponds to $a(\tau)\propto\exp(H\tau)$ in proper time $\tau$, but to $a(t)\propto 1/t$ in conformal time $t$. Thus, the space-time simulated by the trapped ions depends on which time coordinate is identified with the laboratory time.

Vector Field

In contrast to the scalar (spin-zero) field $\Phi$ above, we may also consider a vector (spin-one) field $A^{\mu}$ such as a Proca field as described by the Proca equation of motion $(\nabla_{\mu}\nabla^{\mu}+m^{2}c^{2}/\hbar^{2})A^{\nu}=0$, where we have used the generalized Lorenz gauge $\nabla_{\mu}A^{\mu}=0$. After a mode decomposition A^{\mu}(t,\mbox{\boldmathr})\to A_{\mbox{\boldmath\scriptstyle k},\sigma}(t), where $\sigma$ labels the polarization, we arrive at (in terms of the conformal time)

[TABLE]

We see that the time-dependent scale factor $a(t)$ only enters via the mass term. As a result, photons are not created (as they are massless), which is a direct consequence of the conformal invariance of the electromagnetic field in 3+1 dimensions. However, one could create massive vector bosons such as $W^{\pm}$ or $Z^{0}$. This would require sufficiently large frequency scales in $a(t)$ which reach or even exceed the mass scale set by $m_{W}$ or $m_{Z}$. These large frequency scales in $a(t)$ arise naturally in the early universe, but one should check whether the above description in Eq. (S5) is a valid approximation in this period (symmetry breaking etc.). On the other hand, for Dirac (e.g., spin-1/2) fields $\Psi$, a quite analog mapping is possible. However, these are fermionic fields, while the vibrational excitations of the ions are of bosonic nature, i.e., the quantum statistics is then not directly captured.

Sauter-Schwinger Effect

As the next example, let us consider the Sauter-Schwinger effect, i.e., the creation of charged particle pairs by a strong electric field $E$. We assume a spatially homogeneous but time-dependent electric field \mbox{\boldmathE}(t) which can be represented by the vector potential \mbox{\boldmathA}(t) in temporal gauge via \mbox{\boldmathE}(t)=\mbox{\boldmath\dot{A}}(t). The equation of motion for a complex scalar field $\Phi$ reads (in flat space-time)

[TABLE]

where $q$ denotes the charge of the scalar field. A spatial Fourier transform simplifies this equation to

[TABLE]

Again we find the analogy to a harmonic oscillator with a time-dependent frequency.

In order to observe the non-perturbative Sauter-Schwinger effect associated to the strong-field regime, the electric field should be strong enough and slowly varying. This requirement can be quantified by means of the dimensionless (relativistic) Keldysh parameter

[TABLE]

where $\omega$ is a characteristic frequency scale of the electric field. Small Keldysh parameters $\gamma\ll 1$ correspond to the non-perturbative (strong-field) tunneling regime of the Sauter-Schwinger effect while large Keldysh parameters $\gamma\gg 1$ signify the perturbative (weak-field) multi-photon regime. Comparing Eqs. (S7) and (S8), we see that reaching the non-perturbative strong-field regime $\gamma\ll 1$ of the Sauter-Schwinger effect requires large variations of the harmonic oscillator frequency in Eq. (S7). In contrast, small changes of the harmonic oscillator frequency in Eq. (S7) such as in parametric down conversion in quantum optics (see also  (?)) correspond to the perturbative (weak-field) regime $\gamma\gg 1$.

Hawking Radiation

As our final example, let us consider Hawking radiation, i.e., black-hole evaporation. In terms of the Kruskal light-cone coordinate $U$, the black-hole metric in the (1+1 dimensional) near-horizon region reads

[TABLE]

where $\kappa$ is the surface gravity of the black hole and determines the Hawking temperature via $T_{\rm H}=\hbar\kappa/(2\pi k_{\rm B})$, see, e.g.,  (?). The term $e^{\kappa t}$ describes the tearing apart of waves by the strong tidal forces near the horizon – which is quite analog to particle creation by an expanding universe. Introducing the new time coordinate $T$ (not to be confused with the Hawking temperature $T_{\rm H}$) similar to the conformal time via $e^{-\kappa t}=-\kappa T$, the above metric reads

[TABLE]

where the global pre-factor $1/(\kappa T)^{2}$ now plays the role of the scale parameter $a^{2}(T)$.

The equation of motion for a scalar field $\Phi$ reads

[TABLE]

Finally, a Fourier expansion into modes $e^{iKU}$ and a field redefinition $\varphi_{K}(T)=e^{icKT}\phi_{K}(T)$ yields

[TABLE]

which is formally the same equation as in the manuscript. Thus, adjusting the time-dependence of the mode frequency according to $\omega^{2}_{K}(T)=c^{2}K^{2}+m^{2}c^{4}/(\hbar\kappa T)^{2}$ would mimic the tearing apart of quantum vacuum fluctuations near the black hole horizon by the strong gravitational field (as characterized by the surface gravity $\kappa$).

However, one should bear in mind that a full treatment of particle creation – especially in this scenario – must also include a consideration of the asymptotic regions (not just the vicinity of the horizon). Nevertheless, the tearing apart of quantum vacuum fluctuations is the most vital ingredient for Hawking radiation, see, e.g.,  (?). The entanglement generated in this process explains the thermal nature of the outgoing radiation: If we consider only one partner of the pair forming an entangled state, its density matrix corresponds to a mixed (thermal) state.

Dictionary

Let us discuss the analogy to the ion trap in more detail. Linearizing the radial displacements $\delta q_{i}$ of the ions (labeled by $i$), we obtain the equations of motion

[TABLE]

where $\omega^{2}_{\rm rad}(t)$ encodes the radial confinement controlled by the time-dependent trap potential and the matrix $M_{ij}$ describes the Coulomb interactions between the ions. Diagonalizing this matrix with the eigenvalues $\lambda_{I}$ (not to be confused with the wavelength) labeled by the index $I$ is then analog to the Fourier transform above and we obtain the mode equations

[TABLE]

which have a similar form as in the subsections above. The Coulomb interaction strength $\lambda_{I}$ is then analog to the internal (propagating) dynamics c^{2}\mbox{\boldmathk}^{2}. In both cases, they are responsible for generating entanglement. The trap potential $\omega^{2}_{\rm rad}$ simulates the mass term $m^{2}c^{4}/\hbar^{2}$ plus possible potential terms $\propto\dot{a}^{2}$ or $\propto\ddot{a}$ as in Eq. (S2). Its time-dependence $\omega^{2}_{\rm rad}(t)$ represents the external influence and is then analog to the cosmic expansion. As shown above, the radial displacements $\delta q_{i}$ of the individual ions correspond to the field $\Phi$ (or $A^{\mu}$ etc.) in position representation while the modes $\delta q_{I}$ (such as the center-of-mass or rocking mode) are analog to the Fourier modes \phi_{\mbox{\boldmath\scriptstyle k}}.

Holding the axial positions of the ions fixed, the Coulomb interaction strength $\lambda_{I}$ stays constant. Thus, when identifying $\lambda_{I}$ with the internal dynamics c^{2}\mbox{\boldmathk}^{2}, Eq. (S14) is most similar to the evolution equation (S5) of a vector field in 3+1 dimensions, or Eq. (1) of a scalar field in 1+1 dimensions, where the additional potential terms $\propto\dot{a}^{2}$ or $\propto\ddot{a}$ are absent. This again shows that the explicit details of the equations of motion – albeit they are all qualitatively similar – depend on the field (e.g., $\Phi$ or $A^{\mu}$), the number of dimensions and the used time coordinate etc. For the vector field $A^{\mu}$, we could identify the two transversal polarizations $\sigma$ with the two radial directions.

Identifying the harmonic oscillator potential $\omega^{2}(t)$ with c^{2}\mbox{\boldmathk}^{2}+a^{2}(t)m^{2}c^{4}/\hbar^{2} as in the manuscript and in Eq. (S5), we see that a relative change of $\omega(t)$, such as a certain number of $e$-foldings, implies an even larger relative change of $a(t)$, i.e., more $e$-foldings, because the remaining constant term c^{2}\mbox{\boldmathk}^{2} is positive. Thus, the reported 1.6 $e$-foldings of $\omega(t)$ imply, in principle, more than 1.6 $e$-foldings of space $a(t)$.

Pairs of phonons (as vibrational excitation quanta) of the ions represent pairs of particles as lasting excitations of the field $\Phi$ (or $A^{\mu}$ etc.). Finally, the entanglement of distant space points created during inflation corresponds to the motional entanglement of the two ions. After separating the ions and reading out the motional states individually, the entanglement would manifest itself as a thermal (i.e., mixed) state with a temperature set by the squeezing parameter – in complete analogy to Hawking radiation which appears thermal because one can only observe the one partner of the particle pair escaping to infinity (the other one is falling into the black hole).

2 Experimental methods

Our experimental setup features a linear radio-frequency (rf) ion trap and most of the relevant techniques employed here are described in detail in  (?, ?). We trap and manipulate different isotopes of Mg+ in a common electric potential. Radial confinement is dominated by an rf potential oscillating at $\Omega_{\text{rf}}/(2\pi)\approx$56\text{,}\mathrm{MHz}$$, while fine tuning along the radial directions is accomplished by dedicated control voltages applied to six control electrodes. In addition, we use these electrodes to adjust the axial confinement, to which rf fields contribute effectively less than a few percent. Our rf generator setup driving the two rf electrodes includes a signal generator with fixed amplitude, an rf mixer controlled via an arbitrary waveform generator  (?) (with output amplitude $U_{\text{c}}$), a high power amplifier (output power $\approx 40$ dBm), a bandpass filter ($800\text{\,}\mathrm{kHz}$ bandwidth) and a helical resonator with a loaded quality factor $Q\approx 100$. With our real-time data acquisition system, with a $10\text{\,}\mathrm{ns}$ timing resolution, we tune and ramp the applied rf voltage amplitude $U$ between $0\text{\,}\mathrm{V}$ and $800\text{\,}\mathrm{V}$.

We use a single ${}^{25}\text{Mg}^{+}$ in our trap to calibrate and optimize control parameter settings, i.e., tuning of offset electric fields and curvatures, while for the main sequences (main text), we add one ${}^{26}\text{Mg}^{+}$ and form a strongly coupled hybrid crystal. We choose the ${}^{25}\text{Mg}^{+}$ to implement a pseudo spin system (internal/electronic degree of freedom) on two suitable S1/2 hyperfine ground states, labeled $\ket{\downarrow}$ and $\ket{\uparrow}$. We employ Doppler cooling and subsequent resolved sideband cooling via two-photon stimulated-Raman transitions  (?, ?, ?, ?) applied to the ${}^{25}\text{Mg}^{+}$. In this way, we cool up to six motional modes (external degrees of freedom) and yield low mean thermal occupation numbers. We employ established techniques to map motional state distributions to the spin states, and electron shelving to distinguish $\ket{\downarrow}$ (bright state) from $\ket{\uparrow}$ (dark state), cf.  (?). The mapping technique is based on the coherent coupling of the internal and external degrees of freedom via resonant transitions $\ket{\downarrow,n}\leftrightarrow\ket{\uparrow,n+1}$ (blue sideband) and $\ket{\uparrow,n}\leftrightarrow\ket{\downarrow,n-1}$ (red sideband), where $\ket{n}$ denotes the corresponding number state of the individually addressed mode.

In a first series of calibration measurements, we record the variation of motional frequencies as a function of $U$, results are shown in Fig. S1A. We distinguish two radial modes, the high-frequency (hf) mode $\omega_{\text{hf}}$ and the mid-frequency (mf) mode $\omega_{\text{mf}}$ (blue squares), from the axial (low-frequency, lf) mode $\omega_{\text{lf}}$ (green diamonds), solid lines depict model fits to the data. In addition, we estimate the frequency variation of the rocking mode $\omega_{a}\approx\sqrt{\omega_{\text{mf}}^{2}-\omega_{\text{lf}}^{2}}$ (orange circles). During the course of our experimental runs, we monitor changes of $U$ with a dedicated pickup electrode that is capacitively coupled to the rf electrodes, and the corresponding pickup signal amplitude is denoted by $U_{p}$. Typically, we program a sixth-order smooth-step waveform $U_{c}(t)$ to modulate $U$, with fixed ramp duration $t_{\text{ramp}}=$1\text{,}\mathrm{\SIUnitSymbolMicro s}$$ (cf. Fig. 2B), pulsing between high amplitude $U_{\text{high}}$ and low amplitude $U_{\text{low}}$. We record $U_{p}(t)$ and the input waveform $U_{c}(t)$ to characterize the real-time switching of $U$. In Figure S1B, we depict corresponding results of a single pulse (two ramps), taking into account the calibration measurements in Fig. S1A. Due to the finite bandwidth of our rf generator setup (including all wiring), the waveform $U_{p}$ is shifted and distorted with respect to $U_{c}$. We account for the shift by an appropriate delay $t_{\text{delay}}$ in our experimental sequences. While the distortion can affect the (quantitative) details of the mode dynamics, general (qualitative) features are robust and final squeezing amplitudes vary by a few percent only, see below.

In a second series of calibration runs, we benchmark the real-time switching capabilities of our setup with a spin-motional echo sequence  (?) on the mf mode with an interferometric duration (arm length) $t_{\text{se}}=2\,t_{\text{ramp}}+t_{\text{hold}}+t_{\text{delay}}$, see Fig. S1C. We initialize the state $\ket{\downarrow,n_{\text{mf}}=0}$ and, subsequently, create a coherent superposition $\ket{\downarrow,0}+\ket{\uparrow,1}$ with a $\pi/2$-pulse on the blue sideband. During the first arm, we ramp down $U$, hold at $U_{\text{low}}$ for a duration $t_{\text{hold}}$, and ramp back up to $U_{\text{high}}$. Note, during $t_{\text{se}}$ the spin-motional superposition state accumulates a phase that is proportional to $\int_{0}^{t_{\text{se}}}\omega_{\text{mf}}(t)dt$. Following an echo pulse, i.e., a $\pi$-pulse on the blue sideband, $U$ is kept high (no pulse) during the second arm and the sequence concludes with another $\pi/2$-pulse on the blue sideband, and spin-state detection. Overall, this sequence is sensitive to differential phase accumulations of both arms, i.e., mode frequency differences can be detected. Repeating measurements for variable $t_{\text{hold}}$ and $t_{\text{se}}$, accordingly, we record final $P_{|\downarrow\rangle}$ and show corresponding results in Fig. S1D. From a sinusoidal model fit to the data we extract a mode frequency difference of $2\pi\times$1.53(3)\text{,}\mathrm{MHz}, in agreement with the static calibration values $\omega_{\text{mf}}(U_{\text{high}})-\omega_{\text{mf}}(U_{\text{low}})=2\pi\times$1.48(1)\text{\,}\mathrm{MHz}, cf. Fig. S1A.

In a third series of calibration measurements, we estimate motional heating rates  (?) for the two-ion configuration. In particular, we probe the mid-frequency mode and the accompanying rocking mode for $U_{\text{high}}$. We record rates of $\lesssim$2000\text{,}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{/}\mathrm{s} for the center-of-mass and $\lesssim$20\text{\,}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{/}\mathrm{s} for the rocking mode, while heating rates for the other four modes of the crystal are comparable or smaller. For total sequence lengths of smaller than $200\text{\,}\mathrm{\SIUnitSymbolMicro s}$ in our main experiments, we can neglect incoherent (thermal) mode excitations on the rocking mode. Note, even when conservatively estimating a $60$-fold increased rate for $U_{\text{low}}$, the rocking mode is heated less than $\approx 0.02$ quanta during our sequences. This reflects the fact that out-of-phase modes are intrinsically resilient against homogeneous noise fields.

Finally, we apply our one pulse sequences with variable offset electric fields and curvatures to minimize coherent displacement contributions, i.e., compensating (static) stray electric fields. In our experiments, residual stray fields displace the ions from the potential minimum to an equilibrium position given by the restoring force of the potential $\propto\omega^{2}$. If $\omega$ is switched, the instantaneous equilibrium position is shifted and, thus, the ions experience an effective force that excites (unwanted) coherent motion. While coarse tuning is performed with single ions, fine tuning is done on the hybrid crystal with an axial ion separation of $\approx$5.5\text{,}\mathrm{\SIUnitSymbolMicro m}. From corresponding results, we deduce residual, homogeneous stray fields of $1\text{\,}\mathrm{V}\mathrm{/}\mathrm{m}$ along the mid frequency mode. Stray fields acting differentially on each ion, i.e., stray curvatures, are smaller than $5\cdot 10^{-8}\,$\text{\,}\mathrm{V}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}.

3 Motional state reconstruction

The analysis for motional state reconstruction consists of two crucial steps: (I) Using well-established techniques to determine population distributions (data points in Fig. 3B, C and Fig. 4D) and (II) to yield estimates for squeezing and displacement excitations. For the first step, we employ red or blue sideband couplings to the relevant rocking mode in individual sequences for variable durations $t_{\text{sb}}$, and record $P_{|\downarrow\rangle}$. Figure S2 depicts the experimental data for the states depicted in Figs. 3 and 4 where each data point is the result of 200 repetitions. We reconstruct the motional state distribution $P_{n}$ (population of each number state, data points in Figs. 3B-C and 4D) from combined model fits to the probability $P_{|\downarrow\rangle}(t)$ (solid lines), see  (?). Note, sideband coupling (Rabi) rates depend on $\ket{n}$ and motional state distributions lead to characteristic evolutions of $P_{|\downarrow\rangle}$, and we choose to consider up to $n=8$ in the reconstruction. Independent from any further considerations, we note that, in particular, data points shown in Fig. 4D exhibit direct evidence for substantial squeezing and negligible displacement and thermal excitations. Only even number states are significantly occupied (particle pairs), while populations are in agreement with zero for odd states (unwanted excitations).

In the second step of our analysis, we estimate contributions of squeezing $r$ and displacement $\alpha$ (due to residual stray electric fields), while heating effects are measured to be small in our case (see above) and neglected in the following. To this end we define creation (annihilation) operators $\hat{a}^{\dagger}$ ($\hat{a}$) and combine a squeezing operator $\hat{S}(\xi)=\exp\left(\frac{1}{2}\left(\xi^{*}\hat{a}^{2}-\xi\hat{a}^{\dagger 2}\right)\right)$, with $\xi=re^{i\theta}$, and a displacement operator $\hat{D}(\alpha)=\exp\left(\alpha\hat{a}^{\dagger}-\alpha^{*}\hat{a}\right)$ to create a squeezed coherent state $\rho_{\text{fin}}$ by

[TABLE]

Here, $\rho_{\text{ini}}(\bar{n}_{\text{th}})$ denotes the initial (thermal) density matrix. The diagonal elements of $\rho_{\text{fin}}$ represent the elements $P_{n}^{\text{par}}(r,\alpha,\bar{n}_{\text{th}})$, which we fit to the measured $P_{n}$ to estimate the final experimental state. In particular, we allow for a differential phase between squeezing and displacement to reduce systematic deviations. From $\rho_{\text{fin}}$ the amplitudes $r$ and $|\alpha|$ can be extracted similarly to the numerically simulated states (see below). We note that it is also possible to reconstruct the full density matrix, i.e., Wigner function, experimentally  (?), from which, in turn, squeezing and coherent displacements could be distinguished directly.

4 Numerical simulations

In the following, we perform numerical simulations to gain an additional perspective on the underlying systematic effects in the relevant dynamics, in particular, the effect of finite ramp durations and the accuracy of $P_{n}^{\text{par}}(r,\alpha,\bar{n}_{\text{th}})$. We consider a single harmonic oscillator with time-dependent eigenfrequency $\omega(t)$ and recast the corresponding Hamiltonian  (?):

[TABLE]

Here, $\hat{a}_{\omega}^{\dagger}(t)$ and $\hat{a}_{\omega}(t)$ denote time-dependent creation and annihilation operators, that are transformed into time-independent basis operators $\hat{a}^{\dagger}$ and $\hat{a}$ defined for the initial mode frequency $\omega_{0}=\omega(0)$. Further, we include the effect of stray electric fields by an additional term:

[TABLE]

where $F_{0}(t)$ denotes an effective stray force amplitude and $x_{0}=\sqrt{\hbar/\left(2m\omega_{0}\right)}$ is the width of the ground state wave function of the oscillator with mass $m$. In addition, we define position $\hat{x}=\sqrt{\hbar/\left(2m\omega_{0}\right)}\,\left(\hat{a}^{\dagger}+\hat{a}\right)$ and momentum $\hat{p}=i\sqrt{\hbar m\omega_{0}/2}\,\left(\hat{a}^{\dagger}-\hat{a}\right)$ operators, corresponding expectation values $x(t)=\langle\hat{x}\rangle$ and $p(t)=\langle\hat{p}\rangle$, and variances $\left(\Delta x\right)^{2}(t)=\langle\hat{x}^{2}\rangle-\langle\hat{x}\rangle^{2}$ and $\left(\Delta p\right)^{2}(t)=\langle\hat{p}^{2}\rangle-\langle\hat{p}\rangle^{2}$. The squeezing and displacement amplitudes are then given by $r(t)=\frac{1}{2}\operatorname{arccosh}\left[\left(\Delta x\right)^{2}(t)+\left(\Delta p\right)^{2}(t)\right]$ and $|\alpha(t)|=\sqrt{x^{2}(t)+p^{2}(t)}$, respectively. We numerically simulate the total Hamiltonian $\hat{H}+\hat{H}_{F}$ using QUTIP  (?) and obtain unitary time-evolutions of the density matrix, from which we can calculate all of the defined expectation values, neglecting any decoherence effects.

First, we assess the effect of variable ramp durations and depict final $r$ for the single ramp, with fixed stray field amplitude, from low frequency $\omega_{\text{low}}=\omega_{\text{a}}(U_{\text{low}})$ to high frequency $\omega_{\text{high}}=\omega_{\text{a}}(U_{\text{high}})$ in Fig. S3A. For shorter $t_{\text{ramp}}$, squeezing approaches the ultimate value, corresponding to an instantaneous switching, given by $r_{\text{max}}=\frac{1}{2}\ln\left(\frac{\omega_{\text{high}}}{\omega_{\text{low}}}\right)\approx 0.8$. Assuming a variation of $10\text{\,}\mathrm{\char 37\relax}$ of the ramp duration around $1\text{\,}\mathrm{\SIUnitSymbolMicro s}$, the squeezing amplitude varies by less than $10\text{\,}\mathrm{\char 37\relax}$. In comparison, we illustrate variable contributions of final $|\alpha|$ as a function of $t_{\text{ramp}}$, for fixed stray field amplitudes, in Fig. S3B.

Second, we study the accuracy of $P_{n}^{\text{par}}(r,\alpha,\bar{n}_{\text{th}})$ for the case of our two-pulse sequence. To this end, we insert corresponding experimental parameter settings, determined by the above described calibration measurements, into the numerical simulations. In Figure S4, we show final numerical $r$ and $|\alpha|$ for variable changes of the free evolution duration $\Delta t_{\text{free}}$, including a fixed stray field (amplitude $11\text{\,}\mathrm{m}\mathrm{V}\mathrm{/}\mathrm{m}$). We fit numerical $P_{n}(\Delta t_{\text{free}})$ with $P_{n}^{\text{par}}(r,\alpha,\bar{n}_{\text{th}})$ to extract final $r$ (dashed orange line) and $|\alpha|$ (dashed gray line) and compare it to exact numerical results (solid lines). In Figure S4B we plot corresponding residuals and find overall agreement with residual deviations on the order of a few percent.

5 Spatial entanglement of two ions

We present a summarizing view on the spatial entanglement in our bipartite system, ion A and B, that builds on original work of  (?, ?) and consider the specific condition in our experiment. Ions A and B are confined in three spatial degrees of freedom, interact individually with the common storage field, and with each other due to their mutual Coulomb interaction. We focus on the description of one spatial degree of freedom, e.g., the weaker radial direction, and define individual, single ion oscillator states $\ket{n}_{\text{A, B}}$ with corresponding (uncoupled) eigenfrequencies $\omega_{\text{A}}=\omega_{\text{B}}$. Due to the ion-ion interaction, both harmonic oscillators are coupled with rate $\omega_{\text{coupl}}$. We distinguish short timescales $t\ll 1/\omega_{\text{coupl}}$ from long timescales $t\gg 1/\omega_{\text{coupl}}$, where the motion is typically described by a common in-phase and out-of-phase motion with eigenfrequencies $\omega_{\text{+, -}}$ and corresponding number states $\ket{n}_{\text{+, -}}$, respectively. Further, we assume an initial, vacuum state  (?):

[TABLE]

with $e^{-\beta}=\sqrt{(\kappa-1/2)/(\kappa+1/2)}$ and $\kappa=1/4\left(\sqrt{\omega_{+}/\omega_{-}}+\sqrt{\omega_{-}/\omega_{+}}\right)$. In the basis of $\ket{n}_{\text{A, B}}$, this corresponds to an entangled, squeezed state and, for our experimental parameters $\kappa\approx 0.5008$. In our experiments, we increase the spatial entanglement by the coherent generation of squeezing in the out-of-phase mode, while squeezing of the in-phase mode is, estimated from numerical simulations, negligible. We can separate the residual (classical) coherent displacement from the non-classical dynamics and, for illustration, we assume pure motional squeezing $\hat{S}(\xi)$ to estimate our final state:

[TABLE]

Note, the second term is related to our state illustration $\ket{{}^{\bullet}\,{}_{\bullet}}+\ket{{}_{\bullet}\,{}^{\bullet}}$ given in the main text.

To quantify the amount of entanglement, we, further, consider that the ions are initially prepared in a thermal state described by average phonon numbers $\bar{n}_{\text{+, -}}$ and states remain Gaussian throughout our experimental sequences. The entanglement of formation of the ions is then given by  (?, ?, ?, ?):

[TABLE]

where $\chi$ is related to the so-called symplectic eigenvalues:

[TABLE]

via $\chi=\min(\lambda_{1},\lambda_{2})$. The symplectic eigenvalues and, therefore, $E_{F}$ are functions of the applied squeezing $r$, i.e., $E_{F}(r=0)\approx 10^{-5}$ for our initial state and $E_{F}(r=0.83)\approx 0.41$ for our final state in Fig. 4D. Note, Equations (S24) and (S25) are a generalization of results obtained in  (?), which take the impact of finite coupling strength $\omega_{+}/\omega_{-}$ into account. For pure states and in the absence of squeezing, i.e, $\bar{n}_{-}=\bar{n}_{+}=r=0$, the entanglement of formation Eq. (S21) reduces to the von Neumann entanglement as stated in  (?).