Magic angles

The function {u(x, y, z) = 2z^2 - x^2 - y^2} is notable for the following combination of properties:

  1. It is a harmonic function: {\Delta u = -2 - 2 + 4 = 0}
  2. It vanishes at the origin {(0, 0, 0)} together with its gradient.
  3. It is positive on the cone {C=\{(x, y, z) : z > \sqrt{(x^2+y^2)/2}\}}

The cone C has the opening angle {\theta_3 = \cos^{-1}(1/\sqrt{3}) \approx 57.4^\circ} which is known as the magic angle in the context of NMR spectroscopy. Let us consider the mathematical side of its magic.

If C is replaced by any other cone, the properties 1-2-3 cannot be satisfied by a harmonic function in a neighborhood of the origin. That is, C is the largest cone such that a harmonic function can have a critical point at its vertex which is also a point of its extremum on the cone. Why is that?

Let {u} be a harmonic function in some neighborhood of {(0,0,0)} and suppose {u(0,0,0)=0}, {\nabla u(0,0,0) = 0}, and {u>0} on some cone {C = \{(x, y, z) \colon z>c \sqrt{x^2+y^2+z^2} \}}. Expand {u} into a sum of polynomials {p_2+p_3+\cdots } where each {p_k} is a harmonic homogeneous polynomial of degree {k}. Let {m} be the smallest integer such that {p_m} is not identically zero. Then {p_m} has the same properties as {u} itself, since it dominates the other terms near the origin. We may as well replace {u} by {p_m}: that is, {u} is a spherical harmonic from now on.

Rotating {u} around the {z}-axis preserves all the properties of interest: harmonic, positive on the cone, zero gradient at the origin. Averaging over all these rotations we get a rotationally symmetric function known as a zonal spherical harmonic. Up to a constant factor, such a function is given by {P_m(\cos \phi)} where {\phi} is a spherical coordinate (angle with the {z}-axis) and {P_m} is the Legendre polynomial of degree {m}.

The positivity condition requires {P_m(t) > 0} for {t>c}. In other words, the bound on {\theta} comes from the greatest zero of the Legendre polynomial. As is true for orthogonal polynomials in general, the zeros are interlaced: that is, the zeros of {P_m=0} appear strictly between any two consecutive zeros of {P_{m+1}}. It follows that the value of the greatest zero grows with {m}. Thus, it is smallest when {m=2}. Since {P_2(t) = (3t^2-1)/2}, the zero of interest is {1/\sqrt{3}}, and we conclude that {c \ge  1/\sqrt{3}}. Hence the magic angle.

A piece of magic cone

The magic angle is easy to visualize: it is formed by an edge of a cube and its space diagonal. So, the magic cone with vertex at (0,0,0) is the one that passes through (1, 1, 1), as shown above.

In other dimensions the zonal harmonics are expressed in terms of Gegenbauer polynomials (which reduce to Chebyshev polynomials in dimensions 2 and 4). The above argument applies to them just as well. The relevant Gegenbauer polynomial of degree {2} is {nx^2-1} up to a constant. Thus, in {\mathbb R^n} the magic angle is {\cos^{-1}(1/\sqrt{n})}, illustrated by the harmonic function {u(x)=nx_n^2 - |x|^2}.

This analysis is reminiscent of the Hopf Lemma which asserts that if a positive harmonic function in a smooth domain has boundary value 0 at some point, the normal derivative cannot vanish at that point. The smoothness requirement is typically expressed as the interior ball condition: one can touch every boundary point by a ball contained in the domain. The consideration of magic angles shows that if the function is also harmonic in a larger domain, the interior ball condition can be replaced by the interior cone condition, provided that the opening angle of the cone is greater than {\cos^{-1}(1/\sqrt{n})}.

Rarefaction colliding with two shocks

The Burgers equation {v_t+vv_x =0} is a great illustration of shock formation and propagation, with {v(x,t)} representing velocity at time {t} at position {x}. Let’s consider it with piecewise constant initial condition

