Rotatingframe

Project.toml

@@ -4,6 +4,8 @@ authors = ["Joseph Broz <[email protected]>"]
 version = "0.5.0"
 
 [deps]
+CGcoefficient = "c862aa61-d51e-47d1-b396-b1e789b4e0b6"
+DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
@@ -15,6 +17,7 @@ QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
 QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
 YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
 

src/IonSim.jl

@@ -40,6 +40,7 @@ include("lasers.jl")
 include("iontraps.jl")
 include("chambers.jl")
 include("operators.jl")
+include("rotatingframes.jl")
 include("hamiltonians.jl")
 include("timeevolution.jl")
 include("species/_include_species.jl")

src/hamiltonians.jl

@@ -51,6 +51,7 @@ Constructs the Hamiltonian for `chamber` as a function of time. Return type is a
 """
 function hamiltonian(
     chamber::Chamber;
+    rotatingframe::Union{RotatingFrame,String}="interaction",
     timescale::Real=1,
     lamb_dicke_order::Union{Vector{Int}, Int}=1,
     rwa_cutoff::Real=Inf,
@@ -59,6 +60,7 @@ function hamiltonian(
 )
     return hamiltonian(
         chamber,
+        rotatingframe,
         iontrap(chamber),
         timescale,
         lamb_dicke_order,
@@ -74,8 +76,10 @@ end
 
 # At each time step, this function updates in-place the 2D array describing the full system
 # Hamiltonian.
+
 function hamiltonian(
     chamber::Chamber,
+    rf::Union{RotatingFrame,String},
     iontrap::LinearChain,
     timescale::Real,
     lamb_dicke_order::Union{Vector{Int}, Int},
@@ -85,14 +89,24 @@ function hamiltonian(
 )
     b, indxs, cindxs = _setup_base_hamiltonian(
         chamber,
+        rf,
         timescale,
         lamb_dicke_order,
         rwa_cutoff,
         displacement,
         time_dependent_eta
     )
     aui, gbi, gbs, bfunc, δνi, δνfuncs = _setup_fluctuation_hamiltonian(chamber, timescale)
+    interaction_picture = (typeof(rf)<:String)
+    if interaction_picture
+        @assert rf == "interaction" "argument rotatingframe must be the string \"interaction\" or be of type RotatingFrame"
+    else
+        Sdiag = rotating_eigenenergies(chamber,rf,timescale)
+    end
     S = SparseOperator(basis(chamber))
+
+    #could make rotating_eigenenergies return list and indices instead
+    #try to place this inside of an existing loop
     function f(t, ψ)  # a two argument function is required in the QuantumOptics solvers
         @inbounds begin
             @simd for i in 1:length(indxs)
@@ -116,6 +130,15 @@ function hamiltonian(
                     end
                 end
             end
+
+            #Hit a new loop with an @inbounds begin\end and run a single loop over the nonzero diagonal elements
+            #Have a list of nonzero diagonal elements and a list of the indices corresponding to them
+            #only need to do this loop if we're not in the interaction picture
+            if !interaction_picture
+                for j in 1:length(basis(chamber))
+                    S.data[j,j] = Sdiag.data[j,j]#rotating_eigenenergies(chamber,rf,timescale).data[j,j]#S.data[j,j] + 
+                end
+            end
             if length(gbi) == 0 && length(δνi) == 0
                 return S
             else
@@ -139,6 +162,7 @@ function hamiltonian(
                 end
             end
         end
+
         return S
     end
     return f
@@ -178,24 +202,35 @@ function does not keeps track of only one of these pairs.
 =#
 function _setup_base_hamiltonian(
     chamber,
+    rf,
     timescale,
     lamb_dicke_order,
     rwa_cutoff,
     displacement,
     time_dependent_eta
 )
     rwa_cutoff *= timescale
+
+    interaction_picture = (typeof(rf)<:String)
+
     allmodes = reverse(modes(chamber))
     L = length(allmodes)
-    νlist = Tuple([frequency(mode) for mode in allmodes])
+
+    if interaction_picture
+        νlist = Tuple([frequency(mode) for mode in allmodes])
+    else
+        vibrationalϕlist = [rotatingphase(rf, create(chamber, L-(i-1))) for i in 1:L]        
+        νlist = Tuple([vibrationalϕlist[i] for i in 1:L])
+    end
+
     mode_dims = [modecutoff(mode) + 1 for mode in allmodes]
 
     all_ions = reverse(ions(chamber))
     Q = prod([shape(ion)[1] for ion in all_ions])
     ion_arrays = [spdiagm(0 => [true for _ in 1:shape(ion)[1]]) for ion in all_ions]
 
     ηm, Δm, Ωm =
-        _ηmatrix(chamber), _Δmatrix(chamber, timescale), _Ωmatrix(chamber, timescale)
+        _ηmatrix(chamber), _Δmatrix(chamber, rf, timescale), _Ωmatrix(chamber, timescale)
     lamb_dicke_order = _check_lamb_dicke_order(lamb_dicke_order, L)
     ld_array, rows, vals = _ld_array(mode_dims, lamb_dicke_order, νlist, timescale)
     if displacement == "truncated" && time_dependent_eta
@@ -252,11 +287,13 @@ function _setup_base_hamiltonian(
                     ri = rows[i]
                     ri < j && continue
                     cf = vals[i]
+
                     pflag = abs(Δ_2π + cf) > rwa_cutoff
                     nflag = abs(Δ_2π - cf) > rwa_cutoff
+
                     (pflag && nflag) && continue
                     rev_indxs = false
-                    idxs = _inv_get_kron_indxs((rows[i], j), mode_dims)
+                    idxs = inv_get_kron_indxs((rows[i], j), mode_dims)
                     for l in 1:L
                         (idxs[1][l] ≠ idxs[2][l] && typeof(ηnm[l]) <: Number) && @goto cl
                     end
@@ -404,6 +441,34 @@ function _setup_base_hamiltonian(
     return functions, repeated_indices, conj_repeated_indices
 end
 
+function rotating_eigenenergies(chamber, rf, timescale)
+    S_diag = SparseOperator(basis(chamber))
+    allmodes = (modes(chamber))
+    allions = (ions(chamber))
+    mode_dims = [modecutoff(mode) + 1 for mode in allmodes]
+    ion_dims = [shape(ion)[1] for ion in allions]
+    dims = mode_dims
+    for dim in ion_dims
+        push!(dims,dim)
+    end
+    for j in 1:length(basis(chamber))
+        idxs = inv_get_kron_indxs((j,j), dims)[1]
+        Energy = 0
+
+        for (m,mode) in enumerate(allmodes)
+            Energy = Energy + frequency(mode)*(idxs[m]-1)
+        end
+
+        for (n,ion) in enumerate(allions)
+            Energy = Energy + energy(ion,sublevels(ion)[idxs[n+length(allmodes)]],B=chamber.B)
+        end
+        S_diag.data[j,j] = timescale*(Energy - unitary(rf)[j])
+    end
+    return S_diag
+end
+
+
+
 # δν(t) × aᵀa terms for Hamiltonian. This function returns an array of functions
 # δν_functions = [2π×ν.δν(t)×timescale for ν in modes]. It also returns an array of arrays
 # of arrays of indices, δν_indices, such that δν_indices[i][j] lists all diagonal elements
@@ -517,23 +582,30 @@ end
 # respectively. For each row/column we have a vector of detunings from the laser frequency
 # for each ion transition. We need to separate this calculation from _Ωmatrix to implement
 # RWA easily.
-function _Δmatrix(chamber, timescale)
+function _Δmatrix(chamber, rf, timescale)
     all_ions = ions(chamber)
     all_lasers = lasers(chamber)
     (N, M) = length(all_ions), length(all_lasers)
     B = bfield(chamber)
     ∇B = bgradient(chamber)
     Δnmkj = Array{Vector}(undef, N, M)
+    interaction_picture = (typeof(rf)<:String)
     for n in 1:N, m in 1:M
         Btot = bfield(chamber, all_ions[n])
         v = Vector{Float64}(undef, 0)
         for transition in subleveltransitions(all_ions[n])
-            ωa = transitionfrequency(all_ions[n], transition, B=Btot)
+            #RF Edit: Instead of adjusting by the atomic frequency, adjust according to the rotating frame.
+            #ωa = transitionfrequency(all_ions[n], transition, B=Btot)
+            if interaction_picture
+                ϕ = transitionfrequency(all_ions[n], transition, B=Btot)
+            else
+                ϕ = rotatingphase(rf, n, transition)
+            end
             push!(
                 v,
                 2π *
                 timescale *
-                ((c / wavelength(all_lasers[m])) + detuning(all_lasers[m]) - ωa)
+                ((c / wavelength(all_lasers[m])) + detuning(all_lasers[m]) - ϕ)
             )
         end
         Δnmkj[n, m] = v
@@ -632,7 +704,7 @@ end
 
 # The inverse of _get_kron_indxs. If T = X₁ ⊗ X₂ ⊗ X₃ and X₁, X₂, X₃ are M×M, N×N and L×L
 # dimension matrices, then we should input dims=(M, N, L).
-function _inv_get_kron_indxs(indxs, dims)
+function inv_get_kron_indxs(indxs, dims)
     row, col = indxs
     N = length(dims)
     ret_rows = Array{Int64}(undef, N)
@@ -703,4 +775,4 @@ function _ld_array(mode_dims, lamb_dicke_order, νlist, timescale)
     end
     length(a) == 1 ? ld_array = a[1] : ld_array = kron(a...)
     return ld_array, rowvals(ld_array), log.(nonzeros(ld_array))
-end
+end;

src/operators.jl

@@ -24,18 +24,91 @@ returns the creation operator for `v` such that: `create(v) * v[i] = √(i+1) *
 """
 create(v::VibrationalMode) = SparseOperator(v, diagm(-1 => sqrt.(1:(modecutoff(v)))))
 
