Entanglement and entropy production in coupled single-mode Bose-Einstein condensates

TL;DR

This paper studies how entanglement entropy evolves over time in coupled Bose-Einstein condensates, revealing entropy oscillations, saturation, and a phase transition linked to self-trapping, with results matching classical predictions.

Contribution

It combines numerical and analytical methods to analyze entropy dynamics in coupled condensates, highlighting the connection to semiclassical behavior and phase transitions.

Findings

Entropy oscillates due to coherent condensate dynamics.

Entropy saturates at a value indicating dephasing without full equilibration.

Stationary entropy aligns with classical microcanonical ensemble predictions.

Abstract

We investigate the time evolution of the entanglement entropy of coupled single-mode Bose-Einstein condensates in a double well potential at $T=0$ temperature, by combining numerical results with analytical approximations. We find that the coherent oscillations of the condensates result in entropy oscillations on the top of a linear entropy generation at short time scales. Due to dephasing, the entropy eventually saturates to a stationary value, in spite of the lack of equilibration. We show that this long time limit of the entropy reflects the semiclassical dynamics of the system, revealing the self-trapping phase transition of the condensates at large interaction strength by a sudden entropy jump. We compare the stationary limit of the entropy to the prediction of a classical microcanonical ensemble, and find surprisingly good agreement in spite of the non-equilibrium state of the system. Our predictions should be experimentally observable on a Bose-Einstein condensate in a double well potential or on a two-component condensate with inter-state coupling.

Equations

$$\hat{H}=-J\left({\hat{a}}_L^{†}{\hat{a}}_R+{\hat{a}}_R^{†}{\hat{a}}_L\right)+\frac{U}{2}\left({\hat{N}}_L^2-{\hat{N}}_L+{\hat{N}}_R^2-{\hat{N}}_R\right).$$

$$\mathcal{S}(t)=-\sum _{n_L=0}^N{P}_t(n_L)\log P_t(n_L),$$

$${\hat{S}}_z=\frac{1}{2}\left({\hat{N}}_L-{\hat{N}}_R\right),{\hat{S}}_x=\frac{1}{2}\left({\hat{a}}_L^{†}{\hat{a}}_R+{\hat{a}}_R^{†}{\hat{a}}_L\right),$$

$$\hat{H}=-2J{\hat{S}}_x+U{\hat{S}}_z^2.$$

$$\mathcal{S}(t)=-\sum _{m=-N/2}^{N/2}{P}_t(m)\log P_t(m),$$

$$\alpha \equiv \frac{NU}{2J},$$

$$\vec{S}=\frac{N}{2}\vec{\Omega }.$$

$${\partial }_t{\Omega }_x=-UN{\Omega }_y{\Omega }_z,$$

$${\partial }_t{\Omega }_y=2J{\Omega }_z+UN{\Omega }_x{\Omega }_z,$$

$${\partial }_t{\Omega }_z=-2J{\Omega }_y,$$

$${\Omega }_x=-\frac{1}{\alpha },{\Omega }_y=0,{\Omega }_z=\pm \sqrt{1-\frac{1}{{\alpha }^2}}.$$

$${\left({\Omega }_x+\frac{1}{\alpha }\right)}^2+{\Omega }_y^2=const.$$

$${\alpha }_{\theta }=\frac{2}{1-\sin \theta }.$$

$${\alpha }_{\theta }^{\mathrm{fix}}=-\frac{1}{\sin \theta },$$

$$|{\Omega }_{\theta }⟩=$$

$$\sum _{m=-N/2}^{N/2}\sqrt{\left(\genfrac{}{}{0ex}{}{N}{m+\frac{N}{2}}\right)}{\left(\cos \frac{\theta }{2}\right)}^{m+\frac{N}{2}}{\left(\sin \frac{\theta }{2}\right)}^{\frac{N}{2}-m}|m⟩,$$

$$|{\Omega }_{\theta }⟩\approx \sqrt{\frac{2}{N\pi {\sin }^2\theta }}\sum _m\exp \left(-\frac{{\left(m-\frac{N}{2}\cos \theta \right)}^2}{N{\sin }^2\theta }\right)|m⟩,$$

$$\mathcal{S}(t=0)\approx \frac{1}{2}\log \frac{\pi eN{\sin }^2\theta }{2}.$$

$$\tilde{\mathcal{S}}(t)=\mathcal{S}(t)-\frac{1}{2}\log \frac{\pi eN{\sin }^2\theta }{2}.$$

$$|\psi (t)⟩=\sum _m{e}^{-i\phi (t)m}{c}_m(t)|m⟩.$$

$$c_m(t)\to \psi (x,t)\equiv$$

$${\left(\frac{2N\mathrm{Re}c(t)}{\pi }\right)}^{1/4}\exp (-c(t)N(x-x_0(t){)}^2).$$

$${\mathcal{S}}_G(t)=-\frac{1}{2}\log \left(4{\sin }^2\theta \mathrm{Re}c(t)\right).$$

$$\delta ⟨{\psi }_G|i{\partial }_t-\hat{H}|{\psi }_G⟩=0,$$

$$⟨{\psi }_G|i{\partial }_t-\hat{H}|{\psi }_G⟩=$$

$${\partial }_t\phi (t)\frac{N}{2}\int \mathrm{d}xx|\psi (x,t){|}^2+i\int \mathrm{d}x{\psi }^{\ast }(x,t){\partial }_t\psi (x,t)$$

