Diffusion Through a Half Space : Equivalence Between Different Formulations of the Unique Solution

Diffusion through a half space involves a classical parabolic partial differential equation that is encountered in many fields of physics and has significant engineering applications, concerning particularly heat and mass transfer. However, in the specialized literature, the solution is usually achieved restricting the problem to particular cases and attaining apparently different formulations, thus a comprehensive overview is hindered. In this paper, the solution of the diffusion equation in a half space with a boundary condition of the first kind is worked out by means of the Fourier’s Transform, the Green’s function and the similarity variable, with a proof of equivalence – not found elsewhere – of these different approaches. The keystone of the proof rests on the square completion method applied to Gaussian-like integrals, widely used in Quantum Field Theory.


Introduction
One of the most important mathematical methods of Quantum Field Theory is square completion to compute Gaussian integrals that arise from the Path Integral Approach pioneered by Feynman [1 -2].As Zee [3] says, "Believe it or not, a significant fraction of the theoretical physics literature consists of performing variations and elaborations of this basic Gaussian integral".Although all books at the advanced undergraduate and most books at the graduate level use the method of canonical quantization (which avoids Gaussian integrals) or defer path integrals to the last chapters, the book by Zee introduces path integrals from the beginning.
The purpose of this paper is to show how the square completion method to compute Gaussian-like integrals allows understanding the equivalence of apparently very different formulations of the solution of the standard parabolic PDE encountered in heat conduction and other diffusion problems that play an important role in many Engineering applications.A paper by Slutsky [4] applied the full-blown machinery of path integrals to diffusion in the context of polymer physics.A similar though shorter treatment of linear polymer molecules as random walks is found in earlier works such as the books by Schulman [5] and Carrà [6].On the other hand, Hall [7] very recently discussed the connection between random walks and path integrals.Our purpose is somewhat more limited and, at the same time, more accessible to a broader audience.We want to show that an important integral, which lies at the core of the path integral approach to Quantum Field Theory, emerges naturally from the juxtaposition of classical methods to solve the diffusion PDE and highlights the hidden connections among them.

