FeynmanAthermal (orig 1955 version) vw numeric solution:

Also an example to show comparison with perturbative stuff,
Examples/FeynmanAthermalAsymptoticComparison.jl

Test cross-checks against Schultz numeric solutions:
test/FeynmanAthermal.jl

All this work was previously done in a Jupyter notebook ahead of my APS
March talk; see here:
https://github.com/jarvist/ijulia-notebooks/blob/master/2018-02-Athermal%20Feynman%20polarons.ipynb

Examples/FeynmanAthermalAsymptoticComparison.jl

@@ -0,0 +1,68 @@
+# FeynmanAthermalAsymptoticComparison.jl - reproduce some 1950s numeric results
+
+push!(LOAD_PATH,"../src/") # load module from local directory
+using PolaronMobility 
+
+using Plots
+
+# Here we collect all the data we need for plotting
+αrange=0.2:0.1:10 # unstable below α=0.2
+vs=[ feynmanvw(α)[1]  for α in αrange] # takes a while; computes all vw
+ws=[ feynmanvw(α)[2]  for α in αrange] 
+# Urgh - can't seem to do this more elegantly; recalculating all the params just to strip the second component of the tuple
+
+Es=map(F,vs,ws,αrange); # calculate polaron energy for these params
+
+p=plot(xlim=αrange,ylim=(0,), 
+    xlabel="α parameter", ylabel="Reduced polaron units", 
+    title="Athermal (numeric, variational) v,w, and E for α" )
+
+plot!(αrange,vs,label="v")
+plot!(αrange,ws,label="w")
+plot!(αrange,Es,label="E")
+hline!([0.0],label="",color=:black)
+
+# I like this plot!
+# You can see how v=w=3 as α-->0.0
+# w then tends down to 1.0 in a nice sigmoid form as α --> ~10+
+# v increases monotonically
+
+savefig("AthermalvwE.pdf")
+
+# Feynman Stat Mech (1972), p. 240
+smallαw(α)=3.0
+P(α)=2/smallαw(α) * (sqrt(1+smallαw(α))-1) # Nb: error in Feynman Stat Mech: sqrt(1+W) , as in Feynman1955
+smallαv(α)=3*(1 + 2*α*(1-P(α))/(3*smallαw(α))) 
+
+c=0.5772 # Euler Mascheroni constant
+
+largeαw(α)=1.0
+largeαv(α)=(4*α^2/9π) - (4*(log(2) + c/2)-1)
+
+p=plot(xlim=αrange,ylim=(0,15), xlabel="α parameter", ylabel="v,w, Reduced polaron units", title="Asymtotic limits for v,w" )
+
+plot!(αrange,α->smallαw(α),label="w, smallαw",style=:dash)
+plot!(αrange,α->smallαv(α),label="v, smallαv",style=:solid)
+
+plot!(αrange,α->largeαw(α),label="w, largeαw",style=:dash)
+plot!(αrange,α->largeαv(α),label="v, largeαv",style=:solid)
+
+plot!(αrange,vs,label="v, numeric",w=3,style=:solid)
+plot!(αrange,ws,label="w, numeric",w=3,style=:dash)
+
+savefig("AthermalvwAsymptotic.pdf")
+
+smallαE(α)=-α-α^2/81
+FeynmanStatMechlargeαE(α)=-α^2/(3π) - 3/2 *(2*log(2)+c)-3/4 # As in Feynman Stat Mech; with 2π corrected to 3π ...
+
+Feynman1955largeαE(α)=-α^2/3π - 3*log(2) # as in Feynman I, 1955
+
+p=plot(xlim=αrange,ylim=(-Inf,0), xlabel="α parameter", ylabel="Energy, Reduced polaron units", title="Asymtotic limits for E(α)" )
+
+plot!(αrange,α->smallαE(α),label="smallαE")
+plot!(αrange,α->FeynmanStatMechlargeαE(α),label="FeynmanStatMechlargeαE")
+plot!(αrange,α->Feynman1955largeαE(α),label="Feynman1955largeαE")
+plot!(αrange,Es,label="E (numeric integration)")
+
+savefig("AthermalAsymptoticEnergy.pdf")
+

src/FeynmanKadanoffOsakaHellwarth.jl

@@ -5,7 +5,7 @@
 # These codes were developed with Julia 0.5.0 - Julia 0.6.2, and require the Optim and Plots packages.
 
 export Polaron # Type to hold the data