$$-U\frac{N^2}{4}\int \mathrm{d}x{x}^2|\psi (x,t){|}^2+\frac{JN}{2}(e^{i\phi (t)}\times$$

$$\int \mathrm{d}x{\psi }^{\ast }(x,t)\psi (x-\frac{2}{N},t)\sqrt{1-x^2+\frac{2}{N}(1+x)}+c.c.).$$

$${\partial }_t{x}_0=-2J\sqrt{1-x_0^2}\sin \phi ,$$

$${\partial }_t\phi =UN{x}_0+2J\frac{x_0}{\sqrt{1-x_0^2}}\cos \phi .$$

$$\vec{\Omega }=(\sqrt{1-x_0^2}\cos \phi ,\sqrt{1-x_0^2}\sin \phi ,x_0).$$

$$i{\partial }_{t}c=-\frac{\alpha J}{2}-\frac{J\cos \phi }{2(1-x_0^2{)}^{3/2}}$$

$$-4J\frac{x_0}{\sqrt{1-x_0^2}}\sin \phi c+8J\cos \phi \sqrt{1-x_0^2}{c}^2.$$

$$\overline{\mathcal{S}}=\frac{1}{T}{\int }_0^T\mathrm{d}t\mathcal{S}(t),$$

$${\mathcal{S}}_{\mathrm{class}}=-\int \mathrm{d}{\Omega }_{z}P({\Omega }_z)\log P({\Omega }_z)+{\mathcal{S}}_0.$$

Full text

I Introduction

Entanglement is a fundamental concept of quantum mechanics, manifesting in strong, non-local correlations between subsystems. Constituting one of the most crucial differences between classical and quantum physics, entanglement is studied in a diverse area of physics, ranging from quantum gravity blackhole to topological order in condensed matter systems topology . In recent years, entanglement generation in non-equilibrium many-body systems received a special attention, due to the intimate connection between entanglement spreading and the equilibration in closed systems NonEquilibriumReview ; NonEquilibrium ; schmiedmayer ; trotzky . Even for globally pure quantum states, the generation of strong entanglement between subsystems allows the thermalization of an isolated quantum system under its own coherent dynamics in a sense that measurements of local observables become indistinguishable from the predictions of an equilibrium thermal ensemble thermalization ; thermalization2 ; thermalization3 . In contrast, the large number of conserved quantities in integrable systems can prevent entanglement spreading, and result in the failure of thermalization. The slow, logarithmic increase of entanglement has been suggested as a fingerprint of non-ergodic many-body localized phases mbllog ; mbllog2 ; mbllog3 , whereas delocalized phases are characterized by a linear, light cone-like propagation of correlations lightcone ; ballistic .

Despite its fundamental importance, the experimental investigation of entanglement in correlated many-body systems remains challenging, since it usually requires information on the full quantum state. However, the swift experimental progress in recent years opened up unprecedented possibilities to study entanglement in ultracold atomic settings ZwergerReview . Site-resolved control of ultracold atoms in optical lattices allowed the direct measurement of Rényi entanglement entropy and mutual information lightcone2 ; greiner1 ; bosonentanglement , as well as the investigation of the intimate relation between the quantum purity of subsystems and the thermalization of an isolated non-equilibrium system greiner2 .

Entanglement is also at the heart of entropy production in closed quantum systems deutsch ; rigol . Taking two coupled quantum systems, even coherent evolution of the whole system produces entanglement entropy for each subsystem and may lead to equilibration. Two coupled single mode condensates provide one of the simplest examples to study this phenomenon in detail. In contrast to the case of small subsystems  lightcone2 ; greiner1 ; bosonentanglement , here the two coupled subsystems are equally large, and we find that equilibration can be understood by approximating the state of the full system by a microcanonical ensemble rather than a thermal Gibbs ensemble. This system can be realized by loading a Bose-Einstein condensate into a double well potential (see Fig. 1). Assuming that the atoms in the left and right wells condense into a single wave function, the dynamics is governed by the Hamiltonian milburn

[TABLE]

Here the bosonic operators $\hat{a}_{L}^{\dagger}$ and $\hat{a}_{R}^{\dagger}$ create particles into the left and right potential wells respectively, and $\hat{N}_{i}=\hat{a}_{i}^{\dagger}\hat{a}_{i}$ for $i=L,R$. The first term in the Hamiltonian describes the tunneling of particles, while the second term takes into account the interaction between the bosons in the same potential well. For given total particle number $N$, the entanglement entropy between the left and right wells is simply given by entanglement

[TABLE]

where $P_{t}(n_{L})$ denotes the probability of state $\hat{N}_{L}=n_{L}$ at time $t$ footnote0 . Here, concentrating on the effect of dephasing during the coherent, unitary time evolution of a closed quantum system, we investigate the time dependence of the entanglement entropy $\mathcal{S}(t)$ at $T=0$ temperature. Importantly, while usually an entropy measurement would require detailed knowledge of a complicated quantum state, here the full time evolution of $\mathcal{S}(t)$ can be investigated experimentally, since it only requires measuring the number of particles in the left and right wells. In contrast to earlier entropy measurements in optical lattices, involving a small sublattice with only a few atoms, coupled single-mode Bose-Einstein condensates would allow to study entanglement in large correlated many-body systems oberthalerentanglement ; riedelentanglement ; olsen .

Let us note that besides the double well experiment illustrated above, the Hamiltonian (1) can also be realized in a two component condensate trapped in a single well. E.g. two atomic hyperfine states forming the condensates may be coupled through microwaves hall , while their interaction may be tuned using a Feshbach resonance feshbachrev .

