Functions

A_j(β::Float64, v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1), n::Float64)

Generalisation of the A function from Equation 62b in Hellwarth and Biaggio[1]. This is
the Helmholtz free energy of the trial model.

 - β is the reduced thermodynamic temperature ħ ω_j / (kB T) associated with the `jth` phonon mode.
 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.
 - n is the total number of phonon modes.

[1] Hellwarth, R. W., Biaggio, I. (1999). Mobility of an electron in a multimode polar
 lattice. Physical Review B, 60(1), 299–307. https://doi.org/10.1103/PhysRevB.60.299

A_j(v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1), n::Float64)

Generalisation of the A function from Equation 62b in Hellwarth Biaggio[1] but to
zero temperature. This is the ground-state energy of the trial model.

 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.
 - n is the total number of phonon modes.

[1] Hellwarth, R. W., Biaggio, I. (1999). Mobility of an electron in a multimode polar
 lattice. Physical Review B, 60(1), 299–307. https://doi.org/10.1103/PhysRevB.60.299

B_j(α::Float64, β::Float64, v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1))

Generalisation of the B function from Equation 62c in Hellwarth and Biaggio[1]. This is
the expected value of the exact action <S_j> taken w.r.t trial action, given for the
`jth` phonon mode.

 - α is the partial dielectric electron-phonon coupling parameter for the `jth` phonon
 mode.  
 - β is the reduced thermodynamic temperature ħ ω_j / (kB T) associated with the
 `jth` phonon mode.  
 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.

 [1] Hellwarth, R. W., Biaggio, I. (1999). Mobility of an electron in a multimode polar
 lattice. Physical Review B, 60(1), 299–307. https://doi.org/10.1103/PhysRevB.60.299

B_j(α::Float64, v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1))

Generalisation of the B function from Equation 62c in Biaggio and Hellwarth [1] to zero
temperature. This is the expected value of the exact action <S_j> taken w.r.t trial
action, given for the `jth` phonon mode.

 - α is the partial dielectric electron-phonon coupling parameter for the `jth` phonon mode.
 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.

 [1] Hellwarth, R. W., Biaggio, I. (1999). Mobility of an electron in a multimode polar
 lattice. Physical Review B, 60(1), 299–307. https://doi.org/10.1103/PhysRevB.60.299

C_ij(i::Integer, j::Integer, v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1))

Calculates the element to the coupling matrix C_ij (a generalisation of Feynman's `C`
coupling variational parameter) between the electron and the `ith' and `jth' fictitious
masses that approximates the exact electron-phonon interaction with a harmonic coupling
to a massive fictitious particle. NB: Not to be confused with the number of physical
phonon branches; many phonon branches could be approximated with one or two etc.
fictitious masses for example. The number of fictitious mass does not necessarily need
to match the number of phonon branches.

 - i, j enumerate the current fictitious masses under focus (also the index of the element in the coupling matrix C)
 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.

C_j(β::Float64, v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1), n::Float64)

Generalisation of the C function from Equation 62e in Hellwarth and Biaggio[1]. This is
the expected value of the trial action <S_0> taken w.r.t trial action.

 - β is the reduced thermodynamic temperature ħ ω_j / (kB T) associated with the `jth`
 phonon mode.  
 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.
 - n is the total number of phonon modes.

 [1] Hellwarth, R. W., Biaggio, I. (1999). Mobility of an electron in a multimode polar
 lattice. Physical Review B, 60(1), 299–307. https://doi.org/10.1103/PhysRevB.60.299

C_j(v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1), n::Float64)

Generalisation of the C function from Equation 62e in Biaggio and Hellwarth [] but to
zero temperaure. This is the expected value of the trial action <S_0> taken w.r.t trial
action.

 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.
 - n is the total number of phonon modes.

[1] Hellwarth, R. W., Biaggio, I. (1999). Mobility of an electron in a multimode polar
 lattice. Physical Review B, 60(1), 299–307. https://doi.org/10.1103/PhysRevB.60.299

D_j(τ::Float64, β::Float64, v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1))

Calculates the recoil function that approximates the exact influence (recoil effects) of
the phonon bath on the electron with the influence of the harmonicly coupled fictitious
masses on the electron. It appears in the exponent of the intermediate scattering
function.

 - τ is the imaginary time variable.
 - β is the reduced thermodynamic temperature ħ ω_j / (kB T) associated with the `jth` phonon mode.
 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.

