Persistence Exponent for the Simple Diffusion Equation: The Exact Solution for Any Integer Dimension

Devashish Sanyal

Department of Theoretical Physics, Institute of Physics, Bhubaneswar, India.

DOI: 10.4236/jmp.2021.1210083                

Abstract

The persistence exponent for the simple diffusion equation , with random Gaussian initial condition, has been calculated exactly using a method known as selective averaging. The probability that the value of the field at a specified spatial coordinate remains positive throughout for a certain time t behaves as for asymptotically large time t. The value of , calculated here for any integer dimension d, is for and 1 otherwise. This exact theoretical result is being reported possibly for the first time and is not in agreement with the accepted values for respectively.

Keywords

Non-Equilibrium Statistical Mechanics, Diffusion Equation, Persistence Exponent, Probability Theory

1. Introduction

The problem in the present paper is to find the persistence exponent for the simple diffusion equation ${\varphi }_t\left(x,t\right)=\Delta \varphi \left(x,t\right)$ . The diffusion equation is an equation that has no stochasticity. In the present problem, the stochasticity is introduced through the random initial conditions. The problem is about evaluating the probability of a certain event. The event is that $\varphi$ at a specified location remains positive throughout the time evolution till a certain time t i.e. the $\varphi$ at the location does not change sign even once. This probability for asymptotically large time is characterised by an exponent ${\theta }_o$ called the persistence exponent. Persistence exponent for the diffusion equation has been a subject of interest to physicists [1] - [7] etc., researchers in mathematics and statistics [8] [9] etc. as well as experimentalists [10] . The interest in the persistence exponent is just not confined to the diffusion equation but to other areas of non-equilibrium physics. Among them random walk [11] , walk in a random environment with or without bias [12] , surface growth [13] , diffusing particle in a random potential with a small concentration of absorbers [14] , behaviour of financial markets [15] etc. are worth mentioning. There are few exact calculations for the persistence exponent in the literature. The case of a simple random walk in one dimension

gives the exponent ${\theta }_o=\frac{1}{2}$ . Even the calculation of persistence exponents for Gaussian processes may not be straight forward.

We revisit the problem of simple diffusion. It is strongly non-Markovian in nature. The problem involves the partial differential equation ${\varphi }_t=\Delta \varphi$ with random Gaussian initial conditions. It appears to remain an unsolved problem even though results [1] [2] and several others have been reported. The problem of diffusion may require a better understanding in the context of persistence. The article tries to find an exact solution to the problem.

2. Simple Diffusion Equation, Random Initial Conditions and Persistence Exponent

The diffusion equation ${\varphi }_t=\Delta \varphi$ is a coarse grained differential equation whose solution is uniquely determined by the initial condition. In the present problem, the initial condition is not fixed but is chosen from a distribution. The initial value of $\varphi$ at every coordinate is chosen from a Gaussian distribution with mean 0, variance k and the initial values of $\varphi$ at any two coordinates are statistically independent.

In order to calculate persistence exponent ${\theta }_o$ we have to calculate the probability that the field $\varphi$ at a specified coordinate does not flip sign even once throughout a time t. This probability ${\mathcal{P}}^{+}\left(t\right)$ of $\varphi$ always remaining +ve behaves in the limit of asymptotically large time as ${\mathcal{P}}^{+}\left(t\right)\sim t^{-{\theta }_o}$ . This is true for a non-stationary process like in the present case. In this article any position x coordinate is a vector quantity in a d dimensional space. The moments of the initial condition distribution described above are given by

$\langle \varphi \left(x\mathrm{,0}\right)\rangle =0$ (1-a)

$\langle \varphi \left(x_1\mathrm{,0}\right)\varphi \left(x_2\mathrm{,0}\right)\rangle =k{\delta }^{\left(d\right)}\left(x_1-x_2\right)$ (1-b)

where k is the variance of the distribution. The solution for the diffusion equation may be written in terms of the initial condition as

$\varphi \left(x\mathrm{,}t\right)={\displaystyle \int {\text{d}}^d{x}^{\prime }G\left(x-x^{\prime }\mathrm{,}t\right)\varphi \left(x^{\prime }\mathrm{,0}\right)}$ (2)

