Harmonic Functions: Why They’re Nifty

Harmonic functions arise in countless engineering and physics applications. This post will review the definition, associated properties and typical applications of these mathematical marvels.

The Problem

In one dimension, Laplace’s Equation is the differential equation $u^{\prime\prime}(x) = 0$, where $u(x)$ is a twice-differentiable scalar function defined on an open interval of the real number line, i.e., $u(x): \, \Omega \subset \mathbb{R} \to \mathbb{R}$. The only solutions are constant functions like $u(x) = \pi$ or linear functions like $u(x) = 2x + \pi$.

A note on notation: functions that are twice continuously differentiable are denoted $C^2$. For more detail, see this Wiki.

More generally, if $u(x)$ is $C^2$ and $u(x) :\, \Omega \subset \mathbb{R}^{n} \to \mathbb{R}$, $x=(x_1, x_2, … , x_n)$, Laplace’s equation is the homogeneous linear partial differential equation:

$$\begin{eqnarray}
\Delta u(x) ={\partial^2 u \over \partial x_1^2} + {\partial^2 u \over \partial x_2^2} + … + {\partial^2 u \over \partial x_n^2}= \sum_{i=1}^n {\partial^2 u \over \partial x_i^2} = 0
\end{eqnarray}$$

The left-hand side symbol above, $\Delta u(x)$, is called the Laplacian and is sometimes denoted $\nabla^2 u$ instead, which is just $\nabla \cdot \nabla u$, the divergence of the gradient.

The solutions of $\Delta u(x) = 0$ are called harmonic functions. The name has its root in the study of harmonic motion (sound waves, for example) but the equation has proven applicable to a whole universe of problems beyond that. As will be shown, harmonic functions only become interesting in the context of specific boundary conditions and singularities.

Note: harmonic functions can also be complex-valued, $u(x) :\, \Omega \subset \mathbb{R}^{n} \to \mathbb{C}$; such a function is harmonic exactly when its real and imaginary parts are.

Properties and Examples

It’s helpful to provide a few examples first. Namely, what kind of functions are harmonic? Constant and linear functions certainly fit the bill since their second derivatives are always zero. The rest are, of course, non-linear. For example,

$e^x \sin y + x + 1$
$x^2 – y^2 + x$
$e^{-x}\, (x\sin y \,-\, y\cos y)$

They can be complex-valued as well:

$f(z) = (x + iy)^2 = x^2 -y^2 + 2xyi$

Harmonic functions are linearly superimposable, which means that sums and scalar multiples of harmonic functions are themselves harmonic. So based on the examples above we can conclude that $e^x \sin y + x + 1 + 4(x^2 – y^2 + x)$ is also harmonic. And so on.

Harmonic functions are also invariant under rigid motions, meaning harmonic functions can be translated and rotated and still retain their harmonicity. For examples of this, see this Wiki page.

On their own these two properties, invariance and linearity, are incredibly powerful and yet it’s only the tip of the iceberg…

The Mean Value Property

The mean value property (MVP) is arguably the most important property of harmonic functions, as it describes how they behave within bounded regions. The outcome may seem strange, but over any ball contained in the region, the average value of a harmonic function equals its value at the ball’s center. Stranger still, this is also the average value over the surface of that ball. This is true in any number of dimensions too, from an interval to a disk to an n-dimensional ball.

The mean value property can be expressed mathematically as shown below.

$$\begin{eqnarray}
\def\avint{\mathop{\,\rlap{-}\!\!\int}\nolimits}
u(x_0)
= \avint_{ B(x_0,r)} u \, dy
= \avint_{\partial B(x_0,r)}u \, dS
\end{eqnarray}$$

Some of the symbols here may be unfamiliar. $B(x_0,r)$ is shorthand for an n-dimensional ball of radius $r$ centered at $x_0$. The symbol $\partial B(x_0,r)$ refers to the boundary of that ball, a sphere (the skin of an apple, if you will). Lastly, the dashed integral symbol denotes an average: the integral of the function over the region divided by the volume of that region:

$$\begin{eqnarray}
\avint_{ B(x_0,r)} u \, dy = {1 \over Vol} \int_{ B(x_0,r)} u \, dy
\end{eqnarray}$$

