free_fermions_exact.py

on github (download).

r"""Provides exact entropies for time evolution of free fermion models."""
# Copyright (C) TeNPy Developers, GNU GPLv3

import numpy as np
import matplotlib.pyplot as plt


def hopping_matrix(L, t=1., mu_staggered=0., boundary_conditions='open'):
    r"""Generate matrix for fermionic hopping in real-space c_i operators.

    Provide hopping matrix h_{ij} for

    .. math ::
        H = \sum_{ij} h_{ij} c_i^\dagger c_j
          = t \sum_{i} (c_i^\dagger c^\j + h.c.)  + \mu_s (-1)^{i} n_i

    for open boundary conditions on L sites.
    """
    assert boundary_conditions in ['open', 'periodic']
    H = np.zeros((L,L))
    for i in range(L - int(boundary_conditions == 'open')):
        H[i, (i+1) % L] = t
        H[(i+1) % L, i] = t
    for i in range(L):
        H[i, i] =  mu_staggered*(-1)**i
    return H


def charge_density_wave(L):
    """Initial state |010101010...>"""
    psi0 = np.mod(np.arange(0,L),2)
    return np.diag(psi0)


def time_evolved_state(H_hop, psi0, time_list):
    """Evolve an initial correlation matrix with H_hop."""
    E, U = np.linalg.eigh(H_hop)
    for time in time_list:
        X = np.dot(np.dot(np.conj(U),np.diag(np.exp(1j*time*E))),U.T)
        psi_t = np.dot(X,np.dot(psi0, np.conj(X.T)))
        yield psi_t  # c^dagger_i c_j


def entanglement_entropy(psi, bond=None):
    """Calculate entanglement entropy for cutting at given bond (default: half-chain)."""
    L = psi.shape[0]
    if bond is None:
        bond = L // 2
    z = np.linalg.eigvalsh(psi[:bond, :bond])
    z = z[np.logical_and(z > 1.e-15, z < 1. - 1.e-15)]
    S = - np.sum(z * np.log(z) + (1. - z) * np.log(1. - z))
    return S


def XX_model_ground_state_energy(L, h_staggered, boundary_conditions='open'):
    """
    """
    H_hop = hopping_matrix(L, 2., 2.*h_staggered, boundary_conditions=boundary_conditions)
    E = np.linalg.eigvalsh(H_hop)
    E_shift = 0. if L % 2 == 0 else - h_staggered
    # the shift stems from the 1/2 in mapping sigmaz = (n - 1/2) for the h_s terms
    # which cancels out for even L due to the alternating sign of h_s
    return np.sum(E[:L//2]) + E_shift

def XX_model_time_evolved_entropies(L, h_staggered, time_list, bond=None,
                                    boundary_conditions='open'):
    r"""Half-chain entanglement entropies for time evolving Neel state with XX chain.

    The XX chain given by the hamiltonian (here, X,Y,Z = Pauli matrices)

    .. math ::
        H = \sum_{i} (X_i X_{i+1} + Y_i Y_{i+1}) - h_s (-1)^i Z_i

    maps to free fermions through the Jordan-Wigner transformation.
    This function returns the entropies for a time evolution starting from the Neel state
    |up down up down ...>.
    """
    psi_0 = charge_density_wave(L)
    H_hop = hopping_matrix(L, 2., h_staggered, boundary_conditions=boundary_conditions)
    S_list = []
    for psi_t in time_evolved_state(H_hop, psi_0, time_list):
        S_list.append(entanglement_entropy(psi_t))
    return np.array(S_list)


if __name__ == "__main__":
    L = 100
    h_staggered = 1.
    time_list = np.linspace(0, 20., 100)
    for h_s in [0., 1., 2.]:
        S_list = XX_model_time_evolved_entropies(L, h_s, time_list)
        plt.plot(time_list, S_list, label="$h_s = {h_s:1.1f}$".format(h_s=h_s))
    h_s = 0.
    S_list = XX_model_time_evolved_entropies(L, h_s, time_list,
                                           boundary_conditions='periodic')
    plt.plot(time_list, S_list, linestyle='--',
             label="$h_s = {h_s:1.1f}$, periodic".format(h_s=h_s))
    S_list = XX_model_time_evolved_entropies(10, h_s, time_list)
    plt.plot(time_list, S_list, linestyle=':',
             label="$h_s = {h_s:1.1f}, L=10$".format(h_s=h_s))

    plt.xlabel('$t$')
    plt.ylabel('$S$')
    plt.legend(loc='best')
    plt.show()