Updated Volkov phases to use "automatic integration"

Author: jagot

src/automatic_integration.jl

@@ -117,6 +117,8 @@ function AutomaticIntegration(f, x, x₀; k=3, init=zero(f(x₀)))
     ai
 end
 
+Base.eltype(ai::AutomaticIntegration) = eltype(ai.∫f)
+
 (ai::AutomaticIntegration)(x; i=find_closest(ai.x, x)) =
     ai.∫f[i] + primitive_integral(ai.f, ai.x[i], x, ai.c, ai.w)
 

src/dyson_expansions.jl

@@ -361,8 +361,9 @@ function evaluate_momenta!(𝐩s, prefactors, system, unique_momenta, indetermin
     for idm in indeterminate_momenta
         uidm = unique_momenta[idm]
         a,b = i[uidm[1]],i[uidm[2]]
-        set_momentum!(𝐩s, stationary_momentum(system.volkov, a, b), idm)
-        τ = system.t[a]-system.t[b]
+        ta,tb = system.t[a], system.t[b]
+        set_momentum!(𝐩s, stationary_momentum(system.volkov, ta, tb), idm)
+        τ = ta - tb
         ζ = (2π/(im*τ + ϵ))^(3/2)
         prefactors[idm] = ζ
     end
@@ -424,9 +425,10 @@ function integrate_diagram(::Type{Amp}, system::System, diagram::Diagram, iref,
         Sₑₗ = zero(ctT)
         for j = 1:order
             a,b = is[j], is[j+1]
-            τ = system.t[a] - system.t[b]
+            ta,tb = system.t[a], system.t[b]
+            τ = ta - tb
             Sᵢₒₙ += Eᵢₒₙₛ[j]*τ
-            Sₑₗ += volkov_phase(𝐩s[j], system.volkov, a, b)
+            Sₑₗ += volkov_phase(𝐩s[j], system.volkov, ta, tb)
         end
         S = Sᵢₒₙ + Sₑₗ
         aₚᵣₒₚ = prod(prefactors)*exp(-im*S)

src/volkov.jl

@@ -1,79 +1,49 @@
-struct VolkovPhases{Time,T,U}
-    t::Time
-    ∫A::Vector{T}
-    ∫A²::Vector{U}
+struct VolkovPhases{ReferenceTime,∫At,∫A²t}
+    T::ReferenceTime
+    ∫A::∫At
+    ∫A²::∫A²t
 end
 
-function VolkovPhases(F::ElectricFields.AbstractField, t::AbstractVector)
-    At = typeof(vector_potential(F, t[1]))
-    Tt = eltype(t)
-    nt = length(t)
-    ∫Av = zeros(At, nt)
-    ∫A²v = zeros(Tt, nt)
-
-    tend = t[end]
+function _VolkovPhases(F::ElectricFields.AbstractField, t, tend)
     Aend = ElectricFields.vector_potential(F, tend)
 
     A = t -> ElectricFields.vector_potential(F, t) - Aend
-    ∫A = AccumulatedIntegral(A, At, Tt; t=tend)
-    ∫A² = AccumulatedIntegral(t -> sum(abs2, A(t)), Tt, Tt; t=tend)
-
-    for i = nt-1:-1:1
-        set!(∫A, t[i])
-        set!(∫A², t[i])
-        ∫Av[i] = ∫A.∫f
-        ∫A²v[i] = ∫A².∫f
-    end
+    ∫A = AutomaticIntegration(A, t, tend)
+    square(z) = z^2
+    ∫A² = AutomaticIntegration(t -> sum(square, A(t)), t, tend)
 
-    VolkovPhases(t .- tend, ∫Av, ∫A²v)
+    VolkovPhases(tend, ∫A, ∫A²)
 end
 
-# This is second-order accurate, assuming Av[1] == Av[end] == 0
-function VolkovPhases(Av::AbstractVector{At}, t::AbstractVector{Tt}) where {At,Tt}
-    nt = length(t)
-    dt = step(t)
-    ∫Av = zeros(At, nt)
-    ∫A²v = zeros(Tt, nt)
-
-    tend = t[end]
-    Aend = Av[end]
+VolkovPhases(F::ElectricFields.AbstractField, t::AbstractVector) =
+    _VolkovPhases(F, t, t[end])
 
-    ∫Av[end] = -Av[end]*dt
-    ∫A²v[end] = -sum(abs2, Av[end])*dt
+VolkovPhases(F::ElectricFields.AbstractField, tre::AbstractVector, tim::AbstractVector) =
+    _VolkovPhases(F, ComplexPlane(tre, tim), tre[end])
 
-    for i = nt-1:-1:1
-        dA = (Av[i] + Av[i+1])/2
-        dA² = sum(abs2, dA)
-        ∫Av[i] = ∫Av[i+1] - dA*dt
-        ∫A²v[i] = ∫A²v[i+1] - dA²*dt
-    end
-
-    VolkovPhases(t .- tend, ∫Av, ∫A²v)
-end
+volkov_phase_k²(k, v::VolkovPhases, t) = norm(k)^2*t
 
-volkov_phase_k²(k, v, i) = norm(k)^2*v.t[i]
+kdotA(k::Number, A::Number) = k*A
+kdotA(k::SVector{3}, A::Number) = k[3]*A
+kdotA(k::SVector{3}, A::SVector{3}) = dot(k, A)
+volkov_phase_2kA(k, v::VolkovPhases, t) = 2kdotA(k, v.∫A(t))
 
-volkov_phase_2kA(k::Number, v::VolkovPhases{<:Any,<:Number,<:Number}, i) =
-    2k*v.∫A[i]
+volkov_phase_A²(v::VolkovPhases, t) = v.∫A²(t)
 
-volkov_phase_2kA(k::SVector{3}, v::VolkovPhases{<:Any,<:Number,<:Number}, i) =
-    2k[3]*v.∫A[i]
-
-volkov_phase_2kA(k::SVector{3}, v::VolkovPhases{<:Any,<:SVector{3},<:Number}, i) =
-    2dot(k, v.∫A[i])
-
-volkov_phase_A²(v, i) = v.∫A²[i]
-
-volkov_phase(k, v::VolkovPhases, i) =
-    -(volkov_phase_k²(k, v, i) +
-      volkov_phase_2kA(k, v, i) +
-      volkov_phase_A²(v, i))/2
+function volkov_phase(k, v::VolkovPhases, t)
+    # We only have to subtract the reference time from the integral
+    # over k², since it has been properly accounted for when setting
+    # up the integrals over A and A².
+    -(volkov_phase_k²(k, v, t - v.T) +
+      volkov_phase_2kA(k, v, t) +
+      volkov_phase_A²(v, t))/2
+end
 
-volkov_phase(k, v::VolkovPhases, i, j) =
-    -(volkov_phase(k, v, i) -
-      volkov_phase(k, v, j))
+volkov_phase(k, v::VolkovPhases, a, b) =
+    -(volkov_phase(k, v, a) -
+      volkov_phase(k, v, b))
 
-function stationary_momentum(v::VolkovPhases, i, j)
-    τ = v.t[i] - v.t[j]
-    -1/τ*(v.∫A[i]-v.∫A[j])
+function stationary_momentum(v::VolkovPhases, a, b)
+    τ = a - b
+    -1/τ*v.∫A(b, a)
 end