Proof. In view of equation (2.2), for each vector $\xi \in {\mathbb {S}}$ , consider the linear map $M_{\xi }$ on ${\mathbb {C}}^3$ defined as

$$\begin{align*}M_{\xi} V:=V-i\xi\times V\,. \end{align*}$$

More explicitly, $M_{\xi }$ is the matrix

$$ \begin{align*}M_{\xi}=\left( \begin{array}{ccc} -1& -i \xi_3 & i \xi_2 \\ i \xi_3 & -1 & -i \xi_1 \\ -i \xi_2& i \xi_1 &-1 \\ \end{array} \right)\,. \end{align*} $$

The determinant of this matrix is $\det M_{\xi }=\xi _1^2+\xi _2^2+\xi _3^2-1$ , and in fact it is easy to see that $M_{\xi }$ has rank 2 for any unit vector $\xi $ . Since $M_{\xi } p(\xi )=0$ for all $\xi \in {\mathbb {S}}$ and $p(\xi )$ only vanishes if $\xi =\xi _{\pm }$ , we then obtain that the kernel of $M_{\xi }$ is spanned by the vector $p(\xi )$ whenever $\xi $ is not one of the poles $\xi _{\pm }$ . In a neighborhood of the poles, the kernel of $M_{\xi }$ can be described as the linear span of $\widetilde p(\xi ):= q(\xi _2)\, (\xi _1\xi _2+i\xi _3,\xi _2^2-1, \xi _2\xi _3-i \xi _1)$ .

Since $M_{\xi } f(\xi )=0$ in the sense of distributions by equation (2.2), it stems from the above analysis that one can write

$$\begin{align*}f(\xi)= {\alpha}(\xi)\, p(\xi) \end{align*}$$

for $\xi $ away from the poles, and

$$\begin{align*}f(\xi)= \beta(\xi)\, \widetilde p(\xi) \end{align*}$$

in a neighborhood of the poles; here, ${\alpha }$ and $\beta $ are complex-valued scalars. As $p(\xi )\times \widetilde p(\xi )=0$ for all $\xi \in {\mathbb {S}}$ , we immediately infer that

$$\begin{align*}p\times f=0\,. \end{align*}$$

Also, as the support of a function is a closed set, p is bounded away from zero on the support of $\chi $ , so we have that

$$\begin{align*}\varphi:= \chi \frac{f\cdot p}{|p|^2}\in H^s({\mathbb{S}})\,. \end{align*}$$

As f is Hermitian, this proves the first part of the proposition. The second statement follows immediately from the fact that

$$\begin{align*}M_{\xi}[\varphi(\xi)p(\xi)]= \varphi(\xi)\,M_{\xi} p(\xi)=0\,. \end{align*}$$

Proof. The first assertion is a straightforward consequence of the first part of Proposition 2.1 and of the estimate (2.3). Indeed, since f is a compactly supported distribution, then $f\in H^s({\mathbb {S}},{\mathbb {C}}^3)$ for some s. Take any $\varepsilon '\in (\varepsilon _{f,k},\varepsilon )$ , and let us consider a function $\Theta $ as above such that $\|\Theta f\|_{H^{-k}({\mathbb {S}})} <\varepsilon '$ . Since $\varepsilon '>\varepsilon _{f,k}$ , it is obvious that we can assume that $\Theta =1$ in a small neighborhood of the poles $\xi _{\pm }$ . Applying Proposition 2.1, we infer that $\chi f=\varphi p$ with $\chi :=1-\Theta $ and some Hermitian scalar function $\varphi \in H^s({\mathbb {S}})$ . In view of the fact that the map $f\mapsto U_f$ is linear and of the bound (2.3), we then have

$$\begin{align*}\sup_{R>0}\frac1R\int_{B_R} \frac{|U_f(x)-U_{\varphi p}(x)|^2}{1+|x|^{2k}}\, dx=\sup_{R>0}\frac1R\int_{B_R} \frac{|U_{\Theta f}(x)|^2}{1+|x|^{2k}}\, dx \leqslant C\|\Theta f\|_{H^{-k}({\mathbb{S}},{\mathbb{C}}^3)}<C\varepsilon'\,. \end{align*}$$

As finite linear combinations of spherical harmonics are dense in $H^s({\mathbb {S}})$ , if $s=-k$ we can approximate $\varphi $ in the $H^{-k}( {\mathbb {S}})$ norm by a Hermitian function $\varphi '$ of this form; then

$$ \begin{align*} \sup_{R>0}\frac1R\int_{B_R} &\frac{|U_f(x)-U_{\varphi' p}(x)|^2}{1+|x|^{2k}}\, dx\\ &\qquad\qquad \leqslant \sup_{R>0}\frac1R\int_{B_R} \frac{|U_{f}(x)-U_{\varphi p}(x)|^2}{1+|x|^{2k}}\, dx+ \sup_{R>0}\frac1R\int_{B_R} \frac{|U_{(\varphi'-\varphi)p}(x)|^2}{1+|x|^{2k}}\, dx <C\varepsilon \end{align*} $$

provided that $\|\varphi -\varphi '\|_{H^{-k}({\mathbb {S}})}<\varepsilon -\varepsilon '$ .

Finally, to see that $\varepsilon _{f,0}=0$ if $f\in L^2({\mathbb {S}},{\mathbb {C}}^3)$ , let us take a smooth function $\Theta :{\mathbb {R}}^3\to [0,1]$ supported in the unit ball and such that $\Theta (0)=1$ . Setting

$$\begin{align*}\Theta_n(\xi):= \Theta(n \xi- n\xi_+) + \Theta(n \xi- n\xi_-)\,, \end{align*}$$

we trivially get that $\|\Theta _n f\|_{L^2({\mathbb {S}})}\leqslant \|f\|_{L^2({\mathbb {S}})}$ for all $n\geqslant 2$ and that $\Theta _n f$ tends to zero almost everywhere in ${\mathbb {S}}$ as $n\to \infty $ . The dominated convergence theorem then shows that $\|\Theta _n f\|_{L^2({\mathbb {S}})}\to 0$ as $n\to \infty $ , thus proving the claim.

