Quantum Simulation of Abelian Lattice Gauge Theories via State-Dependent Hopping

TL;DR

This paper proposes a quantum simulation method for Abelian gauge theories using state-dependent hopping controlled by impurity particles, applicable to systems like quantum electrodynamics, and assesses its feasibility in various experimental setups.

Contribution

It introduces a novel quantum simulator architecture for U(1) gauge theories based on controllable particle hopping via impurity states, avoiding fine-tuning of interactions.

Findings

Feasible experimental scheme for simulating gauge theories.

Effective in double well potentials at single-particle and many-body levels.

Applicable to atom-ion systems with micro-motion considerations.

Abstract

We develop a quantum simulator architecture that is suitable for the simulation of $U(1)$ Abelian gauge theories such as quantum electrodynamics. Our approach relies on the ability to control the hopping of a particle through a barrier by means of the internal quantum states of a neutral or charged impurity-particle sitting at the barrier. This scheme is experimentally feasible, as the correlated hopping does not require fine-tuning of the intra- and inter-species interactions. We investigate the applicability of the scheme in a double well potential, which is the basic building block of the simulator, both at the single-particle and the many-body mean-field level. Moreover, we evaluate its performance for different particle interactions and trapping, and, specifically for atom-ion systems, in the presence of micro-motion.

Equations

$$\begin{array}{rl}{\hat{H}}_{\text{QLM}}=& -J\sum _k\left({\hat{b}}_k^{†}{\hat{U}}_{k,k+1}{\hat{b}}_{k+1}+\text{h.c.}\right)\\ & +m\sum _k{\left(-1\right)}^k{\hat{b}}_k^{†}{\hat{b}}_k+\frac{g^2}{2}\sum _k{\hat{E}}_{k,k+1}^2.\end{array}$$

$${\hat{G}}_k=\sum _{\hat{\mu }=1}^D\left({\hat{E}}_{k,k+\hat{\mu }}-{\hat{E}}_{k-\hat{\mu },k}\right)-{\hat{\rho }}_k$$

$$\begin{array}{rl}{e}^{i{\sum }_{k^{\prime }}{\varphi }_{k^{\prime }}{\hat{G}}_{k^{\prime }}}{\hat{b}}_k{e}^{-i{\sum }_{k^{\prime \prime }}{\varphi }_{k^{\prime \prime }}{\hat{G}}_{k^{\prime \prime }}}& =e^{i{\varphi }_k}{\hat{b}}_k\\ e^{i{\sum }_{k^{\prime }}{\varphi }_{k^{\prime }}{\hat{G}}_{k^{\prime }}}{\hat{U}}_{k,k+1}{e}^{-i{\sum }_{k^{\prime \prime }}{\varphi }_{k^{\prime \prime }}{\hat{G}}_{k^{\prime \prime }}}& =e^{i{\varphi }_k}{\hat{U}}_{k,k+1}{e}^{-i{\varphi }_{k+1}}\\ e^{i{\sum }_{k^{\prime }}{\varphi }_{k^{\prime }}{\hat{G}}_{k^{\prime }}}{\hat{E}}_{k,k+1}{e}^{-i{\sum }_{k^{\prime \prime }}{\varphi }_{k^{\prime \prime }}{\hat{G}}_{k^{\prime \prime }}}& =0.\end{array}$$

$$\hat{H}=-J\frac{N}{2}\sum _k\left(e^{-i{\theta }_k}{\tilde{\sigma }}_{k,k+1}^{+}{e}^{i{\theta }_{k+1}}+\text{H.c.}\right).$$

$$\hat{H}=\sum _k{\hat{b}}_k^{†}{\hat{J}}_{k,k+1}{\hat{b}}_{k+1}+\text{H.c.}+\sum _k\frac{{\hat{U}}_k}{2}{\hat{n}}_k({\hat{n}}_k-1).$$

$${\hat{H}}_{(k,k+1)}^0=E_k^0{\hat{n}}_k+E_{k+1}^0{\hat{n}}_{k+1}+\frac{\hslash {\Omega }_R}{2}{\hat{\sigma }}_{k,k+1}^x.$$

$$\hat{H}=\sum _k{\hat{H}}_{(k,k+1)}+\text{on-site terms}.$$

$${\hat{H}}_{(k,k+1)}^{\text{hopp}}={\hat{b}}_k^{†}(J_{k,k+1}^0+J_{k,k+1}^z{\hat{\sigma }}_{k,k+1}^z){\hat{b}}_{k+1}+\text{H.c.},$$

$$\begin{array}{rl}{\hat{H}}_{\text{rot}}=& e^{i{\hat{H}}_{(k,k+1)}^{0}t}{\hat{H}}_{(k,k+1)}^{\text{hopp}}{e}^{-i{\hat{H}}_{(k,k+1)}^{0}t}\\ =& e^{i\left(E_k^0-E_{k+1}^0\right)t}{\hat{b}}_k^{†}{\hat{b}}_{k+1}\\ & \left(\frac{J_{k,k+1}^0}{2}+J_{k,k+1}^z{e}^{-i\hslash {\Omega }_{R}t}{\tilde{\sigma }}_{k,k+1}^{+}+\text{H.c.}\right)+\text{H.c.}\end{array}$$

$${\hat{H}}_{(k,k+1)}^{\text{t}}=J_{k,k+1}^z{\hat{b}}_k^{†}{\tilde{\sigma }}_{k,k+1}^{+}{\hat{b}}_{k+1}+\text{H.c.}$$

$${\hat{H}}_{\text{i}}^{\text{int}}$$

$${\hat{V}}_{\text{p-i}}$$

$$i\hslash \frac{\partial }{\partial t}|\Psi (t)⟩=\hat{H}(t)|\Psi (t)⟩,$$

$$\Psi (x_{\text{p}},x_{\text{i}},s;t)=\sum _{n,m,p}{C}_{n,m,p}(t){\psi }_n(x_{\text{p}}){\varphi }_m(x_{\text{i}}){\chi }_p.$$

$$\Psi (x_{\text{p}},s;t)=\sum _{n,p}{C}_{n,p}(t){\psi }_n(x_{\text{p}}){\chi }_p.$$

$$\hat{H}=$$

$$+\frac{\hslash {\Omega }_R}{2}{\hat{\sigma }}_x$$

$$+\sum _{s=\uparrow ,\downarrow }{\upsilon }_{1D}^{e,s}(x_{\text{p}})+{\upsilon }_{1D}^{o,s}(x_{\text{p}})$$

$$=$$

$$\left[\begin{array}{cccc}{E}_L^0+J_{LL}^{\downarrow }& J_{LR}^{\downarrow }& \frac{\hslash {\Omega }_R}{2}& 0\\ J_{RL}^{\downarrow }& E_R^0+J_{RR}^{\downarrow }& 0& \frac{\hslash {\Omega }_R}{2}\\ \frac{\hslash {\Omega }_R}{2}& 0& E_L^0+J_{LL}^{\uparrow }& J_{LR}^{\uparrow }\\ 0& \frac{\hslash {\Omega }_R}{2}& J_{RL}^{\uparrow }& E_R^0+J_{RR}^{\uparrow }\end{array}\right],$$

$$J_{NM}^s$$

$$=⟨N|g_{1D}^{e,s}{\hat{\delta }}_{\pm }|M⟩+⟨N|g_{1D}^{o,s}{\delta }^{\prime }(x_{\text{p}}){\hat{\partial }}_{\pm }|M⟩$$

$$=g_{1D}^{e,s}\cdot N(0)\cdot [M(0^{+})+M(0^{-})]/2$$

$$-g_{1D}^{o,s}\cdot N^{\prime }(x_{\text{p}}){|}_{x_{\text{p}}=0}\cdot [M^{\prime }(0^{+})+M^{\prime }(0^{-})]/2$$

$${\hat{H}}_{\text{i}}(t)=-\frac{{\hslash }^2}{2m_{\text{i}}}\frac{{\partial }^2}{\partial x_{\text{i}}^2}+\frac{m_i{\Omega }^2(t)}{8}{x}_{\text{i}}^2$$

$${\hat{H}}_{\text{i}}(t)=-\frac{{\hslash }^2}{2m_{\text{i}}}\frac{{\partial }^2}{\partial x_{\text{i}}^2}+\frac{m_{\text{i}}{\omega }_{\text{i}}^2}{2}{x}_{\text{i}}^2+\frac{m_{\text{i}}{\omega }_{\text{i}}^2}{2}{m}_{\text{i}}{x}_{\text{i}}^2\left[\frac{\Omega (t{)}^2}{4{\omega }_{\text{i}}^2}-1\right]$$

$$\hat{H}=-J\sum _k\left[{\hat{b}}_k^{†}{\tilde{\sigma }}_{k,k+1}^{+}{\hat{b}}_{k+1}+\text{h.c.}\right]+m\sum _k{\left(-1\right)}^k{\hat{b}}_k^{†}{\hat{b}}_k.$$

