DQPT in the Ising model(@id demo_dqpt)

In this tutorial we will try to reproduce the results from this paper. The needed packages are

using MPSKit, MPSKitModels, TensorKit

Dynamical quantum phase transitions (DQPT in short) are signatures of equilibrium phase transitions in a dynamical quantity - the loschmidth echo. This quantity is given by $L(t) = \frac{-2}{N} ln(| < \psi(t) | \psi(0) > |)$ where $N$ is the system size. One typically starts from a groundstate and then quenches the hamiltonian to a different point. Non analycities in the loschmidth echo are called 'dynamical quantum phase transitions'.

In the mentioned paper they work with

$H(g) = - \sum^{N-1}_{i=1} \sigma^z_i \sigma^z_{i+1} + g \sum_{i=1}^N \sigma^x_i$

and show that divergences occur when quenching across the critical point (g₀ → g₁) for $t^*_n = t^*(n+\frac{1}{2})$ with $t^* = \pi/e(g_1,k^*)$, $cos(k^*) = (1+g_0 g_1) / (g_0 + g_1)$, $e(g,k) = \sqrt{(g-cos k)^2 + sin^2 k}$.

The outline of the tutorial is as follows. We will pick $g₀ = 0.5$, $g₁ = 2.0$, and perform the time evolution at different system sizes and compare with the thermodynamic limit. For those $g$ we expect non-analicities to occur at $t_n ≈ 2.35 (n + 1/2)$.

First we construct the hamiltonian in mpo form, and obtain the pre-quenched groundstate:

H₀ = transverse_field_ising(; g=-0.5)

L = 20
ψ₀ = FiniteMPS(rand, ComplexF64, L, ℂ^2, ℂ^10)
ψ₀, _ = find_groundstate(ψ₀, H₀, DMRG());