Proof. Let us fix an open set $V\supset K$ whose closure is contained in the open neighborhood where v is defined, and a large ball $B_R\supset \overline V$ . Since $\mathbb R^3\backslash K$ is connected, it is obvious that we can take V so that $\mathbb R^3\backslash \overline V$ is connected as well. By the approximation theorem with decay for Beltrami fields [Reference Enciso and Peralta-Salas14, Theorem 8.3], there is a Beltrami field w that approximates v as

$$\begin{align*}\|w-v\|_{C^k(V)}<\varepsilon \end{align*}$$

and is bounded as $|w(x)|<C/|x|$ . As the Fourier transform of w is supported on ${\mathbb {S}}$ , Herglotz’s theorem [Reference Hörmander26, Theorem 7.1.28] shows that one can write $w=U_f$ for some vector-valued Hermitian field $f\in L^2({\mathbb {S}},{\mathbb {C}}^3)$ that satisfies the distributional equation (2.2). Proposition 2.4 then shows that there exists some Hermitian scalar function $\varphi \in C^{\infty }({\mathbb {S}})$ such that

$$\begin{align*}\|U_f-U_{\varphi p}\|_{L^2(B_R)}<C\varepsilon\, \end{align*}$$

so that $\|v-U_{\varphi p}\|_{L^2(V)}<C\varepsilon $ . As the difference $v-U_{\varphi p}$ satisfies the Helmholtz equation

$$\begin{align*}{\Delta} (v-U_{\varphi p}) + v-U_{\varphi p}=0 \end{align*}$$

in V, and $K\subset \!\subset V$ , standard elliptic estimates then allow us to promote this bound to

$$\begin{align*}\|v-U_{\varphi p}\|_{C^k(K)}<C\varepsilon\,, \end{align*}$$

as we wished to prove.

Proof. For $l\geqslant 0$ and $-l\leqslant m\leqslant l$ , $a_{lm}^2$ are independent, identically distributed random variables with expected value 1. As the number of these variables with $l\leqslant n$ is

$$\begin{align*}\sum_{l=0}^n\sum_{m=-l}^{l}1=(n+1)^2\,, \end{align*}$$

the strong law of large numbers ensures that the sample average, that is, the random variable

$$\begin{align*}X_n:=\frac1{(n+1)^2}\sum_{l=0}^n\sum_{m=-l}^{l}a_{ l m}^2\,, \end{align*}$$

converges to 1 almost surely as $n\to \infty $ . Now, consider the truncation

$$\begin{align*}\varphi_n(\xi) := \sum_{l=0}^n\sum_{m=-l}^{l} i^l \, a_{ l m}\,Y_{lm}(\xi)\,. \end{align*}$$

As the spherical harmonics $Y_{lm}$ are orthonormal, the $L^2$ norm of $\varphi _n$ is

$$\begin{align*}\|\varphi_n\|_{L^2({\mathbb{S}})}^2=\sum_{l=0}^n\sum_{m=-l}^{l}a_{ l m}^2= (n+1)^2 X_n\,, \end{align*}$$

and $\|\varphi _n\|_{L^2({\mathbb {S}})}^2$ tends to $\|\varphi \|_{L^2({\mathbb {S}})}^2$ (which may be infinite) as $n\to \infty $ . Since $X_n\to 1$ almost surely, we obtain from the above formula that $(n+1)^{-2}\|\varphi _n\|_{L^2({\mathbb {S}})}^2$ tends to 1 almost surely. Therefore, $\varphi $ is not in $L^2({\mathbb {S}})$ with probability 1.

On the other hand, since

$$\begin{align*}\|\varphi\|_{H^{-s}({\mathbb{S}})}^2=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\frac{a_{ l m}^2}{(l+1)^{2s}}\,, \end{align*}$$

it is straightforward to see that the expected value

$$\begin{align*}{\mathbb{E}}\|\varphi\|_{H^{-1-\delta}({\mathbb{S}})}^2=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\frac{{\mathbb{E}} a_{ l m}^2}{(l+1)^{2+2\delta}}= \sum_{l=0}^{\infty}\frac{2l+1}{(l+1)^{2+2\delta}} \end{align*}$$

is finite for all $\delta>0$ . Hence, $\varphi \in H^{-1-\delta }({\mathbb {S}})$ almost surely, so $u:= U_{\varphi p}$ is well-defined with probability 1.

To prove the representation formula for u and its convergence, let us begin by noting that

$$ \begin{align*} U_{i^l Y_{lm} p}(x)&= \int_{\mathbb{S}} i^l p(\xi) \,Y_{lm}(\xi)\, e^{i\xi\cdot x}\, d{\sigma}(\xi)\\ &= p(D)\int_{\mathbb{S}} i^l Y_{lm}(\xi)\, e^{i\xi\cdot x}\, d{\sigma}(\xi)\,. \end{align*} $$

Using either the theory of point pair invariants and zonal spherical functions [Reference Canzani and Sarnak8, Proposition 4] or special function identities [Reference Enciso, Peralta-Salas and Romaniega15, Proposition 2.1], the Fourier transform of $Y_{lm}\, d{\sigma }$ has been shown to be

$$\begin{align*}\int_{\mathbb{S}} i^l Y_{lm}(\xi)\, e^{i\xi\cdot x}\, d{\sigma}(\xi)= (2\pi)^{\frac32} Y_{lm} \left(\frac{x}{|x|}\right) \frac{J_{ l + \frac12}(|x|)}{|x|^{1/2}}\,. \end{align*}$$

This permits to formally write u as equation (3.1). To show that this series converges in $C^k$ on compact sets, for any large n, any $N>n$ and any fixed positive integer k consider the quantity

