Waves on donuts

August 2021

The Laplace eigenequation on a sphere, a torus, or any Riemannian manifold is

Δf+λf=0\Delta f + \lambda f = 0

where Δ=x2+y2\Delta = \partial^2_x + \partial^2_y is the Laplace operator on the manifold. On the torus T=R2/Z2\mathbb{T} = \mathbb{R}^2/\mathbb{Z}^2, this equation is well understood and the goal of this note is simply to plot the typical solutions – they look like the donut above.

The only values λ\lambda for which there are solutions to (1) are multiples of 4π24\pi^2,

λn=4π2n\lambda_n = 4\pi^2 n

where nn is an integer that can itself be written as the sum of two squares, ie n=x2+y2n = x^2 + y^2. The function f(z1,z2)=e2iπ(xz1+yz2)f(z_1, z_2) = e^{2i\pi (xz_1 + yz_2)} is easily seen to be an eigenfunction: indeed,

(z12+z22)f(z1,z2)=(2πi)2x12f(z1,z2)+(2πi)2x22f(z1,z2)=4π2(x12+x22)f(z1,z2)=λnf(z1,z2).\begin{aligned}(\partial^2_{z_1} + \partial^2_{z_2}) f (z_1, z_2) &= (2\pi i )^2 x_1^2 f(z_1, z_2) + (2\pi i )^2 x_2^2 f(z_1, z_2) \\ &= -4\pi^2 (x_1^2 + x_2^2)f(z_1, z_2)\\ &= - \lambda_n f(z_1, z_2). \end{aligned}

Sum of squares

Determining if an integer nn satisfies this property, to be the sum of two squares, is an old problem with a long history - it dates back to Diophante. Fermat gave the first characterization of such prime numbers. The general theorem is as follows.

An integer n1n \geqslant 1 can be written as the sum of two squares if and only if, in its prime decomposition, the primes pp with of the form p=4k+3p = 4k + 3 have an odd power. If it is the case, then the number of ways r(n)r(n) to write nn as a sum of two squares is given by

r(n)=4(d1(n)d3(n))r(n) = 4(d_1(n) - d_3(n))

where dk(n)d_k(n) is the number of divisors of nn, not necessarily primes, which are congruent to kk modulo 44.

For example, 10=2×510 = 2 \times 5; it has exactly four divisors 1,2,5,101,2,5,10, two of which are congruent to 1 modulo 4, two others to 2. Consequently, d1(10)=2d_1(10) = 2 and d3(10)=0d_3(10) = 0 and r(10)=8r(10) = 8. These "sum-of-squares-representations" of 10 are

12+3232+12(1)2+32(3)2+1212+(3)232+(1)2(1)2+(3)2(3)2+(1)2\begin{aligned} & 1^2 + 3^2 &&&&&&& 3^2 + 1^2\\ & (-1)^2 + 3^2 &&&&&&& (-3)^2 + 1^2\\ & 1^2 + (-3)^2 &&&&&&& 3^2 + (-1)^2\\ & (-1)^2 + (-3)^2 &&&&&&& (-3)^2 + (-1)^2\\ \end{aligned}

Here is a little function for generating all these representations (it is by no means optimized in any way):

function SOSrep(n::Int)
    """ Returns the list of Sum-Of-Squares representations.
    Returns an empty list if there are none. """
    square_reps = []
    for i in 1:sqrt(n)-1
        if isinteger(sqrt(n - i^2))
            j = sqrt(n - i^2)
            append!(square_reps,  [[i,j], [i,-j], [j,i], [j,-i]])
        end
    end
    if isinteger(sqrt(n))
        j = sqrt(n)
        append!(square_reps, [[0,j], [0, -j], [j, 0], [-j, 0]])
    end
    return square_reps
end

Random arithmetic waves

Let us go back to the Laplace equation (1). Let nn be an integer that can be written as a sum of squares and let λn\lambda_n be the associated eigenvalue. Then, the subspace En\mathscr{E}_n of solutions of (1) has dimensions r(n)r(n), the number of ways to write nn as a sum of squares discussed above. Let S\mathscr{S} be the set of couples x=(x1,x2)x=(x_1,x_2) with n=x12+x22n= x_1^2 + x_2^2; then, a basis of En\mathscr{E}_n is given by the elementary functions

ex: (z1,z2)Texp{2iπ(x1z1+x2z2)} e_{x} : (z_1, z_2) \in \mathbb{T} \mapsto \exp\{ 2i\pi (x_1 z_1 + x_2 z_2) \}

and the solutions of (1) are all the linear combinations of these r(n)r(n) functions, namely the functions of the form

xSαxex \sum_{x \in \mathscr{S}} \alpha_x e_x