Hamiltonian (1) can also be rewritten in a more convenient form. Using the Schwinger boson representation milburn , we introduce spin operators

[TABLE]

of length $N/2$, with $N$ denoting the total number of particles. Apart from a redundant constant term, $\hat{H}$ can then be expressed as

[TABLE]

In this new representation, the entanglement entropy between the left and right wells corresponds to the entropy associated with $\hat{S}_{z}$ footnote ,

[TABLE]

with $P_{t}(m)$ denoting the probability of state $\hat{S}_{z}=m$ at time $t$.

Let us note that the spin Hamiltonian Eq. (3) is a special case of the Lipkin-Meshkov-Glick model, describing mutually interacting spin-1/2 particles, embedded in a magnetic field  LMG . In this context $\hat{S}_{\alpha}=\sum_{i}\hat{\sigma}_{i}^{\alpha}/2$ is the total spin operator, with $\hat{\sigma}_{i}^{\alpha}$ denoting the Pauli matrices at site $i$ for $\alpha=x,y,z$. Depending on the the strength of the magnetic field, the Lipkin-Meshkov-Glick model shows a second order quantum phase transition. The entanglement properties of the ground state of this system have been analyzed by calculating the von Neumann entropy of a subsystem consisting of $L$ sites. In particular, it has been shown that the entanglement entropy shows a logarithmic divergence at the critical point of the quantum phase transition  LMG1 ; QPT . Similar divergence in the entanglement properties of the ground state at the critical point has also been observed in other systems, like the Dicke model or the transverse field Ising model  dicke ; fisherinfo . Moreover, the dynamics of the von Neumann entropy of a single spin in the Lipkin-Meshkov-Glick model has also been investigated  LMG2 . In this work we concentrate on a different type of entanglement entropy, associated with the spin operator $\hat{S}_{z}$.

Our main purpose here is to analyze the time evolution of entropy (4) for different initial states and interaction strengths, by combining numerical results with analytical calculations. We demonstrate that $\mathcal{S}(t)$ exhibits coherent oscillations, reflecting the quantum mechanical dynamics of the coupled single-mode condensates. At the same time, $\mathcal{S}(t)$ shows a steady increase, and eventually reaches a stationary, ”equilibrium” value, even though this closed system always remains in a pure state.

The dynamics of the system depends crucially on the dimensionless parameter

[TABLE]

characterizing the strength of interactions dynamics . For $\alpha<1$, the average population imbalance between the two potential wells, $N_{L}-N_{R}$, oscillates between positive and negative values. For $\alpha>1$, however, the system undergoes a self-trapping transition oberthaler ; oberthalersemiclass . Here, for large initial particle number imbalance $N_{L}-N_{R}$ and strong enough interactions $\alpha\gg 1$, the interaction energy of the initial state prevents levelling off the number of particles in the two wells, and the amplitude of population imbalance oscillations is suppressed (see Sec. II for more details).

We show a typical example of entropy production, e.g. steady increase of the entropy, in Fig 2. The initial state of the system corresponds to maximal population imbalance, $|\hat{S}_{z}=N/2\rangle$, and the time evolution of $\mathcal{S}$ is calculated numerically by exact diagonalization. The entropy shows oscillations on the top of a steady increase, before saturating to a constant value. Moreover, the curves corresponding to the same $\alpha$, but different total particle numbers can be scaled together. As we discussed earlier, these oscillations during entropy production should be experimentally accessible (for a discussion of experimental parameters see Sec. V).

The long time limit of the entropy also reflects the self-trapping transition by showing a sudden jump at the ”phase boundary”. Interestingly, the computed asymptotic entropy value agrees well with the predictions of a classical microcanonical ensemble, where the normalized spin vector $\vec{\Omega}\equiv 2\vec{S}/N$ is distributed uniformly along a classical trajectory. The remarkable success of classical description can be understood by investigating the overlap between the wave function $|\psi(t)\rangle$ and the spin coherent states $|\Omega_{\theta,\varphi}\rangle$, polarized into direction $(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta)$. We plotted this overlap on the unit sphere in Fig. 3, for a maximally polarized initial state $\hat{S}_{z}=N/2$, with the interaction strength tuned to the boundary of self-trapping transition. For sufficiently large $t$, the dephasing between different energy eigenstates leads to the broadening of the wave function, and the overlap $|\langle\Omega_{\theta,\varphi}|\psi(t)\rangle|^{2}$ traces out precisely the semiclassical trajectory (black line in Fig. 3).

The paper is organized as follows. In Sec. II, we outline the semi-classical dynamics of Hamiltonian (3) phasespace ; semiclass ; semiclass2 . We analyze the entropy oscillations and the entropy production for short times, and compare the exact dynamics to a semi-classical approximation in Sec. III. In Sec. IV we concentrate on the stationary long time limit of the entropy, and show that it is well approximated by the classical entropy of the microcanonical ensemble. We outline the experimental realization of Hamiltonian (3) in microwave measurements with 87Rb atoms in Sec. V. Our conclusions are summarized in Sec. VI.

II Semiclassical dynamics

In the following sections we will investigate the time evolution of the entropy (4), with different spin coherent states taken as initial conditions. More precisely, we consider spin coherent states lying in the $x-z$ plane, polarized in the direction $\vec{\Omega}_{\theta}=(\sin\theta,\,0,\,\cos\theta)$, for different angles $-\pi/2<\theta<\pi/2$. These states are eigenstates of the spin operator $\hat{S}_{\theta}=\hat{S}_{x}\sin\theta+\hat{S}_{z}\cos\theta$, with eigenvalue $\hat{S}_{\theta}=N/2$. For large total particle number $N$, the semiclassical approximation yields a good description semiclass ; oberthalersemiclass ; semiclass2 , and the spin operators in Eq. (3) can be replaced by the components of a classical vector

