A new Julia package for Gaussian Quantum Optics [$1600]

[Krastanov]

Bug bounty logistic details (click to expand)

To claim exclusive time to work on this bounty either post a comment here or message [skrastanov@umass.edu](mailto:skrastanov@umass.edu) with: - your name - github username - **(optional)** a brief list of previous pertinent projects you have engaged in Currently the project is claimed by `no one` until `...`. If you want to, you can work on this project without making a claim, however claims are encouraged to give you and other contributors peace of mind. Whoever has made a claim takes precedence when solutions are considered. You can always propose your own funded project, if you would like to contribute something of value that is not yet covered by an official bounty.

A new Julia package for Gaussian Quantum Optics [$1600]

A package for working with quantum modes and linear quantum optics (quadratic Hamiltonians), implementing the entirety of the Gaussian Quantum Information review paper. The package needs to have some baseline documentation and tests, to conform at least partially to the QuantumInterface.jl API, and to not have architectural limitations that make impossible its future use with GPU arrays or autodifferentiation.

On successful completion of this bounty, future related bounties can/will be proposed.

Required skills: knowledge of Gaussian quantum optics

Reviewer: Stefan Krastanov

Duration: 3 months

Payout procedure:

The Funding for these bounties comes from the National Science Foundation and from the NSF Center for Quantum Networks. The payouts are managed by the NumFOCUS foundation and processed in bulk once every two months. If you live in a country in which NumFOCUS can make payments, you can participate in this bounty program.

Click here for more details about the bug bounty program.

[apkille]

Hi Stefan! I would definitely be interested in completing this bounty. I can start it in the next few weeks, once I finish my SciML integration project with QuantumOptics.jl.

I've been wanting to create a package from the ground up, and I have a research background in quantum optics, so I'm quite enthusiastic about this bounty. Since I've been working with QuantumOptics for the package to have use with GPU arrays, autodiff, sensitivity analysis, etc., I know how to avoid any related architectural limitations in this new package. Plus, I am already familiar with QuantumInterface, having worked on QuantumSymbolics and QuantumOptics.

You and I can also flesh out a lot of the details of this new package at QNumerics, if you would be interested.

[Krastanov]

Hi, Andrew! That would be pretty great, and you have a good track record with this type of projects. My only concern is having you spread too thin. Would the following be agreeable:

1. Before officially reserving this bounty for you, waiting for the completion of your QuantumSymbolics summer project (acknowledging that some of the time sink there is me being slow to do some reviews)
2. Independently of point 1, waiting until the end of the summer school (qnumerics.org for anyone else reading through this conversation) for officially reserving this bounty so that other interested contributors have a chance to engage.

This is a big enough project that we are happy to fund multiple contributors engaging on it together without lowering individual bounty amounts.

[apkille]

Sure, that would be OK. Of course I wouldn't start this project until I finished the QuantumSymbolics project, but just for the record, the QS project is wrapping up pretty soon (only 1 or 2 small PRs left). But yes, I'm happy to be accommodating in any way :)

[Krastanov]

@apkille , in preparation for potentially doing this bounty, could you send me a rough project description, similar to the GSoC style description for your previous project?

[apkille]

Sure, I'll email it to you.

[apkille]

To add context for people in the future: this package was discussed at QNumerics between Christian, Gabe, and myself (sorry I can't ping you two yet because I don't know your GitHub usernames!). Here's a very rough sketch of what an interface could look like for such a package:

import QuantumOpticsBase: coherentstate
import QuantumInterface: AbstractOperator
import StaticArrays: SMatrix, SVector

abstract type BosonicState end
abstract type GaussianOperator <: AbstractOperator end

"""
    GaussianState(means::V, covariance::T)

Arbitrary Gaussian state characterized by mean value vector `mean` and covariance matrix `covariance`.
"""
struct GaussianState{V,T} <: BosonicState
    mean::V
    covariance::T
end

"""
    vacuumstate([Tm = SVector{2}, Tc = SMatrix{2,2}])

Create a vacuum state, a Gaussian state with zero photons.
"""
function vacuumstate(::Type{Tm}, ::Type{Tc}) where {Tm, Tc}
    mean = Tm([0.0, 0.0])
    covar = Tc([1.0 0.0; 0.0 1.0])
    return GaussianState(mean, covar)
end
vacuumstate() = vacuumstate(SVector{2}, SMatrix{2,2})

"""
    coherentstate([Tm = SVector{2}, Tc = SMatrix{2,2},] alpha::N)

Create a coherent state, the single-photon quantum analogue of monochromatic classical electromagnetic field, where `alpha` is the eigenvalue of
the bosonic annihilation operator.
"""
function coherentstate(::Type{Tm}, ::Type{Tc}, alpha::A) where {Tm, Tc, A<:Number}
    q = real(alpha)
    p = imag(alpha)
    mean = Tm([q, p])
    covar = Tc([1.0 0.0; 0.0 1.0])
    return GaussianState(mean, covar)