$$V_{\text{p}}^{\mathrm{ext}}(x)=\{\begin{array}{ll}\frac{1}{2}{m}_{\text{p}}{\omega }_L^2(x+d_L{)}^2,& \text{if }x<0\\ \frac{1}{2}{m}_{\text{p}}{\omega }_R^2(-x+d_R{)}^2+\Delta ,& \text{otherwise}\end{array}$$

$$-\frac{{\hslash }^2}{2m_{\text{p}}}\frac{{\partial }^2\psi (x)}{\partial x^2}+V_{\text{p}}^{\mathrm{ext}}(x)\psi (x)=E\psi (x)$$

$$\frac{{\partial }^2\psi (x)}{\partial x^2}-\left[(x+d_L{)}^2-2E\right]\psi (x)=0.$$

$$\frac{{\partial }^2\psi (x)}{\partial x^2}-\left[{\left(x+d_R\right)}^2-2\left({\nu }_1+\frac{1}{2}\right)\right]\psi (x)=0.$$

$$\frac{{\partial }^2\psi (y)}{\partial y^2}-\left[\frac{y^2}{4}-\left({\nu }_1+\frac{1}{2}\right)\right]\psi (y)=0.$$

$$\frac{{\partial }^{2}w}{\partial z^2}-\left(\frac{z^2}{4}+a\right)w=0,$$

$$\frac{\partial D_{\nu }}{\partial z}$$

Full text

I Introduction

Quantum simulation aims at developing devices that can emulate physical systems and experimental conditions that are still poorly understood (e.g., high-temperature superconductivity and quark confinement). It was Feynman that originally suggested to use single purpose quantum computers to simulate a quantum system of interest, which is hardly controllable Feynman (1982). In this regard, trapped ions and quantum degenerate atomic gases are ideal candidates, particularly because of the excellent controllability of the relevant system parameters, e.g., inter-particle interactions and coupling to external fields. Indeed, nowadays atomic quantum simulation of various lattice models relevant in condensed-matter physics is a well-established research area Bloch et al. (2012); Blatt and Roos (2012). A new frontier is the quantum simulation of field theories Casanova et al. (2011); García-Álvarez et al. (2015) and strongly correlated systems with dynamical gauge fields, e.g., lattice gauge theories, which are non-perturbative formulations of high-energy physics models such as quantum electrodynamics (QED). Numerical methods that simulate these models range from traditional quantum Monte-Carlo algorithms Wilson (1974); Kogut (1983) to tensor networks Tagliacozzo et al. (2013); Haegeman et al. (2010); Rico et al. (2014); Tagliacozzo et al. (2014); Bañuls et al. (2015); Pichler et al. (2016). There are problems, however, such as the real-time simulation of heavy ion collisions, that remain intractable with such methods. Hence, simulating such problems with quantum systems specifically tailored for this purpose is necessary. In this regard, let us note that very recently an experiment with a small trapped ion quantum computer has simulated the Schwinger pair production mechanism Martinez et al. (2016), thus demonstrating that such quantum simulations are within reach.

The degrees of freedom of usual gauge theories are described by matter fields $\hat{b}_{k}$ and gauge fields $\hat{U}_{k,k+1}$, where $k$ points to a vertex in a regular lattice, while for the sake of simplicity, we refer to $k+1$ its nearest neighbour vertex (see also Fig. 1). The former fields can be, for instance, fermionic, if we are characterising quark matter, or bosonic, if we are interested in Higgs physics Kasamatsu et al. (2013); Edmonds et al. (2013); Bazavov et al. (2015); González-Cuadra et al. (2017). The dynamics of matter and gauge fields is invariant under local gauge transformations and, assuming the locality of the interactions, the most relevant term is $\hat{b}^{{\dagger}}_{k}\hat{U}_{k,k+1}\hat{b}_{k+1}$ + H.c., i.e., a correlated hopping of the matter field mediated by the excitation of the gauge field Horn (1981); Orland and Rohrlich (1990); Chandrasekharan and Wiese (1997) (e.g., in QED $\hat{U}_{k,k+1}$ is the exponential of the vector potential Kogut and Stephanov (2003)).

Up to now, there are two strategies to build an analog lattice gauge quantum simulator: (a) The (local) gauge symmetry is imposed as a constraint, that is, the gauge variant interactions are cancelled by a large energy penalty Banerjee et al. (2012); Marcos et al. (2013) in such a way that the unphysical gauge variant states are very unlikely populated; (b) the local symmetry is mapped into a fundamental symmetry of the system Zohar et al. (2016). In the first strategy, the symmetry is not a fundamental one for the simulator as in (b), but it arises in the low-energy sector of the simulated Hamiltonian such that gauge variant interactions are cancelled by a large energy penalty. For instance, in the proposal of Ref. Banerjee et al. (2012) for the quantum simulation of QED, the boson-boson and fermion-boson interactions are chosen to be of the same order such that all gauge variant states are removed from the low-energy sector of the system Hilbert space, and thus fulfilling gauge invariance. Here, we aim at realising a similar gauge invariant Hamiltonian by means of the global symmetry, i.e., strategy (b), which offers interesting advantages, as we will show in the following.

To this end, in the next section we introduce the basic formulation of a lattice gauge theory that we use, the so-called Quantum Link Models and we describe how the Higgs or Meissner effect would appear in this context. Then in Sec. III, we summarise the main ingredients of the atomic physics setup needed for the quantum simulation of this type of models via state-dependent hopping. To provide a fine-grained characterisation of the experimental architecture, we develop the microscopic description of the main interactions between the atoms (i.e., the matter field) and the impurities (i.e., the gauge field) that give as a result a state-dependent hopping. We analyse several cases: from a simplified model with a four-level system, through the neutral atom impurity situation, ending at the ion impurity setup. In Sec. V we discuss a many-body problem where the proposed quantum simulator could be an efficient and useful tool, we follow the time evolution of a quenched experiment in lattice gauge theories. Finally, in Sec. VI we discuss parameters for experimental implementations, while in Sec. VII we summarise our findings and provide an outlook. Technical derivations are provided in an appendix.

II Quantum link formulation of a $U(1)$ lattice gauge model

In a $U(1)$ lattice gauge model, the hopping of the matter field is mediated by the gauge field via the term, $\hat{b}^{\dagger}_{k}\hat{U}_{k,k+1}\hat{b}_{k+1}$, which represents the minimal coupling of the matter and gauge degrees of freedom. The fundamental gauge degrees of freedom $\hat{U}_{k,k+1}$ describe quantum operators associated with the lattice links. The complete Hamiltonian, up to the magnetic term, is given by Kogut and Stephanov (2003)

[TABLE]

Here $\hat{E}_{k,k+1}$ is the electric field operator, $g$ the gauge coupling, $m$ the staggered mass and $J$ the hopping parameter. In the Hamiltonian, the first term corresponds to hopping of the matter field between two adjacent lattice sites, which is associated with a change in the gauge field when the charge moves sites. The second term in the equation, the mass term, defines the staggered matter field, where excitations of these sites corresponds to the creation of a particle-antiparticle pair with a mass gap of $m$. The last term corresponds to the electric field energy.

The electric field and the link operator are conjugate operators, which obey $[\hat{E}_{k,k+1},\hat{U}_{k,k+1}]=\hat{U}_{k,k+1}$. Quantum link models Horn (1981); Orland and Rohrlich (1990); Chandrasekharan and Wiese (1997) is a framework for lattice gauge theories where the dynamical gauge fields are represented by systems of discrete quantum degrees of freedom with only a finite-dimensional Hilbert space per link. The previous commutation relation is realised by a quantum link operator $\hat{U}_{k,k+1}=\hat{S}^{+}_{k,k+1}$, which is a raising operator for the electric flux $\hat{E}_{k,k+1}=\hat{S}^{z}_{k,k+1}$ associated with the link connecting neighbouring lattice sites $k$ and $k+1$. Hence, the local gauge Hilbert space is spanned by just $(2S+1)$ quantum states per link. This is in contrast with the original Hamiltonian formulation Kogut and Stephanov (2003), where $\hat{U}_{k,k+1}$ is defined on an infinite dimensional Hilbert space.

A $U(1)$ gauge invariant model is one that has as a local symmetry the conservation of the electric flux around any vertex in the lattice. In any spatial dimension $D$, this constraint reads $\sum^{D}_{\hat{\mu}=1}\left(\hat{E}_{k,k+\hat{\mu}}-\hat{E}_{k-\hat{\mu},k}\right)=\hat{\rho}_{k}$, which means that the difference between the electric fluxes that enter in a vertex minus the electric fluxes that go out from the vertex are equal the charge at this vertex, where $\hat{\mu}$ is the unit vector in the positive $\mu$ direction. This is nothing but the Gauss law of QED that in the continuum is given by $\vec{\nabla}\vec{E}(x)=\rho(x)$. The generator of the $U(1)$ transformation

[TABLE]

is a constant of motion and, as a symmetry generator, commutes with the Hamiltonian, i.e. $[\hat{H},\hat{G}_{k}]=0,~{}\forall k$. Once we have defined the generator of the local $U(1)$ transformation, it is straightforward to check the action of the unitary operator on the relevant fields, for instance,

[TABLE]

For the implementation that we describe, we use a spin-$1/2$ representation of the quantum link in a rotated basis, where $\hat{E}\mapsto\hat{\sigma}^{x}/2$ and $\hat{U}\mapsto\tilde{\sigma}^{+}=|+\rangle\langle-|$ with $|\pm\rangle$ being the eigenstates of the Pauli matrix $\hat{\sigma}_{x}$, so the interaction term between the matter and gauge degrees of freedom will be of the form, $\hat{b}^{\dagger}_{k}\tilde{\sigma}^{+}_{k,k+1}\hat{b}_{k+1}$, and the electric field energy $\hat{E}^{2}\propto\left[\hat{\sigma}^{x}\right]^{2}=\mathbb{I}$. In this regard, let us note that in actual QED the electric field and the gauge degrees of freedom are indeed continuum variables. However, in recent numerical studies Kühn et al. (2014); Buyens et al. (2017), it has been shown that already at small dimensions (on the order of 5), the continuum limit can be reproduced with good accuracy. Hence, by choosing appropriately a manifold with five Zeeman or hyperfine levels, as shown recently in Ref. Kasper et al. (2017), it is possible, also with our scheme, to quantum simulate the Schwinger model in the continuum.

II.1 Bosonic Josephson junctions: A toy model for the Higgs physics

Due to the bosonic nature of the matter degrees of freedom in the previous $U(1)$ quantum link model, an interesting scenario appears when there is an ensemble of interacting bosonic atoms at every lattice site, for previous implementations of this scenario see Refs. Kasamatsu et al. (2013); Edmonds et al. (2013); Bazavov et al. (2015); González-Cuadra et al. (2017). For a large number of bosons per site, we can assume $b_{k}\sim\sqrt{n_{k}}e^{i\theta_{k}}$ with $n_{k}$ the average number of particles within each condensate, which can be fixed to a constant $\frac{N}{2}$, and $\theta_{k}$ the condensate phase. Thus, the first term in Eq. (1) becomes (see also Appendix B for details):

[TABLE]

This Hamiltonian describes the dynamics of a Bose-condensed matter field interacting with a $U(1)$ gauge field. In such a situation, due to a combination of the spontaneous symmetry breaking of the bosonic field and the gauge interaction, the gauge field acquires a mass, as we will see in the following.

Let us underscore, however, that contrarily to the large particle number of the gauge field discussed in Ref. Kasper et al. (2017), where the goal is to attain the continuum limit of the gauge field, here we aim at considering large ensembles of particles and weak interactions in the matter field. In such a way the quantum simulator that we are going to introduce in the next sections would in principle enable the investigation of different high-energy and condensed matter physics problems. Indeed, the condensate wave function would represent the bosonic Higgs-field, and thus the Higgs model in high-energy physics or equivalently the Meissner effect in superconductors Fradkin (2013) could we explored.

To be specific, a particular way to obtain a gauge theory with non-zero mass is via the so-called Higgs mechanism, a combination of spontaneous symmetry breaking and gauge symmetry, where the interaction of the Goldstone boson and the gauge boson give rise to a massive theory, i.e., short-range correlated. A simplified model of the spontaneous symmetry breaking is given by the Hamiltonian: $\hat{H}=-J\sum_{k}\hat{b}^{\dagger}_{k+1}\hat{b}_{k}+V\left(\hat{\rho}_{k}\right)+\text{h.c.}$, where $\hat{b}_{k}$ are usual bosonic operators at site $k$ and $\hat{\rho}_{k}=\hat{b}^{\dagger}_{k}\hat{b}_{k}$. In our setting the density-dependent potential is given by $V\left(\hat{\rho}_{k}\right)=\sum_{k}\epsilon_{k}\hat{\rho}_{k}+U/2\sum_{k}\hat{\rho}_{k}(\hat{\rho}_{k}-1)$, where the first term represents the energy due to the confinement and the second one the interaction between the bosons. If there is a non-trivial potential that depends only on the density of bosons $V\left(\rho_{k}\right)$, as in our physical setup, we can expand the bosonic operators around the potential minima, i.e., \frac{\partial}{\partial\rho}V\big{|}_{\rho=\rho_{0}}=0 and $b_{k}\sim\sqrt{\rho_{0}}e^{i\theta_{k}}$. Such approximation holds definitely in the large particle limit and for weak inter-particle interactions. Under these assumptions, the Hamiltonian is recast into: $\hat{H}=2\rho_{0}\sum_{k}\cos{\left(\theta_{k+1}-\theta_{k}\right)}\sim\int dx\left[\frac{\partial}{\partial x}\theta(x)\right]^{\dagger}\frac{\partial}{\partial x}\theta(x)$. Here, we assumed that the potential $V$ is around its minimum value and, in the last equality of the Hamiltonian, we assumed that the lattice spacing is small enough as well as the function $\theta(x)$ is continuous and differentiable. When the ground state of the system gets a non-trivial density expectation value, a massless (Goldstone) boson appears due to the spontaneous global symmetry breaking.

Now, if the bosonic degrees of freedom are coupled to a gauge field we have: $\hat{H}=\sum_{k}J_{k}\hat{b}^{\dagger}_{k+1}e^{iga\hat{A}_{k+1,k}}\hat{b}_{k}+V\left(\hat{\rho}_{k}\right)+\text{h.c.}$, where we assume again that the potential $V(\hat{\rho}_{k})$ depends on the boson density only, $a$ is the lattice spacing and $g$ the coupling constant. Before expanding the interacting Hamiltonian, we perform a local gauge transformation such that: $\hat{\tilde{b}}_{k}=\hat{b}_{k}e^{-i\theta_{k}}$ and $ga\tilde{A}_{k+1,k}=ga\hat{A}_{k+1,k}+\theta_{k}-\theta_{k+1}$, then the interacting Hamiltonian reduces to a mass term in the gauge field, $\hat{H}=g^{2}\rho_{0}\int dx~{}\hat{A}(x)^{2}$. This mechanism is known as Higgs effect in particle physics and Meissner effect in condensed matter physics, and it could be quantum simulated with our proposed implementaion.

III System, Target Hamiltonian and Conditions for quantum simulation

The system under investigation is illustrated schematically in Fig. 1 a). A superlattice (blue line) with a periodic sequence of double well-like potentials, each trapping neutral atoms (hereafter called particles), is formed. In another periodic potential (green line), which does not influence the former superlattice, either neutral or charged atoms (hereafter called impurities) with two addressable internal states $|\uparrow\rangle,\,|\downarrow\rangle$ are trapped at the position of each double well potential barrier. The particles represent the matter field, while the impurities play the role of the dynamical gauge field. In the limit where the dynamics of the impurities is frozen (i.e., tight trapping), the particle-impurity interaction can be treated as an external potential for the particles. As shown in Refs. Micheli et al. (2004); Gerritsma et al. (2012); Schurer et al. (2016), the internal state of the impurity can control the hopping, $J$, of a particle confined, e.g., in a double well potential, via the particle-impurity interaction, and therefore be promoted to an operator $\hat{J}=J^{\downarrow}|\downarrow\rangle\langle\downarrow|+J^{\uparrow}|\uparrow\rangle\langle\uparrow|$. We note that the physical origin of such state dependence is due to the reliance of the particle-impurity scattering length on the electronic internal states, especially at the short-range part of the interaction potential when the particle and impurity electronic clouds do overlap in a collision. Such controlled hopping can be generalised to a lattice, as depicted in Fig. 1 a), in such a way that, within the tight binding approximation, the following Hamiltonian for the particles holds