D_j(τ::Float64, v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1))

Calculates the recoil function that approximates the exact influence (recoil effects) of
the phonon bath on the electron with the influence of the harmonicly coupled fictitious
masses on the electron. It appears in the exponent of the intermediate scattering
function. This function works at zero temperature.

 - τ is the imaginary time variable.
 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.

DieletricFromIRmode(IRmode)

Calculate dielectric from an individual mode.

IRmode is a tuple f,S with Frequency in THz; InfraRed activity in e^2 amu^-1

HellwarthAScheme(phonon_modes; T=295, convergence=1e-6)

Multiple phonon mode reduction to a single effective frequency.
Temperature dependent, defaults to T = 295 K.

Solved iteratively by bisection until Δfreq<convergence.

Follows Hellwarth et al. 1999 PRB 'A' scheme, Eqn 50 RHS.

HellwarthBScheme(LO)
LO an array assumed to be of [freq ; absolute ir activity ]

Multiple phonon mode reduction to a single effective frequency.
Hellwarth et al. 1999 PRB, 'B scheme'; the athermal method.  Averaging procedure is
constructed by considering the average effect of the action of multiple branches.

Follows Eqn (58) in this paper, assuming typo on LHS, should actually be W_e.

IRtoDielectric(IRmodes,volume)

From absolute value of IR activities of phonon modes, generate a per-mode contribution to the low-frequency dielectric constant.

IRmodes are tuples f,S with Frequency in THz; InfraRed activity in e^2 amu^-1

IRtoalpha(IR,volume)

Calculates contribution to dielectric constant for each polar phonon mode, and thereby the Frohlich alpha contribution for this mode.

function ImX(nurange,v,w,βred,α,ω,mb)

Impedance in (47a) from Feynman1962, directly solving freq dep without taking Hellwarth1999 limit of v->0 .

Calculates a frequency dependent (over range of nu) susceptibility which can be linked back to mobility.

HERE BE DRAGONS! Not well tested or validated code; the central numeric integration is highly oscillatory and becomes intractable for large nu.

effectivefreqs(freqsandiractivity::Matrix{Float64}, numvarparams::Integer)

Generates a matrix of effective phonon modes with frequencies and infra-red activities
derived from a larger matrix using the Principal Component Analysis (PCA) method.

 - freqs_and_ir_activity: is a matrix containing the phonon mode frequencies (in THz) in
 the first column and the infra-red activities (in e^2 amu^-1) in the second column.
 - num_var_params: is the number of effective modes required (which needs to be less
 than the number of modes in freqs_and_ir_activity)

*** POSSIBLY REDUNDANT ***

feynmanvw(α, β; v = 0.0, w = 0.0)

Calculate v and w variational polaron parameters, for the supplied α Frohlich coupling,
and inverse temperature β.  
This version uses the Osaka thermal action symmetrised for computation (Hellwarth 1999).
Returns v,w.

feynmanvw(α; v = 0.0, w = 0.0)

Calculate v and w variational polaron parameters,
for the supplied α Frohlich coupling.
This version uses the original athermal action (Feynman 1955).
Returns v,w.

frohlichPartial((f, ϵmode); ϵo,ϵ_s,meff)

Calculate a (partial) dielectric electron-phonon coupling element. f - frequency of mode in THz ϵmode - this mode's contribution to dielectric ϵo - optical dielectric ϵ_s - total static dielectric contribution

frohlichalpha(εInf,εS,freq,m_eff)

Calculates the Frohlich alpha parameter, for a given dielectric constant,
frequency (f) of phonon in Hertz, and effective mass (in units of the
bare electron mass).

See Feynman 1955:
http://dx.doi.org/10.1103/PhysRev.97.660

h_i(i::Integer, v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1))

