Simulation of a time-depend 1-D Bose-Hubbard model

ITensor Julia Questions
1

Dear ITensor community,

I’m trying to simulate a time-dependent 1-D Bose-Hubbard model, where the Hamiltonian: H = H_{\text{hop}} + H_{\text{int}} + H_{\text{tilt}}, where

H_{\text{hop}} = \sum_{i=1}^{L-1} J(t) (\hat{a}_i^{\dagger} \hat{a}_{i+1} + \hat{a}_{i+1}^{\dagger} \hat{a}_i)

is the hopping Hamiltonian, L is the number of lattice sites, and J(t) \in (0,1) is the hopping strength;

H_{\text{int}} = \sum_{i=1}^L U(t) \hat{n}_i(\hat{n}_i-1)

is the interaction Hamiltonian, and U(t) \in (-1,1);

H_{\text{tilt}} = \sum_{i=1}^L \Delta(t)(i-S) \hat{a}_i^{\dagger} \hat{a}_i

is the “tilt” Hamiltonian where S = (L+1)/2 is the center of the 1-D lattice, and \Delta(t) \in (-1,1). For this model, J(t), U(t), \Delta(t) can be controlled, and the QFI
F(\psi)=4\left( \langle \hat{S}_z^2\rangle_{\psi}-\langle\hat{S}_z\rangle_{\psi}^2 \right) of the state need to be computed.

using ITensors, ITensorMPS

struct SimulationConfig
    n_sites::Int
    n_particles::Int
    time_step::Float64
    final_time::Float64
    cutoff::Float64
end

const CONFIG = SimulationConfig(
    3,              # n_sites
    100,            # n_particles
    0.1,            # time_step
    10.0,           # final_time
    1E-8,           # cutoff
)

Here I’m showing a model with 3 lattice sites and 100 bosons, first I define on-site Hamiltonian (H_{\text{int}}+H_{\text{tilt}}) and hopping

function build_onsite_hamiltonian(site_index::Int, U::Float64, Δ::Float64, s, config::SimulationConfig)
    center = (config.n_sites + 1) / 2
    h = U * op("N * N", s[site_index]) + U * (-1.0) * op("N", s[site_index])
    h += Δ * (site_index - center) * op("N", s[site_index])
    return h
end


function build_hopping_hamiltonian(site_index::Int, J::Float64, s)
    h = -J * op("Adag", s[site_index]) * op("A", s[site_index + 1])
    h += -J * op("A", s[site_index]) * op("Adag", s[site_index + 1])
    return h
end

Then I use a second order trotterization

function make_trotter_gates_2nd(J::Float64, U::Float64, Δ::Float64, dt::Float64, s, config::SimulationConfig)
    gates = ITensor[]

    # First half of one-site gates
    for j in 1:config.n_sites
        hj = build_onsite_hamiltonian(j, U, Δ, s, config)
        Gj = exp(-im * dt/2 * hj)
        push!(gates, Gj)
    end

    # Two-site hopping gates
    for j in 1:(config.n_sites - 1)
        hj = build_hopping_hamiltonian(j, J, s)
        Gj = exp(-im * dt * hj)
        push!(gates, Gj)
    end

    # Second half of one-site gates
    for j in 1:config.n_sites
        hj = build_onsite_hamiltonian(j, U, Δ, s, config)
        Gj = exp(-im * dt/2 * hj)
        push!(gates, Gj)
    end

    return gates
end

To compute the QFI, for the 3-site model 4(\hat{n}_1^2+\hat{n}_3^2-2\hat{n}_1\hat{n}_3 - (\hat{n}_3-\hat{n}_1)^2) 

function compute_qfi(psi, s, cutoff::Float64)
    n_vals = expect(psi, "N")
    n1, n3 = n_vals[1], n_vals[3]

    N_vals = expect(psi, "N * N")
    N1, N3 = N_vals[1], N_vals[3]

    Op1 = op("N", s[1])
    Op3 = op("N", s[3])
    psi_temp = apply(Op3, psi; cutoff)
    n1n3 = real(inner(apply(Op1, psi; cutoff), psi_temp))

    # QFI formula
    QFI = real(4 * ((N1 + N3 - 2*n1n3) - (-n1 + n3)^2))
    return QFI
end

Then I do the time evolution. In this example, at each time step, J,U,\Delta are random values

