Build the unitary matrix corresponding to one period of a periodic (Floquet) quantum circuit.

FloquetCircuits.jl

module FloquetCircuits

import Base: Matrix
using LinearAlgebra: I
using TensorOperations: tensorcontract

export Gate, FloquetCircuit, Matrix

"Quantum gate acting on some number of qubits."
struct Gate
    # sites the gate acts on
    sites::Vector{Int64}
    # unitary matrix
    matrix::Array{Complex{Float64}, 2}
    # default constructor
    function Gate(sites::Vector{Int64}, matrix::Array{Complex{Float64}, 2})
        # check consistency
        s = size(matrix)
        n = length(sites)
        # on-site dimension
        d = 2
        @assert s[1] == s[2] == d^n "Size of matrix is inconsistent with sites."
        new(sites, matrix)
    end
end

"Quantum circuit consisting of a series of gates applied to L qubits."
struct FloquetCircuit
    # number of sites
    L::Int64
    # gates
    gates::Vector{Gate}
    # default constructor
    function FloquetCircuit(L::Int64, gates::Vector{Gate})
        # check consistency
        for g in gates
            for x in g.sites
                @assert 1 ≤ x ≤ L "Gate sites inconsistent with L."
            end
        end
        new(L, gates)
    end
end

"Convert a quantum circuit to a dense unitary matrix."
function Matrix(circ::FloquetCircuit)
    # for use in einstein summation later
    alph = collect("abcdefghijklmnopqrstuvwxyz")
    ALPH = [uppercase(ltr) for ltr in alph]
    # on-site and full hilbert space dimensions
    d = 2
    L = circ.L
    D = d^L
    # initialize to the identity
    U = Matrix{Complex{Float64}}(I, D, D)
    # split the row index into a product over site indices
    sU = tuple(fill(d, L)..., D)
    vU = reshape(U, sU)
    # multiply by the gates in the circuit
    for g in circ.gates
        # number of sites gate acts on
        x = collect(g.sites)
        l = length(x)
        # split the row and column indices into products over site indices
        su = tuple(fill(d, 2*l)...)
        vu = reshape(g.matrix, su)
        # tensor contraction pattern
        Upat = alph[1:L+1]
        upat = [ALPH[x]; alph[x]]
        uUpat = copy(Upat); uUpat[x] = ALPH[x]
        # contract gate into the unitary
        vU .= tensorcontract(
            vu, tuple(upat...),
            vU, tuple(Upat...),
            tuple(uUpat...)
        )
    end
    U
end

end