Simulating correlation-spreading dynamics
in two-dimensional quantum many-bodysystems by the tensor-network method
Ryui KANEKO

Dept. of Physics, Kindai Univ.
Collaborator: Ippei DANSHITA (Kindai Univ.)
R. Kaneko and I. Danshita, Commun. Phys. 5, 65 (2022)

1 / 32

Outline• Introduction• Quench dynamics by analog quantum simulation• Importance of comparison with numerical simulations• Numerical difficulty in simulating the dynamics of 2D quantum systems• Bose-Hubbard model: quench from a Mott insulating state• Motivation• Lack of reliable 2D methods• How far one can go by tensor-network states in 2D?
• Tensor-network method• Simple update, projected entangled pair states (PEPS)• Results• Good agreement with experimental results• Estimate group and phase velocities for smaller U/J that has not been
investigated in the experiment• Transverse-field Ising model: quench from a disordered state• Motivation• To what extent is PEPS useful?
• Preliminary results
2 / 32

Outline• Introduction• Quench dynamics by analog quantum simulation• Importance of comparison with numerical simulations• Numerical difficulty in simulating the dynamics of 2D quantum systems• Bose-Hubbard model: quench from a Mott insulating state• Motivation• Lack of reliable 2D methods• How far one can go by tensor-network states in 2D?
• Tensor-network method• Simple update, projected entangled pair states (PEPS)• Results• Good agreement with experimental results• Estimate group and phase velocities for smaller U/J that has not been
investigated in the experiment• Transverse-field Ising model: quench from a disordered state• Motivation• To what extent is PEPS useful?
• Preliminary results
3 / 32

Outline• Introduction• Quench dynamics by analog quantum simulation• Importance of comparison with numerical simulations• Numerical difficulty in simulating the dynamics of 2D quantum systems• Bose-Hubbard model: quench from a Mott insulating state• Motivation• Lack of reliable 2D methods• How far one can go by tensor-network states in 2D?
• Tensor-network method• Simple update, projected entangled pair states (PEPS)• Results• Good agreement with experimental results• Estimate group and phase velocities for smaller U/J that has not been
investigated in the experiment• Transverse-field Ising model: quench from a disordered state• Motivation• To what extent is PEPS useful?
• Preliminary results
4 / 32

Outline• Introduction
×10−2
2.0
Γ/J = 9.0
• Quench dynamics by analog quantum
simulation
C zz (r = 1,t)
0.0
• Importance of comparison with numerical simulations
0.5• Numerical difficulty in simulating zzthe dynamics
of 2D quantum systems
C (r = 2, t)
0.0• Bose-Hubbard model: quench from a Mott insulating
state
0.2
C zz (r = 3, t)

0.0• Motivation−0.2
0.2• Lack of reliable 2D methodszz
C (r = 4, t) 0.0• How far one can go by tensor-network
states in 2D?
0.1• Tensor-network method
C zz (r = 5, t)• Simple update, projected entangled
pair0.0states (PEPS)−0.1• Results√
0.5
C zz (r = 2, t)results• Good agreement with experimental
0.0−0.5• Estimate group and phase velocities forsmaller
U/J
that 4has not 6been
0
2time tJc investigated in the experimentquenchΓΓ=∞

PEPS, D = 5, 6, 7, 8
ED, Ns = 5 × 5• Transverse-field Ising model: quench from a disordered state• Motivation• To what extent is PEPS useful?
• Preliminary results
5 / 32

IntroductionAnalog quantum simulators

Let the nature do the quantum simulations using highly controllable experimental devices

Ultracold atoms in optical lattices[I.Bloch,Nature.453.1016(’08);
C.Gross,I.Bloch,Science.357.995(’17); W.Hofstetter,T.Qin,J.Phys.B:At.Mol.Opt.Phys.51.082001(’18)]

Rydberg atoms in
optical tweezer arrays[H.Bernien et al.,Nature.551.579(’17); A.Keesling etal.,Nature.568.207(’19)]
Δ |r〉
1,013 nmVij