\displaystyle u(x,0) = \begin{cases} 1,\quad & x<0 \\ 0,\quad & 0<x<1 \\ 2,\quad & 1<x<2 \\ 0,\quad & x>2 \end{cases}

The equation, rewritten as {(v_t,v_x)\cdot (1,v)=0}, states that the function {v} stays constant on each line of the form {x= x_0 + v(x_0,0)t}. Since we know {v(x_0,0)}, we can go ahead and plot these lines (characteristics). The vertical axis is {t}, horizontal {x}.


Clearly, this picture isn’t complete. The gap near {1} is due to the jump of velocity from {0} to {2}: the particles in front move faster than those in back. This creates a rarefaction wave: the gap gets filled by characteristics emanating from {(1,0)}. They have different slopes, and so the velocity {v} also varies within the rarefaction: it is {v(x,t)=(x-1)/t}, namely the velocity with which one has to travel from {(1,0)} to {(x,t)}.


The intersecting families of characteristics indicate a shock wave. Its propagation speed is the average of two values of {v} to the left and to the right. Let’s draw the shock waves in red.


Characteristics terminate when they run into shock. Truncating the constant-speed characteristics clears up the picture:


This is already accurate up to the time {t=1}. But after that we encounter a complication: the shock wave to the right separates velocity fields of varying intensity, due to rarefaction to the left of the shock. Its propagation is now described by the ODE

\displaystyle \frac{dx}{dt} = \frac12 \left( \frac{x-1}{t} + 0 \right) = \frac{x-1}{2t}

which can be solved easily: {x(t) = 2\sqrt{t} +1}.

Similarly, the shock on the left catches up with rarefaction at {t=2}, and its position after that is found from the ODE

\displaystyle \frac{dx}{dt} = \frac12 \left( 1 + \frac{x-1}{t} \right) = \frac{x-1+t}{2t}

Namely, {x(t) = t-\sqrt{2t}+1} for {t>2}. Let’s plot:


The space-time trajectories of both shocks became curved. Although initially, the shock on the right moved twice as fast as the shock on the left, the rarefaction wave pulls them toward each other. What is this leading to? Yet another collision.

The final collision occurs at the time {t=6+4\sqrt{2} \approx 11.5} when the two shock waves meet. At this moment, rarefaction ceases to exist. The single shock wave forms, separating two constant velocity fields (velocities {1} and {0}). It propagates to the right at constant speed {1/2}. Here is the complete picture of the process:


I don’t know where the holes in the spacetime came from; they aren’t predicted by the Burgers equation.

Nodal lines

Wikipedia article on nodes offers this 1D illustration: a node is an interior point at which a standing wave does not move.

Standing wave and its nodes
Standing wave and its nodes

(At the endpoints the wave is forced to stay put, so I would not count them as nodes despite being marked on the plot.)

A standing wave in one dimension is described by the equation {f''+\omega^2 f=0}, where {\omega} is its (angular) frequency. The function {u(x,t) = f(x)\cos \omega t} solves the wave equation {u_{tt}=u_{xx}}: the wave vibrates without moving, hence the name. In mathematics, these are the (Dirichlet) eigenfunctions of the Laplacian.

Subject to boundary conditions {f(0)=0 = f(\pi)} (fixed ends), all standing waves on the interval {(0,\pi)} are of the form {\sin nx} for {n=1,2,3,\dots}. Their eigenvalues are exactly the perfect squares, and the nodes are equally spaced on the interval.

Things get more interesting in two dimensions. For simplicity consider the square {Q=(0,\pi)\times (0,\pi)}. Eigenfunctions with zero value on the boundary are of the form {f(x,y) = \sin mx \sin ny} for positive integers {m,n}. The set of eigenvalues has richer structure, it consists of the integers that can be expressed as the sum of two positive squares: 2, 5, 8, 10, 13, 17,…