$$\begin{align*}q_{n,N}(x):= \sum_{|{\alpha}|\leqslant k} \left|\sum_{l=n}^N \sum_{m=-l}^l a_{lm} D^{\alpha} p(D)\left[Y_{lm} \left(\frac{x}{|x|}\right) \frac{J_{ l + \frac12}(|x|)}{|x|^{1/2}}\right]\right|\,, \end{align*}$$

where we are using the standard multi-index notation. Since $p(D)$ is a third-order operator, for all $|x|<R$ we obviously have

$$ \begin{align*} q_{n,N}&(x)\leqslant C_k \sum_{l=n}^N \sum_{m=-l}^l |a_{lm}|\|Y_{lm}\|_{C^{k+3}({\mathbb{S}})}\left\| \frac{J_{ l + \frac12}(r)}{r^{1/2}}\right\|_{C^{k+3}((0,R))}\\ &\leqslant C_k \left(\sum_{l=n}^N \sum_{m=-l}^l \frac{a_{lm}^2}{(l+1)^{2+2\delta}}\right)^{\frac12}\left(\sum_{l=n}^N \sum_{m=-l}^l(l+1)^{2+2\delta} \|Y_{lm}\|_{C^{k+3}({\mathbb{S}})}^2 \left\| \frac{J_{ l + \frac12}(r)}{r^{1/2}}\right\|^2_{C^{k+3}((0,R))}\right)^{\frac12}, \end{align*} $$

where here $r:=|x|$ and we have used the Cauchy–Schwartz inequality to pass to the second line. The Sobolev inequality immediately gives

$$\begin{align*}\|Y_{lm}\|_{C^{k+3}({\mathbb{S}})}\leqslant C\|Y_{lm}\|_{H^{k+5}({\mathbb{S}})}\leqslant C(l+1)^{k+5}\,. \end{align*}$$

To estimate the Bessel function, recall the large-degree asymptotics

$$\begin{align*}J_{\nu}(r)\sim (2\pi\nu)^{-\frac12} \left(\frac{er}{2\nu}\right)^{\nu}\,, \end{align*}$$

which holds as $\nu \to \infty $ for uniformly bounded r. As the derivative of a Bessel function can be written in terms of Bessel functions via the recurrence relation

$$\begin{align*}\frac d{dr}J_{\nu}(r)= -J_{\nu+1}(r) + \frac\nu r J_{\nu}(r)\,, \end{align*}$$

it follows that the $C^{k+3}$ norm of $J_{l+\frac 12}(r)/r^{1/2}$ tends to 0 exponentially as $l\to \infty $ on compact sets:

$$\begin{align*}\left\| \frac{J_{ l + \frac12}(r)}{r^{1/2}}\right\|_{C^{k+3}((0,R))}\leqslant \left(\frac{CR}{l}\right)^{l-k-3}\,. \end{align*}$$

Since we have proven that

$$\begin{align*}\sum_{l=0}^{\infty}\sum_{m=-l}^l \frac{a_{lm}^2}{(l+1)^{2+2\delta}}<\infty \end{align*}$$

almost surely, now one only has to put together the estimates above to see that, almost surely, $q_{n,N}(x)$ tends to 0 as $n\to \infty $ uniformly for all $N>n$ and for all x in a compact subset of ${\mathbb {R}}^3$ . This establishes the convergence of the series and completes the proof of the proposition.

Proof. As $a_{lm}$ are independent standard Gaussian variables, ${\mathbb {E}}(a_{lm}a_{l'm'})=\delta _{ll'}\delta _{mm'}$ , so the covariance matrix is

$$ \begin{align*} \kappa_{jk}&(x,y)= {\mathbb{E}}[u_j(x) u_k(y)]={\mathbb{E}}[u_j(x) \overline{u_k(y)}]\\ &= \sum_{l=0}^{\infty}\sum_{m=-l}^l \sum_{l'=0}^{\infty}\sum_{m=-l'}^{l'} i^{l-l'}{\mathbb{E}}(a_{lm}a_{l'm'})\int_{\mathbb{S}}\int_{\mathbb{S}} e^{ix\cdot \xi-iy\cdot \eta} p_j(\xi)\, \overline{p_k(\eta)}\, Y_{lm}(\xi)\,{Y_{l'm'}(\eta)}\, d{\sigma}(\xi)\, d{\sigma}(\eta)\\ &= \sum_{l=0}^{\infty}\sum_{m=-l}^l \int_{\mathbb{S}}\int_{\mathbb{S}} e^{ix\cdot \xi-iy\cdot \eta} p_j(\xi)\, \overline{p_k(\eta)}\, Y_{lm}(\xi)\,{Y_{lm}(\eta)}\, d{\sigma}(\xi)\, d{\sigma}(\eta)\,. \end{align*} $$

Here, we have used that u and the spherical harmonics $Y_{lm}$ are real-valued. Since $Y_{lm}$ is an orthonormal basis, one has that

$$\begin{align*}\sum_{l=0}^{\infty}\sum_{m=-l}^l \int_{\mathbb{S}}\int_{\mathbb{S}} \psi(\xi)\,\phi(\eta)\, Y_{lm}(\xi)\, Y_{lm}(\eta)\, d{\sigma}(\xi)\, d{\sigma}(\eta)=\int_{\mathbb{S}} \psi(\xi)\,\phi(\xi)\, d{\sigma}(\xi) \end{align*}$$

for any functions $\psi ,\phi \in L^2({\mathbb {S}})$ . Hence, we can get rid of the sums in the above formula and write

(3.3) $$ \begin{align} \kappa_{jk}(x,y)& = \int_{\mathbb{S}} e^{i(x-y)\cdot\xi}\,p_j(\xi)\, \overline{p_k(\xi)} \, d{\sigma}(\xi)\,, \end{align} $$

which yields the formula for the spectral measure of u. Using now that p is Hermitian (i.e., $\overline { p(\xi )}=p(-\xi )$ ) and a well-known representation formula for the Bessel function $J_{1/2}$ , the above integral can be equivalently written as

