Impurity and soliton dynamics in a Fermi gas with nearest-neighbor interactions

TL;DR

This paper investigates how impurities and local quenches affect soliton-antisoliton pair formation and dynamics in a one-dimensional Fermi gas with nearest-neighbor interactions, combining numerical and analytical approaches.

Contribution

It introduces a combined numerical and analytical study of impurity-induced soliton dynamics in a Fermi gas with nearest-neighbor interactions, revealing novel confinement phenomena.

Findings

Impurities induce soliton-antisoliton pairs in the Fermi gas.

Antisolitons can become confined near their origin due to existing solitons.

The dynamics differ significantly depending on the initial soliton presence.

Abstract

We study spinless fermions with repulsive nearest-neighbor interactions perturbed by an impurity particle or a local potential quench. Using the numerical time-evolving block decimation method and a simplified analytic model, we show that the pertubations create a soliton-antisoliton pair. If solitons are already present in the bath, the two excitations have a drastically different dynamics: The antisoliton does not annihilate with the solitons and is therefore confined close to its origin while the soliton excitation propagates. We discuss the consequences for experiments with ultracold gases.

Tables

Table 1. Table 1: The number of antisolitons in the initial state, the difference Δ ​ N AS = N AS ​ ( t = 8 / J ) − N AS ​ ( 0 ) Δ subscript 𝑁 AS subscript 𝑁 AS 𝑡 8 𝐽 subscript 𝑁 AS 0 \Delta N_{\text{AS}}=N_{\text{AS}}(t=8/J)-N_{\text{AS}}(0) , and the initial occupation probability of site j 0 subscript 𝑗 0 j_{0} . The numbers are shown for the MI and incommensurate states and interactions U = V = 10 ​ J 𝑈 𝑉 10 𝐽 U=V=10J and U = V = 50 ​ J 𝑈 𝑉 50 𝐽 U=V=50J . The number at time t = 8 / J 𝑡 8 𝐽 t=8/J corresponds to the time evolution with the local potential quench. For the impurity case, the numbers are the same within the accuracy shown here.

	$U/J=V/J$	$N_{\text{AS}}\left(0\right)$	$\Delta N_{\text{AS}}$	$⟨n_{j_0}^b\left(0\right)⟩$
MI	50	0.02	1	1
	10	0.6	0.9	1
incomm.	50	0.02	0.8	0.8
	10	0.5	0.7	0.8

Equations

$$H_{\text{b}}=-J\sum _{⟨i,j⟩}{b}_i^{†}{b}_j+V\sum _j\left(n_j^b-\frac{1}{2}\right)\left(n_{j+1}^b-\frac{1}{2}\right),$$

$$H^{\prime }$$

$$⟨n_j^b⟩=\frac{1}{2}\sum _{l_1,\cdots ,l_N}|{\phi }_{k_1,\cdots ,k_N}^{l_1,\cdots ,l_N}{|}^2\times \prod _{d=1}^{j-1}\left[2\sum _{\alpha =1}^N{\delta }_{d,l_{\alpha }}-1\right]+\frac{1}{2}.$$

$$|\psi (0)⟩=$$

$$\to$$

$$H_{\text{b}}=J_{xy}\left(\frac{1}{2}\sum _{⟨i,j⟩}{S}_i^{+}{S}_j^{-}+\Delta \sum _j{S}_j^z{S}_{j+1}^z\right)$$

$$\begin{array}{rl}{J}_{xy}& =2J,\\ \Delta & =\frac{V}{J_{xy}}.\end{array}$$

$$H_b=\frac{J_{xy}}{4}\sum _{j=1}^{L-1}\left({\sigma }_j^x{\sigma }_{j+1}^x+{\sigma }_j^y{\sigma }_{j+1}^y+\Delta {\sigma }_j^z{\sigma }_{j+1}^z\right).$$

$$n_j^b=S_j^z+\frac{1}{2}=\frac{1}{2}{\sigma }_j^z+\frac{1}{2}.$$

$$H^{\prime }=\frac{U}{2}{\sigma }_{j_0}^z,$$

$$\begin{array}{rl}{\tau }_{j+\frac{1}{2}}^z& ={\sigma }_j^z{\sigma }_{j+1}^z,\\ {\tau }_{j+\frac{1}{2}}^y& =\prod _{i=1}^j{\sigma }_i^x.\end{array}$$

$$H_b=\frac{J_{xy}}{4}\sum _{j=2}^{L-2}\left(1-{\tau }_{j+\frac{1}{2}}^z\right){\tau }_{j-\frac{1}{2}}^y{\tau }_{j+\frac{3}{2}}^y+\frac{J_{xy}\Delta }{4}\sum _{j=1}^{L-1}{\tau }_{j+\frac{1}{2}}^z$$

$${\sigma }_{j_0}^z={\sigma }_1^z\prod _{j=1}^{j_0-1}{\tau }_{j+\frac{1}{2}}^z.$$

$${\tau }_l^z$$

$${\tau }_l^{+}$$

$$\begin{array}{rl}{H}_b& =\frac{J_{xy}}{2}\sum _{l=2}^{L-2}[\left(1-c_l^{†}{c}_l\right)c_{l-1}^{†}{c}_{l+1}\\ & -\left(1-c_l^{†}{c}_l\right)c_{l-1}^{†}{c}_{l+1}^{†}+\text{H.c.}]+\frac{J_{xy}\Delta }{2}\sum _{l=1}^{L-1}{c}_l^{†}{c}_l\end{array}$$

$${\sigma }_{j_0}^z=\prod _{l=1}^{j_0-1}\left(2c_l^{†}{c}_l-1\right).$$

$$\left(1-c_l^{†}{c}_l\right)\approx 1$$

