Random Phase Product Sate for Canonical Ensemble
Revised 7/1/2020 00:19:00

Toshiaki IITAKA (飯高敏晃)
RIKEN Center for Computational Science (R-CCS),
2-1 Hirosawa, Wako, Saitama, 351-0198, Japan.
e-mail address: tiitaka@riken.jp
URL: http://www.iitaka.org/rpps/
Method of random phase product state (RPPS) is proposed to calculate canonical ensemble average of quantum
systems described with matrix product states and also with tensor network states in general. The RPPS method is anextension of the method of random phase state for the conventional full Hilbert space representation. The validity of themethod is confirmed by comparing the average energy of N-site spin-1/2 antiferromagnetic Heisenberg chain model with
open boundary conditions with the result of direct method (for up to N=14) and minimally entangled typical thermal state(METTS) method (for N=100). Numerical advantages of the RPPS method such as parallelization, combined calculationof thermal averages at different temperatures, parameters for controlling error are discussed. View point of self-averaging
for the super-convergence of random state method is emphasized in addition to that of typicality.

According to the textbook of statistical mechanics[1],physical quantities of the system in thermal equilibrium
at inverse temperature β = 1/ k BT are calculated basedon the canonical ensemble, in which each microscopic
quantum state φn with energy En is realized according
to the Gibbs probability Pn ∝ e− β En.Calculating

physical quantities of quantum many body system atfinite temperature is one of the most important problems
in computational physics, which requires, however,
formidable computational resources if direct
diagonalization method is applied to large systems.
There are two sources of the difficulty. First, thedimension of the state vector increases exponentially as
the system size increases. For example, state vector of thespin one-half Heisenberg chain has dimension of
M = 2 N where N is the number of spins [2]. Thisdifficulty may be mitigated by using tensor network
representation [3, 4] of quantum states, which efficiently
represents physically meaningful states near the groundstates using much less parameters than the conventional
full Hilbert space representation.
Second, the number of excited states to be included in
the calculation increases exponentially as temperature
increases. [1]. In contrast, at very low temperature, only
few states close to the ground states contribute to thethermal average and these low lying states are efficiently
calculated with Lanczos methods [5, 6]. To cope with thesecond difficulty, method of random state has beensuccessfully applied to many problems (see [6] andreferences therein) and also discussed in terms oftypicality[7-12]. Among them, method of random phasestate [13-18] has superior property of smaller statisticalfluctuation. The essence of random (phase) state method
is sampling the entire Hilbert space with much smaller
number of random states than the full basis set. So far,the method of random state has been used mostly with thefull Hilbert space representation and only few authors
used it with tensor network representation[7-10].

In the following, the method of random state isreviewed and extended to tensor network representation
using, as an example, matrix product state (MPS) [3, 4]of N-site spin-1/2 antiferromagnetic Heisenberg chain
model with open boundary condition. However,extension to tensor network in general is straight forward.
The Hamiltonian of N-site spin-1/2 antiferromagnetic
Heisenberg chain with open boundary condition isdefined as
N −1

 

∑ S ⋅S

H
=

i

i =1



(1)

i +1

where Si is spin-1/2 operator at site i.
A quantum state of this system is expressed in the full
Hilbert space representation as

φ = ∑ cσ σ σ 1 σ N = ∑ cσ σ1Nσ1σ Nσ(2)where the tensor of rank N, cσ

= cσ1σ N , is thecomponents of φ projected onto the complete
orthonormal basis set

σ = σ 1 σ N = σ 1 ⊗ σ 2 ⊗  σ N(3)constructed as the direct product of eigenstates of local
spin operatorS zi σ i = σ i σ i2with eigenvalues σ i = ±1 .
In order to sample the entire Hilbert space, random
state is defined [16] as

Φ =

1
ξσ σ
∑
M σ

(5).

(4)

σ

where ξ is a set of independent and identically
distributed (i.i.d.) complex random numbers that have thefollowing statistical relations

ξσ

=0

(6)

