WIP: Envelope Plasma

examples/simple_interface/PlasmaSSFBS.jl

@@ -0,0 +1,19 @@
+using Luna
+
+radius = 9e-6 # HCF core radius
+flength = 25e-2 # HCF length
+
+gas = :Ar
+pressure = 5.0 # gas pressure in bar
+
+λ0 = 1500e-9 # central wavelength of the pump pulse
+τfwhm = 30e-15 # FWHM duration of the pump pulse
+energy = 1.7e-6 # energy in the pump pulse
+
+pssfbs = prop_capillary(radius, flength, gas, pressure; λ0, τfwhm, energy,
+                        trange=1.6e-12, λlims=(300e-9, 1.8e-6), loss=false, envelope=true, plasma=true)
+
+Plotting.prop_2D(pssfbs, :λ; modes=:sum, trange=(-700e-15, 100e-15), λrange=(300e-9, 1800e-9),
+                 dBmin=-30)
+Plotting.time_1D(pssfbs, range(0,flength,length=5), modes=:sum, trange=(-700e-15, 100e-15))
+Plotting.spec_1D(pssfbs, range(0,flength,length=5), modes=:sum, λrange=(300e-9, 1800e-9))

src/Interface.jl

@@ -509,8 +509,20 @@ function makeplasma!(out, grid, gas, plasma::Symbol, pol)
     push!(out, Nonlinear.PlasmaCumtrapz(grid.to, Et, ionrate, ionpot))
 end
 
+function makeplasma!(out, grid::Grid.EnvGrid, gas, plasma::Symbol, pol)
+    ionpot = PhysData.ionisation_potential(gas)
+    if plasma == :ADK
+        error("ADK not supported for envelopes at present")
+    elseif plasma == :PPT
+        ionrate = Ionisation.ionrate_fun!_PPTcached(gas, grid.referenceλ, rcycle=false)
+    else
+        throw(DomainError(plasma, "Unknown ionisation rate $plasma."))
+    end
+    Et = pol ? Array{Complex{Float64}}(undef, length(grid.to), 2) : Vector{Complex{Float64}}(undef, length(grid.to))
+    push!(out, Nonlinear.PlasmaFourier(grid.ωo, Et, ionrate, ionpot, grid.to[2] - grid.to[1]))
+end
+
 function makeresponse(grid::Grid.EnvGrid, gas, raman, kerr, plasma, thg, pol)
-    plasma && error("Plasma response for envelope fields has not been implemented yet.")
     isnothing(thg) && (thg = false) 
     out = Any[]
     if kerr
@@ -523,6 +535,7 @@ function makeresponse(grid::Grid.EnvGrid, gas, raman, kerr, plasma, thg, pol)
         end
     end
     raman && push!(out, Nonlinear.RamanPolarEnv(grid.to, Raman.raman_response(grid.to, gas)))
+    makeplasma!(out, grid, gas, plasma, pol)
     Tuple(out)
 end
 

src/Nonlinear.jl

@@ -83,9 +83,11 @@ function Kerr_env_thg(γ3, ω0, t)
     end
 end
 
+abstract type PlasmaPolar end
+
 "Response type for cumtrapz-based plasma polarisation, adapted from:
 M. Geissler, G. Tempea, A. Scrinzi, M. Schnürer, F. Krausz, and T. Brabec, Physical Review Letters 83, 2930 (1999)."
-struct PlasmaCumtrapz{R, EType, tType}
+struct PlasmaCumtrapz{R, EType, tType} <: PlasmaPolar
     ratefunc::R # the ionization rate function
     ionpot::Float64 # the ionization potential (for calculation of ionization loss)
     rate::tType # buffer to hold the rate
@@ -113,7 +115,7 @@ function PlasmaCumtrapz(t, E, ratefunc, ionpot)
 end
 
 "The plasma response for a scalar electric field"
-function PlasmaScalar!(Plas::PlasmaCumtrapz, E)
+function PlasmaScalar!(Plas::PlasmaCumtrapz, E::Vector{Float64})
     Plas.ratefunc(Plas.rate, E)
     Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt)
     @. Plas.fraction = 1-exp(-Plas.fraction)
