Further development of Fock states (#69)

* fock rules

* Fock ops

* docs and tests

* reorganize predef file and docstrings

* express and tests

* update metadata

* add fock.md to make.jl

* review changes

* rm space

* some minor stylistic changes and and also breaking your code with a non-working example

* update changelog and kwdef QuantumOpticsRepr

---------

Co-authored-by: Stefan Krastanov <[email protected]>

CHANGELOG.md

@@ -1,8 +1,12 @@
 # News
 
-## v0.3.5 - dev
+## v0.4.0 - 2024-08-03
 
 - Cleaned up metadata decoration of struct definitions.
+- Added documentation for quantum harmonic oscillators.
+- Added phase-shift and displacement operators `DisplaceOp` and `PhaseShiftOp`.
+- Simplification rules for Fock objects.
+- **(breaking)** `FockBasisState` was renamed to `FockState`.
 
 ## v0.3.4 - 2024-07-22
 

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumSymbolics"
 uuid = "efa7fd63-0460-4890-beb7-be1bbdfbaeae"
 authors = ["QuantumSymbolics.jl contributors"]
-version = "0.3.5-dev"
+version = "0.4.0"
 
 [deps]
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"

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,15 +1,16 @@
 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
+    qubit_basis
 import QuantumSymbolics: express, express_nolookup
 using TermInterface
 using TermInterface: isexpr, head, operation, arguments, metadata
@@ -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.cutoff),s.idx)
+express_nolookup(s::CoherentState, r::QuantumOpticsRepr) = coherentstate(FockBasis(r.cutoff),s.alpha)
+express_nolookup(o::NumberOp, r::QuantumOpticsRepr) = number(FockBasis(r.cutoff))
+express_nolookup(o::CreateOp, r::QuantumOpticsRepr) = create(FockBasis(r.cutoff))
+express_nolookup(o::DestroyOp, r::QuantumOpticsRepr) = destroy(FockBasis(r.cutoff))
+express_nolookup(o::DisplaceOp, r::QuantumOpticsRepr) = displace(FockBasis(r.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.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
@@ -28,7 +28,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,
@@ -40,10 +40,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
@@ -119,7 +119,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))
@@ -142,6 +141,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,9 @@ end
 ##
 
 """Representation using kets, bras, density matrices, and superoperators governed by `QuantumOptics.jl`."""
-struct QuantumOpticsRepr <: AbstractRepresentation end
+@kwdef struct QuantumOpticsRepr <: AbstractRepresentation 
+    cutoff::Int = 2
+end
 """Similar to `QuantumOpticsRepr`, but using trajectories instead of superoperators."""
 struct QuantumMCRepr <: AbstractRepresentation end
 """Representation using tableaux governed by `QuantumClifford.jl`"""