end
coherentstate(alpha::A) where {A<:Number} = coherentstate(SVector{2}, SMatrix{2,2}, alpha)

"""
    DisplaceOperator(alpha::N)

Displacement operator, a Gaussian unitary defined by the relation `D(α)|0⟩ = |α⟩`.
"""
struct DisplaceOperator{N} <: GaussianOperator
    alpha::N
end

function Base.:(*)(x::DisplaceOperator{<:Number}, y::GaussianState{V,T}) where {V,T}
    mean = y.mean
    if length(mean) > 2
        throw(ArgumentError("You are applying a single-mode displacement operator to a multi-mode Gaussian state."))
    else
        covar = y.covariance
        q = y.mean[1]
        p = y.mean[2]
        return GaussianState(V([q+sqrt(2)*real(x.alpha),p+sqrt(2)*imag(x.alpha)]), covar)
    end
end

function Base.:(*)(x::DisplaceOperator{AbstractVector}, y::GaussianState{V,T}) where {V,T}
    # insert code for applying multi-mode displacement operator to multi-mode Gaussian state
end

[apkille]

Several things here to note:

• We should avoid code reuse as much as possible. Here, there are many overlapping capabilities with QuantumOptics.jl (a state vector package). QO.jl has a base functionality package called QOBase.jl, as it is a very large package and convenient to keep basic definitions somewhere else (similar to DiffEqBase.jl and DifferentialEquations.jl). We can reuse a lot of their functions.
• We should use QuantumInterface.jl (as mentioned by @krastanov), as we want easy extension capabilities with other QuantumSavory packages and QuantumOptics. Edit: IMO, it might be good to add some abstract type definitions to QuantumInterface to handle this specific package. For instance, the only abstract base class for quantum states at the moment is StateVector, which isn't useful for our purposes, as it is intended for types that wrap around a basis and abstract array.

Edit: Regarding the last bullet, StateVector probably can be used in our case, as it should be considered as a vector in a Hilbert space, rather than a state vector representation in QO.jl

[apkille]

Here's a very rough sketch of what an interface could look like for such a package:

Just to explain my thoughts more about the sketch:

• IMO, we should have a basic concrete type for Gaussian states (here it is called GaussianState) that wraps around its two statistical moments, which fully characterize such an object: the mean vector and covariance matrix. We can create specific Gaussian states with functions, e.g., create a vacuum state by calling vacuumstate() or vacuumstate(::Type, ::Type), the latter if we care about what type of array we are using; both commands return a GaussianState object with its statistical moments corresponding to a vacuum state:

  julia> vacuumstate()
  GaussianState{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}([0.0, 0.0], [1.0 0.0; 0.0 1.0])

• Then, we create specific Gaussian operators as subtypes of AbstractOperator from QuantumInterface. For instance, we have above DisplaceOperator, which can hold, say a single complex number, which corresponds to an operator that's applied to a single mode. Or, a DisplaceOperator can hold an array of complex numbers, which corresponds to an operator applied to a multimode system:

  julia> D = DisplaceOperator(0.5 + 0.5im)
  DisplaceOperator{ComplexF64}(0.5 + 0.5im)

• We can also, for instance, handle Gaussian operations on Gaussian states with *:

  julia> D*v
  GaussianState{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}([0.7071067811865476, 0.7071067811865476], [1.0 0.0; 0.0 1.0])

[apkille]

I'm happy to hear anyone's thoughts on this approach and discuss this further. Some more food for thought:

• We could have nice plotting extensions with Makie.jl and Plots.jl, for instance, to easily create phase-space diagrams of quantum states.
• I may be biased, but it would be really nice (and actually fairly easy) to have express from QuantumSymbolics be able to convert symbolic representations to this package.
• Conversion capabilities to the QuantumOptics state vector representations is also quite important IMO.

[jgr-rgb]

Hi Andrew! I completely agree with your thinking, and the way you have laid it out is really helpful as someone with minimal julia experience. You can find the short overleaf document that we started last night here which sets the groundwork for general functionality that the toolbox needs (at least initially). It seems that we are thinking along the same lines structurally

[jgr-rgb]

Since any gaussian operator can be expressed as a symplectic matrix which applies to the covariance matrix, I think that we can just wrap these symplectic matrices in a julia function, which would then be applied to the gaussian state object that you described above. Then, once we have a final covariance matrix, we can do routines like you mentioned for plotting. I completely agree that symbolic representations would be VERY helpful, so I think it would be good if we plan on that functionality

[jgr-rgb]

We are going to try and take a crack right now on building some of the functions based on your starting code above, and will report back with any results (maybe tomorrow in person?)

[apkille]

