Universal quantum computation and quantum error correction with ultracold atomic mixtures

In this section, we discuss the fundamental unit of information of our proposal: a collective spin with controllable length. In order to realize such a collective spin we consider bosonic atoms of type A, which are localized at the local minima y of the optical potential, e.g. arrays of optical tweezer or an optical lattice (see figure 1(a)). The annihilation (creation) operators of the atoms on site y are a_m (y) and ${a}_{m}^{{\dagger}}(\mathbf{y})$, where m indicates the two internal states (0, 1) of the atom. The local confining potential is chosen such that there is no hopping. Hence, the atom number per site ${N}_{\mathrm{A}}(\mathbf{y})={a}_{1}^{{\dagger}}(\mathbf{y}){a}_{1}(\mathbf{y})+{a}_{0}^{{\dagger}}(\mathbf{y}){a}_{0}(\mathbf{y})$ is conserved [47] and the atoms form a collective spin via the Schwinger representation:

$$\begin{equation}{L}_{z}(\mathbf{y})=\frac{1}{2}[{a}_{1}^{{\dagger}}(\mathbf{y}){a}_{1}(\mathbf{y})-{a}_{0}^{{\dagger}}(\mathbf{y}){a}_{0}(\mathbf{y})],\end{equation} \tag{ 1a }$$

$$\begin{equation}{L}_{+}(\mathbf{y})={a}_{1}^{{\dagger}}(\mathbf{y}){a}_{0}(\mathbf{y}),\end{equation} \tag{ 1b }$$

$$\begin{equation}{L}_{-}(\mathbf{y})={a}_{0}^{{\dagger}}(\mathbf{y}){a}_{1}(\mathbf{y}),\end{equation} \tag{ 1c }$$

with L_x (y) = [L₊(y) + L₋(y)]/2 and L_y (y) = (−i)[L₊(y) − L₋(y)]/2. The conservation of atom number per site determines also the magnitude of the angular momentum ℓ = N_A/2. The eigenstates of the L_z operator are denoted by |m_ℓ ⟩ with m_ℓ being an integer or half-integer ranging from m_ℓ = −ℓ, ..., ℓ. By defining the computational basis |j⟩ ≡ |−ℓ + j⟩ with j = 0, ..., N_A a collective spin can be interpreted as a qudit with dimension d = N_A + 1. Moreover, we can access the qubit (N_A = 1), qudit (N_A > 1), and continuous variable (N_A → ∞) regime by tuning the atom number per site. In particular, the Hilbert space dimension of a single collective spins scales linearly with the numbers of atoms.

In order to realize unitary operations acting on a single collective spin, we consider the time evolution generated by the Hamiltonian

$$\begin{equation}{H}_{\mathrm{A}} = \sum\limits _{\mathbf{y}}\left[\chi (\mathbf{y}){L}_{z}^{2}(\mathbf{y})+{\Delta}(\mathbf{y}){L}_{z}(\mathbf{y})+{\Omega}(\mathbf{y}){L}_{x}(\mathbf{y})\right] .\end{equation} \tag{ 2 }$$

As detailed in appendix A, this Hamiltonian can be implemented by a two component Bose gas localized in optical tweezers, where the first term corresponds to the interaction between the atoms, the second is due to the presence of the magnetic field, and the third term to a Rabi coupling between the hyperfine states. The interactions between the atoms are described by ${L}_{z}^{2}(\mathbf{y})=\frac{1}{4}\left[{n}_{1}^{2}(\mathbf{y})+{n}_{0}^{2}(\mathbf{y})-2{n}_{0}(\mathbf{y}){n}_{1}(\mathbf{y})\right]$, which involve both spin-states, but no spin-changing collisions.

The Hamiltonian in equation (2) can generate all unitary operations on a single collective spin since all couplings χ(y), Δ(y), and Ω(y) can be switched on and off independently on each site in ultracold atom systems [48]. For the F = 1 hyperfine manifold, the coupling χ(y) = 0 can be achieved by hosting the atoms in the magnetic substates m_F = 1 and m_F = −1, while χ(y) ≠ 0 when the atoms are in the magnetic substates m_F = 1 and m_F = 0. The switching of χ(y) was illustrated in the context of mixture experiments in a recent work using Raman transitions [49]. A great advantage of using laser addressing is that the magnetic field does not have to be controlled at the μm scale, which is the typical separation of optical tweezers. Similarly, the couplings Ω(y) and Δ(y) can be tuned via Raman coupling in the same way as in the high-fidelity implementations with trapped ions [50].

For the qubit case (d = N_A + 1 = 2), the operators L_x and L_z together with the identity 1 suffice to obtain all U(2) transformations [51]. For qudits (d = N_A + 1 > 2) the additional non-linear operation ${L}_{z}^{2}\left.\right)$ is required to approximately synthesize any unitary U(d) via Trotterization. Namely, using ${\text{e}}^{-\text{i}\hat{A}\delta t}\enspace {\text{e}}^{\text{i}\hat{B}\delta t}\enspace {\text{e}}^{\text{i}\hat{A}\delta t}\enspace {\text{e}}^{-\text{i}\hat{B}\delta t}={\text{e}}^{[\hat{A},\hat{B}]\delta {t}^{2}}+O\left(\delta {t}^{3}\right)$, we can engineer all possible commutators of L_x , L_z , and ${L}_{z}^{2}$ and higher order commutators. These commutators span the d-dimensional Hermitian matrices with the basis ${\left\{{\mathcal{M}}_{i}\right\}}_{i=1}^{{d}^{2}}$, which in turn allows to synthesize all unitary matrices. In appendix B we give an explicit construction of the enveloping algebra of L_x , L_z , and ${L}_{z}^{2}$ for d = 3, which allows one to generate each element of U(3). As detailed in reference [52], the construction given in the appendix can be generalized to d > 3. In this way, the operators L_x , L_z , and ${L}_{z}^{2}$ provide universal control over each single collective spin. In the next section, we consider the entanglement of several collective spins to achieve exponential growth of the Hilbert space.

The entanglement of different collective spins can be achieved by the exchange of delocalized phonons. The phonons are the Bogoliubov excitations of a weakly interacting condensate of atomic species B in n dimensions, where we consider n = 1, 2 (see figure 1). The purely phononic part of the Hamiltonian is given by

$$\begin{equation}{H}_{\mathrm{B}}=\sum\limits _{\mathbf{k}}\;\hslash {\omega }_{\mathbf{k}}{b}_{\mathbf{k}}^{{\dagger}}{b}_{\mathbf{k}},\end{equation} \tag{ 3 }$$

where b_k is the annihilation operator of a phonon mode at wave number k with components ${k}_{i}=\frac{\pi }{L}{m}_{i}$ and m_i is a positive integer. The phonon dispersion at low momenta is linear ω_k ≈ c|k|, where $c=\sqrt{{\tilde{g}}_{\mathrm{B}}{n}_{\mathrm{B}}/{M}_{\mathrm{B}}}$ is the speed of sound determined by the interaction strength ${\tilde{g}}_{\mathrm{B}}$ of the B atoms with mass M_B and the density n_B, see appendix C for more details.

By immersing the collective spins into the phonon bath, we induce interactions between the phonons and the spins rooted in the contact interactions of the two atomic species A and B. The spin–phonon interaction can be described by the Hamiltonian

