diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..056043b --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,70 @@ +# from: https://discourse.julialang.org/t/easy-workflow-file-for-setting-up-github-actions-ci-for-your-julia-package/49765 +name: CI +on: + pull_request: + branches: + - main + push: + branches: + - main + tags: '*' +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.7' # The minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'. + - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia. + - 'nightly' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: actions/cache@v1 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v1 + with: + file: lcov.info + docs: + name: Documentation + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: '1' + - run: | + julia --project=docs -e ' + using Pkg + Pkg.develop(PackageSpec(path=pwd())) + Pkg.instantiate()' + - run: | + julia --project=docs -e ' + using Documenter: doctest, DocMeta + using MoleculeSpectrum + DocMeta.setdocmeta!(MoleculeSpectrum, :DocTestSetup, :(using MoleculeSpectrum;); recursive=true) + doctest(MoleculeSpectrum)' + - run: julia --project=docs docs/make.jl + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..95f1421 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2022 JILA KRb Team + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md index 2c4dc8d..8f52f07 100644 --- a/README.md +++ b/README.md @@ -1,14 +1,20 @@ # MoleculeSpectrum.jl +[![CI Actions Status](https://github.com/cal-miller-harvard/MoleculeSpectrum/actions/workflows/ci.yml/badge.svg)](https://github.com/cal-miller-harvard/MoleculeSpectrum/actions) + Compute the energy levels of diatomic ${}^{1} \Sigma$ molecules in magnetic, electric, and optical fields. +Documentation is available [here](https://cal-miller-harvard.github.io/MoleculeSpectrum). + # Todo - [x] Add transition strengths - [ ] Dump to file, and/or dataframes? - [ ] Graphing - [x] Improve code organization -- [ ] Documentation - - [ ] Make it build! +- [x] Documentation + - [x] Make it build! + - [ ] Finish documenting all functions + - [ ] Document units of all numeric quantities - [ ] Remove unnecessary types in function arguments (?) - [ ] API ergonomics - [x] ExternalFields diff --git a/docs/Project.toml b/docs/Project.toml index 6c894be..fa45f32 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,3 @@ [deps] -HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" MoleculeSpectrum = "35893526-f1dd-48d2-8180-97f6707ed43b" -WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" diff --git a/docs/make.jl b/docs/make.jl index 5456c13..7a0dbfc 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,5 +1,18 @@ -push!(LOAD_PATH,"../src/") - -using Documenter, MoleculeSpectrum - -makedocs(sitename="MoleculeSpectrum documentation") \ No newline at end of file +push!(LOAD_PATH,"../src/") + +using Documenter, MoleculeSpectrum + +makedocs( + sitename="MoleculeSpectrum Documentation", + format = Documenter.HTML( + prettyurls = get(ENV, "CI", nothing) == "true" + ), + pages = [ + "MoleculeSpectrum" => "index.md", + "API Documentation" => "api.md" + ] +) + +deploydocs( + repo = "github.com/cal-miller-harvard/MoleculeSpectrum.jl.git", #TODO: Change to your username +) \ No newline at end of file diff --git a/docs/src/api.md b/docs/src/api.md new file mode 100644 index 0000000..e98bf73 --- /dev/null +++ b/docs/src/api.md @@ -0,0 +1,109 @@ +# MoleculeSpectrum.jl API Documentation + +```@meta + CurrentModule = MoleculeSpectrum + DocTestSetup = quote + using MoleculeSpectrum + end +``` + +## Contents +```@contents + Pages = ["api.md"] + Depth = 4 +``` + +## Molecular Parameters +### Molecular Parameters +```@docs +MolecularParameters +KRb_Parameters_Neyenhuis +KRb_Parameters_Ospelkaus +DEFAULT_MOLECULAR_PARAMETERS +TOY_MOLECULE_PARAMETERS +``` +### Polarizability Parameters +```@docs +Polarizability +KRb_Polarizability +``` +### Zeeman Parameters +```@docs +ZeemanParameters +KRb_Zeeman +``` +### Nuclear Parameters +```@docs +NuclearParameters +KRb_Nuclear_Neyenhuis +KRb_Nuclear_Ospelkaus +``` + +## Vectors and Tensors +### Spherical vectors +```@docs +SphericalVector +VectorX +VectorY +VectorZ +``` +### Unit vectors +```@docs +SphericalUnitVector +UnitVectorX +UnitVectorY +UnitVectorZ +``` +### Tensors +```@docs +T⁽¹⁾ +T⁽²⁾ +get_tensor_component +tensor_dot +``` + +## States +### Types +```@docs +State +KRbState +``` +### Functions +```@docs +index_to_state +state_to_index +order_by_overlap_with +max_overlap_with +find_closest_basis_state +decompose_to_basis_states +``` + +## Fields +```@docs +ExternalFields +DEFAULT_FIELDS +TEST_FIELDS +``` + +## Hamiltonians +```@docs +hamiltonian +HamiltonianParts +make_hamiltonian_parts +make_krb_hamiltonian_parts +``` + +## Spectra and Analysis +```@docs +Spectrum +calculate_spectrum +get_energy +get_energy_difference +find_transition_strengths +plot_transition_strengths +calculate_dipolar_interaction +``` + +## Index +```@index +``` \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 96f72db..7c5d0e1 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,15 +1,21 @@ -```@meta -CurrentModule = MoleculeSpectrum -``` - -```@docs -State(N, mₙ, I₁, mᵢ₁, I₂, mᵢ₂) -``` - -```@docs - plot_transition_strengths(spectrum::Spectrum, g::State, frequency_range; polarization::SphericalUnitVector=Unpolarized()) -``` - -```@docs - transition_strengths(spectrum::Spectrum, g::State, frequency_range; polarization::SphericalUnitVector=Unpolarized()) -``` \ No newline at end of file +# MoleculeSpectrum.jl + +Computes the energy levels of diatomic ${}^{1} \Sigma$ molecules in magnetic, electric, and optical fields. + +# Todo +- [x] Add transition strengths +- [ ] Dump to file, and/or dataframes? +- [ ] Graphing +- [x] Improve code organization +- [ ] Documentation + - [ ] Make it build! +- [ ] Remove unnecessary types in function arguments (?) +- [ ] API ergonomics + - [x] ExternalFields + - [ ] Get dipole matrix elements + - [ ] More generally... inspect any matrix element? + - [ ] Make more `calculate` functions polymorphic over `State`, eigenstate index, etc + - [ ] ... +- [x] Switch to `get_eigenstates` function +- [x] Explicitly type fields of all structs +- [ ] Many more tests!!! \ No newline at end of file diff --git a/src/MoleculeSpectrum.jl b/src/MoleculeSpectrum.jl index 68cd9d8..4d9ee07 100644 --- a/src/MoleculeSpectrum.jl +++ b/src/MoleculeSpectrum.jl @@ -1,277 +1,314 @@ -module MoleculeSpectrum - -import WignerSymbols: wigner3j -import HalfIntegers: HalfInt -import Gadfly -using LinearAlgebra, SparseArrays, StaticArrays, Test - -export ZeemanParameters, NuclearParameters, Polarizability, MolecularParameters -export KRb_Zeeman, KRb_Nuclear_Neyenhuis, KRb_Nuclear_Ospelkaus, KRb_Polarizability -export KRb_Parameters_Neyenhuis, - KRb_Parameters_Ospelkaus, DEFAULT_MOLECULAR_PARAMETERS, TOY_MOLECULE_PARAMETERS - -export SphericalVector, VectorX, VectorY, VectorZ -export SphericalUnitVector, UnitVectorX, UnitVectorY, UnitVectorZ, Unpolarized -export T⁽¹⁾, T⁽²⁾, get_tensor_component, tensor_dot -export ExternalFields, DEFAULT_FIELDS, TEST_FIELDS - -export State, KRbState, index_to_state, state_to_index -export order_by_overlap_with, max_overlap_with, find_closest_basis_state, decompose_to_basis_states -export get_energy, get_energy_difference - -export HamiltonianParts, make_hamiltonian_parts, hamiltonian, make_krb_hamiltonian_parts - -export Spectrum, calculate_spectrum -export find_transition_strengths, plot_transition_strengths -export calculate_dipolar_interaction - -module Constants -"Nuclear magneton in MHz/G\n" -const μN = 7.622593285e-4 -"Factor to convert from Debye to C m\n" -const DToSI = 3.33564e-30 -"Planck's constant (SI)\n" -const h = 6.62607004e-34 -"Convert from D*(V/cm) to MHz\n" -const DVcm⁻¹ToMHz = (DToSI / h) * 1e-4 -end # module - -include("molecular_parameters.jl") -include("fields.jl") - -""" - State - -Represents a molecular state in the uncoupled basis. -""" -struct State - N::Int - mₙ::Int - I::SVector{2,HalfInt} # [K, Rb] - mᵢ::SVector{2,HalfInt} -end - -""" - State(N, mₙ, I₁, mᵢ₁, I₂, mᵢ₂) - -Creates a basis state ``|N, mₙ, I₁, mᵢ₁, I₂, mᵢ₂⟩``. - -""" -State(N, mₙ, I₁, mᵢ₁, I₂, mᵢ₂) = State(N, mₙ, SVector(HalfInt(I₁), HalfInt(I₂)), SVector(HalfInt(mᵢ₁), HalfInt(mᵢ₂))) - -""" - KRbState(N, mₙ, mK, mRb) - -Creates a basis state ``|N, m_n, m_{\\text{K}}, m_{\\text{Rb}}⟩`` for ``{}^{40}\\text{K}^{87}\\text{Rb}``. - -This is a wrapper around [`State`](@ref) to avoid having to specify the nuclear spins ``I_k`` each time. - -See also [`State`](@ref). -""" -KRbState(N, mₙ, mK, mRb) = - State(N, mₙ, KRb_Parameters_Neyenhuis.I, [HalfInt(mK) HalfInt(mRb)]) - -include("matrix_elements.jl") -include("hamiltonian.jl") - -struct Spectrum - hamiltonian_parts::HamiltonianParts - energies::Vector{Float64} - eigenstates::Vector{Vector{ComplexF64}} -end - -get_eigenstate(spectrum, k) = spectrum.eigenstates[k] -get_energies(spectrum) = spectrum.energies -get_energies(spectrum, range) = filter(x -> (x >= range[1] && x <= range[2]), spectrum.energies) - -include("utility.jl") - -""" - calculate_spectrum(hamiltonian_parts, external_fields) - -Compute the energies and eigenstates under the external fields. - -To avoid reconstructing the Hamiltonian each time, `hamiltonian_parts` can be reused over calls -to `calculate_spectrum`. The output [`Spectrum`](@ref) object is used as an input for -further analysis, for example in [`transition_strengths`](@ref). - -See also [`make_hamiltonian_parts`](@ref), [`make_krb_hamiltonian_parts`](@ref), -[`Spectrum`](@ref). -""" -function calculate_spectrum( - hamiltonian_parts::HamiltonianParts, - external_fields::ExternalFields, -)::Spectrum - h = hamiltonian(hamiltonian_parts, external_fields) - energies = eigvals(h) - eigenstates = [c for c in eachcol(eigvecs(h))] - return Spectrum(hamiltonian_parts, energies, eigenstates) -end - -""" - find_transition_strengths(spectrum::Spectrum, g::State, frequency_range; polarization::Union{Int, SphericalUnitVector, Nothing}=nothing) - -Compute electric dipole transitions out of `g` with energy between `frequency_range[1]` and `frequency_range[2]`. - -The output is a `Vector` of tuples `(frequency, strength, closest_basis_state, eigenstate_index)`, produced in -decreasing order of transition strength. The `strength` is the absolute value of the dipole matrix element, -normalized by ``D/\\sqrt{3}`` (the maximum transition dipole between ``N = 0`` and ``N = 1``). We use -the absolute value of the matrix element, rather than its square, so the results are proportional to -Rabi frequency ``Ω``. - -The `polarization` keyword argument can be used to choose a specific microwave polarization. This defaults to `nothing`, -which returns the incoherent sum over ``σ-``, ``π``, and ``σ+``. If `polarization` is an `Int`, then it is interpreted as -the spherical component `p = -1:1` of the dipole operator (`p == -1` corresponds to ``σ-`` polarization). If -`polarization` is a `SphericalUnitVector`, then the polarization is interpreted as linear along that axis. - -There is a convenience method [`plot_transition_strengths`](@ref) that immediately produces a plot from the result. - -See also [`plot_transition_strengths`](@ref), [`make_hamiltonian_parts`](@ref), [`make_krb_hamiltonian_parts`](@ref), -[`Spectrum`](@ref). -""" -function find_transition_strengths( - spectrum::Spectrum, - g::State, - frequency_range; - polarization::Union{Int, SphericalUnitVector, Nothing} = nothing -) - energies = get_energies(spectrum) - - (overlap, index_g) = max_overlap_with(spectrum, g) - if overlap < 0.5 - @warn "The best overlap with your requested ground state is < 0.5." - end - E_g = energies[index_g] - g_state = get_eigenstate(spectrum, index_g) - - state_range = - searchsortedfirst(energies, frequency_range[1] + E_g):searchsortedlast(energies, frequency_range[2] + E_g) - states = [get_eigenstate(spectrum, k) for k in state_range] - frequencies = [energies[k] - E_g for k in state_range] - - if polarization === nothing - strengths = [calculate_transition_strength_incoherent(spectrum, g_state, e) for e in states] - else - strengths = [abs(calculate_transition_strength_coherent(spectrum, g_state, e, polarization)) for e in states] - end - closest_basis_states = map(e -> find_closest_basis_state(spectrum, e), state_range) - - out = [x for x in zip(frequencies, strengths, closest_basis_states, state_range)] - return sort!(out, by = t -> t[2], rev = true) -end - - -""" - calculate_transition_strength_incoherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}) - -Compute the transition strength from `g` to `e` for an even incoherent mixture of polarizations. - -The incoherent sum is formed by calculating the squared matrix elements of ``|⟨g|H_p|e⟩|^2`` for each -spherical component `p = -1:1` of the dipole Hamiltonian, summing the three values, and then taking the square root. -""" -function calculate_transition_strength_incoherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}) - h_dipole = spectrum.hamiltonian_parts.dipole_relative - intensities = [abs2(e' * h_dipole[p] * g) for p in 1:3] - - # Normalize to d/sqrt(3), which is the largest transition dipole (between |0,0> and |1,0>) - strength = sqrt(reduce(+, intensities)) / (1 / sqrt(3)) - return strength -end - -""" - calculate_transition_strength_coherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}, polarization::Int) - calculate_transition_strength_coherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}, polarization::SphericalUnitVector) - -Compute the transition strength from `g` to `e`, driven by a field with the given `polarization`. - -If `polarization` is an `Int`, then it is interpreted as the spherical component `p = -1:1` of the dipole operator (`p == -1` corresponds to -``σ-`` polarization). If `polarization` is a `SphericalUnitVector`, then the polarization is interpreted as linear along that axis. - -""" -function calculate_transition_strength_coherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}, polarization::Int)::ComplexF64 - @assert polarization <= 1 && polarization >= -1 - - index = polarization + 2 # components are p = -1, 0, 1 - h_dipole = spectrum.hamiltonian_parts.dipole_relative[index] - return e' * h_dipole * g / (1 / sqrt(3)) -end - -function calculate_transition_strength_coherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}, polarization::SphericalUnitVector)::ComplexF64 - h_dipole = tensor_dot(T⁽¹⁾(polarization), spectrum.hamiltonian_parts.dipole_relative) - return e' * h_dipole * g / (1 / sqrt(3)) -end - -""" - plot_transition_strengths(spectrum::Spectrum, g::State, frequency_range; polarization::Union{Int, SphericalUnitVector, Nothing}=nothing) - -Plot the frequencies and strengths of electric dipole transitions out of `g`, -with energy between `frequency_range[1]` and `frequency_range[2]`. - -The `polarization` keyword argument can be used to choose a specific microwave polarization. This defaults to `nothing`, -which returns the incoherent sum over ``σ-``, ``π``, and ``σ+``. If `polarization` is an `Int`, then it is interpreted as -the spherical component `p = -1:1` of the dipole operator (`p == -1` corresponds to ``σ-`` polarization). If -`polarization` is a `SphericalUnitVector`, then the polarization is interpreted as linear along that axis. - -This method calls [`find_transition_strengths`](@ref) internally. - -See also [`find_transition_strengths`](@ref), [`make_hamiltonian_parts`](@ref), [`make_krb_hamiltonian_parts`](@ref), -[`Spectrum`](@ref). -""" -function plot_transition_strengths( - spectrum::Spectrum, - g::State, - frequency_range; - polarization::Union{Int, SphericalUnitVector, Nothing}=nothing -) - transitions = - find_transition_strengths(spectrum, g, frequency_range; polarization = polarization) - - freqs = [t[1] for t in transitions] - strengths = [abs(t[2]) for t in transitions] - - Gadfly.plot( - x = freqs, - y = strengths, - Gadfly.Geom.hair, - Gadfly.Geom.point, - Gadfly.Guide.xlabel("Frequency (MHz)"), - Gadfly.Guide.ylabel("Transition strength"), - ) -end - -function calculate_dipolar_interaction( - spectrum::Spectrum, - g::State, - e::State; - p::Int = 0 -) - (overlap, index_g) = max_overlap_with(spectrum, g) - if overlap < 0.5 - @warn "The best overlap with your requested ground state is < 0.5." - end - - (overlap, index_e) = max_overlap_with(spectrum, e) - if overlap < 0.5 - @warn "The best overlap with your requested excited state is < 0.5." - end - - return calculate_dipolar_interaction( - spectrum, - get_eigenstate(spectrum, index_g), - get_eigenstate(spectrum, index_e); - p=p - ) -end - -function calculate_dipolar_interaction( - spectrum::Spectrum, - g::Vector{ComplexF64}, - e::Vector{ComplexF64}; - p::Int = 0 -) - d_1 = [Complex(calculate_transition_strength_coherent(spectrum, g, e, pol)) for pol=-1:1] - d_2 = [Complex(calculate_transition_strength_coherent(spectrum, e, g, pol)) for pol=-1:1] - return get_tensor_component(p, T⁽²⁾(d_1, d_2)) * sqrt(6) / 2 -end - -end # module +module MoleculeSpectrum + +import WignerSymbols: wigner3j +import HalfIntegers: HalfInt +import Gadfly +using LinearAlgebra, SparseArrays, StaticArrays, Test + +export ZeemanParameters, NuclearParameters, Polarizability, MolecularParameters +export KRb_Zeeman, KRb_Nuclear_Neyenhuis, KRb_Nuclear_Ospelkaus, KRb_Polarizability +export KRb_Parameters_Neyenhuis, + KRb_Parameters_Ospelkaus, DEFAULT_MOLECULAR_PARAMETERS, TOY_MOLECULE_PARAMETERS + +export SphericalVector, VectorX, VectorY, VectorZ +export SphericalUnitVector, UnitVectorX, UnitVectorY, UnitVectorZ +export T⁽¹⁾, T⁽²⁾, get_tensor_component, tensor_dot +export ExternalFields, DEFAULT_FIELDS, TEST_FIELDS + +export State, KRbState, index_to_state, state_to_index +export order_by_overlap_with, max_overlap_with, find_closest_basis_state, decompose_to_basis_states +export get_energy, get_energy_difference + +export HamiltonianParts, make_hamiltonian_parts, hamiltonian, make_krb_hamiltonian_parts + +export Spectrum, calculate_spectrum +export find_transition_strengths, plot_transition_strengths +export calculate_dipolar_interaction + +module Constants +"Nuclear magneton in MHz/G\n" +const μN = 7.622593285e-4 +"Factor to convert from Debye to C m\n" +const DToSI = 3.33564e-30 +"Planck's constant (SI)\n" +const h = 6.62607004e-34 +"Convert from D*(V/cm) to MHz\n" +const DVcm⁻¹ToMHz = (DToSI / h) * 1e-4 +end # module + +include("molecular_parameters.jl") +include("fields.jl") + +""" +``` +struct State + N::Int + mₙ::Int + I::SVector{2,HalfInt} + mᵢ::SVector{2,HalfInt} +end +``` + +Represents a molecular state in the uncoupled basis. +""" +struct State + N::Int + mₙ::Int + I::SVector{2,HalfInt} # [K, Rb] + mᵢ::SVector{2,HalfInt} +end + +""" + State(N, mₙ, I₁, mᵢ₁, I₂, mᵢ₂) + +Creates a basis state ``|N, mₙ, I₁, mᵢ₁, I₂, mᵢ₂⟩``. + +""" +State(N, mₙ, I₁, mᵢ₁, I₂, mᵢ₂) = State(N, mₙ, SVector(HalfInt(I₁), HalfInt(I₂)), SVector(HalfInt(mᵢ₁), HalfInt(mᵢ₂))) + +""" + KRbState(N, mₙ, mK, mRb) + +Creates a basis state ``|N, m_n, m_{\\text{K}}, m_{\\text{Rb}}⟩`` for ``{}^{40}\\text{K}^{87}\\text{Rb}``. + +This is a wrapper around [`State`](@ref) to avoid having to specify the nuclear spins ``I_k`` each time. + +See also [`State`](@ref). +""" +KRbState(N, mₙ, mK, mRb) = + State(N, mₙ, KRb_Parameters_Neyenhuis.I, [HalfInt(mK) HalfInt(mRb)]) + +include("matrix_elements.jl") +include("hamiltonian.jl") + +""" + Spectrum(hamiltonian_parts, energies, eigenstates) + +Represents a molecular spectrum, including eigenvalues and eigenstates, +given [`HamiltonianParts`](@ref) and [`ExternalFields`](@ref). + +Created from [`calculate_spectrum`](@ref). See [`decompose_to_basis_states`](@ref) and +[`get_energy`](@ref)for access to eigenstates and energies. +""" +struct Spectrum + hamiltonian_parts::HamiltonianParts + energies::Vector{Float64} + eigenstates::Vector{Vector{ComplexF64}} +end + +""" + get_eigenstate(spectrum, k) + +Returns the vector of coefficients corresponding to the `k`th lowest energy state. +""" +get_eigenstate(spectrum, k) = spectrum.eigenstates[k] + +""" + get_energies(spectrum) + get_energies(spectrum, range) + +Returns a list of the energies, in increasing order. + +If the optional argument `range` is included, only return energies between `range[1]` and `range[2]`. +""" +get_energies(spectrum) = spectrum.energies +get_energies(spectrum, range) = filter(x -> (x >= range[1] && x <= range[2]), spectrum.energies) + +include("utility.jl") + +""" + calculate_spectrum(hamiltonian_parts, external_fields) + +Compute the energies and eigenstates under the external fields. + +To avoid reconstructing the Hamiltonian each time, `hamiltonian_parts` can be reused over calls +to `calculate_spectrum`. The output [`Spectrum`](@ref) object is used as an input for +further analysis, for example in [`find_transition_strengths`](@ref). + +See also [`make_hamiltonian_parts`](@ref), [`make_krb_hamiltonian_parts`](@ref), +[`Spectrum`](@ref). +""" +function calculate_spectrum( + hamiltonian_parts::HamiltonianParts, + external_fields::ExternalFields, +)::Spectrum + h = hamiltonian(hamiltonian_parts, external_fields) + e = eigen(h) + energies = e.values + eigenstates = [c for c in eachcol(e.vectors)] + return Spectrum(hamiltonian_parts, energies, eigenstates) +end + +""" + find_transition_strengths(spectrum::Spectrum, g::State, frequency_range; polarization::Union{Int, SphericalUnitVector, Nothing}=nothing) + +Compute electric dipole transitions out of `g` with energy between `frequency_range[1]` and `frequency_range[2]`. + +The output is a `Vector` of tuples `(frequency, strength, closest_basis_state, eigenstate_index)`, produced in +decreasing order of transition strength. The `strength` is the absolute value of the dipole matrix element, +normalized by ``D/\\sqrt{3}`` (the maximum transition dipole between ``N = 0`` and ``N = 1``). We use +the absolute value of the matrix element, rather than its square, so the results are proportional to +Rabi frequency ``Ω``. + +The `polarization` keyword argument can be used to choose a specific microwave polarization. This defaults to `nothing`, +which returns the incoherent sum over ``σ-``, ``π``, and ``σ+``. If `polarization` is an `Int`, then it is interpreted as +the spherical component `p = -1:1` of the dipole operator (`p == -1` corresponds to ``σ-`` polarization). If +`polarization` is a `SphericalUnitVector`, then the polarization is interpreted as linear along that axis. + +There is a convenience method [`plot_transition_strengths`](@ref) that immediately produces a plot from the result. + +See also [`plot_transition_strengths`](@ref), [`make_hamiltonian_parts`](@ref), [`make_krb_hamiltonian_parts`](@ref), +[`Spectrum`](@ref). +""" +function find_transition_strengths( + spectrum::Spectrum, + g::State, + frequency_range; + polarization::Union{Int, SphericalUnitVector, Nothing} = nothing +) + energies = get_energies(spectrum) + + (overlap, index_g) = max_overlap_with(spectrum, g) + if overlap < 0.5 + @warn "The best overlap with your requested ground state is < 0.5." + end + E_g = energies[index_g] + g_state = get_eigenstate(spectrum, index_g) + + state_range = + searchsortedfirst(energies, frequency_range[1] + E_g):searchsortedlast(energies, frequency_range[2] + E_g) + states = [get_eigenstate(spectrum, k) for k in state_range] + frequencies = [energies[k] - E_g for k in state_range] + + if polarization === nothing + strengths = [calculate_transition_strength_incoherent(spectrum, g_state, e) for e in states] + else + strengths = [abs(calculate_transition_strength_coherent(spectrum, g_state, e, polarization)) for e in states] + end + closest_basis_states = map(e -> find_closest_basis_state(spectrum, e), state_range) + + out = [x for x in zip(frequencies, strengths, closest_basis_states, state_range)] + return sort!(out, by = t -> t[2], rev = true) +end + + +""" + calculate_transition_strength_incoherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}) + +Compute the transition strength from `g` to `e` for an even incoherent mixture of polarizations. + +The incoherent sum is formed by calculating the squared matrix elements of ``|⟨g|H_p|e⟩|^2`` for each +spherical component `p = -1:1` of the dipole Hamiltonian, summing the three values, and then taking the square root. +""" +function calculate_transition_strength_incoherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}) + h_dipole = spectrum.hamiltonian_parts.dipole_relative + intensities = [abs2(e' * h_dipole[p] * g) for p in 1:3] + + # Normalize to d/sqrt(3), which is the largest transition dipole (between |0,0> and |1,0>) + strength = sqrt(reduce(+, intensities)) / (1 / sqrt(3)) + return strength +end + +""" + calculate_transition_strength_coherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}, polarization::Int) + calculate_transition_strength_coherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}, polarization::SphericalUnitVector) + +Compute the transition strength from `g` to `e`, driven by a field with the given `polarization`. + +If `polarization` is an `Int`, then it is interpreted as the spherical component `p = -1:1` of the dipole operator (`p == -1` corresponds to +``σ-`` polarization). If `polarization` is a `SphericalUnitVector`, then the polarization is interpreted as linear along that axis. + +""" +function calculate_transition_strength_coherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}, polarization::Int)::ComplexF64 + @assert polarization <= 1 && polarization >= -1 + + index = polarization + 2 # components are p = -1, 0, 1 + h_dipole = spectrum.hamiltonian_parts.dipole_relative[index] + return e' * h_dipole * g / (1 / sqrt(3)) +end + +function calculate_transition_strength_coherent(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}, polarization::SphericalUnitVector)::ComplexF64 + h_dipole = tensor_dot(T⁽¹⁾(polarization), spectrum.hamiltonian_parts.dipole_relative) + return e' * h_dipole * g / (1 / sqrt(3)) +end + +""" + plot_transition_strengths(spectrum::Spectrum, g::State, frequency_range; polarization::Union{Int, SphericalUnitVector, Nothing}=nothing) + +Plot the frequencies and strengths of electric dipole transitions out of `g`, +with energy between `frequency_range[1]` and `frequency_range[2]`. + +The `polarization` keyword argument can be used to choose a specific microwave polarization. This defaults to `nothing`, +which returns the incoherent sum over ``σ-``, ``π``, and ``σ+``. If `polarization` is an `Int`, then it is interpreted as +the spherical component `p = -1:1` of the dipole operator (`p == -1` corresponds to ``σ-`` polarization). If +`polarization` is a `SphericalUnitVector`, then the polarization is interpreted as linear along that axis. + +This method calls [`find_transition_strengths`](@ref) internally. + +See also [`find_transition_strengths`](@ref), [`make_hamiltonian_parts`](@ref), [`make_krb_hamiltonian_parts`](@ref), +[`Spectrum`](@ref). +""" +function plot_transition_strengths( + spectrum::Spectrum, + g::State, + frequency_range; + polarization::Union{Int, SphericalUnitVector, Nothing}=nothing +) + transitions = + find_transition_strengths(spectrum, g, frequency_range; polarization = polarization) + + freqs = [t[1] for t in transitions] + strengths = [abs(t[2]) for t in transitions] + + Gadfly.plot( + x = freqs, + y = strengths, + Gadfly.Geom.hair, + Gadfly.Geom.point, + Gadfly.Guide.xlabel("Frequency (MHz)"), + Gadfly.Guide.ylabel("Transition strength"), + ) +end + +""" + calculate_dipolar_interaction(spectrum::Spectrum, g::State, e::State; p::Int = 0) + calculate_dipolar_interaction(spectrum::Spectrum, g::Vector{ComplexF64}, e::Vector{ComplexF64}; p::Int = 0) + +#TODO: Document `calculate_dipolar_interaction` +""" +function calculate_dipolar_interaction( + spectrum::Spectrum, + g::State, + e::State; + p::Int = 0 +) + (overlap, index_g) = max_overlap_with(spectrum, g) + if overlap < 0.5 + @warn "The best overlap with your requested ground state is < 0.5." + end + + (overlap, index_e) = max_overlap_with(spectrum, e) + if overlap < 0.5 + @warn "The best overlap with your requested excited state is < 0.5." + end + + return calculate_dipolar_interaction( + spectrum, + get_eigenstate(spectrum, index_g), + get_eigenstate(spectrum, index_e); + p=p + ) +end + +function calculate_dipolar_interaction( + spectrum::Spectrum, + g::Vector{ComplexF64}, + e::Vector{ComplexF64}; + p::Int = 0 +) + d_1 = [Complex(calculate_transition_strength_coherent(spectrum, g, e, pol)) for pol=-1:1] + d_2 = [Complex(calculate_transition_strength_coherent(spectrum, e, g, pol)) for pol=-1:1] + return get_tensor_component(p, T⁽²⁾(d_1, d_2)) * sqrt(6) / 2 +end + +end # module diff --git a/src/fields.jl b/src/fields.jl index 6b8dfc4..85b1d89 100644 --- a/src/fields.jl +++ b/src/fields.jl @@ -1,338 +1,358 @@ -""" - SphericalVector(magnitude, θ, φ) - -Construct a vector with `magnitude`, polar angle `θ`, and azimuthal angle `φ`. - -Represents an external field vector in spherical coordinates used -to construct [`ExternalFields`](@ref). Currently, `SphericalVector`s -may be negated but no other mathematical operations are implemented. - -Vectors along ``x``, ``y``, or ``z`` can be quickly constructed -using [`VectorX`](@ref), [`VectorY`](@ref), and [`VectorZ`](@ref), -respectively. - -See also [`SphericalUnitVector`](@ref), [`ExternalFields`](@ref). -""" -struct SphericalVector - "Magnitude\n" - magnitude::Float64 - "Polar angle (rad)\n" - θ::Float64 - "Azimuthal angle (rad)\n" - φ::Float64 - - function SphericalVector(magnitude, θ, φ) - if magnitude < 0 - error("Magnitude must be nonnegative") - elseif θ < 0 || θ > π - error("Polar angle must be in [0, π]") - else - φr = rem2pi(φ, RoundDown) - if φr != φ - @warn "φ was provided outside of [0, 2π)" - end - return new(magnitude, θ, φr) - end - end -end - -function Base.:-(sv::SphericalVector) - if sv.magnitude == 0 - return sv - else - return SphericalVector(sv.magnitude, π - sv.θ, sv.φ + π) - end -end - -""" - VectorX(magnitude) - -Construct a [`SphericalVector`](@ref) with `magnitude` along `x`. - -# Examples -```jldoctest -julia> VectorX(5.0) -SphericalVector(5.0, 1.5707963267948966, 0.0) -``` -""" -VectorX(magnitude) = SphericalVector(magnitude, π / 2, 0) - -""" - VectorY(magnitude) - -Construct a [`SphericalVector`](@ref) with `magnitude` along `y`. - -# Examples -```jldoctest -julia> VectorY(2.0) -SphericalVector(2.0, 1.5707963267948966, 1.5707963267948966) -``` -""" -VectorY(magnitude) = SphericalVector(magnitude, π / 2, π / 2) - -""" - VectorZ(magnitude) - -Construct a [`SphericalVector`](@ref) with `magnitude` along `z`. - -# Examples -```jldoctest -julia> VectorZ(10.0) -SphericalVector(10.0, 0.0, 0.0) -``` -""" -VectorZ(magnitude) = SphericalVector(magnitude, 0, 0) - -""" - SphericalUnitVector(magnitude, θ, φ) - SphericalUnitVector(v::SphericalVector) - -Construct a unit vector with polar angle `θ`, and azimuthal angle `φ`. - -Represents the direction of an external field vector in spherical coordinates. -Currently, `SphericalUnitVector`s may be negated but no other -mathematical operations are implemented. - -Vectors along ``x``, ``y``, or ``z`` can be quickly constructed -using [`UnitVectorX`](@ref), [`UnitVectorY`](@ref), and [`UnitVectorZ`](@ref), -respectively. - -See also [`SphericalVector`](@ref), [`transition_strengths`](@ref), -[`T⁽¹⁾`](@ref), [`T⁽²⁾`](@ref). -""" -struct SphericalUnitVector - "Polar angle (rad)\n" - θ::Float64 - "Azimuthal angle (rad)\n" - φ::Float64 - - function SphericalUnitVector(θ, φ) - if θ < 0 || θ > π - error("Polar angle must be in [0, π]") - else - φr = rem2pi(φ, RoundDown) - if φr != φ - @warn "φ was provided outside of [0, 2π)" - end - return new(θ, φr) - end - end - - SphericalUnitVector(v::SphericalVector) = new(v.θ, v.φ) -end - -Base.:-(uv::SphericalUnitVector) = SphericalUnitVector(π - uv.θ, uv.φ + π) - -""" - UnitVectorX() - -Construct a [`SphericalUnitVector`](@ref) along `x`. - -# Examples -```jldoctest -julia> UnitVectorX() -SphericalUnitVector(1.5707963267948966, 0.0) -``` -""" -UnitVectorX() = SphericalUnitVector(π / 2, 0) - -""" - UnitVectorY() - -Construct a [`SphericalUnitVector`](@ref) along `y`. - -# Examples -```jldoctest -julia> UnitVectorY() -SphericalUnitVector(1.5707963267948966, 1.5707963267948966) -``` -""" -UnitVectorY() = SphericalUnitVector(π / 2, π / 2) - -""" - UnitVectorZ() - -Construct a [`SphericalUnitVector`](@ref) along `z`. - -# Examples -```jldoctest -julia> UnitVectorZ() -SphericalUnitVector(0.0, 0.0) -``` -""" -UnitVectorZ() = SphericalUnitVector(0, 0) - -""" - T⁽¹⁾(v) - -Construct the components of the rank 1 spherical tensor ``T⁽¹⁾(v)``. - -# Examples -```jldoctest -julia> T⁽¹⁾(UnitVectorX()) -3-element StaticArrays.SVector{3, ComplexF64} with indices SOneTo(3): - 0.7071067811865475 - 0.0im - 6.123233995736766e-17 + 0.0im - -0.7071067811865475 - 0.0im -``` - -```jldoctest -julia> T⁽¹⁾(UnitVectorY()) -3-element StaticArrays.SVector{3, ComplexF64} with indices SOneTo(3): - 4.329780281177466e-17 - 0.7071067811865475im - 6.123233995736766e-17 + 0.0im - -4.329780281177466e-17 - 0.7071067811865475im -``` - -```jldoctest -julia> T⁽¹⁾(UnitVectorZ()) -3-element StaticArrays.SVector{3, ComplexF64} with indices SOneTo(3): - 0.0 - 0.0im - 1.0 + 0.0im - -0.0 - 0.0im -``` -""" -function T⁽¹⁾(v::SphericalUnitVector)::SVector{3,ComplexF64} - θ = v.θ - φ = v.φ - - x = sin(θ) * cos(φ) - y = sin(θ) * sin(φ) - z = cos(θ) - - T11 = -(1 / sqrt(2)) * (x + im * y) - T10 = z - - return SVector(-conj(T11), T10, T11) -end - -""" - T⁽²⁾(v) - -Construct the components of the rank 2 spherical tensor ``T⁽²⁾(v, v)``. - -# Examples -```jldoctest -julia> T⁽²⁾(UnitVectorX()) -5-element StaticArrays.SVector{5, ComplexF64} with indices SOneTo(5): - 0.5 - 0.0im - 6.123233995736766e-17 - 0.0im - -0.4082482904638631 + 0.0im - -6.123233995736766e-17 - 0.0im - 0.5 + 0.0im -``` - -```jldoctest -julia> T⁽²⁾(UnitVectorY()) -5-element StaticArrays.SVector{5, ComplexF64} with indices SOneTo(5): - -0.5 - 6.123233995736766e-17im - 3.749399456654644e-33 - 6.123233995736766e-17im - -0.4082482904638631 + 0.0im - -3.749399456654644e-33 - 6.123233995736766e-17im - -0.5 + 6.123233995736766e-17im -``` - -```jldoctest -julia> T⁽²⁾(UnitVectorZ()) -5-element StaticArrays.SVector{5, ComplexF64} with indices SOneTo(5): - 0.0 - 0.0im - 0.0 - 0.0im - 0.8164965809277261 + 0.0im - -0.0 - 0.0im - 0.0 + 0.0im -``` -""" -function T⁽²⁾(v::SphericalUnitVector)::SVector{5,ComplexF64} - θ = v.θ - φ = v.φ - - x = sin(θ) * cos(φ) - y = sin(θ) * sin(φ) - z = cos(θ) - - T20 = (2 * z^2 - x^2 - y^2) / sqrt(6) - T21 = -(1 / 2) * (x * z + z * x + im * (y * z + z * y)) - T22 = (1 / 2) * (x * x - y * y + im * (x * y + y * x)) - - return SVector(conj(T22), -conj(T21), T20, T21, T22) -end - -function T⁽²⁾(T1a, T1b)::SVector{5,ComplexF64} - (T1am1, T1a0, T1ap1) = T1a - (T1bm1, T1b0, T1bp1) = T1b - - T22 = T1ap1 * T1bp1 - T21 = (1/sqrt(2)) * (T1ap1 * T1b0 + T1a0 * T1bp1) - T20 = (1/sqrt(6)) * (2*T1a0*T1b0 + T1ap1 * T1bm1 + T1am1 * T1bp1) - T2m1 = (1/sqrt(2)) * (T1am1 * T1b0 + T1a0 * T1bm1) - T2m2 = T1am1 * T1bm1 - return SVector(T2m2, T2m1, T20, T21, T22) -end - -function get_tensor_component(p::Int, tensor) - rank::Int = (length(tensor) - 1) // 2 # Length should be 2*k + 1 - return tensor[1+(p+rank)] -end - -""" - tensor_dot(a, b) - -Contract two spherical tensors `a` and `b`. -""" -function tensor_dot(a, b) - @assert size(a, 1) == size(b, 1) - @assert isodd(size(a, 1)) - - mapreduce(p -> conj(a[p]) .* b[p], +, eachindex(a)) -end - -""" - ExternalFields(B::SphericalVector, E::SphericalVector, Optical::Vector{SphericalVector}) - ExternalFields(B::Float64, E::Float64) - -External magnetic, electric, optical fields to use in constructing the Hamiltonian. - -If `B` and `E` are provided as `Float64`s, then the fields are assumed to be along `z`. -The `Optical` argument can also be left as an empty vector `[]`. - -See also [`calculate_spectrum`](@ref), [`hamiltonian`](@ref), [`SphericalVector`](@ref). - -# Examples -```jldoctest -julia> ExternalFields(VectorZ(545.9), VectorX(1020.0), []) -ExternalFields(SphericalVector(545.9, 0.0, 0.0), SphericalVector(1020.0, 1.5707963267948966, 0.0), SphericalVector[]) -``` - -```jldoctest -julia> ExternalFields(545.9, 1020.0) -ExternalFields(SphericalVector(545.9, 0.0, 0.0), SphericalVector(1020.0, 0.0, 0.0), SphericalVector[]) -``` - -```jldoctest -julia> ExternalFields(VectorZ(545.9), VectorX(1020.0), [VectorY(2300.), SphericalVector(2300., π/2, π/4)]) -[...] -``` -""" -struct ExternalFields - "Magnetic field (G)\n" - B::SphericalVector - "Electric field (V/cm)\n" - E::SphericalVector - "Laser fields (W/cm^2)\n" - Optical::Vector{SphericalVector} - - ExternalFields(B::SphericalVector, E::SphericalVector, Optical) = new(B, E, Optical) -end - -ExternalFields() = ExternalFields(0.0, 0.0, []) -ExternalFields(B, E) = ExternalFields(B, E, []) -ExternalFields(B::Float64, E, Optical) = ExternalFields(VectorZ(B), E, Optical) -ExternalFields(B, E::Float64, Optical) = ExternalFields(B, VectorZ(E), Optical) -ExternalFields(B::Float64, E::Float64, Optical) = ExternalFields(VectorZ(B), VectorZ(E), Optical) - -const DEFAULT_FIELDS = ExternalFields(545.9, 0.0) -const TEST_FIELDS = ExternalFields( - SphericalVector(545.9, π / 4, π / 4), - SphericalVector(1000.0, 3π / 4, 7π / 4), - [SphericalVector(2350.0, 0.0, 0.0)], -) +""" + SphericalVector(magnitude, θ, φ) + +Construct a vector with `magnitude`, polar angle `θ`, and azimuthal angle `φ`. + +Represents an external field vector in spherical coordinates used +to construct [`ExternalFields`](@ref). Currently, `SphericalVector`s +may be negated but no other mathematical operations are implemented. + +Vectors along ``x``, ``y``, or ``z`` can be quickly constructed +using [`VectorX`](@ref), [`VectorY`](@ref), and [`VectorZ`](@ref), +respectively. + +See also [`SphericalUnitVector`](@ref), [`ExternalFields`](@ref). +""" +struct SphericalVector + "Magnitude\n" + magnitude::Float64 + "Polar angle (rad)\n" + θ::Float64 + "Azimuthal angle (rad)\n" + φ::Float64 + + function SphericalVector(magnitude, θ, φ) + if magnitude < 0 + error("Magnitude must be nonnegative") + elseif θ < 0 || θ > π + error("Polar angle must be in [0, π]") + else + φr = rem2pi(φ, RoundDown) + if φr != φ + @warn "φ was provided outside of [0, 2π)" + end + return new(magnitude, θ, φr) + end + end +end + +function Base.:-(sv::SphericalVector) + if sv.magnitude == 0 + return sv + else + return SphericalVector(sv.magnitude, π - sv.θ, sv.φ + π) + end +end + +""" + VectorX(magnitude) + +Construct a [`SphericalVector`](@ref) with `magnitude` along `x`. + +# Examples +```jldoctest +julia> VectorX(5.0) +SphericalVector(5.0, 1.5707963267948966, 0.0) +``` +""" +VectorX(magnitude) = SphericalVector(magnitude, π / 2, 0) + +""" + VectorY(magnitude) + +Construct a [`SphericalVector`](@ref) with `magnitude` along `y`. + +# Examples +```jldoctest +julia> VectorY(2.0) +SphericalVector(2.0, 1.5707963267948966, 1.5707963267948966) +``` +""" +VectorY(magnitude) = SphericalVector(magnitude, π / 2, π / 2) + +""" + VectorZ(magnitude) + +Construct a [`SphericalVector`](@ref) with `magnitude` along `z`. + +# Examples +```jldoctest +julia> VectorZ(10.0) +SphericalVector(10.0, 0.0, 0.0) +``` +""" +VectorZ(magnitude) = SphericalVector(magnitude, 0, 0) + +""" + SphericalUnitVector(magnitude, θ, φ) + SphericalUnitVector(v::SphericalVector) + +Construct a unit vector with polar angle `θ`, and azimuthal angle `φ`. + +Represents the direction of an external field vector in spherical coordinates. +Currently, `SphericalUnitVector`s may be negated but no other +mathematical operations are implemented. + +Vectors along ``x``, ``y``, or ``z`` can be quickly constructed +using [`UnitVectorX`](@ref), [`UnitVectorY`](@ref), and [`UnitVectorZ`](@ref), +respectively. + +See also [`SphericalVector`](@ref), [`find_transition_strengths`](@ref), +[`T⁽¹⁾`](@ref), [`T⁽²⁾`](@ref). +""" +struct SphericalUnitVector + "Polar angle (rad)\n" + θ::Float64 + "Azimuthal angle (rad)\n" + φ::Float64 + + function SphericalUnitVector(θ, φ) + if θ < 0 || θ > π + error("Polar angle must be in [0, π]") + else + φr = rem2pi(φ, RoundDown) + if φr != φ + @warn "φ was provided outside of [0, 2π)" + end + return new(θ, φr) + end + end + + SphericalUnitVector(v::SphericalVector) = new(v.θ, v.φ) +end + +Base.:-(uv::SphericalUnitVector) = SphericalUnitVector(π - uv.θ, uv.φ + π) + +""" + UnitVectorX() + +Construct a [`SphericalUnitVector`](@ref) along `x`. + +# Examples +```jldoctest +julia> UnitVectorX() +SphericalUnitVector(1.5707963267948966, 0.0) +``` +""" +UnitVectorX() = SphericalUnitVector(π / 2, 0) + +""" + UnitVectorY() + +Construct a [`SphericalUnitVector`](@ref) along `y`. + +# Examples +```jldoctest +julia> UnitVectorY() +SphericalUnitVector(1.5707963267948966, 1.5707963267948966) +``` +""" +UnitVectorY() = SphericalUnitVector(π / 2, π / 2) + +""" + UnitVectorZ() + +Construct a [`SphericalUnitVector`](@ref) along `z`. + +# Examples +```jldoctest +julia> UnitVectorZ() +SphericalUnitVector(0.0, 0.0) +``` +""" +UnitVectorZ() = SphericalUnitVector(0, 0) + +""" + T⁽¹⁾(v) + +Construct the components of the rank 1 spherical tensor ``T⁽¹⁾(v)``. + +# Examples +```jldoctest +julia> T⁽¹⁾(UnitVectorX()) +3-element StaticArrays.SVector{3, ComplexF64} with indices SOneTo(3): + 0.7071067811865475 - 0.0im + 6.123233995736766e-17 + 0.0im + -0.7071067811865475 - 0.0im +``` + +```jldoctest +julia> T⁽¹⁾(UnitVectorY()) +3-element StaticArrays.SVector{3, ComplexF64} with indices SOneTo(3): + 4.329780281177466e-17 - 0.7071067811865475im + 6.123233995736766e-17 + 0.0im + -4.329780281177466e-17 - 0.7071067811865475im +``` + +```jldoctest +julia> T⁽¹⁾(UnitVectorZ()) +3-element StaticArrays.SVector{3, ComplexF64} with indices SOneTo(3): + 0.0 - 0.0im + 1.0 + 0.0im + -0.0 - 0.0im +``` +""" +function T⁽¹⁾(v::SphericalUnitVector)::SVector{3,ComplexF64} + θ = v.θ + φ = v.φ + + x = sin(θ) * cos(φ) + y = sin(θ) * sin(φ) + z = cos(θ) + + T11 = -(1 / sqrt(2)) * (x + im * y) + T10 = z + + return SVector(-conj(T11), T10, T11) +end + +""" + T⁽²⁾(v) + +Construct the components of the rank 2 spherical tensor ``T⁽²⁾(v, v)``. + +# Examples +```jldoctest +julia> T⁽²⁾(UnitVectorX()) +5-element StaticArrays.SVector{5, ComplexF64} with indices SOneTo(5): + 0.5 - 0.0im + 6.123233995736766e-17 - 0.0im + -0.4082482904638631 + 0.0im + -6.123233995736766e-17 - 0.0im + 0.5 + 0.0im +``` + +```jldoctest +julia> T⁽²⁾(UnitVectorY()) +5-element StaticArrays.SVector{5, ComplexF64} with indices SOneTo(5): + -0.5 - 6.123233995736766e-17im + 3.749399456654644e-33 - 6.123233995736766e-17im + -0.4082482904638631 + 0.0im + -3.749399456654644e-33 - 6.123233995736766e-17im + -0.5 + 6.123233995736766e-17im +``` + +```jldoctest +julia> T⁽²⁾(UnitVectorZ()) +5-element StaticArrays.SVector{5, ComplexF64} with indices SOneTo(5): + 0.0 - 0.0im + 0.0 - 0.0im + 0.8164965809277261 + 0.0im + -0.0 - 0.0im + 0.0 + 0.0im +``` +""" +function T⁽²⁾(v::SphericalUnitVector)::SVector{5,ComplexF64} + θ = v.θ + φ = v.φ + + x = sin(θ) * cos(φ) + y = sin(θ) * sin(φ) + z = cos(θ) + + T20 = (2 * z^2 - x^2 - y^2) / sqrt(6) + T21 = -(1 / 2) * (x * z + z * x + im * (y * z + z * y)) + T22 = (1 / 2) * (x * x - y * y + im * (x * y + y * x)) + + return SVector(conj(T22), -conj(T21), T20, T21, T22) +end + +function T⁽²⁾(T1a, T1b)::SVector{5,ComplexF64} + (T1am1, T1a0, T1ap1) = T1a + (T1bm1, T1b0, T1bp1) = T1b + + T22 = T1ap1 * T1bp1 + T21 = (1/sqrt(2)) * (T1ap1 * T1b0 + T1a0 * T1bp1) + T20 = (1/sqrt(6)) * (2*T1a0*T1b0 + T1ap1 * T1bm1 + T1am1 * T1bp1) + T2m1 = (1/sqrt(2)) * (T1am1 * T1b0 + T1a0 * T1bm1) + T2m2 = T1am1 * T1bm1 + return SVector(T2m2, T2m1, T20, T21, T22) +end + +""" + get_tensor_component(p::Int, tensor) + +#TODO: Document what `get_tensor_component` does. +""" +function get_tensor_component(p::Int, tensor) + rank::Int = (length(tensor) - 1) // 2 # Length should be 2*k + 1 + return tensor[1+(p+rank)] +end + +""" + tensor_dot(a, b) + +Contract two spherical tensors `a` and `b`. +""" +function tensor_dot(a, b) + @assert size(a, 1) == size(b, 1) + @assert isodd(size(a, 1)) + + mapreduce(p -> conj(a[p]) .* b[p], +, eachindex(a)) +end + +""" + ExternalFields(B::SphericalVector, E::SphericalVector, Optical::Vector{SphericalVector}) + ExternalFields(B::Float64, E::Float64) + +External magnetic, electric, optical fields to use in constructing the Hamiltonian. + +If `B` and `E` are provided as `Float64`s, then the fields are assumed to be along `z`. +The `Optical` argument can also be left as an empty vector `[]`. + +See also [`calculate_spectrum`](@ref), [`hamiltonian`](@ref), [`SphericalVector`](@ref). + +# Examples +```jldoctest; output = false +julia> ExternalFields(VectorZ(545.9), VectorX(1020.0), []) +ExternalFields(SphericalVector(545.9, 0.0, 0.0), SphericalVector(1020.0, 1.5707963267948966, 0.0), SphericalVector[]) +``` + +```jldoctest; output = false +julia> ExternalFields(545.9, 1020.0) +ExternalFields(SphericalVector(545.9, 0.0, 0.0), SphericalVector(1020.0, 0.0, 0.0), SphericalVector[]) +``` + +```jldoctest; output = false +julia> ExternalFields(VectorZ(545.9), VectorX(1020.0), [VectorY(2300.), SphericalVector(2300., π/2, π/4)]) +ExternalFields(SphericalVector(545.9, 0.0, 0.0), SphericalVector(1020.0, 1.5707963267948966, 0.0), SphericalVector[SphericalVector(2300.0, 1.5707963267948966, 1.5707963267948966), SphericalVector(2300.0, 1.5707963267948966, 0.7853981633974483)]) +``` +""" +struct ExternalFields + "Magnetic field (G)\n" + B::SphericalVector + "Electric field (V/cm)\n" + E::SphericalVector + "Laser fields (W/cm^2)\n" + Optical::Vector{SphericalVector} + + ExternalFields(B::SphericalVector, E::SphericalVector, Optical) = new(B, E, Optical) +end + +ExternalFields() = ExternalFields(0.0, 0.0, []) +ExternalFields(B, E) = ExternalFields(B, E, []) +ExternalFields(B::Float64, E, Optical) = ExternalFields(VectorZ(B), E, Optical) +ExternalFields(B, E::Float64, Optical) = ExternalFields(B, VectorZ(E), Optical) +ExternalFields(B::Float64, E::Float64, Optical) = ExternalFields(VectorZ(B), VectorZ(E), Optical) + +""" + DEFAULT_FIELDS = ExternalFields(545.9, 0.0) + +Defines an [`ExternalFields`](@ref) with a magnetic field of 545.9 G along the `z` axis. +""" +const DEFAULT_FIELDS = ExternalFields(545.9, 0.0) + +""" + TEST_FIELDS = ExternalFields( + SphericalVector(545.9, π / 4, π / 4), + SphericalVector(1000.0, 3π / 4, 7π / 4), + [SphericalVector(2350.0, 0.0, 0.0)], + ) + +Defines an [`ExternalFields`](@ref) with a magnetic, electric, and optical fields for tests. +""" +const TEST_FIELDS = ExternalFields( + SphericalVector(545.9, π / 4, π / 4), + SphericalVector(1000.0, 3π / 4, 7π / 4), + [SphericalVector(2350.0, 0.0, 0.0)], +) diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index cac1882..9f76c36 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -128,7 +128,18 @@ end const SparseHamiltonian = SparseMatrixCSC{ComplexF64, Int64} """ - HamiltonianParts +``` +struct HamiltonianParts + basis::Vector{State} + rotation::SparseHamiltonian + dipole::SVector{3, SparseHamiltonian} + dipole_relative::SVector{3, SparseHamiltonian} # used for transition strengths + hyperfine::SparseHamiltonian + zeeman::SVector{3, SparseHamiltonian} + ac_scalar::SparseHamiltonian + ac_tensor::SVector{5, SparseHamiltonian} +end +``` Contains all parts of the Hamiltonian except external fields. @@ -217,6 +228,8 @@ end Construct the full Hamiltonian including magnetic, electric, and optical fields. +Returns a dense Hermitian matrix. + The field-independent building blocks in `parts` can be reused over calls to `hamiltonian` to avoid recalculating the matrix elements each time. diff --git a/src/molecular_parameters.jl b/src/molecular_parameters.jl index ef3383f..b2d5dda 100644 --- a/src/molecular_parameters.jl +++ b/src/molecular_parameters.jl @@ -1,124 +1,201 @@ -""" - ZeemanParameters - -Contains the g-factors and nuclear shielding factors -for computing Zeeman shifts. -""" -struct ZeemanParameters - "Rotational g factor\n" - gᵣ::Float64 - "Nuclear g factor\n" - gᵢ::SVector{2,Float64} - "Nuclear shielding factor\n" - σᵢ::SVector{2,Float64} -end - -""" - NuclearParameters - -Contains the nuclear electric quadrupole moments, -nuclear spin-rotation couplings, and nuclear -spin-spin coupling for calculating the hyperfine -Hamiltonian. -""" -struct NuclearParameters - "Nuclear electric quadrupole (MHz)\n" - eqQᵢ::SVector{2,Float64} - "Nuclear spin-rotation interaction (MHz)\n" - cᵢ::SVector{2,Float64} - "Nuclear spin-spin scalar interaction (MHz)\n" - c₄::Float64 -end - -""" - Polarizability - -Contains the parallel and perpendicular ac -polarizabilities at a particular optical wavelength λ. - -So far, it is assumed that λ is far-detuned from -any electronic transitions, such that the polarizability -does not depend on rotational state ``N``. -""" -struct Polarizability - "Parallel ac polarizability (MHz / (W / cm^2))\n" - α_par::Float64 - "Perpendicular ac polarizability (MHz / (W / cm^2))\n" - α_perp::Float64 -end - -""" - MolecularParameters - -Contains the coupling constants for the molecular -Hamiltonian and the nuclear angular momenta (needed to -construct the basis states). -""" -struct MolecularParameters - "Permanent dipole moment (Debye)\n" - dₚ::Float64 - "Rotational constant (MHz)\n" - Bᵣ::Float64 - "Nuclear angular momenta\n" - I::SVector{2,HalfInt} - "Zeeman parameters\n" - zeeman::ZeemanParameters - "Nuclear Parameters\n" - nuclear::NuclearParameters - "Molecular polarizability at the trapping wavelength\n" - α::Polarizability -end - -"Theoretical values from Aldegunde et al. PRA (2008)\n" # TODO: Put exact citation -const KRb_Zeeman = ZeemanParameters(0.014, [-0.324, 1.834], [1321e-6, 3469e-6]) - -"Experimental values from Neyenhuis et al., PRL 109, 230403 (2012)\n" -const KRb_Polarizability = Polarizability(10.0e-5, 3.3e-5) - -"Experimental values from Neyenhuis et al., PRL 109, 230403 (2012)\n" -const KRb_Nuclear_Neyenhuis = - NuclearParameters([0.45, -1.308], [-24.1e-6, 420.1e-6], -2030.4e-6) - -"Experimental values from Ospelkaus et al., PRL 104, 030402 (2010)\n" -const KRb_Nuclear_Ospelkaus = - NuclearParameters([0.45, -1.41], [-24.1e-6, 420.1e-6], -2030.4e-6) - -"Experimental values from Neyenhuis et al., PRL 109, 230403 (2012)\n" -const KRb_Parameters_Neyenhuis = MolecularParameters( - 0.574, - 1113.9514, - [HalfInt(4), HalfInt(3 / 2)], - KRb_Zeeman, - KRb_Nuclear_Neyenhuis, - KRb_Polarizability, -) - -"Experimental values from Ospelkaus et al., PRL 104, 030402 (2010)\n" -const KRb_Parameters_Ospelkaus = MolecularParameters( - 0.574, - 1113.950, - [HalfInt(4), HalfInt(3 / 2)], - KRb_Zeeman, - KRb_Nuclear_Ospelkaus, - KRb_Polarizability, -) - -const DEFAULT_MOLECULAR_PARAMETERS = KRb_Parameters_Neyenhuis - -""" - TOY_MOLECULE_PARAMETERS - -Toy model values with dipole = 1 D, no hyperfine structure, etc. -Intended for testing. - -Note that the formulas break down for `I = 0`, which is why we use -`I = 1` here. -""" -const TOY_MOLECULE_PARAMETERS = MolecularParameters( - 1.0, - 1000.0, - [HalfInt(1), HalfInt(1)], # Some of the matrix elements don't make sense for I = 0 - ZeemanParameters(0.0, [0.0, 0.0], [0.0, 0.0]), - NuclearParameters([0.0, 0.0], [0.0, 0.0], 0.0), - Polarizability(0.0, 0.0), -) +""" +``` +struct ZeemanParameters + "Rotational g factor" + gᵣ::Float64 + "Nuclear g factor" + gᵢ::SVector{2,Float64} + "Nuclear shielding factor" + σᵢ::SVector{2,Float64} +end +``` + +Contains the g-factors and nuclear shielding factors +for computing Zeeman shifts. +""" +struct ZeemanParameters + "Rotational g factor\n" + gᵣ::Float64 + "Nuclear g factor\n" + gᵢ::SVector{2,Float64} + "Nuclear shielding factor\n" + σᵢ::SVector{2,Float64} +end + +""" +``` +struct NuclearParameters + "Nuclear electric quadrupole (MHz)" + eqQᵢ::SVector{2,Float64} + "Nuclear spin-rotation interaction (MHz)" + cᵢ::SVector{2,Float64} + "Nuclear spin-spin scalar interaction (MHz)" + c₄::Float64 +end +``` + +Contains the nuclear electric quadrupole moments, +nuclear spin-rotation couplings, and nuclear +spin-spin coupling for calculating the hyperfine +Hamiltonian. +""" +struct NuclearParameters + "Nuclear electric quadrupole (MHz)\n" + eqQᵢ::SVector{2,Float64} + "Nuclear spin-rotation interaction (MHz)\n" + cᵢ::SVector{2,Float64} + "Nuclear spin-spin scalar interaction (MHz)\n" + c₄::Float64 +end + +""" +``` +struct Polarizability + "Parallel ac polarizability (MHz / (W / cm^2))" + α_par::Float64 + "Perpendicular ac polarizability (MHz / (W / cm^2))" + α_perp::Float64 +end +``` + +Contains the parallel and perpendicular ac +polarizabilities at a particular optical wavelength λ. + +So far, it is assumed that λ is far-detuned from +any electronic transitions, such that the polarizability +does not depend on rotational state ``N``. +""" +struct Polarizability + "Parallel ac polarizability (MHz / (W / cm^2))\n" + α_par::Float64 + "Perpendicular ac polarizability (MHz / (W / cm^2))\n" + α_perp::Float64 +end + +""" +``` +struct MolecularParameters + "Permanent dipole moment (Debye)" + dₚ::Float64 + "Rotational constant (MHz)" + Bᵣ::Float64 + "Nuclear angular momenta" + I::SVector{2,HalfInt} + "Zeeman parameters" + zeeman::ZeemanParameters + "Nuclear Parameters" + nuclear::NuclearParameters + "Molecular polarizability at the trapping wavelength" + α::Polarizability +end +``` + +Contains the coupling constants for the molecular +Hamiltonian and the nuclear angular momenta (needed to +construct the basis states). +""" +struct MolecularParameters + "Permanent dipole moment (Debye)\n" + dₚ::Float64 + "Rotational constant (MHz)\n" + Bᵣ::Float64 + "Nuclear angular momenta\n" + I::SVector{2,HalfInt} + "Zeeman parameters\n" + zeeman::ZeemanParameters + "Nuclear Parameters\n" + nuclear::NuclearParameters + "Molecular polarizability at the trapping wavelength\n" + α::Polarizability +end + +" + KRb_Zeeman = ZeemanParameters(0.014, [-0.324, 1.834], [1321e-6, 3469e-6]) + +[`ZeemanParameters`](@ref) with theoretical values from [Aldegunde et al., PRA 78, 033434 (2008)](https://doi.org/10.1103/PhysRevA.78.033434)\n" +const KRb_Zeeman = ZeemanParameters(0.014, [-0.324, 1.834], [1321e-6, 3469e-6]) + +" + KRb_Polarizability = Polarizability(10.0e-5, 3.3e-5) + +[`Polarizability`](@ref) with experimental values from [Neyenhuis et al., PRL 109, 230403 (2012)](https://doi.org/10.1103/PhysRevLett.109.230403)\n" +const KRb_Polarizability = Polarizability(10.0e-5, 3.3e-5) + +" + KRb_Nuclear_Neyenhuis = NuclearParameters([0.45, -1.308], [-24.1e-6, 420.1e-6], -2030.4e-6) + +[`NuclearParameters`](@ref) with experimental values from [Neyenhuis et al., PRL 109, 230403 (2012)](https://doi.org/10.1103/PhysRevLett.109.230403)\n" +const KRb_Nuclear_Neyenhuis = + NuclearParameters([0.45, -1.308], [-24.1e-6, 420.1e-6], -2030.4e-6) + +" + KRb_Nuclear_Ospelkaus = NuclearParameters([0.45, -1.41], [-24.1e-6, 420.1e-6], -2030.4e-6) + +[`NuclearParameters`](@ref) with experimental values from [Ospelkaus et al., PRL 104, 030402 (2010)](https://doi.org/10.1103/PhysRevLett.104.030402)\n" +const KRb_Nuclear_Ospelkaus = + NuclearParameters([0.45, -1.41], [-24.1e-6, 420.1e-6], -2030.4e-6) + +" + KRb_Parameters_Neyenhuis = MolecularParameters( + 0.574, + 1113.9514, + [HalfInt(4), HalfInt(3 / 2)], + KRb_Zeeman, + KRb_Nuclear_Neyenhuis, + KRb_Polarizability, + ) + +[`MolecularParameters`](@ref) with experimental values from [Neyenhuis et al., PRL 109, 230403 (2012)](https://doi.org/10.1103/PhysRevLett.109.230403)\n" +const KRb_Parameters_Neyenhuis = MolecularParameters( + 0.574, + 1113.9514, + [HalfInt(4), HalfInt(3 / 2)], + KRb_Zeeman, + KRb_Nuclear_Neyenhuis, + KRb_Polarizability, +) + +" + KRb_Parameters_Ospelkaus = MolecularParameters( + 0.574, + 1113.950, + [HalfInt(4), HalfInt(3 / 2)], + KRb_Zeeman, + KRb_Nuclear_Ospelkaus, + KRb_Polarizability, + ) + +[`MolecularParameters`](@ref) with experimental values from [Ospelkaus et al., PRL 104, 030402 (2010)](https://doi.org/10.1103/PhysRevLett.104.030402)\n" +const KRb_Parameters_Ospelkaus = MolecularParameters( + 0.574, + 1113.950, + [HalfInt(4), HalfInt(3 / 2)], + KRb_Zeeman, + KRb_Nuclear_Ospelkaus, + KRb_Polarizability, +) + +""" + DEFAULT_MOLECULAR_PARAMETERS + +Default molecular parameters; alias for [`KRb_Parameters_Neyenhuis`](@ref) +""" +const DEFAULT_MOLECULAR_PARAMETERS = KRb_Parameters_Neyenhuis + +""" + TOY_MOLECULE_PARAMETERS + +Toy model values with dipole = 1 D, no hyperfine structure, etc. +Intended for testing. + +Note that the formulas break down for `I = 0`, which is why we use +`I = 1` here. +""" +const TOY_MOLECULE_PARAMETERS = MolecularParameters( + 1.0, + 1000.0, + [HalfInt(1), HalfInt(1)], # Some of the matrix elements don't make sense for I = 0 + ZeemanParameters(0.0, [0.0, 0.0], [0.0, 0.0]), + NuclearParameters([0.0, 0.0], [0.0, 0.0], 0.0), + Polarizability(0.0, 0.0), +)