[TABLE]

Here $\hat{N}=\sum_{k}\hat{n}_{k}=\sum_{k}\hat{b}^{{\dagger}}_{k}\hat{b}_{k}$ is the particle number, $\hat{J}_{k,k+1}$ is the state-dependent hopping, $\hat{U}_{k}$ denotes the on-site energy, whose state-dependence is very weak Gerritsma et al. (2012); Negretti et al. (2014).

Our goal is to show that the first sum in Eq. (4) can be mapped into $\propto\sum_{k}\hat{b}^{\dagger}_{k}\,\tilde{\sigma}^{+}_{k,k+1}\hat{b}_{k+1}+$ H.c. with $\tilde{\sigma}^{+}=|+\rangle\langle-|$, where the eigenstates of the Pauli matrix $\hat{\sigma}^{x}$ are given by $|\pm\rangle=(|\uparrow\rangle\pm|\downarrow\rangle)/\sqrt{2}$, while the ones of $\hat{\sigma}^{z}$ by $|\uparrow\rangle$ and $|\downarrow\rangle$, respectively. Hence, if a particle tunnels from the site indexed by $k$ to the site indexed by $k+1$, then the internal state of the impurity is flipped [see also panel c) of Fig. 1]. This interaction is gauge-invariant and of the form $\hat{b}^{{\dagger}}_{k}\,\hat{U}_{k,k+1}\hat{b}_{k+1}$ + H.c.

