Mixed fermion-spin-boson example

This example combines several features at once: fermionic, spin, and bosonic Hilbert spaces; custom constraints; and reduced density matrices.

Consider a spinful fermionic mode $f$, a localized spin-half $S$, and a bosonic mode $a$ with the Hamiltonian

\[H = \omega a^\dagger a + \Delta_f s_z + \Delta_S S_z + J\vec{s} \cdot \vec{S} + g\left(a^\dagger f_\downarrow^\dagger f_\uparrow + a f_\uparrow^\dagger f_\downarrow\right),\]

where $s_\alpha = \tfrac{1}{2} f^\dagger \sigma_\alpha f$ is the spin operator of the fermion. This model conserves $Q = a^\dagger a + s_z + S_z$ and the fermion number $N_f$. We will calculate the ground state in the $N_f = 1$, $Q = 0$ sector, and then the mutual information between different subsystems.

using FermionicHilbertSpaces, LinearAlgebraimport FermionicHilbertSpaces: particle_number, von_neumann_entropy@fermions f@spin S 1//2@boson aH_up = hilbert_space(f, [:↑])H_down = hilbert_space(f, [:↓])H_spin = hilbert_space(S)H_boson = hilbert_space(a, 4)  # Truncated boson with 4 levelsspaces = [H_up, H_down, H_spin, H_boson]Q = FilterConstraint(iszero ∘ sum;    maps = [state -> 1//2 * particle_number(state),    state -> -1//2 * particle_number(state),    state -> state.m,    state -> particle_number(state)],    spaces)Nf = NumberConservation(1; spaces = [H_up, H_down])H = tensor_product(spaces; constraint = Q * Nf)basisstates(H)

3-element basis of ConstrainedSpace(ProductSpace(32-dim, 3 factors), 3-dim):
 (01, ↓, 1)
 (01, ↑, 0)
 (10, ↓, 0)

The sector of interest has three states. Let's define the symbolic operators, build the Hamiltonian, and compute the ground state.

sz = (f[:↑]'*f[:↑] - f[:↓]'*f[:↓]) / 2sx = (f[:↑]'*f[:↓] + f[:↓]'*f[:↑]) / 2sy = (f[:↑]'*f[:↓] - f[:↓]'*f[:↑]) / (2im)s = [sx, sy, sz]ω, Δf, ΔS, J, g = 1.0, 0.8, 1.1, 0.25, 0.15ham = ω*a'*a + Δf*sz + ΔS*S[:z] + J*s'*S[:] +    g * (a'*f[:↓]'*f[:↑] + a*f[:↑]'*f[:↓])ψ = eigvecs(Hermitian(representation(ham, H, :dense)))[:, 1]

3-element Vector{ComplexF64}:
 -0.32327569786942045 + 0.0im
  -0.2865128353741618 + 0.0im
   0.9018886950909727 + 0.0im

Now compute the reduced density matrices and the mutual information between different subsystem pairs.

function mutual_info(ψ, HA, HB, H)    HAB = tensor_product(HA, HB)    S = von_neumann_entropy    S(partial_trace(ψ, H => HA)) + S(partial_trace(ψ, H => HB)) - S(partial_trace(ψ, H => HAB))endH_fermions = tensor_product(H_up, H_down; constraint = NumberConservation(1))subsystems = [(H_fermions, H_spin), (H_fermions, H_boson), (H_spin, H_boson)][mutual_info(ψ, HA, HB, H) for (HA, HB) in subsystems]

3-element Vector{Float64}:
 0.4302190755235034
 0.5322818481366896
 0.1374678871212185

The three entries are the mutual information of the fermion with the spin, of the fermion with the boson, and of the spin with the boson.