$$H_b\approx \frac{J_{xy}}{2}\sum _{l=2}^{L-2}\left(c_{l-1}^{†}{c}_{l+1}+\text{H.c.}\right).$$

$${\phi }_k^l=\sqrt{\frac{2}{L}}\sin (kl),$$

$$c_l^{†}=\sum _k{\phi }_k^{l\ast }{c}_k^{†},$$

$$H_b=\sum _k{ϵ}_k{c}_k^{†}{c}_k,$$

$$c_l^{†}\to e^{-i\frac{\pi }{2}l}{c}_l^{†},$$

$${\phi }_k^l=e^{-i\frac{\pi }{2}l}\sqrt{\frac{2}{L}}\sin (kl).$$

$${ϵ}_k=-J_{xy}\cos (2k).$$

$$\begin{array}{r}|{\Psi }_{1,2,\cdots ,N}⟩=\sum _{l_1,\cdots ,l_N}{\phi }_{k_1,\cdots ,k_N}^{l_1,\cdots ,l_N}|l_1⟩\otimes |l_2⟩\otimes \cdots |l_N⟩.\end{array}$$

$$|l⟩=|0,\cdots ,0,1_l,0,\cdots ,0⟩.$$

$${\phi }_{k_1,\cdots ,k_N}^{l_1,\cdots ,l_N}=\frac{1}{\sqrt{N!}}\left|\begin{array}{cccc}{\phi }_{k_1}^{l_1}& {\phi }_{k_2}^{l_1}& \cdots & {\phi }_{k_N}^{l_1}\\ {\phi }_{k_1}^{l_2}& {\phi }_{k_2}^{l_2}& \cdots & {\phi }_{k_N}^{l_2}\\ ⋮& \ddots & & ⋮\\ {\phi }_{k_1}^{l_N}& \cdots & & {\phi }_{k_N}^{l_N}\end{array}\right|.$$

$$\begin{array}{rl}& ⟨{\Psi }_{1,\cdots ,N}|S_j^z|{\Psi }_{1,\cdots ,N}⟩\\ & =\frac{1}{2}⟨{\Psi }_{1,\cdots ,N}|\prod _{d=1}^{j-1}\left(2c_d^{†}{c}_d-1\right)|{\Psi }_{1,\cdots ,N}⟩\\ & =\frac{1}{2}\sum _{l_1,\cdots ,l_N}|{\phi }_{k_1,\cdots ,k_N}^{l_1,\cdots ,l_N}{|}^2\prod _{d=1}^{j-1}\left[2\sum _{\alpha =1}^N{\delta }_{d,l_{\alpha }}-1\right],\end{array}$$

$$⟨n_j^b⟩=⟨S_j^z⟩+\frac{1}{2}.$$

Full text

I Introduction

The properties of an impurity coupled to a bath are a paradigmatic problem of many body physics. For classical baths, the diffusion of massive particles is described by Brownian motion. A quantum bath leads to more complex physics. In two and three dimensions, the problem is described in terms of quasiparticles in which the impurity is surrounded by excitations of the bath. This is the case for polarons Mahan_Many-particle_2000 ; Grusdt_New_2016 , where the bath is made of phonons, bosonic particles, or fermions such as in a Fermi liquid. Novel effects exists in one dimension where the excitations of the bath are drastically affected by interactions. In particular, if the one-dimensional quantum bath is a Tomonaga-Luttinger liquid Giamarchi , the massless excitations of the bath lead to new diffusion properties for the impurity, such as a subdiffusion Zvonarev_Spin_2007 ; Lamacraft_Dispersion_2009 , and behavior ranging from polaronic to Anderson’s orthogonality catastrophe Kantian .

Cold atomic gases provide an ideal testing ground for such nonequilibrium quantum many-body phenomena due to their weak coupling to the environment and tunable parameters. The motion of initially localized impurities of a different spin state Kohl and atom species Catani was measured in one-dimensional tubes in the continuum. Recent developments in experimental setups allow addressing and imaging atoms in optical lattices at the resolution of a single lattice site. This has made it possible to measure the dynamics of specific nonequilibrium many-body states Hung_In_2014 ; Ott_Single_2016 ; Kuhr_Quantum-gas_2016 . One can for instance create a local energy shift by focusing a laser beam on a selected site and flip the spin of the atom at that site by a microwave pulse. These techniques were used for recording the time evolution of initially localized spin impurities Fukuhara_Quantum_2013 and magnons Fukuhara_Microscopic_2013 in bosonic rubidium. The quantum gas microscope technique has also been extended to fermionic atoms Cheuk_Quantum-Gas_2015 ; Parsons_Site-Resolved_2015 ; haller_single-atom_2015 ; Edge_Imaging_2015 ; Boll_Spin_2016 ; Drewes_Antiferromagnetic_2016 .

On the theoretical front, mobile impurities in homogeneous baths have been largely studied Jaksch ; Kamenev ; Demler ; Massel ; Kantian ; Grusdt ; Sirker . The situation is more complicated if the bath has a structure. In particular, long-range interactions can lead to a periodic arrangement of the bath particles. New behaviors are possible in such systems, such as the localization of the impurity Horovitz_Phase_2013 . For a fermionic bath with nearest-neighbor interactions in the Mott insulator (MI) state, a diffusive motion of the impurity was predicted to be connected to soliton excitations in the bath Visuri2 . This raises the question of how the impurity and the excitations behave when the ground state of the bath already contains solitons. Soliton excitations occur in a variety of one- and quasi-one-dimensional compounds, such as quantum spin chains Nagler_Propagating_1982 ; mourigal_fractional_2013 ; braun_emergence_2005 ; Umegaki_Spinon_2015 . Quantum gas experiments provide complementary systems and offer the possibility to measure spatially the motion of such excitations and the impurity. These phenomena can thus be relevant in experiments with polar KRb molecules Yan_Observation_2013 ; Hazzard_Many-body_2014 , Rydberg atoms Barredo_Coherent_2015 ; labuhn_tunable_2016 , and magnetic atoms dePaz_Nonequilibrium_2013 ; Baier_Extended_2016 in optical lattices, as well as ions in rf traps Neyenhuis_Observation_2016 , where effects of long-range interactions have been observed.