where $G\left(x\mathrm{,}t\right)={\left(4\pi t\right)}^{-d/2}\exp \left(-x^2/4t\right)$ . The plan for the evaluation of the exponent is as follows. First, we have to calculate the probability of $\varphi$ attaining a specific final value $\beta$ at a certain $x=x_o$ starting from a definite initial value $\alpha$ of $\varphi$ at $x=x_o$ . In order to evaluate it we use the method of selective averaging. The paths that take the initial $\alpha$ to the final value $\beta$ also comprise those where $\varphi \left(x_o\right)$ flips sign at least once during time evolution. The probability of such paths is to be subtracted out. Finally, there has to be an integration over the final $\beta$ from 0 to $\infty$ , followed by an integration over $\alpha$ from 0 to $\infty$ .

Selective averaging means averaging over the initial field $\varphi \left(x\mathrm{,0}\right)$ , except when $x=x_o$ . In other words, the averaging is done over all the initial configurations such that $\varphi$ at $x=x_o$ is kept fixed at $\alpha$ (say) i.e. $\varphi \left(x_o,0\right)=\alpha$ while for $x\ne x_o$ $\varphi$ varies according to Gaussian distribution. In this paper the selective distribution, denoted by subscript s, is characterized by the moments,

${\langle \varphi \left(x\mathrm{,0}\right)\rangle }_s=\alpha {\delta }^{\left(d\right)}\left(x-x_o\right)$ (3-a)

${\langle \varphi \left(x_1\mathrm{,0}\right)\varphi \left(x_2\mathrm{,0}\right)\rangle }_s=\left\{k+\left[{\alpha }^2-k\right]{\delta }^{\left(d\right)}\left(x_1-x_o\right)\right\}{\delta }^{\left(d\right)}\left(x_1-x_2\right)$ (3-b)

It may be verified from (3-a), (3-b) that if $x\ne x_o$ , $x_1\ne x_o$ , $x_2\ne x_o$ , we get (1-a), (1-b) and for $x=x_1=x_2=x_o$ , (3-a), (3-b) give $\alpha$ , ${\alpha }^2$ as expected. Using (3-a), (3-b), we can calculate the moments of the random variable $\varphi \left(x_o\mathrm{,}t\right)$ ,

${\langle \varphi \left(x_o\mathrm{,}t\right)\rangle }_s={\left(4\pi t\right)}^{-d/2}\alpha$ (4)

$\begin{array}{c}{\langle {\varphi }^2\left(x_o\mathrm{,}t\right)\rangle }_s={\displaystyle \int {\text{d}}^d{x^{\prime }}_1{\text{d}}^d{x^{\prime }}_2{\left(4\pi t\right)}^{-d}\exp \left[-\frac{{\left(x_o-{x^{\prime }}_1\right)}^2}{4t}\right]}\\ \times \exp \left[-\frac{{\left(x_o-{x^{\prime }}_2\right)}^2}{4t}\right]{\langle \varphi \left({x^{\prime }}_1\mathrm{,0}\right)\varphi \left({x^{\prime }}_2\mathrm{,0}\right)\rangle }_s\\ =k{\displaystyle \int {\text{d}}^d{x^{\prime }}_1{\left(4\pi t\right)}^{-d}\exp \left[-\frac{{\left(x_o-{x^{\prime }}_1\right)}^2}{2t}\right]}-\frac{k}{{\left(4\pi t\right)}^d}+\frac{{\alpha }^2}{{\left(4\pi t\right)}^d}\end{array}$ (5)

While evaluating the second order moment, we have used the relation in (3-b). Hence the mean and the variance of the distribution for $\varphi \left(x_o\mathrm{,}t\right)$ , represented by $\mu$ and ${\sigma }^2$ respectively, are

$\mu ={\langle \varphi \left(x_o\mathrm{,}t\right)\rangle }_s={\left(4\pi t\right)}^{-d/2}\alpha$ (6-a)

$\begin{array}{c}{\sigma }^2={\langle {\varphi }^2\left(x_o\mathrm{,}t\right)\rangle }_s-{\langle \varphi \left(x_o\mathrm{,}t\right)\rangle }_s^2\\ =k{\left(4\pi \right)}^{-d}{2}^{\left(d/2-1\right)}{k}_d\Gamma \left(d/2\right)t^{-d/2}-k{\left(4\pi t\right)}^{-d}\end{array}$ (6-b)

