M1 transitions

src/constants.jl

@@ -2,7 +2,7 @@ module PhysicalConstants
 
 import Base.sqrt
 
-export μB, ħ, c, e, ϵ₀, α, kB, eye3, c_rank1, c_rank2
+export μB, ħ, c, e, ϵ₀, α, kB, gₛ, eye3, c_rank1, c_rank2
 
 """
     PhysicalConstant(x::Real)
@@ -31,6 +31,8 @@ const ϵ₀ = PhysicalConstant(8.85418782e-12, "(s^4A^2) / (m^3 kg)")
 const α = PhysicalConstant(0.007297352557920479, "")
 """`kB` = 1.38064852e-23 ``m^2kg/(s^2K)``"""
 const kB = PhysicalConstant(1.38064852e-23, "m^2kg/(s^2K)")
+"""`gₛ` = -2.00231930436256"""
+const gₛ = PhysicalConstant(-2.00231930436256, "")
 
 #############################################################################################
 # 3D real-space tensors

src/ions.jl

@@ -1,6 +1,6 @@
 using WignerSymbols: wigner3j, wigner6j
 using LinearAlgebra: cross
-#using .PhysicalConstants: e, ħ, α, μB, e, eye3, c_rank1, c_rank2
+#using .PhysicalConstants: e, ħ, α, μB, e, gₛ, eye3, c_rank1, c_rank2
 using IonSim.PhysicalConstants
 
 export Ion,
@@ -396,6 +396,7 @@ function transitionwavelength(I::Ion, transition::Tuple; B = 0, ignore_starkshif
            transitionfrequency(I, transition, B = B, ignore_starkshift = ignore_starkshift)
 end
 
+#TODO: Add M1s here!!
 """
     leveltransitions(I::Ion)
 Returns all allowed transitions between levels of `I` as a vector of `Tuple{String,String}`.
