Add symplectic representation feature

CHANGELOG.md

@@ -1,8 +1,9 @@
 # News
 
-## v1.1.1 - dev
+## v1.2.0 - dev
 
 - Add a benchmark suite as a part of the Github workflows.
+- **(breaking)** Require a `SymplecticBasis` subtype to be defined when creating a Gaussian object.
 
 ## v1.1.0 - 2024-11-18
 

benchmark/benchmarks.jl

@@ -16,8 +16,9 @@ for op in op_list
 end
 
 for nmodes in (2, 10, 50, 100, 200)
-    SUITE["operations"]["unitary product"][string(nmodes)] = @benchmarkable op * state setup=(op = randunitary($nmodes); state = randstate($nmodes))
-    SUITE["operations"]["channel product"][string(nmodes)] = @benchmarkable ch * state setup=(ch = randchannel($nmodes); state = randstate($nmodes))
-    SUITE["operations"]["tensor product"][string(nmodes)] = @benchmarkable state1 ⊗ state2 setup=(state1 = randstate($nmodes); state2 = randstate($nmodes))
-    SUITE["operations"]["partial trace"][string(nmodes)] = @benchmarkable ptrace(state, indices) setup=(state = randstate($nmodes); indices = collect(2:$nmodes))
+    basis = QuadPairBasis(nmodes)
+    SUITE["operations"]["unitary product"][string(nmodes)] = @benchmarkable op * state setup=(op = randunitary($basis); state = randstate($basis))
+    SUITE["operations"]["channel product"][string(nmodes)] = @benchmarkable ch * state setup=(ch = randchannel($basis); state = randstate($basis))
+    SUITE["operations"]["tensor product"][string(nmodes)] = @benchmarkable state1 ⊗ state2 setup=(state1 = randstate($basis); state2 = randstate($basis))
+    SUITE["operations"]["partial trace"][string(nmodes)] = @benchmarkable ptrace(state, indices) setup=(state = randstate($basis); indices = collect(2:$nmodes))
 end

docs/src/manual.md