Calculates the normal-mode (the eigenmodes) frequency of the coupling between the
electron and the `ith' fictitious mass that approximates the exact electron-phonon
interaction with a harmonic coupling to a massive fictitious particle. NB: Not to be
confused with the number of physical phonon branches; many phonon branches could be
approximated with one or two etc. fictitious masses for example. The number of
fictitious mass does not necessarily need to match the number of phonon branches.

 - i enumerates the current fictitious mass.
 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.

load_polaron(p::NewPolaron, prefix)

Loads data from file "polaronfilepath" into a NewPolaron type.

makepolaron(ϵoptic, ϵstatic, phononfreq, meff; temp = 300.0, efieldfreq = 0.0, volume = nothing, iractivity = nothing, Nparams = 1, rtol = 1e-3, verbose = true)

Solves the Feynman polaron problem variationally with finite temperature
Osaka energies. From the resulting v and w parameters, calculates polaron
structure (wave function size, etc.).  Uses FHIP to predict
a temperature- and frequency-dependent mobility, complex conductivity and impedance for the material system.
Can evaluate polaron properties for materials with multiple phonon branches using infrared activities and phonon branch frequencies.

Inputs are:
• temperature range (e.g. 10:50:1000),
• reduced dielectric constants (e.g. 5, 20),
• characteristic dielectric phonon frequency (e.g. 2.25) - units TeraHertz
• bare-band effective-mass (e.g. 0.12) - units electron mass
• electric field frequency range - units TeraHertz
• unit cell volume for the material - units m³
• infrared activities - units ?
• phonon mode frequencies - units TeraHertz
• number of variational parameters to minimise the polaron energy
• relative tolerance for the accuracy of any quadrature integrations

Returns a structure of type NewPolaron, containing arrays of useful
information.  Also prints a lot of information to the standard out - which
may be more useful if you're just inquiring as to a particular data point,
rather than plotting a temperature-dependent parameter.

As an example, to calculate the electron polaron in MAPI, at temperatures 0:100:400 K and electric field frequencies 0.0:0.1:5.0 THz, and inclusive of 15 optical phonon modes:
make_polaron(
    4.5, 
    24.1, 
    [4.016471586720514, 3.887605410774121, 3.5313112232401513, 2.755392921480459, 2.4380741812443247, 2.2490917637719408, 2.079632190634424, 2.0336707697261187, 1.5673011873879714, 1.0188379384951798, 1.0022960504442775, 0.9970130778462072, 0.9201781906386209, 0.800604081794174, 0.5738689505255512], 
    0.12; 
    temp = 0.0:100.0:400.0, 
    efield_freq = 0.0:0.1:5.0, 
    volume = (6.29e-10)^3, 
    ir_activity = [0.08168931020200264, 0.006311654262282101, 0.05353548710183397, 0.021303020776321225, 0.23162784335484837, 0.2622203718355982, 0.23382298607799906, 0.0623239656843172, 0.0367465760261409, 0.0126328938653956, 0.006817361620021601, 0.0103757951973341, 0.01095811116040592, 0.0016830270365341532, 0.00646428491253749], 
    N_params = 1)

makepolaron(ϵoptic, ϵstatic, phononfreq, meff; temp = 300.0, efieldfreq = 0.0, volume = nothing, iractivity = nothing, Nparams = 1, rtol = 1e-3, verbose = true)

Same as above but from a model system with a specified alpha value rather than from material properties.

Inputs are:
• temperature range (e.g. 10:50:1000),
• Frohlich alpha parameter,
• electric field frequency range - units TeraHertz
• phonon mode frequency weight (choose either 1 or 2π)
• number of variational parameters to minimise the polaron energy
• relative tolerance for the accuracy of any quadrature integrations

Returns a structure of type NewPolaron, containing arrays of useful
information.  Also prints a lot of information to the standard out - which
may be more useful if you're just inquiring as to a particular data point,
rather than plotting a temperature-dependent parameter.

As an example, to calculate the polaron of a α = 6 system, at temperatures 0:100:400 K and electric field frequencies 0.0:0.1:5.0 THz:
make_polaron(6.0; temp = 0.0:100.0:400.0, efield_freq = 0.0:0.1:5.0)

multi_F(v, w, α, β; ω = 1.0)

Calculates the Helmholtz free energy of the polaron for a material with multiple phonon
branches. This generalises Osaka's free energy expression (below Equation (22) in [1]).

 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.
 - α is the Frohlich coupling constant.
 - β is the thermodynamic temperature.

[1] Osaka, Y. (1959). Polaron State at a Finite Temperature. Progress of Theoretical
Physics, 22(3), 437–446. https://doi.org/10.1143/ptp.22.437

multi_F(v, w, α; ω)

Calculates the Helmholtz free energy of the polaron for a material with multiple phonon
branches. This generalises Osaka's free energy expression (below Equation (22) in [1]).

    - α is the Frohlich alpha parameter.
    - ω is a vector containing the phonon mode frequencies (in THz).
    - v and w determines if the function should start with a random initial set of
    variational parameters (= 0.0) or a given set of variational parameter values.  

[1] Osaka, Y. (1959). Polaron State at a Finite Temperature. Progress of Theoretical
Physics, 22(3), 437–446. https://doi.org/10.1143/ptp.22.437

multifrohlichalpha(ϵoptic::Float64, ϵionic::Float64, ϵtotal::Float64, phononmodefreq::Float64, m_eff::Float64)

Calculates the partial dielectric electron-phonon coupling parameter for a given
longitudinal optical phonon mode. This decomposes the original Frohlich alpha coupling
parameter (defined for a single phonon branch) into contributions from multiple phonon
branches.

 - ϵ_optic is the optical dielectric constant of the material.
 - ϵ_ionic is the ionic dielectric contribution from the phonon mode.
 - ϵ_total is the total ionic dielectric contribution from all phonon modes of the material.
 - phonon_mode_freq is the frequency of the phonon mode (THz).
 - m_eff is the band mass of the electron (in units of electron mass m_e)

optical_absorption(Ω::Float64, β::Float64, α::Float64, v::Float64, w::Float64; rtol = 1e-3)

Calculate the absorption coefficient Γ(Ω) for the polaron at at finite temperatures
(equation (11a) in [1]) for a given frequency Ω.  β is thermodynamic beta. v and w are
the variational Polaron parameters that minimise the free energy, for the supplied
α Frohlich coupling. rtol specifies the relative error tolerance for the QuadGK integral
in the memory function.  

[1] Devreese, J., De Sitter, J., & Goovaerts, M. (1972). Optical Absorption of Polarons
in the Feynman-Hellwarth-Iddings-Platzman Approximation. Physical Review B, 5(6),
2367–2381. doi:10.1103/physrevb.5.2367

polaroncomplexconductivity(Ω::Float64, β::Float64, α::Float64, v::Float64, w::Float64)

Calculate the complex conductivity σ(Ω) of the polaron at finite temperatures for
a given frequency Ω (equal to 1 / Z(Ω) with Z the complex impedence). β is the
thermodynamic beta. v and w are the variational polaron parameters that minimise the
free energy, for the supplied α Frohlich coupling. rtol specifies the relative error
tolerance for the QuadGK integral in the memory function.

polaroncompleximpedence(Ω::Float64, β::Float64, α::Float64, v::Float64, w::Float64)

Calculate the complex impedence Z(Ω) of the polaron at finite temperatures for a given
frequency Ω (equation (41) in FHIP 1962 [1]). β is the thermodynamic beta. v and w are
the variational polaron parameters that minimise the free energy, for the supplied
α Frohlich coupling. rtol specifies the relative error tolerance for the QuadGK
integral in the memory function.

polaronmemoryfunction(Ω::Float64, β::Float64, α::Float64, v::Float64, w::Float64)

Calculate the memory function χ(Ω) of the polaron at finite temperatures (equation (35)
in FHIP 1962 [1]) for a given frequency Ω. This function includes the zero-temperature
and DC limits. β is the thermodynamic beta. v and w are the variational polaron
parameters that minimise the free energy, for the supplied α Frohlich coupling. rtol
specifies the relative error tolerance for the QuadGK integral and corresponds to the
error of the entire function. 
    
Finite temperature and finite frequency memory function, including the limits to zero
frequency Ω → 0 or zero temperature β → ∞.  

[1] R. P. Feynman, R. W. Hellwarth, C. K. Iddings, and P. M. Platzman, Mobility of slow
electrons in a polar crystal, PhysicalReview127, 1004 (1962)

function multimemoryfunction(Ω::Float64, β::Array{Float64}(undef, 1), α::Array{Float64}(undef, 1), v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1), ω::Array{Float64}(undef, 1), m_eff::Float64)

Calculate polaron complex memory function inclusive of multiple phonon branches j, each
with angular frequency ω[j] (rad THz).

 - Ω is the frequency (THz) of applied electric field.
 - β is an array of reduced thermodynamic betas, one for each phonon frequency ω[j]. 
 - α is an array of decomposed Frohlich alphas, one for each phonon frequency ω[j]. 
 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.
 - m_eff is the is the conduction band mass of the particle (typically electron / hole, in units of electron mass m_e).

polaron_mobility(β::Float64, α::Float64, v::Float64, w::Float64; rtol = 1e-3)

Calculate the dc mobility μ of the polaron at finite temperatues (equation (11.5) in
[3]) for a given frequency Ω. β is the thermodynamic beta. v and w are the variational
polaron parameters that minimise the free energy, for the supplied α Frohlich coupling.
rtol specifies the relative error for the integral to reach.

polaronmobility(Trange, εInf, εS, freq, effectivemass; verbose::Bool=false)

Solves the Feynman polaron problem variationally with finite temperature
Osaka energies.  From the resulting v, and w parameters, calculates polaron
structure (wave function size, etc.).  Uses FHIP, Kadanoff (Boltzmann
relaxation time) and Hellwarth direct contour integration to predict
a temperature-dependent mobility for the material system.
Input is a temperature range (e.g. 10:50:1000),
reduced dielectric constants (e.g. 5, 20),
characteristic dielectric phonon frequency (e.g. 2.25E12) - units Hertz
bare-band effective-mass (e.g. 012) - units electron mass.

Returns a structure of type Polaron, containing arrays of useful
information.  Also prints a lot of information to the standard out - which
may be more useful if you're just inquiring as to a particular data point,
rather than plotting a temperature-dependent parameter.

As an example, to calculate the electron polaron in MAPI at 300 K:
polaronmobility(300, 4.5, 24.1, 2.25E12, 0.12)

save_polaron(p::NewPolaron, prefix)

Saves data from 'polaron' into file "prefix". This is a .jdl file for storing the polaron data whilst preserving types. Allows for saving multidimensional arrays that sometimes arise in the polaron data. Each parameter in the NewPolaron type is saved as a dictionary entry. E.g. NewPolaron.α is saved under JLD.load("prefix.jld")["alpha"].

savepolaron(fileprefix, p::Polaron)

Saves data from polaron 'p' into file "fileprefix". This is a simple space-delimited text file, with each entry a separate temperature, for plotting with Gnuplot or similar.

Structure of file is written to the header:

Ts, βreds, Kμs, Hμs, FHIPμs, vs, ws, ks, Ms, As, Bs, Cs, Fs, Taus, rfsis

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

variation(α::Vector{Real}, β::Vector{Real}; v::Real, w::Real, ω::Vector{Real}, N::Integer)

Minimises the multiple phonon mode free energy function for a set of v_p and w_p
variational parameters.  The variational parameters follow the inequality: v_1 > w_1
> v_2 > w_2 > ... > v_N > w_N.

 - β is the thermodynamic temperature.
 - α is the Frohlich alpha parameter.
 - ω is a vecotr containing the phonon mode frequencies (in THz).
 - v and w determines if the function should start with a random initial set of variational parameters (= 0.0) or a given set of variational parameter values.
 - N specifies the number of variational parameter pairs, v_p and w_p, to use in minimising the free energy.

κ_i(i::Integer, v::Array{Float64}(undef, 1), w::Array{Float64}(undef, 1))

Calculates the spring-constant coupling the electron to the `ith' fictitious mass that
approximates the exact electron-phonon interaction with a harmonic coupling to a massive
fictitious particle. NB: Not to be confused with the number of physical phonon branches;
many phonon branches could be approximated with one or two etc. fictitious masses for
example. The number of fictitious mass does not necessarily need to match the number of
phonon branches.

 - i enumerates the current fictitious mass.
 - v is an one-dimensional array of the v variational parameters.
 - w is an one-dimensional array of the w variational parameters.

ϵionicmode(phononmodefreq::Float64, ir_activity::Float64, volume::Float64)

Calculate the ionic contribution to the dielectric function for a given phonon mode.
phonon_mode_freq is the frequency of the mode in THz.

 - ir_activity is the infra-red activity of the mode in e^2 amu^-1.
 - volume is the volume of the unit cell of the material in m^3.

ϵtotal(freqsandiractivity::Matrix{Float64}, volume::Float64)

Calculate the total ionic contribution to the dielectric function from all phonon modes.

 - freqs_and_ir_activity is a matrix containexing the phonon mode frequencies (in THz)
 in the first column and the infra-red activities (in e^2 amu^-1) in the second column.
 - volume is the volume of the unit cell of the material in m^3.