Dynamics of trapped one-dimensional bosons for intermediate repulsive interactions

TL;DR

This paper investigates the dynamics of a few one-dimensional bosons under harmonic confinement with varying interaction strengths, revealing crossovers in motion period and changes in Bragg-scattering peaks across different regimes.

Contribution

It provides a numerical analysis of the time-evolution of few-body 1D bosons for arbitrary interaction strengths, highlighting new dynamical behaviors and structural features.

Findings

The motion period exhibits two crossovers as interaction strength varies.

Bragg-scattering peaks depend on interaction strength, from weak to Tonks-Girardeau regime.

Numerical methods successfully model dynamics for arbitrary two-body interactions.

Abstract

The time-evolution of few number of interacting, harmonically confined one-dimensional bosons is numerically obtained for arbitrary two-body $\delta-$potential interaction strengths. It is demonstrated that the period of the motion in a Newton's cradle configuration undergoes two crossovers as a function of interactions. Furthermore, through the evaluation of the structure factor, the dependence of Bragg-scattering peaks on the interaction strength ranging from the weak coupling regime to the impenetrable Tonks-Girardeau case is illustrated.

Equations

$$\mathcal{H}(\mathbf{r})=\sum _{i=1}^N\left[-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x_i^2}+\frac{1}{2}m{\omega }_0^2{x}_i^2\right]+c\sum _{i\le j}^N\delta (x_i-x_j),$$

$${\psi }_1(x_1)=A\exp \left[-\frac{(x_1+\alpha {)}^2}{2{\sigma }^2}\right],$$

$${\psi }_i(x_i)=A\exp \left[-\frac{[x_i+2(N/2-i)ϵ{]}^2}{2{\sigma }^2}\right],$$

$$\begin{array}{rl}U(\Delta t)& =e^{-\frac{i}{\hslash }({\mathcal{H}}_{\mathrm{k}}+{\mathcal{H}}_{\mathrm{p}})\Delta t}\\ & \simeq U_{\mathrm{TS}}(\Delta t)\equiv e^{-\frac{i}{2\hslash }{\mathcal{H}}_{\mathrm{p}}\Delta t}{e}^{-\frac{i}{\hslash }{\mathcal{H}}_{\mathrm{k}}\Delta t}{e}^{-\frac{i}{2\hslash }{\mathcal{H}}_{\mathrm{p}}\Delta t}.\end{array}$$

$$\begin{array}{rl}& \rho (x,x^{\prime },t)\\ & =\int {\Psi }^{\ast }(x,\cdots ,x_N,t)\Psi (x^{\prime },\cdots ,x_N,t)d{x}_2\cdots d{x}_N.\end{array}$$

$$\mathbf{J}(\mathbf{r},t)=\sum _{i=1}^N{\hat{\mathbf{x}}}_i{J}_i(\mathbf{r},t)=\frac{\hslash }{2m}\Im \{{\Psi }^{\ast }(\mathbf{r},t)\nabla \Psi (\mathbf{r},t)\},$$

$$\begin{array}{rl}{I}_{i\to j}(t)=& \frac{\hslash }{2m}\int dx\cdots d{x}_{j-1}d{x}_{j+1}\cdots d{x}_N\\ & \times \Im {\left\{{\Psi }^{\ast }(\mathbf{r},t)\frac{\partial }{\partial x_i}\Psi (\mathbf{r},t)\right\}}_{x_i=x_j},\end{array}$$

$$S(k,\omega )=\int e^{i(\omega t-kx)}⟨\rho (x,t)\rho (0,0)⟩dxdt.$$

$$f(\Delta k,t)\propto \int e^{-ix\Delta k}\rho (x,t)dx,$$

Full text

Experiments on one-dimensional (1D) bose gases Görlitz et al. (2001) demonstrated peculiar properties that are missing in their higher dimensional counterparts. In the absence of an external trap, the 1D bose gas with two-body contact interactions of any strength is described by the exactly integrable Lieb-Liniger (LL) model Lieb and Liniger (1963); *LL2. Through the work of Olshanii Olshanii (1998), it has became possible to relate the interaction strength parameter $c$ of LL model to the real experimental parameters of ultracold gases in 1D traps. In practice, $c$ can be tuned via Feshbach resonances Chin et al. (2010). Since then, the infinitely strong repulsive interaction limit of the trapped LL gas, also known as the Tonks-Girardeau (TG) regime Girardeau and Wright (2002); *fermi-bose-review, has been experimentally realized in optical traps Paredes et al. (2004); Kinoshita et al. (2006); *Kinoshita1125. Recently, an experimental implementation of the quantum Newton’s cradle of tunable interaction strength has been achieved using the dipole-diple interaction between highly magnetic dysprosium atoms Tang et al. (2018).