$$\begin{equation}{H}_{\mathrm{A}\mathrm{B}}=\sum\limits _{\mathbf{y},\mathbf{k}}\left[{\bar{g}}_{\mathbf{k}}(\mathbf{y})+\delta {g}_{\mathbf{k}}(\mathbf{y}){L}_{z}(\mathbf{y})\right]({b}_{\mathbf{k}}+\text{H.c.}),\end{equation} \tag{ 4 }$$

with the explicit forms of the coupling constants ${\bar{g}}_{\mathbf{k}}(\mathbf{y})$ and δg_k (y). A detailed derivation of equation (4) is given in appendix D. This interaction is similar to the phonon–ion interactions in trapped ion systems [16] or photon–atom interactions [53]. The first term of equation (4) leads to a constant polaronic shift [45, 54], which can be absorbed by redefining the phonon operators. Consequently, we focus here on the second term, which can be used to generate entanglement between the spins.

Since L_z (y) is a conserved quantity for Ω(y) = 0 it is possible to decouple the phonons and spins in equation (4) by shifting the phonon operators (see appendix E). Eliminating the phonons leads to an effective long-range spin–spin interaction

$$\begin{equation}{H}_{I}=-\sum\limits _{\mathbf{x},\mathbf{y}}\;g(\mathbf{x},\mathbf{y}){L}_{z}(\mathbf{x}){L}_{z}(\mathbf{y}),\end{equation} \tag{ 5 }$$

where we introduced the coupling g(x, y) between the spins after a proper redefinition of χ(y) and Δ(y).

Explicitly, the coupling between the spins is given by

$$\begin{equation} g(\mathbf{x},\mathbf{y})\;=\;g\sum\limits _{\mathbf{k}} \frac{{({u}_{\mathbf{k}}+{v}_{\mathbf{k}})}^{2}}{\hslash {\omega }_{\mathbf{k}}}\enspace {\mathrm{e}}^{-\frac{1}{4}\vert \mathbf{k}{\vert }^{2}{\sigma }_{ A}^{2}} \prod\limits _{i=1}^{n}\mathrm{sin}({k}_{i}{x}_{i})\mathrm{sin}({k}_{i}{y}_{i}).\end{equation} \tag{ 6 }$$

The overall prefactor

$$\begin{equation}g={n}_{\mathrm{B}}{\left(2/L\right)}^{n}{({\tilde{g}}_{\mathrm{A}\mathrm{B}}^{1}-{\tilde{g}}_{\mathrm{A}\mathrm{B}}^{0})}^{2}\end{equation} \tag{ 7 }$$

is determined by the inter-species interaction ${\tilde{g}}_{\mathrm{A}\mathrm{B}}^{0}$ and ${\tilde{g}}_{\mathrm{A}\mathrm{B}}^{1}$, where the ˜ indicates the renormalization according to the optical potential. Further, we introduce the Bogoliubov amplitudes u_k and v_k , which are determined by performing the Bogoliubov approximation in a box potential of length L [55, 56]. We further approximate the Bogoliubov mode functions by sin(k_i x_i ). The length scale σ_A is the harmonic oscillator length of the tweezer confining the A atoms, which provides a cutoff for the momentum sum in equation (6). For more details we refer to appendix D.

