Feature: Floquet spectrum

[oameye]

Is there some interest that I add an example in docs or even a function to construct the quasienergy operator in Floquet theory. Once diagonalized it provides the Floquet spectrum of the time-periodic Hamiltonian.

Here is some code example to create a the basis

struct FloquetBasis <: Basis
    shape::Vector{Int}
    m::Int
    FloquetBasis(m::Int) = new(Int[2*m+1], m) # Constructor
end

function QuantumOptics.number(b::FloquetBasis)
    N = first(b.shape)
    diag = ComplexF64.(N:-1:1)
    data = spdiagm(0 => diag)
    SparseOperator(b, data)
end

function floquet_cell(v::Ket)
    fock, floquet = v.basis.bases
    m̄ = identityoperator(fock) ⊗ number(floquet)
    real(expect(m̄, v))
end

function photon_number(v::Ket)
    fock, floquet = v.basis.bases
    n̄ = number(fock) ⊗ identityoperator(floquet)
    real(expect(n̄, v))
end

function parity(v::Ket)
    fock, floquet = v.basis.bases
    diag_vec = Complex[(-1)^(i-1) for i in 1:first(fock.shape)]
    p = Operator(fock, Diagonal(diag_vec));    
    p̄ = p ⊗ identityoperator(floquet)
    round(Int,real(expect(p̄, v)))
end
function get_quasienergy_operator(b; ω, ω₀=1.0, F=0.0, α=0.0, f::Function=isinteger)
    dim_floquet = b.shape[2]
    m = (dim_floquet-1)/2

    # Fundamental operators
    a = destroy(b.bases[1])
    at = create(b.bases[1])
    n = number(b.bases[1])
    I = identityoperator(b.bases[1])
    zero = spzeros(ComplexF64, size(I))

    # Hamiltonians
    H₀ᵃ(ω) = -((ω - ω₀) - 3*α/(4*ω₀))*n - (F/√(8*ω₀))*(a + at) + (3*α/(8*ω₀^2))*at*at*a*a
    H₋₂ᵃ(ω) = - (F/√(8*ω₀))*a + (3*α/(4*ω₀^2))*a*a + (α/(4*ω₀^2))*at*a*a*a
    H₂ᵃ(ω) =  - (F/√(8*ω₀))*at + (3*α/(4*ω₀^2))*at*at + (α/(4*ω₀^2))*at*at*at*a
    H₋₄ᵃ(ω) = (α/(16*ω₀^2))*a*a*a*a
    H₄ᵃ(ω) = (α/(16*ω₀^2))*at*at*at*at

    matrix_of_matrix = [zero for i in 1:dim_floquet, j in 1:dim_floquet]

    filter_f(lin_idxs, M=matrix_of_matrix) = LinearIndices(M)[filter(idx -> !any(iszero.(M[get_diag_idxs(idx)])), CartesianIndices(M)[lin_idxs])]

    matrix_of_matrix[diagind(matrix_of_matrix)] = [f(i) ? getfield(H₀ᵃ(ω) + i*ω*I, :data) : zero for i in -m:m]
    matrix_of_matrix[filter_f(diagind(matrix_of_matrix,-2), matrix_of_matrix)] .= [H₂ᵃ(ω).data]
    matrix_of_matrix[filter_f(diagind(matrix_of_matrix,2), matrix_of_matrix)] .= [H₋₂ᵃ(ω).data]
    matrix_of_matrix[filter_f(diagind(matrix_of_matrix,-4), matrix_of_matrix)] .= [H₄ᵃ(ω).data]
    matrix_of_matrix[filter_f(diagind(matrix_of_matrix,4), matrix_of_matrix)] .= [H₋₄ᵃ(ω).data]
    matrix = reduce(vcat, [reduce(hcat, matrix_of_matrix[i, :]) for i in 1:dim_floquet])

    return SparseOperator(b, matrix)
end

Qᵃ =  get_quasienergy_operator_a(basis; ω=1.02, ω₀=1, F=0.0001, α=0.01)
Qᵃ.data

qutip floquet

[Krastanov]

Apologies for the slow answer. There is certainly a lot of interest in this being added. Please feel free to create a PR with this, even if it is in an unfinished state (maybe another community member will be able to help finish it).