On the Widder inversion method in problems of statistical mechanics

An alternative way of reconstructing a function from its Laplace transform using the Widder inversion method was shown to be useful in treating some problems of nonequilibrium statistical mechanics. As an example of a successful application of the method, a decay was investigated of the excited atomic state in a simple, but nonetheless physically relevant, model featuring a two-level atom interacting with a continuum of field modes.


Widder inversion formula
It is well known that many physical problems, being treated by methods of nonequilibrium statistical mechanics, can be reduced to a system of differential or integro-differential equations.Sometimes this system can be simplified by means of the direct Laplace transformation which results in the corresponding system of algebraic equations.If this algebraic system permits an explicit solution, then the corresponding solution of the initial system of equations can be reconstructed in principle with the aid of the inverse Laplace transformation.
It is generally believed that the inverse Laplace transformation has something to do exclusively with the theory of complex numbers and the analysis of the functions of complex variables.Thus if f (z) is a Laplace transformation of some function φ(t) where the integral converges for some complex value of z, then the inverse Laplace transformation can be carried out in the form where c is a sufficiently large real constant and the integration runs along a vertical line.
An alternative approach to calculating the inverse Laplace transformations was worked out in detail by Widder [1].Among numerous statements related to the problem in question he proved, in particular, two important theorems.The first theorem, being constructive by nature, states that a convergent sequence of subsequent approximations, which represents an alternative way of reconstructing the function φ(t), can be built explicitly: Theorem 1.If the function φ(t) is integrable in the interval (0, R) for every positive R and if the integral for almost all positive values of t.Here the operator L t is defined as the limit where and f (k) (x) stands for the derivative of the k-th order.It is seen from this theorem that both the direct and the inverse Laplace transformations can rely on real variables t and x only.
The second theorem reveals how this sequence of approximations L k,t [f (x)] converges to the function φ(t): Theorem 2. If the function φ(t) is integrable in the interval (0, R) for every positive R, if it is continuous in the interval 0 a t b, and if the integral converges for some value of x, then uniformly in the interval a ′ t b ′ , where a < a ′ < b ′ < b.
The usage of the conventional inverse transformation (1) assumes the Bromwhich contour integration in the complex z-plane as a rule.Consequently, the function f (z) must be treated as a function of the complex variable and its properties must be thoroughly investigated by means of the corresponding methods of complex analysis.Such an analysis can be very complicated and laborious, depending on the complexity of a particular Laplace transform f (z).
Direct numerical methods, when applied to the integral (1), also come across some difficulties because of the highly oscillatory behaviour of the integrand, especially for large t.Moreover, any numerical method provides us with anything but a set of separate values of φ(t) calculated for several singular points t and for a given set of values of other relevant parameters only.
Contrariwise, the algorithm (2), wherever applicable, permits to circumvent both difficulties.From the analytical point of view, its implementation does not require, as a rule, a detailed knowledge of the functional properties of f (z).From the numerical point of view, some explicit analytical expression exists at each step k.This expression approximates not only singular values but the function φ(t) as a whole and thus may be used as a substitute for this function in a variety of necessities.
As to the convergence speed, it can be shown under the conditions of the theorem 1, that, for any given t, the sequence (3) converges not slower than 1/k at least starting from some k = k 0 .And this k 0 can be chosen one and the same for all t in the interval a ′ t b ′ indicated in the conditions of the theorem 2.
Technically, calculation of the approximation L k,t [f (x)] is uniform at any given order k and comprises, as its essential part, a subsequent k-fold differentiation of the function f (x) only.The mere need to calculate this derivative in case of large enough k could pose a considerable challenge even not so long ago -and it seems to be the only reason why the Widder inverse transformation has been neglected for a long time -but the situation has changed by now and keeps on changing for the better with the appearance and further rapid development of programming tools allowing to unite analytical and numerical computing within one integrated environment.And all this is accompanied by the growing raw computational power of every new generation of computers.

Atom-field interaction model
To show the efficiency of the Widder inversion method in solving actual physical problems we choose, as an example, a well-known model of a two-level atom interacting with a continuum of radiation field modes.The Hamiltonian of this model is where and 55 × 10 16 s −1 is the energy separation between 2P and 1S atomic states, λ = γω 0 (γ is the Einstein coefficient for spontaneous Lyman-α transition), Ω = 3/(2a 0 ) = (3/2) αm e c 2 = 8.498 × 10 18 s −1 (a 0 is the Bohr radius) is the natural cutoff frequency, α is the fine-structure constant, m e is the electron mass and I A and I R are the unit operators in the corresponding subspaces of the atom and the radiation field.Here we used as the excited and as the ground state of an hydrogenic atom.The following normalization rules hold for one-photon states: where |V is the vacuum state and the notation of [3] has been used.It is well known that if we are only interested in the quantum mechanical evolution of the excited state |1 ⊗ |V then the model ( 6) permits a further reduction to the one-photon subspace.In this case only this excited state and the one-photon states |2 ⊗ |ω should be taken into account which results in the following effective Hamiltonian: Later on we deal with just this effective Hamiltonian as with an independent sustainable self-consistent quantum model without invoking the original model (6).