$$ \begin{align*} \int_{\mathbb{S}} e^{ix\cdot\xi}\, p_j(\xi)\, \overline{p_k(\xi)}\, d{\sigma}(\xi)&= p_j(D)\, {p_k(-D)}\int_{\mathbb{S}} e^{ix\cdot\xi} d{\sigma}(\xi)\\ &= (2\pi)^{\frac32}p_j(D)\, {p_k(-D)}\frac{J_{1/2}(|x|)}{|x|^{1/2}}\,. \end{align*} $$

The proposition then follows.

Proof. The formula for the spectral measure computed in Proposition 3.4 implies that

$$\begin{align*}\kappa_{jk}(x,x)=\int_{\mathbb{S}} p_j(\xi)\, \overline{p_k(\xi)}\, d{\sigma}(\xi)\,. \end{align*}$$

As p is a polynomial, the computation then boils down to evaluating integrals of the form $\int _{\mathbb {S}} \xi ^{\alpha }\, d{\sigma }(\xi )$ , where ${\alpha }=({\alpha }_1,{\alpha }_2,{\alpha }_3)$ is a multi-index and $\xi ^{\alpha }:=\xi _1^{{\alpha }_1}\xi _2^{{\alpha }_2}\xi _3^{{\alpha }_3}$ . These integrals can be computed in closed form [Reference Folland18]:

(3.4) $$ \begin{align} \int_{\mathbb{S}} \xi^{\alpha}\, d{\sigma}(\xi)= \begin{cases} 2\big[\prod_{j=1}^3\Gamma(\frac{{\alpha}_j+1}2)\big] /\Gamma(\frac{|{\alpha}|+3}2) & \text{if } {\alpha}_1,{\alpha}_2,{\alpha}_3 \text{ are even,}\\[1mm] 0 & \text{otherwise.} \end{cases} \end{align} $$

Here, $\Gamma $ denotes the gamma function.

Armed with this formula and taking into account the explicit expression of the polynomial $p(\xi )$ (cf. equation (2.4)), a tedious but straightforward computation shows

$$\begin{align*}\int_{\mathbb{S}} p_j(\xi)\, \overline{p_k(\xi)}\, d{\sigma}(\xi)=\delta_{jk}\,. \end{align*}$$

The result then follows.

Proof. Since the covariance kernel $\kappa (x,y)$ only depends on $x-y$ , the probability measure ${\mu _u}$ is translationally invariant. Also, note that $(y,w)\mapsto \tau _yw$ defines a continuous map

$$\begin{align*}{\mathbb{R}}^3\times C^k({\mathbb{R}}^3, {\mathbb{R}}^3)\to C^k({\mathbb{R}}^3, {\mathbb{R}}^3)\,, \end{align*}$$

so the map $(y,w)\mapsto \Phi (\tau _yw)$ is measurable on the product space ${\mathbb {R}}^3 \times C^k({\mathbb {R}}^3, {\mathbb {R}}^3)$ . Wiener’s ergodic theorem [Reference Nazarov and Sodin30, Reference Becker5] then ensures that, for $\Phi $ as in the statement, there is a random variable $\Phi ^*\in L^1(C^k({\mathbb {R}}^3\times {\mathbb {R}}^3),{\mu _u})$ such that

as $R\to \infty $ . Furthermore, $\Phi ^*$ is translationally invariant (i.e., $\Phi ^*\circ \tau _y=\Phi ^*$ for all $y\in {\mathbb {R}}^3$ almost surely) and ${\mathbb {E}}\Phi ^* = {\mathbb {E}}\Phi $ .

Also, as the spectral measure (computed in Proposition 3.4 above) has no atoms, a theorem of Grenander, Fomin and Maruyama (see, e.g., [Reference Nazarov and Sodin30, Appendix B] or [Reference Grenander21] and note that the proof carries over to the multivariate and vector-valued case) ensures that the action of the translations $\{\tau _y: y\in {\mathbb {R}}^3\}$ on the probability space $(C^k({\mathbb {R}}^3,{\mathbb {R}}^3),{\mathfrak {S}},{\mu _u})$ is ergodic. As the measurable function $\Phi ^*$ is translationally invariant, one then infers that $\Phi ^*$ is constant ${\mu _u}$ -almost surely. As $\Phi $ and $\Phi ^*$ have the same expectation, then $\Phi ^*={\mathbb {E}} \Phi $ almost surely. The proposition then follows.

Proof. By Proposition 2.5, there exists a Hermitian finite linear combination of spherical harmonics,

$$\begin{align*}\varphi=\sum_{l=0}^n\sum_{m=-l}^l i^l {\alpha}_{lm} Y_{lm}\,, \end{align*}$$

where ${\alpha }_{lm}$ are real numbers (not random variables), such that $\|v-U_{\varphi p}\|_{C^k(K)}<\varepsilon /4$ . Hence

$$\begin{align*}{\mu_u}\big (\big\{ w\in C^k({\mathbb{R}}^3,{\mathbb{R}}^3): \|w-v\|_{C^k(K)}<\varepsilon\big\}\big) \geqslant {\mathbb{P}}\bigg (\bigg\{\|u-U_{\varphi p}\|_{C^k(K)}<\frac\varepsilon4\bigg\}\bigg)\,, \end{align*}$$

where ${\mathbb {P}}$ denotes the natural Gaussian probability measure on the space of sequences $(a_{lm})$ .

Proposition 3.2 shows that the series

$$\begin{align*}\sum_{l=0}^{\infty}\sum_{m=-l}^l a_{lm} U_{i^lY_{lm}p} \end{align*}$$

converges in $C^k(K)$ almost surely, so for any fixed $\delta>0$ there exists some number N (which one can assume larger than n) such that

$$\begin{align*}{\mathbb{P}}\bigg(\bigg\{\bigg\| \sum_{l=N+1}^{\infty}\sum_{m=-l}^l a_{lm} U_{i^l Y_{lm}p}\bigg\|_{C^k(K)}<\frac\varepsilon8\bigg\}\bigg)>1-\delta\,. \end{align*}$$