Since any gaussian operator can be expressed as a symplectic matrix which applies to the covariance matrix, I think that we can just wrap these symplectic matrices in a julia function, which would then be applied to the gaussian state object that you described above. Then, once we have a final covariance matrix, we can do routines like you mentioned for plotting. I completely agree that symbolic representations would be VERY helpful, so I think it would be good if we plan on that functionality

Oh that's a great idea. So maybe we could be more general and define a GaussianOperator concrete type as

struct GaussianOperator{T} <: AbstractOperator
    displace::V
    symplectic::T
end

Then create, say a displacement operator with the function displace(alpha::T, modes::N) -> GaussianOperator (which we can import from QOBase.jl).

We are going to try and take a crack right now on building some of the functions based on your starting code above, and will report back with any results (maybe tomorrow in person?)

Yes, definitely. Feel free to ask questions on here in the meantime. If you encounter errors, it could be that I messed something up in that code sketch.

Edit: just to be clear: the displace field in GaussianOperator corresponds to d∈R²ᴺ in https://arxiv.org/abs/1110.3234, not the displacement operator QuantumOpticsBase.displace directly. Maybe there's a better name we could choose for the field.

[apkille]

Maybe something like this:

import QuantumOpticsBase: coherentstate, displace
import StaticArrays: SMatrix, SVector
import LinearAlgebra: I

"""
    GaussianState(means::V, covariance::T)

Arbitrary Gaussian state characterized by mean value vector `mean` and covariance matrix `covariance`.
"""
struct GaussianState{V,T}
    mean::V
    covariance::T
end

"""
    vacuumstate([Tm = SVector{2}, Tc = SMatrix{2,2}])

Create a vacuum state, a Gaussian state with zero photons.
"""
function vacuumstate(::Type{Tm}, ::Type{Tc}) where {Tm, Tc}
    mean = Tm([0.0, 0.0])
    covar = Tc([1.0 0.0; 0.0 1.0])
    return GaussianState(mean, covar)
end
vacuumstate() = vacuumstate(SVector{2}, SMatrix{2,2})

"""
    coherentstate([Tm = SVector{2}, Tc = SMatrix{2,2},] alpha::N)

Create a coherent state, the single-photon quantum analogue of monochromatic classical electromagnetic field, where `alpha` is the eigenvalue of
the bosonic annihilation operator.
"""
function coherentstate(::Type{Tm}, ::Type{Tc}, alpha::A) where {Tm, Tc, A<:Number}
    q = real(alpha)
    p = imag(alpha)
    mean = Tm([q, p])
    covar = Tc([1.0 0.0; 0.0 1.0])
    return GaussianState(mean, covar)
end
coherentstate(alpha::A) where {A<:Number} = coherentstate(SVector{2}, SMatrix{2,2}, alpha)

"""
    GaussianOperator(displacement::V, symplectic::T)

Arbitrary Gaussian operator characterized by displacement vector `displacement` and symplectic form `symplectic`.
"""
struct GaussianOperator{V,T}
    displacement::V
    symplectic::T
end

"""
    displace(alpha<:Number, modes<:Int)

Create a displacement operator on N modes, a Gaussian unitary defined by the relation `D(α)|0⟩ = |α⟩`.
"""
function displace(alpha::AbstractVector{Complex{T}}) where {T}
    modes = length(alpha)
    d = sqrt(2)*reinterpret(T, alpha)
    return GaussianOperator(d, Matrix(I, 2*modes, 2*modes))
end

Base.:(*)(x::GaussianOperator, y::GaussianState) = GaussianState(x.symplectic*y.mean + x.displacement, x.symplectic*y.covariance*transpose(x.symplectic))

so you have

julia> d = displace([0.5+0.5im])
GaussianOperator{Vector{Float64}, Matrix{Bool}}([0.7071067811865476, 0.7071067811865476], Bool[1 0; 0 1])

julia> v = vacuumstate()
GaussianState{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}([0.0, 0.0], [1.0 0.0; 0.0 1.0])

julia> d*v
GaussianState{Vector{Float64}, Matrix{Float64}}([0.7071067811865476, 0.7071067811865476], [1.0 0.0; 0.0 1.0])

[jgr-rgb]

Exactly! We spent the end of the hack-a-thon discussing how we would have the user define an initial multimode state, and the idea that we landed on is they would do so in the same convention as the QuantumOptics.jl package by writing something like: vacuumState \otimes vacuumstate \otimes coherentState ... but on the backend the \otimes would be performing the direct sum of the covariance matrices and appending the mean vectors. I think we probably want to generalize these initial states in some way since we would also want to user to be able to use squeezed vacuum or squeezed coherent states. Then, once the initial "GaussianState" object is defined by the user, symplectic matrix operations can just be applied in the same way as you defined above. With that functionality out of the way we simple remain with writing operations for representing the gaussian state, say via the Wigner function via a plot, or performing detection on some of the modes. I think that it would also be valuable to consider symbolics right from the start so that the user can define symbolic parameters for the initial state and the operators, so maybe you could give us a crash course in what that would take?