420 nm

ΩR
δ

|e〉

ΩB

|r〉

Ω
|g〉

Trapped ion quantum computers[R.Blatt,C.F.Roos,Nat.Phys.8.277(’12); E.A.Martinez etal.,Nature.534.516(’16); M.Gärtner etal.,Nat.Phys.13.781(’17);
https://physicsworld.com/wpcontent/uploads/2018/12/IonQ-chip.png]|g〉Superconducting quantum circuits[R.Ma et al.,Nature.566.51(’19); Y.Ye et
al.,PRL.123.050502(’19)]

7 / 32

What do we want to do using analog quantum simulators?
• Solve problems that are hard to tackle by classical computers• Prepare the Hamiltonian corresponding to the problem andobtain the equilibrium state (e.g. the ground state)• Simulate Schrödinger equation→ Simulations of isolated quantum many-body systems haveattracted much interest∗ In experiments, quench is realized by very fast sweep• In general, simulating time evolution requires all the information of
eigenstates on classical computers→ It is much harder than the ground-state calculation
8 / 32

In the case of ultracold atoms on optical lattices...

9 / 32

In the case of ultracold atoms on optical lattices...

10 / 32

What do we want to clarify by simulating time evolution?
• How do isolated quantum many-bodysystems thermalize?
• What is the upper limit of the information
propagation (= Lieb-Robinson bound)?
cf. In relativistic system:
Upper limit = speed of light
Theoretical investigation is active recently
cf. Light-cone-like behavior in Bose-Hubbard
models[Y.Takasu et al., sciadv.aba9255(’20)][T.Kuwahara, K.Saito, PRL.127.070403(’21)]

Desirable to simulate the dynamics of correlation spreading
to answer these questions→ Longer-time experimental and numerical simulations are important
11 / 32

Comparisons between experimental and numerical simulations are desired

Propagation velocities can be obtained from equal-time correlations• Two characteristic velocities• e.g. 1D Bose-Hubbard simulator
Correlations after a quench• Phase velocity• Group velocity (≤ Lieb-Robinson bound)wave packet

envelope of
the wave packet

phase velocity

group velocity[M.Cheneau et al.,Nature.481.484(’11)]• In 1D, tensor-network
simulations with matrix productstates (MPS) are popular
12 / 32

Comparisons between experimental and numerical simulations are desired

Propagation velocities can be obtained from equal-time correlations• Two characteristic velocities• e.g. 1D Bose-Hubbard simulator
Correlations after a quench• Phase velocity• Group velocity (≤ Lieb-Robinson bound)[M.Cheneau et al.,Nature.481.484(’11)]

Quasiparticle pictureJ
wave packet

envelope of
the wave packet

phase velocity

group velocity• In 1D, tensor-network
simulations with matrix productstates (MPS) are popular

U

Dynamics of doublons and holons

13 / 32

Numerical simulations in 2D are extremely hard

Δ=(0,2)Δ=(1,1)

Δ=(0,1)[K.Nagao et al., PRR.3.043091(’21)]• e.g. Quench dynamics in the 2D Bose-Hubbard model• Semiclassical approach (truncated Wigner approximation) is not
powerful enough to reproduce the intensity of correlations• Extend the 1D MPS wave functions to 2D
Examine the accuracy of the 2D tensor-network method
14 / 32

Motivation• Numerical simulations of time evolution on classical computers• Crosscheck and predict experimental results• Numerical simulations in 2D are extremely hard so far• Focus on• 2D Bose-Hubbard model• 2D transverse-field Ising model
to examine the accuracy of the 2D tensor-network method

15 / 32

Tensor-network method

Tensor-network states: MPS and PEPS
• Wave function for quantum spin systems:X|ψi =Cs1 ,s2 ,...,sN |s1 , s2 , . . . , sN i#elements = O(eN ){si }• In 1D: Matrix product state (MPS)

Cs1,s2,...,sN

s1 s2

sN

≃s s
1

2