In the above equation $k_d$ denotes the angular integration in d dimensional space while $\Gamma$ represents the usual Gamma function. It may be mentioned that $\varphi \left(x\mathrm{,}t\right)$ in (2) is Gaussian irrespective of whether $\varphi \left(x^{\prime }\mathrm{,0}\right)$ , the initial Gaussian field, is correlated or not. In the present case, though, the initial field is uncorrelated and $\varphi \left(x\mathrm{,}t\right)$ can be proved to be Gaussian using characteristic functions in probability theory [16] . It may be noted that the $\delta$ function distribution is the limiting case of a Gaussian distribution. The expression for the conditional probability for starting at $\alpha$ and being between $\beta$ and $\beta +\text{d}\beta$ at time $t_1$ is

$P\left(\beta \mathrm{<mo>\; |\; </mo>}\alpha \right)\text{d}\beta =\frac{1}{\sqrt{2\pi }\sigma }\exp \left[\frac{-{\left(\beta -\mu \right)}^2}{2{\sigma }^2}\right]\text{d}\beta$ (7)

where $\mu =\mu \left(\alpha \mathrm{,}{t}_1\right)$ and $\sigma =\sigma \left(t_1\right)$ . This probability considers all the paths that start from $\alpha$ to be between $\beta$ and $\beta +\text{d}\beta$ at time $t_1$ including ones that flip en route $\beta$ as depicted in Figure 1. Figure 1 is the projection of the trajectory of the system in the infinite dimensional $\Phi -t$ space on to the $\varphi \left(x_o\right)-t$ plane.

$A\left(\mathrm{0,}\alpha \right)$ represents the starting point and $B\left(t_1\mathrm{,}\beta \right)$ , the destination. AB represents a path along which $\varphi \left(x_o\right)$ does not flip and ADB (blue curve) is a typical path along which $\varphi \left(x_o\right)$ flips. Such paths have to be excluded. The probability of reaching from A to the neighborhood B at asymptotically large time $t_1$ without flipping is given by,

$P^{+}\left(\beta \mathrm{<mo>\; |\; </mo>}\alpha \right)\text{d}\beta =P\left(\beta \mathrm{<mo>\; |\; </mo>}\alpha \right)\text{d}\beta -P\left(\beta \mathrm{<mo>\; |\; </mo>}-\alpha \right)\left(1+O\left(t^{-1}\right)\right)\text{d}\beta$ (8)

The second term represents the probability of paths such as $\bar{A}DB$ originating from $\bar{A}\left(\mathrm{0,}-\alpha \right)$ and terminating in the neighborhood of B at $t_1$ . (8) is not be confused with the method of images in [17] . (8) follows a very different logic in the present case and holds good asymptotically. To prove (8) we will show that there is a one to one mapping from a path $A\to B$ to a path $\bar{A}\to B$ and that the probability of two such paths converges asymptotically. This part is explained in 1) in what follows. Further, to justify (8) we have to show that the “number” of paths $A\to B$ that flip and the “number” of paths $\bar{A}\to B$ converge asymptomatically. This is done in 2). In the subsequent analysis we will consider a d dimensional lattice—lattice spacing being infinitesimally small— instead of continuum for the sake of notational convenience only. The reason for (8) follows.

Figure 1. Projection of $\Phi -t$ trajectory onto the $\varphi \left(x_o\right)-t$ plane.

1) An initial configuration at A of Figure 1 given by $X_{AB}=\left\{\cdots ,{\alpha }_1,\alpha ,{\alpha }_2,\cdots \right\}$ , is considered, where ${\alpha }_1,{\alpha }_2,\cdots$ are the initial values of $\varphi$ at coordinates

$x\ne x_o$ . The corresponding path takes initial $\varphi \left(x_o\right)=\alpha$ to B, then it may be concluded from (2) that $X_{\bar{A}B}=\left\{\cdots ,f{\alpha }_1\mathrm{,}-\alpha \mathrm{,}f{\alpha }_2\mathrm{,}\cdots \right\}$ ( $f=\frac{\beta +{\left(4\pi t\right)}^{-d/2}\alpha }{\beta -{\left(4\pi t\right)}^{-d/2}\alpha }$ ) is a

configuration at $\bar{A}$ which takes initial $\varphi \left(x_o\right)=-\alpha$ to B. Hence there is a one to one mapping of paths from $A\to B$ to those from $\bar{A}\to B$ . It may be underlined here that $f\to 1$ as $t\to \infty$ . This implies that the probability of the two paths approach each other asymptotically.