Performing the momentum sum numerically, the resulting interaction between two spins for one and two dimensions is illustrated in figures 2 and 3. Positioning the spins appropriately within the bath allows one to tune the strength and sign of the spin–spin interaction. The effective interaction scales as n_B independent of the dimension, as long as only the linear part of the dispersion ω_k contributes to the sum of equation (6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to entangle two specific spins, the ability to deliberately switch on and off the interactions g(x, y) is necessary. A first possibility is to physically move the optical tweezers such that there is no overlap of the collective spins with the phonons. This approach is similar to the shuttling approach in trapped ions [57].

The second approach, is similar to the optical shelving used in trapped ions [58]. The interaction between the collective spins is proportional to the scattering length difference, $g(\mathbf{x},\mathbf{y})\propto {\left({\tilde{g}}_{\mathrm{A}\mathrm{B}}^{1}-{\tilde{g}}_{\mathrm{A}\mathrm{B}}^{0}\right)}^{2}$. The first term describes the interaction strength between A atoms in m_F = 1 and B atoms in m_F = 0, while the second term describes the interaction between atoms of both species being in m_F = 0. To shelve spins, we can coherently transfer the m_F = 0 component of A atoms into the m_F = −1 state, while the atoms in the m_F = 1 component are left unaltered. The resulting interaction is strictly zero since $g(\mathbf{x},\mathbf{y})\propto {\left({\tilde{g}}_{\mathrm{A}\mathrm{B}}^{-1}-{\tilde{g}}_{\mathrm{A}\mathrm{B}}^{1}\right)}^{2}=0$.

By tuning the interaction time between the spins, using two additional single qudit gates and adding a phase allow for synthesizing a controlled-Z gate [59] between two collective spins

$$\begin{equation}CZ(\mathbf{x},\mathbf{y})\;=\;{\text{e}}^{\text{i}\frac{2\pi }{d}{L}_{z}(\mathbf{x}){L}_{z}(\mathbf{y})}\enspace {\text{e}}^{\text{i}\frac{2\pi \ell }{d}{L}_{z}(\mathbf{x})}\enspace {\text{e}}^{\text{i}\frac{2\pi \ell }{d}{L}_{z}(\mathbf{y})}\enspace {\text{e}}^{\text{i}\frac{2\pi {\ell }^{2}}{d}},\end{equation} \tag{ 8 }$$

which completes the universal gate set in the multi-spin system [60, 61]. The quality of the local spin addressing and the magnetic field stability [62] will determine the fidelity of the gates.

A universal quantum computer requires a reliable state-preparation and readout. Using the external magnetic field one can prepare all A atoms in one hyperfine component, which corresponds to a fully polarized collective spin along the quantization axis |ψ⟩ = |d⟩ ⊗ ...|d⟩. In case the particle number per site is probabilistic, one can perform post-selection to fix N_A [63–65]. If one uses instead an optical lattice to create the collective spins, one can deterministically control the atom number by preparing first a Mott-insulating state, which nowadays can be prepared with almost unit filling [26] and a subsequent merging of a fixed number of wells. For the preparation as well as for the detection one has to ensure single counting statistics per site, which can be achieved through fluorescence imaging [34].

When working in the large collective spin regime, one can efficiently prepare condensates through evaporative cooling in the optical dipole trap [66]. This results in an atomic cloud of a few hundred to thousand atoms of low entropy in each dipole trap. In case the A atoms are not in the motional ground state, one can perform Raman sideband cooling [67]. For larger collective spins, one can use homodyne detection to map out large collective spins [68].

Having established all the necessary ingredients for a quantum information platform, we discuss details of a possible experimental realization as well as main sources of errors. The experiments can rely on standard experimental tools employed for cold atoms [20]. For concreteness, we focus on a mixture of ³⁹K (species A) forming the collective spins and ²³Na (species B) forming the phonon bath.

The collective spins can be realized through the trapping of a few atoms in optical tweezers distanced by a few micrometers and with motional ground state extension σ_A ≈ 100 nm. With this confinement and working with ≈10 atoms per tweezer leads to three body loss between the species [69] on a time scale of ∼1 Hz. The typical single qubit gates can then be performed with a Rabi frequency of up to a few kHz with above 99% fidelity [38]. In order to ensure the coherent evolution one has to stabilize the magnetic field, which for state of the art experiments is possible up to ∼10 Hz. Readout reliability through fluorescence imaging differ between schemes but can reach >99% fidelity [70–72].

Assuming a box potential for species B [55] in a tube or a slab geometry with linear dimension L ≈ 100 μm we expect approximately 50 collective spins in a one-dimensional geometry or 2500 collective spins in a two dimensional geometry. In both cases, all-to-all coupling can be achieved through phonons of a (quasi-)condensate of approximately N_B ≈ 300 × 10³ atoms for the 1D bath and N_B ≈ 3 × 10⁶ for the 2D bath, with a transverse confinement of ω_⊥ ≈ 2π × 440 Hz, such that the condensate has a chemical potential of μ_B/h ≈ 7.7 kHz in 1D and μ_B/h ≈ 1.9 kHz in 2D [22]. This implies typical lifetimes for single atom losses that can be of several seconds up to a few minutes, setting an upper limit on the length of the quantum circuit [25].

The collective spins are coupled to the phonons through contact interaction according to the scattering length ${a}_{0}^{\mathrm{A}\mathrm{B}}=756{a}_{0}$, ${a}_{1}^{\mathrm{A}\mathrm{B}}=2542{a}_{0}$ and ${a}_{2}^{\mathrm{A}\mathrm{B}}=-437{a}_{0}$. The couplings ${g}_{\mathrm{A}\mathrm{B}}^{m}$ can then be obtained by calculating the corresponding Clebsch–Gordon coefficients [73]. The speed of the entangling gates is then determined by the strength of the interaction between species A and B given by equation (6), as we illustrated in figures 2 and 3. For example, the typical gate speed of the CZ gate connecting two qudits at a distance of 5 μm can be estimated to be of the order of 50 Hz in the one-dimensional case, see figure 2. Using an optical lattice for the collective spins, 20 × 20 lattice sites are within the current experimental feasibilities and allow all qudits to interact with the phonons. The coupling constants and time scales of the main decoherence channels indicate, that the proposed setup is a realistic route for large-scale quantum information processing with neutral atoms, exploiting already available state-of-the-art technology.

QEC allows one to mitigate the effects of a noisy environment and faulty operations on the information stored in quantum states. The main idea of QEC is to embed quantum information in a Hilbert space of larger dimension, enabling a distribution of information that leads to resilience against noise. For this purpose, one can use the combined Hilbert space of multiple qubits, while an alternative approach is to use a single system that has a larger Hilbert space. A way to realize the latter goes back to a seminal work of Gottesmann, Kitaev, and Preskill (GKP) [74], which proposed encoding a qubit into a harmonic oscillator.

The finite-dimensional version of the GKP code encodes a qubit into a collective spin. This code is based on a set of commuting operators, the stabilizer set, which can be measured simultaneously. The quantum information is then encoded into the joint (+1)-eigenspace of the stabilizer set. Errors acting on the encoded information may lead to a change in one or few stabilizer measurements, which helps in detecting and correcting errors.

The formalism of the finite-dimensional GKP code rests on the generalized Pauli operators X and Z defined by

$$\begin{equation}X\vert {j\rangle }_{d}=\vert (j+1)\;\text{mod}\enspace {d\rangle }_{d},\end{equation} \tag{ 9a }$$

$$\begin{equation}Z\vert {j\rangle }_{d}={\omega }^{j}\vert {j\rangle }_{d}.\end{equation} \tag{ 9b }$$

These operators obey the relation ZX = ωXZ, with ω = exp(2πi/d). In particular, two operators of the form X^α Z^β and X^γ Z^δ can commute with each other, if their exponents α, β, γ, and δ are chosen appropriately. Such operators can then be used to form the stabilizer group of a code.

To illustrate the idea behind the finite dimensional GKP code, we choose d = 8 corresponding to seven atoms in an optical tweezer. In order to generate the stabilizer group, we choose the operators

$$\begin{equation}{S}_{1}={X}^{4},\end{equation} \tag{ 10a }$$

$$\begin{equation}{S}_{2}={Z}^{4}.\end{equation} \tag{ 10b }$$

Notice that S₁ and S₂ commute. Their joint (+1)-eigenspace defines a two-dimensional subspace, which can be used to encode a logical qubit. A basis for this code space is

$$\begin{equation}\vert \bar{0}\rangle =\frac{1}{\sqrt{2}}(\vert {0\rangle }_{8}+\vert {4\rangle }_{8}),\end{equation} \tag{ 11a }$$

$$\begin{equation}\vert \bar{1}\rangle =\frac{1}{\sqrt{2}}(\vert {2\rangle }_{8}+\vert {6\rangle }_{8}).\end{equation} \tag{ 11b }$$

A state is then encoded as $\vert {\psi }_{L}\rangle =\alpha \vert \bar{0}\rangle +\beta \vert \bar{1}\rangle$ and one can detect all errors X^a Z^b that have |a|, |b| ⩽ 1. However, from the fact that $X\vert \bar{0}\rangle ={X}^{-1}\vert \bar{1}\rangle$ we infer that some of these errors can only be detected, but not corrected, because the action of X and X⁻¹ on the logical states cannot be distinguished. More explicitly, this can be seen from the conditions for QEC [75].

Experimentally, the two logical states $\vert \bar{0}\rangle ,\vert \bar{1}\rangle$ are prepared by using a control qubit as shown in figure 4(b): the collective spin and the qubit are initially prepared in the product state |0⟩₂|0⟩₈. Applying a Hadamard gate onto the control qubit leads to $\frac{1}{\sqrt{2}}(\vert {0\rangle }_{2}+\vert {1\rangle }_{2})\vert {0\rangle }_{8}$. Next, we employ a conditional CX⁴ gate defined as

$$\begin{equation}{CX}^{4}\vert {j\rangle }_{2}\vert {k\rangle }_{8}=(\mathbb{1}\otimes {X}^{4j})\vert {j\rangle }_{2}\vert {k\rangle }_{8},\end{equation} \tag{ 12 }$$

which results in $\frac{1}{\sqrt{2}}(\vert {0\rangle }_{2}\vert {0\rangle }_{8}+\vert {1\rangle }_{2}\vert {4\rangle }_{8})$. A second Hadamard gate on the control qubit yields

$$\begin{equation}\frac{1}{\sqrt{2}}\left[\vert {0\rangle }_{2}\vert \bar{0}\rangle +\vert {1\rangle }_{2}(Z\vert \bar{0}\rangle )\right].\end{equation} \tag{ 13 }$$

After measuring the Z eigenvalue of the control qubit, the collective spins is projected onto either $\vert \bar{0}\rangle$ or $Z\vert \bar{0}\rangle$. In the latter case, an additional Z gate is applied. In this way, the system is deterministically prepared in $\vert \bar{0}\rangle$. In an experiment the circuit depicted in figure 4(b) can be realized as follows: the circuit contains Hadamard gates, which acts on a qubit, a Z-gate, which acts on a qudit, a CX⁴ gate acting on a qudit and qubit, and a local readout. The Hadamard gate on the qubit is realized by performing a π/2 pulse locally. The Z-gate on the qudit is performed by changing the detuning of the microwave. The CX⁴ is not a native gate described in section 2. However, since the gates generated by L_z , L_x , ${L}_{z}^{2}$, and the CZ gate are universal, we can decompose the CX⁴ gate into a judiciously chosen sequence of local gates and two-body gates. The local gates correspond to local optical control, microwave pulses, and local collision of atoms. The two-body CZ-gate corresponds to phonon entanglement, which can be switched off via optical shelving.

The encoding (11) allows one to detect certain errors affecting the logical qubit. Consider errors of the form X^a Z^b acting on |ψ_L ⟩. From the commutation relation S₁ X^a Z^b = ω^−4b X^a Z^b S₁ and S₂ X^a Z^b = ω^4a X^a Z^b S₂, one obtains

$$\begin{equation}{S}_{1}\left({X}^{a}{Z}^{b}\vert {\psi }_{L}\rangle \right)={\omega }^{-4b}\left({X}^{a}{Z}^{b}\vert {\psi }_{L}\rangle \right),\end{equation} \tag{ 14a }$$

$$\begin{equation}{S}_{2}\left({X}^{a}{Z}^{b}\vert {\psi }_{L}\rangle \right)={\omega }^{4a}\left({X}^{a}{Z}^{b}\vert {\psi }_{L}\rangle \right).\end{equation} \tag{ 14b }$$

Thus an error may move the logical qubit out of the (+1)-eigenspace of S₁ and S₂. In particular, this happens for all non-trivial X^a Z^b with |a|, |b| ⩽ 1. The circuit in figure 4(c) can then be used to detect when this happens. Likewise, it can be shown that the same circuit also detects any linear combinations of such errors.

Zoom In Zoom Out Reset image size

To also correct errors, one has to use a larger collective spin. Using d = 18, the stabilizers S₁ = X⁶ and S₂ = Z⁶ yield a code with basis [74]

$$\begin{equation}\vert \bar{0}\rangle =\frac{1}{\sqrt{3}}\left(\vert {0\rangle }_{18}+\vert {6\rangle }_{18}+\vert 1{2\rangle }_{18}\right),\end{equation} \tag{ 15a }$$

$$\begin{equation}\vert \bar{1}\rangle =\frac{1}{\sqrt{3}}\left(\vert {3\rangle }_{18}+\vert {9\rangle }_{18}+\vert 1{5\rangle }_{18}\right).\end{equation} \tag{ 15b }$$

The circuit given in figure 4(c) can now not only detect, but also correct all errors X^a Z^b with |a|, |b| ⩽ 1. Then, the encoded states can also be manipulated through a set of universal gates. To be specific, we return to the case of d = 8. The logical gates $\bar{X}={X}^{2}$ and $\bar{Z}={Z}^{2}$ act like the usual spin 1/2 Pauli matrices on $\vert \bar{0}\rangle$ and $\vert \bar{1}\rangle$, i.e.

$$\begin{equation}\bar{X}\vert \bar{0}\rangle =\vert \bar{1}\rangle ,\qquad \bar{Z}\vert \bar{0}\rangle =\vert \bar{0}\rangle ,\end{equation} \tag{ 16a }$$

$$\begin{equation}\bar{X}\vert \bar{1}\rangle =\vert \bar{0}\rangle ,\qquad \bar{Z}\vert \bar{1}\rangle =-\vert \bar{1}\rangle ,\end{equation} \tag{ 16b }$$

such that ${\bar{X}}^{2}={\bar{Z}}^{2}=\mathbf{1}$ and $\bar{X}\bar{Z}+\bar{Z}\bar{X}=0$.

As explained in section 1, we can implement all single-qubit logical gates through Trotterization, which allows to synthesize

$$\begin{equation}\bar{\mathcal{R}}(\alpha ,\beta ,\gamma ,\delta )={\text{e}}^{\text{i}\alpha }\enspace {\text{e}}^{\text{i}\beta \bar{Z}}\enspace {\text{e}}^{\text{i}\gamma \bar{X}}\enspace {\text{e}}^{\text{i}\delta \bar{Z}}.\end{equation} \tag{ 17 }$$

Moreover, two logical qubits can be entangled by using a logical controlled Z gate,

$$\begin{equation}C\bar{Z}(\mathbf{x},\mathbf{y})={\text{e}}^{\text{i}\pi \bar{Z}(\mathbf{x})\bar{Z}(\mathbf{y})},\end{equation} \tag{ 18 }$$

which can be also synthesized via Trotterzation using the gates generated by L_x , L_z and ${L}_{z}^{2}$, and the entangling gate CZ(x, y) given in equation (8). Together, the operations $\bar{\mathcal{R}}(\alpha ,\beta ,\gamma ,\delta )$ and $C\bar{Z}$ then generate a universal set of logical gates [60].

In summary, we have constructed an explicit example of a quantum error-detecting code for d = 8, described the experimental preparation of its logical states, and provided a universal gate set to manipulate the stored quantum information. The presented formalism will also work in larger dimensions, allowing for the detection and correction of a larger set of errors [74]. The here proposed platform with a suitable quantum error-correcting scheme, such as finite GKP codes, may allow one to go beyond the abilities of noisy intermediate-scale quantum (NISQ) technologies.

In this work, we have proposed a mixture of two ultracold atomic species as a platform for universal quantum computation. Our proposed implementation uses long-range entangling gates and, in addition, allows for the realization of QEC. The idea of using phonon mediated interactions for quantum information processing is also discussed in the solid state context and made substantial progress [76, 77]. The presented spin–phonon system is not only useful for the processing of quantum information but is also interesting from a quantum many-body perspective. The many body systems realized by our platform is similar to gauge theories coupled to a Higgs field [78–80] and hence is a promising candidate for the investigation of topological matter. Further, the spin–phonon system proposed here may give access to the physics of the Peierls transition [81], the Jahn–Teller effect [82] in a many-body context as well as the study of frustrated spin models [83]. The versatility of this platform also makes it an ideal candidate for a fully programmable quantum simulation, whose large potential has been demonstrated in Rydberg and trapped ion systems [4, 5, 84, 85]. For example, atomic mixtures have been proposed for the quantum simulation of quantum chemistry problems [86].

In future work, the building blocks of our proposed platform could be exchanged ad libitum exploiting the versatility of atomic, molecular, and atomic physics. The basic unit of information can be stored, e.g. with spinful fermions in an optical tweezer [71] or higher dimensional spins [38, 87, 88]. The entanglement bus may be substituted by more complex many-body excitations like magnons or rotons [89, 90]. This exchange of the basic unit of information and the entanglement mechanism has the potential to lead to new implementations of quantum algorithms and QEC schemes, see for example [91–94]. Thus, ultracold atom mixtures present a promising platform to implement fault-tolerant quantum computation in a scalable platform, paving a way beyond the abilities of NISQ technology.

The authors are grateful for fruitful discussions with J Eisert, M Gaerttner, T Gasenzer, M Gluza, S Jochim, M Oberthaler, A P Orioli, P Preiss, H Strobel and all the members of the SynQs seminar. ICFO group acknowledges support from ERC AdG NOQIA, State Research Agency AEI (Severo Ochoa Center of Excellence CEX2019-000910-S, Plan National FIDEUA PID2019-106901GB-I00/10.13039/501100011033, FPI, QUANTERA MAQS PCI2019-111828-2/ 10.13039/501100011033), Fundació Privada Cellex, Fundació Mir-Puig, Generalitat de Catalunya (AGAUR Grant No. 2017 SGR 1341, CERCA program, QuantumCAT / U16-011424, co-funded by ERDF Operational Program of Catalonia 2014-2020), EU Horizon 2020 FET-OPEN OPTOLogic (Grant No 899794), and the National Science Centre, Poland (Symfonia Grant No. 2016/20/W/ST4/00314), Marie Sklodowska-Curie grant STREDCH No.101029393, La Caixa Junior Leaders fellowships (ID100010434), and EU Horizon 2020 under Marie Sklodowska-Curie Grant Agreement No. 847648 (LCF/BQ/PI19/11690013, LCF/BQ/PI20/11760031, LCF/BQ/PR20/11770012). FH acknowledges the Government of Spain (FIS2020-TRANQI and Severo Ochoa CEX2019-000910-S), Fundació Cellex, Fundació Mir-Puig, Generalitat de Catalunya (CERCA, AGAUR SGR 1381), and the Foundation for Polish Science through TEAM-NET (POIR.04.04.00-00-17C1/18-00). PH acknowledges support by Q@TN—Quantum Science and Technologies at Trento, the Provincia Autonoma di Trento, and the ERC Starting Grant StrEnQTh (Project-ID 804305). This work is part of and supported by the DFG Collaborative Research Center 'SFB 1225 (ISOQUANT)'. FJ acknowledges the DFG support through the Project FOR 2724, the Emmy-Noether Grant (Project-ID 377616843) and support by the Bundesministerium für Wirtschaft und Energie through the project 'EnerQuant' (Project-ID 03EI1025C).

The data that support the findings of this study are available upon reasonable request from the authors.

The type A atoms are placed in an external magnetic field, which allows to prepare the atoms in magnetic substates m ∈ {0, 1}. Further, the atoms of type A are localized at specific sites via optical potentials (see figure 1). The single particle physics is described by the Hamiltonian

$$\begin{equation}{H}_{\mathrm{A}}^{m}(\mathbf{x})=-\frac{{\hslash }^{2}{\nabla }_{\mathbf{x}}^{2}}{2{M}_{\mathrm{A}}}+{V}_{\mathrm{A}}(\mathbf{x})+{E}_{\mathrm{A}}^{m}(B),\end{equation} \tag{ A1 }$$

where M_A is the mass of the atoms and ${V}_{\mathrm{A}}\left(\mathbf{x}\right)$ is the optical potential confining the atoms. The energy shift ${E}_{m}^{\mathrm{A}}(B)$ determined by the magnetic field B is given by

$$\begin{equation}{E}_{m}^{\mathrm{A}}(B)={p}_{m}(B)m+{q}_{m}(B){m}^{2},\end{equation} \tag{ A2 }$$

where p_m (B) parametrizes the linear Zeeman shift and q_m (B) the quadratic Zeeman shift. The field operators A_m (x) and ${A}_{m}^{{\dagger}}(\mathbf{x})$ annihilate and create a particle at position $\mathbf{x}={({x}_{1},\dots ,{x}_{d})}^{\text{T}}$ respectively and fulfill canonical commutation relations. The many-body Hamiltonian is given by

$$\begin{align}{H}_{\mathrm{A}}& =\sum\limits _{m}{\int }_{\mathbf{x}}{A}_{m}^{{\dagger}}(\mathbf{x}){H}_{\mathrm{A}}^{m}(\mathbf{x}){A}_{m}(\mathbf{x})+\sum\limits _{m,n}{\int }_{\mathbf{x}}\frac{{g}_{\mathrm{A}}^{mn}}{2}{A}_{m}^{{\dagger}}(\mathbf{x}){A}_{n}^{{\dagger}}(\mathbf{x}){A}_{n}(\mathbf{x}){A}_{m}(\mathbf{x}).\end{align} \tag{ A3 }$$

The atoms are localized at different sites y, which allows one to expand the field operators as

$$\begin{equation}{A}_{m}(\mathbf{x})=\sum\limits _{\mathbf{y}}\;{\varphi }_{\mathrm{A}}(\mathbf{x}-\mathbf{y}){a}_{m}(\mathbf{y})\end{equation} \tag{ A4 }$$

with a wave function φ_A localized at 0 and which does not depend on the magnetic substate. Inserting the expansion (A4) into H_A and neglecting the tunneling terms, we obtain the Hamiltonian

$$\begin{equation}{H}_{\mathrm{A}}=\frac{1}{2}\sum\limits _{m,n,\mathbf{y}}\;{\tilde{g}}_{\mathrm{A}}^{mn}{a}_{m}^{{\dagger}}(\mathbf{y}){a}_{n}^{{\dagger}}(\mathbf{y}){a}_{m}(\mathbf{y}){a}_{n}(\mathbf{y})\enspace \end{equation} \tag{ A5 }$$

with the overlap integrals

$$\begin{equation}{\tilde{g}}_{\mathrm{A}}^{mn}={g}_{\mathrm{A}}^{mn}{\int }_{\mathbf{x}}\vert {\varphi }_{\mathrm{A}}(\mathbf{x}-\mathbf{y}){\vert }^{4}.\end{equation} \tag{ A6 }$$

To be explicit, we approximate φ_A(x) by the ground state wavefunction of an isotropic harmonic oscillator φ_A(x) = φ_A(x₁)φ_A(x₂)φ_A(x₃) with

$$\begin{equation}{\varphi }_{\mathrm{A}}({x}_{i})={(\sqrt{\pi }{\sigma }_{\mathrm{A}})}^{-1/2}\enspace {\mathrm{e}}^{-\frac{1}{2}{\left(\frac{{x}_{i}}{{\sigma }_{A}}\right)}^{2}}\end{equation} \tag{ A7 }$$

with the characteristic length σ_A. Calculating the overlap integral leads to the dimensional reduced coupling constant

$$\begin{equation}{\tilde{g}}_{\mathrm{A}}^{mn}={(2\pi )}^{-3/2}({g}_{\mathrm{A}}^{mn}/{\sigma }_{\mathrm{A}}^{3}).\end{equation} \tag{ A8 }$$

Inserting the Schwinger representation of the angular momentum into equation (A5) leads to equation (2) with the coupling constants

$$\begin{equation}\chi =\frac{1}{2}\left({\tilde{g}}_{\mathrm{A}}^{00}+{\tilde{g}}_{\mathrm{A}}^{11}-2{\tilde{g}}_{\mathrm{A}}^{10}\right),\end{equation} \tag{ A9a }$$

$$\begin{equation}{\Delta}=\frac{1}{2}({N}_{\mathrm{A}}-1)({\tilde{g}}_{\mathrm{A}}^{11}-{\tilde{g}}_{\mathrm{A}}^{00})+{E}_{\mathrm{A}}^{1}-{E}_{\mathrm{A}}^{0},\end{equation} \tag{ A9b }$$

where N_A is the particle number on each site.

In this appendix, we demonstrate for ℓ = 1 that L_z , L_x , and ${L}_{z}^{2}$ are sufficient to generate all U(3) matrices. We use the spin matrices

$$\begin{equation*}{L}_{x}=\frac{1}{\sqrt{2}}\left(\begin{matrix}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \end{matrix}\right),\qquad {L}_{y}=\frac{1}{\sqrt{2}i}\left(\begin{matrix}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \end{matrix}\right),\qquad {L}_{z}=\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \end{matrix}\right).\end{equation*}$$