To demonstrate this, let us for the moment ignore the on-site terms and focus on a single double well only. We start from a microscopic Hamiltonian $\hat{H}_{(k,k+1)}=\hat{H}^{0}_{(k,k+1)}+\hat{H}^{\text{hopp}}_{(k,k+1)}$ that is split in a local part and an interactive part, where

[TABLE]

We note that with such a notation the entire Hamiltonian on the lattice reads

[TABLE]

The natural energy scales of each mode are $E_{k}^{0}$ for the particle at the site indexed by $k$ (i.e., the lowest trapping energy of the particle in the $k$-th site), whereas $\Omega_{R}$ is the Rabi frequency of the coupling between the impurity internal states (see also Fig. 1). The hopping part of Eq. (4) can be rewritten as

[TABLE]

where $J_{k,k+1}^{0}=(J^{\uparrow}_{k,k+1}+J^{\downarrow}_{k,k+1})/2$ and $J_{k,k+1}^{z}=(J^{\uparrow}_{k,k+1}-J^{\downarrow}_{k,k+1})/2$. Moving to the interaction picture we obtain:

[TABLE]

By assuming resonant coupling $\Delta=E_{k}^{0}-E_{k+1}^{0}=\hbar\Omega_{R}$ and by applying the rotating-wave approximation, $\{|\Delta|,\,\hbar\Omega_{R}\}\gg\{J_{k,k+1}^{0},\,J_{k,k+1}^{z}\}$, we arrive at the desired target Hamiltonian

[TABLE]

Gauge variant term are negligible in the limit where the rotating-wave approximation is valid, for instance, terms like $J_{k,k+1}^{0}e^{i\Delta t}\hat{b}^{\dagger}_{k}\,\hat{b}_{k+1}$ or $J_{k,k+1}^{z}e^{i2\Delta t}\hat{b}^{\dagger}_{k}\,\tilde{\sigma}^{-}_{k,k+1}\hat{b}_{k+1}$ oscillate with a high frequency $\Delta$ compared with the energy scales $J_{k,k+1}^{0}$ or $J_{k,k+1}^{z}$, which give a negligible average contribution in a long-term limit. Finally, let us note that since the onsite interaction is of a density-density type, such interaction does not affect the gauge invariance; and in addition, its energy scale and state dependence are typically small Gerritsma et al. (2012); Negretti et al. (2014).

In order to prove the quantum simulation of our scheme, two features have to be checked: (i) the robustness of the local gauge symmetry, and (ii) the correlated hopping of the impurity and the particle modes. To prove the former, we have to evaluate the local generator of the gauge transformation. Being a symmetry, such a generator commutes with $\hat{H}^{\text{t}}_{(k,k+1)}$. From the general expression Eq. (2), it can be easily verified for a single link $(k,k+1)$ that the operators $\hat{G}_{k}=\hat{n}_{k}-\hat{\sigma}^{x}_{k,k+1}/2$ and $\hat{G}_{k+1}=\hat{n}_{k+1}+\hat{\sigma}^{x}_{k,k+1}/2$ are constants of motion and local generators of the $U(1)$ (phase) transformation of the local operators. Indeed, they are the terms that act non-trivially on the link $(k,k+1)$ from Eq. (2). We note that in the QED context $\hat{\sigma}^{x}_{k,k+1}/2$ plays the role of the electric field in-between two sites. Hence, by keeping track of the expectation values of $\hat{G}_{k}$ and $\hat{G}_{k+1}$ we can assess the robustness of the gauge invariance, i.e., to which extent the expectation values of these operators evolve in time. Since the effective dynamics of the system corresponds to correlated hopping of the impurity and the particle modes, the feature (ii) can be assessed by looking at correlation functions, for instance, $\mathcal{C}_{k,k+1}(t)=\langle(\hat{b}^{\dagger}_{k}\hat{b}_{k}-\hat{b}^{\dagger}_{k+1}\hat{b}_{k+1})\hat{\sigma}_{k,k+1}^{x}\rangle-\langle\hat{b}^{\dagger}_{k}\hat{b}_{k}-\hat{b}^{\dagger}_{k+1}\hat{b}_{k+1}\rangle\langle\hat{\sigma}_{k,k+1}^{x}\rangle$. Indeed, this observable can distinguish a correlated particle-impurity hopping from the independent particle hopping and impurity internal state flip.

IV Model system: The double-well

In this section, we consider the microscopic description of the proposed quantum simulator, describing the dynamics of the atom and the impurity and the tools to solve the time-dependent Schrödinger equation. At the beginning, we analyse a simplified four-level system, where we find all the ingredients needed for the quantum simulation of the abelian lattice gauge model via state-dependent hopping. Then, we consider the full dynamics for both a neutral particle and impurity, where the former is trapped in a double well and the latter in a harmonic trap. Finally, we consider the case of a charged impurity, in particular the case of atom-ion quantum systems, and demonstrate the feasibility of our scheme with such compound system as well.

IV.1 The Hamiltonian

We shall verify to which extent $\hat{H}^{\text{t}}_{(k,k+1)}$ can be accomplished by simulating the particle dynamics in a double well interacting with an impurity. We note that the double well represents the basic building block of our quantum simulator (see also Fig. 1). Specifically, we consider the one-dimensional Hamiltonian $\hat{H}_{\text{p-i}}=\hat{H}_{\text{p}}+\hat{H}_{\text{i}}+\hat{V}_{\text{p-i}}+\hat{H}_{\text{i}}^{\text{int}}$, where

[TABLE]

Here $V_{\text{i}}^{\mathrm{ext}}(x_{\text{i}})=\frac{m_{\text{i}}\omega_{\text{i}}^{2}}{2}x_{\text{i}}^{2}$ is the confining potential for the impurity, $V_{\text{p}}^{\mathrm{ext}}(x_{\text{p}})=m_{\text{p}}\omega_{L}^{2}(x_{\text{p}}+d_{L})^{2}[1-\Theta(x_{\text{p}})]/2+[m_{\text{p}}\omega_{R}^{2}(-x_{\text{p}}+d_{R})^{2}/2+\Delta]\Theta(x_{\text{p}})$ is the double well potential for the bosons with $\Theta(x_{\text{p}})$ the Heaviside function, $m_{\text{p}}$ the particle mass, $m_{\text{i}}$ the impurity mass, and $\Delta$ the energy offset between the two wells (see also Fig. 1). Furthermore, the particle-impurity interaction is modelled with a contact potential 111We underscore that such contact potentials are applicable when the effective range of the particle-impurity interaction is smaller than any other length scale involved in the problem. with even $\upsilon_{\mathrm{1D}}^{e,s}(x)=g_{\mathrm{1D}}^{e,s}\hat{\delta}_{\pm}(x)$ and odd $\upsilon_{\mathrm{1D}}^{o,s}(x)=g_{\mathrm{1D}}^{o,s}\delta^{\prime}(x)\hat{\partial}_{\pm}$ terms Girardeau and Olshanii (2004); Valiente and Zinner (2015), whereas the operators $\hat{\delta}_{\pm}(x)$ and $\hat{\partial}_{\pm}$ are defined in the Appendix A. Here $g_{\mathrm{1D}}^{e,s}=-\hbar^{2}/\mu a_{\mathrm{1D}}^{e,s}$ and $g_{\mathrm{1D}}^{o,s}=-\hbar^{2}a_{\mathrm{1D}}^{o,s}/\mu$ with $\mu$ the reduced particle-impurity mass, and $a_{\mathrm{1D}}^{e,s},\,a_{\mathrm{1D}}^{o,s}$ are the 1D spin-dependent scattering lengths for even- and odd-waves, respectively. While the even term has the familiar form of the contact interaction between two ultra-cold bosons, the odd term appears when the interacting particles are spin-polarised fermions or distinguishable particles. As a consequence, the general solution to the scattering problem is a linear combination of even and odd parity solutions.

IV.2 Solution to the time-dependent Schrödinger equation for the double well

In general, the time evolution of a state vector is described by the time-dependent Schrödinger equation:

[TABLE]

which has the solution $|\Psi(t)\rangle=\hat{U}(t)|\Psi(0)\rangle$ with $|\Psi(0)\rangle$ and $\hat{U}(t)$ being the initial state and the time evolution operator, respectively. In order to solve the corresponding system dynamics practically, we expand the particle-impurity wave function as:

[TABLE]

Here $C_{n,m,p}(t)$ are time-dependent coefficients, $\psi_{n}(x_{\text{p}})$ are time-independent basis functions for the particle Hilbert space, $\phi_{m}(x_{\text{i}})$ are time-independent basis functions for the impurity Hilbert space, and $\chi_{p}$ denotes the spinor component of the impurity internal state. By inserting the ansatz (11) into the Schrödinger equation (10) and projecting onto the state $\langle\chi_{p},\phi_{m},\psi_{n}|$, we obtain a set of coupled differential equations for the coefficients $C_{n,m,p}(t)$ that have then to be solved numerically.