where αx\alpha_x are complex numbers. What we call random arithmetic waves are random elements of En\mathscr{E}_n obtained by choosing the αx\alpha_x to be iid standard complex Gaussian variables[1]. The resulting random function fnf_n is a random Gaussian field; it can be viewed as a typical element of En\mathscr{E}_n, and there is a vast litterature on its behaviour. Usually, it is more natural to study real Gaussian fields, and for this we have to ensure that fnf_n is real by enforcing some minor constraints on the αx\alpha_x: we simply ask that αx=αx\alpha_{-x} = \overline{\alpha_x} (complex conjugate).

We generate a random arithmetic wave with the following function:

function generate_RAW(n::Int)
    Λ = SOSrep(n)
    if isempty(Λ)
        return nothing
    else
        E = [(x,y) -> exp(2im*pi*(λ[1]*x + λ[2]*y)) for λ in Λ] #exponentials
        X = [randn() + 1im*randn() for λ in Λ]
        f(x,y) = sum(real(X.*[e(x,y) for e in E])) / sqrt(length(Λ))
        return f
    end
end

Let us generate some arithmetic wave; most numbers have very few sum-of-squares representations (0,4 or 8), but 85=5×1785 = 5 \times 17 has 16 of them.

f = generate_RAW(85)

Also, 71825 has 68 square representations, and 801125 has 124, which is a lot.

Some code

A torus can be parametrized with two coordinates (θ,φ)[0,2π)2(\theta, \varphi) \in [0,2\pi)^2 by the following formulas:

X(θ, φ) = (2 + cos(θ)) * cos(φ)
Y(θ, φ) = (2 + cos(θ)) * sin(φ)
Z(θ, φ) = sin(θ)

And now, let's use Julia's GLMakie, which is Makie's backend for GPU 2d and 3d plotting, to see how ff looks like. We're going to draw a mesh of [0,2π)2[0, 2\pi)^2 with a mesh-size of ε\varepsilon.

ε = 0.01
t = 0:ε:2π+ε ; u = 0:ε:2π+ε
x = [X(θ, φ) for θ in t, φ in u]
y = [Y(θ, φ) for θ in t, φ in u]
z = [Z(θ, φ) for θ in t, φ in u]

For each mesh point (θ,φ)(\theta, \varphi), we'll color the corresponding point of the torus, (X(θ,φ),Y(θ,φ),Z(θ,φ))(X(\theta, \varphi), Y(\theta, \varphi), Z(\theta, \varphi)), with a color representing the field fn(θ,φ)f_n(\theta, \varphi).

using GLMakie

field = [f(θ, φ) for θ in t, φ in u]
fig, ax, pltobj = surface(x, y, z, color = field, 
        colormap = :vik,
        lightposition = Vec3f0(0, 0, 0.8), 
        ambient = Vec3f0(0.6, 0.6, 0.6),
        backlight = 5f0, 
        show_axis = false) 

cbar = Colorbar(fig, pltobj,
        height = Relative(0.4), width = 10 )

fig[1,2] = cbar #hide
set_theme!(figure_padding = 0)#hide

Instead of plotting each color, we can only plot the sign of fnf_n, black if positive and white if negative. The set of points of the torus where fnf_n is zero is called the nodal line while the sets {fn>0}\{f_n >0\} and {fn<0}\{f_n < 0 \} are called nodal sets.

signs = [x > 0 ? 0 : 1 for x in field]
fig2, ax2, pltobj2 = surface(x, y, z, color = signs, 
        colormap = :ice,
        lightposition = Vec3f0(0, 0, 0.8), 
        ambient = Vec3f0(0.6, 0.6, 0.6),
        backlight = 5f0, 
        show_axis = false)

With a few processing and computational power, one can get much detailed pictures !

Here are some pictures of realizations of various arithmetic waves associated with different values of nn. The left picture is the nodal set ({f>0}\{ f > 0 \} is dark), while the right one is a heatmap: white is close to zero, red is negative and green is positive. Note that the left and right picture do not correspond to the same realization of the wave (I'll update this soon to see the same realization).

You can use these pictures if you want, but if so please mention this page.

n=41n = 41, dimension 8

n=100n=100, dimension 12

n=985n=985, dimension 16

n=1105n=1105, dimension 28

n=71825n = 71825, dimension 64

Notes and references

The paper Nodal length fluctuations for RAW by Krishnapur et al. studies the length of the nodal lines of random arithmetic waves, while this paper by Rudnick and Wigman studies the volume of nodal sets. See the references in these papers for more litterature on the topic. For spherical random waves or plane random waves, Dimitry Belyaev has some nice pictures on his website.

Also, a gallery of Makie plots which turned out to be quite useful.

Extra pictures on request.

[1] Usually there is also a normalization by (2r(n))1/2(2r(n))^{-1/2}, but it's not important.