function run_simulation(config::SimulationConfig)
    s = siteinds("Boson", config.n_sites; dim=config.n_particles + 1, conserve_qns=true)

    psi = MPS(s, ["$(config.n_particles)", "0", "0"])

    # Storage for observables and parameters
    QFI_history = Float64[]
    J_history = Float64[]
    U_history = Float64[]
    Δ_history = Float64[]

    # Time evolution loop
    for t in 0:config.time_step:config.final_time
        # Compute and store QFI
        QFI = compute_qfi(psi, s, config.cutoff)
        push!(QFI_history, QFI)

        # Sample random parameters for this time step
        J_t = rand()
        U_t = 2 * rand() - 1
        Δ_t = 2 * rand() - 1

        push!(J_history, J_t)
        push!(U_history, U_t)
        push!(Δ_history, Δ_t)

        # Build and apply 2nd order Trotter gates
        gates = make_trotter_gates_2nd(J_t, U_t, Δ_t, config.time_step, s, config)
        psi = apply(gates, psi; cutoff=config.cutoff)
        normalize!(psi)
    end

    return QFI_history, J_history, U_history, Δ_history
end

QFI_data, J_data, U_data, Δ_data = run_simulation(CONFIG)

However, this implementation is not very efficient; for 3-site, 100 bosons, T=10 and dt=0.1, it took around 100 minutes to run on a cluster; I extend T=20, then it took more than 5 hours to finish. I would like to know how to improve the performance, thank you.

2

Have you tried using TDVP?

3

Thank you for the suggestion. I tried to implement TDVP, but it’s slower than trotterization; even for a small system (3-site, 20 bosons), the TDVP method took around 1000 sec (on my laptop) to finish the time evolution (T=10, dt=0.02). Below is the code:

using ITensors, ITensorMPS
using DelimitedFiles


struct SimulationConfig
    n_sites::Int
    n_particles::Int
    time_step::Float64
    final_time::Float64
    cutoff::Float64
    maxdim::Int
end

# Initialize simulation configuration
const CONFIG = SimulationConfig(
    3,              # n_sites
    20,             # n_particles
    0.02,           # time_step
    10.0,           # final_time
    1E-8,           # cutoff
    50,            # maxdim for TDVP
)


function build_hamiltonian_mpo(J::Float64, U::Float64, Δ::Float64, s, config::SimulationConfig)
    center = (config.n_sites + 1) / 2

    os = OpSum()

    for j in 1:config.n_sites
        # Interaction
        os += U, "N * N", j
        os += -U, "N", j
        # Tilt 
        os += Δ * (j - center), "N", j
    end

    # Hopping 
    for j in 1:(config.n_sites - 1)
        os += -J, "Adag", j, "A", j + 1
        os += -J, "A", j, "Adag", j + 1
    end

    H = MPO(os, s)
    return H
end


function compute_qfi(psi, s, cutoff::Float64)
    n_vals = expect(psi, "N")
    n1, n3 = n_vals[1], n_vals[3]

    N_vals = expect(psi, "N * N")
    N1, N3 = N_vals[1], N_vals[3]

    # Compute ⟨n₁*n₃⟩ directly using inner product
    Op1 = op("N", s[1])
    Op3 = op("N", s[3])
    psi_temp = apply(Op3, psi; cutoff)
    n1n3 = real(inner(apply(Op1, psi; cutoff), psi_temp))

    # QFI formula
    QFI = real(4 * ((N1 + N3 - 2*n1n3) - (-n1 + n3)^2))
    return QFI
end


function run_simulation(config::SimulationConfig)
    s = siteinds("Boson", config.n_sites; dim=config.n_particles + 1, conserve_qns=true)

    # Initialize state: all particles in first site
    psi = MPS(s, ["$(config.n_particles)", "0", "0"])

    # Storage for observables and parameters
    QFI_history = Float64[]
    J_history = Float64[]
    U_history = Float64[]
    Δ_history = Float64[]

    # Time evolution loop
    for t in 0:config.time_step:config.final_time
        # Compute and store QFI
        QFI = compute_qfi(psi, s, config.cutoff)
        push!(QFI_history, QFI)

        # Sample random parameters for this time step
        J_t = rand()
        U_t = 2 * rand() - 1
        Δ_t = 2 * rand() - 1

        push!(J_history, J_t)
        push!(U_history, U_t)
        push!(Δ_history, Δ_t)


        # Build Hamiltonian MPO for current parameters
        H = build_hamiltonian_mpo(J_t, U_t, Δ_t, s, config)

        # Time evolve using TDVP
        psi = tdvp(H, -im * config.time_step, psi;
                   cutoff=config.cutoff,
                   maxdim=config.maxdim,
                   normalize=true,
                   outputlevel=0)
    end

    return QFI_history, J_history, U_history, Δ_history
end


QFI_data, J_data, U_data, Δ_data = run_simulation(CONFIG)
4

Usually with TDVP you can take much larger steps than with trotterization. Especially if you use the expand method beforehand to enlarge the support. If you play around I’m sure you can improve the performance quite a bit, but whether or not it will beat trotterization…