In the case of a surface integral, integrate over the surface and divide by the surface area instead.

$$\begin{eqnarray}
\avint_{ \partial B(x_0,r)} u \, dS = {1 \over Area} \int_{ \partial B(x_0,r)} u \, dS
\end{eqnarray}$$

Hence, the MVP equation asserts that these two averages are both equal to the value of the function at the ball’s center. Pretty cool, eh? The proof of the MVP is widely available [EV] (§2.2.2) [AXL] [SG].

Two consequences of this. First, the maximum principle: on a bounded region, a harmonic function attains its maximum and minimum on the boundary, and if the region is connected, it cannot attain them in the interior unless it is constant. Second, and related, if a harmonic function is bounded on all of $\mathbb{R}^n$ it is necessarily constant (Liouville’s theorem).

In the 1881 second edition of his Treatise on Electricity and Magnetism (first published in 1873), James Clerk Maxwell made the following remark about the Laplacian:

I propose therefore to call $\nabla^2q$ the concentration of q at the point P, because it indicates the excess of the value of q at that point over its mean value in the neighbourhood of the point.

One caution when reading this: Maxwell’s $\nabla^2$ carries the opposite sign from the modern Laplacian, which is why his “concentration” is positive when q exceeds its neighbors. In modern terms, the Laplacian at a point measures how much a function differs from the average of its neighbors. If it’s positive, the value at the point is less than that average; if it’s negative, it’s greater. And when it’s zero, the value is exactly the average of those neighboring points [LRX]. This has an interesting application in the problem of edge detection in digital images—finding sudden changes in color intensity that occur where the content changes, such as a face against a background. For this reason, the Laplacian is often used to detect these types of intensity changes [AIS].

Real Analytic and Complex Analysis

Harmonic functions are real-analytic. Being real-analytic on an open set means that near any point in that set, the function can be expressed as a convergent power series. As it turns out, that power series is a Taylor series. As a consequence, harmonic functions are also infinitely differentiable, a.k.a., smooth or regular.

Note: The reverse is not true: a smooth function isn’t necessarily analytic. See this example.

In two dimensions, harmonic functions have a symbiotic relationship with complex analysis.  By definition, a complex function is analytic (holomorphic) if it is complex differentiable on an open set. This is far more demanding than real differentiability, since the derivative at a point must be the same from every direction in the plane, not just from the right and left. This leads to a number of interesting outcomes.

First, if a function is analytic, then its real and imaginary parts are harmonic. This follows from the Cauchy-Riemann equations. For example, the real and imaginary parts of $z^2$ are the harmonic functions $x^2 – y^2$ and $2xy$.

Conversely, on a simply connected region, every real harmonic function is the real part of some analytic function: it has a harmonic conjugate $v$ such that $u + iv$ is analytic. (A real-valued harmonic function is not itself analytic in the complex sense, since a real-valued analytic function must be constant.)

Dirichlet Problem

When the Laplacian yields something other than zero, the resulting inhomogeneous equation is called Poisson’s equation. Following [EV], we write it as $-\Delta u(x) = f(x)$ and use this sign convention throughout. It should be noted that adding a harmonic function to a solution of Poisson’s equation produces another solution. More on Poisson’s equation shortly.

The problem of solving a PDE with a specific boundary function is called the Dirichlet problem. On a bounded region, the solution of the Dirichlet problem for Laplace’s equation, if it exists, is unique (a consequence of the maximum principle). As a consequence, the solution inside the region is uniquely determined by its values on the boundary.

$$\begin{eqnarray}
\begin{cases}
\phantom{-}
\mathop{{}\bigtriangleup}\nolimits
u = 0 & x \in \Omega
\\\phantom{-}
u(x) = g(x) & x \in \partial \Omega
\end{cases}
\end{eqnarray}$$

For reasonably shaped regions, this is an example of a well-posed problem: a solution exists, it is unique, and it depends continuously on the data, so small changes to the boundary function $g$ produce only small changes in the solution. The solution also continuously approaches the boundary function value as $x$ nears the boundary $\partial \Omega$ (this is the notation for the boundary of $\Omega$).