@@ -554,6 +555,28 @@ function matrix_element(
             units_factor = abs(e * Efield / (2ħ) * sqrt(15 * A12 / (α * c * k^3)))
             return units_factor * hyperfine_factor * geometric_factor / 2π
         end
+    elseif multipole == "M1"
+        # as implemented, this will error as it requires an API change to include l and s.
+        # I'd like to consider adding "term" as a data structure. Almost always, these arguments are
+        # being passed together, and I think it would be useful to have them in a structure that's
+        # easy to unpack, even if just to ease the syntax burden on the user
+        if abs(q) > 1
+            return 0
+        else
+            hyperfine_factor = abs(
+                sqrt((2 * f1 + 1) * (2 * f2 + 1) * (2 * j1 + 1) * (2 * j2 + 1)) * wigner6j(j2, I, f2, f1, 1, j1) *
+                (
+                wigner6j(l, j1, s, j2, l, 1) sqrt((2 * l + 1) * (l + 1) * l) +
+                gₛ * wigner6j(s, j1, l, j2, s, 1) sqrt((2 * s + 1) * (s + 1) * s)
+                ))
+            geometric_factor = abs(
+                sqrt(2j2 + 1) *
+                wigner3j(f2, 1, f1, -m2, q, m1) *
+                (transpose(c_rank1[q + 2, :]) * cross(khat_rotated, ϵhat_rotated))
+            )
+            units_factor = e * mu_B * Efield / c
+            return units_factor * hyperfine_factor * geometric_factor / 2π
+        end
     else
         @error "calculation of atomic transition matrix element for transition type $multipole not currently supported"
     end

src/lightfields.jl

@@ -0,0 +1,187 @@
+using .PhysicalConstants: c, ϵ₀
+
+export Laser, Microwave
+
+abstract type LightField end
+
+function Base.:(==)(LF1::LightField, LF2::LightField)
+    if typeof(LF1) != typeof(LF2)
+        return false
+    end
+    for field in fieldnames(typeof(LF1))
+        if getfield(LF1, field) != getfield(LF2, field)
+            return false
+        end
+    end
+    return true
+end
+
+function Base.print(LF::LightField)
+    println("λ: ", LF.λ, " m")
+    println("Δ: ", LF.Δ, " Hz")
+    println("ϵ̂: ", "(x=$(LF.ϵ.x), y=$(LF.ϵ.y), z=$(LF.ϵ.z))")
+    println("k̂: ", "(z=$(LF.k.x), y=$(LF.k.y), z=$(LF.k.z))")
+    println("E(t=0): ", "$(LF.E(0.0)) V/m")
+    return println("ϕ(t=0): ", "$(LF.ϕ(0.0)) ⋅ 2π")
+end
+
+
+"""
+    Laser(;λ=nothing, E=0, Δ=0, ϵ=(x̂+ŷ)/√2, k=ẑ, ϕ=0, pointing::Array{Tuple{Int,Real}})
+
+The physical parameters defining laser light.
+**args**
+* `λ::Union{Real,Nothing}`: the wavelength of the laser in meters
+* `E::Union{Function,Real}`: magnitude of the E-field in V/m
+* `Δ`: static detuning from f = c/λ in [Hz]
+* `ϵ::NamedTuple`: (ϵ.x, ϵ.y, ϵ.z), polarization direction, requires norm of 1
+* `k::NamedTuple`: (k.x, k.y, k.z), propagation direction, requires norm of 1
+* `ϕ::Union{Function,Real}`: time-dependent phase. of course, this can also be used to model a 
+    time-dependent detuning. Units are in radians. Note: if this is set to a function of time,
+    then when constructing a Hamiltonian with the `hamiltonian` function, the units of time
+    will be as specified by the `timescale` keyword argument.
+* `pointing`: an array of `Tuple{Int,Real}` for describing ion-laser pointing configuration.
+    (first element of the tuple is the index for an ion and the second element is the scaling
+    factor for the laser's Efield which must be between 0 and 1).
+"""
+mutable struct Laser <: LightField
+    λ::Union{Real, Nothing}
+    E::Function
+    Δ::Real
+    ϵ::NamedTuple{(:x, :y, :z)}
+    k::NamedTuple{(:x, :y, :z)}
+    ϕ::Function
+    pointing::Vector
+    function Laser(;
+        λ = nothing,
+        E::TE = 0,
+        Δ = 0,
+        ϵ = (x̂ + ŷ) / √2,
+        k = ẑ,
+        ϕ::Tϕ = 0,
+        pointing = Array{Tuple{Int, <:Real}}(undef, 0)
+    ) where {TE, Tϕ}
+        rtol = 1e-6
+        @assert isapprox(norm(ϵ), 1, rtol = rtol) "!(|ϵ| = 1)"
+        @assert isapprox(norm(k), 1, rtol = rtol) "!(|k| = 1)"
+        # @assert isapprox(ndot(ϵ, k), 0, rtol=rtol) "!(ϵ ⟂ k)"
+        # Above commented out until we figure out a better place to put this warning
+        a = pointing
+        (ion_num, scaling) = map(x -> getfield.(a, x), fieldnames(eltype(a)))
+        @assert length(ion_num) == length(unique(ion_num)) (
+            "a laser is pointing at the same ion twice"
+        )
+        for s in scaling
+            @assert 0 <= s <= 1 "must have s ∈ [0,1]"
+        end
+        TE <: Number ? Et(t) = E : Et = E
+        Tϕ <: Number ? ϕt(t) = ϕ : ϕt = ϕ
+        return new(λ, Et, Δ, ϵ, k, ϕt, pointing)
+    end
+    # for copying
+    Laser(λ, E, Δ, ϵ, k, ϕ, pointing) = new(λ, E, Δ, ϵ, k, ϕ, pointing)
+end
+
+function Base.setproperty!(L::Laser, s::Symbol, v::Tv) where {Tv}
+    rtol = 1e-6
+    if s == :ϵ
+        @assert isapprox(norm(v), 1, rtol = rtol) "!(|ϵ| = 1)"
+        # if ! isapprox(ndot(L.k, v), 0, rtol=rtol)
+        #     @warn "!(ϵ ⟂ k)"
+        # end 
+    elseif s == :k
+        @assert isapprox(norm(v), 1, rtol = rtol) "!(|k| = 1)"
+        # if ! isapprox(ndot(v, L.ϵ), 0, rtol=rtol)
+        #     @warn "!(ϵ ⟂ k)"
+        # end 
+    elseif s == :pointing
+        b = Tv <: Vector{Tuple{Int64, Float64}} || Tv <: Vector{Tuple{Int64, Int64}}
+        @assert b "type != Vector{Tuple{Int,Real}}"
+        (ion_num, scaling) = map(x -> getfield.(v, x), fieldnames(eltype(v)))
+        @assert length(ion_num) == length(unique(ion_num)) (
+            "a laser is pointing at the same ion twice"
+        )
+        for s in scaling
+            @assert 0 <= s <= 1 "must have s ∈ [0,1]"
+        end
+    elseif s == :E || s == :ϕ
+        Tv <: Number ? vt(t) = v : vt = v
+        Core.setproperty!(L, s, vt)
+        return
+    end
+    return Core.setproperty!(L, s, v)
+end
+
+# This could probably be made into one function with Unitfuls?
+laser_from_intensity(λ, i, Δ, ϵ, k, ϕ, pointing) = Laser(λ, sqrt(2i / (c*ϵ₀)), Δ, ϵ, k, ϕ, pointing)
+laser_from_waist_power(λ, P, Δ, ϵ, k, ϕ, pointing, w₀) = laser_from_intensity(λ, 2P/(π*w₀^2), Δ, ϵ, k, ϕ, pointing)
+
+"""
+    Microwave(;λ=nothing, E=0, Δ=0, ϵ=(x̂+ŷ)/√2, k=ẑ, ϕ=0)
+
+The physical parameters defining microwave radiation.
+**args**
+* `λ::Union{Real,Nothing}`: the wavelength of the microwave in meters
+* `E::Union{Function,Real}`: magnitude of the E-field in V/m
+* `Δ`: static detuning from f = c/λ in [Hz]
+* `ϵ::NamedTuple`: (ϵ.x, ϵ.y, ϵ.z), polarization direction, requires norm of 1
+* `k::NamedTuple`: (k.x, k.y, k.z), propagation direction, requires norm of 1
+* `ϕ::Union{Function,Real}`: time-dependent phase. of course, this can also be used to model a 
+    time-dependent detuning. Units are in radians. Note: if this is set to a function of time,
+    then when constructing a Hamiltonian with the `hamiltonian` function, the units of time
+    will be as specified by the `timescale` keyword argument.
+"""
+mutable struct Microwave <: LightField
+    λ::Union{Real, Nothing}
+    E::Function
+    Δ::Real
+    ϵ::NamedTuple{(:x, :y, :z)}
+    k::NamedTuple{(:x, :y, :z)}
+    ϕ::Function
+    pointing::Array{Tuple{Int, <:Real}}(undef, 0)
+    function Microwave(;
+        λ = nothing,
+        E::TE = 0,
+        Δ = 0,
+        ϵ = (x̂ + ŷ) / √2,
+        k = ẑ,
+        ϕ::Tϕ = 0,
+    ) where {TE, Tϕ}
+        rtol = 1e-6
+        @assert isapprox(norm(ϵ), 1, rtol = rtol) "!(|ϵ| = 1)"
+        @assert isapprox(norm(k), 1, rtol = rtol) "!(|k| = 1)"
+        # @assert isapprox(ndot(ϵ, k), 0, rtol=rtol) "!(ϵ ⟂ k)"
+        # Above commented out until we figure out a better place to put this warning
+        TE <: Number ? Et(t) = E : Et = E
+        Tϕ <: Number ? ϕt(t) = ϕ : ϕt = ϕ
+
+        return new(λ, Et, Δ, ϵ, k, ϕt, [])
+    end
+    # for copying
+    Microwave(λ, E, Δ, ϵ, k, ϕ) = new(λ, E, Δ, ϵ, k, ϕ)
+end
+
+
+function Base.setproperty!(MW::Microwave, s::Symbol, v::Tv) where {Tv}
+    rtol = 1e-6
+    if s == :ϵ
+        @assert isapprox(norm(v), 1, rtol = rtol) "!(|ϵ| = 1)"
+        # if ! isapprox(ndot(L.k, v), 0, rtol=rtol)
+        #     @warn "!(ϵ ⟂ k)"
+        # end 
+    elseif s == :k
+        @assert isapprox(norm(v), 1, rtol = rtol) "!(|k| = 1)"
+        # if ! isapprox(ndot(v, L.ϵ), 0, rtol=rtol)
+        #     @warn "!(ϵ ⟂ k)"
+        # end 
+    elseif s == :E || s == :ϕ
+        Tv <: Number ? vt(t) = v : vt = v
+        Core.setproperty!(MW, s, vt)
+        return
+    end
+    return Core.setproperty!(MW, s, v)
+end
+
+# This could probably be made into one function with Unitfuls?
+microwave_from_B(λ, B, Δ, ϵ, k, ϕ) = Microwave(λ, B*c, Δ, ϵ, k, ϕ)
+microwave_from_intensity(λ, i, Δ, ϵ, k, ϕ, pointing) = Microwave(λ, sqrt(2i / (c*ϵ₀)), Δ, ϵ, k, ϕ)