Ai[si]
sisN• In 2D: Projected entangled pair state (PEPS), tensor product state• Dphys = 2S + 1 for spin S(chosen to be sufficiently large forsoft-core bosons)• D = 1: direct product state• D ≥ 2: entangled state• Wave functions are moredim. of local :Dphys
accurate for larger D
Hilbert space• Translational invariant PEPS
(spin, boson, ...)
can treat infinite systems

virtual :D
bond[T.Nishino et al.,PTP.105.409(’01); F.Verstraete, J.Cirac, arXiv:cond-mat/0407066]

17 / 32

Simulating real-time evolution by infinite PEPS
• Real-time evolution of infinite PEPS: |ψ(t)i = e−itH |ψ(0)it=
e−iδtHjj

Time-evolving block decimation in 2D(= simple update) [comp. cost: O(D 5 )][H.C.Jiang,Z.Y.Weng,T.Xiang(’08); P.Corboz et al.(’10)]• Calculation of expectation values for infinite PEPS:
cornerenvironment
Corner transfer matrix renormalizationgroup method [comp. cost O(D 10 )]1edge[R.J.Baxter(’68); T.Nishino,K.Okunishi(’96,’97);-3
R.Orus,G.Vidal(’09)]
10

1
0.95
0.90
0.85
0.80

∆Ε/Ε0

<X>

<X>

∆Ε/Ε0

0.98
• Previous
studies on 2D quench dynamics:
00.96 transverse-field Ising model (tr.-field: hx = ∞ → hx )
e.g.
c
-3
-10 bond dimension D
0.94 . ~/J accessible by increasing
Time
0
0.2
0.4
0
0.2
0.4
10

-3

0
-3

∆Ε/Ε0

<X>

-10
0
1
0
1
0.5
0.5
1
-3
time
time
10
D=2
0.80
D=3
[A.Kshetrimayum et al.,Nat.Commun.8.1291(’17);
P.Czarnik et al.,PRB.99.035115(’19);
D=4
0.60
C.Hubig,J.I.Cirac,SciPost.Phys.6,031(’19)]
D=5
0
D=6
0.40
D=7
0.20
D=8
-3

18 / 32

Quench dynamics in the Bose-Hubbard model
Motivation:• Reproduce experimental results• Examine the parameter region that has not been exploredNumerical setup: Wish to calculate |ψ(t)i = e−iĤt |ψ0 i• Square Bose-Hubbard model: Ĥ =

X

Ĥij

hiji
†

†

Ĥij = −J(âi âj + âj âi ) +U2z[n̂i (n̂i − 1) + n̂j (n̂j − 1)] −µz(n̂i + n̂j ) (z = 4)

[V.Murg et al.,PRA(’07); J.Jordan et al.,PRB(’09); A.Kshetrimayum et al.,PRL(’19); S.S.Jahromi and R.Orus,PRB(’19); P.Schmoll et al.,PRL(’20);
W.-L.Tu et al.,JPCM(’20); H.-K.Wu et al.,PRA(’20); P.C.G.Vlaar and P.Corboz et al.,PRB(’21)]• Simple update by e.g. e−idtĤ/~ ∼Y−idtĤij /~

e

hiji

(use second-order Suzuki-Trotter decomposition in practice)• Very fast (τQ > 0) and sudden (τQ = 0) quenches from Mott insulator ⊗i |ni = 1i
• Experimental setup: U/J ∼ 100 → 19.6 in τQ = 0.1ms(U/J = 19.6 > 16.74 = Uc /J: Mott insulating region)

μ/U

Mott
Insulator

Superfluid

n=1