[ Info: DMRG init:	obj = +9.799964091770e+00	err = 1.5223e-01
[ Info: DMRG   1:	obj = -2.040021714743e+01	err = 2.4038839149e-02	time = 0.09 sec
[ Info: DMRG   2:	obj = -2.040021715170e+01	err = 6.0313575856e-07	time = 0.07 sec
[ Info: DMRG   3:	obj = -2.040021773534e+01	err = 1.6799456960e-05	time = 0.14 sec
[ Info: DMRG   4:	obj = -2.040021786694e+01	err = 1.9058246307e-06	time = 0.09 sec
[ Info: DMRG   5:	obj = -2.040021786703e+01	err = 1.1474711603e-06	time = 0.05 sec
[ Info: DMRG   6:	obj = -2.040021786703e+01	err = 4.3837579221e-10	time = 0.02 sec
[ Info: DMRG conv 7:	obj = -2.040021786703e+01	err = 1.9834477158e-11	time = 0.50 sec

Finite MPS quenching

We can define a helper function that measures the loschmith echo

echo(ψ₀::FiniteMPS, ψₜ::FiniteMPS) = -2 * log(abs(dot(ψ₀, ψₜ))) / length(ψ₀)
@assert isapprox(echo(ψ₀, ψ₀), 0, atol=1e-10)

We will initially use a two-site TDVP scheme to dynamically increase the bond dimension while time evolving, and later on switch to a faster one-site scheme. A single timestep can be done using

H₁ = transverse_field_ising(; g=-2.0)
ψₜ = deepcopy(ψ₀)
dt = 0.01
ψₜ, envs = timestep(ψₜ, H₁, 0, dt, TDVP2(; trscheme=truncdim(20)));

"envs" is a kind of cache object that keeps track of all environments in ψ. It is often advantageous to re-use the environment, so that mpskit doesn't need to recalculate everything.

Putting it all together, we get

function finite_sim(L; dt=0.05, finaltime=5.0)
    ψ₀ = FiniteMPS(rand, ComplexF64, L, ℂ^2, ℂ^10)
    ψ₀, _ = find_groundstate(ψ₀, H₀, DMRG())

    ψₜ = deepcopy(ψ₀)
    envs = environments(ψₜ, H₁)

    echos = [echo(ψₜ, ψ₀)]
    times = collect(0:dt:finaltime)

    for t in times[2:end]
        alg = t > 3 * dt ? TDVP() : TDVP2(; trscheme=truncdim(50))
        ψₜ, envs = timestep(ψₜ, H₁, 0, dt, alg, envs)
        push!(echos, echo(ψₜ, ψ₀))
    end

    return times, echos
end

finite_sim (generic function with 1 method)

Infinite MPS quenching

Similarly we could start with an initial infinite state and find the pre-quench groundstate:

ψ₀ = InfiniteMPS([ℂ^2], [ℂ^10])
ψ₀, _ = find_groundstate(ψ₀, H₀, VUMPS());

[ Info: VUMPS init:	obj = +5.000419382862e-01	err = 3.8507e-01
[ Info: VUMPS   1:	obj = -1.062780898337e+00	err = 2.4151374798e-02	time = 1.31 sec
┌ Warning: ignoring imaginary component -3.817802691569172e-6 from total weight 2.4343214394239743: operator might not be hermitian?
│   α = 1.7081089443203243 - 3.817802691569172e-6im
│   β₁ = 0.14004438155194618
│   β₂ = 1.7287776826281838
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component -4.33258309393697e-6 from total weight 3.290563788912538: operator might not be hermitian?
│   α = 2.784222669687663 - 4.33258309393697e-6im
│   β₁ = 0.13141322180957496
│   β₂ = 1.7488981501547187
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
[ Info: VUMPS   2:	obj = -1.063544409753e+00	err = 1.4337063091e-05	time = 0.08 sec
[ Info: VUMPS   3:	obj = -1.063544409973e+00	err = 1.7048015038e-07	time = 0.01 sec
[ Info: VUMPS   4:	obj = -1.063544409973e+00	err = 1.7183182449e-08	time = 0.01 sec
[ Info: VUMPS   5:	obj = -1.063544409973e+00	err = 5.8865202016e-10	time = 0.01 sec
[ Info: VUMPS conv 6:	obj = -1.063544409973e+00	err = 6.8020365686e-11	time = 1.42 sec

The dot product of two infinite matrix product states scales as $\alpha ^N$ where $α$ is the dominant eigenvalue of the transfer matrix. It is this $α$ that is returned when calling

dot(ψ₀, ψ₀)

1.0000000000000018 - 7.06078012856404e-16im

so the loschmidth echo takes on the pleasant form

echo(ψ₀::InfiniteMPS, ψₜ::InfiniteMPS) = -2 * log(abs(dot(ψ₀, ψₜ)))
@assert isapprox(echo(ψ₀, ψ₀), 0, atol=1e-10)

This time we cannot use a two-site scheme to grow the bond dimension, as this isn't implemented (yet). Instead, we have to make use of the changebonds machinery. Multiple algorithms are available, but we will only focus on OptimalEpand(). Growing the bond dimension by $5$ can be done by calling:

ψₜ = deepcopy(ψ₀)
ψₜ, envs = changebonds(ψₜ, H₁, OptimalExpand(; trscheme=truncdim(5)));

a single timestep is easy

dt = 0.01
ψₜ, envs = timestep(ψₜ, H₁, 0, dt, TDVP(), envs);

With performance in mind we should once again try to re-use these "envs" cache objects. The final code is

function infinite_sim(dt=0.05, finaltime=5.0)
    ψ₀ = InfiniteMPS([ℂ^2], [ℂ^10])
    ψ₀, _ = find_groundstate(ψ₀, H₀, VUMPS())

    ψₜ = deepcopy(ψ₀)
    envs = environments(ψₜ, H₁)

    echos = [echo(ψₜ, ψ₀)]
    times = collect(0:dt:finaltime)

    for t in times[2:end]
        if t < 50dt # if t is sufficiently small, we increase the bond dimension
            ψₜ, envs = changebonds(ψₜ, H₁, OptimalExpand(; trscheme=truncdim(1)), envs)
        end
        (ψₜ, envs) = timestep(ψₜ, H₁, 0, dt, TDVP(), envs)
        push!(echos, echo(ψₜ, ψ₀))
    end

    return times, echos
end

infinite_sim (generic function with 3 methods)

This page was generated using Literate.jl.