ξ σ '*ξ σ= δ σ 'σξ σ 'ξ σ= 0 (8)(7)ξ σ = exp iθ σ σ(14)
and thermal average of physical quantity A withcanonical ensemble is calculated as(9)numbers in the range [ 0, 2π ) where M=2N is thedimension of the full Hilbert space. Random statecomposed of random phase numbers is called random
phase state [14, 16].
Then, normalization and completeness of random state
are expressed as

=1

(10)

Φ(β ) A Φ(β )Tr[e − β H A]=
Tr[e − β H ]=A

are M independent uniform real random

Φ|Φ

Random state can be regarded as the quantum state at
infinite temperature or β = 0 because all state is
occupied with equal probability according to (10)-(13).
Thermal state at inverse temperature β is defined asΦ ( β ) = e − β H /2 Φ

where the double bracket and the asterisk indicate
statistical average and complex conjugate, respectively.
As a special case of such random numbers, random
phase number is defined aswhere θ(13). The double brackets indicate statistical average.
Normalization factor of 1/M is omitted for clarity.

Φ(β ) Φ(β )(15)The essence of the random state method is that the
double summation in expectation value (13) is reduced tosingle summation in trace by using statistical relation(7)as illustrated in Fig.1. When statistical averageΦ X Φ

is calculated as sample average with

N sample realizations of random states, the fluctuation of

δX

the average decreases asΦ Φ =1(11)of large numbers. The statistical fluctuation ofΦ X Φ for each realization of random state is given

(12)

by

in the special case of random phase state, and

Φ Φ

I
=
M

according to the law

N sampleorrespectively. Matrix element of operator X gives thetrace of X ,1
Φ X Φ = ∑ ξ σ'*ξ σ X σ 'σ
M σ'σ
.
1
1
= =
∑ X σσ M tr [ X ]
M σ
(13)

δX = Φ X Φ −

Φ X Φ

(

)

=

1
ξ σ1*ξ σ2 − δ σ1σ2 X σ1σ2∑
M σ1σ2(16)
Therefore the magnitude of the fluctuation is expected as2(ξ ξ − δ ) X
σ3 *

=

(ξ ξ

1
∑
M 2 σ1σ2σ3σ4

δX
=
σ4

σ 3σ 4

σ1

σ2 *

)

− δ σ1σ2 X σ*1σ2 ×

σ 3σ 42
1
1X+ 2 
∑σ1σ 2
2
M σ1 ≠σ2M ξσ− 1 ∑ X σσ
 σ

4

2

.
(17)
According to the inequality

ξ

σ 4

≥

ξ

σ 2

2

=
1

the second term of (17) becomes zero if and only if

ξ σ = 1 for all σ . Since this condition is satisfied by

Fig. 1. Graphic representation of random state method:(a) normalization (10); (b) completeness (12); (c) trace

random phase state, it will give the smallest fluctuation
among the random states in the given basis set. In manycases, where X is such as combinations of local
operators, X σ1σ 2 is sparse (i.e. most of the matrix

elements are zero) and the right hand side of (17)decreases rapidly as

1

1
2− N
= =
N rand M

δX ∝
2

(18)

Now let us extend the method of random state to tensor
network representation. Here, Matrix Product State(MPS), a special form of tensor network decompositionσof c , useful for one dimensional system, is used as an
example. MPS is defined as

∑c σ
σ

=

σ

∑  A   A    A σ   A σ  σ=
∑ µ ∑µ Aµ µσ Aµ µσ  Aµ µσ Aµ σµ σ
[1]σ1

=

[2]σ 2[ N −1] N −1

[N ] Nσσ

1 N −1= δ σ 'ση σ 'η σ= 0 (27)(26)σ

as system size N and the number of random numbers per
sample N rand increases. This is self-averaging[16, 19].

φ

η σ '*η σ

[1] 1

[2] 2

[ N −1] N −1

[N ] N

0 1

1 2

N −2 N −1

N −1 Nµ=0 µ=
N 1

although M = 2 N random numbers of η are notindependent of each other because they are composed ofonly 2 N χ