IV.3 A simplified model: The four-level system

In the following we consider the case for which the Hamiltonian is time-independent and the impurity is treated statically (i.e., no impurity motion, but only spin dynamics). Thus, the above outlined ansatz reduces to

[TABLE]

In particular, we shall focus on the special case for which only two motional states of the particle in the double well are considered: One particle-state in the left well (the ground state), $|L\rangle$, one particle-state in the right well (the first excited state of the double well), $|R\rangle$ (see also Fig. 2). In addition to these states, we have the two impurity spin states $\{|\uparrow\rangle,|\downarrow\rangle\}$ (see also Fig. 2). Under these assumptions, we have in total four possible state configurations. Moreover, we note that the particle-impurity interaction contributes to the particle Hamiltonian as an external contact potential given by Eq. (IV.1) centred in $x_{\text{i}}=0$, where the impurity is assumed to be located (see also Fig. 2). Thus, the particle-impurity Hamiltonian is given by

[TABLE]

Hence, we form the basis $\{|L,\downarrow\rangle,|R,\downarrow\rangle,|L,\uparrow\rangle,|R,\uparrow\rangle\}$, where $\hat{H}_{\text{p}}|L\rangle=E_{L}^{0}|L\rangle$ and $\hat{H}_{\text{p}}|R\rangle=E_{R}^{0}|R\rangle$. Hence, in this basis, the particle-impurity Hamiltonian has the following matrix form:

[TABLE]

with

[TABLE]

for $s=\uparrow,\downarrow$ and $N,M=L,R$. Here the apex ′ denotes the spatial derivative. Exemplary values of the parameters previously introduced are: $g_{\mathrm{1D}}^{o,\uparrow}=\hbar\omega_{L}\ell^{3}$ and $g_{\mathrm{1D}}^{e,\uparrow}=\hbar\omega_{L}\ell$ imply $J^{\uparrow}_{LR}=J^{\uparrow}_{RL}=-0.011\hbar\omega_{L}$, whereas $g_{\mathrm{1D}}^{o,\downarrow}=-\hbar\omega_{L}\ell^{3}$ and $g_{\mathrm{1D}}^{e,\downarrow}=-\hbar\omega_{L}\ell$ imply $J^{\downarrow}_{LR}=J^{\downarrow}_{RL}=-0.011\hbar\omega_{L}$. Note that if the wells are too separated, then the hopping rates will be exponentially suppressed, and therefore no tunnelling occurs.

Now, we have simulated the corresponding dynamics by assuming that for $t=0$ the initial state is $|\Psi(t=0)\rangle=\frac{|R,{\uparrow}\rangle-|R,{\downarrow}\rangle}{\sqrt{2}}\equiv|R,-\rangle$. Hence, the only non-vanishing coefficients are: $C_{R,{\downarrow}}=-1/\sqrt{2}$ and $C_{R,{\uparrow}}=1/\sqrt{2}$. By solving the corresponding coupled differential equations for the coefficients $C_{n,p}(t)$, we can extract the probability of finding the particle in the left well with the flipped spin state for the impurity as $\mathcal{O}^{+}_{L}(t)=|C_{L,\uparrow}(t)+C_{R,\downarrow}(t)|^{2}/2$ [similarly for $\mathcal{O}^{-}_{R}(t)$]. Examples of such calculations for different interaction strengths are illustrated in Fig. 3. As one can see, in the case the couplings for the two spin states have opposite sign, but have the same strength (green lines in Fig. 3), the particle-impurity state oscillates purely between the states $|R,-\rangle$ and $|L,+\rangle$. This situation represents the ideal scenario we are aiming for in our study.

Furthermore, the overall picture is maintained if the strengths are different (not shown). On the other hand, if the couplings have the same sign (dashed blue lines), the process is significantly modified with respect to the desired ideal scenario (green solid lines). This is because the resonance condition has to be replaced with $\hbar\Omega_{R}=E^{0}_{R}-E^{0}_{L}+(J^{\downarrow}_{RR}+J^{\uparrow}_{RR}-J^{\downarrow}_{LL}-J^{\uparrow}_{LL})/2$. Then, the process works again very well as in the other cases (see dash-dotted magenta lines). Hence, the desired process can take place regardless the values of the coupling constants for the even and odd wave and for the two spin-states of the impurity. Thus, this proves that our schemes does not require fine tuning of the inter-particle interactions.

IV.4 Simulation for neutral impurity-particle systems

Upon the previous findings on the four-level system, we investigate now the full dynamics for a neutral impurity and a particle in the double well potential. We begin by consider a tightly trapped impurity such that $\hat{H}_{\text{i}}\equiv 0$ (i.e., static impurity at $x_{\text{i}}=0$). In this case the eigenvalue problem for the particle is analytically solvable, as shown in the Appendix A. We assume a well separation $d_{L}+d_{R}$ such that $J_{\uparrow,\downarrow}$ is not too small, but still smaller than $\Delta$ and $\hbar\Omega_{R}$. This means that in such a configuration the particle and the impurity are interacting. As initial condition for the time-dependent Schrödinger equation we choose the particle state $|R\rangle$, i.e., the first excited eigenstate of $\hat{H}_{\text{p}}$, and for the impurity the state $|-\rangle$ (see Fig. 1). We remind that the goal is to obtain, after a certain time, the target state $|+\rangle$ for the impurity and the particle in the left well, but not necessarily in the ground state $|L\rangle$ of $\hat{H}_{\text{p}}$. The result of this analysis is illustrated in Fig. 4, where the overlap $\mathcal{O}_{L}^{+}(t)$ is shown in the upper panel. Here $\mathcal{O}_{y}^{\alpha}(t)=\int_{\mathbb{R}^{y}}dx_{\text{p}}\int_{\mathbb{R}}dx_{\text{i}}|\langle x_{\text{p}},x_{\text{i}},\alpha_{\text{i}}|\psi_{\text{p-i}}(t)\rangle|^{2}$ with $\alpha=\pm$, $y=L,R$, $\mathbb{R}^{R,L}\equiv\mathbb{R}^{\pm}$, and where $\psi_{\text{p-i}}(t)$ is the time evolved particle-impurity state. As the solid lines show, after a time $t_{\max}\approx 104/\omega_{L}$ the desired target state is reached [$\mathcal{O}_{L}^{+}(t_{\max})\simeq 0.94$]. Further, in the lower panel the time evolution of the correlation $\mathcal{C}(t)\equiv\mathcal{C}_{L,R}(t)$ is displayed with revivals in agreement with the evolution of the overlaps. This clearly demonstrates that the particle hopping is correlated with the internal state flip of the impurity. The inset illustrates the expectation value of $\hat{G}_{L}$ (similarly for $\hat{G}_{R}$). It shows that the generator of the gauge transformation is, to a good approximation, a constant of motion. Hence, the microscopic dynamics of such a setup corresponds to our target Hamiltonian. In reference to the four-level system studied previously, we see that by including more states in the total wave function expansion (12), more excited states participate in the dynamics. Thus, the sine-like behaviour of the populations of the four-level scheme is a bit hampered and additional oscillations with smaller amplitudes on the top of the main oscillation appear.

The desired state transfer is also obtained when the static impurity assumption is relaxed (i.e., $\hat{H}_{\text{i}}\neq 0$). We obtain essentially the same result as for the static particle for $\omega_{\text{i}}=10^{2}\omega_{L}$ (not shown). We note that this situation can be realised experimentally by starting at $t=0$ with a height of the double well barrier such that no particle-impurity interaction takes place and then suddenly quenching it to a value for which the two systems interact and such that the conditions $|\Delta|=\hbar\Omega_{R}$ and $\{|\Delta|,\,\hbar\Omega_{R}\}\gg\{J^{0},\,J^{z}\}$ are fulfilled. This is illustrated in Fig. 4, where the desired behaviour is observed. Hence, the impurity trap does not play any major role in the dynamics (e.g., 6% overlap reduction if $\omega_{\text{i}}=\omega_{L}$).

Finally, let us note that we have also considered the feasibility of the scheme, particularly of the correlated hopping, in the case the particle system in the double well is replaced by a condensate in contact with a two-level spin impurity. We found that even within mean-field theory for the bosonic ensemble and by means of optimal control techniques such a process is not possible and we attribute this failure mainly to nonlinear effects of inter-particle interactions. A detailed analysis together with the derivation of the corresponding equations of motion is provided in the Appendices C and D.

IV.5 Simulation for the atom-ion system

The situations described so far are representatives of an ultra-cold neutral particle-impurity system, such as an atomic boson and fermion. We have also considered the case in which the impurity is replaced by an ion. The system then resembles the ion-controlled double well of Ref. Joger et al. (2014). Here, however, we assume that the long-range atom-ion polarisation potential can be approximated by the contact interaction (IV.1), which is reasonable if the separation between the ions is larger than the range of the atom-ion interaction Negretti et al. (2014). Additionally, we have included the effect of the time-dependent radio-frequency (rf) fields, since most of the ion experiments use the so-called Paul-trap Leibfried et al. (2003). In this way we can assess the impact of the micro-motion on the dynamics. More precisely, we have considered the Hamiltonian for an ion in one spatial dimension subjected to static and time-dependent electric fields given by Cook et al. (1985)