The canonical example of a Laplace Dirichlet problem is a steady-state temperature distribution in a region with one or more heated edges. The image below is an example of this (via Matlab). There are also interesting applications in Brownian Motion.

There are other types of boundary problems too, such as Neumann and Robin, depending on how the boundary constraints are defined. For steady-state heat problems, Neumann conditions are often used to represent adiabatic (insulated) edges.

Operators

The notion of a linear operator $L$ is useful when working with these equations. A linear operator is a mapping from one vector space to another that preserves addition and scalar multiplication. In the realm of differential equations, these vector spaces are typically infinite dimensional vector spaces of real or complex-valued functions subject to constraints.

For example, the space of real-valued continuous functions on an open set or the space of twice-differentiable complex-valued functions of a real variable. $L^2$ is one that often arises in Quantum Mechanics. In all cases, the vectors in these vector spaces are functions. These vector spaces can also be endowed with a norm, which provides a notion of length, and an inner product, which provides a notion of angles; the canonical example is a Hilbert space, a complete inner product space.

The linear part of “linear operator” means that the operation on a sum of terms is the sum of those operations: $L(\alpha x + \beta y) = \alpha L(x) + \beta L(y)$. For harmonic functions, this is consistent with the property of linear superposition already mentioned.

Familiar calculus expressions like the derivative can be expressed as linear operators, for example, $L=\frac{d}{dx}$ or in terms of the function on which it operates, $L(u) = u^{\prime}$.

This is an example of a linear differential operator, one that maps a differentiable function $u(x)$ to a linear combination of its derivatives. For example, $L(u) = u^{\prime\prime} + 2u^{\prime} + u$. Relevant to this discussion, the Laplacian is a linear differential operator, and with our sign convention we take $L = -\Delta$.

There are many vector spaces on which these differential operators can “operate.” The vector space of twice-differentiable functions, $C^2(\mathbb{R})$, for example. Or the space of infinitely differentiable (or smooth) functions $C^{\infty}(\mathbb{R})$ and so on. The point being that a linear operator must have both a linear operation (e.g., the derivative) and a vector space on which it operates.

There are many other types of linear operators. In a system of linear equations, the matrix can be treated as an operator, as can its inverse. And on the topic of inverses, the antiderivative undoes the first derivative, though only up to a constant: since $\frac{d}{dx}$ sends every constant to zero, the antiderivative is only a one-sided inverse unless a boundary condition such as $u(a) = 0$ fixes the constant.

Operators have null spaces, also called kernels (not to be confused with the integral kernel introduced below). These are the functions that the operator maps to zero: $Lu=0$, i.e., $\Delta u =0$. In particular, harmonic functions are the null space of the Laplacian operator.

In operator form, Poisson’s equation takes on the following elegance.

$ L = -\Delta$ (the operator)
$ L u = f $ (Poisson’s equation)
$ u = L^{-1}f $ (its solution!)

With this in mind, the inverse of this operator is the following integral operator

$$\begin{eqnarray}
u(x) = L^{-1}f = \int_{\Omega} G(x, \xi) \, f(\xi) \, d\xi
\end{eqnarray}$$

The new term $G$ is called the kernel of the integral operator. For a linear differential operator, this is just the Green’s function, discussed below. Thus knowing the kernel one can (theoretically, at least) construct a general solution to the differential equation.

Now you may be staring at this thinking, “what if $f = 0$?” as it is with harmonic functions. As written above, this would mean that $u$ is always the zero function, which would make this a very dull topic. The reason this is wrong is because we omitted the boundary conditions, which (as you probably know) are required when solving a differential equation. The operator therefore must also account for those and in doing so, the integral equation will have to include additional terms to represent the boundary conditions. Note that with these types of operators, the boundary conditions are actually embedded into their definition.

The boundary of $\Omega$ is denoted $\partial\Omega$ and so the revised solution will be of this form:

$$\begin{eqnarray}
u = BC(\partial\Omega) + \int_{\Omega} G(x, \xi) \, f(\xi) \, d\xi
\end{eqnarray}$$

This is demonstrated below in the discussion of Poisson Dirichlet problems.