Consider the following matrices generated by commutators of L_x , L_z , and ${L}_{z}^{2}$

$$\begin{align}{\mathcal{M}}_{1}& ={L}_{x},& {\mathcal{M}}_{2}& ={L}_{z},& {\mathcal{M}}_{3}& ={L}_{z}^{2},\\ {\mathcal{M}}_{4}& =i[{\mathcal{M}}_{1},{\mathcal{M}}_{2}],& {\mathcal{M}}_{5}& =i[{\mathcal{M}}_{3},{\mathcal{M}}_{1}],& {\mathcal{M}}_{6}& =i[{\mathcal{M}}_{3},{\mathcal{M}}_{4}],\\ {\mathcal{M}}_{7}& =i[{\mathcal{M}}_{5},{\mathcal{M}}_{1}],& {\mathcal{M}}_{8}& =i[{\mathcal{M}}_{5},{\mathcal{M}}_{4}],& {\mathcal{M}}_{9}& =i[{\mathcal{M}}_{6},{\mathcal{M}}_{4}].\end{align} \tag{ B1 }$$

These commutators form a basis for the Lie algebra of U(3). This can be explicitly checked by constructing a change of basis from the ${\left\{{\mathcal{M}}_{i}\right\}}_{i=1}^{9}$ to the canonical basis of Hermitian matrices ${\left\{{M}_{i}\right\}}_{i=1}^{9}$ given by