[TABLE]

where $\Omega(t)^{2}=\Omega_{\text{rf}}^{2}\cdot[a+2q\cos(\Omega_{\text{rf}}t)]$ with $\Omega_{\rm rf}$ the driving frequency, and $|a|\ll|q|<1$ being some trapping geometrical factors (typical values used in experiments are $|a|\ll|q|\approx 0.2$ for a linear Paul trap). For the numerical simulations, however, it is more convenient to recast it as:

[TABLE]

where $\omega_{\text{i}}=\frac{\Omega_{\text{rf}}}{2}\sqrt{a+q^{2}/2}$ is the so-called secular frequency, that is, the effective harmonic trap frequency of the ion, if the micro-motion would be neglected. Thus, we chose the basis functions $\phi_{m}(x_{\text{i}})$ for the expansion of the particle-ion wave function (11) as the eigenstates of the harmonic oscillator with frequency $\omega_{\text{i}}$.

The result is shown in Fig. 4 (dotted lines) for a driving frequency $\Omega_{\text{rf}}/\omega_{L}=2500$. Similarly to the quenched barrier height case, more pronounced wiggles are observed, but the overall impact is in no way detrimental. Hence, the analysis demonstrates that also with the compound atom-ion system the scheme can be implemented, especially for pairs with a small atom-ion mass ratio (e.g., Li-Yb+), as the influence of the micro-motion Cetina et al. (2012); Joger et al. (2014) and the spin-orbit coupling Tscherbul et al. (2016) is expected to be small. We note that very recently experimental measurements on charge exchange show that the corresponding rate is significantly smaller than the Langevin collision rate, thus confirming that for Li-Yb+ the required coherence for achieving state-dependent tunnelling is indeed possible Joger et al. (2017).

V False vacuum simulation

An instance for which the proposed quantum simulator becomes relevant is in out-of-equilibrium dynamics, where the system is initially prepared in a quantum state that is not an eigenstate of the target Hamiltonian, e.g., due to a collision or a sudden quench of the parameters that define the dynamics. The notion of false vacuum decay is linked to out-of-equilibrium physics, in the sense, that the system is prepared in a state (the so-called false vacuum) that is not an eigenstate of the Hamiltonian with which the system evolves.

As we have seen, the precise Hamiltonian that characterises the dynamics is given by

[TABLE]

Such a model describes a bosonic matter field $\hat{b}_{k}$ with a staggered mass $m$ coupled to a gauge field $\tilde{\sigma}^{+}_{k,k+1}=|+\rangle_{k,k+1}\langle-|$. Note that the mass term in Eq. (16) is realised in our quantum simulator with the double well potential (or superlattice in a lattice system). Due to the mass term, at half filling, the vacuum of the decoupled model corresponds to a staggered distribution of particles depending on the sign of the mass (i.e., in our simulator, deeper wells correspond to a negative sign, whereas energetically higher wells to a positive sign). At this point, it is worth emphasise the role of the superlattice: In a first place, we have seen that it is needed to implement the desired Hamiltonian in an interacting picture, whereby the difference of the energy between the double wells $\Delta$ matches the local Rabi frequency $\Omega_{R}$ of the impurity that is placed on the links between the sites. Then, if the resonance condition is not completely fulfilled, the detuning between these energy levels takes the role of a staggered mass $m$ that characterises two different type of sites: Particle-sites, when the mass term is positive, $m>0$; antiparticle-sites, when the mass term is negative, $m<0$. Moreover, at half-filling, the sites with negative mass will be filled, defining the reference state of the matter field.

Further, because of the gauge invariance when we consider the gauge field configuration on these staggered particle distributions, there are three possible (ground) states:

• •

Positive mass, $(m\gg 0)$:

$|g_{0}\rangle=|\cdots 1_{2k-1},+_{2k-1,2k},0_{2k},-_{2k,2k+1},1_{2k+1},\cdots\rangle$

• •

Negative mass, $(m\ll 0)$:

$|g_{+}\rangle=|\cdots 0_{2k-1},+_{2k-1,2k},1_{2k},+_{2k,2k+1},0_{2k+1},\cdots\rangle$

$|g_{-}\rangle=|\cdots 0_{2k-1},-_{2k-1,2k},1_{2k},-_{2k,2k+1},0_{2k+1},\cdots\rangle$

The first thermodynamical property that characterises all these states is the electric flux density of the system, defined as $\frac{1}{N}\sum_{k}\hat{\sigma}^{x}_{k}$, i.e., the average value of the configurations of the state of the impurity. The first state $|g_{0}\rangle$ have zero net electric flux density, while the other two has positive, $|g_{+}\rangle$, and negative, $|g_{-}\rangle$, electric flux, respectively.

This cartoon structure of the vacuum of the system remains with the full interacting Hamiltonian (16), i.e., the vacuum of the full Hamiltonian with positive mass is unique with zero electric flux, while with negative mass the ground state manifold is two-fold degenerate with positive and negative value of the electric flux (see Fig. 5). These two behaviours, in fact, define two phases separated by a quantum phase transition at a finite negative $m$. The quantum phase transition goes from a charge and parity ordered phase with non-zero electric flux to a disordered one with a net zero electric flux configuration.

Once we have identified the different phases and different vacua identified by the charge and parity symmetries Rico et al. (2014), we could perform a quenched experiment where, the system is prepared in one of the ground states that break the symmetries and we let it evolve with a Hamiltonian from the disorder phase. Then we expect non-trivial oscillations of the order parameter, i.e., the net electric flux, and oscillations between these states. Figure 6 illustrates how the electric flux density evolves from a non-zero value to a vanishing one. This type of toy experiments are examples of out-of-equilibrium dynamics, where the initial state of a many-body system is trapped in a local energy minima after a collision or a sudden change of the thermodynamic conditions, and then, the system evolves from the highly excited state.

VI Experimental considerations

Typically, time scales of ultra-cold atom experiments do not exceed a few hundred of milliseconds. This means that in order to attain the desired state in a reasonable time, we need a trap frequency of a few kHz and similarly for the Rabi frequency. For instance, for $\omega_{L}/(2\pi)=$ 1 kHz we would have $t_{\max}\approx 104/\omega_{L}\simeq 17$ ms. In the case of a quasi-1D setting, a transverse frequency for the particle $\sim 20\omega_{L}$ is necessary to keep the transverse motion frozen to its ground state. These frequencies are attainable in optical lattice or atom-chip experiments and coherence times of a few hundred milliseconds are easily obtainable. For the atom-ion system, the condition is even more easily satisfied: For a driving frequency $\Omega_{\text{rf}}/(2\pi)=2.5$ MHz, we obtain a secular frequency $\omega_{\text{i}}/(2\pi)\simeq 254$ kHz $\gg\omega_{L}$. As far as the odd-wave coupling constants $g_{\mathrm{1D}}^{o,\uparrow,\downarrow}$ are concerned, we note that they do not play an essential role for the successful realisation of the process. Since odd-wave contributions at low temperatures are typically small, we used small values for $g_{\mathrm{1D}}^{o,\uparrow,\downarrow}$ and, in principle, they may be negligible on short time-scales. Instead, the even coupling constants $g_{\mathrm{1D}}^{e,\uparrow,\downarrow}$ play a crucial role for the implementation of our scheme. Note, however, that any kind of combination $g_{\mathrm{1D}}^{e,\uparrow,\downarrow}$ will work, a particular choice of such couplings relies essentially on how fast the process has to be with respect to the decoherence times involved in the specific atomic system used. In fact, the correlated hopping vanishes only in the special case where the couplings have the same strength and sign, since here $J_{k}^{z}$ becomes zero. We underscore, however, that the proposed scheme can be generalised to higher dimensions and suitable parameters can be determined as well.

VII Conclusions and outlook

In this work we have proposed an alternative route to the quantum simulation of lattice gauge theories by means of compound atomic quantum systems like Bose-Fermi or atom-ion systems. Our approach does not require any precise control of the particle-impurity interactions, and therefore it offers advantages compared to previous proposals Banerjee et al. (2012); Kasper et al. (2017). Furthermore, we investigated in great detail the basic building block of the simulator and demonstrated that the correlated hopping for the single particle case can be realised efficiently when the resonance condition $|\Delta|=\hbar\Omega_{R}$ is fulfilled. Interestingly, we have found that the impurity trap plays only a marginal role for the correlated hopping process and that for atom-ion systems the scheme is resilient to micro-motion. Furthermore, we have shown that in the limit of a large number of bosons the system can be used to explore the physics of the Higgs model. In the future it would be interesting to investigate more closely the impurity-induced correlated many-body quantum dynamics and/or when the initial superfluid state in the Bose-Hubbard is suddenly perturbed by the presence of the impurities via state-dependent hopping terms as well as gauge plaquette terms that appear in higher spatial dimensions.

