diff --git a/src/PhysData.jl b/src/PhysData.jl index f05372f4..12497f95 100644 --- a/src/PhysData.jl +++ b/src/PhysData.jl @@ -47,15 +47,19 @@ const m_u = ustrip(CODATA2014.m_u) "Atomic unit of electric polarisability" const au_polarisability = electron^2*ustrip(CODATA2014.a_0)^2/au_energy -const gas = (:Air, :He, :HeJ, :HeB, :Ne, :Ar, :ArB, :Kr, :Xe, :N2, :H2, :O2, :CH4, :SF6, :N2O, :D2) +const gas = (:Air, :He, :HeJ, :HeB, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrB, :KrBideauMehu, :Xe, :N2, :H2, :O2, :CH4, :SF6, :N2O, :D2) const gas_str = Dict( :He => "He", :HeB => "He", :HeJ => "He", :Ar => "Ar", :ArB => "Ar", + :ArBideauMehu => "Ar", :Ne => "Neon", + :NeBideauMehu => "Neon", :Kr => "Krypton", + :KrB => "Krypton", + :KrBideauMehu => "Krypton", :Xe => "Xenon", :Air => "Air", :N2 => "Nitrogen", @@ -138,6 +142,24 @@ function γ_QuanfuHe(A, B, C, dens) return μm -> complex((1 + 1e-8*(A + B/(C - (1e4/μm)^2))))/dens end +""" + γ_BideauMehu2(B1, C1, B2, C2, dens) + +2-term Sellmeier expression for Ne from https://doi.org/10.1016/0022-4073(81)90057-1 +""" +function γ_BideauMehu2(B1, C1, B2, C2, dens) + return μm -> (((B1/(C1-1/μm^2) + B2/(C2-1/μm^2))+1)^2 - 1)/dens +end + +""" + γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) + +3-term Sellmeier expression for Ar and Kr from https://doi.org/10.1016/0022-4073(81)90057-1 +""" +function γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) + return μm -> (((B1/(C1-1/μm^2) + B2/(C2-1/μm^2) + B3/(C3-1/μm^2))+1)^2 - 1)/dens +end + """ sellmeier_gas(material::Symbol) @@ -166,6 +188,12 @@ function sellmeier_gas(material::Symbol) B2 = 4018.63e-8 C2 = 5.728e-3 return γ_Börzsönyi(B1/dens, C1, B2/dens, C2) + elseif material == :NeBideauMehu + B1 = 0.00128145 + C1 = 184.6661 + B2 = 0.0220486 + C2 = 376.84 + return γ_BideauMehu2(B1, C1, B2, C2, dens) elseif material == :Ar B1 = 0.00032323117217767093 C1 = 0.0045416501944977915 @@ -180,12 +208,36 @@ function sellmeier_gas(material::Symbol) B2 = 34458.31e-8 C2 = 8.066e-3 return γ_Börzsönyi(B1/dens, C1, B2/dens, C2) + elseif material == :ArBideauMehu + B1 = 2.50141e-3 + C1 = 91.012 + B2 = 5.00283e-4 + C2 = 87.892 + B3 = 5.22343e-2 + C3 = 214.02 + return γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) elseif material == :Kr + B1 = 7.156207185088921e-6 + C1 = 0.016126479267946525 + B2 = 0.00023020614937540427 + C2 = 0.01348822574907029 + B3 = 0.000604323109684193 + C3 = 0.004204161467601559 + return γ_JCT(B1/dens, C1, B2/dens, C2, B3/dens, C3) + elseif material == :KrB B1 = 26102.88e-8 C1 = 2.01e-6 B2 = 56946.82e-8 C2 = 10.043e-3 return γ_Börzsönyi(B1/dens, C1, B2/dens, C2) + elseif material == :KrBideauMehu + B1 = 0.00253637 + C1 = 65.4742 + B2 = 0.00273649 + C2 = 73.698 + B3 = 0.0620802 + C3 = 181.08 + return γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) elseif material == :Xe B1 = 103701.61e-8 C1 = 12750e-6 @@ -588,7 +640,7 @@ References: function γ3_gas(material::Symbol; source=nothing) # TODO: More Bishop/Shelton; Wahlstrand updated values. if source === nothing - if material in (:He, :HeB, :HeJ, :Ne, :Ar, :ArB, :Kr, :Xe, :N2) + if material in (:He, :HeB, :HeJ, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrB, :KrBideauMehu,:Xe, :N2) source = :Lehmeier elseif material in (:H2, :CH4, :SF6, :D2) source = :Shelton @@ -605,11 +657,11 @@ function γ3_gas(material::Symbol; source=nothing) # Table 1 in [3] if material in (:He, :HeB, :HeJ) fac = 1 - elseif material == :Ne + elseif material in (:Ne, :NeBideauMehu) fac = 1.8 - elseif material in (:Ar, :ArB) + elseif material in (:Ar, :ArB, :ArBideauMehu) fac = 23.5 - elseif material == :Kr + elseif material in (:Kr, :KrB, :KrBideauMehu) fac = 64.0 elseif material == :Xe fac = 188.2 @@ -731,11 +783,11 @@ Possible units are `:SI`, `:atomic` and `:eV`. function ionisation_potential(material; unit=:SI) if material in (:He, :HeB, :HeJ) Ip = 0.9036 - elseif material == :Ne + elseif material in (:Ne, :NeBideauMehu) Ip = 0.7925 - elseif material in (:Ar, :ArB) + elseif material in (:Ar, :ArB, :ArBideauMehu) Ip = 0.5792 - elseif material == :Kr + elseif material in (:Kr, :KrB, :KrBideauMehu) Ip = 0.5142 elseif material == :Xe Ip = 0.4458 @@ -777,11 +829,11 @@ Return the quantum numbers of the `material` for use in the PPT ionisation rate. """ function quantum_numbers(material) # Returns n, l, ion Z - if material in (:Ar, :ArB) + if material in (:Ar, :ArB, :ArBideauMehu) return 3, 1, 1 - elseif material == :Ne + elseif material in (:Ne, :NeBideauMehu) return 2, 1, 1; - elseif material == :Kr + elseif material in (:Kr, :KrB, :KrBideauMehu) return 4, 1, 1 elseif material == :Xe return 5, 1, 1 @@ -822,9 +874,9 @@ function polarisability(material, ion=false; unit=:SI) end if material in (:He, :HeB, :HeJ) return (ion ? 0.2811 : 1.3207)*factor - elseif material == :Ne + elseif material in (:Ne, :NeBideauMehu) return (ion ? 1.2417 : 2.376)*factor - elseif material in (:Ar, :ArB) + elseif material in (:Ar, :ArB, :ArBideauMehu) return (ion ? 6.807 : 10.762)*factor else return missing @@ -861,11 +913,11 @@ Ammosov, M. V., Delone, N. B. & Krainov, V. P. Tunnel Ionization Of Complex Atom function Cnl_ADK(material) if material in (:He, :HeB, :HeJ) return 1.99 - elseif material == :Ne + elseif material in (:Ne, :NeBideauMehu) return 1.31 - elseif material in (:Ar, :ArB) + elseif material in (:Ar, :ArB, :ArBideauMehu) return 1.9 - elseif material == :Kr + elseif material in (:Kr, :KrB, :KrBideauMehu) return 2.17 elseif material == :Xe return 2.27 @@ -889,6 +941,17 @@ function lookup_glass(material::Symbol) return spl end + +function Ith(material; ionLevel=1) + + Ip = ionisation_potential(material; unit=:eV) + Ith = 3.8e9*(Ip^4)/(ionLevel^2) + + return Ith*1e4 # convert from W/cm^2 to W/m^2 + +end + + """ lookup_metal(material::Symbol) diff --git a/src/Tools.jl b/src/Tools.jl index 04c2e9eb..9fa07afc 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -235,6 +235,296 @@ function pressureZDW(a::Number, gas::Symbol, λzd; Pmax=100, clad=:SiO2, kwargs. end end + + +""" + Calculates the critical electron density for a given wavelength, which is the density at which the plasma frequency equals the frequency of the light + Ref: https://doi.org/10.1016/j.physrep.2006.12.005, top of page 56 + Arguments: + λ: the wavelength of light +""" +function critical_electron_density(λ) + ω = PhysData.wlfreq(λ) + return (PhysData.ε_0*PhysData.m_e*ω^2)/PhysData.electron^2 +end + + +""" + Calculates the energy required for a certain soliton order + Arguments: + a: radius of the capillary core + gas: type of gas in the capillary + pressure: pressure of the gas in the capillary + τFWHM: pulse duration (FWHM) + λp: central wavelength of the pump (soliton) pulse + N: soliton order + energy_lims: limits for the root finding algorithm that looks for the energy; by default these are set to 1 nJ to 1 J +""" +function energyN(a, gas, pressure, τFWHM, λp, N; energy_lims=(1e-9, 1.0)) + + EforN = find_zero(energy_lims) do energy + params = Tools.capillary_params(energy, τFWHM, λp, a, gas; P=pressure) + params.N - N + end + + return EforN + +end + +""" + Helper function to calculate the linear contribution to the phase-mismatch Δβ, which is used in the phase-matching condition for RDW emission +""" +function Δβlin(a, gas, pressure, λp) + + mode = Capillary.MarcatiliMode(a, gas, pressure) + ωsol = PhysData.wlfreq(λp) + β1 = Modes.dispersion(mode, 1, ωsol) + β0 = Modes.β(mode, ωsol) + Δβlin(ω) = Modes.β(mode, ω) - β1*(ω - ωsol) - β0 + + return Δβlin + +end + +""" + Calculate the estimated compression factor and quality factor of the soliton self-compression + Schade et al. "Scaling rules for high quality soliton self-compression in hollow-core fibers" https://doi.org/10.1364/OE.426307 + Arguments: + N: soliton order + Returns compression factor and compression quality (Fc, Q) +""" +function compression_estimation(a, gas, pressure, λp, τFWHM, N; include_ξ=false) + + A(N) = (1/(2*N)+1.7/(N^2)) + Q(N) = 3.7/(N + 2.2) + + ωsol = PhysData.wlfreq(λp) + mode = Capillary.MarcatiliMode(a, gas, pressure) + β2 = Modes.dispersion(mode, 2, ωsol) + β3 = Modes.dispersion(mode, 3, ωsol) + ξ = β3/(τFWHM*abs(β2)) + + Fc = 3/(A(N)) + + if include_ξ + Fc *= (1-N*ξ) + end + + return Fc, Q(N) +end + +function Ppeak(τFWHM, energy; pulse_shape=:sech) + + if pulse_shape == :sech + return 0.88*(energy/τFWHM) + elseif pulse_shape == :gauss + return 0.94*(energy/τFWHM) + else + error("pulse_shape must be one of: :sech, :gauss") + end + +end + +""" + Helper function to calculate the nonlinear contribution to the phase-mismatch Δβ, which is used in the phase-matching condition for RDW emission +""" +function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=false, input_pulse_shape=:gauss) + + energy = energyN(a, gas, pressure, τFWHM, λp, soliton_order) + params = Tools.capillary_params(energy, τFWHM, λp, a, gas; P=pressure) + ωsol = PhysData.wlfreq(λp) + + if includeLoss # if include loss is true, decrease the energy for which the peak power is calculated + Lfiss = params.Lfiss + mode = Capillary.MarcatiliMode(a, gas, pressure) + energy = (1.0 - Modes.α(mode, PhysData.wlfreq(λp); z=Lfiss))*energy + end + + if input_pulse_shape == :sech + Pp = Ppeak(τFWHM, energy, pulse_shape=:sech) + elseif input_pulse_shape == :gauss + Pp = 0.936*Ppeak(τFWHM, energy, pulse_shape=:gauss) # 0.88/0.94 + else + error("input_pulse_shape must be one of: :sech, :gauss") + end + + soliton_factor = ((2*soliton_order-1)/soliton_order)^2 # ((2N-1)/N)^2 + + return params.γ*(soliton_factor*Pp)*(1/ωsol) + +end + +""" + Helper function to calculate the ionisation contribution to the phase-mismatch Δβ, which is used in the phase-matching condition for RDW emission +""" +function Δβion(a, gas, pressure, λp, ionisation_fraction) + neutral_density = PhysData.density(gas, pressure) + electron_density = ionisation_fraction*neutral_density + critical_density = critical_electron_density(λp) + ωsol = PhysData.wlfreq(λp) + mode = Capillary.MarcatiliMode(a, gas, pressure) + nlin = real.(Capillary.neff(mode, ωsol)) # taking the real part here, otherwise the ionisation term becomes complex + return (1/(2*PhysData.c*nlin))*(electron_density/critical_density)*(ωsol^2) +end + +""" + A helper type to store parameters for phase-matching +""" +struct PhaseMatching + mode + ωsol + β0 + β1 + ΔβnonlinCoeff + ΔβionCoeff +end + +""" + A helper function to calculate the phase-matching curve from pre-calculated parameters before passing to the root-finding algorithm; this speeds up the calculation a lot compared to when just passing Δβ as a full function, that has to be evaluated at every iteration of the root-finding +""" +function (pm::PhaseMatching)(ω) + + return Modes.β(pm.mode, ω) - pm.β1*(ω - pm.ωsol) - pm.β0 - pm.ΔβnonlinCoeff*ω + pm.ΔβionCoeff*(1/ω) + +end + +""" + A helper function to pre-compute parameters for phase-matching +""" +function make_PhaseMatching(a, gas, pressure, λp; + include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, includeLoss=false, input_pulse_shape=:gauss, + include_Δβion=false, ionisation_fraction=nothing) + + mode = Capillary.MarcatiliMode(a, gas, pressure) + ωsol = PhysData.wlfreq(λp) + β1 = Modes.dispersion(mode, 1, ωsol) + β0 = Modes.β(mode, ωsol) + + if include_Δβnonlin # perform the assertions only if the nonlinear contribution is included, otherwise these parameters are not needed and can be left as nothing + @assert !isnothing(soliton_order) "If include_Δβnonlin is true, soliton_order must be given" + @assert soliton_order > 0 "Soliton order must be positive" + @assert !isnothing(τFWHM) "If include_Δβnonlin is true, τFWHM must be given" + end + + if include_Δβion # perform the assertions only if the ionisation contribution is included, otherwise these parameters are not needed and can be left as nothing + @assert !isnothing(ionisation_fraction) "If include_Δβion is true, ionisation_fraction must be given" + @assert ionisation_fraction >= 0 && ionisation_fraction <= 1 "Ionisation fraction must be between 0 and 1" + end + + ΔβnonlinCoeff = include_Δβnonlin ? Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=includeLoss,input_pulse_shape=input_pulse_shape) : 0.0 + ΔβionCoeff = include_Δβion ? Δβion(a, gas, pressure, λp, ionisation_fraction) : 0.0 + + return PhaseMatching(mode, ωsol, β0, β1, ΔβnonlinCoeff, ΔβionCoeff) +end + +""" + Function for calculating the phase-matched resonant dispersive wave (RDW) wavelength in a gas-filled capillary, taking into account the nonlinear and ionisation contributions to phase-matching. + Ref: PRL 115, 033901 (2015), https://doi.org/10.1103/PhysRevLett.115.033901 + Arguments: + a: radius of the capillary core + gas: type of gas in the capillary + pressure: pressure of the gas in the capillary + λp: central wavelength of the pump (soliton) pulse + kwargs: additional keyword arguments for the mode (e.g. m, n, ϕ, kind, clad; used in calling MarcatiliMode() line 95 in Luna.Capillary) + soliton_order: used to calculate the energy/peak power + τFWHM: pulse duration (FWHM) + ionisation_fraction: electron fraction + λlims: limits for the root finding algorithm that looks for the phase-matching; by default these are set to the minimum wavelength just above the first resonance in the Sellmeier equation for the given gas and the maximum for the pump wavelength minus 1 nm +""" +function λRDWfull(a, gas, pressure, λp; + include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, includeLoss=false, input_pulse_shape=:gauss, + include_Δβion=false, ionisation_fraction=nothing, + λlims=nothing) + + if isnothing(λlims) + λUVlim = Dict( + :He => 90e-9, + :HeB => 90e-9, + :HeJ => 90e-9, + :Ne => 100e-9, + :NeBideauMehu => 100e-9, + :Ar => 110e-9, + :ArB => 110e-9, + :ArBideauMehu => 110e-9, + :Kr => 130e-9, + :KrB => 130e-9, + :KrBideauMehu => 130e-9 + ) + λlims = (λUVlim[gas], λp-10e-9) # pump wavelength minus 10 nm for safety + end + + phase_matching = make_PhaseMatching(a, gas, pressure, λp; + include_Δβnonlin=include_Δβnonlin, soliton_order=soliton_order, τFWHM=τFWHM, includeLoss=includeLoss, input_pulse_shape=input_pulse_shape, + include_Δβion=include_Δβion, ionisation_fraction=ionisation_fraction) + + # if ionization effects are included, shrink λlims, because the dispersion curve gets lifted in the IR and a new phase-matched point appears near the pump + if include_Δβion + + # first find the linear phase-matching wavelength + λlin = Tools.λRDW(a, gas, pressure, λp; λlims=λlims) # use the OG Luna function to calculate the linear phase-matched RDW + + # find the root close to the pump wavelength + ωmax = NaN + try + @assert !ismissing(λlin) "λlin should not be NaN here, the first root-finding has failed" + ωmax = find_zero(phase_matching, extrema(PhysData.wlfreq.((λlin+50e-9, λp-10e-9)))) # linear RDW + 50 nm for safety + λlims = (minimum(λlims), PhysData.wlfreq(ωmax) - 10e-9) + catch err + @warn "A phase-matching wavelength couldn't be found for Δβ..." exception=(err, catch_backtrace()) + end + + end + + # find RDW phase-matched wavelength with the full Δβ and potentially shrunk λlims + ωRDW = NaN + try + ωRDW = find_zero(phase_matching, extrema(PhysData.wlfreq.(λlims))) + catch err + @warn "A phase-matching wavelength couldn't be found for the given parameters..." exception=(err, catch_backtrace()) + end + + return PhysData.wlfreq(ωRDW) + +end + +function fλ(λ, material, pressure; temperature=PhysData.roomtemp) + + χ1function = PhysData.χ1_fun(material, pressure, temperature) + + return Maths.derivative(χ1function, λ, 2) + +end + + +function δ(λ, zdw, material, pressure; kind=:HE, n=1, m=1) + + unm = Capillary.get_unm(n, m, kind) + + return ((unm^2*λ^3)/(8*π^3*PhysData.c^2))*((fλ(λ, material, pressure))/(fλ(zdw, material, pressure)) - 1) + +end + +function NmaxSelfFocusing(λp, zdw, gas, τFWHM, a; safetyFactor=10, pulseShape=:gauss, kind=:HE, n=1, m=1) + + pressure = pressureZDW(a, gas, zdw) + τ0 = τfw_to_τ0(τFWHM, pulseShape) + + return sqrt((τ0^2*λp)/(safetyFactor*abs(δ(λp, zdw, gas, pressure; kind=kind, n=n, m=m)))) + +end + +function NmaxIonisation(λp, zdw, gas, τFWHM, a; safetyFactor=10, pulseShape=:gauss, kind=:HE, n=1, m=1) + + pressure = pressureZDW(a, gas, zdw) + τ0 = τfw_to_τ0(τFWHM, pulseShape) + ωpump = PhysData.wlfreq(λp) + _, _, n2 = getN0n0n2(ωpump, gas; P=pressure) + unm = Capillary.get_unm(n, m, kind) + + return sqrt((τ0^2*n2*PhysData.PhysData.Ith(gas)*unm^2)/(safetyFactor*π*λp*abs(δ(λp, zdw, gas, pressure; kind=kind, n=n, m=m))*fλ(zdw, gas, pressure))) + +end + function aperture_filter(a, dist, radius) w0 = 0.64*a function filter(λ) diff --git a/test/test_tools.jl b/test/test_tools.jl index 85e5aaf3..ec9db806 100644 --- a/test/test_tools.jl +++ b/test/test_tools.jl @@ -10,7 +10,7 @@ p = Tools.params(300e-6, 10e-15, 800e-9, m, :HeB, P=0.4) @test isapprox(p.zdw, 378.8e-9, rtol=1e-2) @test isapprox(p.P0/p.Pcr, 0.0398, rtol=2e-2) # test backup zdw method -p = Tools.capillary_params(6e-9, 20e-15, 800e-9, 14e-6, :Kr, P=15.0) +p = Tools.capillary_params(6e-9, 20e-15, 800e-9, 14e-6, :KrB, P=15.0) @test isapprox(p.zdw, 7.693023014958748e-7, rtol=1e-7) end