2[ i ]σ 2 N independent random numbers, ξ µν.

Therefore random MPS also satisfies all equations for thefull Hilbert space representation including randomvariables up to quadratic order.
Random Phase Product State (RPPS) is defined as
random phase MPS with χ = 1 and regarded as quantum
states at infinite temperature where all state realizes withequal probability and entanglement vanishes. Since it is a
product state, RPPS is defined independent of the
topology of tensor network. Topology of tensor network
matters when the operators such as A and eapplied to the state. RPPS is expressed as−τ 0 H /2

are

(19)
σ

where c is decomposed into the product of N matrices[ i ]σ A of dimension χ or bond dimension[3, 4]. The

bond index µi runs from 1 to χ with constraintµ=µ=
1 due to the open boundary conditions. The
0Nσ
M = 2 N coefficients of c in the full Hilbert space

representation is reduced to ~2 N χ coefficients in MPS.
2

If χ were exponentially large (i.e. 2 N χ ≈ 2 ), MPS
would have sufficient flexibility to express the exact state[4].
Random MPS is defined as
2

1
ησ σ
∑
M σ

Ψ =
where coefficient ηη σ = χ − ( N −1)/2

σ

N

(20).

is defined byξ µ µσ ξ µ µσ ξ µ µσ ξ µ σµ∑µ µ1 N −1µ=0 µ=
N 1

[1] 1

[2] 2

[ N −1] N −1

[N ] N

0 1

1 2

N −2 N −1