The zero sets of eigenfunctions in two dimensions are called nodal lines. At a first glance it may appear that we have nothing interesting: the zero set of {\sin mx \sin ny} is a union of {n-1} equally spaced horizontal lines, and {m-1} equally spaced vertical lines:

Boring nodal lines
This is a square, not a tall rectangle

But there is much more, because a sum of two eigenfunctions with the same eigenvalue is also an eigenfunction. To begin with, we can form linear combinations of {\sin mx \sin ny} and {\sin nx \sin my}. Here are two examples from Partial Differential Equations by Walter Strauss:

When {f(x,y) = \sin 12x \sin y+\sin x \sin 12y }, the square is divided by nodal lines into 12 nodal domains:

Frequency 145, twelve nodal domains
Eigenvalue 145, twelve nodal domains

After slight perturbation {f(x,y) = \sin 12x \sin y+0.9\sin x \sin 12y } there is a single nodal line dividing the square into two regions of intricate geometry:

Also frequency 145, but two  nodal domains
Also eigenvalue 145, but two nodal domains

And then there are numbers that can be written as sums of squares in two different ways. The smallest is {50=1^2+7^2 = 5^2+5^2}, with eigenfunctions such as

\displaystyle    f(x,y) = \sin x\sin 7y +2\sin 5x \sin 5y+\sin 7x\sin y

pictured below.

Frequency 50
Frequency 50

This is too good not to replicate: the eigenfunctions naturally extend as doubly periodic functions with anti-period {\pi}.

Periodic extension
Periodic extension

3 calculus 3 examples

The function {f(x,y)=\dfrac{xy}{x^2+y^2}} might be the world’s most popular example demonstrating that the existence of partial derivatives does not imply differentiability.


But in my opinion, it is somewhat extreme and potentially confusing, with discontinuity added to the mix. I prefer

\displaystyle  f(x,y)=\frac{xy}{\sqrt{x^2+y^2}}

pictured below.


This one is continuous. In fact, it is Lipschitz continuous because the first-order partials {f_x} and {f_y} are bounded. The restriction of {f} to the line {y=x} is {f(x,y)=x^2/\sqrt{2x^2} = |x|/\sqrt{2}}, which is a familiar single-variable example of a nondifferentiable function.

To unify the analysis of such examples, let {f(x,y)=xy\,g(x^2+y^2)}. Then

\displaystyle    f_x = y g+ 2x^2yg'

With {g(t)=t^{-1/2}}, where {t=x^2+y^2}, we get

\displaystyle    f_x = O(t^{1/2}) t^{-1/2} + O(t^{3/2})t^{-3/2} = O(1),\quad t\rightarrow 0

By symmetry, {f_y} is bounded as well.

My favorite example from this family is more subtle, with a deceptively smooth graph:

Looks like xy
Looks like xy

The formula is

\displaystyle    f(x,y)=xy\sqrt{-\log(x^2+y^2)}

Since {f} decays almost quadratically near the origin, it is differentiable at {(0,0)}. Indeed, the first order derivatives {f_x} and {f_y} are continuous, as one may observe using {g(t)=\sqrt{-\log t}} above.

And the second-order partials {f_{xx}} and {f_{yy}} are also continuous, if just barely. Indeed,

\displaystyle    f_{xx} = 6xy g'+ 4x^3yg''