-export feynmanalpha, polaronmobility, savepolaron, plotpolaron
+export feynmanalpha, feynmanvw, F, polaronmobility, savepolaron, plotpolaron
 export HellwarthBScheme, HellwarthAScheme
 export ImX
 
@@ -112,6 +112,7 @@ end
 # Set up equations for the polaron free energy, which we will variationally improve upon
 
 # Hellwarth et al. 1999 PRB - Part IV; T-dep of the Feynman variation parameter
+
 # Originally a Friday afternoon of hacking to try and implement the T-dep electron-phonon coupling from the above PRB
 # Which was unusually successful! And more or less reproduced Table III
 
@@ -135,6 +136,40 @@ F(v,w,β,α)=-(A(v,w,β)+B(v,w,β,α)+C(v,w,β))
 # BUT - this is just the objective function! (The temperature-dependent free-energy.) Not the optimised parameters.
 # Also there's a scary numeric integration (quadgk) buried within...
 
+# In Julia we have 'Multiple dispatch', so let's just construct the Feynman (athermal) 
+# energy as the same signature, but without the thermodynamic beta
+
+# Integrand of (31) in Feynman I (Feynman 1955, Physical Review, "Slow electrons...")
+fF(τ,v,w)=(w^2 * τ + (v^2-w^2)/v*(1-exp(-v*τ)))^-0.5 * exp(-τ)
+# (31) in Feynman I
+AF(v,w,α)=π^(-0.5) * α*v * quadgk(τ->fF(τ,v,w),0,Inf)[1]
+# (33) in Feynman I
+F(v,w,α)=(3/(4*v))*(v-w)^2-AF(v,w,α)
+
+# Let's wrap this in a simple function
+"""
+    feynmanvw(α)
+
+    Calculate v and w variational polaron parameters (Feynman original athermal action), 
+    for the supplied α Frohlich coupling.
+	Returns v,w. 
+"""
+function feynmanvw(α)
+    initial=[7.0,6.0]
+    # Main use of these bounds is stopping v or w going negative, at which you get a NaN error as you are evaluating log(-ve Real)
+    lower=[0.1,0.1]
+    upper=[1000.0,1000.0]
+
+    myf(x) = F(x[1],x[2],α) # Wraps the function so just the two variational params are exposed
+
+    res=optimize(OnceDifferentiable(myf, initial; autodiff = :forward), initial, lower, upper, Fminbox(); optimizer=BFGS)
+    # allow_f_increases=true,  - increases stability of algorith, but makes it more likely to crash as it steps outside Fminbox(), le sigh.
+
+    v,w=Optim.minimizer(res)
+
+    return v,w
+end
+
 # Structure to store data of polaron solution + other parameters, for each temperature
 struct Polaron
     T

test/FeynmanAthermal.jl

@@ -0,0 +1,25 @@
+println("FeynmanAthermal Tests - check athermal feynmanvw(α) -> v,w literature values.")
+
+# Results from Feynman & Hibbs, Emended Edition, p 319.
+Schultz=[
+3.00 3.44 2.55 -3.1333
+5.00 4.02 2.13 -5.4401
+7.00 5.81 1.60 -8.1127 
+9.00 9.85 1.28 -11.486
+11.0 15.5 1.15 -15.710
+]
+
+rows(M::Matrix) = map(x->reshape(getindex(M, x, :), :, size(M)[2]), 1:size(M)[1])
+
+for r in rows(Schultz)
+    α,vSchultz,wSchultz,ESchultz=r # Unpack row of results
+    v,w=feynmanvw(α) # performs the optimisation
+    E=F(v,w,α) # energy at the optimised parameters
+
+    println("α=$α v=$v w=$w E=$E  | Schultz: v=$vSchultz w=$wSchultz E=$ESchultz")
+
+    @test v ≈ vSchultz atol=0.1 # Strangely these need more tolerance than the Energies
+    @test w ≈ wSchultz atol=0.1
+    @test E ≈ ESchultz atol=0.001
+end
+

test/runtests.jl

@@ -4,8 +4,9 @@ using PolaronMobility
 using Base.Test
 
 include("FeynmanAlpha.jl") # Simple test of Frohlich alpha vs. literature values
+include("FeynmanAthermal.jl") # Athermal Feynman tests
 include("HellwarthMultipleBranches.jl") # Test Hellwarth et al. 1999 PRB 'B' multiple branch reduction scheme
 include("FrostPolaronMobility2017.jl") # Reproduce values published in Frost 2017 PRB 
 
-println("That's me! If I finished without interupting, let's assume the tests were A.OK.")
+println("\nThat's me! If I finished without interupting, all tests have passed.")