Searching for a shell

Green's function (part I)

Math  DFT 

Green’s function method is a way to decompose a complex inhomogenous source using Dirac delta functions and then combine individually solutions to obtain the true answer.

Preconditions

Green’s function is a solution to a linear differential equation with a Dirac delta function as the inhomogeneous source and the additional requirement of the boundary condition being homogenous.

To clarify, by homogenous/inhomogeneous, we mean $\textit{something}=0$ . For example, a differential equation is said to be inhomogeneous if $g(x)$ in the following equation does not equals to $0$, homogeneous otherwise:

D^f(x)=g(x)\hat{D} f(x) = g(x)

Similarly, the boundary condition is said to be homogeneous only if it specify the value (or the value of the derivative) of the function on the boundary to be $0$. i.e.

f(x=L)=0 or f′(x=L)=0f(x=L) = 0 \text{ or } f'(x=L) = 0

On the other hand, by linear, we mean no coupling between operators with different derivative order. i.e. no product of function and its derivatives, such as $f(x) \cdot f’(x)$.

Overall, the differential equation should look like:

a1(x)f(x)+a2(x)f′(x)+⋯+an(x)fn(x)=g(x)=δ(x−x′),(1)a_1(x)f(x) + a_2(x)f'(x) + \cdots + a_n(x)f^{n}(x) = g(x) = \delta (x-x'), \tag{1}

or, we can put all operation into an operator and write Eq. 1 as:

D^f(x)=g(x)=δ(x−x′)(2)\hat{D} f(x) = g(x) = \delta(x-x') \tag{2}

where $\hat{D} = a_1(x)+a_2(x)\frac{d}{dx} + \cdots + a_n(x)(\frac{d}{dx})^n$. Since there are no coupling between different order of differentiation, we call this operator a linear operator. The linearity is essential to the superposition theory that we are so used to in quantum mechanics that we’ve just taken it for granted.

Deriving of Green’s function method

Now, how do we apply this particular solution of a particular differential equation to solve all those general questions? More importantly, what kinds of general differential equations can be solved using this method?

Like I said in the beginning, the essence of this method is that it breaks down the complex inhomogeneous source using Dirac delta functions. This means that we can have any continuous complex inhomogeneous source in our general differential equations, but we still need the operator to be linear. i.e.:

D^f(x)=g(x)(3)\hat{D} f(x) = g(x) \tag{3}

where, $g(x)$ is a general continuous function.

$g(x)$ can be expressed (decomposed) using a infinite set of Dirac delta functions with weights as:

g(x)=∫−∞∞δ(x′−x)g(x′)dx′,(4)g(x) = \int_{-\infty}^{\infty} \delta (x'-x) g(x') dx' \tag{4},

where, each Dirac delta function (more appropriately, unit pulse) takes one value from $g(x’)$. Or this can be interpreted as each unit pulse (located at $x$) has a weight of $g(x)$.

By plugging Eq. 4 into Eq. 3, we get:

D^f(x)=g(x)=∫−∞∞δ(x′−x)g(x′)dx′.(5)\hat{D} f(x) = g(x) = \int_{-\infty}^{\infty} \delta (x'-x) g(x') dx'. \tag{5}

Note that Green’s function is the solution to the equation:

D^G(x,x′)=δ(x−x′).(6)\begin{aligned} \hat{D} G(x,x') &= \delta (x-x') \tag{6}\\ \end{aligned}.

Plug Eq. 6 to Eq. 5 (note that we’ve swapped $x$ and $x’$), we get:

D^f(x)=∫−∞∞D^G(x,x′)g(x′)dx′,(7)\hat{D} f(x) = \int_{-\infty}^{\infty} \hat{D} G(x,x') g(x') dx', \tag{7}

and since $\hat{D}$ is a linear operator, the final solution to Eq. 3 is:

f(x)=∫−∞∞G(x,x′)g(x′)dx′.(8)f(x) = \int_{-\infty}^{\infty}G(x,x') g(x') dx'. \tag{8}

Note that to get the particular solution, we also need to add the general solution of homogeneous equation:

D^f(x)=0\hat{D} f(x) = 0

to the solution in Eq. 8.

Now, things are much clearer: The green’s function is to be viewed as the building block of the particular solution since it decomposes the inhomogeneous source using Dirac functions.

Practical example

To get a taste of how the whole process works, let’s consider solving the Poisson equation with Green’s function method! The Poisson equation relates the electrical potential to the charge density through:

∇2ψ(r⃗)=−ρ(r⃗)ϵ0(9)\nabla^2 \psi(\vec r) = - \frac{\rho(\vec r)}{\epsilon_0} \tag{9}

in a non-periodic system, we have the boundary condition:

ψ(r⃗→∞)=0\psi(\vec r \rightarrow \infty) = 0

substituting $g(x)$ in Eq. 8 with $-\frac{\rho(\vec r)}{\epsilon_0}$, we get:

ψ(r⃗)=∫G(r⃗,r⃗′)ρ(r⃗)dr⃗′,(10)\psi(\vec r) = \int G(\vec r, \vec r') \rho(\vec r) d \vec r', \tag{10}

where $G(\vec r, \vec r’)$ satisfies:

∇2G(r⃗,r⃗′)=−1ϵ0δ(r⃗−r⃗′)(11)\nabla^2 G(\vec r, \vec r') = - \frac{1}{\epsilon_0} \delta (\vec r - \vec r') \tag{11}

To solve this equation and get the Greens function needed to solve Eq. 9, first, we set $\vec r’ = 0$ since the delta function only cares about the difference between $\vec r$ and $\vec r’$. This yields $G(\vec r)$ and $\delta(\vec r)$ Then we Fourier transform the Green’s function and the delta function:

  • Green’s function:

    G(k⃗)=F[G(r⃗)](k⃗)=1(2π)3/2∫G(r⃗)e−ik⃗r⃗dr⃗G(r⃗)=F−1[G(k⃗)](r⃗)=1(2π)3/2∫G(k⃗)eik⃗r⃗dk⃗\begin{aligned} G(\vec k) = \mathcal{F}[G(\vec r)](\vec k) &= \frac{1}{(2\pi)^{3/2}} \int G(\vec r) e^{-i\vec k \vec r} d\vec r\\ G(\vec r) = \mathcal{F}^{-1}[G(\vec k)](\vec r) &= \frac{1}{(2\pi)^{3/2}} \int G(\vec k) e^{i\vec k \vec r} d\vec k \end{aligned}

  • delta function:

    δ(k⃗)=F[δ(r⃗)](k⃗)=1(2π)3/2∫δ(r⃗)e−ik⃗r⃗dr⃗δ(r⃗)=F−1[δ(k⃗)](r⃗)=1(2π)3/2∫δ(k⃗)eik⃗r⃗dk⃗\begin{aligned} \delta(\vec k) = \mathcal{F}[\delta(\vec r)](\vec k) &= \frac{1}{(2\pi)^{3/2}} \int \delta(\vec r) e^{-i\vec k \vec r} d\vec r\\ \delta(\vec r) = \mathcal{F}^{-1}[\delta(\vec k)](\vec r) &= \frac{1}{(2\pi)^{3/2}} \int \delta(\vec k) e^{i\vec k \vec r} d\vec k \end{aligned}

Now, Fourier transforming both sides of Eq. 11 gives:

∇2G(r⃗,r⃗′)=−1ϵ0δ(r⃗−r⃗′)∇21(2π)3/2∫G(k⃗)eik⃗r⃗dk⃗=−1ϵ01(2π)3/2∫δ(k⃗)eik⃗r⃗dk⃗1(2π)3/2∫−k2G(k⃗)eik⃗r⃗dk⃗=1(2π)3/2∫−1ϵ0δ(k⃗)eik⃗r⃗dk⃗,(12)\begin{aligned} \nabla^2 G(\vec r, \vec r') &= - \frac{1}{\epsilon_0} \delta (\vec r - \vec r')\\ \nabla^2 \frac{1}{(2\pi)^{3/2}} \int G(\vec k) e^{i\vec k \vec r} d\vec k &= - \frac{1}{\epsilon_0} \frac{1}{(2\pi)^{3/2}} \int \delta(\vec k) e^{i\vec k \vec r} d\vec k\\ \frac{1}{(2\pi)^{3/2}} \int -k^2 G(\vec k) e^{i\vec k \vec r} d\vec k &= \frac{1}{(2\pi)^{3/2}} \int - \frac{1}{\epsilon_0} \delta(\vec k) e^{i\vec k \vec r} d\vec k\\ \end{aligned}, \tag{12}

and we can observe that:

k2G(k⃗)=−δ(k⃗)ϵ0G(k⃗)=−1k2δ(k⃗)ϵ0(13)\begin{aligned} k^2 G(\vec k) &= - \frac{\delta(\vec k)}{\epsilon_0} \\ G(\vec k) &= - \frac{1}{k^2} \frac{\delta(\vec k)}{\epsilon_0} \tag{13} \end{aligned}

To get back to the real-space, inverse FT $G(\vec k)$ and use the inverse FT of $\delta(\vec k)$:

G(r⃗)=1(2π)3/2∫G(k⃗)eik⃗r⃗dk⃗=1(2π)3/2∫−1k2δ(k⃗)ϵ0eik⃗r⃗dk⃗=1(2π)3/2∫−1k21ϵ0(1(2π)3/2∫δ(r⃗′)e−ik⃗r⃗′dr⃗′)eik⃗r⃗dk⃗=1(2π)3∫eik⃗(r⃗−r⃗′)k2dk⃗∫δ(r⃗′)ϵ0dr⃗′=1ϵ01(2π)3∫eik⃗(r⃗)k2dk⃗,\begin{aligned} G(\vec r) &= \frac{1}{(2\pi)^{3/2}} \int G(\vec k) e^{i\vec k \vec r} d\vec k\\ &= \frac{1}{(2\pi)^{3/2}} \int - \frac{1}{k^2} \frac{\delta(\vec k)}{\epsilon_0} e^{i\vec k \vec r} d\vec k\\ &= \frac{1}{(2\pi)^{3/2}} \int - \frac{1}{k^2} \frac{1}{\epsilon_0} \left ( \frac{1}{(2\pi)^{3/2}} \int \delta(\vec r') e^{-i\vec k \vec r'} d\vec r' \right ) e^{i\vec k \vec r} d\vec k\\ &= \frac{1}{(2\pi)^{3}} \int \frac{e^{i\vec k (\vec r - \vec r')}}{k^2} d\vec k \int \frac{\delta(\vec r')}{\epsilon_0} d \vec r'\\ &= \frac{1}{\epsilon_0} \frac{1}{(2\pi)^{3}} \int \frac{e^{i\vec k (\vec r)}}{k^2} d\vec k \end{aligned},

assuming spherical symmetry:;

G(r⃗)=1ϵ01(2π)3∫eik⃗(r⃗)k2dk⃗=1ϵ01(2π)3∫02π∫0π∫0∞k2sin⁡(θ)eikcos⁡(θ)∣r⃗∣k2dϕdθdk=1ϵ01(2π)2∫0π∫0∞k2sin⁡(θ)eikcos⁡(θ)∣r⃗∣k2dθdk,\begin{aligned} G(\vec r) &= \frac{1}{\epsilon_0} \frac{1}{(2\pi)^{3}} \int \frac{e^{i\vec k (\vec r)}}{k^2} d\vec k \\ &= \frac{1}{\epsilon_0} \frac{1}{(2\pi)^{3}} \int_0^{2\pi} \int_{0}^{\pi} \int_{0}^{\infty} k^2 \sin(\theta) \frac{e^{i k \cos(\theta) |\vec r|}}{k^2} d\phi d\theta dk \\ &= \frac{1}{\epsilon_0} \frac{1}{(2\pi)^{2}} \int_{0}^{\pi} \int_{0}^{\infty} k^2 \sin(\theta) \frac{e^{i k \cos(\theta) |\vec r|}}{k^2} d\theta dk \end{aligned},

substituting $\cos(\theta)$ with $u$ and change the limits of the angular integral:

G(r⃗)=1ϵ01(2π)2∫0π∫0∞k2sin⁡(θ)eikcos⁡(θ)∣r⃗∣k2dθdk=1ϵ01(2π)2∫−11∫0∞eiku∣r⃗∣dudk=1ϵ01(2π)2∫0∞1ik∣r⃗∣(eik∣r⃗∣−e−ik∣r⃗∣)dk.\begin{aligned} G(\vec r) &= \frac{1}{\epsilon_0} \frac{1}{(2\pi)^{2}} \int_{0}^{\pi} \int_{0}^{\infty} k^2 \sin(\theta) \frac{e^{i k \cos(\theta) |\vec r|}}{k^2} d\theta dk \\ &= \frac{1}{\epsilon_0} \frac{1}{(2\pi)^{2}} \int_{-1}^{1} \int_{0}^{\infty} e^{i k u |\vec r|} du dk \\ &= \frac{1}{\epsilon_0} \frac{1}{(2\pi)^{2}} \int_{0}^{\infty} \frac{1}{ik |\vec r|} \left( e^{i k |\vec r|} - e^{ -i k |\vec r|} \right) dk \end{aligned}.

We know that:

sin⁡(a)=eia−e−ia2i\sin(a) = \frac{e^{ia} - e^{-ia}}{2i}

so that:

G(r⃗)=1ϵ01(2π)2∫0∞1ik∣r⃗∣(eik∣r⃗∣−e−ik∣r⃗∣)dk=1ϵ012π2∫0∞1k∣r⃗∣sin⁡(k∣r⃗∣)dk\begin{aligned} G(\vec r) &= \frac{1}{\epsilon_0} \frac{1}{(2\pi)^{2}} \int_{0}^{\infty} \frac{1}{ik |\vec r|} \left( e^{i k |\vec r|} - e^{ -i k |\vec r|} \right) dk \\ &= \frac{1}{\epsilon_0} \frac{1}{2\pi^{2}} \int_{0}^{\infty} \frac{1}{k |\vec r|} \sin (k |\vec r|) dk \end{aligned}

substituting $v$ with $(k \mid \vec r \mid)$:

G(r⃗)=1ϵ012π2∫0∞1k∣r⃗∣sin⁡(k∣r⃗∣)dk=1ϵ012π21∣r⃗∣∫0∞sin⁡(v)vdv=14πϵ0∣r⃗∣(14)\begin{aligned} G(\vec r) &= \frac{1}{\epsilon_0} \frac{1}{2\pi^{2}} \int_{0}^{\infty} \frac{1}{k |\vec r|} \sin (k |\vec r|) dk\\ &= \frac{1}{\epsilon_0} \frac{1}{2\pi^{2}} \frac{1}{|\vec r|} \int_{0}^{\infty} \frac{\sin(v)}{v} dv \\ &= \frac{1}{4 \pi \epsilon_0 |\vec r|} \end{aligned} \tag{14}

which, is exactly the potential generated by a single point charge located at the origin point.

If we don’t specify $\vec r’ = 0$, then:

G(r⃗,r⃗′)=14πϵ0∣r⃗−r⃗′∣(15)\begin{aligned} G(\vec r, \vec r') &= \frac{1}{4 \pi \epsilon_0 |\vec r - \vec r'|} \end{aligned} \tag{15}

Finally, substituting Eq. 14 to Eq. 10, we get:

ψ(r⃗)=14πϵ0∫ρ(r⃗′)∣r⃗−r⃗′∣dr⃗′(16)\psi(\vec r) = \frac{1}{4 \pi \epsilon_0} \int \frac{\rho(\vec r')}{|\vec r - \vec r'|} d \vec r' \tag{16}

Which, is exactly the solution to our Poisson equation.

For more information:

  • This article if based on :link: this paper

  • Watch :link: this video from Andrew Dotson for a more “matrix” like explanation 😉.

Next up, :link: the quantum version.

Author | Chengcheng Xiao

Currently a postdoctoral research associate at Imperial College London. Predicting electron behaviour since 2016.