+"""
+    create(c::Chamber, mode_idx::Int64; onlymodes=false)
+Returns the creation operator for the mode specified by mode_idx, tensored with the identity for all other modes.
+If onlymodes=false, then the output is also tensored with the identity for all ions.
+"""
+function create(c::Chamber, mode_idx::Int64; onlymodes=false)
+    @assert (mode_idx <= length(modes(c))) & (mode_idx > 0) "Invalid mode index."
+    v = modes(c)[mode_idx]
+    if mode_idx == 1
+        if length(modes(c)) > 1
+            adag = create(v) ⊗ tensor([one(modes(c)[j]) for j in 2:length(modes(c))]...)
+        else
+            adag = create(v)
+        end
+    else
+        adag = one(modes(c)[1])
+        for j in 2:length(modes(c))
+            if j == mode_idx
+                adag = adag ⊗ create(v)
+            else
+                adag = adag ⊗ one(modes(c)[j])
+            end
+        end
+    end
+    if onlymodes
+        return adag
+    end
+    return tensor([one(ion) for ion in ions(c)]...) ⊗ adag
+end
+
+
+
+
 """
     destroy(v::VibrationalMode)
 Returns the destruction operator for `v` such that: `destroy(v) * v[i] = √i * v[i-1]`.
 """
 destroy(v::VibrationalMode) = create(v)'
 
