diff --git a/src/Ionisation.jl b/src/Ionisation.jl index 1e483933..e4d2cb78 100644 --- a/src/Ionisation.jl +++ b/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) diff --git a/test/test_ionisation.jl b/test/test_ionisation.jl index a28e4ace..e0ca9f6c 100644 --- a/test/test_ionisation.jl +++ b/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