In spite of the rich history backing the subject Cazalilla et al. (2011), most of the studies on the thermodynamic and spectral properties of interacting 1D bose gases have focused on the limiting cases of impenetrable TG regime Atas et al. (2017); De Rosi et al. (2017); Pezer and Buljan (2007); Lang et al. (2015) and the weakly interacting regime which is well described by Bogoliubov theory. As shown in Fig. 1, the strongly interacting regime broadens the Bogoliubov spectrum into a continuum between two types of modes. Type I modes are bosonic quasiparticle modes and the Type II modes are fermionic quasihole modes Lieb (1963).

As far as the present authors can ascertain, it has not yet been resolved whether the finite mutual interaction strength compromises the integrability of a system of bosons in the presence of an external potential. On the other hand, numerical approaches proved powerful in shedding light into the crossover between weak interaction and TG regimes Zöllner et al. (2006); *PhysRevLett.100.040401. In particular, the intermediate window of momentum distributions has only been investigated numerically Astrakharchik and Giorgini (2003); Schmidt and Fleischhauer (2007), despite not being able to describe all physics from infrared to the ultraviolet range. In this Letter, three main results are presented: an analysis of the real-time dynamics reduced density matrix of a trapped 1D few-boson system, the dependence of the motion period on the interaction strength $c$, and the full range structure factor at finite $c$.

Model and the numerical method.

In 1D, $N$ indistinguishable bosons of mass $m$ with repulsive two-body contact interaction $c\geq 0$ inside a harmonic trap of characteristic frequency $\omega_{0}$ can be modelled by the Hamiltonian

[TABLE]

where $\mathbf{r}=(x,x_{2},\cdots,x_{N})$. It is convenient to define a unitless interaction parameter $\gamma=2c\left(\frac{\hbar^{2}}{2m\Delta x^{2}}\right)^{-1}$, where $\Delta x$ is the spatial width that is used in numerical discretization. In the numerical implementation, the $\delta-$potential is alleviated by using a more realistic normalized interaction with a small, but finite range. Normalization ensures that in the zero-width limit it approaches to the $\delta-$distribution.

The symmetrized initial state $\Psi(\mathbf{r},t=0)$ is composed of the Slater permanent of $N$ localized wavepackets $\psi_{j}$, with $j\in[1,N]$. In particular, the Newton’s cradle configuration comprises one swinging wavepacket displaced by $\alpha$ with respect to the equilibrium position

[TABLE]

corresponding to the lifted ball in the cradle. The remaining $N-1$ wavepackets are nearly stationary and centered around the equilibrium, each separated by small distance $\epsilon\ll\alpha$

[TABLE]

with $i\in[2,N]$. Here, $A$ is the normalization, and $\sigma$ is the wavepacket width. Moreover it is taken that $\epsilon/\sigma\approx 1$, because $\epsilon$ is related to an effective radius of the bosons which tends to $\sigma$ as $\gamma\rightarrow\infty$. The goal is to study the dynamics of this system for different values of the interaction parameter $\gamma$ by calculating the time evolution of $\Psi(\mathbf{r},t)$ as governed by the time dependent Schrödinger equation (TDSE).

Diffusion Monte Carlo techniques Zhang et al. (2014); Kalos et al. (1974); Astrakharchik and Giorgini (2003); Gudyma et al. (2015) and density matrix renormalization group Schmidt and Fleischhauer (2007) are common methods to extract the ground state and low-energy excitations in interacting 1D bosonic systems. In addition, the few-body case with a double-well trap has been investigated using the Multi-Configuration Time-Dependent Hartree Meyer et al. (1990) method by Zöllner et al. Zöllner et al. (2006); *PhysRevLett.100.040401. In this study, on the other hand, the $N$-particle TDSE $\mathcal{H}(\mathbf{r},t)\Psi(\mathbf{r},t)=i\hbar\frac{\partial}{\partial t}\Psi(\mathbf{r},t)$ is numerically integrated for low $N$ via a split operator method. Here, the time evolution according to $\mathcal{H}$ of Eq. 1 is obtained by alternated advancement in real and Fourier spaces for $\Psi(\mathbf{r},t=0)=\sum_{\mathcal{P}}\prod_{j=1}^{N}\,\psi_{j}(x_{i})$ via Equations 2, 3. Here, $\mathcal{P}$ denotes the particle permutations. The premise of this method is using the fast Fourier transform to diagonalize the propagator $U(\Delta t)$ by breaking it up into the kinetic ($\mathcal{H}_{\rm k}$) and the potential ($\mathcal{H}_{\rm p}$) parts using Baker-Campbell-Hausdorff identity