J/U
• Use tensor-network library TeNeS[Y.Motoyama et al., Comp.Phys.Commun.279.108437(’22); https://github.com/issp-center-dev/TeNeS, https://github.com/TsuyoshiOkubo/pTNS]

20 / 32

Numerical results: Comparison with the experiment at U/J = 19.6
Crsp (t) =X1††hâ (t)âj (t) + âj (t)âi (t)i2Ns r −r =r iij• Consider finite quench time as in the experiment• Nearly conserved energy for 0 ≤ tJ/~ . 0.4• Physical quantities are likely to be converged forthis short time
21 / 32

Numerical results: Comparison with the experiment at U/J = 19.6
Crsp (t) =X1††hâ (t)âj (t) + âj (t)âi (t)i2Ns r −r =r iijΔ=(0,2)Δ=(1,1)

Δ=(0,1)•••••

Consider finite quench time as in the experiment
Nearly conserved energy for 0 ≤ tJ/~ . 0.4
Single-particle correlations agree very well
How about other parameter regions?
How does the propagation velocity behave?
22 / 32

Numerical results: Estimate propagation velocities from ha0 ar i and hn0 nr i†

Consider a sudden quench
Energy is conserved for longer time
tJ/~ . 0.9 when U/J ∼ 5
Can capture peaks up to |r| ≤ 3
X1††hâ (t)âj (t) + âj (t)âi (t)i2Ns r −r =r ii
Xj
1
Crdd (t) =(hn̂i (t)n̂j (t)i − 1)
Ns r −r =r
Crsp (t) =ij• vphase : captured by single-particle correlation ha†0 ar i
• vgroup : captured by density-density correlation hn0 nr i
23 / 32

Numerical results: Estimate propagation velocities from ha0 ar i and hn0 nr i
†

X
1
†
†
hâ (t)âj (t) + âj (t)âi (t)i2Ns r −r =r ii
Xj
1
Crdd (t) =(hn̂i (t)n̂j (t)i − 1)
Ns r −r =r
Crsp (t) =ij• vphase : captured by single-particle correlation ha†0 ar i
• vgroup : captured by density-density correlation hn0 nr i
24 / 32

Numerical results: U dependence of velocity• For U . zJ (z = 4), single-particle picture (mean-field-like picture) holds
vgroup ∼ 4J/~ [K.Nagao et al.,PRA.99.023622(’19)]• For U  J, quasi-particle picture holds
vgroup ∼ 6J/~ × [1 + O(J 2 /U 2 )] [M.Cheneau et al.,Nature.481.484(’11)]• vgroup estimated from hn0 nr i consistent with• single-particle group velocity deep in superfluid region• strong-coupling result near criticality• vphase and vgroup gradually converge to the same value as U/J is decreased
25 / 32

Conclusions: Bose-Hubbard case
• Quench dynamics from Mott insulator in 2D Bose-Hubbard model• Simulation by infinite PEPS using simple update• Compare PEPS simulations with experiments→ Good agreement for tJ/~ . 0.4 at U/J = 19.6
μ/U

J/U

Mott
Insulator

Superfluid

quench
time

n=1

J/U

• Estimate velocity of correlation spreading for smaller U/Jvphase and vgroup gradually converge
to the same value as U/J is decreased
Might be helpful for future experiments

R. Kaneko and I. Danshita, Commun. Phys. 5, 65 (2022)
26 / 32

Quench dynamics in the 2D transverse-field Ising model
Motivation:• How good is the 2D tensor-network method in this case?
• How does the group velocity for spin correlations look?
(Compare with the recently updated Lieb-Robinson bound)Analog quantum simulations of the quantum Ising model
by Rydberg-atom arrays[H.Bernien et al.,Nat.551.579(’17); A.Keesling et al.,Nat.568.207(’19);
E.Guardado-Sanchez et al.,PRX.8.021069(’18); V.Lienhard et al.,PRX.8.021070(’18);
D.Bluvstein et al.,Science.371.1355(’21); P.Scholl et al.,Nat.595.233(’21); S.Ebadi et al.,Nat.595.227(’21); …]

1, D : dim.
2
i
i
hiji

X
X
X
X V D∆=ΩSix − (∆ − V D)
Siz + V
Siz Sjz +
−
4
2
i
i
i

H =Ω

X

Six − ∆

X

ni + V

X



ni nj

ni = Siz +

hiji

Δ |r〉
1,013 nmVij

420 nm

ΩR
δ

|e〉

ΩB

|r〉

Ω
|g〉

bases
Ω
∆
V

|g〉

Rydberg atoms

quantum Ising model|gi, |ri
Rabi frequency
detuning
van der Waals interaction| ↓i, | ↑ itransverse fieldlongitudinal field
Ising spin interaction

Δ/V

Ω/V

Very recently, real-time dynamics for # qubits > 200
28 / 32

Numerical setup: 2D transverse-field Ising model• Ground-state phase diagram
Disorder

h/J

2

1

0

[R.Kaneko et al.,JPSJ.90.073001(’21)]

H = +J

AF order00.5

XSiz Sjz − Γ

hiji
1

X

Six − h

X

i

Siz

i

1.5

Γ/J

• For simplicity, focus on h = 0 case
→ Map to ferromagnetic model by appropriate unitary transformation
XX
H = −J
Siz Sjz − Γ
Six
hiji

i

• Sudden quench from the Γ = ∞ ground state | →→ · · · →iquench|ψ(t)i = e−iHt/~ | →→ · · · →iΓcΓ=∞
29 / 32

Numerical results: Extract group velocity from spin correlations
2.0

C zz (r = 3, t)

C zz (r = 4, t)

C zz (r = 5, t)
√

0.0
0.5
0.0

Γ/J = 9.0
D=842

longitudinaltransverse

0.2
0.0
−0.2
0.2

0

0.0

0

1

2

5

6

1.0
0.5

longitudinal, r = 5
longitudinal, r = 4
longitudinal, r = 3

2, t)