[TABLE]

The time evolution of the unit vector $\vec{\Omega}$ is governed by the differential equations semiclass

[TABLE]

with initial condition $\vec{\Omega}(t=0)=\vec{\Omega}_{\theta}$. These classical trajectories lie on the unit sphere, and their shape is determined by the parameter $\alpha$ in Eq. (5) semiclass .

Typical trajectories are depicted in Fig. 4. Below the critical value $\alpha_{c}=1$, all trajectories intersect the equator of the sphere, and no self-trapping occurs. Here the equations of motion (II) have two stable fixed points at $\Omega_{x}=\pm 1$, $\Omega_{y}=\Omega_{z}=0$. For $\alpha>\alpha_{c}$, however, a self-trapped regime appears on the unit sphere (see Fig. 4). Here the fixed point at $\Omega_{x}=-1$ becomes unstable, and bifurcates into two new, stable fixed points at semiclass ; oberthalersemiclass

[TABLE]

Trapped trajectories around these fixed points can not cross the equator of the sphere, but are constrained to the $\Omega_{z}>0$ or $\Omega_{z}<0$ hemisphere. Non-trapped trajectories, however, reach both positive and negative $\Omega_{z}$ values. The separatrix, forming the boundary of self-trapping, touches the equator at the unstable fixed point, $\Omega_{x}=-1$ (see Fig. 4).

The semiclassical trajectories can be visualized more easily by noting that their projections on the $x-y$ plane form circles centered at $(-1/\alpha,0)$ (see Fig. 4),

[TABLE]

The trajectory determined by the initial condition $\vec{\Omega}(t=0)=(\sin\theta,0,\cos\theta)$ will coincide with the separatrix at interaction strength

[TABLE]

Another special case occurs, when the initial condition satisfies

[TABLE]

and $\vec{\Omega}(t=0)$ is a stable fixed point of the classical equations of motion (II).

As we will show later, the long time limit of the entropy (4) reflects this semiclassical dynamics (see Sec. IV). The trapping transition at $\alpha=\alpha_{\theta}$, Eq. (8), is revealed by a sudden jump of size $\log 2$ in the entropy, related to the rapid change by a factor of $2$ in the length of the classical trajectory. The classical fixed point (9) corresponds to a local minimum in $\mathcal{S}$ due to the strong confinement of trajectories around this point.

III Entropy production

Now we concentrate on the entropy production at short times, during the first few oscillations of entropy (4). We consider the spin coherent initial state $|\Omega_{\theta}\rangle\equiv|\hat{S}_{\theta}=N/2\rangle$, with $\hat{S}_{\theta}=\hat{S}_{x}\sin\theta+\hat{S}_{z}\cos\theta$. To gain more insight into the structure of the wave function, let us expand this state in the eigenbasis of $\hat{S}_{z}$ Auerbach ,

[TABLE]

with $|m\rangle$ denoting the eigenstate $\hat{S}_{z}=m$. This expression shows that the shifted spin operator $\hat{S}_{z}+N/2$ follows a binomial distribution $B(n,p)$, with number of trials $n=N$ and probability $p=\cos^{2}(\theta/2)$. This binomial distribution yields an expectation value $\langle\hat{S}_{z}\rangle=np-N/2=N/2\cos\theta$ and a variance ${\rm Var}(\hat{S}_{z})=np(1-p)=N\sin^{2}\theta/4$. In the semiclassical limit of large total particle number $N$, this initial state can be approximated by a Gaussian wave function

[TABLE]

excepting the vicinity of $\theta=0$.

The entropy of a Gaussian distribution is known exactly Renyi , yielding the approximation

[TABLE]

Based on this expression, we scale together the curves corresponding to different total particle numbers $N$ by introducing the rescaled entropy

[TABLE]

The time evolution of the rescaled entropy (11), obtained by exact diagonalization, is shown for different interactions strengths $\alpha$ and particle numbers $N$ in Fig. 5, for two different initial angles $\theta$. Different rows correspond to the non-trapped regime $\alpha<\alpha_{\theta}$ (top), lie on the phase boundary $\alpha=\alpha_{\theta}$ (center), and fall in the self-trapped regime $\alpha>\alpha_{\theta}$ (bottom), respectively, with $\alpha_{\theta}$ given by Eq. (8). As expected, for fixed $\alpha$ but different total particle number $N$, the rescaled entropies $\widetilde{\mathcal{S}}$ follow the same curve for short times. The long time limit of $\widetilde{\mathcal{S}}$, however, is proportional to $\log N$ instead of $\log\sqrt{N}$, and is different for each $N$ (see also Fig. 2). The entropy oscillates both in the non-trapped and self-trapped regimes, while increasing steadily towards a stationary value. These entropy oscillations reveal the coherent oscillations of the single mode condensates, while the dephasing between different energy eigenstates is responsible for the steady increase of the entropy. At the phase boundary, the oscillations become washed out, and the entropy increases approximately linearly, until saturating to the long time limit.

