Implementation of laser-dressed ion states

Project.toml

@@ -4,19 +4,24 @@ authors = ["Stefanos Carlström <[email protected]>"]
 version = "0.1.0"
 
 [deps]
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 ElectricFields = "2f84ce32-9bb1-11e8-0d9f-3dce90a4beca"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
 [compat]
 ElectricFields = "0.1, 0.2"
 FastGaussQuadrature = "0.4.7, 0.5"
+FillArrays = "1"
 HCubature = "1.5"
 ProgressMeter = "1.5,1.6"
 SpecialFunctions = "1.3,2"
@@ -26,7 +31,8 @@ UnitfulAtomic = "1.0"
 julia = "1.6"
 
 [extras]
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["DelimitedFiles", "Test"]

docs/Project.toml

@@ -10,8 +10,8 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
 [compat]
 CondaPkg = "0.2"
-Documenter = "0.26"
-ElectricFields = "0.1"
+Documenter = "0.27"
+ElectricFields = "0.1, 0.2"
 FFTW = "1.4.0"
 PythonCall = "0.9"
 PythonPlot = "1"

docs/make.jl

@@ -32,6 +32,7 @@ makedocs(;
                     :Hamiltonian => "\\operator{H}",
                     :hamiltonian => "\\operator{h}",
                     :Lagrangian => "\\operator{L}",
+                    :propU => "\\mathcal{U}",
                     :fock => "\\operator{f}",
                     :lagrange => ["\\epsilon_{#1}", 1],
                     :vary => ["\\delta_{#1}", 1],

docs/plots.jl

@@ -67,10 +67,12 @@ function hhg_example()
 
     # d will be a vector of scalars; by limiting the "memory" of the
     # integrals, we can include only the short trajectory.
-    d = induced_dipole(Iₚ, F, ndt, memory=floor(Int, 0.65ndt));
+    d = induced_dipole(Iₚ, F, ndt, window=flat_window(floor(Int, 0.65ndt)));
+    # d_all includes all trajectories.
+    d_all = induced_dipole(Iₚ, F, ndt);
 
     # d2 will be a vector of 3d vectors
-    d2 = induced_dipole(Iₚ, F2, ndt, memory=floor(Int, 0.65ndt));
+    d2 = induced_dipole(Iₚ, F2, ndt, window=flat_window(floor(Int, 0.65ndt)));
     d2x = [e[1] for e in d2]
     d2y = [e[2] for e in d2]
     d2z = [e[3] for e in d2]
@@ -80,6 +82,7 @@ function hhg_example()
     q = fftshift(fftfreq(length(t), ndt))
     qsel = ind(q,0):ind(q,100)
     D = spectrum(d)
+    D_all = spectrum(d_all)
     D2 = spectrum(d2)
     cutoff = 3.17austrip(ponderomotive_potential(F)) + Iₚ
 
@@ -98,13 +101,15 @@ function hhg_example()
             plot(tplot, d2x, "--")
             plot(tplot, d2y, "--")
             plot(tplot, d2z, "--")
-            legend(["1d", "3d x", "3d y", "3d z"])
+            plot(tplot, d_all, ":")
+            legend(["1d", "3d x", "3d y", "3d z", "1d all traj."], loc=1)
             ylabel(L"\langle\mathbf{r}\rangle(t)")
         end
         csubplot(313) do
             semilogy(q[qsel], abs2.(D[qsel]))
             semilogy(q[qsel], abs2.(D2[qsel,:]), "--")
-            legend(["1d", "3d x", "3d y", "3d z"])
+            semilogy(q[qsel], abs2.(D_all[qsel]), ":")
+            legend(["1d", "3d x", "3d y", "3d z", "1d all traj."])
             axvline(cutoff/photon_energy(F), linestyle=":", color="black")
             xlabel(L"Harmonic order of 800 nm [$q$]")
             ylabel(L"|\mathbf{r}(q)|^2")

docs/src/index.md

@@ -12,7 +12,7 @@ throughout, unless otherwise specified.
 
 At the moment, two kind of observables are supported: induced dipole
 moment and photoelectron spectra. Time integrals are performed
-numerically, with recursive time integrals limited by a `memory`,
+numerically, with recursive time integrals limited by a `window`,
 i.e. how many time steps are considered (default is from the beginning
 of the pulse). Integrals over intermediate photoelectron momenta are
 performed using the saddle-point method, i.e. given two times, the
@@ -75,7 +75,7 @@ julia> Iₚ = 0.5 # Hydrogen
 
 julia> # d will be a vector of scalars; by limiting the "memory" of the
        # integrals, we can include only the short trajectory.
-       d = induced_dipole(Iₚ, F, ndt, memory=floor(Int, 0.65ndt));
+       d = induced_dipole(Iₚ, F, ndt, window=flat_window(floor(Int, 0.65ndt)));
 ┌ Info: Induced dipole calculation
 │   system =
 │    1-channel System:
@@ -91,8 +91,25 @@ julia> # d will be a vector of scalars; by limiting the "memory" of the
 └
 Progress: 100%|████████████████████████████████████████████████████████████████| Time: 0:00:00
 
+julia> # d_all includes all trajectories.
+       d_all = induced_dipole(Iₚ, F, ndt);
+┌ Info: Induced dipole calculation
+│   system =
+│    1-channel SFA System:
+│     1. IonizationChannel: Iₚ = 0.5 Ha = 13.6055 eV
+│
+│   diagram =
+│    Goldstone Diagram:
+│       |0⟩
+│       ╱ ╲⇜
+│     1┃   │𝐩
+│       ╲ ╱⇝
+│       |0⟩
+└
+Progress: 100%|████████████████████████████████████████████████████████████████| Time: 0:00:08
+
 julia> # d2 will be a vector of 3d vectors
-       d2 = induced_dipole(Iₚ, F2, ndt, memory=floor(Int, 0.65ndt));
+       d2 = induced_dipole(Iₚ, F2, ndt, window=flat_window(floor(Int, 0.65ndt)));
 ┌ Info: Induced dipole calculation
 │   system =
 │    1-channel System:

src/StrongFieldApproximation.jl

@@ -6,22 +6,36 @@ using UnitfulAtomic
 using FastGaussQuadrature
 
 using LinearAlgebra
+using SparseArrays
 using StaticArrays
+using FillArrays
+
+import DataStructures: Queue, enqueue!, dequeue!
 
 using ProgressMeter
+using TimerOutputs
 
 include("threading.jl")
 include("telescope_iterators.jl")
 
+include("find_blocks.jl")
+
 include("atom_models.jl")
 include("momentum_grid.jl")
 
 include("units.jl")
 include("field.jl")
+include("ion_propagators.jl")
 include("volkov.jl")
+include("windows.jl")
 
 include("channels.jl")
+include("system.jl")
+include("diagrams.jl")
+include("momenta.jl")
+
 include("dyson_expansions.jl")
+include("multi_channel_sfa.jl")
 
 include("ionization_rates.jl")
 

src/channels.jl

@@ -21,15 +21,26 @@ struct DipoleSourceTerm{Dipole<:Function,Field<:Union{ElectricFields.AbstractFie
         new{Dipole,FieldAmplitude}(d, Fv)
 end
 
+Base.show(io::IO, d::DipoleSourceTerm) =
+    write(io, "DipoleSourceTerm")
+
+function Base.show(io::IO, mime::MIME"text/plain", d::DipoleSourceTerm)
+    show(io, d)
+    print(io, "\nDipole: ")
+    show(io, mime, d.d)
+    print(io, "\nField: ")
+    show(io, d.F)
+end
+
 dipole(F::Number, d::Number) = F*d
-dipole(F::SVector{3}, d::SVector{3}) = dot(F, d)
+dipole(F::SVector{3}, d::SVector{3}) = F[1]*d[1] + F[2]*d[2] + F[3]*d[3]
 dipole(F::Number, d::SVector{3}) = F*d[3]
 
-source_term(d::DipoleSourceTerm{<:Any,<:ElectricFields.AbstractField}, t, p) =
-    dipole(field_amplitude(d.F, t), d.d(p))
+_field_amplitude(F::ElectricFields.AbstractField, t) =
+    field_amplitude(F, t)
+_field_amplitude(F::AbstractVector, i::Integer) = F[i]
 
-source_term(d::DipoleSourceTerm{<:Any,<:AbstractVector}, i::Integer, p) =
-    dipole(d.F[i], d.d(p))
+source_term(d::DipoleSourceTerm, t, p) = dipole(_field_amplitude(d.F, t), d.d(p))
 
 # * Ionizations channels
 
@@ -43,6 +54,8 @@ Represents an ionization channel with energy `E` above the neutral
 struct IonizationChannel{T,SourceTerm<:AbstractSourceTerm}
     E::T
     st::SourceTerm
+    IonizationChannel(E::T, st::SourceTerm) where {T,SourceTerm} =
+        new{T,SourceTerm}(E, st)
 end
 
 source_term(ic::IonizationChannel, args...) =
@@ -54,11 +67,17 @@ IonizationChannel(E, args...) =
 Base.show(io::IO, ic::IonizationChannel) =
     write(io, "IonizationChannel: Iₚ = $(ic.E) Ha = $(27.211ic.E) eV")
 
+function Base.show(io::IO, mime::MIME"text/plain", ic::IonizationChannel)
+    show(io, ic)
+    print(io, "\nSource term: ")
+    show(io, mime, ic.st)
+end
+
 # * Couplings
 
 abstract type AbstractCanonicalMomentumConservation end
-struct CanonicalMomentumConservation <:  AbstractCanonicalMomentumConservation end
-struct NoCanonicalMomentumConservation <:  AbstractCanonicalMomentumConservation end
+struct CanonicalMomentumConservation <: AbstractCanonicalMomentumConservation end
+struct NoCanonicalMomentumConservation <: AbstractCanonicalMomentumConservation end
 
 abstract type AbstractCoupling end
 Base.iszero(::AbstractCoupling) = false
@@ -68,6 +87,9 @@ canonical_momentum_conservation(::Type{<:AbstractCoupling}) = NoCanonicalMomentu
 struct NoCoupling <: AbstractCoupling end
 Base.iszero(::NoCoupling) = true
 Base.zero(::AbstractCoupling) = NoCoupling()
+Base.zero(::Type{<:AbstractCoupling}) = NoCoupling()
+
+Base.show(io::IO, ::NoCoupling) = write(io, "0")
 
 # ** Dipole couplings
 
@@ -79,6 +101,8 @@ end
 canonical_momentum_conservation(::DipoleCoupling) = CanonicalMomentumConservation()
 canonical_momentum_conservation(::Type{<:DipoleCoupling}) = CanonicalMomentumConservation()
 
+Base.show(io::IO, ::DipoleCoupling) = write(io, "𝐅⋅𝐫")
+
 # ** Coulomb couplings
 
 momentum_transfer(k::Number, p::Number) = k-p
@@ -90,3 +114,14 @@ struct CoulombCoupling{Coupling} <: AbstractCoupling
 end
 
 (cc::CoulombCoupling)(𝐤, 𝐩, _) = cc.coupling(𝐤, 𝐩)
+
+Base.show(io::IO, ::CoulombCoupling) = write(io, "ĝ")
+
+function Base.show(io::IO, mime::MIME"text/plain", g::CoulombCoupling)
+    write(io, "CoulombCoupling: ")
+    show(io, mime, g.coupling)
+end
+
+# * Exports
+
+export IonizationChannel