$$\begin{align}{M}_{1}& =\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),& {M}_{2}& =\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),& {M}_{3}& =\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right),& \\ {M}_{4}& =\left(\begin{matrix}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),& {M}_{5}& =\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),& {M}_{6}& =\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \end{matrix}\right),& \\ {M}_{7}& =\left(\begin{matrix}\hfill 0\hfill & \hfill i\hfill & \hfill 0\hfill \\ \hfill -i\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),& {M}_{8}& =\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill i\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -i\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),& {M}_{9}& =\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill i\hfill \\ \hfill 0\hfill & \hfill -i\hfill & \hfill 0\hfill \end{matrix}\right).& \end{align} \tag{ B2 }$$

The Hamiltonian involving only the atomic species B is given by

$$\begin{equation}{H}_{\mathrm{B}}={\int }_{\mathbf{x}}{B}^{{\dagger}}(\mathbf{x}){H}_{\mathrm{B}}^{0}(\mathbf{x})B(\mathbf{x})+\frac{{g}_{\mathrm{B}}}{2}{\int }_{\mathbf{x}}{B}^{{\dagger}}(\mathbf{x}){B}^{{\dagger}}(\mathbf{x})B(\mathbf{x})B(\mathbf{x})\end{equation} \tag{ C1 }$$

