Green's Functions for Partial Differential Equations

A linear partial differential equation $Lu=f$ on a domain $\Omega$ , with prescribed boundary or initial data, is at heart a request to invert the operator $L$ . The Green's function is precisely that inverse written as an integral kernel: it is the response felt at a field point $\mathbf{x}$ when a single unit impulse is placed at a source point $\mathbf{x}'$ . Because $L$ is linear, an arbitrary source $f$ is nothing but a continuous crowd of point sources, and the field it produces is the very same crowd of point responses added up. That one idea threads electrostatics, heat flow, wave propagation, and quantum scattering onto a single string. The construction then splits cleanly by the character of the operator. For elliptic operators such as $-\nabla^2$ (Poisson) or $-\nabla^2-k^2$ (Helmholtz), $G$ is an equilibrium potential or a steady scattered wave, pinned down by spatial boundary conditions alone. For parabolic operators such as $\partial_t-D\nabla^2$ (diffusion), $G$ is the heat kernel — a Gaussian that smooths instantly and carries information at formally infinite speed but with exponentially small tails. For hyperbolic operators such as $\partial_t^2-c^2\nabla^2$ (waves), $G$ is a propagator supported on or inside the light cone, so that causality is written directly into where the kernel is allowed to be nonzero. Two technical pillars make all of this operational. Green's second identity trades a volume PDE for a representation in boundary data and the kernel; the spectral resolution $G=\sum_n\phi_n(\mathbf{x})\phi_n^*(\mathbf{x}')/\lambda_n$ exposes $G$ as $L^{-1}$ acting mode by mode through the spectrum, while the $i\epsilon$ prescription on the poles decides retarded from advanced. In practice one reads $G$ off a Fourier transform in free space and then bends it to the boundary with images, conformal maps, or eigenfunction sums — and the far field of those integrals is exactly where Green's functions hand off to the asymptotic methods this chapter builds toward.

Key takeaways