In this article, we use a combination of numerical and analytic techniques to study the dynamics of an impurity coupled to a bath of fermions with repulsive nearest-neighbor interactions. We compare it to the simpler case of a local quench in potential. Although both perturbations create a soliton and an antisoliton excitation, we show that the antisoliton can form a bound state with the impurity, which drastically affects its motion. We consider baths in the MI state and at an incommensurate filling for which solitons are present in the ground state. In the latter case, energy conservation prevents the antisoliton excitation from annihilating with the solitons and further constrains its motion. We discuss the possibility of observing these phenomena in experiments with ultracold gases.

II The model

We consider a bath described by

[TABLE]

where $b_{j}$ ($b_{j}^{\dagger}$) annihilates (creates) a bath fermion at site $j$, $n_{j}^{b}=b_{j}^{\dagger}b_{j}$ is the number operator of the bath fermions, $J$ is the tunneling amplitude, and $V$ is the nearest-neighbor interaction energy. We consider baths both in the MI state at commensurate (half) filling and at an incommensurate filling slightly above one half. The initial state is shown in the schematic drawing of Fig. 1(a).

Our results also apply for the XXZ spin model, to which model (1) can be mapped Jordan_Uber_1928 . The dynamics of excitations in the XXZ model have been studied by the Bethe Ansatz Vlijm_Spinon_2016 , time-evolving block decimation (TEBD) Ganahl_Observation_2012 ; Visuri2 , and time-dependent DMRG Langer_Real-time_2009 ; Langer_Real-time_2011 methods. The previous studies considered antiferromagnetic Vlijm_Spinon_2016 ; Visuri2 or fully polarized Ganahl_Observation_2012 initial states with one or a few flipped spins at the center. Ballistic and diffusive transport regimes were studied using different types of quenches Langer_Real-time_2009 ; Langer_Real-time_2011 .

Here, we study the dynamics after either an impurity particle or a static potential barrier is introduced at the center site $j_{0}$ at time $t=0$. The time evolution is described by the Hamiltonian $H~{}=~{}H_{\text{b}}+H^{\prime}$, where

[TABLE]

respectively for the impurity and the potential. The repulsive on-site interaction between the impurity and the bath fermions is denoted by $U>0$, the annihilation operator of the impurity by $c_{j}$, and $n_{j}=c_{j}^{\dagger}c_{j}$. The tunneling energies of the impurity and the bath fermions are equal. For the static potential barrier, the barrier height is equal to the on-site interaction energy.

III Results and discussion

We compute both ground-state and time-dependent observables with the numerical TEBD method Daley_Time-dependent_2004 ; Vidal_Efficient_2004 . An odd number of lattice sites $L=61$ and open boundary conditions make the ground state nondegenerate for commensurate filling with $N_{b}=31$ particles. For the incommensurate filling, $N_{b}=33$. A Schmidt number $\chi=120$ is used in the imaginary time evolution. In the real time evolution, we use $\chi=160$ and the second-order Trotter decomposition with time step $\delta t=0.01/J$. The interaction energies are $U=V=50J$ unless mentioned otherwise. We focus on strong interactions to suppress the creation of solitons and antisolitons due to quantum fluctuations and to distinguish clearly the effects of the quenches. We show that essentially the same phenomena are observed also for interactions $U=V=10J$, which are closer to experimentally realizable values.

III.1 Ground state density distribution

In the MI state, for $V\gg J$, there is a particle every two sites and the density oscillates between [math] and $1$. For incommensurate filling, the amplitude of the density oscillation decreases away from the boundaries. A low density of excess particles or holes with respect to half filling leads to soliton or antisoliton excitations Giamarchi where the phase of the density oscillation changes by $\pi$. The maxima of the soliton distribution due to the excess particles can be seen as domain walls (DWs) – pairs of neighboring maxima or minima – in the density distribution of Fig. 1(b). For open boundary conditions, the solitons are not completely delocalized weiss_finite_2008 as they would be for periodic boundary conditions.

When the density of solitons is low and $V\gg J$, the solitons behave as free fermions Giamarchi . In this case, the density profile can be accurately predicted by a DW model, where one uses a bond representation to map the solitons to free fermions in an otherwise empty lattice (see Appendix A). The wave function of $N$ free fermions can be written as $\ket{\Psi_{1,2,\cdots,N}}=\sum_{l_{1},\cdots,l_{N}}\varphi_{k_{1},k_{2},\cdots,k_{N}}^{l_{1},l_{2},\cdots,l_{N}}\times\ket{l_{1}}\otimes\ket{l_{2}}\otimes\cdots\ket{l_{N}}$, where $l_{\alpha}$ is the coordinate of fermion $\alpha$. The single-particle state where site $l$ is occupied and the other sites are empty is denoted by $\ket{l}=\ket{0,\cdots,0,1_{l},0,\cdots,0}$. The coefficient $\varphi_{k_{1},\cdots,k_{N}}^{l_{1},\cdots,l_{N}}$ is given by the Slater determinant formed of the single-particle wave functions $\varphi_{k_{m}}^{l}=\sqrt{\frac{2}{L}}\sin(k_{m}l)$, where $k_{m}=\frac{m\pi}{L}$ and $m=1,\cdots,N$. One can thus calculate the expectation value of the density $\langle n_{j}^{b}\rangle=\bra{\Psi_{1,\cdots,N}}n_{j}^{b}\ket{\Psi_{1,\cdots,N}}$ as

[TABLE]