with

$$\begin{equation}{H}_{\mathrm{B}}^{0}(\mathbf{x})=-\frac{{\hslash }^{2}{\nabla }_{\mathbf{x}}^{2}}{2{M}_{\mathrm{B}}}+{V}_{\mathrm{B}}(\mathbf{x}).\end{equation} \tag{ C2 }$$

In order to create a one- or two-dimensional phononic bath we apply the following harmonic and isotropic confinement

$$\begin{equation}{V}_{\mathrm{B}}(\mathbf{x})=\frac{1}{2}{M}_{\mathrm{B}}{\omega }_{\mathrm{B}}^{2}\sum\limits _{i=n+1}^{d}{x}_{i}^{2}.\end{equation} \tag{ C3 }$$

For sufficiently low temperatures, the transversal degrees of freedom will not be excited, which will confine the particles effectively to one (n = 1) or two (n = 2) dimensions respectively. This freezing out of the transversal directions allows one to write the field operator as

$$\begin{equation}B(\mathbf{x})=B({x}_{1},\dots ,{x}_{n}){\varphi }_{\mathrm{B}}({x}_{n+1},\dots ,{x}_{d})\end{equation} \tag{ C4 }$$

with the transversal wave function φ_B(x_n+1, ..., x_d ) and B(x₁, ..., x_n ) the annihilation operator in n dimensions, and in the following we will use $\tilde{\mathbf{x}}={({x}_{1},\dots ,{x}_{n})}^{\text{T}}$. The stationary Gross–Pitaevskii equation is given by

$$\begin{equation}0=\left[-\frac{{\hslash }^{2}{\nabla }_{\tilde{\mathbf{x}}}^{2}}{2{M}_{\mathrm{B}}}-{\mu }_{\mathrm{B}}+{\tilde{g}}_{\mathrm{B}}\vert {\phi }_{\mathrm{B}}(\tilde{\mathbf{x}}){\vert }^{2}\right]{\phi }_{\mathrm{B}}(\tilde{\mathbf{x}})\end{equation} \tag{ C5 }$$

with the condensate ${\phi }_{\mathrm{B}}(\tilde{\mathbf{x}})$ fulfilling the Dirichlet boundary conditions ${\phi }_{\mathrm{B}}(\tilde{\mathbf{x}})=0$ for $\tilde{\mathbf{x}}\in \partial D$ with D = [0, L]ⁿ and chemical potential μ_B. The coupling constant is given by

$$\begin{equation}{\tilde{g}}_{\mathrm{B}}={g}_{\mathrm{B}}{\int }_{{x}_{n+1},\dots ,{x}_{d}}\vert {\varphi }_{\mathrm{B}}({x}_{n+1},\dots ,{x}_{d}){\vert }^{4},\end{equation} \tag{ C6 }$$

which becomes ${\tilde{g}}_{\mathrm{B}}={(\sqrt{2\pi }{\sigma }_{\mathrm{B}})}^{-(d-n)}{g}_{\mathrm{B}}$ for harmonic confinement. The bulk solution of the Gross–Pitaevskii equation can be approximated by the homogeneous function

$$\begin{equation}{\phi }_{\mathrm{B}}(\tilde{\mathbf{x}})\approx \sqrt{\frac{{\mu }_{\mathrm{B}}}{{\tilde{g}}_{\mathrm{B}}}}\end{equation} \tag{ C7 }$$

leading to the density ${n}_{\mathrm{B}}=\frac{{\mu }_{\mathrm{B}}}{{\tilde{g}}_{\mathrm{B}}}$. In order to study the excitations of the bulk solution, we perform the Bogoliubov approximation

$$\begin{equation}B(\tilde{\mathbf{x}})={\phi }_{\mathrm{B}}(\tilde{\mathbf{x}})+\delta B(\tilde{\mathbf{x}}),\end{equation} \tag{ C8 }$$

where δB and δB^† fulfill canonical commutation relations. The Bogoliubov Hamiltonian approximation of H_B is given by

$$\begin{equation}{H}_{\mathrm{B}}={\int }_{D}\left[\frac{{\hslash }^{2}}{2{M}_{\text{B}}}\vert {\nabla }_{\tilde{\mathbf{x}}}\delta B{\vert }^{2}-{\mu }_{\mathrm{B}}\vert \delta B{\vert }^{2}+2{\tilde{g}}_{\mathrm{B}}\vert \delta B{\vert }^{2}\vert {\phi }_{\mathrm{B}}{\vert }^{2}\right.\left.+\frac{{\tilde{g}}_{\mathrm{B}}}{2}\left[{(\delta {B}^{{\dagger}})}^{2}{({\phi }_{\mathrm{B}})}^{2}+{(\delta B)}^{2}{({\phi }_{\mathrm{B}}^{\ast })}^{2}\right]\right].\end{equation} \tag{ C9 }$$

The equations of motions can be solved by using the mode expansion

$$\begin{equation}\delta B(\tilde{\mathbf{x}},t)=\sum\limits _{\mathbf{k}}\left[{b}_{\mathbf{k}}{u}_{\mathbf{k}}(\tilde{\mathbf{x}}){\text{e}}^{-\text{i}{\omega }_{\mathbf{k}}t}+{b}_{\mathbf{k}}^{{\dagger}}{v}_{\mathbf{k}}^{\ast }(\tilde{\mathbf{x}}){\text{e}}^{\text{i}{\omega }_{\mathbf{k}}t}\right],\end{equation} \tag{ C10 }$$

