Cycle average for ADK rate

src/Ionisation.jl

@@ -15,9 +15,10 @@ import Luna: @hlock
 
 Return a closure `ionrate!(out, E)` which calculates the ADK ionisation rate for the electric
 field `E` and places the result in `out`. If `threshold` is true, use [`ADK_threshold`](@ref)
-to avoid calculation below floating-point precision.
+to avoid calculation below floating-point precision. If `cycle_average` is `true`, calculate
+the cycle-averaged ADK ionisation rate instead.
 """
-function ionrate_fun!_ADK(ionpot::Float64, threshold=true)
+function ionrate_fun!_ADK(ionpot::Float64, threshold=true; cycle_average=false)
     nstar = sqrt(0.5/(ionpot/au_energy))
     cn_sq = 2^(2*nstar)/(nstar*gamma(nstar+1)*gamma(nstar))
     ω_p = ionpot/ħ
@@ -29,14 +30,27 @@ function ionrate_fun!_ADK(ionpot::Float64, threshold=true)
         thr = 0
     end
 
+    # Zenghu Chang: Fundamentals of Attosecond Optics (2011) p. 184
+    # Section 4.2.3.1 Cycle-Averaged Rate
+    # ̄w_ADK(Fₐ) = √(3/π) √(Fₐ/F₀) w_ADK(Fₐ) where Fₐ is the field amplitude 
+    Ip_au = ionpot / au_energy
+    F0_au = (2Ip_au)^(3/2)
+    F0 = F0_au*au_Efield
+    avfac = sqrt.(3/(π*F0))
+
+
     ionrate! = let nstar=nstar, cn_sq=cn_sq, ω_p=ω_p, ω_t_prefac=ω_t_prefac, thr=thr
         function ir(E)
             if abs(E) >= thr
-                (ω_p*cn_sq*
-                (4*ω_p/(ω_t_prefac*abs(E)))^(2*nstar-1)
-                *exp(-4/3*ω_p/(ω_t_prefac*abs(E))))
+                r = (ω_p*cn_sq*
+                    (4*ω_p/(ω_t_prefac*abs(E)))^(2*nstar-1)
+                    *exp(-4/3*ω_p/(ω_t_prefac*abs(E))))
+                if cycle_average
+                    r *= avfac*sqrt(abs(E))
+                end
+                return r
             else
-                zero(E)
+                return zero(E)
             end
         end
         function ionrate!(out, E)
@@ -119,9 +133,9 @@ function ionrate_fun!_PPTcached(material::Symbol, λ0; kwargs...)
 end
 
 function ionrate_fun!_PPTcached(ionpot::Float64, λ0, Z, l;
-                                sum_tol=1e-4, rcycle=true, N=2^16, Emax=nothing,
+                                sum_tol=1e-4, cycle_average=false, N=2^16, Emax=nothing,
                                 cachedir=joinpath(homedir(), ".luna", "pptcache"))
-    h = hash((ionpot, λ0, Z, l, sum_tol, rcycle, N, Emax))
+    h = hash((ionpot, λ0, Z, l, sum_tol, cycle_average, N, Emax))
     fname = string(h, base=16)*".h5"
     fpath = joinpath(cachedir, fname)
     lockpath = joinpath(cachedir, "pptlock")
@@ -134,7 +148,7 @@ function ionrate_fun!_PPTcached(ionpot::Float64, λ0, Z, l;
         return rate
     else
         E, rate = makePPTcache(ionpot::Float64, λ0, Z, l;
-                               sum_tol=sum_tol, rcycle=rcycle, N=N, Emax=Emax)
+                               sum_tol=sum_tol, cycle_average, N=N, Emax=Emax)
         @info "Saving PPT rate cache for $(ionpot/electron) eV, $(λ0*1e9) nm in $cachedir"
         pidlock = mkpidlock(lockpath)
         if isfile(fpath) # makePPTcache takes a while - has another process saved first?
@@ -160,7 +174,7 @@ function loadPPTaccel(fpath)
 end
 
 function makePPTcache(ionpot::Float64, λ0, Z, l;
-                      sum_tol=1e-4, rcycle=true, N=2^16, Emax=nothing)
+                      sum_tol=1e-4, cycle_average=false, N=2^16, Emax=nothing)
     Emax = isnothing(Emax) ? 2*barrier_suppression(ionpot, Z) : Emax
 
     # ω0 = 2π*c/λ0
@@ -169,7 +183,7 @@ function makePPTcache(ionpot::Float64, λ0, Z, l;
 
     E = collect(range(Emin, stop=Emax, length=N));
     @info "Pre-calculating PPT rate for $(ionpot/electron) eV, $(λ0*1e9) nm"
-    rate = ionrate_PPT(ionpot, λ0, Z, l, E; sum_tol=sum_tol, rcycle=rcycle);
+    rate = ionrate_PPT(ionpot, λ0, Z, l, E; sum_tol=sum_tol, cycle_average);
     @info "PPT pre-calcuation done"
     return E, rate
 end
@@ -235,13 +249,13 @@ function ionrate_fun!_PPT(args...)
 end
 
 """