Finally, note that just as with matrices, this inverse exists only if the null space is trivial: $Lu = 0$ must have only the solution $u = 0$. That may seem to contradict calling harmonic functions the null space of the Laplacian, but recall that the boundary conditions are built into the operator’s domain. With homogeneous Dirichlet conditions ($u = 0$ on the boundary), the maximum principle shows the only harmonic function is the zero function, so the inverse exists. In matrix language, zero must not be an eigenvalue: the $\lambda$ of $L u = \lambda u$ must all be non-zero. (For operators on infinite-dimensional spaces the picture is more subtle, but the intuition carries over.)

The Fundamental Solution

The fundamental solution of Laplace’s equation is a function $\Phi(x| \xi)$, which is harmonic everywhere except at a singularity $\xi$. In operator form this is written $-\Delta \Phi(x| \xi)= \delta(x – \xi)$, where $\delta(x – \xi)$ is a Dirac distribution centered at $\xi$. A classic example is the potential field generated by a unit point charge.

The form of the fundamental solution for Laplace’s equation depends on the dimensionality of the problem. With the normalization $-\Delta \Phi = \delta$:

• If $n=1$ it is piecewise linear: $-\tfrac{1}{2}|x – a|$
• For $n=2$ it’s logarithmic: $-\tfrac{1}{2\pi}\ln |x – a| \,, x \ne a$
• For $n \ge 3$ it’s $\dfrac{1}{n(n-2)\alpha(n)}|x – a|^{2-n} \,, x \ne a$

Here $\alpha(n)$ is the volume of the unit ball in $\mathbb{R}^n$.

A few observations around this. The first is that these represent a family of solutions, not a specific one. This is because one can add harmonic functions to this solution without changing the outcome.

The second observation is that for $n>1$, the solutions above are not defined at the point $a$; rather, they are only harmonic on the so-called “punctured region,” $\Omega \backslash \{a\}$. And none of them, including the $n=1$ case, is differentiable at $a$. The singularity that occurs there is of tremendous importance.

Third, these functions are radially symmetric. This again aligns with intuition. Both heat and electric fields radiate radially outward, as does the pull of gravity.

Free space Poisson equation

The solution of Poisson’s equation on all of $\mathbb{R}^n$, the so-called free space solution, can be expressed in terms of the fundamental solution. Again, in the absence of boundary conditions, non-trivial solutions arise as the result of a singularity. This is captured by the Dirac distribution. Since $\Phi(x|\xi)$ depends only on the difference $x – \xi$, we write it as $\Phi(x – \xi)$. We start with a property of the Dirac distribution called the “sifting property”:

$$\begin{eqnarray}
\int_{\mathbb{R}^n}\delta (x-\xi) \, f(\xi)\,d\xi
=
f(x)
\end{eqnarray}$$

The left-hand side is a convolution, so this can be written more compactly as:

$$\begin{eqnarray}
(\delta * f)(x)
=\int_{\mathbb{R}^n}\delta (x-\xi) \, f(\xi)\,d\xi = f(x)
\end{eqnarray}$$

Convolutions have a useful property with respect to differentiation, shown here in the one-dimensional case but true in general too.

$$\begin{eqnarray}
\frac{d}{dx} (f * g)
=\frac{df}{dx} * g =f * \frac{dg}{dx}
\end{eqnarray}$$

Convolution is also commutative, which allows swapping the order of the convolution

$$\begin{eqnarray}
(f * g)(x)
=\int_{\mathbb{R}^n}f(x-\xi) \, g(\xi)\,d\xi
=\int_{\mathbb{R}^n}f(\xi) \, g(x- \xi)\,d\xi
\end{eqnarray}$$

Putting these together

$-\mathop{{}\bigtriangleup}\nolimits \Phi = \delta
\implies$

$-\mathop{{}\bigtriangleup}\nolimits( \Phi * f ) =( -\mathop{{}\bigtriangleup}\nolimits \Phi ) * f = \delta * f = f$

This means that we can express the solution to the free space Poisson equation as a convolution of the fundamental solution with the driving function, $f$.

$$\begin{eqnarray}
u(x) =(\Phi * f)(x) =\int_{\mathbb{R}^n}\Phi(x – \xi) \, f(\xi)\,d\xi
\end{eqnarray}$$