More details are given in Appendix A. The result of the DW model shown in Fig. 1(b) agrees very well with the numerical solution.

III.2 Time evolution in the Mott insulator state

We now turn to the time evolution. In the MI state, at $t=0$, the central site $j_{0}$ is occupied by a bath fermion with a high probability. A local potential quench [Eq. (3)] causes the fermion at $j_{0}$ to tunnel to the neighboring site. A soliton and an antisoliton excitation are created in this process and propagate symmetrically in opposite directions. The density profiles in Fig. 2

show an excitation as an additional DW. We can verify that the two excitations propagate symmetrically by studying the correlation of density on both sides of $j_{0}$, as is done in Appendix B. The symmetric propagation results in the complete inversion of the density profile between the excitations, as shown in Fig. 2(a). The analytic result in Fig. 2(a) is based on the DW model and agrees very well with the numerical solution.

In the case of an impurity, the time evolution is richer. The impurity can oscillate between the two empty sites which form the antisoliton, and thus form a bound state with the antisoliton Visuri2 . In this case, only the soliton excitation propagates and the density is inverted on one side of the lattice. The superposition of evolutions where both excitations propagate or only the soliton propagates in either direction results in a reduced amplitude of the density oscillation with respect to the initial state. In addition, the impurity can move past a neighboring occupied site only in a second-order process with velocity $\frac{4J^{2}}{U}$ Georges_Strongly_2013 , which for $U\gg J$ results in the very slow motion of the impurity observed in Fig. 2. For the incommensurate filling, the numerical results in Fig. 2(b) show a similar population reversal in the regions between the DWs as in the MI. The DWs in the density remain in the final state after the excitation has passed. Using the DW model when solitons are present in the ground state is more delicate since the excitations can interact with these solitons (see Appendix A) and a mapping to free fermions would be inaccurate.

In order to understand how the solitons in the initial state affect the new soliton and antisoliton excitations, we compute the distributions of neighboring filled sites $\langle n_{j}^{b}n_{j+1}^{b}\rangle$ and empty sites or holes $\langle n_{j}^{h}n_{j+1}^{h}\rangle=\langle(1-n_{j}^{b})(1-n_{j+1}^{b})\rangle$ as functions of time. They correspond respectively to the distributions of solitons and antisolitons. Note that the measurement of such correlations is well within reach of experiments with quantum gas microscopes Preiss_Strongly_2015 . The soliton distribution for the MI initial state is shown in panels (a) and (b) of Fig. 3. The distribution is zero at $j_{0}-1$ and $j_{0}$, and obtains maxima at $j_{0}-2$ and $j_{0}+1$ when the bath fermion at $j_{0}$ tunnels either to $j_{0}-1$ or $j_{0}+1$. The distributions result from the superposition of these two configurations. Correspondingly, the antisoliton distributions in panels (c)–(f) have initially a maximum at $j_{0}-1$ and $j_{0}$.

The different perturbations lead to distinctly different dynamics. In Figs. 3(a) and (c), where the potential is quenched at $j_{0}$ in the MI state, the distributions of both solitons and antisolitons become zero at the center, indicating that the two excitations indeed propagate symmetrically in opposite directions. The soliton distributions in Figs. 3(a) and (b) are almost identical for the static barrier and the impurity. For the antisolitons, on the contrary, creating an impurity at $j_{0}$ produces an interference pattern and a maximum remains at the center, as seen in Fig. 3(d). The interference pattern indicates that the antisoliton can be in several momentum states. Whereas in the potential quench, the soliton and antisoliton excitations have a high probability to obtain momenta of equal magnitude in opposite directions, in the case of the impurity, part of the momentum in the direction of the antisoliton can be absorbed by the impurity. The antisoliton can thus have a smaller momentum than the soliton. The impurity can form a bound state with the antisoliton, as shown by the remaining maximum at the center of the lattice. These differences are also seen in the line profiles of Fig. 4(c) taken at time $t=6/J$. The soliton distributions at $t=6/J$ are very similar to the antisoliton distribution for the static barrier in Fig. 4(c).

III.3 Time evolution at incommensurate filling

For incommensurate filling, an additional effect is observed. Figure 4(b) shows that the soliton excitation does not propagate as a free particle but is slowed down by the interactions with the other solitons in the system. Interestingly, the antisoliton excitation moves even less. After the initial propagation, it stays confined in the central region of the lattice, as shown in Figs. 3(e) and (f) and Fig. 4(d). This confinement can be explained by energy conservation: the motion of an antisoliton excitation past a soliton would mean that they annihilate. The time evolution can be drawn schematically as

[TABLE]

where the sites with (without) bath particles are denoted by o (x). The impurity is not drawn but is considered to stay at the underlined central site. When $U=V\gg J$ and kinetic energy is not taken into account, the total energy in the first configuration is $E=2V+U=3V$. On the second line, $E=3V$, whereas on the last line, where the antisoliton has annihilated on the right side, $E=2V$. The annihilation should have a very low probability for $V\gg J$ since the released energy cannot be absorbed as kinetic energy. Only a soliton can move past another soliton and conserve energy. Since the solitons in the initial state are not completely localized, the antisoliton distribution does not go to zero abruptly in Fig. 3 but rather diminishes smoothly.