2) In this part we will address the fact that in the asymptotically large time limit it is a very good approximation to say that there is a one to one correspondence between the paths from A that flip to those from $\bar{A}\to B$ . This may be used as it is a controlled approximation for it improves with increasing t. In order to see this point let us consider a point $C\left(t_2\mathrm{,}\beta \right)$ (not shown in the Figure 1) where $t_2>t_1$ . Let $Y_{AB}=\left\{\cdots ,{\gamma }_1,\alpha ,{\gamma }_2,\cdots \right\}$ be the initial configuration corresponding to path ADB ( the path in blue in Figure 1) where ${\gamma }_1,{\gamma }_2,\cdots$ are the initial values of $\varphi$ at coordinates $x\ne x_o$ . This path crosses zero while reaching B. It can be

shown that $Y_{AC}=\left\{\cdots ,f_1{\gamma }_1,\alpha ,f_1{\gamma }_2,\cdots \right\}$ ( $f_1={\left(\frac{t_2}{t_1}\right)}^{d/2}\frac{\beta -{\left(4\pi t_2\right)}^{-d/2}\alpha }{\beta -{\left(4\pi t_1\right)}^{-d/2}\alpha }$ ) is the

corresponding initial configuration for a path $A\to C$ . The exact expression for $f_1$ contains a coordinate dependent term whose leading order behavior for large t is 1. Since $t_2>t_1$ , we have $f_1>1$ for sufficiently large $t_1$ . Let the time coordinate at D be $t_D$ , then $\varphi \left(x_o,t_D\right)=0$ for the path ADB. Then one may arrive from (2) that $\varphi \left(x_o,t_D\right)<0$ for the initial configuration $Y_{AC}$ . Hence, one can conclude that the path corresponding to $Y_{AC}$ must have flipped at an earlier time than $t_D$ . Therefore, if a path from $A\to B$ flips, the corresponding path from $A\to C$ flips at an earlier time. Since $t_2>t_1$ , the “number” of paths flipping while going from $A\to C$ is more than those from $A\to B$ . Thus the “number” of paths from $A\to B$ that flip is a fraction $f_2$ of those from $\bar{A}\to B$ where $f_2=1-O\left(t_1^{-a}\right)$ for large $t_1$ , a being some positive number.

On account of 1), 2) we say that the probability of the paths (like ADB in Figure 1) that flip while reaching B in the large time limit is given by $P\left(\beta \mathrm{<mo>\; |\; </mo>}-\alpha \right)h_{\text{correction}}$ , where $h_{\text{correction}}=1+O\left(t^{-b}\right)$ , $b=1$ , being Taylor expansion in $t^{-1}$ . In principle, the coefficient of $t^{-b}$ may be a function of $\beta$ . When integrating over $\beta$ —as will be done later—the contribution to the integral comes from the vicinity of $\beta =-\mu \sim t^{-d/2}$ which is vanishingly small in the asymptotic limit. Also, $\text{d}\beta \sim \sigma \sim t^{-d/4}$ . The coefficient is Taylor expanded about $\beta =0$ and only the zeroth order term or the term independent of $\beta$ is retained. So the probability of $\varphi \left(x_o\right)$ not changing sign when reaching the neighborhood ( $\text{d}\beta$ ) of B is, for asymtotically large time,

$\begin{array}{l}\underset{t\to \infty }{lim}\left[P\left(\beta \mathrm{<mo>\; |\; </mo>}\alpha \right)\text{d}\beta -P\left(\beta \mathrm{<mo>\; |\; </mo>}-\alpha \right)h_{\text{correction}}\text{d}\beta \right]\\ =P\left(\beta \mathrm{<mo>\; |\; </mo>}\alpha \right)\text{d}\beta -P\left(\beta \mathrm{<mo>\; |\; </mo>}-\alpha \right)\left(1+O\left(t^{-1}\right)\right)\text{d}\beta \end{array}$ (9)

This leads us to (8). The final $\beta$ may have any value as long as it remains positive. The probability of $\varphi \left(x_o\right)$ starting from $\alpha$ and reaching a final positive value without ever changing sign is

${\mathbb{P}}^{+}\left(\alpha \right)={\displaystyle {\int }_0^{\infty }\text{d}\beta P^{+}\left(\beta \mathrm{<mo>\; |\; </mo>}\alpha \right)}$ (10)