ED, Ns = 5 × 5

3
4
peak time tJ

1.5 D = 8

0.1
0.0
−0.1

0.5
0.0
PEPS, D = 5, 6, 7, 8 −0.5

C zz (r =distance r

C zz (r = 2, t)

6v group(r)/JC zz (r = 1, t)

×10−2
Γ/J = 9.0

0.0
0

2

4
time tJ

6

0

2

4

6
Γ/J

transverse, r = 5transverse, r = 4
transverse, r = 3

8

10• Distances longer than exact diag are calculable• Estimated group velocity for Γ  J: vspin /J = 1.07 ± 0.20
• Best Lieb-Robinson bound for any correlations in 2D TFIsing model for Γ  J:vLR /J = 7.55
[Z.Wang, K.R.A.Hazzard, PRXQuantum.1.010303(’20)]• vspin  vLR
→ Our data is more meaningful when we need to compare the spin correlations

30 / 32

Conclusions: Ising case
• Quench dynamics from the disordered state in 2D transverse-field Isingmodel• Simulation by infinite PEPS using simple updatev group(r)/J

1.5 D = 8

quench

1.0
0.5
0.0

Γc

Γ=∞

longitudinal, r = 5
longitudinal, r = 4
longitudinal, r = 3

0

2

4

6
Γ/J

transverse, r = 5transverse, r = 4
transverse, r = 3

8

10• Our estimate of the group velocity: vspin /J ∼ 1• This is much smaller than the current best Lieb-Robinson bound:vLR /J = 7.55
• Our group velocity and spin correlations are helpful for crosschecking
experimental data
31 / 32

Conclusions• Simulating the dynamics of 2D systems by the tensor-network method with iPEPS• Focus on the quench in the 2D Bose-Hubbard and transverse-field Ising models
Bose-Hubbard case
• Good agreement with the experiment

Ising case
• Group velocity satisfies the
Lieb-Robinson bound(but the value is much smaller thanthe bound)v group(r)/J

1.5 D = 8

• Examine the parameter region that has
not been explored

1.0
0.5
0.0

longitudinal, r = 5
longitudinal, r = 4
longitudinal, r = 3

0

2

4

6
Γ/J

transverse, r = 5transverse, r = 4
transverse, r = 3

810• vspin  vLR is intrinsic?
• Recent vLR is still loose?
• Provide numerical data that can becompared with future experiments
32 / 32