With the convention that ${\alpha }_{lm}:=0$ for $l>n$ , note that

$$\begin{align*}\|u-U_{\varphi p}\|_{C^k(K)}\leqslant \sum_{l=0}^N\sum_{m=-l}^l |a_{lm}-{\alpha}_{lm}| \|U_{i^lY_{lm}p}\|_{C^k(K)} + \left\| \sum_{l=N+1}^{\infty}\sum_{m=-l}^l a_{lm} U_{i^lY_{lm}p}\right\|_{C^k(K)}\,. \end{align*}$$

Therefore, if we set $M:= 8(N+1)^2\max _{l\leqslant N}\max _{-l\leqslant m\leqslant l} \|U_{i^lY_{lm}p}\|_{C^k(K)}$ , it follows that

$$ \begin{align*} {\mathbb{P}} &\left(\left\{\|u-U_{\varphi p}\|_{C^k(K)}<\frac\varepsilon4\right\}\right) \\ &\qquad\qquad\qquad\geqslant {\mathbb{P}}\bigg (\bigg\{\bigg\| \sum_{l=N+1}^{\infty}\sum_{m=-l}^l a_{lm} U_{i^lY_{lm}p}\bigg\|_{C^k(K)}<\frac\varepsilon 8 \bigg\}\bigg)\, \prod_{l=0}^N \prod_{m=-l}^l {\mathbb{P}}\left(\left\{ |a_{lm}-{\alpha}_{lm}|<\frac\varepsilon M\right\}\right)\,, \end{align*} $$

which is strictly positive. The proposition then follows.

Proof. Let us take cylindrical coordinates $(z,r,\theta )\in {\mathbb {R}}\times {\mathbb {R}}^+\times {\mathbb {T}}$ , with ${\mathbb {T}}:={\mathbb {R}}/2\pi \mathbb Z$ , defined as

$$\begin{align*}z:= x_3 \,,\qquad (r\cos\theta,r\sin\theta):=(x_1,x_2)\,. \end{align*}$$

We now consider the axisymmetric vector field v in ${\mathbb {R}}^3$ given by

(5.1) $$ \begin{align} v:=\frac1r\Big(\partial_r\psi\, E_z-\partial_z\psi\,E_r+\frac{\psi}{r}\,E_{\theta}\Big)\,. \end{align} $$

Here,

$$\begin{align*}\psi:=\cos z+3rJ_1(r) \end{align*}$$

with $J_1$ being the Bessel function of the first kind and order 1, and the vector fields

$$\begin{align*}E_z:= (0,0,1)\,,\qquad E_r:= \frac1r(x_1,x_2,0)\,,\qquad E_{\theta}:= (-x_2,x_1,0)\,, \end{align*}$$

which are often denoted by ${\partial }_z$ , ${\partial }_r$ , ${\partial }_{\theta }$ in the dynamical systems literature, have been chosen so that

$$\begin{align*}E_z\cdot\nabla\phi={\partial}_z\phi\,,\qquad E_r\cdot\nabla\phi={\partial}_r\phi\,,\qquad E_{\theta}\cdot\nabla\phi={\partial}_{\theta}\phi \end{align*}$$

for any function $\phi $ . Notice that $v\cdot \nabla \psi =0$ , so the scalar function $\psi $ is a first integral of v. This means that the trajectories of the field v are tangent to the level sets of $\psi $ .

The vector field v is not defined on the z-axis, so we shall consider the domain in Euclidean 3-space

$$ \begin{align*}\Omega:=\left\{(z,r,\theta):\, (z,r)\in{\mathcal{D}}\,,\; \theta\in{\mathbb{T}}\right\}\,, \end{align*} $$

where ${\mathcal {D}}$ is the domain in the $(z,r)$ -plane given by

$$ \begin{align*}{\mathcal{D}}:=\left\{(z,r):\,-10<z<10,\, \frac{9}{10}<r<\frac{18}{5}\right\}\,. \end{align*} $$

The reason for choosing this particular domain of ${\mathbb {R}}^3$ will become clear later in the proof; for the time being, let us just note that $\psi (z,r)>0$ if $(z,r)\in {\mathcal {D}}$ .

Also, observe that, away from the axis $r=0$ , the vector field v is smooth and satisfies the Beltrami field equation $\operatorname {\mathrm {curl}} v=v$ .

We claim that, in $\Omega $ , v has two hyperbolic periodic orbits joined by a heteroclinic cycle. Indeed, noticing that

$$\begin{align*}(\partial_z\psi,\partial_r\psi)=(-\sin z,3rJ_0(r))\,, \end{align*}$$

where we have used the identity ${\partial }_r[rJ_1(r)]=rJ_0(r)$ , it follows that the points $p_{\pm }:=(\pm \pi ,j_{0,1})\in {\mathcal {D}}$ are critical points of $\psi $ . Here, $j_{0,1}=2.4048\dots $ is the first zero of the Bessel function $J_0$ . Plugging this fact in equation (5.1), this implies that, on the circles in 3-space

$$\begin{align*}\gamma_{\pm}:= \{(z,r,\theta): (z,r)=p_{\pm}\,,\; \theta\in{\mathbb{T}}\}\,, \end{align*}$$

the field v takes the form

$$\begin{align*}v(p_{\pm},\theta)=\frac{c_0}{j_{0,1}^2}\,E_{\theta} \end{align*}$$

with $c_0:=3j_{0,1}J_1(j_{0,1})-1>0$ . Therefore, we conclude that the circles $\gamma _{\pm }$ are periodic orbits of v contained in $\Omega $ .

It is standard that the stability of these periodic orbits can be analyzed using the associated normal variational equation. Denoting by $(v_z,v_r,v_{\theta })$ the components of the field v in the basis $\{E_z,E_r,E_{\theta }\}$ , this is the linear ordinary differential equation (ODE)

$$\begin{align*}\dot\eta = A\eta\,, \end{align*}$$

where $\eta $ takes values in ${\mathbb {R}}^2$ and A is the constant matrix

