Ground States and Area Laws

Area Laws

When a 1D Hamiltonian $H$ has a gap $\Delta = E_0 - E_1 > 0$ between the ground state energy $E_0$ and the first excited state $E_1$ then the entanglement entropy $S_A$ of the ground state satisfies:

$$S_A \leq C , \ C \in \mathbb{R}$$

for a bipartition of the system $A|B$. This is proven in reference and is an example of more general area laws of entanglement entropy:

$$S_A \propto V^{d-1}$$

where $V$ is the volume of the system and $d$ its physical dimension. The idea is that for these special low-lying states, the only contributions to their entanglement entropy come from the area of contact between the subregion $A$ and its complement $B$, not from the size of $A$ itself. In this tutorial we will explicitly investigate this bound in 1D systems.

Gapped Systems

A rich playground for studying gapped systems is spin chains, chains of localized spins that interact through some magnetic interactions. In this tutorial we will examine the 1D Heisenberg model:

$$H = \sum_{i=1}^{N-1} = J ( S^{+}_i S^{-}_{i+1} + S^{-}_i S^{+}_{i+1} + Z_i Z_{i+1} ) + h ( X_i + Z_i)$$

with open boundary conditions. In the code below, we will consider $J=3.2$, $h = 0.2$ . For such a choice, the system has a gap which is a necessary ingredient to test the 1D scaling law of entanglement entropy. The energy gap $\Delta$ also plays a crucial role in quantum phase transitions, as at the critical point of the transition the gap closes, i.e. $\Delta = 0$. As such, extracting the gap information can help in identifying when a change of phase has occurred.

Make the Hamiltonian

We can construct the Hamiltonian with an operator_sum:

def make_hamiltonian(real j, real h, integer N)
{
    var H = operator_sum(as_real)

    for(var s = 0 ; s < N-1 ; ++s)
    {
        H += j * Splus(s)*Sminus(s+1) + j * Sminus(s)*Splus(s+1) + j * Z(s)*Z(s+1)
    }
    for(var s = 0 ; s < N ; ++s)
    {
        H += h * X(s) + h * Z(s)
    }
    
    return H
}

This function takes in our choice $J$ for the heisenberg interaction strength, the $[101]$ magnetic field strength $h$, and the number of sites in the chain and returns $H$.

Obtaining the Ground State with Exact Diagonalization

Once we have built our Hamiltonian, we can pass it our iterative eigensolvers to obtain a number of desired eigenvector and eigenvalue pairs

var options = ["is_hermitian" : true, "is_real" : true, "krylov_dimension": 6, "tolerance": 1.0e-6]
var results = iterative_eigensolver(H, 2**N, neig, options)

The first line builds the options map to pass to our solver, it helps in determining which solver to use to get the best performance. The Hamiltonian is Hermitian and comprised of only real valued operators and thus the solver will select the Lanczos algorithm to use. Specifying that the operator is real can provide a speed up in the computation as the storage of complex numbers is double that of real numbers. Moreover, as this algorithm uses a Krylov subspace method, we can also specify the number of Krylov vectors to use as well as the numerical precision we want to use. The error of this iterative method is controlled by the residual norm $||Hv - \lambda v||_2$, the input precision establishes when this norm is considered small enough.

The second line calls the solver and begins the computation, taking the operator_sum, the total Hilbert space size, the number of eigenvalues to target and the options map. The results are stored in a list containing the energies and the eigenstates of $H$ that we will use for our study.

Getting the Entanglement Entropy

Next we need to calculate the entanglement entropy from the ground state we have just obtained. First we must perform the partial trace to part of our system to obtain the reduced density matrix $\rho_A$ for the sites we are interested in

$$\rho_A = \text{Tr}_B[\rho] = \sum_{m_b \in \mathcal{H}_B} \langle m_b | \rho | m_b \rangle.$$

This operation sums over the states of the Hilbert space corresponding the physical system $B$ we will be tracing out. The entanglement entropy of the subsystem $A$ is given by

$$S(\rho_A) = -\text{Tr}[\rho_A \log(\rho_A)] = -\sum_n \lambda_n \log(\lambda_n)$$

where $\lambda_n$ are the eigenvalues of $\rho_A$. We can perform both of these calculations in aleph with two functions

var evec = results[1]
var rdm = reduced_density_matrix(evec.col(0),[0..5])
var EE = von_neumann_entropy(rdm)

The first line obtains $\rho_A$ by passing our ground state and the sites we are interested in, the non-included sites are traced over. In this case we are only studying the interval of sites $[0,4]$. The second line calculates the entanglement entropy from the input density matrix.

Putting It All Together

We can now study the scaling of entanglement entropy of the ground state with system size by using the above code and making a for loop.

var j = 3.2
var h = 0.2
var Nmax = 20

for(N: [10..Nmax+1][::2])
{
    var H = make_hamiltonian(j,h,N)
    var options = ["is_hermitian" : true, "is_real" : true, "krylov_dimension": 6, "tolerance": 1.0e-6]
    var results = iterative_eigensolver(H, 2**N, 2, options)
    var eval = results[0]
    var evec = results[1]

    var rdm = reduced_density_matrix(evec.col(0),[0..5])
    var ee = von_neumann_entropy(rdm)

    var out = "N = ${N} , Gap = ${eval[1]-eval[0]}, S(rho) = ${ee}"
    print(out)
}

What's Next?

• Can we conclude that the ground state satisfies an area law?
• Can we instead keep N fixed and grow our subregion size?
• What is the effect of changing the number of krylov vectors used?
• How does our choice of where we cut our subsystem affect the results?
• How big a system can be used for these calculations?
• How stable is this gap with system size? Should it remain in the thermodynamic limit?

Documentation Contributors

Sebastien J Avakian