Boundary-driven spin transport

A closed spin chain merely redistributes magnetisation already present in the system. To observe sustained transport, we instead attach two reservoirs that favour opposite edge polarisations. The left reservoir tries to polarise the first spin upward, the right reservoir tries to polarise the final spin downward, and the chain must continuously carry spin between them.

This boundary-driven setup is a standard nonequilibrium open-system model. It produces three readable signatures:

• a spin current that grows from zero,
• a magnetisation profile across the chain,
• an approximately uniform current once the bulk stops accumulating magnetisation.

The implementation is correspondingly compact:

1. define the XXZ Hamiltonian as a physical OpSum;
2. encode the two reservoirs with four boundary jump_ops;
3. build one Liouvillian MPO;
4. evolve the vectorized density matrix with two-site TDVP.

Model and physical scales

We use the spin-$1/2$ XXZ Hamiltonian

\[H = J\sum_{j=1}^{N-1} \left( S_j^x S_{j+1}^x + S_j^y S_{j+1}^y + \Delta S_j^z S_{j+1}^z \right).\]

The $XY$ terms move spin between neighbouring sites, while $\Delta$ sets the interaction anisotropy. The density matrix obeys

\[\frac{d\rho}{dt} = -i[H,\rho] + \sum_k \mathcal{D}[L_k]\rho ,\]

with opposing boundary reservoirs

\[L_{1,+}=\sqrt{\Gamma(1+\mu)}\,S_1^+, \qquad L_{1,-}=\sqrt{\Gamma(1-\mu)}\,S_1^-,\]

\[L_{N,+}=\sqrt{\Gamma(1-\mu)}\,S_N^+, \qquad L_{N,-}=\sqrt{\Gamma(1+\mu)}\,S_N^-.\]

Here, the parameter $\mu$ controls the bias between the left and right reservoirs, with its range restricted to $-1 \leq \mu \leq 1$ . An isolated left spin would approach $\langle S_1^z\rangle=\mu/2$; an isolated right spin would approach $\langle S_N^z\rangle=-\mu/2$. Their disagreement drives the chain away from equilibrium.

Because the bulk Hamiltonian conserves total $S^z$, the local magnetisation $m_j=\langle S_j^z\rangle$ obeys a lattice continuity equation. The bond current consistent with $H$ is

\[\mathcal{J}_j = J\left( S_j^x S_{j+1}^y - S_j^y S_{j+1}^x \right).\]

We start from $|\mathrm{Dn}\cdots\mathrm{Dn}\rangle$ and report the final magnetisation profile together with the mean and standard deviation of the bond currents.

Parameters and operators

using Printf
using Statistics: mean, std
using ITensors
using ProcessTensors

const N = 4
const J = 1.0
const Δ = 0.5
const Γ = 1.0
const μ = 0.4
const dt = 0.1
const final_time = 8.0
const maxdim = 64
const cutoff = 1e-10

physical_sites = siteinds("S=1/2", N)
liouville_sites = liouv_sites(physical_sites)

hamiltonian = let H = OpSum()
    for j in 1:(N - 1)
        H += J, "Sx", j, "Sx", j + 1
        H += J, "Sy", j, "Sy", j + 1
        H += J * Δ, "Sz", j, "Sz", j + 1
    end
    H
end

sum(
  1.0 Sx(1,) Sx(2,)
  1.0 Sy(1,) Sy(2,)
  0.5 Sz(1,) Sz(2,)
  1.0 Sx(2,) Sx(3,)
  1.0 Sy(2,) Sy(3,)
  0.5 Sz(2,) Sz(3,)
  1.0 Sx(3,) Sx(4,)
  1.0 Sy(3,) Sy(4,)
  0.5 Sz(3,) Sz(4,)
)

Package jump tuples use the dissipative rate itself, not its square root, so the factors written under the square roots above enter as Γ * (1 ± μ).

jump_operators = [
    (Γ * (1 + μ), "S+", 1),
    (Γ * (1 - μ), "S-", 1),
    (Γ * (1 - μ), "S+", N),
    (Γ * (1 + μ), "S-", N),
]

initial_state = MPS(physical_sites, fill("Dn", N))
initial_density = to_dm(initial_state)
initial_density_liouville =
    to_liouville(initial_density; sites=liouville_sites)

liouvillian = liouvillian_mpo(
    hamiltonian,
    liouville_sites;
    jump_ops=jump_operators,
)

magnetisation_ops = [
    let observable = OpSum()
        observable += 1.0, "Sz", j
        to_liouville(MPO(observable, physical_sites); sites=liouville_sites)
    end for j in 1:N
]

current_ops = [
    let observable = OpSum()
        observable += J, "Sx", j, "Sy", j + 1
        observable += -J, "Sy", j, "Sx", j + 1
        to_liouville(MPO(observable, physical_sites); sites=liouville_sites)
    end for j in 1:(N - 1)
]