$$\begin{align*}A:= \frac{{\partial}(v_z,v_r)}{{\partial}(z,r)}\bigg|_{(z,r)=p_{\pm}}= \left( \begin{array}{cc} 0 & 3J^{\prime}_0(j_{0,1}) \\ -{1}/{j_{0,1}} &0 \\ \end{array} \right)\,. \end{align*}$$

The Lyapunov exponents of the periodic orbit $\gamma _{\pm }$ are the eigenvalues of the matrix A. Therefore, since $J^{\prime }_0(j_{0,1})<0$ , these periodic orbits have a positive and a negative Lyapunov exponent, so they are hyperbolic periodic orbits of saddle type.

Since $\psi $ is a first integral of v and $\psi (p_{\pm })=c_0$ , the set

$$\begin{align*}\{(z,r,\theta): \psi(z,r)=c_0\} \end{align*}$$

is an invariant singular surface of the vector field v. This set contains two regular surfaces $\Gamma _1$ and $\Gamma _2$ diffeomorphic to a cylinder. We label them so that $\Gamma _1$ is contained in the half space $\{r\leqslant j_{0,1}\}$ and $\Gamma _2$ in $\{r\geqslant j_{0,1}\}$ . The boundaries of these cylinders are the periodic orbits $\gamma _{\pm }$ . The surface $\Gamma _1$ is the stable manifold of $\gamma _+$ that coincides with an unstable manifold of $\gamma _-$ , while $\Gamma _2$ is the unstable manifold of $\gamma _+$ that coincides with a stable manifold of $\gamma _-$ . Hence, the union $\Gamma _1\cup \Gamma _2$ of both cylinders then form an heteroclinic cycle between the periodic orbits $\gamma _+$ and $\gamma _-$ , and one can see that it is contained in $\Omega $ .

Let us now perturb the Beltrami field v in $\Omega $ by adding a vector field w (to be fixed later) that also satisfies the Beltrami field equation $\operatorname {\mathrm {curl}} w=w$ . Our goal is to break the heteroclinic cycle $\Gamma _1\cup \Gamma _2$ in order to produce transverse intersections of the stable and unstable manifolds of $\gamma _+^{\varepsilon }$ and $\gamma _-^{\varepsilon }$ , where $\gamma _{\pm }^{\varepsilon }$ denote the hyperbolic periodic orbits of the perturbed vector field

$$\begin{align*}X:=v+\varepsilon w=\left(\frac{\partial_r\psi}{r}+\varepsilon w_z\right)\, E_z+\left (-\frac{\partial_z\psi}{r}+\varepsilon w_r\right)\,E_r+\left (\frac{\psi}{r^2}+\varepsilon w_{\theta}\right)\,E_{\theta}\,. \end{align*}$$

As before, $(w_z,w_r,w_{\theta })$ denote the components of the vector field w in the basis $\{E_z,E_r,E_{\theta }\}$ , which are functions of all three cylindrical coordinates $(z,r,\theta )$ . If $\varepsilon>0$ is small enough, the $\theta $ -component of X is positive on the domain $\Omega $ , so we can divide X by the factor $X_{\theta }:=\frac {\psi }{r^2}+\varepsilon w_{\theta }>0$ to obtain another vector field Y that has the same integral curves up to a reparametrization:

(5.2) $$ \begin{align} Y:=\frac{X}{X_{\theta}}=\frac{r\partial_r\psi+\varepsilon r^2w_z}{\psi+\varepsilon r^2w_{\theta}}\, E_z+\frac{-r\partial_z\psi+\varepsilon r^2w_r}{\psi+\varepsilon r^2w_{\theta}}\,E_r+E_{\theta}\,. \end{align} $$

Substituting the expression of $\psi (z,r)$ and expanding in the small parameter $\varepsilon $ , the analysis of the integral curves of Y reduces to that of the following nonautonomous system of ODEs in the planar domain ${\mathcal {D}}$ :

(5.3) $$ \begin{align} \frac{d z}{dt}&=\frac{3r^2J_0(r)}{\psi(z,r)}+\varepsilon \left(\frac{r^2w_z(z,r,t)}{\psi(z,r)}-\frac{3r^4J_0(r)w_{\theta}(z,r,t)}{\psi(z,r)^2}\right)+O(\varepsilon^2)\,,\end{align} $$

(5.4) $$ \begin{align}\frac{d r}{dt}&=\frac{r\sin z}{\psi(z,r)}+\varepsilon\left (\frac{r^2w_r(z,r,t)}{\psi(z,r)}-\frac{r^3\sin z \, w_{\theta}(z,r,t)}{\psi(z,r)^2}\right)+O(\varepsilon^2)\,. \end{align} $$

Notice that the dependence on t is $2\pi $ -periodic and that we have replaced $\theta $ by t in the function $w_z(z,r,\theta )$ (and similarly $w_r,w_{\theta }$ ) because the $\theta $ -component of the vector field Y is 1. When $\varepsilon =0$ , one has

(5.5) $$ \begin{align} &\dot z=\frac{3r^2J_0(r)}{\psi(z,r)}\,, \end{align} $$

(5.6) $$ \begin{align} &\hspace{-5pt}\dot r=\frac{r\sin z}{\psi(z,r)}\,. \end{align} $$

Hence, the unperturbed system is Hamiltonian with symplectic form $\omega :=r^{-1}dz\wedge dr$ and Hamiltonian function $H(z,r):=\log \psi (z,r)$ . The periodic orbits $\gamma _{\pm }$ of v and their heteroclinic cycle $\Gamma _1\cup \Gamma _2$ correspond to the (hyperbolic) fixed points $p_{\pm }$ of the unperturbed system joined by two heteroclinic connections $\widetilde \Gamma _k:=\Gamma _k\cap \{\theta =0\}$ , $k=1,2$ . These are precisely the two pieces of the level curve $\{H(z,r)=\log c_0\}$ that are contained in ${\mathcal {D}}$ . Let us denote by

