Dynamics of entanglement entropy of interacting fermions in a 1D driven harmonic trap

Following up on a recent analysis of two cold atoms in a time-dependent harmonic trap in one dimension, we explore the entanglement entropy of two and three fermions in the same situation when driven through a parametric resonance. We find that the presence of such a resonance in the two-particle system leaves a clear imprint on the entanglement entropy. We show how the signal is modified by attractive and repulsive contact interactions, and how it remains present for the three-particle system. Additionaly, we extend the work of recent experiments to demonstrate how restricting observation to a limited subsystem gives rise to locally thermal behavior.


Introduction
The link between the unitary evolution of microscopic quantum systems and the ergodic dynamics of macroscopic classical systems has been under investigation since the earliest days of quantum statistical mechanics. The fundamental question can be stated as, "How does a quantum state, which remains pure under action of the time-evolution operator, give rise to thermal behavior?" A prominent theory to address this issue is the eigenstate thermalization hypothesis, which is closely related to quantum chaos and entanglement (see e.g. [1][2][3]. More recently, experiments have directly shown local thermal behavior in globally pure states [4].
In this work, we investigate local thermalization by exploring the time-dependent dynamics of entanglement in small systems of spin-1/2 particles in a one-dimensional harmonic trap. Specifically, we consider the two-body problem with one particle per species, which we will call the 1+1 case, and the three-body problem of two particles of one species and one of the other, i.e. 2+1 particles. We focus on a situation recently analyzed in Ref. [5], where it was found that the interacting and noninteracting two-body problems can display dramatically different behaviors when driven through a parametric resonance. Specifically, we consider dynamics governed by a Hamiltonian where the kinetic and interaction energy operators are time-independent and given bŷ Speaker, e-mail: joshmck@unc.edu respectively, whereas the external potential iŝ where ω 2 (t) = ω 2 0 (1 + α sin (t/T )), for a fixed time-independent frequency ω 0 and parameters α, T to be explored; note α is dimensionless, but T has units of time, i.e. inverse energy. This form was chosen because it follows the protocol of Ref. [5]. As usual, the operatorsψ † s ,ψ s used above are the creation and annihilation operators for particles of spin s =↑, ↓ andn s is the corresponding particle density operator. The coupling g is directly related to the scattering length a 0 via g = 2/a 0 , and we set m = = 1.
The work of Ref. [5] compared the behavior of different cases with g ≤ 0, i.e. either noninteracting or repulsively interacting (note our sign convention). Here, we reproduce part of that analysis as a way to check our calculations and include the attractively interacting case as well. We further extend that work by studying the behavior of the Rényi entanglement entropy (defined next) as a function of subsystem size and for different coupling strengths and driving parameters. To provide a more intuitive measure of thermal behavior, we present our results as the fraction of the Hilbert space required to describe the system.

Entropy
Since the seminal work of Shannon in the 1940s [6], entropy has been widely used as a precise means of quantifying information. Defined for a discrete random variable X with possible values x, the Shannon entropy S S takes the form where the summation becomes an integral for continuous X. The essential question that this quantity answers is roughly, "How much additional information is needed to unambiguously specify a particular configuration?" Thus for states that are fully known, the entropy is zero, while for states with high multiplicity, the entropy is large. In this sense, the Shannon entropy parallels the combinatorial entropy of the microcanonical ensemble, with its macrostates and microstates. For uniform probability distributions p(x), the Shannon entropy reduces to the Boltzmann entropy; on the other hand, if p(x) is defined by Boltzmann weights and the partition function, the Shannon entropy gives the same result as the canonical ensemble for the entropy of the ideal gas (under the usual assumptions of indistinguishability and phase space in units of Planck's constant h).
The quantum mechanical analog of the Shannon entropy is the von Neumann entropy, whereρ is the density operator. The "variable" now under consideration is the space of states accessible to a quantum system, and the missing information quantified by this entropy is the amount needed to distinguish which particular state the system is in. The von Neumann entropy may be considered as a special case (n → 1) of the more general n-th order Rényi entropies, Of particular physical interest is the second order Rényi entropy, as Trρ 2 reflects the purity of a quantum state, equal to 1 for a pure state and less than 1 otherwise. Correspondingly, S 2 grows with increased mixing of states. For a maximally mixed state, where V H is the volume of the Hilbert space. Accordingly, we identify exp(S 2 ) as the weighted number of states contributing toρ and introduce the fraction of the Hilbert space occupied, exp(S 2 )/V H , as a normalized measure valid for any state. When this quantity approaches unity, all states become equally likely to be occupied, in analogy to the infinite-temperature limit.

Entropy of entanglement and reduced density matrices
To characterize spatial entanglement, we divide the 1D region on which the full Hilbert space H exists into two subregions, A and its complement,Ā. The full space H may be written as a direct product space of the Hilbert spaces of the subregions as Since the subsystem A is our focus and its complementĀ is the "environment," we have access only to the density matrix of A,ρ A , which is the reduced density matrix of the full system, where the trace is over the states of HĀ in the position-space basis. Onceρ A is known, the Rényi and von Neumann entanglement entropies can be computed via and Both entropies naturally vanish when A is an empty set or the whole system, as thenρ A corresponds to a single state (which is trivially factorizable).

Computational method, lattice, scales
To study the small systems we are concerned with here, we employ the non-stochastic method of Ref. [7] to build the transfer matrixT = exp −τĤ .
In Ref. [7], where no external potential was present, the transfer matrix was approximated for small timestep τ using a symmetric Trotter-Suzuki decomposition aŝ In the present work, we extend the above to includeV ext by settinĝ  With this new transfer matrix, we evolve a trial state in imaginary time to project onto the ground state of the system by maintaining a constantV ext . From that point on, we implement a Wick rotation and set τ → iδt to study the real-time dynamics of the ground state under time variations inV ext (as specified in the Introduction; see below for more details).
In our calculations, we used a discretized version of the Hamiltonian of Eq. (1) with the kinetic energy piece implemented in momentum space and the interaction and external potential pieces applied in coordinate space. Our calculations were carried out in a spatial lattice of N x points where the lattice spacing was defined as = 1 to set the scale for every other quantity. The evolution in imaginary time was performed with a lattice spacing τ = 0.005 2 , whereas the real-time evolution used a smaller spacing δt = 0.001 2 . The dimensionful quantities defining the Hamiltonian of Eq. (1) are g, ω 0 and T . From those, we form two dimensionless parameters g 2 /ω 0 and ωT , which are held constant to their desired physical values when taking the continuum limit. To take that limit, we must ensure the following ordering of scales to ensure that the trapping potential is insensitive to the lattice spacing and volume. In practice, we reduced the numerical value of ω 0 until no differences were observed. To keep the physics constant, this required reducing g 2 and increasing T . This is consistent with the idea that we must have g 2 1/ 2 to ensure that the internal length scales of the two-body interaction are insensitive to the lattice spacing. Below, whenever a specific value of g is quoted, it should be understood as corresponding to ω 0 = 1. We used lattice sizes N x = 20, 40 and found the finite-volume effects on the (1+1)-and (2+1)-particle systems to be negligible.
It is worth noting that for the onedimensional systems considered here, other methods such as conventional exact diagonalization are available. However, the approach presented here does not rely on the dense-matrix form of the quantum operators, nor on any particular symmetries of the system, and therefore it is applicable "as is" to the two-and three-dimensional counterpart systems, and in principle to higher particle number as well (memory constraints permitting; see also our comment below on attempting the 3-body problem).

Crosschecks with results of Ebert et al.
As a cross-check for correctness of our lattice approach, we reproduced the results of Ref. [5], comparing the two methods for the following quantities: (1) the root-mean-square separation of particles under an exponential driving potential with finite cutoff, which show excellent agreement with Ref. [5], as shown at left in Fig. 1; (2) the energy in a harmonic trap external potential with periodic frequency changes for interacting and noninteracting cases (right of Fig. 1), which agrees with the corresponding plot of Ref. [5] once the center-of-mass energy is taken into account; (3) the stability of the system under variation of the harmonic trap parameters α and T (Fig. 2). As our approach did not involve an analytical solution, which would clearly indicate regions of instability, we used the maximum energy E max /E 0 (average energy for the attractive case) attained over a trajectory to provide a numerical indicator for resonant behavior. Our results are quantitatively consistent with those of Ref. [5].

Local thermal behavior
The trajectory represented in Fig. 3 is the same as that in Fig. 1 (right). As described in Ref. [5], the driving frequency is tuned to the gap between the ground and second excited states, such that at t ≈ 225 and t ≈ 700, the system is nearly entirely in the excited state, while at t = 0 and t ≈ 480, it is in the ground state. Viewed from top to bottom, it is evident that the fraction of the Hilbert space required to described the state decreases as the subsystem size increases, approaching 1/V H A as L A approaches the full system size, L. As the system is excited, the smallest subsystem's (topmost curve, blue) Hilbert space occupation decreases because particles are driven away from the center of the trap. For larger subsystems, however, the opposite occurs because prior to the excitation, the particles were mostly contained within the boundaries of L A . Thus, we see that the size of the subsystem observed relative to its exterior can cause a large variation in the description of the quantum state. Additionally, we see that the local basis of very small subsystems is inefficient at describing states in the larger basis of the full system, leading to apparent thermal behavior in the subsystem. While in agreement with Ref. [4], we extend that work by observing the system as a function of time as the external driving continues.

Entanglement entropy
In this section we present the calculation of the entanglement entropy across a parametric resonance, focusing for definiteness on the case α = 0.5 (see Fig. 2). In the absence of interactions, the trap oscillations drive the system across a parametric resonance which, for α = 0.5, is present around T 0.5, 1.0, and 1.5 (and traces of another one can be seen in Fig. 2 for T 2.1). In Fig. 4 we show the impact of the driving potential on the time evolution of the second Rényi entanglement entropy S 2 . While it is a non-trivial problem to understand the behavior of S 2 analytically (because it is a complicated quantity that varies with the subsystem size L A ), certain features can be gleaned which, as we will see below, are a clear imprints of the interactions, while other features are essentially unchanged relative to the noninteracting case. Specifically, we see that the resonance yields a seemingly chaotic signal when plotted in the t, T plane (to be contrasted with the more regular structure observed away from the resonance in those plots) and changes the entanglement entropy by a large factor. In Fig. 4 we show the second Rényi entanglement entropy S 2 plotted as exp(S 2 )/V H A , where V H A is the number of states in the Hilbert space where S 2 is computed (which depends on the subsystem size L A ). Since this work started, we have noticed that the time-evolution operator did not remain exactly unitary for interacting (2+1)particle systems. While this is likely a numerical issue, we have yet to determine the exact cause. This non-conservation of probability may very well account for the unexpected increase in entropy seen in Fig. 4 over time, as the evolution becomes essentially random.

Summary and Conclusions
In this work we studied the time evolution of systems of 1+1 and 2+1 spin-1/2 particles in a onedimensional time-varying external trapping potential. Expanding on the work of Ref. [5], which we reproduced as a crosscheck for the (1+1)-particle case, we determined the impact on the entanglement entropy of driving the system through a resonance. It is worth noting that we used very different techniques from Ref. [5]: we put the particles on a spatial lattice and carried out the time evolution by first projecting onto the ground state (thus preparing the system) using imaginary-time evolution (i.e. applying the transfer matrix), followed by real-time evolution to study the dynamics.
The objective of Ref. [5] was to investigate the interaction effects on the two-body problem by considering two different scenarios: a rapid trap change and a periodic trap change. While the former behaved essentially as the non-interacting case, the latter, which is the main one we considered here, can enhance small fluctuations by virtue of parametric resonance. Such a situation is interesting for several reasons: small interaction effects may have a large impact on the overall behavior of the system; the driving protocol is similar to the quantum-control protocol used in experiments; and resonances can be used for heating or cooling. To add insight into this problem, and inspired by the work of Ref. [4], we examined the emergence of local thermal behavior and showed how it varied over time under the influence of external driving. To this end, we put forward a non-stochastic lattice technique that only incurs uncertainties in lattice spacing. As expected, the Rényi entanglement entropy showed thermal behavior diminishing as the subsystem size grew to the full system. This work is therefore a very first step towards building a more robust technique that can address larger few-body problems (current memory limitations should allow for at least 10 particles in three spatial dimensions, but using supercomputers and aggressive optimization) and a much wider range of observables and time-varying situations.
This material is based upon work supported by the National Science Foundation under Grant No. PHY1452635 (Computational Physics Program).