We would now calculate (10) for asymptotically large value of t. Under the circumstances the second term on the R.H.S of (6-b) can be neglected. Further $\frac{{\mu }^2}{{\sigma }^2}\sim \alpha t^{-d/2}$ . Hence for $\alpha \ll t^{d/2}$ , ${\mu }^2/{\sigma }^2\ll 1$ . The expression (10) is evaluated using the identity [18]

$\boxed{{\displaystyle {\int }_0^{\infty }\text{d}x\exp \left(\frac{-x^2}{4\beta }-\gamma x\right)}=\sqrt{\pi \beta }\exp \left(\beta {\gamma }^2\right)\left[1-erf\left(\gamma \sqrt{\beta }\right)\right]}$ (11)

Evaluation leads to a sum of two terms—one is proportional to $t^{-d/4}$ and the other is proportional to $t^{-1}$ . Hence we obtain

${\mathbb{P}}^{+}\left(\alpha \right)\sim \alpha t^{-d/4}$ (12)

for $d\le 4$ . In arriving at the above result the asymptotic expansion of “error function” erf has been used for small argument. Finally, the expression for ${\mathcal{P}}^{+}\left(t\right)$ is obtained by integrating $\alpha$ over a Gaussian distribution.

${\mathcal{P}}^{+}\left(t\right)={\displaystyle {\int }_0^{\infty }\text{d}\alpha {\mathbb{P}}^{+}\left(\alpha \right)Q\left(\alpha \right)}$ (13)

where $Q\left(\alpha \right)$ is the Gaussian distribution for initial $\varphi \left(x_o\mathrm{,}o\right)=\alpha$ with variance k as mentioned at the beginning. If $k\ll t^{d/2}$ , it may be concluded from (12) and (13) that ${\mathcal{P}}^{+}\left(t\right)\sim t^{-d/4}$ or $t^{-1}$ depending on whether $d\le 4$ or not. This gives ${\theta }_o=d/4$ or 1.

3. Result and Conclusion

In the previous section, exact calculation has been carried out to determine the probability ${\mathcal{P}}^{+}\left(t\right)$ of the sign of the field $\varphi$ remaining positive throughout an asymtotically large time t. The probability is ${\mathcal{P}}^{+}\left(t\right)\sim t^{-d/4}$ . Hence, the persistence exponent is ${\theta }_o=d/4$ valid for any arbitrary integer dimension $d\le 4$ . The exponents for $d=1,2,3$ are 0.25, 0.50, 0.75 respectively.

The result may be experimentally verified for a system initially at thermal equilibrium defined by a temperature T. The equilibrium is then disturbed in a suitable manner. The time evolution of the coarse grained temperature at any point satisfies the simple diffusion equation; hence, this time evolution can be studied to find the persistence exponent.

The answer for the exponent ${\theta }_o$ obtained in this paper is in disagreement with all the papers cited in the beginning. The first results for the persistence exponent in the case of the diffusion problem were published in [1] [2] back to back. The papers used a two time correlation function and explicitly applied the approximation (IIA), Independent Interval Approximation, to get to the answer. Application of the two time correlation function is not suitable here and so is IIA which is a Markovian approximation. Further, the papers use Monte Carlo simulation to confirm the result. Monte Carlo method appears to be unsuitable for this problem. Hence all the papers that reproduce the results of [1] [2] are not expected to give the correct answer. In [3] , the authors have defined a correlation function $C\left(T\right)$ , just like [1] [2] , to carry the calculation forward. Let us also consider [6] where the authors have used Kac Polynomials [19] to obtain the “exact exponent” in 2d. The answer obtained agrees perfectly with [1] [2] . In the course of the calculation, they have used that the zero crossing property is

governed by the covariance $c\left(T\right)=\text{sech}\left(\frac{T}{2}\right)$ [6] of the stationary Gaussian