Exact non-decay probability amplitude
As shown in [2,3], an exact equation for the probability amplitude b 1 (t ) for finding a two-level hydrogenic atom in the excited state and no photons in the radiation field is where F (ω) is the natural smooth cutoff function To solve equation (9) the Laplace transformation b1 (z) = where b 1 (0) is an initial value of b 1 (t).Following [3], we substitute u = iz + ω 0 and, integrating over ω, introduce a new function where the notation log(u) stands for the multivalued natural logarithm of a complex variable.The Laplace inversion of b1 (z) yields (after substitution z = −iu where C is some contour in the complex u-plane (see figure 1 in [3]).Further simplification of equation ( 10) requires a laborious analysis based on the theory of complex variables and lengthy transformations on some Riemann surface which were fulfilled in [3].Here, we only outline the final result consisting of the residue term and the so-called integral cut contribution: where and the value u 0 is the Weisskopf-Wigner pole of B1 (u) which lies on the lower Riemann sheet of log(u).Parametrization u = (1 − i)Ωs is used in the integral term and where the subscripts 0 and −1 denote different branches of the functions N(u) and I(u) on the Riemann surface of log(u) (see [3] for details).Usually all frequencies are expressed in units of the characteristic frequency Ω for convenience.Thus in the present case the characteristic time unit would be 1/Ω.In terms of the scaled variables the probability amplitude b 1 (t) reads where ū = (1 − i)s.Correspondingly where ω0 = ω 0 /Ω, ū = u/Ω, ū0 = u 0 /Ω, t = Ωt and χ = λ/2π.Integrating (13) two times by parts, we obtain an identity Here .
The right-hand side of equation ( 14) can be transformed into where s 0 is some arbitrary positive parameter.For a given t one can choose this parameter in such a way that the second integral term in equation ( 15) is negligible while the integrand in the first integral term does not oscillate very fast at the same time.The first term in equation ( 15) represents asymptotic nonexponential behaviour of the probability amplitude for t → +∞ predicted in [5][6][7] and thoroughly investigated in [3].

Calculations with the Widder inversion formula
In [4] the probability P ( t ) = |b 1 ( t )| 2 for finding a two-level hydrogenic atom in the excited state and no photons in the radiation field was investigated in the time interval 0 t 10000 based on equation ( 12).Here, the same calculations were carried out by means of the Widder inversion formula (2,4) for k = 15, 20.Results are presented as corresponding graphs plotted in solid lines against the same results obtained from equation ( 12 The main feature of all these graphs is that the approximate solution (2) converges uniformly to the exact one (12) from below.
It is also seen that at the beginning of the time interval 0 t 1 both graphs nearly coincide already for k = 20.Thus the approximate solution (2) is able to reproduce deviation from the exponential Weisskopf-Wigner law at the initial time interval correctly.This deviation may be characterized by a special parameter -the so-called Zeno time τ z -if one approximates the actual behaviour of the probability by the polynomial function P ( t ) = 1 − t 2 /τ 2 z as it was done in [8].
All calculations were carried out with Mathematica 3.0 computing package on a more than humble hardware -Pentium-100 processor with 128 Mb RAM.Higher orders of approximation k > 20 are accessible even for this computer but on the considerable time expense.For the sake of better precision, Mathematica's internal system variables $MinP recision and $MaxExtraP recision were put to 60.This measure probably allowed to avoid possible divergencies when calculating the Widder inversion.For the same purpose, all numerical parameters, including the scaled time t, entered in calculations as rational numbers making Mathematica to assume they were exact ones and thus forcing high-precision calculations.

∞ 0 b 1
(t)e −tz dt can be applied and the resulting algebraic equation for the function b1 (z) can be solved: b1 ) which are plotted in dashed lines.Deviations from the exact solution ∆( t ) = |b 1 ( t )| 2 − |b Widder 1 ( t )| 2 are plotted separately for the same time intervals.

Figure 1 .
Figure 1.Non-decay probability P ( t ) for k = 15 and k = 20 for various time intervals.