N −1 N(21)[ i ]σwith i.i.d. complex random numbers, ξ µν , that satisfy
the statistical relations,[ i ]σξ µνξ=0[ i ]σ '* [ j ]σ
µ 'ν 'µνξ[ j ]σξ µ[i'ν]σ' 'ξ µν(22)= δ ijδσ 'σ δ µ ' µ δν 'ν

=0

(23)

(24).

Random MPS composed of random phase numbers iscalled random phase MPS.
σUsing (21)-(24), it is shown that η satisfies statistical
relations corresponding to (6)-(8),

ησ

=0

(25)

Ψ =

(

(28).
Random phase number is used so that statistical fluctuationfrom diagonal matrix element is suppressed. However,other i.i.d. random numbers such as Gaussian random
number will work as well.
The fluctuation of the matrix element (13) with RPPS
decreases as

δX ∝
2

1

1
=
N rand 2 N

as system size N increases, which is much slower than(18) for the full Hilbert space representation. RPPS canbe also regarded as the random matrix product state(RMPS) in Ref. [7] with χ = 1 and is analogous toclassical product state (CPS) in minimally entangled
typical thermal state method (METTS) [20, 21].
Thermal MPS state is defined with RPPS as(29)
Ψ ( β ) = e − β H /2 Ψand the thermal average of physical quantity A isexpressed as

=
A

)

1
exp ( iθ [1]σ1 ) exp ( iθ [2]σ 2 ) exp iθ [ N ]σ N σ∑
M σ

Ψ(β ) A Ψ(β )Tr[e − β H A]=
Tr[e − β H ]Ψ(β ) Ψ(β )(30)where A and H are matrix product operators (MPO)[4].
In numerical calculation, the canonical ensemble
average (30) is calculated as shown in Fig. 3 (a).
Manipulation of MPS is taken care of by numericallibrary ITensor [22]. First, an initial RPPS is generated.
Second, the thermal stateLΨ (i ) ( β=) e −τ 0 H /2  Ψ (i ) is calculated by−τ H /2
repeatedly applying Boltzmann MPO [23], e 0 , with

small imaginary time step τ 0 = β / L , where L is a large

H a( H − c)=(31)
is used in the program for numerical convenience. The
constants a>0 and c are determined so that all
eigenvalues of H lie in the range [0, −1] . Third,
expectation value of the observable

(a)
0
-0.1

energy per site

integer. Entanglement may be increased by theapplication of operators and the resulting MPS may
require larger bond dimensions to satisfy the required
cutoff error condition. Such automatic control of bonddimension is taken care of by iTensor. The shifted andscaled Hamiltonian

-0.2
-0.3
-0.4
-0.5

0

5

10

(i )

(i )

Z (i ) =Ψ (i ) ( β ) Ψ (i ) ( β ) are calculated. AfterAβ==∑ Z (i )N sample

Z

0

.

i =1

The accuracy of RPPS method is controlled by the
number of initial RPPS, N sample , scaled imaginary time
step τ 0 , cutoff parameter of SVD, ε cutoff [19]. The
accuracy of RPPS method for the Heisenberg chain up to
N spin = 14 with N sample = 10000 , τ 0 = 0.02 ,ε cutoff = 10
error of 10

−13

−4

25

χmax

A(i )
(i )

20

5

N sample
i =1

25

(b)

canonical average of the observable at inverse
temperature β is obtained as
(i )

20

10

averaging with N sample initial random states, the

∑A

15

β

A =
Ψ ( β ) A Ψ ( β ) and partition function(i )has been confirmed within statisticalin comparison with the results of the direct

method[24].The accuracy of N spin = 100 calculation
with control parameters N sample = 100 , τ 0 = 0.02 ,ε cutoff = 10−10 was confirmed within statistical error of
10−4 in comparison with the result of METTS method
using the same control parameters. In Fig. 2, the averageenergy and the maximum bond dimension are shown as afunction of inverse temperature. It is noted that relatively
small bond dimension of χ  10 is sufficient for
converged results.

0

5

10

β

15

Fig. 2. Temperature dependence of (a) average energy,and (b) maximum bond dimension for 100 site spin-1/2
antiferromagnetic Heisenberg chain with open boundary
condition and control parameters N sample = 100 ,τ 0 = 0.02 , ε cutoff = 10−10 .
Many other numerical methods are known for thermal
calculation of tensor network state. Among them, Matrix
Product Purification method [19, 25, 26] utilizes
augmented sites (or ancilla) and poorly scale to large
systems at low temperature. RMPS studied by Garnerone[7] reduces to RPPS when χ = 1 . METTS method [20,21] scales to large systems even at low temperature. The
similarity of algorithm of RPPS method and METTS
method are demonstrated in Fig. 3. In this sense, RPPS
method is algorithm of new generation that has inheritedgene from both RMPS method and METTS method.
RPPS method has better parallelization efficiency than
METTS method because it has no sequential loop. In
addition, RPPS method can calculate thermal averages atintermediate inverse temperatures β l = lτ 0 , where=
l 1, 2,3, , L − 1 , as the byproducts of the calculation
at β = Lτ 0 . Unlike METTS method, RPPS method does
not have problems related to the ergodicity of sampling
because it calculates the trace over the entire Hilbert
space as shown in (13). As METTS method does, RPPS
method starts with an initial product state and the bonddimension is increased during imaginary time evolution

as cutoff accuracy requires, which makes RPPS method
numerically more efficient than RMPS method.
FORALL initial RPPS(a) − β H  (i )
(β )
exp 
Ψ ( i )=
 Ψ
 2 

Ai =
Ψ (i ) ( β ) A Ψ (i ) ( β )Zi =Ψ (i ) ( β ) Ψ (i ) ( β )END FORALL
N sample

=
A β

∑ A

−β H

Tr[e
A]
=
Tr[e − β H ]i =1
N sample

∑ Z
i =1

i

i

FORALL initial CPS(b) −β H Ψ=exp i (β )
 Ψi
 2 

Ai =

Ψ (i ) ( β ) A Ψ (i ) ( β )Ψ (i ) ( β ) Ψ (i ) ( β )Project Ψ i ( β ) onto CPS Ψ i+1 .

END FORALL

A β =

1

N sample

N sample

i =1

∑ A

i

Fig. 3. Flow chart of (a) Random Phase Product State(RPPS) method and (b) Minimally Entangled Typical
Thermal State (METTS) method [20, 21]. Sequential
loop is designated with bold line.

Since the program for RPPS method with Heisenberg
chain model is written in C++ language using hybridparallelization with OpenMP thread parallelization for
ITensor library and MPI process parallelization for the
initial RPPSs, the program will run extremely efficient on
PC clusters and supercomputers such as FUGAKU. The
source code will be submitted to ITensor website [22]and MateriApps website [27] after the publication of thispaper.
For large and almost homogeneous systems,
macroscopic quantities calculated with random statemethod in the full Hilbert space representation sometimes
converge extremely rapidly and calculation with only one
random state provides a good approximation. This superconvergence has been discussed in terms of typicality [712]. Derivation of random state methods in this article
also emphasizes the self-averaging [14, 19] due to therandom phase on the off-diagonal matrix elements in(13).
In summary, we have introduced a numerical methodthat is useful for calculating canonical ensemble averagein tensor network representation.
This research was supported by MEXT as“Exploratory Challenge on Post-K computer” (Challengeof Basic Science – Exploring Extremes through MultiPhysics and Multi-Scale Simulations). This research used
computational resources of the Hokusai system provided
by RIKEN and the HPCI system provided by JCAHPC
through the HPCI System Research Project hp190167
and hp190196.
References[1]
F. Reif, Fundamentals of Statistical and
Thermal Physics (Waveland Press, 2009).
[2]
J. J. Sakurai, and J. Napolitano, Modern
Quantum Mechanics (Pearson, San Francisco, 2013).
[3]
R. Orus, Annals of Physics 349, 117 (2014).
[4]
U. Schollwock, Annals of Physics 326, 96(2011).
[5]
J. K. Cullum, and R. A. Willoughby, Lanczos
Algorithms for Large Symmetric Eigenvalue
Computations: Vol. 1: Theory (Siam, 2002), Vol. 41.
[6]
J. Schnack, J. Richter, and R. Steinigeweg,
Physical Review Research 2, 013186 (2020).
[7]
S. Garnerone, T. R. de Oliveira, and P. Zanardi,
Phys. Rev. A 81, 032336 (2010).
[8]
S. Garnerone et al., Phys. Rev. A 82, 052312(2010).
[9]
S. Garnerone, Phys. Rev. B 88, 165140 (2013).
[10]
S. Garnerone, and T. R. de Oliveira, Phys. Rev.
B 87, 214426 (2013).
[11]
S. Sugiura, and A. Shimizu, Phys. Rev. Lett.
108, 240401 (2012).
[12]
S. Sugiura, and A. Shimizu, Phys. Rev. Lett.
111, 010401 (2013).
[13]
R. Alben et al., Phys. Rev. B 12, 4090 (1975).
[14]
T. Iitaka et al., Phys. Rev. E 56, 1222 (1997).
[15]
T. Iitaka, and T. Ebisuzaki, Phys. Rev. Lett. 90,
047203 (2003).
[16]
T. Iitaka, and T. Ebisuzaki, Phys. Rev. E 69,
057701 (2004).
[17]
M. Machida, T. Iitaka, and S. Miyashita, J Phys
Soc Jpn 74, 107 (2005).
[18]
K. Seki, and S. Yunoki, Phys. Rev. B 101,
235115 (2020).
[19]
M. Binder, and T. Barthel, Phys. Rev. B 92,
125119 (2015).
[20]
S. R. White, Phys. Rev. Lett. 102, 190601(2009).
[21]
E. M. Stoudenmire, and S. R. White, New
Journal of Physics 12, 055026 (2010).
[22]
https://itensor.org.
[23]
M. P. Zaletel et al., Phys. Rev. B 91, 165112(2015).
[24]
http://gtoth.eu/qubit4matlab.html.
[25]
F. Verstraete, J. J. García-Ripoll, and J. I. Cirac,
Phys. Rev. Lett. 93, 207204 (2004).

[26]M. Zwolak, and G. Vidal, Phys. Rev. Lett. 93,
207205 (2004).
[27]
https://ma.issp.u-tokyo.ac.jp/.