-    ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, rcycle=true)
+    ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, cycle_average=false)
 
 Create closure to calculate PPT ionisation rate.
 
 # Keyword arguments
 - `sum_tol::Number`: Relative tolerance used to truncate the infinite sum.
-- `rcycle::Bool`: If `true` (default), remove the cycle-average factor.
+- `cycle_average::Bool`: If `false` (default), calculate the cycle-averaged rate
 
 # References
 [1] Ilkov, F. A., Decker, J. E. & Chin, S. L.
@@ -254,7 +268,7 @@ Ultrashort filaments of light in weakly ionized, optically transparent
 media. Rep. Prog. Phys. 70, 1633–1713 (2007)
 (Appendix A)
 """
-function ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, rcycle=true)
+function ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, cycle_average=false)
     Ip_au = ionpot / au_energy
     ns = Z/sqrt(2Ip_au)
     ls = ns-1
@@ -294,7 +308,7 @@ function ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, rcycle=true)
                 lret *= exp(-2v*(asinh(γ) - γ*sqrt(1+γ2)/(1+2γ2)))
                 lret *= Ip_au * γ2/(1+γ2)
                 # Remove cycle average factor, see eq. (2) of [1]
-                if rcycle
+                if !cycle_average
                     lret *= sqrt(π*E0_au/(3E_au))
                 end
                 k = ceil(v)

test/test_ionisation.jl

@@ -1,5 +1,6 @@
 import Test: @test, @testset
-import Luna: Ionisation, Maths, PhysData
+import Luna: Ionisation, Maths, PhysData, Tools
+import NumericalIntegration: integrate, SimpsonEven
 
 @test Ionisation.ionrate_ADK(:He, 1e10) ≈ 1.2416371415312408e-18
 @test Ionisation.ionrate_ADK(:He, 2e10) ≈ 1.0772390893742478
@@ -44,6 +45,36 @@ outneg = similar(out)
 ratefun!(outneg, -E)
 @test out == outneg
 
+##
+@testset "cycle-averaging $gas" for gas in (:He, :Ar, :Xe)
+
+    bsi = Tools.field_to_intensity(
+        Ionisation.barrier_suppression(
+            PhysData.ionisation_potential(gas),
+            1))
+    isy = 10 .^ collect(16:0.01:log10(bsi))
+    E0 = Tools.intensity_to_field.(isy)
+
+    Ip_au = PhysData.ionisation_potential(gas; unit=:atomic)
+    F0_au = (2Ip_au)^(3/2)
+    F0 = F0_au*PhysData.au_Efield
+
+    t = collect(range(0, 2π; length=2^12))
+    Et = sin.(t)
+
+    adk_avg = zero(E0)
+    rf! = Ionisation.ionrate_fun!_ADK(gas)
+    out = zero(Et)
+    for (idx, E0i) in enumerate(E0)
+        rf!(out, E0i*Et)
+        adk_avg[idx] = 1/2π * integrate(t, out, SimpsonEven())
+    end
+
+    adk = Ionisation.ionrate_ADK.(gas, E0)
+    adk_avg_kw = Ionisation.ionrate_ADK.(gas, E0; cycle_average=true)
+    @test all(isapprox.(adk_avg, adk_avg_kw; rtol=0.05))
+end
+##
 # this_folder = dirname(@__FILE__)
 # import CSV
 # import PyPlot: plt, pygui