One caveat: $f$ must be well behaved at infinity for this integral to make sense. The standard assumption [EV] is that $f$ is twice continuously differentiable and has compact support, meaning it is zero outside some finite ball. This is written $f \in C_c^2(\mathbb{R}^n)$. (Sufficiently rapid decay at infinity also works.)

Green’s Functions

Heads up, we do a deep-dive of Green’s functions in this post. Check it out.

The Green’s function $G(x, \xi)$ for Poisson’s equation is a fundamental solution, corrected by a harmonic function so that it satisfies homogeneous boundary conditions. It is sometimes called the influence function. As mentioned above, these represent the kernel of the inverse differential operator, $ u = L^{-1}f $.

The second variable in $G(x, \xi)$, $\xi$, indicates that we define a Green’s function at every point in the region, thus providing a mechanism to superimpose them across that region. We set $G(x, \xi)$ to zero on the boundary so that the only component of $u(x)$ on the boundary is $g(x)$, the boundary value function. Green, in his seminal 1828 essay in which this concept was introduced, suggested that we think of the zero potential on the boundary as the result of grounding the conducting surface. But this technique proved to extend far beyond problems in electrostatics.

Collectively, these ideas are expressed by the following accessory equation:

$$\begin{eqnarray}
\begin{cases}
\phantom{-}
-\mathop{{}\bigtriangleup}\nolimits
G(x, \xi) = \delta(x – \xi) & x \in \Omega
\\\phantom{-}
G(x, \xi) = 0 & x \in \partial \Omega
\end{cases}
\end{eqnarray}$$

Green’s functions are specific to the boundary geometry too, meaning they can be reused when solving the Dirichlet Problem for the same region regardless of the forcing and boundary functions.

One side note. The Green’s function here is symmetric: $G(x, \xi) = G(\xi, x)$. This follows from the fact that the Laplacian with Dirichlet boundary conditions is self-adjoint (see [EV], §2.2.4, for a short proof). The upshot is that we can switch variables when doing certain integrations. This comes in handy when proving some of the results shown in this section.

It’s also worth stressing that Green’s functions exist for many different differential equations, not just Poisson’s equation. It is but one of the many techniques available for cracking these often intractable problems. We now demonstrate how it is used in the context of the Dirichlet Problem.

Poisson Dirichlet Problem

In terms of the Green’s function, the general solution for the Poisson equation $-\Delta u = f$ with Dirichlet conditions $u = g$ on $\partial\Omega$ is shown below. The proof is shown in [EV], [SG] and others. But it should look familiar based on the discussion of inverse differential operators above.

$$\begin{eqnarray}
u(x) =
– \int_{\partial \Omega} g(\xi)\frac {\partial G(x, \xi)}{\partial \nu} dS(\xi)
+ \int_{\Omega} f(\xi) G(x, \xi) d\xi
\end{eqnarray}$$

Here $\nu$ is the outward unit normal to the boundary. The term $-\frac {\partial G(x, \xi)}{\partial \nu}$ is called the boundary influence function, which for a disk of radius $R$ becomes the well-known Poisson Kernel. Note that this is just the derivative of the kernel of the inverse operator in the direction of the unit normal to the boundary.

In particular, when $f(x)=0$ as it is for Laplace’s equation, the second term vanishes, leaving the following closed-form solution on a disk of radius $R$:

$$\begin{eqnarray}
u(r, \theta) = \frac{1}{2\pi}\int_0^{2\pi}
\frac
{R^2 – r^2}
{ r^2 – 2Rr\cos(\theta – \psi) + R^2} g(\psi)
d\psi
\end{eqnarray}$$

As a sanity check, at the center ($r = 0$) the fraction equals 1 and the formula reduces to the average of $g$ over the boundary, exactly as the mean value property requires. You should also notice that the denominator vanishes when $r = R$ and $\theta = \psi$, i.e., at the boundary point $(R, \psi)$. The Poisson Kernel therefore represents the value that would be generated by a delta distribution on the boundary at that point. This means that the function value $u(r, \theta)$ is the weighted sum of all such distributions along the boundary. Compare this to how the solution of Poisson’s equation is a sum of Green’s functions across the region, weighted by the forcing function.