@@ -155,8 +157,110 @@ function PlasmaVector!(Plas::PlasmaCumtrapz, E)
     Maths.cumtrapz!(Plas.P, Plas.J, Plas.δt)
 end
 
+struct PlasmaFourier{R, EType, tType, FTType} <: PlasmaPolar
+    ratefunc::R # the ionization rate function
+    ionpot::Float64 # the ionization potential (for calculation of ionization loss)
+    rate::tType # buffer to hold the rate
+    fraction::tType # buffer to hold the ionization fraction
+    phase::EType # buffer to hold the plasma induced (mostly) phase modulation
+    J::EType # buffer to hold the plasma current
+    P::EType # buffer to hold the plasma polarisation
+    ω::tType # buffer to hold the frequency grid
+    FT::FTType # Fourier transform to use for integrations
+    δt::Float64
+end
+
+"""
+    make_fft(E::Vector)
+
+Plan a suitable FFT for the electric field type `E`.
+
+"""
+function make_fft(E::Array{T,N}) where T<:Real where N
+    Utils.loadFFTwisdom()
+    FT = FFTW.plan_rfft(E, 1, flags=Luna.settings["fftw_flag"])
+    inv(FT)
+    Utils.saveFFTwisdom()
+    FT
+end
+
+function make_fft(E::Array{Complex{T},N}) where T<:Real where N
+    Utils.loadFFTwisdom()
+    FT = FFTW.plan_fft(E, 1, flags=Luna.settings["fftw_flag"])
+    inv(FT)
+    Utils.saveFFTwisdom()
+    FT
+end
+
+function PlasmaFourier(ω, E::Array{T,N}, ratefunc, ionpot, δt) where T<:Real where N
+    rate = similar(E)
+    fraction = similar(E)
+    phase = similar(E)
+    J = similar(E)
+    P = similar(E)
+    FT = make_fft(E)
+    PlasmaFourier(ratefunc, ionpot, rate, fraction, phase, J, P, ω, FT, δt)
+end
+
+function PlasmaFourier(ω, E::Array{Complex{T},N}, ratefunc, ionpot, δt) where T<:Real where N
+    rate = similar(ω)
+    fraction = similar(ω)
+    phase = similar(E)
+    J = similar(E)
+    P = similar(E)
+    FT = make_fft(E)
+    PlasmaFourier(ratefunc, ionpot, rate, fraction, phase, J, P, ω, FT, δt)
+end
+
+"The plasma response for a scalar electric field"
+function PlasmaScalar!(Plas::PlasmaFourier, E::Vector{Complex{Float64}})
+    # The equations used here for envelopes are very similar to the real field cases above.
+    # So far I have not been able to derive them directly. The modified loss term (Jabs)
+    # is taken from [1]. The plasma term is obtained by inspection. Essentially, you can start
+    # from the definition of the plasma refractive index, insert it into the general relation
+    # between polarisation and refractive index, and obtain this version (in the frequency domain)
+    # but with the electron density in the time domain(!). I also think this term could be obtained
+    # from Geissler's model [2], by making use of the Bedrosian identity for Hilbert transforms of
+    # products of functions, but I have not yet achieved it.
+    Plas.ratefunc(Plas.rate, E)
+    Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt)
+    @. Plas.fraction = 1-exp(-Plas.fraction)
+    @. Plas.phase = Plas.fraction * e_ratio * E
+    Plas.J .= Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E
+    # here we perform cumulative integration via the Fourier transform.
+    Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Plas.J) ./ Plas.ω)
+end
+
+function PlasmaVector!(Plas::PlasmaFourier, E::Array{Complex{Float64},2})
+    Ex = E[:,1]
+    Ey = E[:,2]
+    Em = @. hypot.(Ex, Ey)
+    Plas.ratefunc(Plas.rate, Em)
+    Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt)
+    @. Plas.fraction = 1-exp(-Plas.fraction)
+    @. Plas.phase = Plas.fraction * e_ratio * E
+    Plas.J .= Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ Em.^2 .* E
+    Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Plas.J) ./ Plas.ω)
+end
+
+"The plasma response for a scalar electric field"
+function PlasmaScalar!(Plas::PlasmaFourier, E::Vector{Float64})
+    Plas.ratefunc(Plas.rate, E)
+    Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt)
+    @. Plas.fraction = 1-exp(-Plas.fraction)
+    @. Plas.phase = Plas.fraction * e_ratio * E
+    fill!(Plas.J, 0.0)
+    for ii in eachindex(E)
+        if abs(E[ii]) > 0
+            Plas.J[ii] += Plas.ionpot * Plas.rate[ii] * (1-Plas.fraction[ii])/E[ii]
+        end
+    end
+    # here we perform cumulative integration via the Fourier transform.
+    Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Plas.J) ./ Plas.ω)
+end
+
 "Handle plasma polarisation routing to `PlasmaVector` or `PlasmaScalar`."
