5. Green’s function of Laplacian (ongoing)

v2: Added section 3.

Some prerequisites: vector calculus, Heaviside step function, Dirac delta function

1. Introduction

In electrostatics, the presence of electric charge influences the value of electric field through Gauss’ law

\displaystyle  \vec{{\nabla}}\cdot\vec{E} = \frac{{\rho}}{{\epsilon}_0}, \ \ \ \ \ (1)

where {\vec{E}} is electric field, {{\rho}} is charge density, and {{\epsilon}_0} is permittivity of free space. It is often convenient to introduce a scalar potential {{\phi}} via

\displaystyle  \vec{E} = -\vec{{\nabla}} {\phi}. \ \ \ \ \ (2)

The Gauss’ law then becomes

\displaystyle  {\nabla}^2{\phi} = -\frac{{\rho}}{{\epsilon}_0}. \ \ \ \ \ (3)

In order to solve this equation, one considers Green’s function {G(\vec{r})} such that

\displaystyle  {\nabla}^2 G(\vec{r}) = {\delta}^{(3)}(\vec{r}), \ \ \ \ \ (4)

where {\vec{r} = (x,y,z).} A solution [In fact, in order to solve for Green’s function, one has to also specify boundary conditions. So in our case, it is possible that there are also other solutions.] is

\displaystyle  G(\vec{r}) = -\frac{1}{4{\pi} r}, \ \ \ \ \ (5)

where {r = \sqrt{x^2 + y^2 + z^2}.} This then allows us to obtain the scalar potential