The main features of the time evolution of $\mathcal{S}$ can be understood in terms of the classical trajectories in the semiclassicl limit of large $N$. As supported by the detailed analysis below, the state of the system can be visualized as an extended packet on the unit sphere around the classical unit vector $\vec{\Omega}(t)$. For an initial state $\hat{S}_{\theta}=N/2$, this packet has a Gaussian shape around $\vec{\Omega}_{\theta}$, with variance $\sim 1/N$. The center of the packet, $\vec{\Omega}(t)$, follows the classical equations of motion (II), while the typical width of the packet increases. This broadening occurs because the packet gets more elongated along the classical trajectory due to the dephasing between different energy eigenstates. At the same time, the width perpendicular to the trajectory decreases to keep the volume of the packet constant.

The distribution of $\Omega_{z}$ and the corresponding entropy $\mathcal{S}$ can be determined by projecting this packet to the $z$ axis. The oscillations of the center of the packet, $\vec{\Omega}(t)$, result in entropy oscillations, and their period is given by the period of the classical trajectory. For the separatrix this period is infinity, explaining the vanishing entropy oscillations at the boundary of self-trapping. For a more detailed analysis, notice that the entropy is proportional to $\log\sigma$, with $\sigma$ denoting the typical width of the distribution of $\Omega_{z}$. In spite of the broadening of the packet, $\sigma$ can display a very different behavior depending on the position along the trajectory. At the upper and lower turning points, where the tangent vector of the trajectory is perpendicular to the axis $\hat{z}$, the projection yields a sharp distribution for $\Omega_{z}$, resulting in local minima for the entropy. Since the width of the packet perpendicular to the trajectory decreases, $\sigma$ can even decrease compared to the width of the initial state, resulting in decreasing local minima (see the first row in Fig.5). On the other hand, at the horizontal turning points, where the tangent vector is parallel to $\hat{z}$, $\sigma$ is maximal. This maximal value increases with time as the wave packet gets more elongated along the trajectory, yielding increasing entropy maxima after each oscillation.

To substantiate these arguments and to get a quantitative description for the time evolution of entropy $\widetilde{\mathcal{S}}$, we applied a Gaussian Ansatz for the wave function. This approximation relies on the observation that the initial state is well described by the Gaussian expression (10). We assume that the wave function keeps this Gaussian form during the time evolution. As a first step, we expand the wave function according to the eigenstates of $\hat{S}_{z}$,

[TABLE]

Having separated a rapidly oscillating phase factor $e^{-i\varphi(t)m}$ - corresponding to the rotation of the state around the $z$ axis - we can assume that the coefficients $c_{m}(t)$ are slowly varying functions of $m$.

Let us introduce a new variable $x=2m/N$. In the limit of large $N$, $x$ can be treated as a continuous variable WKB . We can replace the discrete, slowly varying coefficients $c_{m}(t)$ by a continuous function $\psi(x,t)$, and assume a Gaussian form,

[TABLE]

This Ansatz yields a Gaussian distribution for the normalized spin operator $\hat{\Omega}_{z}\equiv 2\hat{S}_{z}/N$, with expectation value $\langle\hat{\Omega}_{z}\rangle(t)=x_{0}(t)$ and variance $1/(N{\rm Re}\,c(t))$. Moreover, for the Gaussian wave function given by Eq. (III), the rescaled entropy (11) can be expressed as

[TABLE]

The optimal parameters of the Gaussian wave function, $|\psi_{G}\rangle$, are determined from the variational condition

[TABLE]

where

[TABLE]

For large total particle number $N$, Eq. (14) can be expanded systematically according to the powers of $N$. The leading order contributions result in the semiclassical equations of motion

[TABLE]

These equations determine the same trajectories as Eqs. (II), with the unit vector $\vec{\Omega}$ given by

[TABLE]

The next order of the expansion yields the time evolution of $c(t)$,

[TABLE]

Notice that $c(t)$ only depends on the dimensionless time $tJ$, the parameter $\alpha$ and the initial condition $\theta$, but not on the particle number $N$. Concentrating on the semiclassical limit of large $N$, we neglect the remaining $O(1/N)$ corrections.

The Gaussian entropy $\eqref{eq:SG}$ is plotted together with exact numerical results in Fig. 5. As noted above, $\mathcal{S}_{G}$ is independent of $N$ up to corrections of the order $1/N$, neglected in our semiclassical approximation. We find that Eq. (13) yields a surprisingly good approximation for the dynamics at short times.

Let us emphasize that the Gaussian ansatz relied on the observation that the spin coherent initial state results in a Gaussian distribution for $\hat{S}_{z}$, thus this description remains valid only on time scales shorter than the time scale of the entropy saturation. Indeed, in Fig. 5 we have found that a Gaussian wave function remains a reasonable approximation on such short time scales. However, as we show in Sec. IV below, this ansatz breaks down as the entropy saturates to the stationary long time limit (see also the middle row of Fig. 5). On such long time scales it has to be replaced by a non-Gaussian semiclassical approximation, derived in Sec. IV and Appendix A.

IV Long time limit of entropy and equilibration

We now turn to the long time behavior of the entropy $\mathcal{S}$, and show how it reflects the semiclassical dynamics discussed in Sec. II. Due to the discrete spectrum of Hamiltonian (3), the entropy shows several revivals and, strictly speaking, it never reaches a stationary value. However, the period of these revivals is typically very long compared to experimentally relevant time scales, and it is still meaningful to consider the steady state at intermediate times revival .

We defined the long time limit of the entropy (4) as the time average

[TABLE]

with $T$ chosen large enough to reach a stationary value. As before, we used spin coherent initial states $|\hat{S}_{\theta}=N/2\rangle$. The numerical results from exact diagonalization are shown in Fiq. 6 as a function of the angle $\theta$, for two different parameters $\alpha$.