[TABLE]

The approximate expression $U_{\rm TS}(\Delta t)$ is also known as the Trotter-Suzuki expansion Trotter (1959); *suzuki1976generalized. The wavefunction is iteratively advanced in time by $\Psi(\mathbf{r},\Delta t)=U_{\rm TS}(\Delta t)\Psi(\mathbf{r},t=0)$ at the cost of an error proportional to $\Delta t^{3}$. Here, it is chosen that $\Delta t=t/n$ for the total number $n$ of time steps of the simulation. The expansion approaches the exact expression as $n\rightarrow\infty$ due to the Lie product formula.

The $N$-particle wave-function $\Psi(\mathbf{r},t)$ is computed over several cycles of motion for specific values of $\gamma$. Subsequently, the reduced real-space density matrix $\rho(x,x^{\prime},t)$ is given by

[TABLE]

In what follows, the case where $x=x^{\prime}$ will be considered, such that $\rho(x,x,t):=\rho(x,t)$ which describes the local density profile of the system.

Time evolution of density profiles.

The simulated space-time profile of $\rho(x,t)$ for $N=3$ is plotted in Fig. 2, and in the left panel of Fig. 2.b for $\gamma=0.3$ and for $\gamma=3$, which is one of the main results of this Letter. Here, the time and the space axes are respectively scaled in terms of the motion period $\tau_{0}=2\pi/\omega_{0}$ and the oscillator length $l_{0}=\sqrt{\hbar/m\omega_{0}}$, both of which calculated for $\gamma=0$ case. It is seen that in the limit $\gamma\to\infty$, for the aformentioned inital conditions, $\rho(x,t)$ assumes a periodic profile reminiscent of the worldlines of a classical Newton’s cradle. On the other hand, tuning $\gamma$ results in alterations in the general form of the real-space density, as well as the motion period.

First observation is a qualitative one; the interference between the spacetime densities of the particles reduces with increasing $\gamma$. This can be seen for $N=3$ in the left panel of Fig. 2. It is seen that for $\gamma=3$ the real-space densities of the two central bosons cease to overlap. The depletion of the overlap region leads bosons to isolate each other, as $N=3$ humps emerge in $\rho(x,t)$. This process, referred to as fragmentation, also signals that bosons are more likely to get reflected upon collisions. In fact, the $\gamma\to\infty$ limit imitates the Pauli exclusion principle, as fragmentation saturates into a so-called fermionized state Girardeau (1960). On the other hand, for $\gamma=0.3$, such an overlap occurs at each event of collision, thereby the $N=3$ humps in $\rho(x,t)$ merge into a single one (see Fig. 2), indicating a finite transmission amplitude.

In addition to fragmentation and the alterations in transmission rates, the motion period deviates from $\tau_{0}$ as a function of $\gamma$ as shown in the lower panel of Fig. 3.b. This is the second main result of the present work. Here, two crossovers for the frequency deviations between zero to weak interaction, and weak to strong interaction are demonstrated. The results indicate that a maximal deviation at an intermediate value $\gamma\approx 0.4$ is followed by a return towards the initial value. Moreover, a minor reduction of $\tau$ at strong interactions signals a scattering resonance. The greatest deviation from $\tau_{0}$ is determined to be 5.8%. Note that $\tau_{0}$ is equivalent to the oscillation period of a single particle, which is readily attained from the trap parameters in experiments. The measurement of deviations should reveal the empirical connection of the numerical parameter $\gamma$. The present results are in accordance with the previous experimental Moritz et al. (2003); Haller et al. (2009); Fang et al. (2014) results along with the numerical Schmitz et al. (2013) and theoretical Menotti and Stringari (2002); Gudyma et al. (2015) studies focusing on either the lower or the higher crossover for breathing modes.

Tunneling rates.

The time dependent amplitude profiles given in Fig. 2 by themselves do not yield quantitative information about the transmission rates of particles at a collision event, which occurs when their wavepackets overlap. To obtain this quantity, a convenient starting point is to calculate the rate of change of the probability of finding a particle in state $\Psi(\mathbf{r},t)$, i.e. the probability current. From the continuity equation, the total probability current density $\mathbf{J}(\mathbf{r},t)$ can be found to be

[TABLE]

