Description
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
endIn 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
endSimilarly to the mpn functions I proposed, get_coherence can also be extended to multiple subsystems.
Alternate Solutions
No response
Additional Context
No response