Problem Statement
We consider an important class of partial differential equations in the general form belonging to the category of parabolic equations.They are used to represent in different contexts a kind of transport referred to as diffusion [8].For instance, setting  T u  , temperature [K] and , thermal diffusivity [m 2 s -1 ], Eq. ( 1) describes heat conduction in a homogeneous isotropic continuum with constant properties and without heat sources.This equation was first derived by Fourier [9]. A c u  , molar concentration of the chemical species A [mol m -3 ] and classical problem in Metallurgy is the estimate of the decarburization depth in steel [11].A variant of the diffusion PDE is also used in Nuclear Reactor Physics to model neutron density in the one-speed approximation [12].
, probability density function of the velocity of a particle and D , diffusion coefficient.This is known as the Fokker-Planck equation with zero drift coefficient [13].It is interesting to note that this equation has been recently applied beyond Physics to study the volatility in financial markets [14].A classical problem is the determination of , initially at a uniform value 0 u with the interface subjected to a first kind boundary condition for 0  t (Dirichelet's problem).The differential problem is stated as is a solution as well) and 0 0  t conventionally.

General Solution by Means of Fourier Analysis
The problem is approached by the Fourier analysis [15].At first, we consider the Fourier transform of Eq. (2c) with respect to time where  is the angular frequency [rad s -1 ].On the other hand,   Notice that Eqs.(3) and ( 4) have the same coefficient   2 1 2   .However, different choices are possible.The reader is referred to Appendix A for a brief discussion on this subject.
The variable separation method is applied to Eq. (2a) looking for particular solutions in the form as suggested by the integrand in Eq. ( 4).Replacing in Eq. (2a) and dividing both members by The characteristic equation associated to Eq. ( 6) is giving which is often called wavenumber [16] [rad m -1 ].Hence, the general solution of Eq. ( 6) turns out to be Replacing in Eq. ( 5) which is often called a thermal wave even though the second-order derivative with respect to time, characteristic of the wave equation, does not appear in Eq. (1).A thorough discussion about the concept of wave and thermal waves is given by Salazar [17].
The boundary conditions Eqs. (2c) and (2d) are applied to calculate the coefficients Eq. ( 13) is applied under the assumption that D is real and positive.This requirement is actually a consequence of the second principle of thermodynamics.On the other hand, if D were imaginary Eq. ( 1) would turn into the well-known Schrödinger equation which does not describe diffusive transport, hence is not treated here.
Eq. ( 10) becomes Passing to the limit, Eq. ( 13), as u must be finite Repeating the same procedure for , Eq. ( 15) reduces to a constant that can be neglected as Eq.(2a) only contains derivatives of u.
A unique representation is obtained introducing the sign function, strictly related to the Heaviside step function as it will be shown in Section 5: By integration over the angular frequency it is obtained which, for 0  x , reduces to Eq. ( 4).Finally, to eliminate the transformed function U it is convenient to express the boundary condition Eq. (2c) from Eq. (3) as follows Handbooks usually report particular cases of Eq. ( 20).The reader should address in particular the book of Prestini [18] where the presented approach is developed in a less general way but with very interesting practical applications.Restricting the attention to the heat conduction problem, many authors deal with the cases of constant and periodic heating [19 -26] though the general problem is not discussed in detail.

The Similarity Solution
Most of the cited bibliography directly refer to the similarity solution of (2a).Actually, dimensional analysis shows that the dependence on x and t is condensed in the combinations The former is sometimes called the Boltzmann number, whereas the latter is simply known as the similarity variable.The physical meaning of B is discussed in Appendix B.
Generally speaking, similarity solutions are only a subset of the existing solutions.In this case, however, it can be shown that all the solutions are self-similar.
Adopting the Boltzmann number, the following identities hold Hence, replacing in Eq. (2a), an ordinary differential equation is obtained: (26) which is written as a first order equation setting and performing a second integration taking into account Eq. ( 21) The integral in Eq. ( 29) cannot be evaluated as a combination of elementary functions.It is a transcendental function as it is shown by the Liouville's theory [27].
It is customary to define the error function (31) It is then obtained from Eq. ( 29) where multiplicative factors and additive constants have been lumped in 1 C and 2 C respectively.A particularly useful case study is obtained if Eq. ( 2c) is written as This is the response of the half-space to a step variation of u on its interface.

Figures 1a and 1b report  
2  u and   t x u , , respectively, to clarify the meaning of the term similarity.It is evident that each spatial distribution of u at a certain time instant is self-similar because, when reported in terms of the similarity variable  , all the distributions collapse into a unique curve.It is also evident that the diffusive transport described by (1) occurs instantaneously in contrast with the basic tenets of Special Relativity.Actually, since , the propagation speed turns out to be infinite, meaning that the effect of a perturbation at the interface 0  x is immediately felt at any distance from the interface.This is a theoretical problem arising from the constitutive equations relating the diffusive flux to the gradient of u , such as Fourier's law and Fick's first law.However, this effect is quite small in the most common situations and it is usually neglected [28].
To show that all the solutions of the problem defined by Eqs.(2a) to (2d) are self-similar, it is convenient to switch to a dimensionless formulation.Recalling the definition of the Boltzmann number, Eq. ( 21), characteristic time Replacing in Eqs.(2a) to (2d) the dimensionless problem results , the latter showing that a double infinity of solutions is derived by choosing arbitrarily At this point it is interesting to seek a general solution of problem defined by Eqs.(2a) to (2d) in the form of infinite series of particular solutions like Eq. (33) where self-similarity is evident rather than Eq.(20).In the following it is discussed the method for representing the new form of the solution and the equivalence with Eq. (20).

Integral Representations of the Dirac's  Function and Heaviside's Step Applied to the Diffusion Problem
The Dirac's delta function is defined as [29]  Accordingly, an important property is that any function There are different representations [30] of such a function that today mathematicians prefer to call more properly a distribution.A useful one for the purpose of this paper is found in Mandl [31]  The following developments justify this choice.Actually, if Eq. ( 3) is replaced in Eq. ( 4), it yields On the other hand, according to Eq. ( 38) Comparing Eq. ( 40) and Eq. ( 41) it is seen that which is equivalent to Eq. ( 39).
The Heaviside's step function is defined as [3]  The relation between H and the sign function used in Eq. ( 20) is formally expressed as   The application of the step function usually requires, as seen for the Dirac's delta, suitable representations.For the purpose of this paper, it is convenient to use the following [32]: Considering Eq. ( 39) it can be shown that which is more easily understood if H is thought as the limit of a ramp that rises from 0 to 1 about 0  t .This relation is useful to transform Eq. (38) in another useful representation of any continuous function.Integrating by parts The geometrical meaning of Eq. ( 47) is a representation of as the superposition of elementary steps of height df (Duhamel's formula) [20].Equation ( 47) is then applied to the solution of problem defined by Eqs.(2a) to (2d) as follows: where   0 , 0   u since no perturbation is applied at the interface before the initial time.In any case, the value of    , 0 u would be only an additive constant.
The same holds for    , 0 u since any physical perturbation has finite duration.The solution is then built as a continuous linear combination of the particular solution, Eq. ( 33), corresponding to the response of the half space to a constant perturbation at the interface, that is It is worthwhile noting that the argument of the error function is prevented from assuming imaginary values when . Equation ( 49) clearly shows that the diffusion process obeys to the principle of delayed causality.Actually, a generic boundary condition at the interface is decomposed as the superposition of elementary steps, the response to which is delayed by the time interval t t   .Some attempts to modify the Fourier's law in order to prevent an infinite speed of propagation of the perturbations, as observed in the previous section, happened to violate the delayed causality principle [33].

Solution by Means of the Green's Function
The solution of Eq. (2a) in the same form as Eq. ( 20) is also obtained by means of the Green's function [34].Since the derivation is less straightforward than making use of the Fourier's transform, only an outline will be given in the following.For this purpose, it is convenient to consider the inhomogeneous diffusion equation: where   physically represents a source term.
The Green's function   is the solution of Eq. ( 50) if under appropriate initial and boundary conditions, that is The relation between u and G is found as The Green's function is determined as follows.
At first,   Then the partial derivatives that appear in Eq. (52) result respectively Replacing Eqs.(55), ( 56) and (57) in Eq. (52) the Fourier's back-transform of the Green's function results Hence, from Eq. ( 54) which is more conveniently rewritten as with Integration is performed by means of the residuals theorem and the Jordan's lemma in the complex plane along the loop depicted in Figure 2. The integration loop is suitably selected so that, when , only the contribution along the diameter (which extends to the whole real axis) is different from zero, as a consequence of the Jordan's lemma.
The integrand poles are: The residual is then Hence, the residuals theorem allows writing Finally, replacing Eq. (65) in Eq. ( 60) The partial derivative of Eq. (66) with respect to x coincides with Eq. ( 20) if the boundary condition is set as Actually, replacing Eq. ( 68) in Eq. ( 20) it follows Despite the different sing in the square brackets, Eqs.(67) and (69) describe the same scalar field: it is sufficient to replace simultaneously  with   and t with t  .
In summary, it has been shown that the solution of Eq. (2a) with the boundary condition Eq. ( 68) is equivalent to the solution of Eq. ( 52).In other words, through the Green's function the solution of the homogeneous diffusion equation is derived from the inhomogeneous one endowed with a suitable source term.

Equivalence of the General Solutions by Means of the Square Completion Method
Apparently, the two general solutions developed in Sections 3 and 5, respectively, are quite different.In particular, Eq. ( 20) involves two improper integrals whereas only one appears in Eq. (49); Eq. ( 20) involves complex functions, whereas Eq. ( 49) is restricted to the real field; Eq. ( 20) involves the boundary condition   t u  , 0 , whereas its derivative with respect to time appears in Eq. ( 49).Nevertheless, the boundary value problems involving Eq. ( 1) do have a unique solution [19] so that Eq. ( 20) and Eq. ( 49) must be equivalent.The proof of equivalence is given in the following.
Equation ( 20) is written as Replacing in Eq. ( 70) (72) Comparison between Eq. (49) and Eq.(72) shows that their equivalence would imply apart from an integration constant that is equal to zero according to the initial condition Eq. (2b).
From the Leibnitz formula The latter expression is manipulated as first by applying to the argument of the exponential the square completion method [35], widely used in Quantum Field Theory to compute path integrals [5], [36].Incidentally, it is worthwhile mentioning that the first mathematician who studied the so-called Gaussian integrals was not Gauss but De Moivre in 1733 [37].
Replacing Eqs. ( 62) to (64) in Eq. (80) yields, finally The equivalence of Eq. ( 20) and Eq. ( 49) is then proven thanks to the decisive resort to the square completion method, just as in many gaussian-like integrals found in the path integral approach to Quantum Field Theory.This result is written in a more interesting form as This heuristic argument is also useful to understand the microscopic origin of irreversibility: at each interaction of the particle, any information about the previous one is lost, hence diffusion cannot be inverted.A recent paper by Brazzle [41] discusses a pedagogical method to simulate diffusion by means of spreadsheet computation, which is widely available.Actual random walks are built on a variety of lattices.However the easiest method remains the one followed by Gautreau [40].

Figure 1 -
Figure 1 -Response of the half-space to a step variation of u on its interface in terms of the similarity variable (a) and of the natural variables (b).

u
as a characteristic value of u, the following set of dimensionless quantities is identified .e. the set of all the solutions is split into two equivalence classes.As each class includes only self-similar solutions, all the possible solutions are self-similar.

Figure 2 -
Figure 2 -Integration loop and poles of the integrand for Eq.(61).Black dots represent the poles for 0   hand term shows the same form of the Boltzmann number Eq.