Implementation of an additional ionization rate function for bulk

[Gvaz98]

Hello,

I am trying to add the Keldysh ionization rate to include plasma effects for bulk samples. I have already written the equation in Ionization.jl copying the structure of the ADK rate https://github.com/Gvaz98/Luna.jl/blob/08338aff8473233f8de79af3c2d12074923a1cd0/src/Ionisation.jl#L389) , but have a few doubts:

1. The equations include a typical reduced electron-hole mass. But from what I see from the other rate equations, the variables are normalized to atomic energies. So in this case it needs to be divided it by PhysData.m_u? I ask because the ADK rate uses the SI mass of an electron.
2. The equation includes a summation of a function, f, for now I write it as sum(f, UnitRange(0, Nsum)), where Nsum is an optional parameter with default 1e3. Is there a better/faster option?

Additionally, for a first try, I changed one of the free space examples to include plasma effects https://github.com/Gvaz98/Luna.jl/blob/b9fb68ee9c13e49b0daf5a809aadd7cfcdeea975/examples/low_level_interface/freespace/radial_hcf_keldysh_plasma.jl However, when propagation starts the rate function receives NaN for the electric field, not sure why.

Thanks in advance

[chrisbrahms]

Hi @Gvaz98, thanks a lot for putting this together! Unfortunately I don't have a huge amount of time right now as I'm travelling. I will be able to look at the error in more detail in August, but just a quick thought--it may be easier to get it to work with the other freespace example, as the one you modified actually simulates HCF propagation, but using the Hankel transform (this is formally different but accesses the same physics).

Regarding your two questions:

1. The ADK equation we use is written in SI units as a whole including all the constants. This is not the case for PPT and your direct translation of the Keldysh equations, so yes you should convert to atomic units (and then you need to convert the rate to SI units at the end).
2. You can have a look at the converge_series function in the Maths module: https://github.com/LupoLab/Luna.jl/blob/675f0362b316585d114d11fa03519de59ce11797/src/Maths.jl#L716 Rather than defining the number of steps, it defines a tolerance. If the Keldysh rate turns out to be slow, you can also use the same pre-calculation approach we use for the PPT rate (though this would very much be a second step).

[chrisbrahms]

@jtravs if you have a moment you could also look at this before I'm back?

[Gvaz98]

Hello again,

Small update. Managed to correct the NaN error. Currently, I am testing with the following example. Still radial to reduce the debug time. https://github.com/Gvaz98/Luna.jl/blob/changes/examples/low_level_interface/freespace/radial_env_keldysh_plasma.jl

Now I am having issues with an “InexactError” e.g. ERROR: LoadError: InexactError: Float64(-7.86401876814954e-20 - 1.6661659225482432e-18im) It derives from the calculation of the phase in https://github.com/Gvaz98/Luna.jl/blob/bc5da7bbb0d1f0e54c663494f44d9abcc97f5d9e/src/Nonlinear.jl#L146

@. Plas.phase = Plas.fraction * e_ratio * E

In particular because the phase is Float64 and the result on the right-hand side is a ComplexF64 because of E. But I found some strange behaviours:

• When I run with PPT or ADK the type of E is a SubArray{ComplexF64, 1, Matrix{ComplexF64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true} but with Keldysh is a 256-element view(::Matrix{ComplexF64}, :, 1) with eltype ComplexF64. PPT and ADk also give the same error at higher energies.
• I ran plasma_ssfbs_modeAvg.jl and the type of E is a simple Vector{Float64}.

Here is the link to what I added so far. https://github.com/LupoLab/Luna.jl/compare/master...Gvaz98:Luna.jl:changes#

EDIT: From running other examples, I assume the phase should also be a complex, correct? For now I am using the following as a plasma response to assure the phase is complex.:

Nonlinear.PlasmaCumtrapz(grid.to, convert(Array{ComplexF64}, grid.to), ionrate, ionpot)

Probably not the most elegant, but at least the code runs. Still need to benchmark it specially because I am still using only SI units in Keldysh. I don't see a reference to atomic units in the paper.

[chrisbrahms]

Hi!

Sorry for the slow response. I've been away for several weeks but I'm back now.

The main problem you're running into here is that the PlasmaCumtrapz model for the plasma only works for real fields. In principle you could adopt #241 to extend it to envelope fields, but as you can from the discussion there, it's not fully worked through yet. The easiest thing is probably just to switch to real fields instead.

[Gvaz98]

Hello,

Sorry, I have not had much time and could not benchmark yet.

Just a few questions. How much does changing from envelope to real field change the physics? Or is just the equations? Just to be sure, could you verify if the equations that I am using from 10.1016/j.physrep.2006.12.005 can be done fully in SI? The equations include constants that should be 1 in atomic units and would normally be omitted. Plus, I tried atomic resulting rates several orders lower than PPT.

Meanwhile, I think I managed to make a cached version of the rate function. Although I am not sure in how to define the Emax and Emin for it. For now, I am using:

function makeKeldyshcache(ionpot::Float64, λ0;
                        rtol = 1e-6, maxiter = 10000, N=2^16, Emax=nothing)
    Imax=ionpot/λ0^3*c #Intensity in W/m^2 where the energy matches the ionisation potential, the duration λ/c and area of λ^2.
    #i.e. a λ^3 laser
    Emax = isnothing(Emax) ? 1e4*sqrt(μ_0*c*Imax) : Emax
    Emin = Emax/50000

    E = collect(range(Emin, stop=Emax, length=N));
    @info @sprintf("Pre-calculating Keldysh rate rate for %.2f eV, %.1f nm...", ionpot/electron, 1e9λ0)
    rate = ionrate_Keldysh(ionpot, λ0, E; rtol = rtol, maxiter = maxiter);
    @info "...Keldysh pre-calculation done"
    return E, rate
end

Where the “safety factor” 1e4 is arbitrary.