The semiclassical fixed point, Eq. (9), appears as a sharp local minimum in the time averaged entropy. Since the entropy is related to the width of the distribution of $\hat{S}_{z}$, this entropy minimum follows from the strong confinement of classical trajectories around the stable fixed point, leading to sharp distributions for $\Omega_{z}$. The separatrix of the self-trapped phase, Eq. (8), is revealed by a sudden jump of size $\log 2$ in entropy (16). This sudden entropy gain is related to the doubling of the length of classical trajectories at the self-trapping transition, doubling the phase space available for $\vec{\Omega}$ (see Fig. 4).

The long time limit of the entanglement entropy, Eq. (16), can be understood in terms of a semiclassical microcanonical ensemble. In the semiclassical approximation, the trajectory is determined by Eqs. (II), which conserves the energy of the classical Hamiltonian. In a microcanonical description, the spin vector $\vec{\Omega}$ is randomly distributed along this trajectory (the surface of constant energy in general), resulting in a uniform distribution on the classical trajectory. This classical trajectory amounts in a continuous distribution for the $z$-component of the spin, $\Omega_{z}$ (see Appendix A). Denoting the corresponding probability density by $P(\Omega_{z})$, the classical entropy is given by

[TABLE]

Here $\mathcal{S}_{0}$ denotes an arbitrary constant entropy shift, accounting for some unknown box size $\Delta\Omega_{z}$. Eq. (17) yields good agreement with the numerical results by using a single fitting parameter $\mathcal{S}_{0}=7.0$ (see Fig. 6). The deviation between the semiclassical approximation and $\overline{\mathcal{S}}$ gets larger only in the immediate vicinity of the semiclassical fixed point, Eq. (9). Here $\mathcal{S}_{\rm class}$ diverges, because the variance of the continuous classical distribution $P(\Omega_{z})$ approaches zero. However, the time averaged entropy $\overline{\mathcal{S}}$ remains non-negative even at the fixed point, and its minimal value is determined by the width of the spin coherent initial state in the eigenbasis of $\hat{S}_{z}$. For the semiclassical case of large $N$ this yields a minimal variance ${\rm Var}(\hat{S}_{z})\approx N/(4\,\alpha^{2})$ and minimal entropy

[TABLE]

A classical microcanonical equilibrium ensemble thus yields a good approximation for the time averaged entropy, in spite of the unitary time evolution of the system. To see how dephasing alone can lead to equilibration, we plotted the overlap between the wave function $|\psi(t)\rangle$ and the spin coherent states polarized into direction $(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta)$, $|\Omega_{\theta,\varphi}\rangle$, for long times in Fig. 7. For the initial state we chose the maximally polarized state $\hat{S}_{z}=N/2$, and selected two different interaction strengths, corresponding the non-trapped and self-trapped regimes, respectively. In accordance with the classical picture presented above, the wave function broadens, and extends to the vicinity of the classical trajectory on both sides of the self-trapping transition. Note that in spite of the success of this semiclassical description, Fig. 7 still reflects the non-classical nature of the exact wave function, by displaying several very sharp minima on the unit sphere, originating from quantum interference. Notice also that the state of the system never becomes stationary, and at any time the density displays several maxima along the trajectory. The position of these maxima depends on time, leading to revival effects. However, the time averaged entropy is always well approximated by a state spreaded uniformly along the classical trajectory.

As a final remark, let us note that in order to observe self-trapping, the upper bound of the time averaging Eq. (16), $T$, should not be too large. Strictly speaking, a real self-trapping transition only occurs in the classical system. The quantum system would eventually tunnel to reversed population imbalances tunneling , and the entropy would increase by $\log 2$, compensating the entropy loss of self-trapped regime. We have chosen $T$ much smaller than the time scale of this tunneling process, $T_{\rm tunnel}$. The period $T_{\rm tunnel}$ depends on the small energy difference $\Delta E_{\rm doublet}$ between quasi-degenerate doublets $(|u\rangle\pm|d\rangle)/\sqrt{2}$, with the states $|u\rangle$ and $|d\rangle$ confined to the $\hat{S}_{z}>0$ and $\hat{S}_{z}<0$ hemispheres, respectively. The splitting $\Delta E_{\rm doublet}$ is tiny even for moderate particle numbers  footnote2 , resulting in an exponentially long tunneling time $T_{\rm tunnel}$, far beyond the time scales available for experiments or numerical simulations. However, for small systems consisting only of a few particles, this tunneling process can be observed in simulations. An example of this population inversion is shown in Fig. 8, for a small particle number $N=10$  footnote4 .

V Microwave experiments

As mentioned in the introduction, Hamiltonian (3) can also be realized in microwave experiments with ultracold atoms, making use of different hyperfine states of the atoms. Below we outline how Hamiltonian (3) arises in this setup, and we briefly discuss the optimal experimental parameters for the observation of entropy oscillations for the specific case of 87Rb atoms.

In 87Rb experiments, one can tune interactions by utilizing the two hyperfine states $|0\rangle\equiv|F=1,\,m_{F}=1\rangle$ and $|1\rangle\equiv|F=2,\,m_{F}=-1\rangle$, with $F$ and $m_{F}$ denoting the total spin of the atom and its projection to the quantization axis, respectively feshbachrev . These states can be trapped in optical dipolar traps, while they can be coupled by microwave pulses. Below we concentrate on an experimentally relevant setting, where these atoms are trapped in a spherically symmetric harmonic trap with trapping frequency $\nu_{0}=50$Hz. For weak enough interactions, all atoms occupy the ground state of this harmonic potential,