+"""
+    destroy(c::Chamber, mode_idx::Int64; onlymodes=false)
+Returns the destruction operator for the mode specified by mode_idx, tensored with the identity for all other modes.
+If onlymodes=false, then the output is also tensored with the identity for all ions.
+"""
+function destroy(c::Chamber, mode_idx::Int64; onlymodes=false)
+    return SparseOperator(create(c,mode_idx,onlymodes=onlymodes)')
+end
+
 """
     number(v::VibrationalMode)
 Returns the number operator for `v` such that:  `number(v) * v[i] = i * v[i]`.
 """
 number(v::VibrationalMode) = SparseOperator(v, diagm(0 => 0:(modecutoff(v))))
 
+"""
+    number(c::Chamber, mode_idx::Int64; onlymodes=false)
+Returns the number operator for the mode specified by mode_idx, tensored with the identity for all other modes.
+If onlymodes=false, then the output is also tensored with the identity for all ions.
+"""
+
+function number(c::Chamber, mode_idx::Int64; onlymodes=false)
+    @assert (mode_idx <= length(modes(c))) & (mode_idx > 0) "Invalid mode index."
+    v = modes(c)[mode_idx]
+    if mode_idx == 1
+        if length(modes(c))>1
+            num = number(v) ⊗ tensor([one(modes(c)[j]) for j in 2:length(modes(c))]...)
+        else
+            num = number(v)
+        end
+    else
+        num = one(modes(c)[1])
+        for j in 2:length(modes(c))
+            if j == mode_idx
+                num = num ⊗ number(v)
+            else
+                num = num ⊗ one(modes(c)[j])
+            end
+        end
+    end
+    if onlymodes
+        return num
+    end
+    return SparseOperator(tensor([one(ion) for ion in ions(c)]...) ⊗ num)
+end
+
 """
     displace(v::VibrationalMode, α::Number; method="truncated")
 Returns the displacement operator ``D(α)`` corresponding to `v`.
@@ -68,6 +141,33 @@ function displace(v::VibrationalMode, α::Number; method="truncated")
     end
 end
 
+"""
+    displace(c::Chamber, mode_idx::Int64, α::Number; method="truncated", onlymodes=false)
+Returns the displacement operator for the mode specified by mode_idx, D(α), tensored with the identity for all other modes.
+If onlymodes=false, then the output is also tensored with the identity for all ions.
+"""
+function displace(c::Chamber, mode_idx::Int64, α::Number; method="truncated", onlymodes=false)
+    @assert (mode_idx <= length(modes(c))) & (mode_idx > 0) "Invalid mode index."
+    v = modes(c)[mode_idx]
+    if mode_idx == 1
+        dis = displace(v,α,method=method) ⊗ tensor([one(modes(c)[j]) for j in 2:length(modes(c))]...)
+    else
+        dis = one(modes(c)[1])
+        for j in 2:length(modes(c))
+            if j == mode_idx
+                dis = dis ⊗ displace(v,α,method=method)
+            else
+                dis = dis ⊗ one(modes(c)[j])
+            end
+        end
+    end
+    if onlymodes
+        return dis
+    end
+    return SparseOperator(tensor([one(ion) for ion in ions(c)]...) ⊗ dis)
+end
+
+
 """
     thermalstate(v::VibrationalMode, n̄::Real; method="truncated")
 Returns a thermal density matrix with ``⟨a^†a⟩ ≈ n̄``. Note: approximate because we are
@@ -183,6 +283,32 @@ sigma(ion::Ion, ψ1::T, ψ2::T) where {T <: Union{Tuple{String, Real}, String, I
     sparse(projector(ion[ψ1], dagger(ion[ψ2])))
 sigma(ion::Ion, ψ1::Union{Tuple{String, Real}, String, Int}) = sigma(ion, ψ1, ψ1)
 
+"""
+    sigma(c::Chamber, ion_idx::Int64, ψ1::sublevel, ψ2::sublevel, onlyions=false)
+Returns the transition operator for the ion specified by ion_idx undergoing the transition ψ2 --> ψ1, tensored with the identity for all other ions.
+If onlyions=false, then the output is also tensored with the identity for all modes.
+"""
+function sigma(c::Chamber, ion_idx::Int64, ψ1::T, ψ2::T, onlyions=false) where {T <: Union{Tuple{String, Real}, String, Int}}
+    @assert (ion_idx <= length(ions(c))) & (ion_idx > 0) "Invalid ion index."
+    ion = ions(c)[ion_idx]
+    if ion_idx == 1
+        sig = sigma(ion,ψ1,ψ2) ⊗ tensor([one(ions(c)[j]) for j in 2:length(ions(c))]...)
+    else
+        sig = one(ions(c)[1])
+        for j in 2:length(ions(c))
+            if j == ion_idx
+                sig = sig ⊗ sigma(ion,ψ1,ψ2)
+            else
+                sig = sig ⊗ one(ions(c)[j])
+            end
+        end
+    end
+    if onlyions
+        return SparseOperator(sig)
+    end
+    return SparseOperator(sig ⊗ tensor([one(mode) for mode in modes(c)]...))
+end
+
 """
     ionprojector(obj, sublevels...; only_ions=false)