process i.e. the diffusion equation with time redefined. Similarly in [4] , correlator in time $F_{\epsilon }\left(\tau -{\tau }^{\prime }\right)$ has been used in the calculation. The point is that the covariance/correlator/correlation function is a misleading quantity for the problem for reasons mentioned below. The model presented in the paper has randomness only in the initial condition. Once the system starts evolving, there is no further randomness. It evolves in accordance with the kernel in (2). It is encoded in the initial condition when and where the $\varphi$ will flip. The probability of each path is uniquely determined by the probability of initial condition, hence the problem with covariance/correlator. The covariance function imposes stochasticity on the present problem throughout the entire time evolution. We now have a different model with the same correlation function but no unique dependence of the probability of the path on initial condition. It also makes the problem Markovian. Hence, all the previous results are in perfect agreement though the calculated exponent will be different from the actual value. The value of the exponent does not depend on only correlation function, but it depends on other details of the model too. Further, there even appears to be experimental proof [10] for the results of [1] [2] . The experimental setup of [10] does not represent the diffusion model described in this paper. The setup satisfies the approximations of the previous papers and hence the agreement with their result.

Acknowledgements

The author would like to thank CSIR, India for Fellowship during the course of the work (2004) at IACS, India.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Majumdar, S.N., Sire, C., Bray, A.J. and Cornell, S.J. (1996) Physical Review Letters, 77, 2867-2870. https://doi.org/10.1103/PhysRevLett.77.2867
[2]	Derrida, B., Hakim, V. and Zeitak, R. (1996) Physical Review Letters, 77, 2871-2874. https://doi.org/10.1103/PhysRevLett.77.2871
[3]	Ehrhardt, G.C.M.A. and Bray, A.J. (2002) Physical Review Letters, 88, Article ID: 070601. https://doi.org/10.1103/PhysRevLett.88.070601
[4]	Hilhorst, H.J. (2000) Physica A: Statistical Mechanics and Its Applications, 277, 124-126. https://doi.org/10.1016/S0378-4371(99)00509-9
[5]	Schehr, G. and Majumdar, S.N. (2008) Journal of Statistical Physics, 132, 235-273. https://doi.org/10.1007/s10955-008-9574-3
[6]	Poplavskyi, M. and Schehr, G. (2018) Physical Review Letters, 121, Article ID: 150601. https://doi.org/10.1103/PhysRevLett.121.150601
[7]	Barbier-Chebbah, A., Benichou, O. and Voituriez, R. (2020) Physical Review E, 102, Article ID: 062115. https://doi.org/10.1103/PhysRevE.102.062115
[8]	Aurzada, F. and Simon, T. (2015) Persistence Probabilities and Exponents. In: Lévy Matters V, Springer, Berlin, 183-224. https://doi.org/10.1007/978-3-319-23138-9_3
[9]	Dembo, A. and Mukherjee, S. (2015) The Annals of Probability, 43, 85-118. https://doi.org/10.1214/13-AOP852
[10]	Wong, G.P., Mair, R.W., Walsworth, R.L. and Cory, D.G. (2001) Physical Review Letters, 86, 4156-4159. https://doi.org/10.1103/PhysRevLett.86.4156
[11]	Schwarz, J.M. and Maimon, R. (2001) Physical Review E, 64, Article ID: 016120. https://doi.org/10.1103/PhysRevE.64.016120
[12]	Le Doussal, P., Monthus, C. and Fisher, D.S. (1999) Physical Review E, 59, 4795-4840. https://doi.org/10.1103/PhysRevE.59.4795
[13]	Krug, J., Kallabis, H., Majumdar, S.N., Cornell, S.J., Bray, A.J. and Sire, C. (1997) Physical Review E, 56, 2702-2712. https://doi.org/10.1103/PhysRevE.56.2702
[14]	Le Doussal, P. (2009) Journal of Statistical Mechanics: Theory and Experiment, No. 7, P07032. https://doi.org/10.1088/1742-5468/2009/07/P07032
[15]	Constantin, M. and Das Sarma, S. (2005) Physical Review E, 72, Article ID: 051106. https://doi.org/10.1103/PhysRevE.72.051106
[16]	Mathews, J. and Walker, R.L. (1970) Mathematical Methods of Physics. Volume 501, WA Benjamin, New York.
[17]	Chadrasekhar, S. (1989) Stochastic, Statistical and Hydromagnetic Problems in Physics and Astronomy. Selected Papers Vol. 3.
[18]	Gradshteyn, I.S. and Ryzhik, I.M. (2014) Table of Integrals, Series, and Products. Academic Press, Cambridge.
[19]	Kac, M. (1943) Bulletin of the American Mathematical Society, 49, 314-320. https://doi.org/10.1090/S0002-9904-1943-07912-8