@assert isapprox(real(tr(initial_density)), 1.0; atol=1e-12)

Liouville-space TDVP

Vectorization turns the master equation into

\[\frac{d}{dt}|\rho(t)\rangle\rangle = \mathcal L|\rho(t)\rangle\rangle .\]

The Liouvillian is time-independent, so we build it once and advance with two-site TDVP. The evolution argument is the real interval dt because liouvillian_mpo already includes the Hamiltonian factor $-i$.

nsteps = round(Int, final_time / dt)
times = collect(range(0.0; step=dt, length=nsteps + 1))
@assert isapprox(nsteps * dt, final_time; atol=100eps(Float64))

Time evolution

At each stored time we record the mean bond current. The full magnetisation and current profiles are evaluated once at the end.

trajectory = let
    density = copy(initial_density_liouville)
    mean_current = Float64[]
    trace_errors = Float64[]
    bond_dimensions = Int[]

    for step in eachindex(times)
        density_trace = tr(to_hilbert(density))
        bond_currents = [
            real(inner(observable, density) / density_trace)
            for observable in current_ops
        ]
        push!(mean_current, mean(bond_currents))
        push!(trace_errors, abs(density_trace - 1))
        push!(bond_dimensions, maxlinkdim(density))

        step == length(times) && continue

        density = tdvp(
            liouvillian,
            dt,
            density;
            time_step=dt,
            nsite=2,
            maxdim=maxdim,
            cutoff=cutoff,
            outputlevel=0,
        )
    end

    density_trace = tr(to_hilbert(density))
    magnetisation_profile = [
        real(inner(observable, density) / density_trace)
        for observable in magnetisation_ops
    ]
    current_profile = [
        real(inner(observable, density) / density_trace)
        for observable in current_ops
    ]

    (
        mean_current=mean_current,
        magnetisation_profile=magnetisation_profile,
        current_profile=current_profile,
        trace_errors=trace_errors,
        bond_dimensions=bond_dimensions,
    )
end

(mean_current = [0.0, 0.0006034540455180815, 0.00218673324875952, 0.004462003760058623, 0.007202161656387689, 0.01023036661911011, 0.013411980365508407, 0.016647154372545783, 0.019864072095275232, 0.023012540129862025, 0.026060432435469693, 0.028989365684385048, 0.031791276814156284, 0.03446570596551502, 0.037017483130036145, 0.03945477981082928, 0.041787616877850244, 0.04402667097362447, 0.04618236124535408, 0.04826417744965376, 0.050280214267809244, 0.05223715924925327, 0.05413959442014465, 0.055989934426488246, 0.05778908215115816, 0.059536301665408686, 0.06122941806386099, 0.06286505594346854, 0.06443890181020474, 0.06594597766894793, 0.06738091491982577, 0.06873821944665687, 0.07001252043730896, 0.0711988174408537, 0.07229285004759502, 0.07329052395244941, 0.07418889747061108, 0.07498620053856632, 0.07568154195346775, 0.07627510752910749, 0.07676817062593284, 0.07716307443324015, 0.07746318860810165, 0.07767284323455009, 0.07779724332298561, 0.07784236723016807, 0.07781485245176258, 0.0777218722301836, 0.07757100633889615, 0.0773701092592063, 0.07712717876597552, 0.07685022769373524, 0.07654716137342961, 0.07622566292132453, 0.07589308823404099, 0.07555637220530351, 0.07522194733842782, 0.07489567559077485, 0.07458279395867824, 0.07428787399926542, 0.07401479519395131, 0.07376673179118863, 0.0735461525265459, 0.07335483240876244, 0.07319387558278223, 0.07306374813579958, 0.07296431960031506, 0.07289491182866638, 0.07285435386549313, 0.07284104142658221, 0.07285299960259206, 0.07288794744188452, 0.07294336312551186, 0.07301654852639268, 0.07310469204086285, 0.07320492869096545, 0.07331439661690596, 0.07343028920793011, 0.07354990225346807, 0.07367067563183764, 0.07379022918843257], magnetisation_profile = [0.15573004469954754, -0.023890177453871505, -0.04634814523298402, -0.17127744430778194], current_profile = [0.08535304082618868, 0.07614104505452253, 0.05987660168458651], trace_errors = [0.0, 2.529159297992578e-6, 9.328098711014476e-6, 1.0359841089879751e-5, 1.3565096696810386e-5, 1.595507901241966e-5, 2.11203533544424e-5, 2.5193928733235893e-5, 2.5209124003167105e-5, 2.525457484737059e-5, 2.5356507114370623e-5, 2.5547902431424836e-5, 2.5864620007769102e-5, 2.6341226087045712e-5, 2.700208318054349e-5, 2.785921085143811e-5, 2.8911045958130817e-5, 3.0142065651464954e-5, 3.15240870029961e-5, 3.3018857580513325e-5, 3.4581398736427814e-5, 3.5373062948601255e-5, 3.537306294926738e-5, 3.5373062949267376e-5, 3.537306295037759e-5, 3.537306295015555e-5, 3.537306295082168e-5, 3.5373062950155554e-5, 3.5373062950599636e-5, 3.5373062949267376e-5, 3.5373062949045345e-5, 3.537306294704695e-5, 3.5373062947046936e-5, 3.5340125748817606e-5, 3.487274482337774e-5, 3.460024920578903e-5, 3.460024920556699e-5, 3.4600249204234716e-5, 3.4600249204012665e-5, 3.46002492046788e-5, 3.460024920512288e-5, 3.460024920467879e-5, 3.46002492062331e-5, 3.460024920601105e-5, 3.460024920578901e-5, 3.460024920689923e-5, 3.460024920689923e-5, 3.460024920712128e-5, 3.460024920645515e-5, 3.460024920623311e-5, 3.4600249206011066e-5, 3.460024920645516e-5, 3.4600249206899244e-5, 3.460024920667719e-5, 3.460024920645515e-5, 3.460024920667719e-5, 3.460024920645515e-5, 3.460024920689924e-5, 3.46002492066772e-5, 3.460024920645516e-5, 3.460024920601107e-5, 3.460024920667721e-5, 3.460024920778744e-5, 3.460024920778745e-5, 3.460024920800949e-5, 3.460024920867563e-5, 3.460024921000791e-5, 3.460024921134018e-5, 3.4600249210674037e-5, 3.460024921089608e-5, 3.4600249211118126e-5, 3.460024921134017e-5, 3.460024921222834e-5, 3.4600249212894476e-5, 3.46002492140047e-5, 3.4600249214226743e-5, 3.460024921267243e-5, 3.460024921311653e-5, 3.460024921311652e-5, 3.460024921245038e-5, 3.460024921222834e-5], bond_dimensions = [1, 5, 5, 5, 7, 8, 10, 11, 11, 12, 13, 14, 15, 15, 15, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 15, 15, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16])