Acknowledgements

We gratefully acknowledge A. Alberti, M. Dalmonte, R. Gerritsma, S. Montangero, and M. Valiente for valuable feedback on the manuscript, A. B. Michelsen for performing the optimal control calculation, and J. M. Schurer for discussions (AN). This work was supported by The Hamburg Centre for Ultrafast Imaging (AN), the UPV/EHU grant EHUA15/17, the Spanish Ministerio de Economía y Competitividad through the project MINECO/FEDER FIS2015-69983-P, Basque Government IT986-16 (ER), and the Danish Council for Independent Research under the Sapere Aude program (ASD and NTZ).

Appendix A Analytical solutions to the Double Harmonic Potential

We consider a particle in a double harmonic potential as illustrated in Fig. 7 (thick black line), whose mathematical expression is given by:

[TABLE]

where $d_{L}$ and $d_{R}$ are the centres of the double harmonic potential one in left and the other in right of $x=0$, respectively. Here $\omega_{L}$ and $\omega_{R}$ are the trapping frequencies, $m_{\text{p}}$ is the mass of the trapped particle, and $\Delta$ is the height difference between the two harmonic potentials (i.e., the energy offset). Hereafter we shall assume that the potential is a continuous function. Thus, we must have $\frac{1}{2}m_{\text{p}}\omega_{L}^{2}d_{L}^{2}=\frac{1}{2}m_{\text{p}}\omega_{R}^{2}d_{R}^{2}+\Delta$.

Our goal is to solve the following eigenvalue problem

[TABLE]

For the sake of convenience, we introduce the length $\ell=\sqrt{\hbar/(m_{\text{p}}\omega_{L})}$ and energy $\epsilon=\hbar\omega_{L}$ scales. By transforming all the lengths with respect to $\ell$ and the energies with respect to $\epsilon$, and by rewriting the energy offset as $\Delta\equiv\hbar\omega_{L}\delta$, the continuity condition becomes $d_{L}^{2}=r^{2}\cdot d_{R}^{2}+2\delta$ where $r\equiv\omega_{R}/\omega_{L}$. In the rescaled units the Schrödinger equation for $x<0$ reads:

[TABLE]

By defining $E=\nu_{1}+\frac{1}{2}=r\cdot(\nu_{2}+\frac{1}{2})+\delta$, where $\nu_{1}$ and $\nu_{2}$ are real numbers (not necessary integers) to be determined, the above equation turns into:

[TABLE]

Hereafter, we introduce the new variable $y=-\sqrt{2}(x+d_{L})$, which turns the above equation into:

[TABLE]

However, note that the introduced variable, $y$, could also be defined as $y=\sqrt{2}(x+d_{L})$, which is a solution to the differential equation, too. Nevertheless, only $y=-\sqrt{2}(x+d_{L})$ gives a normalised and convergent wave function within the interval $(-\infty,0)$, whereas the former does not. We note that Eq. (21) has the form of a parabolic cylinder differential equation

[TABLE]

whose solutions are named parabolic cylinder functions and denoted as $D_{-a-\frac{1}{2}}(z)$ Abramowitz and Stegun (1992). In order to solve our eigenvalue problem, we shall make use of the following property of such functions

[TABLE]

which enable us to obtain the following solutions for the differential equation (21):

[TABLE]

Here $c_{1}$ and $c_{2}$ are some constants, to be determined later in order to normalise the wave function. The derivative of such wave function yields

[TABLE]

where $z_{1}\equiv-\sqrt{2}(x+d_{L})$ and $z_{2}\equiv\sqrt{2r}(x-d_{R})$. In order to find the unknowns, $c_{1}$, $c_{2}$, $\nu_{1}$ and $\nu_{2}$, one must impose the following conditions:

$\psi(x)|_{x=0}$ must be continuous. 2. 2.

$\psi(x)^{\prime}|_{x=0}$ must be continuous. 3. 3.

Energy is the same: $\nu_{1}+1/2=r(\nu_{2}+1/2)+\delta$. 4. 4.

Normalisation: $1=\int|\psi(x)|^{2}dx$.

If the potentials are symmetric around $x=0$, one can use the fact that the even solutions must have $\psi(x)^{\prime}|_{x=0}=0$ and the odd solutions $\psi(x)|_{x=0}=0$. In general, however, one must derive additional (transcendental) equations in order to determine the unknowns. Special solutions are exemplarily displayed in Fig. 7.

When the particle in the double well interacts with the (static) impurity located in $x=0$, we have to incorporate in the Scrödinger equation (18) the following additional (external) potential

[TABLE]

where

[TABLE]

Here $g_{\mathrm{1D}}^{e,s}=-\hbar^{2}/\mu a_{\mathrm{1D}}^{e,s}$ and $g_{\mathrm{1D}}^{o,s}=-\hbar^{2}a_{\mathrm{1D}}^{o,s}/\mu$ with $\mu$ being the reduced particle-impurity mass, and $a_{\mathrm{1D}}^{e,s},\,a_{\mathrm{1D}}^{o,s}$ are the 1D spin-dependent scattering lengths for even- and odd-waves, respectively. Besides, the action of the two operators appearing in Eq. (26) on a wave-function $\psi(x)$ is given by:

[TABLE]

where $\psi(0^{\pm})=\lim_{x\rightarrow 0^{\pm}}\psi(x)$. With these conditions we end up with:

[TABLE]

for even and odd scattering, respectively.

By using these conditions instead of the above outlined conditions 1. and 2. of the non-interacting particle problem, we end up with the following transcendental equation:

[TABLE]

where $z_{1}=-\sqrt{2}d_{L}$ and $z_{2}=-\sqrt{2r}d_{R}$, and from condition 3. we arrive at $\nu_{1}=r(\nu_{2}+1/2)+\delta-1/2$. Hence, we are basically left with one equation with one unknown, $\nu_{2}$. Therefore we can define the left-hand-side of Eq. (30) as being a function of $\nu_{2}$, $f(\nu_{2})$, which we have to solve as $f(\nu_{2})=0$ as a final step. This is done by finding the zeros of such a function numerically. Once the sequence of real numbers $\nu_{2}$ are found, the corresponding $\nu_{1}$ are obtained accordingly. Finally, by using the normalisation condition 4., we can obtain the constants $c_{1}$ and $c_{2}$.

A.1 Limiting cases

Let us note that in the case for which $r=1$, $\delta=0$, $d_{L}=d_{R}=0$, we have $\nu_{1}=\nu_{2}=\nu$. Moreover, if $a_{\text{1D}}^{o}=0$ and $-1/a_{\text{1D}}^{e}\equiv g^{e}_{\text{1D}}$, we obtain the non-trivial result calculated for two interacting atoms in a harmonic trap Busch et al. (1998):

[TABLE]

In the other limit, that is when $a_{\text{1D}}^{o}\equiv-g^{o}_{\text{1D}}$ and $1/a_{\text{1D}}^{e}=0$, we have

[TABLE]

Note that in both limits $E=\nu+1/2$ holds.

Appendix B Bosonic Josephson junction in the large particle limit

In this section we would like to show in a more rigorous way that, within the two-mode approximation, when the particle number per well is sufficiently large, the ensemble of bosons can be treated as a Bose-Einstein condensate, and therefore the Hamiltonian term

[TABLE]

can be replaced by Eq. (3). To this end, let us consider $N$ particles in both wells with a well-defined phase difference $\theta=\theta_{L}-\theta_{R}$, with $N\gg\Delta N$, i.e., the total number of particles is much greater than any variation due to the dynamics. In this scenario, the initial state of the compound bosonic ensemble and impurity spin system is given by

[TABLE]

where $|N,\theta\rangle$ denotes the state of the bosonic ensemble with $N$ particles and a given phase difference $\theta$, $|\text{vac}\rangle$ the vacuum state, and $|\sigma\rangle$ the spin state of the impurity. The distribution $\frac{1}{n!\left(N-n\right)!}$ is exponentially picked at $n=N/2$. When $N\gg 1$, the state characterises two condensates of $\frac{N}{2}$ particles in each well, i.e.,

[TABLE]

By exploiting the identity $[\hat{b},(\hat{b}^{\dagger})^{n}]=n(\hat{b}^{\dagger})^{n-1}$, the action of the Hamiltonian $\hat{H}_{\text{2m}}$ on this state is given by

[TABLE]

Hence, this shows that the state $|N,\theta\rangle\otimes|\sigma\rangle$ is almost an eigenstate of the Hamiltonian $\hat{H}_{\text{2m}}$. In particular, because of the last equality in the equation above, the operators $b^{\dagger}_{L,R}$ become proportional to $\sqrt{N}e^{i\theta_{L,R}}$ in the large-$N$ limit, thus demonstrating that in such a case the effective Hamiltonian is given by Eq. (3).

Appendix C Effective equation of motion for a condensate coupled to a spin impurity