The Poisson Kernel for the unit disk can be expressed in rectangular coordinates using the law of cosines:

$$\begin{eqnarray}
u(x) =
{1 -|x|^2 \over 2 \pi}
\int_{\partial B(0,1)} \frac {g(\xi)}{|x-\xi |^2} dS(\xi)
\end{eqnarray}$$

More generally, for the unit ball in n dimensions, where $\alpha(n)$ is again the volume of the unit ball,

$$\begin{eqnarray}
u(x) =
{1 -|x|^2 \over n \alpha(n)}
\int_{\partial B(0,1)} \frac {g(\xi)}{|x-\xi |^n} dS(\xi)
\end{eqnarray}$$

In general, finding the Green’s functions and thus solving Poisson’s equation is not easy since these integrals generally do not have closed-form solutions. But solutions do exist for a few simple geometries (i.e., the shape of the boundary). The solution on a ball is well-known, as is the solution in the half-space, both found using the method of images (see [EV], for example).

Fortunately in two dimensions there are many more options thanks to the correspondence between harmonic functions and complex analysis. In particular, we have the wonderful technique of conformal mappings and by extension, the technique of Schwarz-Christoffel. For example, if $w(z)$ is a conformal map of a simply connected region onto the unit disk, the Green’s function for that region is obtained by carrying the unit disk’s Green’s function back through the map:

$$\begin{eqnarray}
G(z, \xi) =
{1 \over 2\pi} \ln{
|1 – \overline{w(\xi)}\, w(z) | \over |w(z)-w(\xi)|
}
\end{eqnarray}$$

Exterior Problems

Exterior solutions are also interesting. In two dimensions, the bounded solution to the exterior Dirichlet problem for a disk tends toward a finite limit at infinity that happens to be the mean value of the boundary data. That is, the solution at infinity for the region outside a disk is the same as the value at the center of the interior solution. This can be shown fairly easily by taking the limit toward infinity of the Poisson Kernel in the exterior Dirichlet Problem.

$$\begin{eqnarray}
\lim_{r\to\infty}
u(r, \theta) = \frac{1}{2\pi }\int_0^{2\pi}
g(\psi)
d\psi
\end{eqnarray}$$

In higher dimensions, a limiting value at infinity must be prescribed (conventionally zero) to ensure a unique solution. Boundedness alone is not enough: outside the unit sphere in $\mathbb{R}^n$, $n \ge 3$, both $u = 0$ and $u = 1 – |x|^{2-n}$ are bounded, harmonic, and zero on the sphere. For a proof see [AXL] or this paper by Giles Auchmuty and Qi Han. This result is interesting given the previous result for the two-dimensional problem and arguably physically consistent too. See [SG].

Bibliography and Further Reading

[AIS] Sinha, U. The Sobel and Laplacian Edge Detectors. AI Shack.

[ALF] Ahlfors, L., 1979, Complex Analysis, 3rd ed., McGraw-Hill Book Company.

[AXL] Axler, S., Bourdon, P., Ramey, W., 2001, Harmonic Function Theory, 2nd ed., Springer.

[EV] Evans, L., 2010, Partial Differential Equations, 2nd ed. (Graduate Studies in Mathematics), American Mathematical Society.

[EV2] Evans, L., n.d. Overview article on partial differential equations, UC Berkeley.

[HT] Hunter, J., Nachtergaele, B., 2001. Applied Analysis, World Scientific.

[HT2] Hunter, J., 2014. Notes on Partial Differential Equations, University of California at Davis.

[LRX] Lamoureux, M., 2004. The mathematics of PDEs and the wave equation. University of Calgary.

[LV] Levandosky, J., 2003. Math 220B Lectures, Stanford University.

[NH] Nehari, Z., 1952, Conformal Mapping, McGraw-Hill Book Company.

[OL2] Olver, P., 2018, Complex Analysis and Conformal Mapping, University of Minnesota.

[STR] Strauss, W. 2007. Partial Differential Equations: An Introduction, 2nd ed., Wiley.

[SG] Stakgold, I., 1967, Boundary Value Problems of Mathematical Physics, Macmillan.