Utility functions for bosonic states

[labay11]

Problem Description

Define and correct some utility functions to get the mean photon number and coherence for bosonic systems.

In arithmetic_and_attributes.jl there is already a defined function to get the coherence for a quantum state. However, the current implementation is only valid for single systems. For composite systems the function computes an annihilation operator which is the product over all the dimensions, such definition is not correct as the total annihilation operator should be the kronecker product of a in each subsysems.

Proposed Solution

Here are some functions I use to calculate the mean photon number without the need to construct the operator. This is particularly faster for density matrix as you can skip the matrix product.

function mpn(ρ::QuantumObject{T,OperatorKetQuantumObject}) where {T}
    if length(ρ.dims) > 1
        throw(ArgumentError("Mean photon number not implemented for composite OPeratorKetQuantumObject"))
    end
    r = ρ.dims[1]
    return real(sum(ρ[(k-1)*r+k] * (k - 1) for k in 1:r))
end


function mpn(ψ::QuantumObject{T,KetQuantumObject}) where {T}
    if length(ψ.dims) == 1
        return real(mapreduce(k -> abs(ψ[k])^2 * (k - 1), +, 1:ρ.dims[1]))
    else
        R = CartesianIndices((ψ.dims...,))
        off = circshift(ψ.dims, 1)
        off[end] = 1

        x = 0.0im
        for j in R
            J = sum((Tuple(j) .- 1) .* off) + 1
            x += abs(ψ[J])^2 * prod(Tuple(j) .- 1)
        end
        return real(x)
    end
end

function mpn(ρ::QuantumObject{T,OperatorQuantumObject}) where {T}
    if length(ρ.dims) == 1
        return real(mapreduce(k -> abs(ρ[k, k]) * (k - 1), +, 1:ρ.dims[1]))
    else
        R = CartesianIndices((ρ.dims...,))
        off = circshift(ψ.dims, 1)
        off[end] = 1

        J = 1
        x = 0.0im
        for j in R
            # convert the cartesian index j into an array index J
            J = sum((Tuple(j) .- 1) .* off) + 1
            x += ρ[J, J] * prod(Tuple(j) .- 1)
        end
        return real(x)
    end
end

In the case of the get_coherent functions I would propose to divide into two functions: get_coherence which returns expect(a, *) and remove_coherence which returns the state with D' * ψ applied.

function get_coherence(ψ::QuantumObject{<:AbstractArray,KetQuantumObject,1})
    return mapreduce(n -> sqrt(n-1) * ψ.data[n] * conj(ψ.data[n-1]), +, 2:ψ.dims[1])
end

function get_coherence(ρ::QuantumObject{T,OperatorQuantumObject,1})
    return mapreduce(n -> sqrt(n-1) * ρ.data[n, n-1], +, 2:ψ.dims[1])
end

function remove_coherence(ψ::QuantumObject{<:AbstractArray,KetQuantumObject})
    α = get_coherence(ψ)
    a = mapreduce(dim -> destroy(dim), kron, ψ.dims)
    D = exp(α * a' - conj(α) * a)

    return D' * ρ * D
end

Similarly to the mpn functions I proposed, get_coherence can also be extended to multiple subsystems.

Alternate Solutions

No response

Additional Context

No response

[albertomercurio]

Hi @labay11,

Thank you for suggesting new features.

I like the idea of avoiding to define the bosonic operators when computing the coherence of the state, by directly iterating over the elements of the array. There are a few comments about the functions you defined.

Perhaps it is better to change the name of the mpn function to something more specific and readable. Maybe mean_occupation, or mean_population. I would avoid to use the photon word, because they represent generic bosons as phonons, magnons and many others. If you have other suggestions about the name, please write yours here.

Then, I would include the options to also choose for which subsystem we want to compute the mean occupation or the coherence. Immagine to have many bosons, and I only want to compute the coherence of the second one only, or the third and the fourth.

The definition of your mpn functions are very smart, even if I think we can improve them a bit. For example, we can use the @inbounds macro to avoid bounds checking when indexing over the array. Also, the following code should be a bit more efficient

x = sum(R) do j
    j_tuple = Tuple(j) .- 1
    J = dot(j_tuple, off) + 1
    abs2(ψ[J]) * prod(j_tuple)
end

and also the other implementations.

Overall, I think this is a good idea. Do you want to contribute with a Pull Request?

[ytdHuang]

Does Python qutip also have a similar function ?

If it does, we can make the function name same as theirs.

[labay11]

@albertomercurio I like the idea of being able to calculate the mean value in just one subsystem. Also the name mean_occupation sounds good to me.

@ytdHuang I think it is not implemented in QuTiP.