$$ \begin{align*}\gamma_k(t)=(Z_k(t;0,r_k),R_k(t;0,r_k))\end{align*} $$

the integral curves of the separatrices that solve equations (5.5) and (5.6) with initial conditions $(0,r_k)\in \widetilde \Gamma _k$ . Of course, the closure of the set ${\{\gamma _k(t):t\in {\mathbb {R}}\}}$ is $\widetilde \Gamma _k$ , and the stability analysis of the periodic integral curves $\gamma _{\pm }$ readily implies that $\lim _{t\to \pm (-1)^{k+1}\infty }\gamma _k(t)=p_{\pm }$ .

By the implicit function theorem, the perturbed system (5.3) and (5.4) has exactly two hyperbolic fixed points $p_{\pm }^{\varepsilon }\in {\mathcal {D}}$ so that $p_{\pm }^{\varepsilon } \to p_{\pm }$ as $\varepsilon \to 0$ . The technical tool to prove that the unstable (resp. stable) manifold of $p_+^{\varepsilon }$ and the stable (resp. unstable) manifold of $p_-^{\varepsilon }$ intersect transversely when $\varepsilon>0$ is small is the Melnikov function. We define the vector fields $Y_0$ , $Y_1$ , respectively, as the unperturbed system and the first order in $\varepsilon $ perturbation, that is,

$$ \begin{align*} &Y_0:= \frac{3r^2J_0(r)}{\psi(z,r)}\,E_z+\frac{r\sin z}{\psi(z,r)}\,E_r\,,\\ &Y_1:= \left(\frac{r^2w_z}{\psi(z,r)}-\frac{3r^4J_0(r)w_{\theta}}{\psi(z,r)^2}\right)\,E_z+\left(\frac{r^2w_r}{\psi(z,r)}-\frac{r^3\sin z w_{\theta}}{\psi(z,r)^2}\right)\,E_r\,. \end{align*} $$

Since the unperturbed system is Hamiltonian, we can apply Lemma 5.4 below (which is a variation on known results in Melnikov theory) to conclude that if the Melnikov functions

(5.7) $$ \begin{align} M_k(t_0):=\int_{-\infty}^{\infty} \omega(Y_0,Y_1)|_{\gamma_k(t-t_0)}\, dt\, \end{align} $$

have simple zeros for each $k=1,2$ , then the aforementioned transverse intersections exist, and that actually the heteroclinic connections intersect at infinitely many points. The integrand $\omega (Y_0,Y_1)$ denotes the action of the symplectic 2-form $\omega $ on the vector fields $Y_0,Y_1$ , evaluated on the integral curve $\gamma _k(t-t_0)$ . It is standard that the improper integral in the definition of the Melnikov functions is absolutely convergent because of the hyperbolicity of the fixed points joined by the separatrices (see, e.g., [Reference Guckenheimer and Holmes22, Section 4.5]). Also, notice that although [Reference Guckenheimer and Holmes22, Section 4.5] concerns transverse intersections of homoclinic connections, the analysis applies verbatim to transverse intersections of heteroclinic connections.

More explicitly, the Melnikov functions are given by

$$\begin{align*}M_k(t_0)=\frac1{c_0^2}\int_{-\infty}^{\infty} R_k(t)^2\big[w_z(Z_k(t),R_k(t),t)\sin Z_k(t)-3R_k(t)J_0(R_k(t))w_r(Z_k(t),R_k(t),t)\big]\,dt\,, \end{align*}$$

where $R_k(t)\equiv R_k(t;0,r_k)$ and $Z_k(t)\equiv Z_k(t;0,r_k)$ . It is well known that the existence of transverse intersections is independent of the choice of initial condition.

To analyze these Melnikov integrals, let us now choose the particular perturbation

