Further development of Fock states

docs/make.jl

@@ -26,6 +26,7 @@ function main()
     pages = [
         "QuantumSymbolics.jl" => "index.md",
         "Qubit Basis Choice" => "qubit_basis.md",
+        "Fock States" => "fock.md",
         "Express Functionality" => "express.md",
         "API" => "API.md",
     ]

docs/src/fock.md

@@ -0,0 +1,147 @@
+# Fock States
+
+```@meta
+DocTestSetup = quote
+    using QuantumSymbolics, QuantumOptics
+end
+```
+
+In this section, we describe symbolic representations of Fock states in QuantumSymbolics, which can be numerically translated to `QuantumOptics.jl`. A Fock state is the number state representation $|n\rangle$ of a system with $n$ particles. 
+
+One can define a basis of the Fock space with `FockBasis(N, offset=0)`, where `N` is the highest available Fock state and `offset` is the lowest cutoff state, which is by default zero. In the following example, we create a `FockState` with 3 quanta in an infinite-dimension Fock space:
+
+```jldoctest
+julia> b = FockBasis(Inf, 0.0)
+Fock(cutoff=Inf)
+
+julia> f = FockState(3, b)
+|3⟩
+```
+
+Both vacuum (ground) and single-photon states in an infinite-dimension Fock basis 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 state $|\alpha\rangle$, where $\alpha$ is a complex amplitude, with `CoherentState(α::Number, basis::Basis)`:
+
+```jldoctest
+julia> b = FockBasis(Inf, 0.0);
+
+julia> c = CoherentState(im, b)
+|im⟩
+```
+
+## Operators
+
+Operations on fock states are supported, and can be simplified with `qsimplify` and its rewriter `qsimplify_fock`. For instance, we can apply the raising (creation) $a^{\dagger}$ and lowering (annihilation or destroy) $a$ operators on a fock state as follows:
+
+```jldoctest
+julia> b = FockBasis(Inf, 0.0);
+
+julia> f = FockState(3, b);
+
+julia> raise = Create*f
+a†|3⟩
+
+julia> qsimplify(raise, rewriter=qsimplify_fock)
+2.0|4⟩
+
+julia> lower = Destroy*f
+a|3⟩
+
+julia> qsimplify(lower, rewriter=qsimplify_fock)
+1.7320508075688772|2⟩
+```
+Or, we can apply the number operator $\hat{n}$ to our fock state:
+
+```jldoctest
+julia> b = FockBasis(Inf, 0.0);
+
+julia> f = FockState(3, b);
+
+julia> num = N*f
+n|3⟩
+
+julia> qsimplify(num, rewriter=qsimplify_fock)
+3|3⟩
+```
+
+In the infinite dimension case for Fock states, constants are defined for number and ladder operators:
+
+- `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. Consider the following example:
+
+```jldoctest
+julia> b = FockBasis(Inf, 0.0);
+
+julia> displace = DisplaceOp(im, b)
+D(im)
+
+julia> c = qsimplify(displace*vac, rewriter=qsimplify_fock)
+|im⟩
+
+julia> phase = PhaseShiftOp(pi, b)
+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 operators, which can be defined in any Fock basis.
+
+- Number operator: `NumberOp(basis:Basis)`,
+- Creation operator: `CreateOp(basis::Basis)`,
+- Annihilation operator: `DestroyOp(basis::Basis)`,
+- Phase-shift operator: `PhaseShiftOp(phase::Number, basis:Basis)`,
+- Displacement operator: `DisplaceOp(alpha::Number, basis::Basis)`.
+
+## Numerical Conversions to QuantumOptics.jl
+
+Fock systems can be translated to the ket representation with `express`. For instance:
+
+```jldoctest
+julia> using QuantumOptics
+
+julia> b = FockBasis(3);
+
+julia> f = FockState(1, b);
+
+julia> express(f)
+Ket(dim=4)
+  basis: Fock(cutoff=3)
+ 0.0 + 0.0im
+ 1.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+
+julia> express(CreateOp(b)) |> dense
+Operator(dim=4x4)
+  basis: Fock(cutoff=3)
+ 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  1.41421+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
+
+julia> express(CreateOp(b)*f)
+Ket(dim=4)
+  basis: Fock(cutoff=3)
+                0.0 + 0.0im
+                0.0 + 0.0im
+ 1.4142135623730951 + 0.0im
+                0.0 + 0.0im
+
+julia> express(DestroyOp(b)*f)
+Ket(dim=4)
+  basis: Fock(cutoff=3)
+ 1.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.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,28 +71,67 @@
 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]
+function express_nolookup(s::FockState, r::QuantumOpticsRepr)
+    b = basis(s)
+    i = s.idx
+    if !samebases(b, inf_fock_basis)
+        fockstate(b, i)
+    else
+        @warn "Fock space cutoff is not specified so we default to 2"
+        return fockstate(_bf2, i)
+    end
+end
+function express_nolookup(s::CoherentState, r::QuantumOpticsRepr)
+    b = basis(s)
+    alpha = s.alpha
+    if !samebases(b, inf_fock_basis)
+        coherentstate(b, alpha)
+    else
+        @warn "Fock space cutoff is not specified so we default to 2"
+        return coherentstate(_bf2, alpha)
+    end
 end
 function express_nolookup(o::NumberOp, r::QuantumOpticsRepr)
-    @warn "Fock space cutoff is not specified so we default to 2"
-    return _n₂
+    b = basis(o)
+    if !samebases(b, inf_fock_basis)
+        number(b)
+    else
+        @warn "Fock space cutoff is not specified so we default to 2"
+        return number(_bf2)
+    end
 end
 function express_nolookup(o::CreateOp, r::QuantumOpticsRepr)
-    @warn "Fock space cutoff is not specified so we default to 2"
-    return _ad₂
+    b = basis(o)
+    if !samebases(b, inf_fock_basis)
+        create(b)
+    else
+        @warn "Fock space cutoff is not specified so we default to 2"
+        return create(_bf2)
+    end
 end
 function express_nolookup(o::DestroyOp, r::QuantumOpticsRepr)
-    @warn "Fock space cutoff is not specified so we default to 2"
-    return _a₂
+    b = basis(o)
+    if !samebases(b, inf_fock_basis)
+        destroy(b)
+    else
+        @warn "Fock space cutoff is not specified so we default to 2"
+        return destroy(_bf2)
+    end
+end
+function express_nolookup(o::DisplaceOp, r::QuantumOpticsRepr)
+    b = basis(o)
+    alpha = o.alpha
+    if !samebases(b, inf_fock_basis)
+        displace(b, alpha)
+    else
+        @warn "Fock space cutoff is not specified so we default to 2"
+        return displace(_bf2, alpha)
+    end
 end
-
 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
+    if !samebases(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"

src/QSymbolicsBase/QSymbolicsBase.jl

@@ -38,10 +38,10 @@ export SymQObj,QObj,
        SProjector,MixedState,IdentityOp,SInvOperator,
        SApplyKet,SApplyBra,SMulOperator,SSuperOpApply,SCommutator,SAnticommutator,SDagger,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
 
@@ -142,7 +142,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))
@@ -164,6 +163,7 @@ include("basic_ops_inhomogeneous.jl")
 include("linalg.jl")
 include("predefined.jl")
 include("predefined_CPTP.jl")
+include("predefined_fock.jl")
 
 ##
 # Symbolic and simplification rules

src/QSymbolicsBase/predefined.jl

@@ -25,24 +25,6 @@ symbollabel(x::YBasisState) = "Y$(num_to_sub(x.idx))"
 end
 symbollabel(x::ZBasisState) = "Z$(num_to_sub(x.idx))"
 
-@withmetadata struct FockBasisState <: SpecialKet
-    idx::Int
-    basis::Basis
-end
-symbollabel(x::FockBasisState) = "$(x.idx)"
-
-@withmetadata struct DiscreteCoherentState <: SpecialKet
-    alpha::Number # TODO parameterize
-    basis::Basis
-end
-symbollabel(x::DiscreteCoherentState) = "$(x.alpha)"
-
-@withmetadata struct ContinuousCoherentState <: SpecialKet
-    alpha::Number # TODO parameterize
-    basis::Basis
-end
-symbollabel(x::ContinuousCoherentState) = "$(x.alpha)"
-
 @withmetadata struct MomentumEigenState <: SpecialKet
     p::Number # TODO parameterize
     basis::Basis
@@ -69,13 +51,6 @@ const Z1 = const Z₁ = const L0 = const L₀ = ZBasisState(1, qubit_basis)
 """Basis state of σᶻ"""
 const Z2 = const Z₂ = const L1 = const L₁ = ZBasisState(2, qubit_basis)
 
-const inf_fock_basis = FockBasis(Inf,0.)
-"""Vacuum basis state of n"""
-const vac = const F₀ = const F0 = FockBasisState(0,inf_fock_basis)
-"""Single photon basis state of n"""
-const F₁ = const F1 = FockBasisState(1,inf_fock_basis)
-
-
 ##
 # Gates and Operators on qubits
 ##
@@ -173,29 +148,6 @@ const CNOT = CNOTGate()
 """CPHASE gate"""
 const CPHASE = CPHASEGate()
 
-##
-# Gates and Operators on harmonic oscillators
-##
-
-abstract type AbstractSingleBosonOp <: Symbolic{AbstractOperator} end
-abstract type AbstractSingleBosonGate <: AbstractSingleBosonOp end # TODO maybe an IsUnitaryTrait is a better choice
-isexpr(::AbstractSingleBosonGate) = false
-basis(::AbstractSingleBosonGate) = inf_fock_basis
-
-@withmetadata struct NumberOp <: AbstractSingleBosonOp end
-symbollabel(::NumberOp) = "n"
-@withmetadata struct CreateOp <: AbstractSingleBosonOp end
-symbollabel(::CreateOp) = "a†"
-@withmetadata struct DestroyOp <: AbstractSingleBosonOp end
-symbollabel(::DestroyOp) = "a"
-
-"""Number operator, also available as the constant `n̂`"""
-const N = const n̂ = NumberOp()
-"""Creation operator, also available as the constant `âꜛ` - there is no unicode dagger superscript, so we use the uparrow"""
-const Create = const âꜛ = CreateOp()
-"""Annihilation operator, also available as the constant `â`"""
-const Destroy = const â = DestroyOp()
-
 ##
 # Other special or useful objects
 ##