where the sum does not include the condensate mode [95]. This expansion leads to the following generalized eigenvalue problem

$$\begin{equation}{\omega }_{\mathbf{k}}\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {u}_{\mathbf{k}}\hfill \\ \hfill {v}_{\mathbf{k}}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill h(\tilde{\mathbf{x}})\hfill & \hfill {\tilde{g}}_{\mathrm{B}}{({\phi }_{\mathrm{B}})}^{2}\hfill \\ \hfill {\tilde{g}}_{\mathrm{B}}{({\phi }_{\mathrm{B}}^{\ast })}^{2}\hfill & \hfill h(\tilde{\mathbf{x}})\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {u}_{\mathbf{k}}\hfill \\ \hfill {v}_{\mathbf{k}}\hfill \end{matrix}\right),\end{equation} \tag{ C11 }$$

with

$$\begin{equation}h(\tilde{\mathbf{x}})=-\frac{{\hslash }^{2}{\nabla }_{\tilde{\mathbf{x}}}^{2}}{2{M}_{\mathrm{B}}}-{\mu }_{\mathrm{B}}+2{\tilde{g}}_{\mathrm{B}}\vert {\phi }_{\mathrm{B}}(\tilde{\mathbf{x}}){\vert }^{2}.\end{equation} \tag{ C12 }$$

Solving this generalized eigenvalue problem will lead to the orthonormalization condition

$$\begin{equation}{\int }_{D}\left[{u}_{\mathbf{k}}^{\ast }(\tilde{\mathbf{x}}){u}_{{\mathbf{k}}^{\prime }}(\tilde{\mathbf{x}})-{v}_{\mathbf{k}}^{\ast }(\tilde{\mathbf{x}}){v}_{{\mathbf{k}}^{\prime }}(\tilde{\mathbf{x}})\right]={\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }}.\end{equation} \tag{ C13 }$$

Since the background-field is approximately constant and because of the Dirichlet boundary conditions, we make the ansatz

$$\begin{equation}{u}_{\mathbf{k}}(\tilde{\mathbf{x}})={u}_{\mathbf{k}}{\left(\frac{2}{L}\right)}^{n/2}\prod\limits _{i=1}^{n}\mathrm{sin}({k}_{i}{x}_{i}),\end{equation} \tag{ C14a }$$

$$\begin{equation}{v}_{\mathbf{k}}(\tilde{\mathbf{x}})={v}_{\mathbf{k}}{\left(\frac{2}{L}\right)}^{n/2}\prod\limits _{i=1}^{n}\mathrm{sin}({k}_{i}{x}_{i}),\end{equation} \tag{ C14b }$$

with ${k}_{i}=\frac{{n}_{i}\pi }{L}$ and n_i ⩾ 2 and the amplitudes are given by

$$\begin{equation}{u}_{\mathbf{k}}^{2}=\frac{1}{2}\left(\frac{{\varepsilon }_{\mathbf{k}}}{{\omega }_{\mathbf{k}}}+1\right),\end{equation} \tag{ C15a }$$

$$\begin{equation}{v}_{\mathbf{k}}^{2}=\frac{1}{2}\left(\frac{{\varepsilon }_{\mathbf{k}}}{{\omega }_{\mathbf{k}}}-1\right),\end{equation} \tag{ C15b }$$

with ${\varepsilon }_{\mathbf{k}}=({\hslash }^{2}{\mathbf{k}}^{2})/(2{M}_{\mathrm{B}})+{\tilde{g}}_{\mathrm{B}}\vert {\phi }_{\mathrm{B}}{\vert }^{2}$. The Bogoliubov eigenfrequencies of the excitations are given by

$$\begin{equation}\hslash {\omega }_{\mathbf{k}}=\sqrt{\frac{{\hslash }^{2}{\mathbf{k}}^{2}}{2{M}_{\mathrm{B}}}\left(\frac{{\hslash }^{2}{\mathbf{k}}^{2}}{2{M}_{\mathrm{B}}}+2{\mu }_{\mathrm{B}}\right)}.\end{equation} \tag{ C16 }$$

The Dirichlet boundary conditions can be achieved by box potentials, which lead to a homogeneous Bose–Einstein condensate.

In order to derive the interaction between the phonons and the collective spins, we start from Hamiltonian modeling the interaction between the A and B atoms

$$\begin{equation}{H}_{\mathrm{A}\mathrm{B}}=\sum\limits _{m}{\int }_{\mathbf{x}}\frac{{g}_{\mathrm{A}\mathrm{B}}^{m}}{2}{A}_{m}^{{\dagger}}(\mathbf{x}){A}_{m}(\mathbf{x}){B}^{{\dagger}}(\mathbf{x})B(\mathbf{x}).\end{equation} \tag{ D1 }$$

After reducing the dimension and considering the tight confinement of the spins, i.e. equations (C4) and (A4), we obtain the Hamiltonian

$$\begin{equation}{H}_{\mathrm{A}\mathrm{B}}=\sum\limits _{m,\mathbf{y}}\frac{{\tilde{g}}_{\mathrm{A}\mathrm{B}}^{m}}{2}{a}_{m}^{{\dagger}}(\mathbf{y}){a}_{m}(\mathbf{y}){\int }_{D}\vert {\varphi }_{\mathrm{A}}({x}_{1}-{y}_{1})\dots {\varphi }_{\mathrm{A}}({x}_{n}-{y}_{n}){\vert }^{2}{B}^{{\dagger}}(\tilde{\mathbf{x}})B(\tilde{\mathbf{x}}),\end{equation} \tag{ D2 }$$

where we neglected hopping of the A atoms and introduced the coupling constant

$$\begin{equation}{\tilde{g}}_{\mathrm{A}\mathrm{B}}^{m}={g}_{\mathrm{A}\mathrm{B}}^{m}{\int }_{{x}_{n+1},\dots ,{x}_{d}}\vert {\varphi }_{\mathrm{B}}({x}_{n+1})\dots {\varphi }_{\mathrm{B}}({x}_{d}){\vert }^{2}\vert {\varphi }_{\mathrm{A}}({x}_{n+1}-{y}_{n+1})\dots {\varphi }_{\mathrm{A}}({x}_{d}-{y}_{d}){\vert }^{2},\end{equation} \tag{ D3 }$$

which can be written as

$$\begin{equation}{\tilde{g}}_{\mathrm{A}\mathrm{B}}^{m}={g}_{\mathrm{A}\mathrm{B}}^{m}{\left[\pi ({\sigma }_{\mathrm{A}}^{2}+{\sigma }_{\mathrm{B}}^{2})\right]}^{-(d-n)/2}\end{equation} \tag{ D4 }$$

for harmonic confinement of the A and B atoms with harmonic oscillator length scale σ_A and σ_B respectively. Expanding the field operator B in fluctuations as in equation (C8) and neglecting terms of order $\mathcal{O}(\delta {B}^{2})$ one obtains a Hamiltonian

$$\begin{equation}{H}_{\mathrm{A}\mathrm{B}}={H}_{\mathrm{A}\mathrm{B}}^{(0)}+{H}_{\mathrm{A}\mathrm{B}}^{(1)},\end{equation} \tag{ D5 }$$

where ${H}_{\mathrm{A}\mathrm{B}}^{(0)}$ is independent of the fluctuations and ${H}_{\mathrm{A}\mathrm{B}}^{(1)}$ is linear in the fluctuations. The first contribution is given by