(5.8) $$ \begin{align} w=J_1(r)\sin \theta\, E_z+\frac{J_1(r)}{r}\cos \theta\, E_r-\frac{J_1'(r)\sin\theta}{r}\,E_{\theta}\,. \end{align} $$

It is easy to check that $\operatorname {\mathrm {curl}} w=w$ in ${\mathbb {R}}^3$ ; in fact $w=(\operatorname {\mathrm {curl}}\operatorname {\mathrm {curl}} +\operatorname {\mathrm {curl}})(J_0(r),0,0)$ (or, to put it differently, $w=U_{\varphi ' q(\xi _1)^{-1} p}$ , where the distribution $\varphi '$ on the sphere ${\mathbb {S}}$ is the Lebesgue measure of the equator, normalized to unit mass). With this choice, the Melnikov functions take the form

$$ \begin{align*} c_0^2M_k(t_0)&=\int_{-\infty}^{\infty} R_k(t)^2\big[J_1(R_k(t))\sin Z_k(t)\sin (t+t_0)-3J_0(R_k(t))J_1(R_k(t))\cos (t+t_0)\big]\,dt \\ &=: a_k\sin t_0 +b_k \cos t_0\,, \end{align*} $$

where the constants $a_k,b_k$ are given by the integrals

$$ \begin{align*} a_k&=\int_{-\infty}^{\infty} R_k(t)^2\big[J_1(R_k(t))\sin Z_k(t)\cos t+3J_0(R_k(t))J_1(R_k(t))\sin t\big]\,dt\,,\\ b_k&=\int_{-\infty}^{\infty} R_k(t)^2\big[J_1(R_k(t))\sin Z_k(t)\sin t-3J_0(R_k(t))J_1(R_k(t))\cos t\big]\,dt\,. \end{align*} $$

Since the Hamiltonian function has the symmetry $H(-z,r)=H(z,r)$ , it follows that $R_k(t)=R_k(-t)$ and $Z_k(t)=-Z_k(-t)$ . This immediately yields that $a_1=a_2=0$ . Moreover, it is not hard to compute the constants $b_1$ and $b_2$ numerically:

$$\begin{align*}\quad b_1=3.5508\dots\,,\quad b_2=0.2497\dots \end{align*}$$

Therefore, the function $M_k(t_0)=b_k\cos t_0$ is a nonzero multiple of the cosine, so it obviously has exactly two zeros in the interval $[0,2\pi )$ , which are nondegenerate. It then follows from Lemma 5.4 below that the two heteroclinic connections joining $p_{\pm }^{\varepsilon }$ intersect transversely. In turn, this implies [Reference Wiggins33, Theorem 26.1.3] that each hyperbolic fixed point $p_{\pm }^{\varepsilon }$ has transverse homoclinic intersections, so by the Birkhoff–Smale theorem [Reference Guckenheimer and Holmes22, Theorem 5.3.5] the perturbed system (5.3) and (5.4) (with w given by equation (5.8)) has a compact hyperbolic invariant set on which the dynamics is topologically conjugate to a Bernoulli shift. This set is contained in a neighborhood of the heteroclinic cycle $\widetilde \Gamma _1 \cup \widetilde \Gamma _2$ , and hence in the planar domain ${\mathcal {D}}$ where the system is defined. This immediately implies that the vector field Y defined in equation (5.2), which is the suspension of the nonautonomous planar system (5.3), has a compact normally hyperbolic invariant set K on which its time-T flow is topologically conjugate to a Bernoulli shift, where $T:= 2\pi N$ for some positive integer $N>0$ . The invariant set K is contained in $\Omega $ because it lies in a small neighborhood of the invariant set $\Gamma _1\cup \Gamma _2$ . Since the integral curves of X and Y are the same, up to a reparametrization, K is also a chaotic invariant set of the Beltrami field X in $\Omega $ .

Finally, since ${\mathbb {R}}^3\backslash \overline \Omega $ is connected, and of course the vector field X satisfies the Beltrami equation in an open neighborhood of $\overline \Omega $ , for each $\delta>0$ , Proposition 2.5 shows that there is a Hermitian finite linear combination of spherical harmonics $\varphi $ such that

$$\begin{align*}\|X-U_{\varphi p}\|_{C^1(\Omega)}<\delta\,. \end{align*}$$

If $\delta $ is small enough, the stability of transverse intersections implies that the Beltrami field $U_{\varphi p}$ has a compact chaotic invariant set $K_{\delta }$ in a small neighborhood of K on which a suitable reparametrization of its time-T flow is conjugate to a Bernoulli shift, so the proposition follows.

Proof. Taking $R_0$ so that the horseshoe of $u_0$ is contained in $B_{R_0}$ , the result is a straightforward consequence of the lower semicontinuity of $N^{\mathrm {h}}_u(R)$ , cf. Proposition 5.1.

Proof. If $\varepsilon $ is small enough, the perturbed system has two hyperbolic fixed points $p_{\pm }^{\varepsilon }$ . To analyze how the heteroclinic connection is perturbed, we take a point $p_0\in \widetilde \Gamma $ and we compute the so-called displacement (or distance) function $\Delta (t_0)$ on a section $\Sigma $ based at $p_0$ and transverse to $\widetilde \Gamma $ . Recall that the function $\varepsilon \Delta (t_0)$ gives the distance of the splitting, up to order $O(\varepsilon ^2)$ , between the corresponding stable and unstable manifolds of the perturbed system at the section $\Sigma $ .

A standard analysis, cf. [Reference Holmes25, Equation (22)] or the proof of [Reference Guckenheimer and Holmes22, Theorem 4.5.3], yields the following formula for $\Delta (t_0)$ :

(5.9) $$ \begin{align} \Delta(t_0)=\frac{1}{|Y_0(p_0)|}\,\int_{-\infty}^{\infty} Y_1(\gamma(t-t_0))\times Y_0(\gamma(t-t_0)) e^{-\int_0^{t-t_0}\text{Tr}\, DY_0(\gamma(s))\,ds}\,dt\,, \end{align} $$

where we have omitted the dependence of the integral curve on the initial condition $p_0\in \widetilde \Gamma $ . Here, we are using the notation $X\times Y:= X_1Y_2-X_2Y_1$ for vectors $X,Y\in {\mathbb {R}}^2$ and $\text {Tr}\, DY_0$ is the trace of the Jacobian matrix of the unperturbed field $Y_0$ .

Take coordinates in ${\mathcal {D}}$ , which we will call $(z,r)$ just as in the proof of Proposition 5.2, and write the symplectic form as $\omega =\rho (z,r)\,dz\wedge dr$ , where $\rho (z,r)$ is a smooth function that does not vanish. Let us call here $\{e_z,e_r\}$ the basis of vector fields dual to $\{dz,dr\}$ (which are usually denoted by ${\partial }_z$ and ${\partial }_r$ , as they correspond to the partial derivatives with respect to the coordinates z and r). The Hamiltonian field $Y_0$ reads in these coordinates as

$$ \begin{align*}Y_0=\frac{1}{\rho(z,r)}\Big(\partial_r H\, e_z-\partial_zH\, e_r\Big)\,. \end{align*} $$

Noting that

$$\begin{align*}Y_1(\gamma(t-t_0))\times Y_0(\gamma(t-t_0))=\frac{\omega(Y_0,Y_1)|_{\gamma(t-t_0)}}{\rho(\gamma(t-t_0))} \end{align*}$$

and

(5.10) $$ \begin{align} e^{-\int_0^{t-t_0}\text{Tr}\, DY_0(\gamma(s))\,ds}&=e^{\int_{0}^{t-t_0}Y_0(\gamma(s))\cdot \nabla \log \rho(\gamma(s))\,ds}\end{align} $$

(5.11) $$ \begin{align}&=e^{\int_0^{t-t_0}\frac{d\log \rho(\gamma(s))}{ds}\,ds}=\frac{\rho(\gamma(t-t_0))}{\rho(p_0)}\,, \end{align} $$

equation (5.9) implies that

$$\begin{align*}\Delta(t_0)=\frac{M(t_0)}{|Y_0(p_0)|\rho(p_0)}\,, \end{align*}$$

so the claim follows because $M(t_0)$ coincides with the displacement function up to a constant proportionality factor.