@@ -38,11 +38,13 @@ such state vector simulations, we recommend using the [QuantumOptics.jl](https:/
 which can be typed in the Julia REPL as `\otimes<TAB>`. Take the following example, where we produce a 3-mode Gaussian state that consists of a coherent state, vacuumstate, and squeezed state:
 
 ```jldoctest
-julia> coherentstate(-1.0) ⊗ vacuumstate() ⊗ squeezedstate(0.25, pi/4)
-GaussianState for 3 modes.
+julia> basis = QuadPairBasis(1);
+
+julia> coherentstate(basis, -1.0+im) ⊗ vacuumstate(basis) ⊗ squeezedstate(basis, 0.25, pi/4)
+GaussianState for 3 modes in QuadPairBasis representation.
 mean: 6-element Vector{Float64}:
  -1.4142135623730951
-  0.0
+  1.4142135623730951
   0.0
   0.0
   0.0
@@ -56,6 +58,30 @@ covariance: 6×6 Matrix{Float64}:
  0.0  0.0  0.0  0.0  0.184235  0.748048
 ```
 
+Note that in the above example, we defined the symplectic basis of the Gaussian state to be [`QuadPairBasis`](@ref), which determines the arrangement of our quadrature field operators to be pairwise: $\mathbf{\hat{x}} = (q_1, p_1, q_2, p_2, q_3, p_3)^{\text{T}}$. If we wanted the field operators to be ordered blockwise, i.e.,
+$\mathbf{\hat{x}} = (q_1, q_2, q_3, p_1, p_2, p_3)^{\text{T}}$ then we would call [`QuadBlockBasis`](@ref) instead:
+
+```julia
+julia> basis = QuadBlockBasis(1);
+
+julia> coherentstate(basis, -1.0+im) ⊗ vacuumstate(basis) ⊗ squeezedstate(basis, 0.25, pi/4)
+GaussianState for 3 modes in QuadBlockBasis representation.
+mean: 6-element Vector{Float64}:
+ -1.4142135623730951
+  0.0
+  0.0
+  1.4142135623730951
+  0.0
+  0.0
+covariance: 6×6 Matrix{Float64}:
+ 1.0  0.0  0.0       0.0  0.0  0.0
+ 0.0  1.0  0.0       0.0  0.0  0.0
+ 0.0  0.0  0.379578  0.0  0.0  0.184235
+ 0.0  0.0  0.0       1.0  0.0  0.0
+ 0.0  0.0  0.0       0.0  1.0  0.0
+ 0.0  0.0  0.184235  0.0  0.0  0.748048
+```
+
 ## Gaussian Operators
 
 To transform Gaussian states into Gaussian states, we need Gaussian maps. There are various ways to construct Gaussian transformations, which we will discuss in this section.
@@ -74,16 +100,11 @@ an `N`-mode Gaussian bosonic system. As long as we have a displacement vector of
 !!! note
     A matrix $\mathbf{S}$ of size $2N\times 2N$ is symplectic when it satisfies the following relation:
 
-    $$\mathbf{S} \mathbf{\Omega} \mathbf{S}^{\text{T}} = \mathbf{\Omega}, \qquad \mathbf{\Omega} = \bigoplus_{i=1}^{N} \begin{pmatrix} 0 & -1 \\ -1 & 0 \end{pmatrix}.$$
+    $$\mathbf{S} \mathbf{\Omega} \mathbf{S}^{\text{T}} = \mathbf{\Omega}.
 
     In this library, we define symplectic matrices with respect to $\Omega$, the *symplectic form*, which satisfies the canonical
     commutation relation $[\hat{x}_i, \hat{x}_j] = 2i\Omega_{ij}$, where $\hat{x}_i$ and $\hat{x}_j$
-    are quadrature field operators. However, a different basis can of course be used for a symplectic matrix. 
-    Another common way of expressing symplectic matrices is with the block diagonal symplectic form:
-
-    $$\mathbf{J} = \begin{pmatrix} 0 & \mathbf{I}_N \\ -\mathbf{I}_N & 0 \end{pmatrix},$$
-
-    where $\mathbf{I}_N$ is the identity matrix of size $N \times N$.
+    are quadrature field operators.
 
 This library has a number of predefined Gaussian unitaries, which are listed below:
 
@@ -113,10 +134,12 @@ Listed below are a list of predefined Gaussian channels supported by Gabs:
     method can be called with an additional noise matrix to create a [`GaussianChannel`](@ref) object. For instance, a noisy displacement operator can be called with [`displace`](@ref) as follows:
 
     ```jldoctest
+    julia> basis = QuadPairBasis(1);
+
     julia> noise = [1.0 -2.0; 4.0 -3.0];
 
-    julia> displace(1.0-im, noise)
-    GaussianChannel for 1 mode.
+    julia> displace(basis, 1.0-im, noise)
+    GaussianChannel for 1 mode in QuadPairBasis representation.
     displacement: 2-element Vector{Float64}:
       1.4142135623730951
      -1.4142135623730951

docs/src/tutorials.md

@@ -16,14 +16,16 @@ currently has support for the following distributions, which can be called with
 Below is a code example that plots the wigner distribution of a vacuum state:
 ```@example
 using Gabs, CairoMakie
-state = vacuumstate()
+basis = QuadPairBasis(1)
+state = vacuumstate(basis)
 q, p = collect(-5.0:0.5:5.0), collect(-5.0:0.5:5.0) # phase space coordinates
 heatmap(q, p, state, dist = :wigner)
 ```
 Of course, more plotting sugar can be added to this example with internal Makie attributes:
 ```@example
 using Gabs, CairoMakie
-state = vacuumstate()
+basis = QuadPairBasis(1)
+state = vacuumstate(basis)
 q, p = collect(-5.0:0.5:5.0), collect(-5.0:0.5:5.0) # phase space coordinates
 fig = Figure(fontsize=15, size = (375, 300), fonts = (; regular="CMU Serif"))
 ax = Axis(fig[1,1], xlabel = L"q", ylabel = L"p")
@@ -41,8 +43,8 @@ follow the [`AbstractArray` interface](https://docs.julialang.org/en/v1/manual/i
 Let's explore this feature in detail, using [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) and [SparseArrays.jl](https://github.com/JuliaSparse/SparseArrays.jl) as a case study. To create a coherent state that wraps around pure Julia arrays, the function
 [`coherentstate`](@ref) can be called with a single complex argument:
 ```jldoctest
-julia> coherentstate(1.0-im)
-GaussianState for 1 mode.
+julia> coherentstate(QuadPairBasis(1), 1.0-im)
+GaussianState for 1 mode in QuadPairBasis representation.
 mean: 2-element Vector{Float64}:
   1.4142135623730951
  -1.4142135623730951
@@ -57,8 +59,8 @@ can specify an array type in its first (and second) arguments. Let's see an exam
 ```jldoctest
 julia> using StaticArrays
 
-julia> state = coherentstate(SVector{2}, SMatrix{2,2}, 1.0-im)
-GaussianState for 1 mode.
+julia> state = coherentstate(SVector{2}, SMatrix{2,2}, QuadPairBasis(1), 1.0-im)
+GaussianState for 1 mode in QuadPairBasis representation.
 mean: 2-element SVector{2, Float64} with indices SOneTo(2):
   1.4142135623730951
  -1.4142135623730951
@@ -67,7 +69,7 @@ covariance: 2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
  0.0  1.0
 
 julia> tp = state ⊗ state
-GaussianState for 2 modes.
+GaussianState for 2 modes in QuadPairBasis representation.
 mean: 4-element SVector{4, Float64} with indices SOneTo(4):
   1.4142135623730951
  -1.4142135623730951
@@ -82,7 +84,7 @@ covariance: 4×4 SMatrix{4, 4, Float64, 16} with indices SOneTo(4)×SOneTo(4):
 julia> using SparseArrays
 
 julia> ptrace(SparseVector, SparseMatrixCSC, tp, 1)
-GaussianState for 1 mode.
+GaussianState for 1 mode in QuadPairBasis representation.
 mean: 2-element SparseVector{Float64, Int64} with 2 stored entries:
   [1]  =  1.41421
   [2]  =  -1.41421
@@ -97,8 +99,8 @@ Importantly, methods that create or manipulate a Gaussian state, such as [`tenso
     and matrices, then you can specify the array type with a single argument. For instance, to initialize a state that contains `Array`s holding numbers of type `Float32` rather
     than `Float64`, simply pass `Array{Float32}` to any relevant Gabs.jl method:
     ```jldoctest
-    julia> state = displace(Array{Float32}, 1.0-im)
-    GaussianUnitary for 1 mode.
+    julia> state = displace(Array{Float32}, QuadPairBasis(1), 1.0-im)
+    GaussianUnitary for 1 mode in QuadPairBasis representation.
     displacement: 2-element Vector{Float32}:
       1.4142135
      -1.4142135

src/Gabs.jl

@@ -10,6 +10,8 @@ import QuantumInterface: StateVector, AbstractOperator, apply!, tensor, ⊗
 export 
     # types
     GaussianState, GaussianUnitary, GaussianChannel, Generaldyne,
+    # symplectic representations
+    QuadPairBasis, QuadBlockBasis,
     # operations
     tensor, ⊗, apply!, ptrace, output, prob,
     # predefined Gaussian states
@@ -30,6 +32,8 @@ include("errors.jl")
 
 include("utils.jl")
 
+include("symplectic.jl")
+
 include("types.jl")
 
 include("states.jl")