In this section we derive the equations of motion for a bosonic ensemble in interaction with a quantum spin, thus without motional degrees of freedom for the latter. In particular, we focus on the simplest scenario for which the bosonic ensemble is in the condensate state (i.e., mean field theory), but we take into account all correlations between the ensemble and the spin system.

C.1 Many-body Hamiltonian and effective equations of motion

The Hamiltonian describing $N$ interaction bosons of mass $m$ with a spin-$\frac{1}{2}$ impurity is given by

[TABLE]

with $\hat{p}_{j}=-i\hbar\frac{\partial}{\partial x_{j}}$ being the momentum operator of the $j$-th boson, $V_{\text{ext}}(x_{j})$ the external potential (e.g., the double well), $U_{\alpha}(x_{j})$ a potential that depends on the internal state $|\alpha\rangle$ of the impurity (in our setting the contact potentials), $\Omega_{R}$ the Rabi frequency, and $\hat{\sigma}_{x}=|0\rangle\langle 1|+|1\rangle\langle 0|$ the Pauli matrix. The coupling constant $g$ describes the interaction among the bosons in quasi-1D, which typically depends on the three-dimensional s-wave scattering length and the transverse trap frequency. Note that if $U_{\alpha}(x_{j})=0$ $\forall\alpha$, then no coupling between the spin impurity and the ensemble of bosons exists. This implies that the two subsystems are fully uncorrelated and that the full state is simply given by the tensor product of a state living in the Hilbert space associated to the ensemble of bosons and a state of the impurity in $\mathbb{C}^{2}$.

C.2 Many-body state ansatz

The most general state for such an ensemble of $N$ bosons coupled to a two-level quantum system is given by

[TABLE]

with $|c_{0}|^{2}+|c_{1}|^{2}=1$. Such a state encompasses correlations among the bosons as well as between the bosonic ensemble and the spin impurity. Next, we assume that the $N$-particle state is a product state for each spin component (mean-field approximation), namely $|\Phi_{\alpha}\rangle=\prod_{j=1}^{N}|\varphi_{\alpha}^{(j)}\rangle$. Thus, the above outlined ansatz reduces to

[TABLE]

with $\langle\varphi_{\alpha}|\varphi_{\alpha}\rangle=1$ $\forall\alpha$.

C.3 Condensate-impurity equations of motion:

We aim at determining the equations of motion for the coefficients $c_{\alpha}$ as well as for the condensate wave functions $\varphi_{\alpha}(x)$. To this end, we employ the so-called Dirac-Frankel variational principle, $\langle\delta\Psi|(i\hbar\partial_{t}-\hat{H})|\Psi\rangle=0$, where $\delta\Psi$ denotes the variation of the many-body wave function with respect to the free parameters $c_{\alpha}$ and $\varphi_{\alpha}$.

After a lengthy calculation we arrive at the equations of motion for the spin-impurity coefficients ($\overline{\alpha}=1-\alpha$)

[TABLE]

with $\hat{H}_{0}^{\alpha}=\frac{\hat{p}^{2}}{2m}+V_{\text{ext}}+U_{\alpha}$, whereas for the condensate wave functions we obtain

[TABLE]

with $\langle\varphi_{\alpha}|\delta(x-y)|\varphi_{\alpha}\rangle=|\varphi_{\alpha}(x)|^{2}$, and

[TABLE]

Now, by substituting $\dot{c}_{\alpha}$ in Eq. (C.3) with Eq. (C.3) we finally arrive at

[TABLE]

where $\hat{H}_{\text{gp}}^{\alpha}[\varphi_{\alpha}]=\frac{\hat{p}^{2}}{2m}+V_{\text{ext}}+U_{\alpha}+g(N-1)|\varphi_{\alpha}|^{2}$ and $\langle\hat{H}_{\text{gp}}^{\alpha}[\varphi_{\alpha}]\rangle=\langle\varphi_{\alpha}|\hat{H}_{\text{gp}}[\varphi_{\alpha}]|\varphi_{\alpha}\rangle$. With these definitions, we can rewrite Eq. (C.3) as

[TABLE]

To further simplify the above outlined equations of motion, we first perform a unitary transformation on the coefficients, that is, $\mathbf{c}=\hat{U}\mathbf{C}$ with $\mathbf{c}\equiv(c_{0},c_{1})^{T}$ (similarly for $\mathbf{C}$), where

[TABLE]

By applying this unitary, we get

[TABLE]

Now, let us note that $\langle\hat{H}_{\text{gp}}^{\alpha}[\varphi_{\alpha}/\sqrt{2}]\rangle$ is a real number, as it corresponds to the expectation value of a (nonlinear) hermitian operator, but that also $i\hbar\langle\varphi_{\alpha}|\dot{\varphi}_{\alpha}\rangle$ is a real number. To see this, let us note that $\frac{d}{dt}\langle\varphi_{\alpha}|\varphi_{\alpha}\rangle=2\Re\{\langle\varphi_{\alpha}|\dot{\varphi}_{\alpha}\rangle\}=0$. This implies that $\langle\varphi_{\alpha}|\dot{\varphi}_{\alpha}\rangle$ is purely imaginary and as a consequence $i\hbar\langle\varphi_{\alpha}|\dot{\varphi}_{\alpha}\rangle$ is purely real. Hence, $\eta_{\alpha}(t)$ is a good phase definition. Thus, by defining $|\phi_{\alpha}\rangle:=|\varphi_{\alpha}e^{-\frac{i}{\hbar}f_{\alpha}(t)}\rangle$, we have

[TABLE]

Now, since $|\varphi_{\alpha}\rangle=e^{\frac{i}{\hbar}f_{\alpha}(t)}|\phi_{\alpha}\rangle$, we can use Eq. (C.3) in order to get a differential equation for $|\phi_{\alpha}\rangle$. After a straightforward calculation, one arrives at

[TABLE]

Alternatively, we can combine Eq. (49) and Eq. (C.3) in a single one by defining $|\psi_{\alpha}\rangle:=C_{\alpha}|\phi_{\alpha}\rangle$. This yields

[TABLE]

with $|\phi_{\alpha}\rangle=|\psi_{\alpha}\rangle/\sqrt{\langle\psi_{\alpha}|\psi_{\alpha}\rangle}$. In this way, one can simulate a single system of coupled differential equations (C.3), and therefore forget about the equation of motion of the coefficients $C_{\alpha}$. Note that $\sum_{\alpha=0}^{1}\int_{\mathbb{R}}dx|\psi_{\alpha}(x,t)|^{2}=1$, and thereby $\int_{\mathbb{R}}dx|\psi_{\alpha}(x,t)|^{2}\neq 1$.

Appendix D Correlated hopping in the many-body mean-field scenario

We have investigated numerically the possibility to realise the correlated hopping, as for the single particle case discussed in Sec. IV.4, when the system is an ensemble of ultra-cold bosons described in the mean-field approximation, as discussed in the previous appendix. To this end, we solved the equations of motion (C.3). In Ref. Brouzos et al. (2015) it has been shown that for a bosonic Josephson junction without contact to a spin impurity, whereby the condensate is prepared initially in one of the two wells, the compensating control pulse technique allows to transfer the entire condensate to the other well. Such a pulse, that represents the energy offset between the two wells, has a simple analytical form: $\Delta(t)=-\Delta U-U\cos(Jt)$. Here $J=(E_{R}^{0}-E_{L}^{0})/\hbar$, $\Delta U=U_{L}-U_{R}$, $U=(U_{L}+U_{R})/2$,

[TABLE]

and similarly for $U_{R}$. Inspired by this work, we modulated the Rabi frequency by using the following compensating control pulse:

[TABLE]

Contrarily to the case without spin impurity, we need to modulate the Rabi frequency over time in order to fulfil the resonance condition that, because of the inter-particle interaction, changes in time as the condensate tunnels through the barrier. The result of such a transfer process is illustrated in Fig. 8. As it is shown, the condensate transfer does not work very well (roughly 50% only) and, importantly, the correlated hopping is not occurring at all, as other (undesired) states are strongly populated.

We then used the pulse (53) as initial guess for the Rabi frequency for an optimal control study and at the same time we allowed control for the energy offset $\Delta$. The outcome of such an investigation is illustrated in Fig. 8. As it can be seen, even by varying two experimentally accessible parameters via optimal control theory we were unable to realise the desired process. Even though this is not a rigorous mathematical proof, it strongly indicates that the correlated hopping in presence of interactions is substantially modified. Moreover, the process seems not attainable within a mean-field description. This might be also due to the fact that our derivation of the target gauge invariant Hamiltonian relies on the single particle tunnelling and that in the many-body scenario high-order tunnelling effects might play a crucial role, while in our description such effects are not allowed. Indeed, as it has been observed for a purely bosonic Josephson junction without spin impurity, at times close to the so-called quantum speed limit the mean-field description breaks down Brouzos et al. (2015). Hence, it would be interesting to understand, whether in the presence of a gauge field, the process can take place by means of a concerted action of inter- and intra-particle correlations, a task that we shall pursue in the future.