To verify that the antisoliton does not annihilate with the solitons, we compute the total number of solitons $N_{\text{S}}(t)=\sum_{j}\bra{\psi(t)}n_{j}^{b}n_{j+1}^{b}\ket{\psi(t)}$ and antisolitons $N_{\text{AS}}(t)=\sum_{j}\bra{\psi(t)}n_{j}^{h}n_{j+1}^{h}\ket{\psi(t)}$. Figure 4(a) shows $N_{\text{AS}}(t)$, which is close to zero at $t=0$ and increases until $t\approx 1/J$, the time scale for the tunneling of the bath particle away from $j_{0}$. The number saturates to a value close to $1$ in the MI state and to a smaller value at incommensurate filling. For $t\gtrsim 1/J$, we have checked that both $N_{\text{S}}(t)$ and $N_{\text{AS}}(t)$ stay constant, confirming that no annihilation takes place. The lower saturation value for incommensurate filling is due to the finite probability of site $j_{0}$ being initially empty, in which case excitations would not be created. Consistently, the value to which $N_{\text{AS}}(t)$ saturates is approximately the initial occupation probability $\langle n_{j_{0}}^{b}(0)\rangle$.

IV Experimental realization with ultracold dipolar gases

Lattice models with long-range interactions have so far been realized with bosonic dipolar atoms dePaz_Nonequilibrium_2013 ; Baier_Extended_2016 and molecules Yan_Observation_2013 ; Hazzard_Many-body_2014 . In one dimension, spinless fermions with nearest-neighbor interactions can be mapped to a spin model Jordan_Uber_1928 . Since spin systems can be mapped to hard-core bosons Giamarchi , bosonic particles with both hard-core and long-range interactions could be used to realize the fermionic bath studied here.

In quantum gas experiments, $\frac{V}{J}$ between zero and approximately $2$ was measured in extended Bose-Hubbard Baier_Extended_2016 and $t-J$-like dePaz_Nonequilibrium_2013 models realized with magnetic atoms. In optical lattices, $U$ is tunable by Feshbach resonances inouye_observation_1998 and $J$ by the lattice spacing and depth, which allows to tune $\frac{V}{J}$ Cugliandolo . For the same atoms and laser wavelengths as in Baier_Extended_2016 , we estimate that a larger lattice depth would allow to reach $V\approx 10J$, where $J\approx 2.7$ Hz Visuri2 . A coherent Bose-Einstein condensate was preserved for up to 1 s Baier_Extended_2016 , which gives a time scale $t\approx 2.7\frac{1}{J}$ sufficiently long to observe the different dynamics resulting from the two different perturbations. The time scale required to observe the confinement of the antisoliton at an incommensurate filling of the bath is of the same order of magnitude. Another possibility to realize the type of bath studied here are polar molecules confined to deep lattices. Spin exchange by dipole-dipole interaction was demonstrated with immobile KRb molecules using different rotational states as pseudospin states Yan_Observation_2013 ; Hazzard_Many-body_2014 .

The results of Sec. III have been obtained for a large repulsion for clarity. Figures 4(e) and (f) show that the same clear differences are observable for $U=V=10J$, although the probability of finding neighboring empty sites in the ground state is larger. The larger effective tunneling energy of the impurity leads to a broader maximum of the antisoliton distribution at the center of the lattice. Figures corresponding to Figs. 2–4(a) are presented in Appendix C for $U=V=10J$.

Excitations to higher bands can be prevented by making the band gap an order of magnitude larger than the tunneling energies and interactions dePaz_Nonequilibrium_2013 ; Baier_Extended_2016 . We therefore expect choosing interactions $U=V\approx 10J$ larger than the bandwidth of the lowest band and smaller than the band gap to be feasible so that the single-band approximation is valid. We consider a system initially in the ground state at zero temperature, whereas in experiments, the temperature of the gas cloud is nonzero. For temperatures close to or larger than $V$, thermally excited soliton-antisoliton pairs could have an effect on the dynamics. To minimize the effects of thermal excitations, post-selection techniques Fukuhara_Quantum_2013 could be used. Note that initial-state preparation can be used to create an excited initial state with alternating occupation Schreiber_observation_2015 .

In experiments with ultracold gases, uniform box potentials have recently been realized Meyrath_Bose-Einstein_2005 ; van_Es_Box_2010 ; Gaunt_Bose-Einstein_2013 ; Mukherjee_Homogeneous_2016 . Most experiments however use a harmonic potential for confining the gas cloud. We show in Appendix D that a Mott insulator ground state can be realized in the central region of a sufficiently shallow harmonic trap. The dynamics in this case is the same as in a box trap. To realize an incommensurate phase, a box trap is required. On the other hand, superimposing a harmonic potential and a box leads to a stronger confinement of the antisoliton than a uniform box potential, as shown in Appendix D.

V Conclusions

In summary, we observe that the presence of solitons in the bath leads to very different dynamics of the soliton and antisoliton excitations created by a local perturbation. The antisoliton does not annihilate with the solitons and is therefore confined close to its origin while the soliton excitation propagates. This is an example of the restrictions imposed by energy conservation on the dynamics: For interactions larger than the bandwidth, the energy released in an annihilation could not be absorbed as kinetic energy. Besides the numerical results, the simplified analytic model developed here offers a basis for understanding the dynamics of soliton excitations which occur in various physical systems.

Acknowledgements.

This work was supported by the Academy of Finland through its Centres of Excellence Programme (2012-2017) and under Project Nos. 263347, 251748, and 272490, and by the European Research Council (ERC-2013-AdG-340748-CODE). Computing resources were provided by CSC–the Finnish IT Centre for Science and the Aalto Science-IT Project. This work was supported in part by the Swiss NSF under Division II and by the ARO-MURI Non-equilibrium Many-body Dynamics grant (W911NF-14-1-0003).

Appendix A Simplified domain wall model

Using a dual bond representation, one can build a simplified model of the soliton excitations. Here, we consider the ground state for the commensurate and incommensurate filling with a low density of excess particles. In this case, the solitons can be assumed to not interact. A system with filling slightly below one half, where antisolitons exist in the ground state, could be treated in the same way. We only consider dynamics in the case of a local potential quench and the MI initial state, where the soliton or antisoliton excitation can be described in terms of a single particle. Extending the model to account for the motion of the impurity or dynamics in the case of solitons in the ground state would require taking the interactions into account.