\displaystyle  \begin{array}{rl} {\phi}(\vec{r}) &\displaystyle=\int d^3\vec{r'} G(\vec{r}-\vec{r'})\frac{{\rho}(\vec{r'})}{{\epsilon}_0}\\ &\\ &\displaystyle=-\frac{1}{4{\pi}{\epsilon}_0}\int d^3\vec{r'} \frac{{\rho}(\vec{r'})}{|\vec{r}-\vec{r'}|}, \end{array} \ \ \ \ \ (6)

whose form looks very familiar to any physicist.

What I want to discuss in this post is not on electrostatics, but on Green’s function of Laplacian operator. I will quote the Green’s functions and discuss a few alternative ways to show that they are indeed Green’s functions of the Laplacian operator. The main purpose is just to do it for fun, not necessarily following any mainstream presentation. I might even discuss some heuristic approaches.

2. Green’s function of 1D Laplacian

In one dimension, let us consider

\displaystyle  \frac{d^2}{dx^2}G(x) = {\delta}(x). \ \ \ \ \ (7)

A solution to this equation is simply

\displaystyle  G(x) = \frac{|x|}{2}. \ \ \ \ \ (8)

Simply by looking at this Green’s function, one can be convinced because it is not a smooth function. While it can be differentiated twice, the second derivative is infinite at {x=0,} and hence is not well-defined. This behaviour seems to get along well with Dirac’s delta function on the RHS of eq.(7).

Let us discuss a few ways to verify that eq.(8) satisfies eq.(7).

2.1. Using Heaviside step function

The behaviour of a discontinuous function can be well captured by Heaviside step function, which is defined as

\displaystyle  {\Theta}(x) = \begin{cases} 1;&x\geq 0\\ 0;&x< 0 \end{cases} . \ \ \ \ \ (9)

Its plot looks like a step:

The Heaviside step function can be used in order to simplify the writing of the definition of a piecewise function. For example, suppose we have a function

\displaystyle  f(x) = \begin{cases} h(x);&x\geq 0\\ g(x);&x< 0 \end{cases} . \ \ \ \ \ (10)

This function can be written using only one line as

\displaystyle  f(x) = h(x){\Theta}(x) + g(x)(1-{\Theta}(x)). \ \ \ \ \ (11)

So we can write

\displaystyle  |x| = \begin{cases} x;&x\geq 0\\ -x;&x< 0 \end{cases} \ \ \ \ \ (12)

as

\displaystyle  |x| = x{\Theta}(x) - x(1-{\Theta}(x)). \ \ \ \ \ (13)

The derivative of Heaviside step function is given by Dirac delta function:

\displaystyle  \frac{d{\Theta}(x)}{d x} = {\delta}(x). \ \ \ \ \ (14)

Therefore,

\displaystyle  \begin{array}{rl} \displaystyle\frac{d|x|}{dx} &\displaystyle={\Theta}(x) - (1-{\Theta}(x))+2x{\delta}(x)\\ &\displaystyle={\Theta}(x) - (1-{\Theta}(x))\\ &\displaystyle=2{\Theta}(x) - 1, \end{array} \ \ \ \ \ (15)

\displaystyle  \begin{array}{rl} \displaystyle\frac{d^2|x|}{dx^2} =2{\delta}(x). \end{array} \ \ \ \ \ (16)

So from this result, we can easily see that eq.(7) is indeed satisfied.

2.2. Replacing {|x|} by the limit of smooth function sequence

The presence of Dirac delta function might not be so friendly for many, especially mathematicians. So instead of dealing directly with Dirac delta function, one can also consider it as a sequence of functions. To make it even more convenient, one can demand that functions in the sequence are all smooth. So in this way, since the RHS of eq.(7) is replaced by a sequence of smooth functions, one has to also replace the LHS of eq.(7) in the similar manner. This amounts to replacing {|x|} by a sequence of smooth functions.

A simple choice is given by

\displaystyle  \{\sqrt{x^2 + {\epsilon}}| {\epsilon}>0\}. \ \ \ \ \ (17)

Clearly, the limit as {{\epsilon}\rightarrow 0} of this sequence is {|x|.} So

\displaystyle  \frac{d^2\sqrt{x^2 + {\epsilon}}}{dx^2} = \frac{{\epsilon}}{(x^2 + {\epsilon})^{3/2}}. \ \ \ \ \ (18)

It can easily be seen that

\displaystyle  \lim_{{\epsilon}\rightarrow 0}\frac{d^2\sqrt{x^2 + {\epsilon}}}{dx^2} = \begin{cases} 0; & x\neq 0\\ \infty; & x=0 \end{cases} . \ \ \ \ \ (19)

Next, let us consider

\displaystyle  \begin{array}{rl} \displaystyle\int_{-\infty}^\infty dx\frac{{\epsilon}}{(x^2 + {\epsilon})^{3/2}} &\displaystyle=\sqrt{{\epsilon}}\int_{-{\pi}/2}^{{\pi}/2}d{\theta} \sec^2{\theta} \frac{{\epsilon}}{{\epsilon}^{3/2}\sec^3{\theta}}\\ &\\ &\displaystyle=\int_{-{\pi}/2}^{{\pi}/2}d{\theta} \cos{\theta}\\ &\\ &\displaystyle=\sin{\theta}|_{-{\pi}/2}^{{\pi}/2}\\ &\\ &\displaystyle=2, \end{array} \ \ \ \ \ (20)

where we have made the substitution {x = \sqrt{{\epsilon}}\tan{\theta}.} By combining the results eq.(19)(20), we see that

\displaystyle  \lim_{{\epsilon}\rightarrow 0}\frac{d^2}{dx^2}\left(\frac{\sqrt{x^2 + {\epsilon}}}{2}\right) ={\delta}(x). \ \ \ \ \ (21)

By moving {\lim_{{\epsilon}\rightarrow 0}} into the parenthesis, the expression in the parenthesis then becomes {|x|/2,} and hence we recover eq.(7).

Note that during the derivation, we have not been careful with moving around the limit {{\epsilon}\rightarrow 0.} In general, the order of the operations matter, especially when infinity is involved. In fact, we should verify that the operations we have made are allowed, or else we should have kept the order of the operations at all times. Well, we can consider ourselves lucky as physicists as we are leaving all these worries to mathematicians… [Actually, I am not yet able to tackle these issues nor see the need for doing this any soon. So readers interested in learning these delicate matters are advised to look elsewhere.]

2.3. Replacing {|x|} by the limit of smooth function sequence (II)

Let us consider an alternative sequence which represents {|x|:}

\displaystyle  \left\{\frac{x^2}{\sqrt{x^2+{\epsilon}}}\ |\ {\epsilon}>0\right\}. \ \ \ \ \ (22)

For convenience, call

\displaystyle  f_{\epsilon}(x)\equiv \frac{x^2}{\sqrt{x^2+{\epsilon}}}. \ \ \ \ \ (23)

Compute

\displaystyle  \frac{df_{\epsilon}(x)}{dx} = \frac{{\left(x^{2} + 2 \, {\epsilon}\right)} x}{{\left(x^{2} + {\epsilon}\right)}^{\frac{3}{2}}}, \ \ \ \ \ (24)

\displaystyle  \frac{d^2f_{\epsilon}(x)}{dx^2} = -\frac{{\left(x^{2} - 2 \, {\epsilon}\right)} {\epsilon}}{{\left(x^{2} + {\epsilon}\right)}^{\frac{5}{2}}}. \ \ \ \ \ (25)

So

\displaystyle  \lim_{{\epsilon}\rightarrow 0}\frac{d^2 f_{\epsilon}(x)}{dx^2} = \begin{cases} 0; & x\neq 0\\ \infty; & x=0 \end{cases} . \ \ \ \ \ (26)

Next, let us consider

\displaystyle  \begin{array}{rl} \displaystyle\int_{-\infty}^\infty dx \frac{d^2f_{\epsilon}(x)}{dx^2} &\displaystyle=\frac{df_{\epsilon}(x)}{dx}\bigg|_{-\infty}^\infty\\ &\\ &\displaystyle=2. \end{array} \ \ \ \ \ (27)

From the results eq.(26)(27), we see that

\displaystyle  \lim_{{\epsilon}\rightarrow 0}\frac{d^2}{dx^2}\left(\frac{f_{\epsilon}(x)}{2}\right) = {\delta}(x). \ \ \ \ \ (28)

But the the limit {{\epsilon}\rightarrow 0} of the quantity in the bracket on the LHS is just {|x|/2.} So we recover eq.(7).

2.4. Replacing {|x|} by the limit of smooth function sequence (III)

The alternative possibilities seem endless. As yet another example, let us consider a sequence

\displaystyle  \left\{f_k(x)\equiv x\tanh kx\ |\ k>0\right\}. \ \ \ \ \ (29)

It reaches {|x|} as {k\rightarrow\infty.} That is

\displaystyle  \lim_{k\rightarrow\infty}f_k(x) = |x|. \ \ \ \ \ (30)

The figure below illustrates that this is really the case.

I leave the verification that this case satisfies eq.(7) as an exercise for interested readers. But to ease off the calculations, I give results of derivatives:

\displaystyle  \frac{d f_k(x)}{dx} = \frac{k x + \cosh\left(k x\right) \sinh\left(k x\right)}{\cosh^2\left(k x\right)}, \ \ \ \ \ (31)

\displaystyle  \frac{d^2 f_k(x)}{dx^2} = -\frac{2k \, {\left(k x \sinh\left(k x\right) - \cosh\left(k x\right)\right)} }{\cosh^3\left(k x\right)}. \ \ \ \ \ (32)

3. Green’s function for 2D Laplacian

Green’s function for 2D Laplacian is given by

\displaystyle  \left(\frac{{\partial}^2}{{\partial} x^2}+\frac{{\partial}^2}{{\partial} y^2}\right) G(x,y) ={\delta}(x,y), \ \ \ \ \ (33)

where {{\delta}(x,y)} is 2D Dirac’s delta function satisfying

\displaystyle  {\delta}(x,y) = \begin{cases} 0,& (x,y)\neq (0,0),\\ \infty,& (x,y) = (0,0), \end{cases} \ \ \ \ \ (34)

and

\displaystyle  \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} dx dy{\delta}(x,y) =1. \ \ \ \ \ (35)

Let us simply write down the Green’s function, which is given by

\displaystyle  G(x,y) = \frac{1}{2{\pi}}\log r, \ \ \ \ \ (36)

where {r = \sqrt{x^2 + y^2}.} A contour plot is given below. The darker the colour the higher the value of Green’s function.

Again, let us show some alternative ways that Green’s function eq.(36) satisfies eq.(33).

3.1. Using divergence theorem

It is clear that Green’s function eq.(36) is well-defined whenever {r\neq 0.} So by excluding the point {r = 0,} it should satisfy

\displaystyle  \left(\frac{{\partial}^2}{{\partial} x^2}+\frac{{\partial}^2}{{\partial} y^2}\right) G(x,y) =0,\qquad r\neq 0. \ \ \ \ \ (37)

This can easily be shown to be the case. For this we make use of the intermediate results

\displaystyle  \frac{{\partial}}{{\partial} x}G = \frac{x}{2{\pi} r^2}, \ \ \ \ \ (38)

\displaystyle  \frac{{\partial}^2}{{\partial} x^2}G = \frac{1}{2{\pi} r^2} - \frac{x^2}{{\pi} r^4}, \ \ \ \ \ (39)

and similar ones for derivatives wrt {y.}

The analysis however is invalid for {r=0.} For this case, let us use divergence theorem. Consider a closed surface {S} which encloses the volume {V.} Divergence theorem states that

\displaystyle  \int_V dV \vec{{\nabla}}\cdot\vec{f} =\int_S dS \hat{n}\cdot\vec{f}, \ \ \ \ \ (40)

where {\hat{n}} is a unit vector pointing out from each point on the surface.

In 2d, a closed surface is in fact a closed loop, whereas the enclosed volume is in fact the region inside the closed loop. The idea is to consider a small loop enclosing the origin. In particular, let us consider a circle of radius {{\epsilon}} centred at {r = 0.} So {\hat n = \hat r.} A sketch is given by the figure below.

Let us consider an integral

\displaystyle  \begin{array}{rl} \displaystyle\lim_{{\epsilon}\rightarrow 0}\int_V dV {\nabla}^2 G &=\displaystyle\lim_{{\epsilon}\rightarrow 0}\int_V dV \vec{{\nabla}}\cdot\vec{{\nabla}}G\\ &\\ &=\displaystyle\lim_{{\epsilon}\rightarrow 0}\int_S dS \hat{n}\cdot\vec{{\nabla}}G. \end{array} \ \ \ \ \ (41)

Let us compute

\displaystyle  \begin{array}{rl} \displaystyle\vec{{\nabla}} G\bigg|_{r={\epsilon}} &=\displaystyle \frac{\hat r}{2{\pi} r}\bigg|_{r={\epsilon}}\\ &\\ &=\displaystyle \frac{\hat r}{2{\pi}{\epsilon}}. \end{array} \ \ \ \ \ (42)

Therefore,

\displaystyle  \hat n\cdot\vec{{\nabla}} G\bigg|_{r={\epsilon}} = \frac{1}{2{\pi}{\epsilon}}, \ \ \ \ \ (43)

and hence

\displaystyle  \begin{array}{rl} \displaystyle\lim_{{\epsilon}\rightarrow 0}\int_V dV {\nabla}^2 G &=\displaystyle\lim_{{\epsilon}\rightarrow 0}\frac{1}{2{\pi}{\epsilon}}\int_S dS\\ &\\ &=\displaystyle\lim_{{\epsilon}\rightarrow 0} 1\\ &=1. \end{array} \ \ \ \ \ (44)

By integrating eq.(33) over the area enclosed by a small circle centred at {r=0,} the LHS agrees with the LHS of eq.(44). By using the property of Dirac’s delta function, the RHS from the two equations can also be shown to agree. So we have finished showing that eq.(36) satisfies eq.(33).

3.2. Working in polar coordinates with cut-off

NB: This method is still not fully justified. Proceed with care.

In polar coordinates, eq.(33) becomes

\displaystyle  {\nabla}^2 G = \left(\frac{{\partial}^2}{{\partial} r^2} + \frac{1}{r}\frac{{\partial}}{{\partial} r} + \frac{1}{r^2}\frac{{\partial}^2}{{\partial}{\theta}^2}\right)G =\frac{1}{2{\pi} r}{\delta}(r). \ \ \ \ \ (45)

Let us cut the solution eq.(36) off for {r<{\epsilon}.} That is, we consider

\displaystyle  G_{\epsilon} = {\Theta}(r-{\epsilon})\frac{1}{2{\pi}}\log r. \ \ \ \ \ (46)

Compute

\displaystyle  \frac{{\partial} G_{\epsilon}}{{\partial} r} =\frac{1}{2{\pi} r}{\Theta}(r-{\epsilon}) + {\delta}(r-{\epsilon})\frac{1}{2{\pi}}\log r, \ \ \ \ \ (47)

\displaystyle  \frac{{\partial}^2 G_{\epsilon}}{{\partial} r^2} =-\frac{1}{2{\pi} r^2}{\Theta}(r-{\epsilon}) + {\delta}(r-{\epsilon})\frac{1}{{\pi} r}+ {\delta}'(r-{\epsilon})\frac{1}{2{\pi}}\log r. \ \ \ \ \ (48)

Therefore,

\displaystyle  {\nabla}^2 G_{\epsilon} = {\delta}(r-{\epsilon})\frac{1}{{\pi} r} + {\delta}'(r-{\epsilon})\frac{1}{2{\pi}}\log r + {\delta}(r-{\epsilon})\frac{1}{2{\pi} r}\log r. \ \ \ \ \ (49)

So by taking the limit {{\epsilon}\rightarrow 0} gives

\displaystyle  {\nabla}^2 G = {\delta}(r)\frac{1}{{\pi} r} + {\delta}'(r)\frac{1}{2{\pi}}\log r + {\delta}(r)\frac{1}{2{\pi} r}\log r. \ \ \ \ \ (50)

However, the RHS of this equation looks different to that of eq.(45). In order to justify this, let us multiply RHS of eq.(45) by {2{\pi} r f(r)} and integrate from {r=0} to {r=\infty,} where {f(r)} is a well-behaved function. This gives

\displaystyle  \int dr 2{\pi} r {\delta}(r) f(r)\frac{1}{2{\pi} r} = f(0). \ \ \ \ \ (51)

On the other hand, the same operation on RHS of eq.(50) gives

\displaystyle  \begin{array}{rl} \displaystyle\int dr 2{\pi} r &\displaystyle f(r) \left({\delta}(r)\frac{1}{{\pi} r} + {\delta}'(r)\frac{1}{2{\pi}}\log r + {\delta}(r)\frac{1}{2{\pi} r}\log r\right)\\ &\\ &=\displaystyle 2f(0) + \int dr r\log r\ f(r) {\delta}'(r) + \int dr {\delta}(r) f(r)\log r\\ &\\ &=\displaystyle 2f(0) - \int dr \left(\log r f(r) + f(r) + r\log r f'(r)\right){\delta}(r) + \int dr {\delta}(r) f(r)\log r\\ &\\ &=\displaystyle f(0) - \int dr\ r\log r f'(r){\delta}(r)\\ &\\ &=\displaystyle f(0) - \lim_{r\rightarrow 0} r\log r f'(r)\\ &\\ &=f(0), \end{array} \ \ \ \ \ (52)

as required. Note that the second term on the penultimate step vanishes because {r\log r\rightarrow 0} as {r\rightarrow 0.}

So we have again shown that eq.(36) satisfies eq.(33). However, the analysis I have shown might still contain some (if not lots of) flaws. So please read with care.

4. Green’s function for 3D Laplacian

[To be continued…]