How to simulate Trotterized time evolution for the Fermi-Hubbard model

In this guide, we show how to use ffsim to simulate Trotterized time evolution for the Fermi-Hubbard model using the split-operator method. In the split-operator method, the Hamiltonian is expressed as a sum of two terms, the one-body part and the two-body part, for the purposes of the Trotter product formula.

First, let’s create a Fermi-Hubbard Hamiltonian and choose the number of electrons, which will remain fixed throughout.

import ffsim

norb_x = 4
norb_y = 4

norb = norb_x * norb_y
nelec = (2, 2)
n_alpha, n_beta = nelec

op = ffsim.fermi_hubbard_2d(
    norb_x=norb_x,
    norb_y=norb_y,
    tunneling=1.0,
    interaction=4.0,
    chemical_potential=0.0,
    periodic_x=True,
    periodic_y=False,
)

Before we can simulate Trotterized time evolution, we need to convert the operator to an instance of the DiagonalCoulombHamiltonian class.

dc_hamiltonian = ffsim.DiagonalCoulombHamiltonian.from_fermion_operator(op)

Now, let’s apply time evolution to an initial state. Here, we’ll use the ground state of the one-body part of the Hamiltonian as the initial state.

import numpy as np

# Construct the initial state.
eigs, orbital_rotation = np.linalg.eigh(dc_hamiltonian.one_body_tensor)
initial_state = ffsim.slater_determinant(
    norb,
    occupied_orbitals=(range(n_alpha), range(n_beta)),
    orbital_rotation=orbital_rotation,
)

# Set the evolution time.
time = 1.0

# Perform the Trotterized time evolution
final_state = ffsim.simulate_trotter_diag_coulomb_split_op(
    initial_state,
    dc_hamiltonian,
    time,
    norb=norb,
    nelec=nelec,
    n_steps=5,
    order=0,
)

To check the error from Trotterization, we can compute the result of exact time evolution and compare it with the final state we got.

import scipy.sparse.linalg

# Convert the Hamiltonian to a LinearOperator
linop = ffsim.linear_operator(dc_hamiltonian, norb=norb, nelec=nelec)

# Compute the exact result of time evolution
exact_state = scipy.sparse.linalg.expm_multiply(
    -1j * time * linop,
    initial_state,
    traceA=-1j * time * ffsim.trace(dc_hamiltonian, norb=norb, nelec=nelec),
)

# Compute fidelity between results from Trotterized evolution and exact evolution
fidelity = abs(np.vdot(final_state, exact_state))

print(f"Fidelity of Trotter-evolved state with exact state: {fidelity}")

Fidelity of Trotter-evolved state with exact state: 0.9725640629741702

To reduce the Trotter error, you can increase the number of steps or use a higher-order formula.

final_state = ffsim.simulate_trotter_diag_coulomb_split_op(
    initial_state,
    dc_hamiltonian,
    time,
    norb=norb,
    nelec=nelec,
    n_steps=10,
    order=1,
)

fidelity = abs(np.vdot(final_state, exact_state))

print(f"Fidelity of Trotter-evolved state with exact state: {fidelity}")

Fidelity of Trotter-evolved state with exact state: 0.9999380169878249