where $\hat{\mathbf{x}}_{i}$ is the unit vector along the coordinate $x_{i}$, and $\Im$ denotes the imaginary part. The components of $\mathbf{J}$ can be associated with each of the $N$ bosons. Therefore, the partial probability density $J_{i}(x,\cdots,x_{i},x_{j}=x_{i}\cdots,x_{N},t)$ can be interpreted as a measure of the probability current density of $i$th particle through $j$th particle upon their head on collision. In other words, by calculating the probability current $I_{i\rightarrow j}(t)$, which is the total flux of $\mathbf{J}$ through the interaction plane $x_{i}=x_{j}$ given by

[TABLE]

one can quantify the likelihood of transmission of one boson through another at a collision. Simulated results for $I_{i\rightarrow j}(t)$ corresponding to different values of $\gamma$ are presented in Fig. 3.a. It is found that the transmission probability $\mathcal{T}(t)\propto|I_{i\rightarrow j}(t)|^{2}$ decreases with increasing $\gamma$ and saturates at zero for $\gamma\gtrsim 3$, corresponding to impenetrable bosons. Being in a low energy scattering regime, this limit, where $\mathcal{T}(t)$ approaches to zero, corresponds to the TG gas, as suggested by OlshaniiOlshanii (1998).

The structure factor.

The dynamical structure factor (DSF) is defined as

[TABLE]

Previously this has been analytically calculated for the LL gas in the absence of an external trap in Caux et al.Caux and Calabrese (2006), and the static structure factor $S(k)=\int\frac{d\omega}{2\pi}S(k,\omega)$ has been obtained using quantum Monte Carlo techniques Astrakharchik and Giorgini (2003). Our numerical result for the DSF at interaction strengths $\gamma=0.003$ and $1000$ is presented in Fig. 4. Here, the $\omega$ axis corresponds to the energy transfer and $k$ corresponds to momentum. According to the results, for weak interactions most of spectral weight of $S(k,\omega)$ is found in the vicinity of a type-I excitation. For strongly interacting bosons, the spectrum is marginally broadened. On the other hand, being far away from the thermodynamic limit ($N=3$) and due to the low resolution, a rigorous comparison with Boguliubov theory would be greatly ambitious.

From an experimental point of view, the dynamics of the system can be characterized by Bragg spectroscopy Sette et al. (1998). One can treat the bosonic system as a perturbative potential that transfers a momentum of $\hbar\Delta k$ to the incoming light. Then the corresponding scattering amplitude can be quantified by the inelastic, time dependent structure factor Kitanine et al. (2011)

[TABLE]

and it can be probed by light scattering experiments. The time evolution of the norm of the structure factor $F(\Delta k,t)=|f(\Delta k,t)|^{2}$ is given in the right panel of Fig. 2.b, which brings us to the third and the final main result presented in this Letter. It is found that the profile of the structure factor retains similar characteristics for different values of $\gamma$. Nonetheless, it is reported that for weaker interactions, the weight of $F(\Delta k,t)$ shifts towards smaller momentum components. This qualitative remark is reflected in the right panel of Fig. 4. It is seen that the most profound contrast between $\gamma=0.3$ and $\gamma=3$ cases occurs at times equal to odd multiples of $\tau_{0}/4$, corresponding to an overlap of the wave-packets. $\mathcal{T}(t)$ is related to the spread of the scattering amplitude at $t=\tau_{0}/4$. This suggests that, despite the indistinguishability of the particles, it is possible to acquire insight into a collision event of bosons through Bragg spectroscopy.

To obtain a quantitative account of the collision dynamics, the structure factor can be investigated for different values of $\gamma$ by focusing on $F(\Delta k,t)$ at specific time slices. In Fig. 5, $F(\Delta k,t)$ is plotted (a) at $t=\tau_{0}/4$, i.e. at the moment of three-body collision, and (b) at $t=\tau_{0}/2$, i.e. when the initially displaced particle is furthest from the equilibrium position of the external trap. Here, for both cases, the basic observation is that high momentum transfer rates monotonically increase with increasing $\gamma$. The tails of the structure factor escalates with increasing $\gamma$ in accordance with Caux et al.Caux et al. (2007) and Minguzzi et al.Minguzzi et al. (2002). On the other hand, given that the infrared behavior is usually obscured by the harmonic trap Paredes et al. (2004), it is also essential to note that these results also shed light into the previously concealed intermediate window of momenta.

Conclusion.

To summarize, the evolution of the reduced real-space density and the one-dimensional structure factor are computed for specific finite values of the two-body interaction strength $c$. The dependence of the motion period and the transmission amplitude of colliding bosons on the numerical interaction parameter $\gamma$ is demonstrated. It is suggested that the period of motion deviates from the characteristic period of the harmonic trap at a maximum about 6%. From the empirical perspective, this constitutes a method for establishing the link between $\gamma$ and the experimental parameters.