Further development of Fock states

docs/make.jl

@@ -27,6 +27,7 @@ function main()
         "Getting Started with QuantumSymbolics.jl" => "introduction.md",
         "Express Functionality" => "express.md",
         "Qubit Basis Choice" => "qubit_basis.md",
+        "Quantum Harmonic Oscillators" => "QHO.md",
         "API" => "API.md",
     ]
     )

docs/src/QHO.md

@@ -0,0 +1,159 @@
+# Quantum Harmonic Oscillators
+
+```@meta
+DocTestSetup = quote
+    using QuantumSymbolics, QuantumOptics
+end
+```
+
+In this section, we describe symbolic representations of bosonic systems in QuantumSymbolics, which can be numerically translated to [`QuantumOptics.jl`](https://github.com/qojulia/QuantumOptics.jl).
+
+## States
+
+A Fock state is a state with well defined number of excitation quanta of a single quantum harmonic oscillator (an eigenstate of the number operator). In the following example, we create a `FockState` with 3 quanta in an infinite-dimension Fock space:
+
+```jldoctest
+julia> f = FockState(3)
+|3⟩
+```
+
+Both vacuum (ground) and single-photon states are defined as constants in both unicode and ASCII for convenience:
+
+- `vac = F₀ = F0` $=|0\rangle$ in the number state representation,
+- `F₁ = F1` $=|1\rangle$ in the number state representation.
+
+To create quantum analogues of a classical harmonic oscillator, or monochromatic electromagnetic waves, we can define a coherent (a.k.a. semi-classical) state $|\alpha\rangle$, where $\alpha$ is a complex amplitude, with `CoherentState(α::Number)`:
+
+```jldoctest
+julia> c = CoherentState(im)
+|im⟩
+```
+!!! note "Naming convention for quantum harmonic oscillator bases"
+    The defined basis for arbitrary symbolic bosonic states is a `FockBasis` object, due to a shared naming interface for Quantum physics packages. For instance, the command `basis(CoherentState(im))` will output `Fock(cutoff=Inf)`. This may lead to confusion, as not all bosonic states are Fock states. However, this is simply a naming convention for the basis, and symbolic and numerical results are not affected by it.
+
+## Operators
+
+Operations on bosonic states are supported, and can be simplified with `qsimplify` and its rewriter `qsimplify_fock`. For instance, we can apply the raising (creation) $\hat{a}^{\dagger}$ and lowering (annihilation or destroy) $\hat{a}$ operators on a Fock state as follows:
+
+```jldoctest
+julia> f = FockState(3);
+
+julia> raise = Create*f
+a†|3⟩
+
+julia> qsimplify(raise, rewriter=qsimplify_fock)
+(sqrt(4))|4⟩
+
+julia> lower = Destroy*f
+a|3⟩
+
+julia> qsimplify(lower, rewriter=qsimplify_fock)
+(sqrt(3))|2⟩
+```
+Or, we can apply the number operator $\hat{n}$ to our Fock state:
+
+```jldoctest
+julia> f = FockState(3);
+
+julia> num = N*f
+n|3⟩
+
+julia> qsimplify(num, rewriter=qsimplify_fock)
+3|3⟩
+```
+
+Constants are defined for number and ladder operators in unicode and ASCII:
+
+- `N = n̂` $=\hat{n}$,
+- `Create = âꜛ` $=\hat{a}^{\dagger}$,
+- `Destroy = â` $=\hat{a}$.
+
+Phase-shift $U(\theta)$ and displacement $D(\alpha)$ operators, defined respectively as 
+$$U(\theta) = \exp\left(-i\theta\hat{n}\right) \quad \text{and} \quad D(\alpha) = \exp\left(\alpha\hat{a}^{\dagger} - \alpha\hat{a}\right),$$ 
+can be defined with usual simplification rules. Consider the following example:
+
+```jldoctest
+julia> displace = DisplaceOp(im)
+D(im)
+
+julia> c = qsimplify(displace*vac, rewriter=qsimplify_fock)
+|im⟩
+
+julia> phase = PhaseShiftOp(pi)
+U(π)
+
+julia> qsimplify(phase*c, rewriter=qsimplify_fock)
+|1.2246467991473532e-16 - 1.0im⟩
+```
+Here, we generated a coherent state $|i\rangle$ from the vacuum state $|0\rangle$ by applying the displacement operator defined by `DisplaceOp`. Then, we shifted its phase by $\pi$ with the phase shift operator (which is called with `PhaseShiftOp`) to get the result $|-i\rangle$.
+
+Summarized below are supported bosonic operators.
+
+- Number operator: `NumberOp()`,
+- Creation operator: `CreateOp()`,
+- Annihilation operator: `DestroyOp()`,
+- Phase-shift operator: `PhaseShiftOp(phase::Number)`,
+- Displacement operator: `DisplaceOp(alpha::Number)`.
+
+## Numerical Conversions to QuantumOptics.jl
+
+Bosonic systems can be translated to the ket representation with `express`. For instance:
+
+```jldoctest
+julia> f = FockState(1);
+
+julia> express(f)
+Ket(dim=3)
+  basis: Fock(cutoff=2)
+ 0.0 + 0.0im
+ 1.0 + 0.0im
+ 0.0 + 0.0im
+
+julia> express(Create) |> dense
+Operator(dim=3x3)
+  basis: Fock(cutoff=2)
+ 0.0+0.0im      0.0+0.0im  0.0+0.0im
+ 1.0+0.0im      0.0+0.0im  0.0+0.0im
+ 0.0+0.0im  1.41421+0.0im  0.0+0.0im
+
+julia> express(Create*f)
+Ket(dim=3)
+  basis: Fock(cutoff=2)
+                0.0 + 0.0im
+                0.0 + 0.0im
+ 1.4142135623730951 + 0.0im
+
+julia> express(Destroy*f)
+Ket(dim=3)
+  basis: Fock(cutoff=2)
+ 1.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+```
+
+!!! warning "Cutoff specifications for numerical representations of quantum harmonic oscillators"
+    Symbolic bosonic states and operators are naturally represented in an infinite dimension basis. For numerical conversions of such quantum objects, a finite cutoff of the highest allowed state must be defined. By default, the basis dimension of numerical conversions is set to 3 (so the number representation cutoff is 2), as demonstrated above. To define a different cutoff, one must customize the `QuantumOpticsRepr` instance, e.g. provide `QuantumOpticsRepr(cutoff=n::Int)` to `express`.
+
+If we wish to specify a different numerical cutoff, say 4, to the previous examples, then we rewrite them as follows:
+
+```jldoctest
+julia> f = FockState(1);
+
+julia> express(f, QuantumOpticsRepr(cutoff=4))
+Ket(dim=5)
+  basis: Fock(cutoff=4)
+ 0.0 + 0.0im
+ 1.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+
+julia> express(Create, QuantumOpticsRepr(4)) |> dense
+Operator(dim=5x5)
+  basis: Fock(cutoff=4)
+ 0.0+0.0im      0.0+0.0im      0.0+0.0im  0.0+0.0im  0.0+0.0im
+ 1.0+0.0im      0.0+0.0im      0.0+0.0im  0.0+0.0im  0.0+0.0im
+ 0.0+0.0im  1.41421+0.0im      0.0+0.0im  0.0+0.0im  0.0+0.0im
+ 0.0+0.0im      0.0+0.0im  1.73205+0.0im  0.0+0.0im  0.0+0.0im
+ 0.0+0.0im      0.0+0.0im      0.0+0.0im  2.0+0.0im  0.0+0.0im
+```

docs/src/express.md

@@ -6,15 +6,16 @@ DocTestSetup = quote
 end
 ```
 
-A principle feature of `QuantumSymbolics` is to numerically represent symbolic quantum expressions in various formalisms using [`express`](@ref). In particular, one can translate symbolic logic to back-end toolboxes such as `QuantumOptics.jl` or `QuantumClifford.jl` for simulating quantum systems with great flexibiity.
+A principle feature of `QuantumSymbolics` is to numerically represent symbolic quantum expressions in various formalisms using [`express`](@ref). In particular, one can translate symbolic logic to back-end toolboxes such as [`QuantumOptics.jl`](https://github.com/qojulia/QuantumOptics.jl) or [`QuantumClifford.jl`](https://github.com/QuantumSavory/QuantumClifford.jl) for simulating quantum systems with great flexibiity.
 
 As a straightforward example, consider the spin-up state $|\uparrow\rangle = |0\rangle$, the eigenstate of the Pauli operator $Z$, which can be expressed in `QuantumSymbolics` as follows:
 
 ```@example 1
 using QuantumSymbolics, QuantumClifford, QuantumOptics # hide
 ψ = Z1
 ```
-Using [`express`](@ref), we can translate this symbolic object into its numerical state vector form in `QuantumOptics.jl`.
+
+Using [`express`](@ref), we can translate this symbolic object into its numerical state vector form in [`QuantumOptics.jl`](https://github.com/qojulia/QuantumOptics.jl).
 
 ```@example 1
 express(ψ)
@@ -26,7 +27,7 @@ By default, [`express`](@ref) converts a quantum object with `QuantumOpticRepr`.
 ψ.metadata
 ```
 
-The caching feature of [`express`](@ref) prevents a specific representation for a symbolic quantum object from being computed more than once. This becomes handy for translations of more complex operations, which can become computationally expensive. We also have the ability to express $|Z_1\rangle$ in the Clifford formalism with `QuantumClifford.jl`:
+The caching feature of [`express`](@ref) prevents a specific representation for a symbolic quantum object from being computed more than once. This becomes handy for translations of more complex operations, which can become computationally expensive. We also have the ability to express $|Z_1\rangle$ in the Clifford formalism with [`QuantumClifford.jl`](https://github.com/QuantumSavory/QuantumClifford.jl):
 
 ```@example 1
 express(ψ, CliffordRepr())
@@ -56,4 +57,24 @@ julia> express(σʸ, CliffordRepr(), UseAsObservable())
 
 julia> express(σʸ, CliffordRepr(), UseAsOperation())
 sY
+```
+
+Another edge case is translations with `QuantumOpticsRepr`, where we can additionally define a finite cutoff for bosonic states and operators, as discussed in the [quantum harmonic oscillators page](@ref Quantum-Harmonic-Oscillators). The default cutoff for such objects is 2, however a different cutoff can be specified by passing an integer to `QuantumOpticsRepr` in an `express` call. Let us see an example with the number operator:
+
+```jldoctest
+julia> express(N) |> dense
+Operator(dim=3x3)
+  basis: Fock(cutoff=2)
+ 0.0+0.0im  0.0+0.0im  0.0+0.0im
+ 0.0+0.0im  1.0+0.0im  0.0+0.0im
+ 0.0+0.0im  0.0+0.0im  2.0+0.0im
+
+julia> express(N, QuantumOpticsRepr(cutoff=4)) |> dense
+Operator(dim=5x5)
+  basis: Fock(cutoff=4)
+ 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
+ 0.0+0.0im  1.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
+ 0.0+0.0im  0.0+0.0im  2.0+0.0im  0.0+0.0im  0.0+0.0im
+ 0.0+0.0im  0.0+0.0im  0.0+0.0im  3.0+0.0im  0.0+0.0im
+ 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  4.0+0.0im
 ```

ext/QuantumOpticsExt/QuantumOpticsExt.jl

@@ -1,13 +1,14 @@
 module QuantumOpticsExt
 
 using QuantumInterface, QuantumOpticsBase
+using QuantumInterface: samebases
 using QuantumSymbolics
 using QuantumSymbolics:
     HGate, XGate, YGate, ZGate, CPHASEGate, CNOTGate, PauliP, PauliM,
     XCXGate, XCYGate, XCZGate, YCXGate, YCYGate, YCZGate, ZCXGate, ZCYGate, ZCZGate,
     XBasisState, YBasisState, ZBasisState,
     NumberOp, CreateOp, DestroyOp,
-    FockBasisState,
+    FockState,
     MixedState, IdentityOp,
     qubit_basis, inf_fock_basis
 import QuantumSymbolics: express, express_nolookup
@@ -70,34 +71,14 @@ express_nolookup(s::XBasisState, ::QuantumOpticsRepr) = (_s₊,_s₋)[s.idx]
 express_nolookup(s::YBasisState, ::QuantumOpticsRepr) = (_i₊,_i₋)[s.idx]
 express_nolookup(s::ZBasisState, ::QuantumOpticsRepr) = (_l0,_l1)[s.idx]
 
-function express_nolookup(o::FockBasisState, r::QuantumOpticsRepr)
-    @warn "Fock space cutoff is not specified so we default to 2"
-    @assert o.idx<2 "without a specified cutoff you can not create states higher than 1 photon"
-    return (_f0₂,_f1₂)[o.idx+1]
-end
-function express_nolookup(o::NumberOp, r::QuantumOpticsRepr)
-    @warn "Fock space cutoff is not specified so we default to 2"
-    return _n₂
-end
-function express_nolookup(o::CreateOp, r::QuantumOpticsRepr)
-    @warn "Fock space cutoff is not specified so we default to 2"
-    return _ad₂
-end
-function express_nolookup(o::DestroyOp, r::QuantumOpticsRepr)
-    @warn "Fock space cutoff is not specified so we default to 2"
-    return _a₂
-end
-
+express_nolookup(s::FockState, r::QuantumOpticsRepr) = fockstate(FockBasis(r.fock_cutoff),s.idx)
+express_nolookup(s::CoherentState, r::QuantumOpticsRepr) = coherentstate(FockBasis(r.fock_cutoff),s.alpha)
+express_nolookup(o::NumberOp, r::QuantumOpticsRepr) = number(FockBasis(r.fock_cutoff))
+express_nolookup(o::CreateOp, r::QuantumOpticsRepr) = create(FockBasis(r.fock_cutoff))
+express_nolookup(o::DestroyOp, r::QuantumOpticsRepr) = destroy(FockBasis(r.fock_cutoff))
+express_nolookup(o::DisplaceOp, r::QuantumOpticsRepr) = displace(FockBasis(r.fock_cutoff), o.alpha)
 express_nolookup(x::MixedState, ::QuantumOpticsRepr) = identityoperator(basis(x))/length(basis(x)) # TODO there is probably a more efficient way to represent it
-function express_nolookup(x::IdentityOp, ::QuantumOpticsRepr)
-    b = basis(x)
-    if b!=inf_fock_basis
-        return identityoperator(basis(x)) # TODO there is probably a more efficient way to represent it
-    else
-        @warn "Fock space cutoff is not specified so we default to 2"
-        return identityoperator(_bf2)
-    end
-end
+express_nolookup(x::IdentityOp, r::QuantumOpticsRepr) = identityoperator(FockBasis(r.fock_cutoff))
 
 express_nolookup(p::PauliNoiseCPTP, ::QuantumOpticsRepr) = LazySuperSum(SpinBasis(1//2), [1-p.px-p.py-p.pz,p.px,p.py,p.pz],
                                                                [LazyPrePost(_id,_id),LazyPrePost(_x,_x),LazyPrePost(_y,_y),LazyPrePost(_z,_z)])

src/QSymbolicsBase/QSymbolicsBase.jl

@@ -1,5 +1,5 @@
 using Symbolics
-import Symbolics: simplify
+import Symbolics: simplify,Term
 using SymbolicUtils
 import SymbolicUtils: Symbolic,_isone,flatten_term,isnotflat,Chain,Fixpoint,Prewalk,sorted_arguments
 using TermInterface
@@ -27,7 +27,7 @@ export SymQObj,QObj,
        I,X,Y,Z,σˣ,σʸ,σᶻ,Pm,Pp,σ₋,σ₊,
        H,CNOT,CPHASE,XCX,XCY,XCZ,YCX,YCY,YCZ,ZCX,ZCY,ZCZ,
        X1,X2,Y1,Y2,Z1,Z2,X₁,X₂,Y₁,Y₂,Z₁,Z₂,L0,L1,Lp,Lm,Lpi,Lmi,L₀,L₁,L₊,L₋,L₊ᵢ,L₋ᵢ,
-       vac,F₀,F0,F₁,F1,
+       vac,F₀,F0,F₁,F1,inf_fock_basis,
        N,n̂,Create,âꜛ,Destroy,â,basis,SpinBasis,FockBasis,
        SBra,SKet,SOperator,SHermitianOperator,SUnitaryOperator,SHermitianUnitaryOperator,SSuperOperator,
        @ket,@bra,@op,@superop,
@@ -39,10 +39,10 @@ export SymQObj,QObj,
        MixedState,IdentityOp,
        SApplyKet,SApplyBra,SMulOperator,SSuperOpApply,SCommutator,SAnticommutator,SBraKet,SOuterKetBra,
        HGate,XGate,YGate,ZGate,CPHASEGate,CNOTGate,
-       XBasisState,YBasisState,ZBasisState,
-       NumberOp,CreateOp,DestroyOp,
+       XBasisState,YBasisState,ZBasisState,FockState,CoherentState,
+       NumberOp,CreateOp,DestroyOp,PhaseShiftOp,DisplaceOp,
        XCXGate,XCYGate,XCZGate,YCXGate,YCYGate,YCZGate,ZCXGate,ZCYGate,ZCZGate,
-       qsimplify,qsimplify_pauli,qsimplify_commutator,qsimplify_anticommutator,
+       qsimplify,qsimplify_pauli,qsimplify_commutator,qsimplify_anticommutator,qsimplify_fock,
        qexpand,
        isunitary,
        KrausRepr,kraus
@@ -144,7 +144,6 @@ end
 Base.isequal(::SymQObj, ::Symbolic{Complex}) = false
 Base.isequal(::Symbolic{Complex}, ::SymQObj) = false
 
-
 # TODO check that this does not cause incredibly bad runtime performance
 # use a macro to provide specializations if that is indeed the case
 propsequal(x,y) = all(n->(n==:metadata || isequal(getproperty(x,n),getproperty(y,n))), propertynames(x))
@@ -167,6 +166,7 @@ include("basic_superops.jl")
 include("linalg.jl")
 include("predefined.jl")
 include("predefined_CPTP.jl")
+include("predefined_fock.jl")
 
 ##
 # Symbolic and simplification rules

src/QSymbolicsBase/express.jl

@@ -78,7 +78,17 @@ end
 ##
 
 """Representation using kets, bras, density matrices, and superoperators governed by `QuantumOptics.jl`."""
-struct QuantumOpticsRepr <: AbstractRepresentation end
+struct QuantumOpticsRepr <: AbstractRepresentation 
+    fock_cutoff::Union{Nothing, Int}
+    function QuantumOpticsRepr(fock_cutoff::Union{Nothing, Int})
+        if isnothing(fock_cutoff)
+            return new(2)
+        else
+            return new(fock_cutoff)
+        end
+    end
+end
+QuantumOpticsRepr() = QuantumOpticsRepr(nothing)
 """Similar to `QuantumOpticsRepr`, but using trajectories instead of superoperators."""
 struct QuantumMCRepr <: AbstractRepresentation end
 """Representation using tableaux governed by `QuantumClifford.jl`"""