Since the growth of {g} is sub-logarithmic, it follows that {g'(t)=o(t^{-1})} and {g''(t)=o(t^{-2})}. Hence,

\displaystyle    f_{xx} = O(t) o(t^{-1}) + O(t^{2}) o(t^{-2}) = o(1),\quad t\rightarrow 0

So, {f_{xx}(x,y)\rightarrow 0 = f_{xx}(0,0)} as {(x,y)\rightarrow (0,0)}. Even though the graph of {f_{xx}} looks quite similar to the first example in this post, this one is continuous. Can’t trust these plots.

Despite its appearance, f_{xx} is continuous
Despite its appearance, f_{xx} is continuous

By symmetry, {f_{yy}} is continuous as well.

But the mixed partial {f_{xy}} does not exist at {(0,0)}, and tends to {+\infty} as {(x,y)\rightarrow (0,0)}. The first claim is obvious once we notice that {f_x(0,y)= y\, g(y^2)} and {g} blows up at {0}. The second one follows from

\displaystyle    f_{xy} = g + 2(x^2+y^2) g' + 4x^2y^2 g''

where {g\rightarrow\infty} while the other two terms tend to zero, as in the estimate for {f_{xx}}. Here is the graph of {f_{xy}}.

Up, up and away
Up, up and away

This example is significant for the theory of partial differential equations, because it shows that a solution of the Poisson equation {f_{xx}+f_{yy} = h } with continuous {h} may fail to be in {C^2} (twice differentiable, with continuous derivatives). The expected gain of two derivatives does not materialize here.

The situation is rectified by upgrading the continuity condition to Hölder continuity. Then {f} indeed gains two derivatives: if {h\in C^\alpha} for some {\alpha\in (0,1)}, then {f\in C^{2,\alpha}}. In particular, the Hölder continuity of {f_{xx} } and {f_{yy} } implies the Hölder continuity of {f_{xy} }.

What difference a difference scheme makes

One of the simplest nonlinear PDE is

{\displaystyle u_t = -uu_x}

inviscid Burgers equation in one dimension. The function {u} represents the velocity of a compressible fluid. The equation can model a shock wave: for discontinuous initial data such as

{\displaystyle u(x,0) = \begin{cases} 1, \quad x\le 0 \\ 0, \quad x>0 \end{cases}}

the solution (understood in appropriate weak sense) has the jump from {1} to {0} at all times, but the position of the jump moves to the right. The equation can also demonstrate the creation of shocks from continuous initial data (nonlinear steepening of waves), but I am not going into that here.

Let’s try to solve this PDE numerically. Given the initial {u(x,0)} as a sequence of numbers {\dots 1 1 1 1 1 0 0 0 0 0 \dots} in the top row of a spreadsheet, I use the difference of adjacent entries as approximation to {u_x}, multiply it by {u}, and use the result to calculate the solution at the next moment of time. The time step is set at {\Delta t=0.1}. On the first try, I use the forward difference scheme for {u_x}, which in spreadsheet terms amounts to formulas such as

F2 = F1 - 0.1*F1*(G1-F1)

The result of running this scheme is nothing short of disastrous:

Forward (downwind) difference
Forward (downwind) difference

Let’s try the backward (upwind) difference (“wind direction” is positive, since we work with positive values of velocity). This means using formulas such as

F2 = F1 - 0.1*F1*(F1-E1)

The numerical solution is less ridiculous, but is still quite wrong: it is constant in time.

Backward (downwind) difference
Backward (downwind) difference

Going back to the PDE, rewrite it as

{\displaystyle u_t = -(u^2/2)_x}

This is the conservation form of the PDE, which is more suitable for interpretation of discontinuous solutions. Indeed, it says exactly that the vector field {V(x,t) = (u^2/2,u)} has zero divergence. The property of having zero divergence can be understood without taking derivatives: it means that the flux integral over every closed curve (e.g., over every rectangle) is zero.

Discretizing the conservation form of the PDE leads to finite differences of squares, such as

F2 = F1 - 0.1*(G1^2-F1^2)/2

(conservative downwind method):

Conservative downwind method
Conservative downwind method

Still not good. But the conservative upwind method, with formulas such as

F2 = F1 - 0.1*(F1^2-E1^2)/2

turns out to work:

Conservative upwind method
Conservative upwind method

Here is the plot of the data:

Numerical shock propagation with the conservative upwind method.
Numerical shock propagation with the conservative upwind method.

The shock front looks a bit diffuse, but it sustains itself: note how the value at “4” converges to {1}. Recalling that every row corresponds to time interval of {0.1}, we can estimate the propagation speed to be about {0.5}, which is what the analysis of PDE would predict.