[jgr-rgb]

Also, Stefan told us that a slightly more formal proposal process would be warranted to reserve this bounty. I think we have all of the pieces for that proposal, but it is just a matter of gathering them, so let's i'll add these details to the overleaf document that I linked above

[apkille]

Exactly! We spent the end of the hack-a-thon discussing how we would have the user define an initial multimode state, and the idea that we landed on is they would do so in the same convention as the QuantumOptics.jl package by writing something like: vacuumState \otimes vacuumstate \otimes coherentState ... but on the backend the \otimes would be performing the direct sum of the covariance matrices and appending the mean vectors.

Sure that's probably fine. Or we could simply use +.

I think we probably want to generalize these initial states in some way since we would also want to user to be able to use squeezed vacuum or squeezed coherent states.

Maybe we could make the GaussianState structs mutable, so we can edit the statistical moments of initial states (if we want to squeeze them, etc).

Then, once the initial "GaussianState" object is defined by the user, symplectic matrix operations can just be applied in the same way as you defined above. With that functionality out of the way we simple remain with writing operations for representing the gaussian state, say via the Wigner function via a plot, or performing detection on some of the modes. I think that it would also be valuable to consider symbolics right from the start so that the user can define symbolic parameters for the initial state and the operators, so maybe you could give us a crash course in what that would take?

Sure I can show you the Fock stuff in QSymbolics during breaks or lunch.

[jgr-rgb]

@apkille @Krastanov

The group of us that has been working on the proposal for this bounty started something of a working GitHub here that describes the philosophy around how the toolbox could work and some basic code that builds upon Andrew's code from above. We also have discussed in separate conversations that a better path forward might be to split this bounty into smaller bounties since we have a number of individuals that are interested in this project (including @11of12 @cjcarver Topher and Ajit). @Krastanov mentioned that one of these smaller bounties could be to write the rules for Gaussian state to enable integration with QuantumSymbolics.jl. What would other smaller bounties look like?

[jgr-rgb]

Oh, and by "working Github" I dont mean function, but rather a space to work on this project.

[jgr-rgb]

@Topher37

[Krastanov]

Thanks folks! I have been getting a bunch of interested folks in this, so I would like to proceed a bit more officially on it:

• if someone is interested to work on this, could you send me a proposal (in the style of the google summer of code program, e.g. see here https://julialang.org/jsoc/ )
• given it is one of the $800+ proposals, which generally involve a lot of software engineering best practices, I require some form of "qualifying" step (see the bounty description page for more on that). A typical thing done in GSoC is to ask the applicant to submit and have merged a pull request on a smaller topic (that is not necessarily related), to prove they know how to use git and have some familiarity with the software. E.g. yesterday at the hackathon we discussed creating some symbolic simplification rules for Gaussian states in quantumsymbolics.jl.

[Krastanov]

also, we can make some of these "qualifying" pull requests small bounties of their own. I do find the qualifying step important, but I also do not want it to end up being just free labor

[akirakyle]

I'll be following this with great enthusiasm! Hopefully I'll have some time to help out since I've spent a lot of time (and frustration) implementing a lot of that Weedbrook et al paper in SymPy and have been really wanting to port it over to julia. A couple of pointers to those trying to implement this, of course all biased by the needs and issues I've had over the years:

• I highly suggest getting a copy of Serafini's textbook Quantum Continuous Variables as it contains a lot more details and derivations than Weedbrook's review paper, which can be immensely helpful when implementing non-trivial functionality.
• It look like your current GaussianOperator corresponds to a unitary channel, however it's quite useful to be able to specify arbitrary gaussian CPTP channels which are specified by a vector of displacements along with two 2N x 2N matrices which Weedbrook calls T and N (see section V.A.). Furthermore there are conditional Gaussian channels (measurements) which are non-trace preserving and which Weedbrook touches on in section II.E., however Serafini covers in much more generality.
• Finding the symplectic eigenvalues and more generally the Williamson decomposition of a covariance matrix gives a lot of insightful information about a gaussian state (e.g. purity, entanglement) so I would encourage you to keep this task in mind as you design the package. Doing this symbolically is unfortunately not possible in general for more than two modes due to the Abel–Ruffini theorem, however I've often found that given the structure of realistic covariance matrices, sympy can in fact find symbolic solutions for the symplectic eigenvalues for up to four modes and seems to be limited by how slow it is. Unfortunately I don't think there exists a pure julia symbolic polynomial factorization library that works with symbolics.jl, something I really wish would be implemented!

Hopefully this is helpful at this stage as you're thinking through the design of such a package.