A.1 Dual bond representation

The model of spinless fermions with nearest-neighbor interactions [Eq. (1) of the main text] can be mapped to the XXZ spin model

[TABLE]

by the Jordan-Wigner transformation Jordan_Uber_1928 . The mapping is illustrated in Fig. 5.

In the first sum of Eq. (4), $\langle i,j\rangle$ denotes neighboring sites. The coupling constant between the spins is denoted by $J_{xy}$ and the anisotropy in the $z$ direction by $\Delta$. The operators $S_{j}^{\pm}~{}=~{}S_{j}^{x}~{}\pm~{}iS_{j}^{y}$ are the raising and lowering operators on site $j$ of a spin-$\frac{1}{2}$ chain, $S_{j}^{a}=\frac{1}{2}\sigma_{j}^{a}$, and $\sigma_{j}^{a}$ are the Pauli matrices with $a=x,y,z$. The parameters in terms of $J$ and $V$ are

[TABLE]

The soliton excitations correspond in the XXZ model to two neighboring spins up, and antisolitons to two neighboring spins down. These excitations can represented by filled bonds, which are then mapped to spinless fermions in an otherwise empty lattice. One starts by replacing the lattice site operators in Hamiltonian (4) by bond operators. In terms of the Pauli matrices $\sigma^{a}$, $H_{b}$ is written as

[TABLE]

In the spin notation, density operator is written as

[TABLE]

The local potential term of Eq. (3) of the main text therefore transforms into

[TABLE]

leaving out the constant term. The matrices $\sigma^{a}$ can be transformed into bond operators by the Kramers-Wannier transformation Kramers_Statistics_1941 ; Essler_Duality_2009

[TABLE]

The indexing of the bonds is illustrated in Fig. 6. When $\sigma^{y}=-i\sigma^{z}\sigma^{x}$ is also transformed, the Hamiltonian (6) becomes

[TABLE]

and $\sigma_{j_{0}}^{z}$ in term (7) becomes

[TABLE]

for $j_{0}\geq 2$. We consider the case where the system is initially in the commensurate Mott insulator phase and $\Delta\gg 1$. When the number of lattice sites is odd and $\sum_{j}\langle S_{j}^{z}\rangle=\frac{1}{2}$, the spins at the edge sites are up and one can fix the boundary condition $\langle\sigma_{1}^{z}\rangle=1$.

If neighboring spins point in opposite directions, the bond operator $\tau^{z}$ gives the value $-1$, and for neighboring spins in the same direction, $+1$. The value $+1$ therefore corresponds to the existence of a domain wall, which can be a soliton (two neighboring spins up) or an antisoliton (two neighboring spins down). We will use here the term domain wall (DW) for both since they have the same bond representation. For the configuration of Fig. 6, the value of $\langle\sigma_{j_{0}}^{z}\rangle$ is determined by the parity of $j_{0}$. When $j_{0}$ is odd, $\langle\sigma_{j_{0}}^{z}\rangle=1$, and when $j_{0}$ is even, $\langle\sigma_{j_{0}}^{z}\rangle=-1$. We focus on the case where $j_{0}$ is odd and there is initially a spin up at $j_{0}$.

Changing into the bond indices $l=j-\frac{1}{2}$, the operators $\tau_{l}^{z}$ and $\tau_{l}^{y}=i\left(\tau_{l}^{-}-\tau_{l}^{+}\right)$ can be transformed into spinless fermion operators using

[TABLE]

In terms of the fermion operators, one can write Essler_Duality_2009

[TABLE]

and

[TABLE]

The creation of a fermion now corresponds to the creation of a DW in the original spin chain, as depicted in Fig. 7. The DW representation can be applied to calculating both ground state quantities and time-dependent ones, as is done in Secs. A.2 and A.3.

A.2 Density distribution in the ground state

The density operator can be written in terms of the $S_{j}^{z}$ spin operator as $n_{j}^{b}=S_{j}^{z}+\frac{1}{2}$. One can directly use the expression (12) for calculating the local magnetization $\langle S_{j}^{z}\rangle$ in a system with only few domain walls. The energy cost of creating a domain wall is large for $\Delta\gg 1$ and one can expect $\langle c_{l-1}^{\dagger}c_{l+1}^{\dagger}\rangle~{}\approx~{}0$. In Hamiltonian (11), the first term on the second line can thus be neglected in this limit. When there is no creation or annihilation of domain walls, the last term of Eq. (11) is a constant shift in energy and can also be left out.

When the density of domain walls is very low, $\langle c_{l}^{\dagger}c_{l}\rangle~{}\approx~{}0$, one can approximate

[TABLE]

on the first line of Eq. (11). This approximation removes the interaction and leads to the free fermion Hamiltonian

[TABLE]

This Hamiltonian is an excellent approximation of (11) as long as the domain walls are far from each other. However, using Hamiltonian (14) directly to compute the ground state density distribution would be inadequate. It separates into two decoupled Hamiltonians for the odd and even bonds, $H_{b}=H_{\text{odd}}+H_{\text{even}}$. The Hamiltonian $H_{b}$ therefore has degenerate pairs of eigenstates corresponding to $H_{\text{odd}}$ and $H_{\text{even}}$, with nearly identical wavefunctions. The wavefunctions corresponding to two eigenstates with the same energy have a similar envelope, but the ones which correspond to the eigenstates of $H_{\text{odd}}$ are only nonzero on the odd bonds, and respectively for $H_{\text{even}}$. In this description, the domain walls in the two degenerate states are discernible, and fermionic statistics do not exclude two domain walls in the same region of space. In Hamiltonian (11), the interaction on the first line prevents the DWs from crossing and ensures the correct fermionic statistics.