[TABLE]

with $l_{0}=\sqrt{\hbar/(2\pi m\nu_{0})}$ denoting the oscillator length. Let us note that for stronger interactions the condensate wave function is better described by a Thomas-Fermi profile instead of $\varphi_{0}$, Eq. (18). However, for not too large interaction strengths, a two mode approximation is still applicable, and the derivation presented below is valid with minor modifications (see also the discussion at the end of the section).

The short range interaction between the Rb atoms is well described by a Dirac-delta potential. Denoting the bosonic creation operators of the hyperfine states $|0\rangle$ and $|1\rangle$ by $\hat{a}_{0}^{\dagger}$ and $\hat{a}_{1}^{\dagger}$, the interaction energy is given by

[TABLE]

Here $U_{\sigma\sigma^{\prime}}=g_{\sigma\sigma^{\prime}}/(2\pi l_{0}^{2})^{3/2}$, and the interaction strength $g_{\sigma\sigma^{\prime}}$ can be expressed with the scattering length of the Rb atoms, $a_{\sigma\sigma^{\prime}}$, as ZwergerReview

[TABLE]

with $m$ denoting the mass of 87Rb.

The bare scattering lengths of 87Rb depend very weakly on the hyperfine states of the atoms, and all interactions are determined by the single length scale $a_{\sigma\sigma^{\prime}}=5.3$nm. The scattering length between hyperfine states $|0\rangle$ and $|1\rangle$, however, can be tuned by a Feshbach resonance, changing $a_{01}$ by as much as $\Delta a_{01}=0.1\,a_{00}=0.53$nm feshbach . Introducing the average interaction strength $\overline{U}=(U_{00}+U_{01})/2$, and the difference $\Delta U=U_{00}-U_{01}$, the interaction energy (V) can be rewritten as

[TABLE]

Here $\hat{N}_{i}=\hat{a}_{i}^{\dagger}\hat{a}_{i}$ for $i=0,1$, $N=\hat{N}_{0}+\hat{N}_{1}$ is the total particle number, and the spin operator is defined as $\hat{S}_{z}=(\hat{N}_{0}-\hat{N}_{1})/2$. For a closed system $\overline{U}N^{2}/2$ is just an irrelevant constant energy shift. Thus the interaction between Rb atoms takes the same form as the interaction term in Hamiltonian (3), with interaction strength $\Delta U$ determined by the difference of scattering lengths $a_{00}-a_{01}$.

Let us mention that instead of controlling the scattering lengths $a_{\sigma\sigma^{\prime}}$ by a Feshbach resonance, the interaction strength $\Delta U$ can also be tuned by applying a microwave trapping potential which depends on the hyperfine state of the atoms riedelentanglement . In this case the atoms occupy state-dependent condensate wave functions, $\varphi_{\sigma}$, and the interaction strength is given by

[TABLE]

Thus $\Delta U=U_{00}-U_{01}$ can be controlled by changing the overlap between the two condensate modes $\varphi_{0}$ and $\varphi_{1}$.

The hyperfine states $|0\rangle$ and $|1\rangle$ can be coupled by a two photon transition, where a detuned microwave pulse couples $|0\rangle$ to an intermediate state $|F=2,m_{F}=0\rangle$, coupled to the final state $|1\rangle$ by a radiofrequency transition hall . This two photon transition gives rise to a hopping term $-J(\hat{a}_{0}^{\dagger}\hat{a}_{1}+\hat{a}_{1}^{\dagger}\hat{a}_{0})$ in the Hamiltonian. In the spin representation this corresponds to the term $-2J\hat{S}_{x}$, thus with the already known form of the interaction, $\Delta U\hat{S}_{z}^{2}$, we recover Hamiltonian (3).

To reach optimal parameters, one needs strong enough interactions and therefore relatively strong confinement. For a trap frequency $\nu_{0}=50{\rm Hz}$, and typical scattering length difference $a_{00}-a_{01}=0.1\,a_{00}=0.53$nm, the interaction strength is $\Delta U/h=0.014$Hz. For atom numbers in the range of $N\sim 2000$, the relevant parameter of the spin model is around $\Delta UN/h=30$Hz. For the typical entropy oscillations plotted in Fig. 2, the parameter $\alpha=\Delta UN/(2J)$ is roughly $\alpha\sim 2$, corresponding to $J/h=7$Hz. With these parameters, the typical time scale of entropy oscillations and entropy generation for the maximally polarized initial state $|\hat{S}_{z}=N/2\rangle$ is expected to be around $t\sim 70$ms,which is much shorter than the lifetime of a condensate, and also much shorter than the coherence time of superposition states microwave ; treutlein ; deutsch2 ; opticaltrap . Therefore the entropy oscillations should be observable on experimentally realistic time scales.

Let us note that for the parameters above $\overline{U}\gg\Delta U$ implies $N\overline{U}/h\gg\nu_{0}$. Since the typical scale of interaction energy is much larger than the trapping frequency $\nu_{0}$, the atoms do not remain in the ground state of the harmonic potential. However, the system can still be described as two coupled single mode condensates with a modified condensate wave function $\varphi_{0}$, because $N\Delta U/J<\nu_{0}$ footnote2 , and the entropy oscillations and entropy generation remain observable with a slightly modified oscillation frequency and entropy production rate.

VI Conclusions

