Table Of ContentFlux quench in a system of interacting spinless fermions in one dimension
Yuya O. Nakagawa,1,∗ Gr´egoire Misguich,2 and Masaki Oshikawa1
1Institute for Solid State Physics, the University of Tokyo,
5-1-5 Kashiwanoha, Kashiwa, Chiba Japan 277-8581
2Institut de Physique Th´eorique, Universit´e Paris Saclay, CEA, CNRS, F-91191 Gif-sur-Yvette, France
(Dated: May 30, 2016)
Westudyaquantumquenchinaone-dimensionalspinlessfermionmodel(equivalenttotheXXZ
spin chain), where a magnetic flux is suddenly switched off. This quench is equivalent to imposing
a pulse of electric field and therefore generates an initial particle current. This current is not a
conservedquantityinpresenceofalatticeandinteractionsandweinvestigatenumericallyitstime-
evolutionafterthequench,usingtheinfinitetime-evolvingblockdecimationmethod. Forrepulsive
interactionsorlargeinitialflux,wefindoscillationsthataregovernedbyexcitationsdeepinsidethe
Fermi sea. At long times we observe that the current remains non-vanishing in the gapless cases,
6
whereas it decays to zero in the gapped cases. Although the linear response theory (valid for a
1
weak flux) predicts the same long-time limit of the current for repulsive and attractive interactions
0
(relationwiththezero-temperatureDrudeweight), largernonlinearitiesareobservedinthecaseof
2
repulsive interactions compared with that of the attractive case.
y
a PACSnumbers: 05.30.-d,05.70.Ln,67.10.Jn,75.10.Jm
M
7 I. INTRODUCTION andlatticeeffectsthatisresponsibleforthenontrivialdy-
2 namics. Second, changing the flux by an integer number
of flux quanta on a periodic chain amounts to a unitary
Thenonequilibriumdynamicsofisolatedquantumsys-
] transformation of the Hamiltonian and therefore leaves
h tems has become a major subject of study in condensed
c matter physics1. Thanks to substantial developments on the energy spectrum unchanged. In other words, for an
e integernumberoffluxquanta,theenergyspectrumisthe
the experimental side, it is now possible to compare the-
m same for the pre-quench and post-quench Hamiltonians.
oretical predictions and experiments with a high accu-
- In that case the dynamics solely comes from a change
racy and controllability, in particular in the field of cold
t
a atoms. Quantum quench, a sudden change of some pa- in the eigenstates, not from their energies. Third, this
t quenchallowstomakecontactwithsometransportprop-
s rameter(s) in a quantum system, is one of the simplest
. ertiesofthesystem. Whenthenumberoffluxquantaper
t protocols to drive systems out of equilibrium. Typically,
a the length of the system is small, the electric field pulse
an initial state is prepared as the ground state of some
m is weak. In this limit, the dynamics may be described
pre-quench Hamiltonian and some external parameter is
- then abruptly changed at t=0. This leads, for t>0, to by using the linear response theory; long-time limit of
d the current should then be directly proportional to the
n a unitary evolution with a different Hamiltonian and to
zero-temperature Drude weight of the model2,3.
o some non-trivial dynamics.
c Quantum quenches in one-dimensional (1d) systems In this study we use the infinite time-evolving block
[ decimation (iTEBD)4 method to monitor the evolution
havealreadybeenintensivelystudiedforseveralreasons.
of the wave function. We focus on the particle current
2 First,theeffectofinteractionsandquantumfluctuations
v are particularly important in 1d. Second, several pow- and analyze its dynamics, including its long-time limit.
7 erful analytical and numerical methods, such as Bethe As an important result, we observe some current oscil-
6 lations at intermediate times. In addition, these current
Ansatz, bosonization, Time-Evolving Block Decimation
1 oscillationsarefoundtobecarriedbyexcitationslocated
(TEBD),andtime-dependentDensityMatrixRenormal-
6
deepinsidetheFermisea. Finally, wefindthatthelong-
0 izationGroup(t-DMRG),areavailableforthesesystems;
time limit of the current depends in a nonlinear way on
. these methods allow us to make predictions concerning
1 the initial flux. These nonlinearities appear to be par-
the dynamics of these quantum many-body systems. In
0
ticularly strong in the case of repulsive interactions be-
the present study we consider a simple quench for inter-
6
tweentheparticles. Atheoreticalunderstandingofthese
1 actingspinlessfermionsin1d,whereanAharonov-Bohm
observations – presently lacking – would require us to
: fluxissuddenlyswitchedoff. Thisisequivalenttoanap-
v go beyond an effective low-energy description, such as
plicationofaninstantaneouspulseofelectricfield,which
i
X generatesaninitialparticlecurrent. Thisquenchhassev- bosonization. Severalotherquantities,likethegrowthof
the entanglement entropy, are also computed.
r eral appealing properties. First, the non-trivial dynam-
a
icscomesspecificallyfromlatticeandinteractionseffects. The remainder of the paper is organized as follows.
Indeed, the current is an exactly conserved quantity for In Sec. II, we introduce the model and the flux quench
free fermions on a lattice, as well as for any model with problem. In Sec. III, we review the numerical method
translation symmetry in the continuum (due to Galilean (iTEBD) and present our numerical results for the dy-
invariance)27. So, it is the combination of interactions namics after the quench. In Sec. IV, we summarize our
2
results and state conclusions. Technical details on nu- of motion in presence of interaction ∆ (cid:54)= 0. The latter,
merical calculations are presented in the Appendix. combined with the presence of lattice, causes Umklapp
scattering as
II. FLUX QUENCH [H ,Jˆ]=− ∆ (cid:88)(cid:16)c†c +c† c (cid:17)(n −n ),
0 2iN i i+1 i+1 i i−1 i+2
i
We consider one of the simplest interacting spinless (3)
fermion systems in one dimension (assuming a periodic where H is the Hamiltonian without flux, H = H(t ≥
0 0
boundary condition): 0). As mentioned in the introduction, the interactions
are essential to produce a nontrivial dynamics. In this
N−1 N−1
1 (cid:88) (cid:16) (cid:17) (cid:88) study we focus on the expectation value of the cur-
H(t)=− e−iθ(t)c†c +h.c. −∆ n˜ n˜ ,
2 i i+1 i i+1 rent, J(t) := (cid:104)ψ(t)|Jˆ|ψ(t)(cid:105). In particular we analyze (i)
i=0 i=0
the time-evolution towards stationary states and (ii) the
(1)
long-time limit of the current. We note that Mierzejeski
where n˜ =c†c −1/2, and N is a total number of sites.
i i i et al.3 recently utilized this flux quench to illustrate the
As is well known, a Jordan-Wigner transformation maps
breakdownofthegeneralizedGibbsensemble6. Also,the
this model to a spin-1/2 XXZ chain5. We focus here
Loschmidtechoassociatedtothisquenchwasconsidered
on the zero chemical potential case, which corresponds
using the Bethe ansatz in Ref. 7, and the flux quench for
tozeroexternalmagneticfieldinthespinlanguage. The
bosons was studied in Refs. 8,9.
phasefactorθ(t)inthehoppingtermsisthevectorpoten-
When θ is small, we expect that the linear response
tial, representing an Aharonov-Bohm flux Φ(t) = Nθ(t) 0
(LR) theory can be applied to obtain J(t) as a response
piercing the ring. In the spin language, it introduces a
to the weak electric pulse as
twist in the xy-plane. In the following, we call θ(t) flux
strength. 1 (cid:90) ∞
InwhatfollowswetakethethermodynamiclimitN → J(t)= 2π dωF(ω)σ(ω)e−iωt+O(θ02), (4)
∞ while keeping the flux strength θ(t) constant. The −∞
model is integrable for any value of ∆, and the phase di- whereF(ω)=θ istheFouriertransformoftheimposed
0
agramhasbeenthoroughlystudied5. For∆>1(∆=1), electric field E(t)=θ δ(t), and σ(ω) is the conductivity
0
the system is gapped (gapless), and there are two de-
generate ground states which are exactly given by the N (cid:90) ∞
σ(ω)= dteiωt(cid:104)[Jˆ(t),Jˆ](cid:105) . (5)
completely empty state and the completely filled state. ω GS,0
0
They correspond to ferromagnetic states in the spin lan-
guage. In these cases, the ground states are completely Here (cid:104)...(cid:105) denotes the expectation in the ground
GS,0
insensitive to the flux and there will be no dynamics as state of H . The conductivity has a zero-frequency com-
0
well. For −1≤∆<1, the system is gapless and its low- ponent,calledtheDrudeweight,aswellasaregularpart:
energyuniversalbehaviorsaredescribedbybosonization
as Tomonaga-Luttinger liquid. For ∆ < −1, the system σ(ω)=2πDδ(ω)+σreg(ω). (6)
is again gapped, and there are two degenerate ground
In the case of flux quench, Eq. (4) gives
states corresponding to an antiferromagnetic long-range
order in the spin language. In contrast to the ferromag- θ (cid:90) ∞
netic case, however, the ground states are still nontrivial J(t)= 0 dωσ(ω)e−iωt+O(θ2) (7)
2π 0
and sensitive to the flux. Therefore, in this paper, we −∞
consider the regime ∆ < 1, where there are nontrivial
and,asnotedinRef.3,thelong-timelimitofthecurrent
effects of the flux.
is proportional to the Drude weight D
Fluxquench–Theproblemwestudyhereisaquantum
quench where the flux θ is varied from θ(t < 0) = θ0 to J(t=∞)=Dθ0, (8)
θ(t ≥ 0) = 0. This sudden change of magnetic flux is
equivalent to imposing an instantaneous pulse of electric whilethefinite-timedynamicsisgovernedbytheregular
field E(t) = −∂ θ(t) = θ δ(t) to the fermions, and it part of the conductivity σ (ω). We will compare our
t 0 reg
induces some particle current at the initial time. In the numerical data on the flux quench with these LR predic-
present setup, the current is always uniform throughout tions later.
the system and thus we can define the current operator Before analyzing this problem in detail, we mention
by the average of local currents as28 somepossibleexperimentalrealizations. Thefluxquench
can be viewed as a sudden momentum shift for the par-
Jˆ= 1 (cid:88)(cid:16)c†c −c† c (cid:17)= 1 (cid:88)sin(q)c˜†c˜ , (2) ticles. It is therefore equivalent to a situation where a
2iN i i+1 i+1 i N q q
moving lattice stops abruptly at t=0 (the lattice veloc-
i q
ity provides the initial momentum shift). This situation
where c˜q := √1N (cid:80)Nr=−01cre−iqr is the annihilation oper- was experimentally realized10 with bosons trapped in an
ator in momentum space. The current is not a constant optical lattice. We may therefore expect that a similar
3
setup could be realized with fermions (1). Besides, a ues of interaction ∆ (note however that the oscillation
quantum quench using an artificial gauge field in optical period for ∆ = ±0.1 is too long to be measured accu-
lattices was also proposed11. This is a direct realization rately). For smaller initial fluxes, θ = π/6 to π/30, we
0
of the flux quench studied here, although bosons were observe a qualitatively different dynamics, depending on
considered. the sign of ∆. For attractive interactions (∆ > 0), the
oscillations (if any) are too slow to be visible within the
simulation time, and the relative decay of the current is
III. NUMERICAL RESULTS small. In that regime J(t) quickly reaches a stationary
value (with the possibility of some short time-scale and
small amplitude oscillations, as visible in the inset of the
Inthissectionwepresentnumericalresultsforthedy-
rightpanelofFig.4). Forrepulsiveinteractions(∆<0),
namics of the current. We employ the iTEBD method4,
some oscillations are visible, although for small initial
which enables to study the system in thermodynamic
flux θ and small |∆| their period can exceed the simula-
limit N →∞. The iTEBD is a numerical scheme based 0
tion time. Besides, the decay of the current is larger and
on the Matrix-Product State (MPS) representation of
the associated relaxation time scale is longer than in the
quantum many-body states in 1d. The MPS can nat-
attractive case.
urally describe a translationally invariant state of an in-
Inordertoquantifythevariousscalesassociatedtothe
finitely long 1d system. In the present problem, the ini-
current dynamics, we use two simple fitting functions,
tial state is translation invariant, and this symmetry is
preserved also at t > 0 by the post-quench Hamiltonian (cid:40)
f(t)=c+(A+Bcos(ωt+φ))e−t/τ
H0. Thus, strictly speaking, there is no finite-size effect , (9)
in our calculation. On the other hand, an exact descrip- g(t)=c+Ae−t/τ
tion of a given quantum state by an MPS generally re-
quires a matrix of infinite dimensions, but we need to where c,A,B,ω,φ,τ are fitting parameters. We use f(t)
approximate it by a finite-dimensional matrix in a prac- when some oscillations are visible within the simulation
tical calculation. The dimension of the matrix is called time scale (t<100), and g(t) otherwise. Some examples
bond dimension, and the use of a finite bond dimension offitsareshowninFig. 2. Amongthefittingparameters,
is a possible source of the error in the calculation. For a we focus on c,ω and τ, which correspond respectively to
groundstateofagapped1dsystem,anMPSwitha(suf- the long-time limit of the current, the frequency of the
ficientlylarge)finitebonddimensionisknowntoprovide oscillations, and the relaxation time.
analmostexactdescriptionofthewavefunction12. How-
ever, forthegroundstateofagaplesssystemwithanin-
finite correlation length, the finite bond dimension of an A. Oscillations
MPS approximation is known to introduce an effective
finite correlation length13. Nevertheless, the MPS de- As shown in Fig. 3, the oscillation frequency ω ex-
scription(andthusiTEBDalgorithm)withafinitebond tractedfromthefitsisapproximatelylinearin|∆|. Note
dimension can provide an accurate result on quantum however, the associated slope appears to depend on the
dynamics up to a certain time14, in particular concern- value of the initial flux. The relation ω ∝ |∆| seems to
ing local observables. Therefore we apply the iTEBD holdin gaplessphase (|∆|≤1). Butit mayalsobevalid
algorithm to the flux quench problem, to obtain the evo- beyond that regime, since ω also appears to be approxi-
lutionofthecurrentforacertainperiodoftimeafterthe mately linear |∆| in the regime −1.5(cid:46)∆≤1, although
quench. this relation clearly breaks down for ∆(cid:46)−1.5.
In practice, first, the ground state of the model (1) We comment on the relation between our numerical
withfluxθ isobtainedbysimulatinganimaginarytime- results and the LR theory. As described in the previous
0
evolution with the iTEBD algorithm. We then compute section, the LR theory (θ (cid:28) 1) relates the real-time
0
the real time-evolution using the Hamiltonian without dynamics of the current to the Fourier transform of the
flux, still with iTEBD. We carefully check the numer- regular part of the conductivity σ (ω). Even though
reg
ical errors by varying the time steps during the time- the applicability of the LR theory is not obvious when
evolution, as well as the bond dimension χ. The results θ is of the order of unity (as for θ = π/6,π/3 and
0 0
shown in the present paper were obtained by using bond π/2), the frequency ω of the observed oscillations may
dimensionsbetweenχ=500andχ=1200. Wecalculate originate from a peak in σ (ω). To our knowledge such
reg
the dynamics for various initial flux strengths θ ranging a structure in σ (ω) has not been explicitly discussed
0 reg
from π/30 to π/2, and for interaction ∆ from −2.0 to in the literature for the XXZ model, but similar results
0.8. Moredetailsonthenumericalcalculationsaregiven have been reported in studies on the finite-temperature
in Appendix B. Drude weight of this model15. Also, it is worthwhile to
The evolution of the expectation value of the current, point out that the current-current correlation function
J(t), is shown in Fig. 1. We summarize the observed dy- (cid:104)Jˆ(t)Jˆ(cid:105) shows a similar oscillatory behavior16 (basically
namics as follows. For large initial flux θ (θ (cid:38) π/3), the integral of this correlation function gives the time-
0 0
thecurrentshowssomedecayandoscillationsforallval- dependence of J(t)).
4
0.35 0.16 0.3
θ = π/2 0.5
0.3 0 0.14 0.8
-0.3
0.25 ∆= -0.1 0.12 -0.5
0.1
-0.3 0.1
0.2
-0.5 0.08
J(t) 0.15 0.5 0.3 J(t) 0.06 -0.8 θ = π/6
0.1 0
-0.8 0.04
0.8
0.05
0.02
-1.0 -1.0
0 -1.5 0 -1.2
-2.0 -1.2 -2.0 -1.5
-0.05 -0.02
0 5 10 15 20 25 30 0 20 40 60 80 100
t t
0.1 0.3 0.1
0.5 θ = π/10
-0.3 0
0.8 0.098
0.08 0.3
-0.5
0.096
0.5
0.06
-0.8 0.094
(t) 0.04 (t) 0.092 -0.3
J θ = π/10 J
0
0.02 0.09
-1.0 0.8
-1.5 0.088
0 -1.2
-2.0 0.086 -0.5
-0.02
0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100
t t
FIG. 1: Dynamics of the current after the quench for θ =π/2 (upper left), θ =π/6 (upper right), and θ =π/10 (bottom
0 0 0
left and right).
θ = π/6, line: fitting, dot: Data 2.5
0 θ = π/2
0.16 θ0 = π/3
0.14 ∆ = 0.5 2 θ00 = π/6
0.8
0.12 -0.5
0.1 1.5
0.08
) ω
(t -0.8
J 0.06 1
0.04
0.02 -1.0 0.5
0
-1.2
-2.0
-0.02 0
0 20 40 60 80 100 -2 -1.5 -1 -0.5 0 0.5 1
t
∆
FIG. 2: Fitting of numerical data for initial flux θ = π/6.
0
Roughly up to first oscillation, the empirical fitting works FIG. 3: Frequency ω from numerical fitting. Error bar is
well. TheblackarrowindicatestheLRpredictionofthelong- estimated by error of fitting.
time limit of the current for ∆=±0.5 (Eqn. (8)).
a population imbalance between q > 0 and q < 0 re-
sults in non-zero current. Figure. 4 shows the dynam-
1. Dynamics of the momentum distribution
ics of the momentum distribution for two cases17: θ =
0
π/3, ∆=−0.5 and θ =π/6, ∆=0.5. At t=0 the mo-
0
To investigate the nature of the current oscillations, mentum distribution is that of the ground state in pres-
wecalculatethemomentumdistributionoftheparticles, ence of a magnetic flux. It corresponds to the ground
n = (cid:104)c˜†c˜ (cid:105). Since we have (cid:104)J(t)(cid:105) = 1 (cid:80) sin(q)n (t), state in zero flux, but shifted by momentum θ (see Ap-
q q q N q q 0
5
θ = π/3, ∆= -0.5 θ = π/6, ∆= 0.5
0 0
1 1
t=0 t=0
0.9 0.9
t=3 t=7.6
0.8 t=5 0.8 t=15
0.7 t=10 0.7
t=15
0.6 0.6
q 0.5 q 0.5
n n
0.4 0.4 0.156
0.28 0.155 J(t)
0.3 00..2246 J(t) 0.3 0.154
0.22 0.153
0.2 0.2 0.2 0.152
0.1 000...111468 0.1 0.151 0 10 20 30 40 50
0 5 10 15 20 25 30
0 0 t
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0t.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
Momentum q /π Momentum q /π
FIG. 4: Time evolution of the momentum distribution n for (left) θ = π/3, ∆ = −0.5 and (right) θ = π/6, ∆ = 0.5. An
q 0 0
oscillatingdip/peakstructureisvisibleintheleftpanel(markedbyarrows),whereassuchstructureisabsentintherightcase.
See also the Supplemental Material17.
pendix A). Thus, at t = 0, n is a quasi (broadened) Finally, we point out that this anomalous dip (peak)
q
Fermi distribution with the shifted Fermi wave vectors structuremightbeobservedinrealexperiments,sincethe
±k(cid:48) = ±π/2+θ . After the flux is quenched to zero, momentum distribution is often accessible in cold atom
F 0
the distribution starts evolving. experiments.
An important observation is the simultaneous appear- To summarize this subsection, we found numerically
anceofa“dip”anda“peak”inn (leftofFigs.4and5). that the oscillation frequency is proportional to the
q
BothstructuresappeartooscillateinphasewithJ(t),as strengthofinteraction|∆|andthattheseoscillationsare
shown in the inset of Fig. 4. We also note that, in situa- governed by excitations located in momentum space far
tionswherecurrentoscillationsareabsent(rightpanelof from the shifted Fermi point.
Fig. 4), no dip/peak is observed. The dip and the peak
correspond to two momenta pdip and ppeak that are sep-
arated by π: ppeak = pdip −π. The momentum pdip is B. The long-time limit
plotted as a function of θ for several values of ∆ in the
0
right of Fig. 5. Here we discuss the long-time limit of the current.
For small θ , pdip approaches the shifted Fermi point From the LR theory in θ0, the long-time limit of the
0
k(cid:48) (→π/2)andthecurrent-carryingmodesbecomelow- current J(t = ∞) is given by J(t = ∞) = Dθ0, where
F
D is the Drude weight of the system. The latter is ex-
energymodes. Thisisexpectedsincetheinitialstateand
actlyknownfortheXXZmodel(equivalenttoourpresent
the ground state are energetically close to each other in
model) at zero temperature19:
thiscase. Howeverforfiniteθ weobservethatthemodes
0
responsible for the current oscillations are located at a
π sinµ
significant distance from the Fermi points, and are not D = , µ=arccos(−∆). (10)
4 µ(π−µ)
low-energy excitations.
The detailed dependence of pdip on the parameters θ0 D is non-zero only in the gapless phase (−1 ≤ ∆ < 1),
and ∆ is not yet understood but we can consider a sim- and vanishes in the gapped phase. Hence the LR theory
plified picture where only two characteristic modes gov- predictsthatJ(t=∞)isnon-zeroonlyingaplessphase.
ern the current dynamics. Located at pdip and ppeak = In addition, D is symmetric under ∆ ↔ −∆ (except for
pdip−π,thesemodesarerelatedthroughsomeUmklapp ∆=±1). So the long-time limit of the current does not
processes induced by the interactions. depend on the sign of ∆ in the LR theory.
As mentioned above, since pdip generically departs We can also estimate the long-time limit of the cur-
from k(cid:48) , the oscillations might not be described by the rent, J(t=∞), by fitting and extrapolating the numeri-
F
Tomonaga-Luttinger liquid (TLL) framework, where the cal data obtained for the finite time after the quench by
physics is entirely described in terms of low-energy exci- the iTEBD method. In the following, we shall compare
tations in the vicinity of the Fermi points5. As a com- J(t = ∞) estimated from the iTEBD calculation and
parison, a different global quench for the XXZ model that predicted by the LR theory. Note however that in
(equivalent to the present Hamiltonian) was considered the case of ∆ < 0 and small initial flux, it is difficult to
inRefs.14and18. There,thestrengthoftheinteraction evaluate the long-time limit from numerical data.
is suddenly changed and several aspects of the dynamics AsshowninFig.1, thenumericalresultsindicatethat
appeared to be well described by the TLL model. J(t = ∞) is non-zero in the gapless phase and zero in
6
∆= -0.5
1 θ =π/2, t=6.0 0.5 ∆= -0.3
0
00..89 θθ00==ππ//34,, tt==91.56.3 0.45 ∆∆== --00..58
θ =π/6, t=31.4
0
0.7 0.4
π
0.6
/
nq 0.5 dip 0.35
0.4 p
0.3
0.3
0.2 0.25
0.1
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
Momentum q /π θ0 /π
FIG. 5: Left: Momentum distribution for ∆=−0.5 with different θ . The time is chosen to match the first minimum of the
0
current after the quench. Right: Momentum of dip, pdip, versus θ .
0
1-J(t=∞)/Dθ
J(t=∞)/D 0
1
0.6
∆=0.3
∆=0.5
0.5 ∆=0.8
0.4 ∆∆== --00..35 0.1
∆= -0.8 ∆=0.3
∆= -1.0 ∆=0.5
0.3 LR theory
∆=0.8
0.2 0.01 ∆= -0.3
∆= -0.5
0.1 ∆= -0.8
∆= -1.0
0 0.001
0 0.1 0.2 0.3 0.4 0.5 0.6 0.1 1
θ θ
0 0
FIG. 6: Numerical data for (left) J(t=∞)/D and (right) normalized deviation from LR theory 1−J(t=∞)/(Dθ ).
0
the gapped phase, which is consistent with the LR re- may attribute this to superfluid correlations in the sys-
sult. On the other hand, it is clear from Fig. 6 that tem. Strictly speaking, superfluidity is absent in one di-
J(t=∞) rapidly deviates from the LR prediction when mension, but a superfluid-like response can be observed
θ increases. This reflects the nonlinearity of current in some dynamical properties of the system20, and these
0
as a function of the initial flux strength. For attrac- are expected to be stronger for attractive interactions
tive interactions (∆>0), the normalized deviation from (∆ > 0) than for repulsive ones (∆ < 0). In the case of
the LR theory (right of Fig. 6) shows power-law decay repulsiveinteractions,thelargernormalcomponentdissi-
at small θ , and these appear to be compatible with patesandthiswouldresultinsmallervaluesofJ(t=∞),
0
J(t = ∞) = Dθ + O(θ3). In contrast, for large re- as observed in Fig. 6.
0 0
pulsive interactions (∆<0), J(t=∞) strongly deviates
from the LR theory even for θ as small as π/30 ≈ 0.1.
0
Wenotethat, despitethegeneraldifficultyindoingsuch 1. Comparing the momentum distributions in the gapless
fitsandextrapolations,thedeviationfromtheLRtheory phase and in the gapped phase
cannot be attributed to some error in the extrapolation.
This is clear by comparing the raw finite-time data and
The dynamics of the momentum distribution shows a
the LR theory prediction, as shown in Fig. 2; any sensi-
qualitativedifferencebetweenthegaplessandthegapped
bleextrapolationwouldgivedifferentJ(t=∞)fromthe
phases, irrespective of initial flux θ . In the gapless
0
LRtheory. Wealsonotethatthisstrongnon-linearityin
phase, the shifted Fermi sea structure of the initial state
presenceofrepulsiveinteractionshasalreadybeennoted
appearstobequalitativelyrobustandsurvivesuptothe
in Ref. 11.
stationary regime (Fig. 4). The imbalance between the
Itisaninterestingfactthatthemagnitudeoftheabove numberofleft(q >0)andright(q <0)movingparticles
nonlinearities strongly depends on the sign of ∆. We in the stationary regime is the source of the persistent
7
currentJ(t=∞)(cid:54)=0. Ontheotherhand,inthegapped andissupportedbyAdvancedLeadingGraduateCourse
phase, the shifted Fermi sea structure of the initial state for Photon Science (ALPS) of JSPS. YON and MO
disappears over a relatively short time scale (Fig. 7). In thank the Yukawa Institute for Theoretical Physics at
that case, the whole momentum distribution moves to- Kyoto University (YITP). Discussions during the YITP
wards the center (q = 0) and the symmetry between workshopYITP-W-14-02on“HiggsModesinCondensed
q >0 and q <0 is restored, leading to J(t=∞)=0. Mater and Quantum gases” were useful to complete this
Thisdifferencebetweenthegaplessandgappedphases work. MO and YON are supported by MEXT/JSPS
might be related to the presence of additional con- KAKENHI Grant Nos. 25103706 and 25400392. GM
served quantities that exist in the gapless phase21–24. acknowledges V. Pasquier for useful discussions, and is
Those additional constants of motion, called quasi-local, supportedbyaJCJCgrantoftheAgenceNationalepour
are responsible for the ballistic transport and the non- la Recherche (Project No. ANR-12-JS04-0010-01). The
zero Drude weight at finite temperature in the gapless computation in this paper has been partially carried out
phase21. Inasimilarmanner, weexpecttheseadditional by using the facilities of the Supercomputer Center, the
conservedquantitiestopreventtherestorationoftheleft- InstituteforSolidStatePhysics,theUniversityofTokyo.
right symmetry in the momentum distribution.
Appendix A: Initial state and twist operator
C. Relaxation time
As mentioned in the introduction, in a ring of length
The relaxation time τ, extracted from the fits, is plot-
N, the zero-flux Hamiltonian H and the Hamiltonian
0
ted in Fig. 8. In general, larger |∆| and larger θ0 result with p∈Z flux quanta are related by a unitary transfor-
in smaller τ (faster decay). This is natural because the
mation:
timederivativeofthecurrentdJ(t)/dtisproportionalto
∆(Eqn.(3)). When∆=0(freefermionpoint), thecur- H =U H U−1, (A1)
rent is conserved and τ must be ∞. In agreement with θ0=2πp/N p 0 p
this fact, τ appears to diverge when |∆|→0.
where U (so-called twist operator) is defined as25
p
(cid:32) N−1 (cid:33)
IV. CONCLUSION U =exp 2ipπ (cid:88) xc†c . (A2)
p N x x
x=0
We have studied numerically a flux quench in an in-
teracting spinless fermion model in one dimension. This So, the ground state |ψ (cid:105) in presence of p flux quanta
p
quenchgeneratessomeparticlecurrentattheinitialtime can be expressed in terms of the zero-flux ground-state
and we monitored and analyzed quantitatively the cur- |ψ (cid:105):
0
rent dynamics that follows. The numerical data reveal
some current oscillations as well as some decay to a sta- |ψ (cid:105)=U |ψ (cid:105). (A3)
p p 0
tionary value. For repulsive interactions and a large ini-
tial flux, the frequency of those oscillations is propor- In other words, when p is an integer, the flux quench
tional to the strength of the interaction in the system. amounts to study the dynamics generated by H when
0
Remarkably,thedynamicsofthemomentumdistribution starting from the initial state U |ψ (cid:105).
p 0
revealsthattheseoscillationsaregovernedbyexcitations U has a simple action on the fermion creation opera-
p
located deep inside the (shifted) Fermi sea. In addition tors:
tothosenoveloscillations,thelong-timelimitofthecur-
(cid:18) (cid:19)
rent exhibits nonlinearities which are particularly strong U c†U−1 =exp 2ipxπ c†, (A4)
in presence of repulsive interactions. p x p N x
As future work, it seems important to understand the
origin of the specific “dip” momentum pdip that gov- which implies that it performs a momentum shift (or
erns the current oscillations. Another interesting direc- boost):
tion of research would be to compute the long-time limit
of the current (beyond the weak-flux regime where the U c˜†U−1 =c˜† . (A5)
p q p q+2πp/N
linear response theory applies) using integrability tech-
niques7,23,24. In the case of a noninteracting fermion problem (∆ = 0
in Eqn. (1)), U maps the Fermi sea |ψ (cid:105) to a “shifted
p 0
Fermi”seaU |ψ (cid:105), withFermipointslocatedat−π/2+
p 0
Acknowledgments 2πp/N and+π/2+2πp/N. Thelatterisanexactexcited
eigenstateoftheHamiltonianH . Therefore,inthenon-
0
YONacknowledgesY.Tada,I.Danshita,S.Furukawa, interacting case, the flux quench does not generate any
andS.C.Furuyaforvaluablecommentsanddiscussions, dynamics.
8
θ = π/2, ∆= -1.2 θ = π/6, ∆= -1.2
0 0
1 1
t=0 t=0
0.9 t=1 0.9 t=1
0.8 t=2 0.8 t=2
t=3
0.7 0.7 t=3
t=10
0.6 0.6 t=5
q 0.5 q 0.5 t=20
n n
0.4 0.4 00..1146
0.3 0.12 J(t)
0.3 0.2 J(t) 0.3 0 .00.81
0.2 0.1 0.2 00..0046
0.1 0 0.1 0.0 02
-0.1 0 5 10 15 20 25
0 0 1 2 3 4 5 6 7 8 9 10 0
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.t4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 t0.2 0.4 0.6 0.8 1
Momentum q /π Momentum q /π
FIG. 7: Time evolution of the momentum distribution in the gapped phase. Left panel: θ = π/2, ∆ = −1.2. Right panel:
0
θ =π/6, ∆=−1.2. Inbothcases,thewholemomentumdistributionmovestowardsthecenter(q=0)andthecurrentdecays
0
to zero.
104 104
θ = π/2
θ0 = π/3
103 θθ00 == ππ//610 103
θ0 = π/20
θ0 = π/30
102 0 102
τ τ
101 101
θ = π/2
θ0 = π/3
100 100 θθ 00= = π /π1/60
θ0 = π/20
θ0 = π/30
10-1 10-1 0
-2 -1.5 -1 -0.5 0 0.5 1 0.1 1
∆ -∆
FIG.8: Relaxationtimeτ versustheinteractionparameter∆. Leftpanel: semi-logplot. Rightpanel: log-logplot. Thedata
points where the decay is too slow to be reliably quantified with our calculations are omitted.
Appendix B: Numerical calculations at δτ = 0.001 converges, we compare the obtained en-
ergy with the exact one26. Our iTEBD energy matches
This appendix provides some details on the numerical the exact value with 5 or 6 digits. Note that the entan-
method. We first prepare the ground state of the Hamil- glement entropy of a half-infinite system should diverge
tonian with flux θ inthegaplessphase,andthatitisapproximatedhereby
0
a finite value (since χ is finite). See Refs. 13 and 14 for
H =−1(cid:88)(cid:16)e−iθ0c†c +h.c.(cid:17)−∆(cid:88)n˜ n˜ , related discussions.
θ0 2 i i+1 i i+1 Next, using iTEBD again, we calculate the real time-
i i
evolution using the Hamiltonian without flux:
(B1)
using an imaginary time-evolution (with iTEBD). A 1(cid:88)(cid:16) (cid:17) (cid:88)
H =− c†c +h.c. −∆ n˜ n˜ . (B2)
second-order Suzuki-Trotter decomposition is used and 0 2 i i+1 i i+1
we take the bond dimension χ between 500 and 1200. i i
The imaginary time step δτ is reduced gradually from We again use a second-order Suzuki-Trotter decomposi-
δτ = 0.1 to δτ = 0.001. δτ is reduced each time the tion and take a real-time step dt = 0.01 or 0.02. One of
imaginary time propagation with a coarser δτ has con- the largest obstacles to calculate the real time-evolution
verged. At each δτ, the convergence is checked by look- of a quantum system by iTEBD is the growth of the en-
ing at the energy and the entanglement entropy of a half tanglemententropy. Thisgrowthisusuallylinearintime
chain between two successive time steps. Our conver- for global quenches, and, as shown in Fig. 9, it appears
gencecriterionis10−8fortheenergyand10−6fortheen- to be the case for the present flux quench. In practice
tanglement entropy. After the imaginary time-evolution this has limited the accessible time scale to t(cid:46)100.
9
6 5.5
∆=0.5
θ = π/2 5 ∆=0.8
y 5 0 y ∆= -0.3
op op 4.5 ∆= -0.5
ntr 4 ntr 4 ∆∆== --01..80
nt e nt e 3.5 ∆= -1.5
e 3 e 3
m m
gle 2 ∆=0.5, χ=1000 gle 2.5
n ∆=0.8, χ=1000 n 2
a a
Ent 1 ∆∆== --01..80,, χχ==11000000 Ent 1 .15 θ0 = π/6, χ=700
∆= -1.5, χ=700
0 0.5
0 5 10 15 20 25 0 20 40 60 80 100
t t
FIG. 9: Growth of entanglement entropy of a half chain for θ = π/2 (Left) and θ = π/6 (Right) and various values of ∆.
0 0
One can check that the (weak) oscillations of the entanglement entropy are in phase with the oscillations of the current.
0.35
In order to check the accuracy of the simulated time-
θ = π/2, ∆ = -0.8 χ=200
0.3 0 χ=400 evolution, we have monitored the accumulated trunca-
χ=700 tionerrors, theconservationofenergy, andtheconserva-
0.25 χ=1000 tionofthenumberofparticles. Inadditiontomonitoring
thesevalues,weconfirmedtheaccuracyofourresultsfor
0.2
(t) a few ∆ and θ0 by comparing the results for several val-
J 0.15 uesofthetimestepdtandthebonddimensionχ. Asan
example, Fig. 10 shows how the expectation value of the
0.1
current is modified when increasing the bond dimension.
0.05
0
0 5 10 15 20 25 30
t
FIG. 10: Expectation value of the current, computed with
severalvaluesofthebonddimensionχ. Inthiscaseweretain
the data for χ = 1000, and only up to t = 20 (indicated by
arrow).
∗ Electronic address: [email protected] (2007).
1 A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat- 11 S. Peotta, C.-C. Chien, and M. Di Ventra, Phys. Rev. A
tore, Rev. Mod. Phys. 83, 863 (2011). 90, 053615 (2014).
2 M. Oshikawa, Phys. Rev. Lett. 90, 236401 (2003). 12 M. B. Hastings, J. Stat. Mech. 2007, P08024 (2007).
3 M.Mierzejewski,P.Prelovˇsek, andT.Prosen,Phys.Rev. 13 L. Tagliacozzo, T. R. de Oliveira, S. Iblisdir, and J. I.
Lett. 113, 020602 (2014). Latorre, Phys. Rev. B 78, 024410 (2008).
4 G. Vidal, Phys. Rev. Lett. 98, 070201 (2007). 14 F. Pollmann, M. Haque, and B. Do´ra, Phys. Rev. B 87,
5 T.Giamarchi,Quantumphysicsinonedimension (Oxford 041109 (2013).
University Press, 2004). 15 J.Herbrych,P.Prelovˇsek, andX.Zotos,Phys.Rev.B84,
6 M.Rigol,V.Dunjko,V.Yurovsky, andM.Olshanii,Phys. 155125 (2011).
Rev. Lett. 98, 050405 (2007). 16 C. Karrasch, J. H. Bardarson, and J. E. Moore, Phys.
7 A. De Luca, Phys. Rev. B 90, 081403 (2014). Rev. Lett. 108, 227206 (2012).
8 I. Danshita and A. Polkovnikov, Phys. Rev. B 82, 094304 17 See ancillary files of this arXiv post for some animations
(2010). showing the dynamics of the momentum distribution .
9 I. Danshita and A. Polkovnikov, Phys. Rev. A 84, 063637 18 C. Karrasch, J. Rentrop, D. Schuricht, and V. Meden,
(2011). Phys. Rev. Lett. 109, 126406 (2012).
10 J.Mun,P.Medley,G.K.Campbell,L.G.Marcassa,D.E. 19 B. Shastry and B. Sutherland, Phys. Rev. Lett. 65, 243
Pritchard, and W. Ketterle, Phys. Rev. Lett. 99, 150604 (1990).
10
20 T.Eggel,M.Cazalilla, andM.Oshikawa,Phys.Rev.Lett. 27 In these two cases the sudden change of flux creates an
107, 275302 (2011). excited eigenstate (a boosted Fermi sea for free fermions)
21 T. Prosen, Phys. Rev. Lett. 106, 217206 (2011). and no dynamics is generated (see also Appendix A). In
22 T. Prosen and E. Ilievski, Phys. Rev. Lett. 111, 057203 presence of a lattice and interactions, Umklapp processes
(2013). can change the particle current while conserving the total
23 R. G. Pereira, V. Pasquier, J. Sirker, and I. Affleck, J. lattice momentum.
Stat. Mech. 2014, P09037 (2014). 28 In presence of a non-zero vector potential θ(t), the ex-
24 T. Prosen, Nuclear Physics B 886, 1177 (2014). pression of the particle current is modified to Jˆ =
25 E. Lieb, T. Schultz, and D. Mattis, Ann. Phys. 16, 407 1 (cid:80) (cid:16)e−iθ(t)c†c −eiθ(t)c† c (cid:17)inordertobegauge-
(1961). 2iN i i i+1 i+1 i
26 M.Takahashi,ThermodynamicsofOne-DimensionalSolv- invariant.
able Models (Cambridge University Press, 1999).