When the number of domain walls is small, their wavelength is very large. We can therefore take the continuum limit and ignore the microscopic lattice, since the envelopes of the wavefunctions vary at a much larger lengthscale. In this limit, we can consider the non-degenerate eigenstates of a particle in a box,

[TABLE]

with the momenta $k=k_{m}=\frac{m\pi}{L}$, $m=1,2,\cdots,L$. We thus consider indistinguishable domain walls which can be located on all bonds. Using these states excludes properly two domain walls from the same region of space by the effect of fermionic statistics even without the interaction term.

The fermion operators have the momentum representation

[TABLE]

Substituting Eq. (16), the Hamiltonian of Eq. (14) can be written as

[TABLE]

where $\epsilon_{k}=J_{xy}\cos(2k)$. For this dispersion, there is a maximum of energy at $k=0$. It is more convenient to have the minimum of energy at $k=0$, and therefore we shift the momentum by $\frac{\pi}{2}$:

[TABLE]

which changes $H_{b}\rightarrow-H_{b}$. The basis functions shift as

[TABLE]

The phase factor in $\varphi_{k}^{l}$ cancels in expectation values which contain $|\varphi_{k}^{l}|^{2}$. In Hamiltonian (17), the dispersion relation becomes

[TABLE]

The ground state of the free Hamiltonian is the many-body state of $N$ noninteracting fermions,

[TABLE]

Here, $l_{\alpha}$ denotes the coordinate of fermion $\alpha$. In the direct product, $\ket{l}$ denotes a single-particle state where site $l$ is occupied and other sites are empty,

[TABLE]

The coefficient $\varphi_{k_{1},\cdots,k_{N}}^{l_{1},\cdots,l_{N}}$ is given by the Slater determinant

[TABLE]

Equation (12) can be used for calculating the expectation value of magnetization $\langle S^{z}_{j}\rangle=\frac{1}{2}\langle\sigma_{j}^{z}\rangle$:

[TABLE]

which gives the expectation value of density as

[TABLE]

The number of fermions $N$ is equal to the number of domain walls $N=N_{\text{DW}}$. In the ground state, the number of DWs is given by the number of bath fermions $N_{b}$,

[TABLE]

In Fig. 9 with $\left(N_{b}-\frac{L}{2}\right)=2.5$, there are four DWs. The most probable locations of the DWs are seen as two neighboring maxima or minima in the density distribution of the upper panel. The result of Eq. (25) and the many-body TEBD solution are shown in Fig. 8

for $N=1$ and $L=60$. There is a very good agreement between the two solutions. When $N>1$, the accuracy of the approximation (14) should reduce since there can be interactions between the DWs. For $N=4$ and $L=61$, we find however a good agreement between Eq. (24) and the numerical TEBD result, as shown in Fig. 1(b) of the main text and in the upper panel of Fig. 9.

One can see the correspondence between solitons and free fermions also by comparing the many-body density distribution in the upper panel of Fig. 9 directly to the density distribution of free fermions

[TABLE]

shown in the lower panel of Fig. 9. The sum is over $k=\frac{\pi}{L+1},\cdots,k_{F}$, where $k_{F}=\frac{N\pi}{L+1}$, and $\ket{FS}=\prod_{k\leq k_{F}}c_{k}^{\dagger}\ket{0}$. The positions of neighboring minima or maxima in the many-body density distribution coincide with the maxima of the free fermion density, which has Friedel oscillations due to the boundaries.

A.3 Time evolution of the density distribution

In the Mott insulator phase, the time evolution of the density distribution after a local perturbation can be calculated in a similar way as the ground state density in Sec. A.2. We limit the discussion to the time evolution of one DW excitation. For incommensurate filling with more than one DW, one could not use the free-fermion description but would have to consider the interaction of the DWs.

As can be seen from Figs. 6 and 7, if the up spin at site $j_{0}$ in the original lattice is exchanged with the down spin at site $j_{0}+1$, the DWs are created at bonds $j_{0}-1$ and $j_{0}+1$. If the up spin at site $j_{0}$ is exchanged with the down spin at site $j_{0}-1$, the DWs are created at bonds $j_{0}-2$ and $j_{0}$. When considering processes in the same energy sector, a DW created at bond $j_{0}-2$ or $j_{0}-1$ can only move on the left, $l<j_{0}$, and a DW created at bond $j_{0}$ or $j_{0}+1$ can only move on the right, $l\geq j_{0}$. One can thus consider the time evolution of only half of the system. We consider an initial state where the DWs have already been created. The DW on the right side can be created at two neighboring bonds, $j_{0}$ or $j_{0}+1$. One can therefore write the initial state as a superposition

[TABLE]

The calculation is presented here for the right side of the lattice. The left side is symmetric.

The time-dependent state is obtained by operating with the time evolution operator,

[TABLE]

where $H_{b}$ is given by eqs. (17) and (20). We operate on the time-evolved state with the operator $\sigma^{z}$ of Eq. (12)

[TABLE]

The local magnetization is now

[TABLE]

Figure 2 of the main text shows the density distribution given by Eq. (31) together with the many-body TEBD result. The initial state is slightly different in the TEBD simulation than in the DW model, which leads to some differences in the time evolution. The motion of the DW in the density is however quite accurately given by the simple model.

Appendix B Density correlation

The complete inversion of density for the static barrier perturbation in Fig. 2 a) of the main text suggests that the soliton and antisoliton excitation propagate in opposite directions with a probability close to 1. The soliton and antisoliton distributions of Figs. 3 and 4 in the main text become zero at the center of the lattice, which is consistent with propagation to opposite directions. We can further investigate the symmetry in the motion of the excitations by computing the correlation

[TABLE]

shown in Fig. 10 (a).

Here, $i=j-j_{0}$ and

[TABLE]