A Green's function $G(\mathbf{x},\mathbf{x}')$ is the response of a linear operator $L$ to a unit point source, $L_xG=\delta(\mathbf{x}-\mathbf{x}')$ , so every solution of $Lu=f$ is assembled by superposition: $u(\mathbf{x})=\int G(\mathbf{x},\mathbf{x}')f(\mathbf{x}')\,\mathrm{d}\mathbf{x}'$ .
Boundary and initial conditions are baked into $G$ : the same operator $-\nabla^2$ yields different kernels for Dirichlet, Neumann, or causal (retarded) requirements — there is no single "the" Green's function until these are fixed.
Operator class sets the character: elliptic (Laplace, Helmholtz) gives static and scattering kernels, parabolic (diffusion) gives the Gaussian heat kernel, hyperbolic (wave) gives light-cone propagators.
Green's second identity converts the volume PDE into a boundary integral, exposing $G$ and its normal derivative on $\partial\Omega$ ; the spectral sum builds $G$ from eigenfunctions divided by eigenvalues.
Dimension dictates the singularity: $-\tfrac{1}{2\pi}\ln r$ in 2D, $1/(4\pi r)$ in 3D, a light-cone $\delta(t-r/c)$ in 1+3D — and this controls how the superposition integral converges.
Key idea: a Green's function is the kernel of the inverse operator $L^{-1}$ subject to specified homogeneous boundary or initial conditions — solve $L_xG(\mathbf{x},\mathbf{x}')=\delta(\mathbf{x}-\mathbf{x}')$ once, and every inhomogeneous problem $Lu=f$ collapses to a single integral plus boundary terms.

Where this lives today

Green's functions are not a historical curiosity; they are the working currency of modern computational and theoretical physics. In electronic-structure theory the many-body Green's function is the central object of the GW approximation and of dynamical mean-field theory (DMFT), where the one-particle propagator's poles are the quasiparticle energies measured in photoemission; computing it accurately for strongly correlated materials — cuprates, heavy-fermion compounds — remains a frontier open problem, and the non-equilibrium (NEGF) version predicts current through a single molecule. Engineers solve Laplace, Helmholtz, and elastodynamic problems with boundary element methods that discretize only the surface, using exactly the free-space kernels above, and antenna designers work with the dyadic Green's function of Maxwell's equations. In gravitational-wave astronomy the templates LIGO matches against data are retarded-potential solutions of the linearized Einstein equations, and in seismology the Earth's Green's function — the ground motion from a point slip — is now reconstructed from ambient-noise correlations, imaging the planet without waiting for earthquakes. Wherever a linear response to a localized cause must be computed, the miller's kernel is still doing the work.

The physical picture

Think of the source $f(\mathbf{x})$ as a crowd of independent point charges, one at every location. Each one, of strength $f(\mathbf{x}')\,\mathrm{d}\mathbf{x}'$ , produces its own little field $G(\mathbf{x},\mathbf{x}')$ centred on $\mathbf{x}'$ ; linearity says the total field is just the sum — the integral — of these elementary echoes. The Green's function is the operator's "unit echo," and the boundary conditions reshape that echo so the sum automatically satisfies the data at the walls. The singularity of $G$ where $\mathbf{x}=\mathbf{x}'$ is not a blemish but a requirement: applying $L$ must return something that integrates to one over an arbitrarily small ball, and in three dimensions the unique locally integrable function whose Laplacian is a negative delta is exactly $1/(4\pi r)$ .

Figure

A unit point source and its 1/r echo — the Green's function of minus the Laplacian in three dimensions. A unit impulse at x′ produces the Coulomb field G = 1/(4π r), which falls off as 1/r: 0.159 at r = 0.5 m, 0.0796 at r = 1 m, 0.0398 at r = 2 m, with G carrying units of inverse length

Boundary conditions act like mirrors. The method of images literally places fictitious charges outside the domain so their fields cancel $G$ on the boundary (Dirichlet) or cancel its normal derivative (Neumann) — you arrange the echoes until the wall looks transparent to the equation while still satisfying the data. Time-dependent operators turn the kernel into a propagator: for the wave equation the echo travels outward at speed $c$ , and the response at $\mathbf{x}$ only "hears" a source event after the light-travel delay $|\mathbf{x}-\mathbf{x}'|/c$ . Causality is then not an extra assumption but the literal statement that $G$ vanishes off the light cone. The second figure lines up the three operator classes so the difference — static, spreading, or sharply on the cone — is visible at a glance.

Figure

The operator class dictates the kernel. Elliptic −∇²G = δ gives the static 1/(4π r) potential; parabolic ∂ₜ − D∇² gives a Gaussian heat kernel whose width √(4Dt) spreads as √t; hyperbolic ∂ₜ² − c²∇² gives a sharp pulse δ(t − r/c)/(4π r) living only on the light cone r = ct

The core equations

Everything below is a special case of one sentence: $G$ is what $L$ turns a point source into, and $u$ is the superposition of those responses. The defining equation states it, and every other relation is a way of solving it or reading a boundary condition off it.

L_x\,G(\mathbf{x},\mathbf{x}') = \delta(\mathbf{x}-\mathbf{x}') \qquad (1)

where: $L_x$ — linear differential operator acting on the unprimed variable; $G$ — Green's function, the kernel of $L^{-1}$ ; $\delta$ — Dirac delta in the relevant dimension; $\mathbf{x},\mathbf{x}'$ — field point and source point. Reading: $G(\cdot,\mathbf{x}')$ is the field produced by a unit spike parked at $\mathbf{x}'$ , and it inherits every boundary condition we impose on $L$ . Once $G$ is known, the inhomogeneous problem is solved by a single integral — no series, no relaxation, just a convolution of the source against the kernel.

u(\mathbf{x}) = \int_\Omega G(\mathbf{x},\mathbf{x}')\,f(\mathbf{x}')\,\mathrm{d}\mathbf{x}' \qquad (2)

where: $u$ — solution of $Lu=f$ obeying homogeneous boundary conditions; $f$ — source or inhomogeneity; $\Omega$ — domain. Reading: apply $L$ under the integral, use (1), and the delta sifts out $f(\mathbf{x})$ — so (2) really does solve the PDE. To handle inhomogeneous boundary data one needs Green's second identity, the integration-by-parts statement that moves derivatives between two fields and leaves a boundary flux behind. It is the engine that turns a volume PDE into a boundary representation.

\int_\Omega\!\big(u\,L^\dagger v - v\,Lu\big)\,\mathrm{d}V = \oint_{\partial\Omega}\! Q[u,v]\cdot \mathrm{d}\mathbf{S} \qquad (3)

where: $L^\dagger$ — formal adjoint of $L$ ; $Q[u,v]$ — bilinear concomitant or boundary flux, equal to $u\nabla v-v\nabla u$ for $L=\nabla^2$ ; $\partial\Omega$ — boundary with outward normal. Reading: the volume mismatch between $L$ and its adjoint is entirely a surface term — the root of "boundary data determines the interior." Setting $v=G$ in (3) and using $-\nabla^2G=\delta$ collapses one volume term to $u(\mathbf{x})$ and yields the representation that every boundary-value technique specializes.

u(\mathbf{x}) = \int_\Omega G\,f\,\mathrm{d}V' + \oint_{\partial\Omega}\!\big(G\,\partial_{n'}u - u\,\partial_{n'}G\big)\,\mathrm{d}S' \qquad (4)

where: $\partial_{n'}$ — outward normal derivative at the source point; $G$ — chosen to satisfy $-\nabla^2G=\delta$ plus a boundary condition that kills one of the two surface terms. Reading: choose $G_D|_{\partial\Omega}=0$ and only the known Dirichlet data $u|_{\partial\Omega}$ survives; choose the Neumann kernel and only $\partial_nu|_{\partial\Omega}$ does. When $L$ is self-adjoint with a complete eigenbasis, $G$ has an explicit spectral form — literally $L^{-1}$ written mode by mode, each eigenfunction weighted by the reciprocal of its eigenvalue.

G(\mathbf{x},\mathbf{x}') = \sum_n \frac{\phi_n(\mathbf{x})\,\phi_n^*(\mathbf{x}')}{\lambda_n} \qquad (5)

where: $\phi_n$ — orthonormal eigenfunctions of $L$ obeying the boundary conditions; $\lambda_n$ — eigenvalue with $L\phi_n=\lambda_n\phi_n$ . Reading: a zero eigenvalue ( $\lambda_n=0$ ) makes the sum blow up — that missing mode is exactly the Neumann solvability condition in disguise. For free space, translation invariance makes the Fourier transform diagonalize $L$ : the kernel is the inverse transform of one over the operator's symbol.

\tilde G(\mathbf{k}) = \frac{1}{\sigma(\mathbf{k})}, \qquad G(\mathbf{x}) = \int \frac{\mathrm{d}^d k}{(2\pi)^d}\,\frac{e^{\,i\mathbf{k}\cdot\mathbf{x}}}{\sigma(\mathbf{k})} \qquad (6)

where: $\sigma(\mathbf{k})$ — symbol of $L$ , defined by $L\,e^{i\mathbf{k}\cdot\mathbf{x}}=\sigma(\mathbf{k})\,e^{i\mathbf{k}\cdot\mathbf{x}}$ , so $\sigma=|\mathbf{k}|^2$ for $-\nabla^2$ ; $d$ — spatial dimension. Reading: dividing by $\sigma$ is the free-space $L^{-1}$ ; the only subtlety is what to do when $\sigma$ vanishes on the contour. That subtlety is the whole story for wave operators. The symbol $\sigma=c^2k^2-\omega^2$ has real zeros, and how one steps around them selects the boundary condition in time.

G_{\mathrm{ret/adv}}(\mathbf{x},t) = \int\!\frac{\mathrm{d}\omega}{2\pi}\frac{\mathrm{d}^3k}{(2\pi)^3}\,\frac{e^{\,i(\mathbf{k}\cdot\mathbf{x}-\omega t)}}{c^2k^2-(\omega\pm i\epsilon)^2} \qquad (7)

where: $\pm i\epsilon$ — infinitesimal pole shift; the $+$ choice, closed in the lower half-plane for $t>0$ , gives the retarded kernel supported on the future cone, the $-$ choice gives the advanced kernel on the past cone. Reading: same operator, same source — only the contour differs, and that alone is the difference between cause-before-effect and its time-reverse.

Derivation AThe free-space Coulomb kernel in three dimensions

The goal is to compute the Green's function of $-\nabla^2$ in three dimensions from scratch, showing that the requirement $-\nabla^2G=\delta^3$ has the unique free-space solution $G=1/(4\pi r)$ — the Coulomb/Newton potential — and thereby to see how the Fourier method turns a differential equation into an algebraic division followed by a single radial integral. We transform, divide by the symbol $k^2$ , and invert, doing the angular integral exactly and the radial integral as a standard Dirichlet integral.

Step 1. Fourier transform the equation with the convention $G=\int\!\frac{\mathrm{d}^3k}{(2\pi)^3}e^{i\mathbf{k}\cdot\mathbf{x}}\tilde G(\mathbf{k})$ . Since $-\nabla^2\to k^2$ and $\delta^3\to1$ , the PDE becomes the algebraic statement $k^2\tilde G(\mathbf{k})=1$ , so $\tilde G(\mathbf{k})=1/k^2$ . Why: translation invariance is what lets plane waves diagonalize the operator, converting a differential equation into a one-line division by the symbol.

Step 2. Invert and do the angular part. Writing $r=|\mathbf{x}|$ and choosing the polar axis of $\mathbf{k}$ along $\mathbf{x}$ , the solid-angle integral gives $\int\mathrm{d}\Omega_k\,e^{ikr\cos\theta}=4\pi\sin(kr)/(kr)$ , so $G(r)=\frac{1}{(2\pi)^2}\int_0^\infty\!\mathrm{d}k\,\frac{2\sin(kr)}{kr}=\frac{1}{2\pi^2 r}\int_0^\infty\!\frac{\sin(kr)}{k}\,\mathrm{d}k$ . Why: the isotropy of $1/k^2$ means only the radial magnitude survives, and the angular integral is elementary.

Step 3. Evaluate the radial integral. The Dirichlet integral $\int_0^\infty\frac{\sin(kr)}{k}\,\mathrm{d}k=\frac{\pi}{2}$ for any $r>0$ , so $G(r)=\frac{1}{2\pi^2 r}\cdot\frac{\pi}{2}=\frac{1}{4\pi r}$ . Why: the integral is independent of the magnitude of $r$ (it depends only on the sign), which is exactly what makes the $1/r$ scaling come out clean.

Step 4. Check by direct action. Away from the origin $\nabla^2(1/r)=0$ ; integrating $-\nabla^2G$ over a small ball and using the divergence theorem gives $-\oint\nabla G\cdot\mathrm{d}\mathbf{S}=-\oint(-1/4\pi r^2)\,r^2\,\mathrm{d}\Omega=1$ , so $-\nabla^2G=\delta^3$ as required. Why: the flux through any small sphere is exactly one regardless of its radius — the operational meaning of "unit point source." *The Coulomb kernel $1/(4\pi|\mathbf{x}|)$ is the standard free-space result, giving the radial value $G=0.0796\,\mathrm{m^{-1}}$ at $r=1\,\mathrm{m}$ , consistent with the figure above.*

Go deeper · Derivation B · 1 minThe heat kernel of the diffusion operator

Derivation BThe heat kernel of the diffusion operator

The goal is to build the Green's function of the diffusion operator $\partial_t-D\nabla^2$ in $d$ dimensions and show it is the normalized Gaussian $G=\Theta(t)(4\pi Dt)^{-d/2}\exp(-|\mathbf{x}|^2/4Dt)$ , the object that is simultaneously the temperature response to a heat pulse and the transition probability of a Brownian particle. We transform in space to reduce the PDE to a first-order ODE in time, solve it with the causal step function, and invert with a Gaussian integral that factorizes over the dimensions.

Step 1. Fourier transform in space only. With $-\nabla^2\to k^2$ , the equation $(\partial_t-D\nabla^2)G=\delta(\mathbf{x})\delta(t)$ becomes $\partial_t\tilde G+Dk^2\tilde G=\delta(t)$ , an ordinary differential equation in $t$ at each wavevector $\mathbf{k}$ . Why: space and time decouple once the Laplacian is replaced by a number, leaving a bank of independent one-dimensional relaxation problems.

Step 2. Solve the ODE causally. For $t>0$ the homogeneous solution decays as $e^{-Dk^2t}$ , and the temporal delta forces the jump $\tilde G(0^+)=1$ , so $\tilde G(\mathbf{k},t)=\Theta(t)\,e^{-Dk^2t}$ . Why: the step function $\Theta(t)$ is where causality enters — the kernel is strictly zero before the pulse, and its unit jump encodes the delta.

Step 3. Invert the transform. The remaining integral $G=\Theta(t)\int\frac{\mathrm{d}^dk}{(2\pi)^d}e^{i\mathbf{k}\cdot\mathbf{x}-Dk^2t}$ factorizes into $d$ identical Gaussian integrals; completing the square in each gives $\int\frac{\mathrm{d}k}{2\pi}e^{ikx_j-Dk^2t}=\frac{1}{\sqrt{4\pi Dt}}e^{-x_j^2/4Dt}$ . Why: a Gaussian stays Gaussian under Fourier transform, so the whole $d$ -dimensional inversion is just the product of one-dimensional results.

Step 4. Collect and check normalization. Multiplying the $d$ factors gives $G=\Theta(t)(4\pi Dt)^{-d/2}e^{-|\mathbf{x}|^2/4Dt}$ , and $\int G\,\mathrm{d}^dx=1$ for every $t>0$ . Why: probability (or heat) is conserved — the point source neither creates nor destroys total content, only spreads it, with width $\sqrt{4Dt}$ that grows as $\sqrt{t}$ . *The prefactor $(4\pi Dt)^{-d/2}$ , the Gaussian exponent, and the unit-normalization integral are all standard results, and the mean-square displacement $\langle x^2\rangle=2Dt$ per axis follows immediately from the kernel.*

Go deeper · Derivation C · 1 minThe retarded propagator of the wave operator

Derivation CThe retarded propagator of the wave operator

The goal is to derive the retarded Green's function of the 3D wave operator, $G_\mathrm{ret}=\delta(t-r/c)/(4\pi r)$ , and to make explicit why the signal lands exactly on the light cone $r=ct$ — Huygens' principle in three space dimensions. We Fourier transform in time to land on a Helmholtz equation, pick the outgoing (retarded) solution set by the pole prescription, and transform back, at which point the frequency integral is nothing but a representation of the delta function.

Step 1. Transform in time. With $\partial_t^2\to-\omega^2$ , the equation $(\partial_t^2-c^2\nabla^2)G=\delta^3(\mathbf{x})\delta(t)$ becomes the Helmholtz equation $(\nabla^2+\omega^2/c^2)\tilde G=-\frac{1}{c^2}\delta^3(\mathbf{x})$ with wavenumber $k=\omega/c$ . Why: at fixed frequency the wave operator is elliptic (Helmholtz), so we can reuse the free-space $1/(4\pi r)$ machinery with a phase.

Step 2. Take the outgoing solution. The Helmholtz Green's function is $\tilde G(\mathbf{x},\omega)=e^{i(\omega/c)|\mathbf{x}|}/(4\pi c^2|\mathbf{x}|)$ ; the $+i(\omega/c)r$ sign is the retarded/outgoing choice, equivalent to displacing the poles by $\omega\to\omega+i\epsilon$ . Why: the Sommerfeld radiation condition — energy flowing out, not in — is what selects this branch and, with it, cause before effect.

Step 3. Transform back. Writing $r=|\mathbf{x}|$ , $G_\mathrm{ret}=\frac{1}{4\pi c^2 r}\int\frac{\mathrm{d}\omega}{2\pi}e^{-i\omega(t-r/c)}$ , and the frequency integral is the standard representation $\int\frac{\mathrm{d}\omega}{2\pi}e^{-i\omega(t-r/c)}=\delta(t-r/c)$ . Why: an unweighted integral of a pure phase over all frequencies is a delta in the conjugate variable — the sharpest possible arrival.

Step 4. Read off causality. Absorbing the constants gives $G_\mathrm{ret}=\delta(t-r/c)/(4\pi r)$ : the kernel is nonzero only when $t=r/c$ , so a source event is felt at distance $r$ precisely after the light-travel time. Why: the delta support is the light cone itself, which is Huygens' principle — in odd space dimensions $\ge 3$ a sharp pulse stays sharp, with no lingering tail. *The retarded propagator $\delta(t-r/c)/(4\pi r)$ is the canonical result, and it reduces to the static $1/(4\pi r)$ Coulomb kernel when the time dependence is dropped.*

Limiting cases

The qualitative form of the free-space Green's function is fixed by two things only: the class of the operator and the spatial dimension. Read across the table to see how the same "invert $L$ " instruction produces radically different kernels.

Case	Result	Meaning
Poisson, 3D	$G=\dfrac{1}{4\pi\lvert\mathbf{x}-\mathbf{x}'\rvert}$	Solves $-\nabla^2G=\delta$ ; the Coulomb/Newton potential, decaying as $1/r$ .
Poisson, 2D	$G=-\dfrac{1}{2\pi}\ln\lvert\mathbf{x}-\mathbf{x}'\rvert$	Logarithmic and not integrable at infinity — 2D electrostatics needs care with the reference point.
Helmholtz, 3D (outgoing)	$G=\dfrac{e^{\,ik\lvert\mathbf{x}-\mathbf{x}'\rvert}}{4\pi\lvert\mathbf{x}-\mathbf{x}'\rvert}$	Sommerfeld radiation condition selects $+ik$ ; the outgoing scattered wave of steady scattering.
Diffusion, $d$ -D	$G=\dfrac{\Theta(t-t')}{[4\pi D(t-t')]^{d/2}}\,e^{-\lvert\mathbf{x}-\mathbf{x}'\rvert^2/4D(t-t')}$	Gaussian heat kernel; smooths instantly, formally infinite signal speed, $\Theta$ enforces time-causality.
Wave, 1+3D (retarded)	$G_\mathrm{ret}=\dfrac{\delta\!\left(t-t'-\lvert\mathbf{x}-\mathbf{x}'\rvert/c\right)}{4\pi\lvert\mathbf{x}-\mathbf{x}'\rvert}$	Sharp light-cone support; Huygens' principle holds in odd space dimension $\ge3$ .

Common trapQuoting "the" Green's function without naming the boundary conditions

A Green's function is not unique until you fix the boundary and initial conditions and, for singular operators, the regularization — and skipping that step is the single most common error. The same operator $-\nabla^2$ has different Green's functions for Dirichlet ( $G|_{\partial\Omega}=0$ ) versus Neumann ( $\partial_nG=-1/A$ ) data; the wave operator has retarded, advanced, and Feynman propagators distinguished only by the $i\epsilon$ prescription — which way the contour steps around the poles. Saying " $G$ for the wave equation is $\delta(t-r/c)/(4\pi r)$ " is meaningful only once you have said "retarded." Two traps ride along: the Neumann problem for $-\nabla^2$ has no ordinary Green's function, because the constant is a zero mode ( $\lambda_0=0$ ) that makes the spectral sum diverge, so it demands the solvability condition $\int_\Omega f\,\mathrm{d}V=-\oint_{\partial\Omega}\partial_nu\,\mathrm{d}S$ ; and in 2D the logarithmic $G$ grows at infinity, so " $G\to0$ as $r\to\infty$ " — automatic in 3D — is simply false.

Worked examples

1 · Force on a charge above a grounded plane (fully worked)

A point charge $q=1\ \mathrm{nC}$ sits at height $d=1\ \mathrm{cm}$ above an infinite grounded conducting plane. Use the Dirichlet Green's function to find the force on it. Take $\varepsilon_0=8.854\times10^{-12}\ \mathrm{F/m}$ .

Step 1. Build the half-space Dirichlet kernel by images. Reflect the source through the plane $z=0$ : for a source at $\mathbf{x}'=(x',y',z')$ the image sits at $\mathbf{x}'^{*}=(x',y',-z')$ , and $G=\frac{1}{4\pi|\mathbf{x}-\mathbf{x}'|}-\frac{1}{4\pi|\mathbf{x}-\mathbf{x}'^{*}|}$ . On $z=0$ the two distances are equal, so $G=0$ — the Dirichlet condition holds — while inside the domain $-\nabla^2G=\delta$ because the image lies outside it. The physical charge behaves as if a mirror charge $-q$ sat at $-d$ .

Step 2. Read off the force. The grounded plane exerts on $q$ exactly the Coulomb attraction of its image $-q$ a distance $2d$ away, so $F=-\dfrac{q^2}{4\pi\varepsilon_0(2d)^2}=-\dfrac{q^2}{16\pi\varepsilon_0 d^2}$ , directed toward the plane. The $1/(2d)^2$ makes the image force four times weaker than two real charges at separation $d$ would suggest.

Step 3. Put in numbers. With $q^2=(10^{-9})^2=10^{-18}\ \mathrm{C^2}$ , $16\pi\varepsilon_0=4.451\times10^{-10}$ , and $d^2=10^{-4}\ \mathrm{m^2}$ , $F=\dfrac{10^{-18}}{4.451\times10^{-10}\times10^{-4}}=2.25\times10^{-5}\ \mathrm{N}$ . The charge is pulled toward the plane with $F\approx 2.2\times10^{-5}\ \mathrm{N}=22\ \mu\mathrm{N}$ , and the induced surface density integrates to exactly $-q$ , so total induced charge is conserved. No infinite series, no relaxation grid — the boundary was solved the moment the image was placed.

2 · How fast heat spreads (the heat kernel)

A point of heat is deposited in a copper bar at $t=0$ ; copper has thermal diffusivity $D=1.1\times10^{-4}\ \mathrm{m^2/s}$ . How long until the disturbance has an rms spread of $1\ \mathrm{cm}$ along the bar, and how does the peak temperature decay?

Step 1. Use the heat kernel's variance. The 1D kernel $G\propto e^{-x^2/4Dt}$ is a Gaussian of variance $\langle x^2\rangle=2Dt$ , so the rms spread is $\sqrt{2Dt}$ — the hallmark $\sqrt{t}$ diffusion law, not the $\propto t$ of ballistic motion. Distance grows only as the square root of time.

Step 2. Solve for the time. Setting $\sqrt{2Dt}=L=10^{-2}\ \mathrm{m}$ gives $t=\dfrac{L^2}{2D}=\dfrac{(10^{-2})^2}{2\times1.1\times10^{-4}}=\dfrac{10^{-4}}{2.2\times10^{-4}}=0.45\ \mathrm{s}$ . Roughly half a second to smear across a centimetre — and to reach $10\ \mathrm{cm}$ would take $100\times$ longer, about $45\ \mathrm{s}$ , the reason thick objects equilibrate so slowly.

Step 3. Track the peak. At $\mathbf{x}=0$ the $d$ -dimensional kernel is $(4\pi Dt)^{-d/2}$ , so in 3D the peak falls as $t^{-3/2}$ ; at $t=1\ \mathrm{s}$ it equals $(4\pi\times1.1\times10^{-4})^{-3/2}\approx1.9\times10^{4}\ \mathrm{m^{-3}}$ . Total heat is conserved ( $\int G\,\mathrm{d}^dx=1$ ) even as the peak collapses — the kernel spreads what it cannot destroy. (For the hyperbolic partner, the retarded kernel $\delta(t-r/c)/(4\pi r)$ would place a pulse from $3\ \mathrm{km}$ away at exactly $t=r/c=10\ \mu\mathrm{s}$ , with the same $1/r$ geometric weight — geometry sets the size, causality the timing.)

Physics of the portrait

The Green's function carries the name of a man who taught himself the mathematics after a day's work in a windmill, and whose central idea sat unread for two decades before it reshaped physics. It is a story of how a single, quietly revolutionary essay became the language of fields.

Origins · 1 minGeorge Green and how this was discovered

George Green

1793–1841 · self-taught miller who founded mathematical potential theory

George Green was a Nottingham miller's son who left school at nine, ran the family windmill, and taught himself mathematics from borrowed books with almost no formal instruction. In 1828 he privately printed An Essay on the Application of Mathematical Analysis to the Theories of Electricity and Magnetism, a work of extraordinary originality: it introduced the word potential for the scalar field whose gradient is the force, proved the identities now called Green's theorem, and — most consequentially — constructed the function that inverts a boundary-value problem by superposing the responses to point sources. Only about fifty subscribers bought it, and it was all but invisible. Green entered Gonville and Caius College, Cambridge, in 1833 at the age of forty, but died in 1841 with his masterwork still unknown. In 1845 the young William Thomson (later Lord Kelvin) stumbled on a copy, grasped its importance, and had it republished; through Thomson and Maxwell the kernel became the standard tool of electrostatics, and in the twentieth century the same object — now called a propagator — became the backbone of quantum field theory. The miller's function is today one of the most-used constructions in all of mathematical physics.

History

1828

— Green's Essay introduces the potential function and the Green's function

1833

— Green matriculates at Gonville and Caius College, Cambridge, aged forty

1845

— Thomson rediscovers Green's Essay and spreads it across Europe

1882

— Kirchhoff casts Huygens' principle as a wave Green's function

1948

— Feynman's propagator recasts the Green's function for quantum fields

Going further

Texts. Stone & Goldbart, Mathematics for Physics, ch. 5–6 (Green functions, PDEs) · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., ch. 10 · Hassani, Mathematical Physics, ch. 21–22 (Green's functions in one and several dimensions) · Byron & Fuller, Mathematics of Classical and Quantum Physics, ch. 7 · Morse & Feshbach, Methods of Theoretical Physics, vol. I, ch. 7 · Jackson, Classical Electrodynamics, 3rd ed., ch. 1–3 and §6.4 (electrostatic and retarded kernels) · Nakahara, Geometry, Topology and Physics, §7 (Green's functions and the Atiyah–Singer setting). Key numerical checks: (1) $-\nabla^2[1/(4\pi r)]=\delta^3$ via the Dirichlet integral $\int_0^\infty\sin(kr)/k\,\mathrm{d}k=\pi/2$ , giving $G(1\,\mathrm{m})=0.0796\,\mathrm{m^{-1}}$ ; (2) the heat-kernel normalization $\int(4\pi Dt)^{-d/2}e^{-x^2/4Dt}\,\mathrm{d}^dx=1$ and $\langle x^2\rangle=2Dt$ , giving $t=L^2/2D=0.45\,\mathrm{s}$ for copper at $L=1\,\mathrm{cm}$ ; (3) the retarded propagator $\delta(t-r/c)/(4\pi r)$ and its $10\,\mu\mathrm{s}$ delay at $3\,\mathrm{km}$ ; (4) the image-charge force $F=q^2/(16\pi\varepsilon_0d^2)=2.2\times10^{-5}\,\mathrm{N}$ for $q=1\,\mathrm{nC}$ , $d=1\,\mathrm{cm}$ .

Builds on: Dirac delta and distributions · Fourier and Laplace transforms · Sturm–Liouville theory and eigenfunction expansions · Contour integration and residues · Vector calculus and Green's identities. Leads to: Method of Stationary Phase · Method of Steepest Descent and Saddle Points · WKB Approximation · Laplace's Method for Integrals · Watson's Lemma.

Ask the tutor: Want to see how the same $1/(4\pi r)$ kernel becomes a Coulomb potential, a scattering amplitude, and a Feynman propagator — just by changing the boundary condition? Ask, and we'll walk one operator through all three.

Is the Green's function the same thing as a propagator?

Yes — "propagator" is the name physicists give a time-dependent Green's function, especially in wave and quantum problems. The retarded propagator of the wave equation propagates a disturbance from a source event to a later field point along the light cone, and Feynman's propagator in quantum field theory is exactly the Green's function of the relevant wave operator (Klein–Gordon, Dirac) with a specific pole prescription. The word emphasizes the dynamical reading — how an impulse travels — while "Green's function" emphasizes the algebraic one — the kernel that inverts $L$ . They are two names for one object, and which you hear depends only on whether the speaker is thinking about fields moving or operators inverting.

Why does the singularity of G depend on the number of dimensions?

Because the divergence theorem ties the flux of $\nabla G$ through a small sphere to a fixed number — one — regardless of the sphere's radius. In $d$ dimensions that sphere's area scales as $r^{d-1}$ , so $\nabla G$ must scale as $r^{-(d-1)}$ , and integrating gives $G\sim r^{-(d-2)}$ for $d\ne2$ : that is $1/r$ in three dimensions and the marginal $\ln r$ in two, where $\nabla G\sim1/r$ integrates to a logarithm. The dimension of space, not any choice of ours, dictates whether the point-source potential falls off, stays flat, or grows — which is why 2D electrostatics behaves so differently from 3D.

What is the difference between retarded and advanced Green's functions, physically?

Both solve the identical wave equation with the identical point source; they differ only in when the response is allowed to be nonzero. The retarded kernel is supported on the future light cone — an effect appears only after its cause, the ordinary arrow of causality radiation theory is built on. The advanced kernel is supported on the past cone — the "response" precedes the source, a mathematically valid but acausal solution. In the frequency domain the distinction is a single sign in the $i\epsilon$ that shifts the poles off the real axis, deciding which way the contour closes. Nature uses the retarded one for radiation; the advanced one is its formal partner, and combinations like the Feynman propagator organize quantum field theory.

Why does the heat kernel imply infinite signal speed, and is that a problem?

The Gaussian $e^{-r^2/4Dt}$ is strictly positive for every $r$ at any $t>0$ , so a pulse deposited at the origin produces an immediate — if exponentially tiny — response arbitrarily far away. Formally that is infinite propagation speed, an artifact of the diffusion equation being parabolic rather than hyperbolic. It is not a physical paradox because the amplitude at large $r$ is unmeasurably small (suppressed by $e^{-r^2/4Dt}$ ), and the diffusion equation is itself only a coarse-grained limit of an underlying finite-speed process (random molecular walks, or the telegrapher's equation, which restores a finite front). For every practical purpose the heat kernel is exact; the "infinite speed" is the price of the continuum approximation.

When should I use the spectral sum for G instead of a Fourier transform?

Use the eigenfunction sum (equation 5) when the domain is bounded and the eigenfunctions are known — a box, a sphere, an interval with fixed ends — because there the continuous spectrum is replaced by a discrete one and the boundary conditions are satisfied term by term. Use the free-space Fourier transform (equation 6) when the domain is all of space and translation invariance makes plane waves the natural basis. The two are one construction, $L^{-1}$ resolved on the spectrum of $L$ , with a discrete versus continuous spectrum. A practical tell: if your problem has walls the wave "sees," expand in modes; if it stretches to infinity, transform.

How do boundary conditions get "baked into" the Green's function?

By construction you demand that $G$ , as a function of the field point, satisfy a homogeneous version of the boundary condition — $G=0$ on the wall for Dirichlet, $\partial_nG$ fixed for Neumann — in addition to $LG=\delta$ . Because the superposition $u=\int Gf$ is a sum of copies of $G$ , whatever homogeneous condition each copy obeys, the sum obeys too, so $u$ automatically meets the boundary requirement without any further adjustment. Inhomogeneous data (a prescribed nonzero $u$ or flux on the wall) is then carried by the surface integral in the boundary representation (equation 4), where the chosen homogeneous condition on $G$ conveniently kills the term you do not know and keeps the one you do. The kernel and the boundary are solved together, once.