-function (Plas::PlasmaCumtrapz)(out, Et, ρ)
+function (Plas::PlasmaPolar)(out, Et, ρ)
     if ndims(Et) > 1
         if size(Et, 2) == 1 # handle scalar case but within modal simulation
             PlasmaScalar!(Plas, reshape(Et, size(Et,1)))
@@ -217,11 +321,8 @@ harmonic generation component of the response.
 """
 function RamanPolarField(t, r; thg=true)
     h = zeros(length(t)*2) # note double grid size, see explanation below
+    FT = make_fft(h)
     ht = view(h, 1:length(t))
-    Utils.loadFFTwisdom()
-    FT = FFTW.plan_rfft(h, 1, flags=Luna.settings["fftw_flag"])
-    inv(FT)
-    Utils.saveFFTwisdom()
     hω = FT * h
     Eω2 = similar(hω)
     Pω = similar(hω)
@@ -242,11 +343,8 @@ using response function `r`.
 """
 function RamanPolarEnv(t, r)
     h = zeros(length(t)*2) # note double grid size, see explanation below
+    FT = make_fft(h)
     ht = view(h, 1:length(t))
-    Utils.loadFFTwisdom()
-    FT = FFTW.plan_fft(h, 1, flags=Luna.settings["fftw_flag"])
-    inv(FT)
-    Utils.saveFFTwisdom()
     hω = FT * h
     Eω2 = similar(hω)
     Pω = similar(hω)

src/Stats.jl

@@ -216,7 +216,7 @@ If oversampling > 1, the field is oversampled before the calculation
     Oversampling can lead to a significant performance hit
 """
 function electrondensity(grid::Grid.RealGrid, ionrate!, dfun, aeff; oversampling=1)
-    to, Eto = Maths.oversample(grid.t, complex(grid.t), factor=oversampling)
+    to, Eto = Maths.oversample(grid.t, grid.t, factor=oversampling)
     δt = to[2] - to[1]
     # ionfrac! stores the time-dependent ionisation fraction in out and returns the max
     # ionisation rate
@@ -238,6 +238,29 @@ function electrondensity(grid::Grid.RealGrid, ionrate!, dfun, aeff; oversampling
     end
 end
 
+function electrondensity(grid::Grid.EnvGrid, ionrate!, dfun, aeff; oversampling=1)
+    to, Eto = Maths.oversample(grid.t, complex(grid.t), factor=oversampling)
+    δt = to[2] - to[1]
+    # ionfrac! stores the time-dependent ionisation fraction in out and returns the max
+    # ionisation rate
+    function ionfrac!(out, Et)
+        ionrate!(out, Et)
+        ratemax = maximum(out)
+        Maths.cumtrapz!(out, δt) # in-place cumulative integration
+        @. out = 1 - exp(-out)
+        return ratemax
+    end
+    frac = similar(to)
+    function addstat!(d, Eω, Et, z, dz)
+        # note: oversampling returns its arguments without any work done if factor==1
+        to, Eto = Maths.oversample(grid.t, Et, factor=oversampling)
+        @. Eto /= sqrt(ε_0*c*aeff(z)/2)
+        ratemax = ionfrac!(frac, abs.(Eto))
+        d["electrondensity"] = frac[end]*dfun(z)
+        d["peak_ionisation_rate"] = ratemax
+    end
+end
+
 """
     electrondensity(grid, ionrate, dfun, modes; oversampling=1)