The value $C_{i,-i}(t)=1$ means that the potential has produced an inversion of the magnetization at distance $i$ in a symmetric way, and the soliton and antisoliton excitations have propagated in opposite directions. A totally asymmetric time evolution where only one of the excitations propagates and the other one stays fixed would lead to $C_{i,-i}(t)=0$. Values between 0 and 1 result from averaging over different possible time evolutions. These cases are detailed with schematic diagrams in Appendix C of Visuri2 .

The correlation in Fig. 10 a) clearly indicates symmetric propagation. This in agreement with the DW model given in Sec. LABEL:sec:DW_model, which indeed leads to the creation of a domain wall on both sides of the center site. In the case of an impurity, the time evolution is more complicated. The impurity can move to the neighboring site and form a bound state with the antisoliton excitation, in which case only the soliton excitation propagates Visuri2 . This results in a reduced amplitude of the magnetization oscillation with respect to the initial state, and the correlation in Fig. 10 (c) has a smaller maximum value than in Fig. 10 (a). In both cases, however, one can see a light-cone propagation leading to zero correlation beyond a certain distance.

The correlation of Eq. (32) shown in Figs. 10(b) and 10(d) has smaller values than in the commensurate case since the particles are initially partly delocalized. The maximum of the correlation is again smaller for the impurity than for the static barrier, since a similar bound state of the impurity and the possible antisoliton excitation can occur here as in the commensurate case. The correlation has an additional maximum between 10 and 20 sites from the center, corresponding to the density profile.

Appendix C Smaller interaction $V=10J$

For smaller interactions $U=V=10J$, we obtain essentially the same results as for $U=V=50J$. Due to a smaller nearest-neighbor repulsion, the particles are more delocalized in the initial state. This leads to a slightly reduced amplitude of the density oscillation in the ground state at $t=0$ in Fig. 11.

Due to a larger effective tunneling $\frac{4J^{2}}{U}$, the impurity moves further, which results in a broader maximum of the antisoliton distribution at the center of the lattice in Fig. 4 of the main text. Essentially the same features as for $U=V=50J$ can be seen in the time evolution of the soliton and antisoliton distributions in Fig. 12:

The excitations propagate away from the center when a static barrier is created, whereas for the impurity, the soliton propagates while the antisoliton distribution has a maximum at the center. For the incommensurate filling shown in panels e) and f), the antisoliton excitations clearly propagate a shorter distance than in the MI.

As seen in Fig. 12, there is a small density of antisolitons in the ground state. This is seen in the total number of antisolitons in Fig. 13, which is larger than for $U=V=50J$. The number of antisolitons at $t=0$ and the increase during the simulation time are shown in Table 1. For $U=V=50J$, $N_{\text{AS}}$ is initially close to zero and the increase is approximately equal to the initial occupation of site $j_{0}$ in both the MI and incommensurate phases. The probability of not creating excitations is thus approximately equal to the probability of site $j_{0}$ being empty. For $U=V=10J$, the increase in the number of antisolitons is slightly smaller than the initial occupation of site $j_{0}$. When the particles are more delocalized, it is possible that both the sites $j_{0}$ and a site next to it are occupied. In this case, an additional soliton-antisoliton pair would not be created.

The presence of solitons can also be seen as DWs in the density distribution. The number of DWs $N_{\text{DW}}(t)$ in the density distribution is equal to the number of neighboring filled sites $N_{\text{S}}(t)$ for $V\rightarrow\infty$. For a finite $V$, more than $N_{\text{DW}}$ pairs of neighboring occupied sites can exist in the ground state due to delocalization, which leads to a smaller amplitude of the density oscillation. Correspondingly, the number of antisolitons in the ground state would be zero for $V\rightarrow\infty$ but can be nonzero for a finite $V$. We have verified that $N_{\text{AS}}(t)\approx N_{\text{S}}(t)$ in the MI state, and $N_{\text{AS}}(t)\approx N_{\text{S}}(t)-N_{\text{DW}}(0)$ for incommensurate filling, with differences of order $10^{-3}$ for $U=V=50J$ and 0.08 for $U=V=10J$. We attribute the slightly larger soliton number to a small asymmetry in the system: a nondegenerate state for odd $L$ at half filling requires one filled site more than the number of empty sites.

Appendix D Harmonic trap

The ground state of the bath with a harmonic potential is described by the Hamiltonian

[TABLE]

where

[TABLE]

and $V_{\text{trap}}(j)=V_{0}(j-j_{0})^{2}$. For a sufficiently small $V_{0}$, i.e. when $V_{0}(j-j_{0})^{2}\lesssim 2V$ at $j=j_{0}+N_{b}-1$, the ground state in the central region of the trap is a Mott insulator, leading to the same dynamics as in a uniform system with open boundary conditions. The ground state density distribution is shown in Fig. 14.

For the ground state to be incommensurate, the particles must be confined by hard boundaries. The lower panel of Fig. 14 shows the ground state density distribution for a stronger harmonic confinement, where the energy is minimized by a higher particle density in the center of the trap surrounded by Mott insulator regions. This distribution is qualitatively different from the incommensurate case of Fig. 9, and leads to different dynamics.

We find that when a shallow harmonic potential is superimposed on a system limited by hard boundaries, the particle density at incommensurate filling is similar to the one in a uniform box but the solitons are distributed closer to the center, as shown in Fig. 15.

At half filling, the ground state density distribution is not modified significantly by the harmonic potential when $V_{\text{trap}}(1)=V_{\text{trap}}(L)\ll 2V$. In the time evolution, $H=H_{\text{b}}+H_{\text{trap}}+H^{\prime}$, where $H_{\text{trap}}$ is modified in the case of the impurity,

[TABLE]

The time-dependent antisoliton distributions in Fig. 16 show that at incommensurate filling, the antisoliton is confined to a smaller region around $j_{0}$ than in Figs. 12(e) and (f).