$$\begin{equation}{H}_{\mathrm{A}\mathrm{B}}^{(0)}=\sum\limits _{m,\mathbf{y}}\;{{\Delta}}_{m}{a}_{m}^{{\dagger}}(\mathbf{y}){a}_{m}(\mathbf{y}),\end{equation} \tag{ D6 }$$

with the coupling constant

$$\begin{equation}{{\Delta}}_{m}=\frac{1}{2}{\tilde{g}}_{\mathrm{A}\mathrm{B}}^{m}{n}_{\mathrm{B}}{\int }_{D}\vert {\varphi }_{\mathrm{A}}({x}_{1}-{y}_{1})\dots {\varphi }_{\mathrm{A}}({x}_{n}-{y}_{n}){\vert }^{2},\end{equation} \tag{ D7 }$$

and since the harmonic oscillator wave functions are normalized we obtain

$$\begin{equation}{{\Delta}}_{m}=\frac{1}{2}{\tilde{g}}_{\mathrm{A}\mathrm{B}}^{m}{n}_{\mathrm{B}}.\end{equation} \tag{ D8 }$$

The contribution linear in the fluctuations is given by

$$\begin{equation}{H}_{\mathrm{A}\mathrm{B}}=\frac{\sqrt{{n}_{\mathrm{B}}}}{2}\sum\limits _{m,\mathbf{y}}\;{\tilde{g}}_{\mathrm{A}\mathrm{B}}^{m}{a}_{m}^{{\dagger}}(\mathbf{y}){a}_{m}(\mathbf{y}){\int }_{D}[\delta B(\tilde{\mathbf{x}})+\text{H.c.}]\vert {\varphi }_{\mathrm{A}}({x}_{1}-{y}_{1})\dots {\varphi }_{\mathrm{A}}({x}_{n}-{y}_{n}){\vert }^{2}.\end{equation} \tag{ D9 }$$

Inserting the mode expansion (C10) into ${H}_{\mathrm{A}\mathrm{B}}^{(1)}$ we obtain

$$\begin{equation}{H}_{\mathrm{A}\mathrm{B}}^{(1)}=\sum\limits _{m,\mathbf{y},\mathbf{k}}\;{\tilde{g}}_{\mathrm{A}\mathrm{B},\mathbf{k}}^{m}(\mathbf{y}){a}_{m}^{{\dagger}}(\mathbf{y}){a}_{m}(\mathbf{y})\left[{b}_{\mathbf{k}}+\text{H.c.}\right],\end{equation} \tag{ D10 }$$

where we introduced the coupling constant

$$\begin{equation}{\tilde{g}}_{\mathrm{A}\mathrm{B},\mathbf{k}}^{m}(\mathbf{y})=\frac{\sqrt{{n}_{\mathrm{B}}}}{2}{\tilde{g}}_{\mathrm{A}\mathrm{B}}^{m}{\int }_{D}\vert {\varphi }_{\mathrm{A}}({x}_{1}-{y}_{1})\dots {\varphi }_{\mathrm{A}}({x}_{n}-{y}_{n}){\vert }^{2}[{u}_{\mathbf{k}}(\tilde{\mathbf{x}})+{v}_{\mathbf{k}}(\tilde{\mathbf{x}})].\end{equation} \tag{ D11 }$$

Inserting (C14) and (A7) we obtain for L ≫ σ_B approximately

$$\begin{align}{\tilde{g}}_{\mathrm{A}\mathrm{B},\mathbf{k}}^{m}(\mathbf{y})& ={\tilde{g}}_{\mathrm{A}\mathrm{B}}^{m}\sqrt{{n}_{\mathrm{B}}}{\left(\frac{2}{L}\right)}^{n/2}({u}_{\mathbf{k}}+{v}_{\mathbf{k}}){\mathrm{e}}^{-\frac{1}{4}{(\mathbf{k}{\sigma }_{ A})}^{2}}\\ & \times \prod\limits _{i=1}^{n}\mathrm{sin}({k}_{i}{y}_{i}).\end{align} \tag{ D12 }$$

Using the Schwinger representation (see equation (1)), we obtain the interaction between the spins and the phonons

$$\begin{equation}{H}_{\mathrm{A}\mathrm{B}}^{(1)}=\sum\limits _{m,\mathbf{y},\mathbf{k}}\;{\tilde{g}}_{\mathrm{A}\mathrm{B},\mathbf{k}}^{m}(\mathbf{y})\left[L(\mathbf{y})-{(-1)}^{m}{L}_{z}(\mathbf{y})\right]\left({b}_{\mathbf{k}}+\text{H.c.}\right),\end{equation} \tag{ D13 }$$

with L(y) being the length of the angular momentum on site y. Reshuffling terms leads to equation (4) with the coupling constants

$$\begin{equation}{\bar{g}}_{\mathbf{k}}(\mathbf{y})=L(\mathbf{y})[{\tilde{g}}_{\mathrm{A}\mathrm{B},\mathbf{k}}^{0}(\mathbf{y})+{\tilde{g}}_{\mathrm{A}\mathrm{B},\mathbf{k}}^{1}(\mathbf{y})],\end{equation} \tag{ D14a }$$

$$\begin{equation}\delta {g}_{\mathbf{k}}(\mathbf{y})={\tilde{g}}_{\mathrm{A}\mathrm{B},\mathbf{k}}^{1}(\mathbf{y})-{\tilde{g}}_{\mathrm{A}\mathrm{B},\mathbf{k}}^{0}(\mathbf{y}).\end{equation} \tag{ D14b }$$

Assuming Ω(y) = 0 and given the approximations of appendices A, C and D the Hamiltonian H = H_A + H_B + H_AB is diagonal in L_z , which allows us to treat phonons and spins separately. The Heisenberg equation of motion for the phonons is

$$\begin{equation}i\hslash {\partial }_{t}{b}_{\mathbf{k}}=\hslash {\omega }_{\mathbf{k}}{b}_{\mathbf{k}}+\delta {b}_{\mathbf{k}}\end{equation} \tag{ E1 }$$

of b_k operators, where we introduced the abbreviation

$$\begin{equation}\delta {b}_{\mathbf{k}}=\sum\limits _{\mathbf{y}}\left[{\bar{g}}_{\mathbf{k}}(\mathbf{y})L(\mathbf{y})+\delta {g}_{\mathbf{k}}(\mathbf{y}){L}_{z}(\mathbf{y})\right].\end{equation} \tag{ E2 }$$

We define a shifted annihilation operator as

$$\begin{equation}{\beta }_{\mathbf{k}}={b}_{\mathbf{k}}+{(\hslash {\omega }_{\mathbf{k}})}^{-1}\delta {b}_{\mathbf{k}}.\end{equation} \tag{ E3 }$$

Inserting the shifted operator in the Hamiltonian equation (4) leads to a spin–spin interaction

$$\begin{equation}{H}_{I}=-\sum\limits _{\mathbf{x},\mathbf{y}}\;g(\mathbf{x},\mathbf{y}){L}_{z}(\mathbf{x}){L}_{z}(\mathbf{y})\end{equation} \tag{ E4 }$$

with

$$\begin{equation}g(\mathbf{x},\mathbf{y})=\sum\limits _{\mathbf{k}}{(\hslash {\omega }_{\mathbf{k}})}^{-1}\delta {g}_{\mathbf{k}}(\mathbf{x})\delta {g}_{\mathbf{k}}(\mathbf{y}).\end{equation} \tag{ E5 }$$

Inserting the explicit expression for δg_k (x), we obtain equation (6).