In this work we analyzed the entropy generation for two coupled single-mode Bose-Einstein condensates, realized by loading a condensate into a double well potential, or by an interstate coupling of a two component Bose-Einstein condensate. This system provides one of the simplest examples to study the entropy production by the coherent time evolution of coupled quantum systems. Even though entanglement measurements in generic correlated many-body systems are challenging, in this setting the entanglement between the two potential wells should be experimentally accessible by measuring the number of atoms in the wells. Besides its experimental relevance, the dynamics of coupled single-mode condensates already shows interesting physics. At large particle number imbalances and sufficiently strong interactions, the system enters a self-trapped regime, where the amplitude of population imbalance oscillations gets suppressed due to the large interaction energy of the initial state.

Concentrating on the entropy production during unitary time evolution, we investigated the time dependence of the entropy at $T=0$ temperature, by combining numerical results with analytical calculations. We found that the coherent oscillations of the single mode condensates manifest in entropy oscillations on the top of a steady entropy generation. These coherent oscillations only vanish in the vicinity of the self-trapping transition, where the entropy increases linearly for short time scales. In this pure quantum state, the entropy production originates from the dephasing between different energy eigenstates, eventually leading to a stationary, saturated entropy. Interestingly, this entropy saturation looks like equilibration, in spite of the coherent time evolution of this closed system. These results should be experimentally observable for realistic parameters in microwave measurements with 87Rb atoms. Here the two modes of the condensate are not spatially separated; instead they correspond to two different hyperfine states of 87Rb.

To gain more insight into the entropy oscillations and entropy production on short times scales, we have shown that the time evolution of the entropy can be understood in terms of the semiclassical trajectories of the system. The wave function can be visualized as a broadening packet on the unit sphere, with its center evolving along the classical trajectory. To obtain a quantitative description, we have shown that a Gaussian Ansatz for the wave function, together with a semiclassical expansion, yields a surprisingly good approximation for the exact time evolution.

We also analyzed how the stationary long time limit of the entropy reflects the semiclassical dynamics of the system. The classical fixed point is revealed by a local minimum in the entropy, related to the strong confinement of trajectories in the vicinity of this point, while the self-trapping transition is accompanied by a sudden entropy jump of size $\log 2$, due to the rapid change by a factor of 2 in the length of trajectories.

In order to investigate the dephasing induced equilibration of the entropy in more detail, we compared the numerical results to the prediction of a classical microcanonical ensemble, where the spin vector is distributed uniformly over the classical trajectory. We found that this ensemble yields a surprisingly accurate description for the stationary limit of the entropy. To gain more insight into the exact time evolution of the system, we calculated the overlap of the wave function with the spin coherent states of different orientations. We have shown that this overlap traces out the classical trajectories on the unit sphere at long times, supporting our picture describing the system in terms of a classical microcanonical ensemble.

In this work we concentrated on the entropy generation in a pure state, and performed all calculations at $T=0$ temperature. Thermal fluctuations are expected to shift the entropy of the system to higher values, while the reduction of coherence starts to wash out the oscillations during the entropy production. However, the entropy oscillations should still remain visible for low enough temperatures of the order of a few hundred nK.

Since entropy generation lies at the heart of equilibration and thermalization in closed systems, the detailed analysis of entropy production in other correlated many-body systems, and entanglement spreading in the presence of conserved quantities, remains a question of fundamental interest.

Acknowledgements.

We acknowledge fruitful discussions with Péter Domokos. This research has been supported by the National Research, Development and Innovation Office - NKFIH Nos. K105149 and SNN118028. ED acknowledges support from Harvard-MIT CUA, NSF Grant No. DMR-1308435, AFOSR Quantum Simulation MURI, AFOSR Grant No. FA9550-16-1-0323.

Appendix A Classical microcanonical distribution of $\Omega_{z}$

In this appendix we outline the calculation of the continuous distribution $P(\Omega_{z})$ used in Eq. (17). In a microcanonical description, we assume that the unit vector $\vec{\Omega}$ is distributed uniformly along the curve of constant energy, selected by the initial state $\vec{\Omega}=(\sin\theta,\,0,\,\cos\theta)$. The projection of this classical trajectory to the $x-y$ plane is a circle, given by

[TABLE]

with $r_{\theta}=|\sin\theta+1/\alpha|$. Based on this expression, the trajectory can be parametrized by an angle $\chi$ as

[TABLE]

The length of a small arc segment, $\Delta s$, can be expressed with the parameter change $\Delta\chi$ as

[TABLE]

with $|\partial_{\chi}\vec{\Omega}(\chi)|$ denoting the length of the tangent vector $\partial_{\chi}\vec{\Omega}(\chi)$. For a uniform distribution along the curve, $\vec{\Omega}$ points to this segment with probability $\Delta s/s_{0}$, where $s_{0}$ is the total length of the trajectory. This results in the following distribution $P(\chi)$ for parameter $\chi$

[TABLE]

The distribution of $\Omega_{z}$ can be expressed as

[TABLE]

where $\chi(\Omega_{z})$ denotes the inverse function of $\Omega_{z}(\chi)$, and the factor 2 arises from the symmetry $\Omega_{z}(\chi)=\Omega_{z}(-\chi)$.

The distribution $P(\Omega_{z})$ as a function of parameter $\alpha$ and initial condition $\theta$ is given by

[TABLE]

with sgn denoting the sign function. The length of the trajectory, $s_{0}$, ensures the correct normalization of the distribution,

[TABLE]

The support of this distribution $(\Omega_{z}^{\rm min},\,\Omega_{z}^{\rm max})$, however, is different on the two sides of the self-trapping transition. In the non-trapped regime

[TABLE]

while for a self-trapped state

[TABLE]