Final diagnostics

mean_bond_current = mean(trajectory.current_profile)
std_bond_current = std(trajectory.current_profile; corrected=false)
left_right_magnetisation =
    first(trajectory.magnetisation_profile) - last(trajectory.magnetisation_profile)
max_trace_error = maximum(trajectory.trace_errors)
max_bond_dimension = maximum(trajectory.bond_dimensions)

@assert all(isfinite, trajectory.mean_current)
@assert all(isfinite, trajectory.magnetisation_profile)
@assert all(isfinite, trajectory.current_profile)
@assert all(value -> -0.5 - 1e-8 ≤ value ≤ 0.5 + 1e-8, trajectory.magnetisation_profile)
@assert max_trace_error < 1e-2

println("Boundary-driven XXZ spin chain")
@printf("  N=%d, J=%.2f, Δ=%.2f, Γ=%.2f, μ=%.2f\n", N, J, Δ, Γ, μ)
@printf("  simulated to t = %.2f\n", final_time)
@printf("  final mean bond current = %.6f\n", mean_bond_current)
@printf("  final bond-current std  = %.3e\n", std_bond_current)
@printf("  final ⟨S₁ᶻ⟩ − ⟨Sₙᶻ⟩   = %.6f\n", left_right_magnetisation)
@printf("  max trace error = %.3e, max bond dim = %d\n", max_trace_error, max_bond_dimension)

Boundary-driven XXZ spin chain
  N=4, J=1.00, Δ=0.50, Γ=1.00, μ=0.40
  simulated to t = 8.00
  final mean bond current = 0.073790
  final bond-current std  = 1.053e-02
  final ⟨S₁ᶻ⟩ − ⟨Sₙᶻ⟩   = 0.327007
  max trace error = 3.537e-05, max bond dim = 16

Transport response and numerical interpretation

Starting from all Dn, every site begins at $\langle S_j^z\rangle=-1/2$. The left reservoir injects upward polarisation, so site $1$ rises toward a positive value, while the right reservoir keeps site $N$ near a negative value. This left-to-right bias drives a nonzero spin current through the XXZ chain. The mean bond current rises from zero during the transient and saturate to different values for each $Delta$ value in the later time, as it nears the steady state.

At late times the magnetisation profile settles and the bond currents become more uniform. Continuity then requires that spin entering one bond leave through the next. The resulting state is a nonequilibrium steady state maintained by the reservoirs, not an equilibrium state of the XXZ Hamiltonian.

Sources of numerical error

• TDVP / truncation error: two-site updates allow the Liouville-MPS bond dimension to grow. cutoff and maxdim control the discarded weight.
• Finite timestep: the real-time TDVP step approximates the short-time propagator generated by $\mathcal L$.
• Finite simulation window: final_time may be shorter than the time needed for a fully uniform bond-current profile.

Trace drift is a cheap diagnostic for the density-matrix evolution. The mean and standard deviation of the final bond currents give a compact check of spatial uniformity at the simulated time.