Discretisation Schemes for Level Sets of Planar Gaussian Fields

Smooth random Gaussian functions play an important role in mathematical physics, a main example being the random plane wave model conjectured by Berry to give a universal description of high-energy eigenfunctions of the Laplacian on generic compact manifolds. Our work is motivated by questions about the geometry of such random functions, in particular relating to the structure of their nodal and level sets. We study four discretisation schemes that extract information about level sets of planar Gaussian fields. Each scheme recovers information up to a different level of precision, and each requires a maximum mesh-size in order to be valid with high probability. The first two schemes are generalisations and enhancements of similar schemes that have appeared in the literature (Beffara and Gayet in Publ Math IHES, 2017. https://doi.org/10.1007/s10240-017-0093-0; Mischaikow and Wanner in Ann Appl Probab 17:980–1018, 2007); these give complete topological information about the level sets on either a local or global scale. As an application, we improve the results in Beffara and Gayet (2017) on Russo–Seymour–Welsh estimates for the nodal set of positively-correlated planar Gaussian fields. The third and fourth schemes are, to the best of our knowledge, completely new. The third scheme is specific to the nodal set of the random plane wave, and provides global topological information about the nodal set up to ‘visible ambiguities’. The fourth scheme gives a way to approximate the mean number of excursion domains of planar Gaussian fields.


Introduction
Let Ψ : R 2 → R be a planar Gaussian field, that is, a random function whose finitedimensional distributions are Gaussian random variables. We shall throughout this paper assume that Ψ is stationary and normalised to have zero mean and unit variance at each point. This implies that Ψ may be defined through its positive-definite correlation kernel κ : R 2 → [−1, 1], satisfying κ(0) = 1 and, for each s, t ∈ R 2 , The main objects of study in this paper are the level sets of Ψ , that is, the random sets N := s ∈ R 2 : Ψ (s) = , ∈ R.
Throughout the paper we assume that κ is C 6 at the origin and κ vv = 0 for each unit vector v ∈ S 1 , where for a function f : R 2 → R and a vector v ∈ R 2 we use f v to denote the derivative of f in the direction v. This ensures that, for a fixed level ∈ R, the level set N almost surely consists of a collection of disjoint simple closed curves. We refer to the components of N as the level lines and the components of R 2 \ N as the excursion domains. In the special case = 0, we refer to N := N 0 as the nodal set, and the level lines and the excursion domains as the nodal lines and nodal domains respectively. A discretisation scheme for a level set is a method of extracting information about the level set through discrete observations of the field. To illustrate, suppose the aim is to assess, for a large fixed box, whether there exists an excursion domain that crosses the box horizontally. By choosing a suitable lattice and a sufficiently fine mesh, we would expect to be able to determine this event with high probability by sampling the value of the random field at the vertices of the lattice.
It is not hard to see that, as the size of the box increases, a progressively finer mesh must be used in order to control the errors that arise in this procedure; this is since any event depending on a fixed scale becomes overwhelmingly likely to occur somewhere inside the box. Our main aim is to quantify the optimum scale at which the mesh-size must decrease as the size of the box grows. As illustrated by our results, the optimum scale depends on the exact property of the level set to be extracted.

Level sets of planar Gaussian fields.
Understanding the geometry of the level sets of planar Gaussian fields has been of great interest to mathematical physicists over the last 30 years. Local properties of the level sets, such as their total length in large boxes, are generally well-understood and easy to compute explicitly by direct integral methods. Global properties, such as the number of excursion domains, or crossing events for excursion domains, are much more subtle to understand. For a general review of properties of level sets of Gaussian random fields, see [13].
In regards to the global properties, a breakthrough result of Nazarov and Sodin [12] has established that, under some weak conditions on the kernel κ, the number of nodal domains of Ψ satisfies a law of large numbers (although not considered in [12], an analysis of the proof shows that this holds also for the excursion domains at any level.) To make this precise, for each bounded domain D let N (D) denote the number of nodal domains in D, i.e., the number of components of D \ N . Then, under certain conditions on κ, there exists a constant c N S = c N S (κ) ≥ 0 such that, for any smooth bounded domain D, N (s D)/Area(s D) → c N S almost surely and in mean, where we use s D to denote the scaled set {sx : x ∈ D} ⊆ R 2 . On the other hand, even in the case that c N S is known to exist, there is no known way to compute its value. Indeed, the value of c N S is not known explicitly for any (non-degenerate) Gaussian field. The connectivity properties of the level sets of planar Gaussian fields are also challenging to understand. One important question is whether the level sets are almost surely bounded. This was confirmed in the special case of the nodal set of positively-correlated planar Gaussian fields-i.e., the case that = 0 and κ ≥ 0-in a result of Alexander [2], but the general case remains open. A recent result of Beffara and Gayet [4] has provided, for the first time, more explicit control over the connectivity of the nodal set of positively-correlated planar Gaussian fields, in the form of a Russo-Seymour-Welsh (RSW) estimate. Conditions under which the nodal set N (and hence also the nodal domains R 2 \ N ) satisfies the RSW estimate in (1) were given in [4]: Theorem 4.9 of [4]. Fix δ > 0. Let Ψ be a stationary planar Gaussian field that is almost surely C 4 and such that the distribution of ∇Ψ (0) is non-degenerate. Suppose further that the correlation kernel κ is invariant under reflection through the horizontal axis and under rotation by π/2, that κ(x) ≥ 0, and that κ(x) = o(|x| −α−δ ) as |x| → ∞ for α = 144 + 128 log 4/3 (3/2) ≈ 325. Then the nodal set N satisfies the RSW estimate in Definition 1.
As an application of our results, we significantly weaken the necessary decay exponent under which the RSW estimate is known to hold, from α ≈ 325 to α = 16 (see Theorem 4 below).
A particularly important planar Gaussian field is the random plane wave, which is the stationary, normalised Gaussian field with the correlation kernel κ(s) := J 0 (k|s|), (2) where J 0 is the zeroth Bessel function, k > 0 is a scale parameter which encodes the frequency or inverse wave-length of the plane wave, and | · | denotes the standard 2 distance. The random plane wave is almost surely smooth, and is the canonical stationary, isotropic Gaussian element in the Hilbert space of planar functions satisfying the Helmholtz equation The random plane wave is of particular importance because it is conjectured by Berry [6] to be a universal model for the high-energy eigenfunctions of the Laplacian in domains with chaotic dynamics. It is also known to be a universal scaling limit of many ensembles of eigenfunctions of the Laplacian on Riemannian manifolds. Solutions of Eq. (3) have been extensively studied, and there are many deterministic results that can be applied to the level sets of the random plane wave. For instance, it is known that there are universal positive constants c 1 and c 2 such that the nodal set N intersects every disc of radius c 1 /k, and each nodal domain contains a disc of radius c 2 /k.
Bogomolny and Schmit conjectured in [7] that the nodal set of the random plane wave is well-approximated by a small perturbation of the square lattice with mesh-size 2π/k, and in particular the nodal domains can be modelled by critical percolation clusters on the square lattice. Based on this idea they made a precise conjecture about the constant c N S in (1) for the random plane wave. Recent computer simulations [5,10] have given very strong evidence that this prediction is slightly inaccurate. On the other hand, these simulations also provided evidence that some global properties of the nodal domains and critical percolation clusters match very well. In particular, this is true for the crossing probabilities discussed above.

Discretisation schemes for level sets of random fields.
Discretisation schemes provide an important tool with which to study the level sets of random fields. To see why, consider that many random fields of interest are extremely rigid, for instance the random plane wave is real analytic. This means that care needs to be taken when working with them, in particular, one cannot condition on the event that the field takes certain values in any open domain without determining the whole field. Conditioning instead on values in a discrete lattice gives a possible way to circumvent this problem.
A second importance of discretisation schemes comes from numerical methods, which are extensively used to study the geometric properties of random fields. These are intrinsically discrete, and hence it is important to have good control on the errors that may arise. As an example of what can go wrong, observe that for any large box there is a positive probability that the random plane wave takes positive values at all lattice points in the box, which suggests a positive probability of having one giant nodal domain; as we just saw, this is deterministically prohibited.
Various discretisation schemes for level sets of planar random fields have previously appeared in the literature. An early work was [11], which developed such a scheme for general planar random fields. This is very similar to our first scheme-assessing the validity of the discretisation on the local scale, see Theorem 2-and is based around controlling the event that the level set N crosses an edge twice (what we call a doublecrossing; see Sect. 2). On the other hand, [11] does not, in general, give a bound on the maximum mesh-size for which the scheme is valid; this is given only in the case of one specific Gaussian field.
More close to our results is the discretisation scheme for planar Gaussian fields in [4,Theorem 1.5]. Again this is very similar to our first scheme in Theorem 2, and is also based around controlling double-crossings. Although [4,Theorem 1.5] does give a general bound on the maximum mesh-size, the provided bound is much more restrictive than the one we give in Theorem 2. In particular, the scheme in [4,Theorem 1.5] requires that the mesh-size ε decays as ε = o(s −8−δ ) for some δ > 0, where s is the scale of the domain in which the level sets are discretised. This can be compared to our scheme in Theorem 2, which is valid if ε = o(s −2−δ ).
In [4], the discretisation scheme was combined with a result of Tassion [15] to establish that the RSW estimate in Definition 1 holds for the nodal set of positively-correlated planar Gaussian fields with sufficiently fast decay of correlations. Since we improve the discretisation scheme in [4, Theorem 1.5], we are able to widen the applicability of the RSW estimate; see Remarks 1 and 2 and Theorem 4.

Main results.
Our main results are a series of four discretisation schemes that extract information about the level sets of planar Gaussian fields. Each scheme recovers information to a different level of precision, and each requires a different maximum mesh-size in order to be valid with high probability. We shall present the schemes in decreasing order of precision. Note that we have chosen to state each of our results as an asymptotic statement. In each case, explicit quantitative bounds on the rate of convergence may be recovered from our proofs (see Sect. 6; note the decay rates of these bounds, up to constants, depend only on the mesh-size ε).
Before we present our results, we state the conditions on the correlation kernel under which our results hold, and give some general definitions.

Conditions on the correlation kernel.
Our results will be valid under certain smoothness and non-degeneracy assumptions on the Gaussian field, expressed through conditions on κ at the origin. These conditions are extremely mild, and will be satisfied in most applications. In particular, it is easy to check that the random plane wave satisfies these conditions. Assumption 1. Suppose that the following hold: 1. (Smoothness) The correlation kernel κ is C 6 at the origin; 2. (Non-degeneracy of first derivatives) For all unit vectors v ∈ S 1 , κ vv (0) < 0.
We briefly discuss the relationship between the conditions on κ in Assumption 1 and the resulting properties of Ψ ; see Sect. 3.2 for precise statements. First, condition (1) ensures that Ψ is three-times differentiable and twice continuously differentiable almost surely. Second, since for each s ∈ R 2 and v 1 , v 2 ∈ S 1 , condition (2) ensures that the distribution of ∇Ψ is non-degenerate. Together these conditions are sufficient to ensure that, for a fixed level ∈ R 2 , the components of N are almost surely simple closed curves.
Alternatively, as is done in [13] for instance, the conditions in Assumption 1 could also be reformulated in terms of the spectral measure ρ = ρ(κ) defined by For instance, (1) could be the replaced by the condition that and (2) could be the replaced with the condition that ρ is not supported on any line.

Level sets and their discretisation.
We begin by introducing the discretisation of the level set N on which our schemes are based. We also make precise the sense in which we shall consider the level set to be well-approximated by its discretisation.
Let L = (V, E, F) be a periodic lattice in R 2 , with vertex set V, edge set E, and face/cell set F; we do not assume all faces in F are equal, but denote by d(L) the largest diameter of any f ∈ F. Let a bounded domain D ⊆ R 2 be called polygonal if its boundary ∂ D consists of a finite number of straight edges. Let a polygonal domain P ⊆ R 2 be called L-compatible if its boundary edges are the union of edges in E. For each ε > 0 and bounded domain D ⊆ R 2 , let P ε (D) be the largest polygonal subdomain P ⊆ D such that P is εL-compatible. This is well-defined by the periodicity of L.
We now introduce the discretisation of the level set N for a fixed ∈ R, based on the lattice L at a mesh-size ε > 0. Define the signed excursion domains and let P ε = (P ε v ) v∈εV ∈ {+1, −1} V be the percolation process on εV induced by these regions (i.e. P ε v = 1 if v ∈ S + , and similarly for S − ); this is well-defined, since almost surely εV ∈ S + ∪ S − . Consider the dual lattice L * , with vertex set F, face set V, and edge set E * . Define the discretised level set N ε ⊆ εL * by prescribing that an edge e * ∈ εE * belongs to N ε if and only if the edge e ∈ εE that is dual to e * has endpoints of opposite sign in P ε ; see Fig. 1.
Finally, we introduce the sense in which we consider the level set N to be wellapproximated by its discretisation N ε (see Fig. 2). For each ε > 0, bounded domain D ⊆ R 2 , and sets M 1 , M 2 ⊆ R 2 that do not intersect the vertices of P ε (D), we say that M 1 and M 2 are ε-homeomorphic in D if there exists a homeomorphism h : P ε (D) → P ε (D) that maps M 1 ∩ P ε (D) onto M 2 ∩ P ε (D), that fixes the vertices V ∩ P ε (D), and such that  Examples of a level set N and its discretisation N ε inside a square domain D for mesh-sizes ε > 0 in which (going left to right) (i) N is ε-homeomorphic to N ε in each face f ∈ D∩εL, (ii) N is ε-homeomorphic to N ε in D but not in the two faces marked with grey circles, and (iii) N is not ε-homeomorphic to N ε in D. Here L is the square lattice Note that the constant 3d(L) is not optimal for our results to hold, but it is convenient as an upper bound.
Our discretisation schemes assess how well the level set N is approximated by its discretisation N ε inside a domain D by applying the notion of an ε-homeomorphism on two scales: the local scale and the global scale. On the local scale, we assess whether N and N ε are ε-homeomorphic in each face f ∈ P ε (D) ∩ εF. On the global scale, we assess whether N and N ε are ε-homeomorphic in D. Whether either of these holds depends, in general, on the fineness of the mesh; see Fig. 3.

1.3.3.
Discretisation on the local scale. The first discretisation scheme gives complete topological information about the level set N within each cell in the lattice; see for instance the left panel of Fig. 3. It requires the finest mesh in order to be valid with high probability.
Theorem 2 (Discretisation on the local scale). Suppose that κ satisfies Assumption 1. Fix a level ∈ R and a bounded domain D ⊆ R 2 . Let ε = ε s > 0 be a sequence such that ε = o(s −2 ) as s → ∞. Then, as s → ∞, Remark 1. The discretisation scheme in Theorem 2 is an enhancement of the discretisation scheme in [4,Theorem 1.5], since it implicitly controls the event that the level set N crosses an edge e ∈ εE twice (i.e. a double-crossing; see Sect. 2), and since it is valid for a significantly coarser mesh. This enhancement may be used to improve the main result in [4], although since we improve this result further using the next discretisation scheme, we postpone this discussion (see Remark 2 and Theorem 4).

Discretisation on the global scale.
The second discretisation scheme gives complete topological information about the level set N on a global, rather than local, scale. This is useful for assessing events involving crossings of level sets and excursion domains on macroscopic scales, where it is only important that N and N ε have the same topology globally, rather than within each cell individually; see for instance the central panel of Fig. 3.

Remark 2.
As a corollary of Theorem 3, we improve the main result of [4] on RSW estimates for the nodal set of positively-correlated planar Gaussian fields. In particular, we significantly weaken the required decay exponent of the correlation kernel, from α ≈ 325 to α = 16 (we also weaken slightly the required non-degeneracy conditions, but this is not as important). So as not to disrupt the exposition of our results, and since the proof is not self-contained (it requires knowledge of [4]), we defer our discussion of the proof of this corollary to Appendix C.
On the other hand, it is very likely that the optimal decay exponent under which the RSW estimate holds in general is actually much lower, perhaps as low as α = 1. This appears challenging to prove, and likely requires new ideas.
Theorem 4 (RSW estimate for the nodal set of positively-correlated planar Gaussian fields). Fix δ > 0. Suppose that κ satisfies Assumption 1, and is also invariant under reflection through the horizontal axis and under rotation by π/2. Suppose further that κ(x) ≥ 0, and κ(x) = o(|x| −16−δ ) as |x| → ∞. Then the nodal set N satisfies the RSW estimate in Definition 1. It follows that the nodal domains R 2 \ N also satisfy this estimate. 1.3.5. Discretisation on the global scale up to visible ambiguities. The third discretisation scheme applies only to the nodal set of the random plane wave, and is based on the observation that, because of the strong rigidity of the plane wave, errors in the discretisation of the nodal set are with high probability contained within local regions which display a certain visible signature in their discretisation (see Remark 4 for a heuristic explanation of this). We call such local regions 'visible ambiguities', and if we are willing to accept that the discretisation is only correct up to these ambiguities, the mesh can be coarsened considerably. This is potentially useful in numerical simulations, since it allows the plane wave to be initially evaluated on a coarse mesh, before a finer mesh is used on a few local regions to resolve the visible ambiguities (see Remark 3).
For clarity of presentation, we only state our result in the case of the regular hexagonal lattice, although we believe that similar results hold, with suitable modification, for any periodic lattice (see Remark 6).
Let us first make precise the concept of a visible ambiguity. Fix the level = 0. For a regular hexagon H , we say that H displays a Type 1 error pattern if its vertices, going clockwise, change sign in the percolation process P ε more than twice. For adjacent regular hexagons H 1 and H 2 with common edge e, we say that H 1 and H 2 display a Type 2 error pattern if: 1. The endpoints of the edge e have the same sign in the percolation process P ε ; and 2. Both H 1 and H 2 have vertices of the opposite sign in the percolation process P ε to the endpoints of e.
Remark that, in the context of the discretised nodal set N ε := N ε 0 , the presence of an error pattern means either that the discretised nodal lines intersect and so cannot be a true representation of the topology of the nodal set (Type 1), or come so close together that there is a high probability that they are not a true representation (Type 2); see Fig. 4. For each ε > 0, we define the set of visible ambiguities π ε ⊆ εF to be the connected components of the union over faces f ∈ εF that display a Type 1 error pattern and adjacent faces f 1 , f 2 ∈ εF that display a Type 2 error pattern (these components will usually consist of at most two adjacent faces, but may contain more). We now make explicit the sense in which we compare the nodal set and its discretisation 'up to visible ambiguities'. For a bounded domain D and a set S such that |S ∩ ∂ D| ∈ 2N, define a setS ⊆ D to be a resolution of S in D if: (i) the components of S that intersect ∂ D are a planar matching of the boundary points S ∩ ∂ D; and (ii) each component ofS that does not intersect ∂ D is a simple closed curve; see Fig. 5. For each componentF ∈ π ε , we define the collection N ε (F) of resolutions of the ambiguity at F to be the set of all possible modification of N ε formed by substituting S = N ε ∩F with a setS ⊆F that is a resolution of S inF. We define the discretised nodal sets with resolved ambiguities M ε to be the collection {N ε (F)}F ∈π ε of all possible such resolutions. For each ε > 0 and bounded domain D ∈ R 2 , we say that the nodal sets N and N ε are ε-homeomorphic in D up to visible ambiguities if there exists a M ∈ M ε such that N and M are ε-homeomorphic in D.
Theorem 5 (Discretisation on the global scale up to visible ambiguities). Suppose that κ satisfies (2) and that L is a regular hexagonal lattice. Fix a smooth bounded domain D ⊆ R 2 . Let ε = ε s > 0 be a sequence such that ε = o(s −1/2 ) as s → ∞. Then, as s → ∞, P N and N ε are ε-homeomorphic in s D up to visible ambiguities → 1.
Remark 3. Theorem 5 demonstrates that the nodal set of the random plane wave can be reconstructed using far fewer observations of (the sign of) the field than are needed in the general case. In particular, if the mesh is allowed to be chosen dynamically, then roughly s 3 observations suffice to reconstruct the nodal set in a domain of radius s (i.e. the s 3+δ vertices of a fixed mesh on the scale ε = s −1/2−δ/2 , then a further s 2 vertices chosen dynamically to resolve, with a mesh on the scale ε = s −1 (c.f. Theorem 3), the roughly s visible ambiguities that arise). This can be contrasted with the roughly s 4 points needed to reconstruct a level set in the general case under the scheme in Theorem 3. It is an interesting question whether a dynamic scheme might be able to lower this bound in the general case.

Remark 4.
We briefly explain why the nodal set of the random plane wave has discretisation errors that are easily identifiable as potential ambiguities. Consider that, for a general smooth function, errors in the discretisation of its nodal set are likely to result from either small nodal domains or low-lying saddle points. For the random plane wave, however, small nodal domains are deterministically ruled out. Moreover, because the wave satisfies the Helmholtz equation (3), any low-lying saddle point has principal axes that intersect at approximately right-angles, unless both eigenvalues of the Hessian ∇ 2 Ψ at the saddle point are very small (which is highly unlikely due to eigenvalue repulsion, see Sect. 3.3). This ensures that a low-lying saddle, and hence a possible discretisation error, is very likely to be easily identifiable by observing the field at nearby lattice vertices.
Observe that the situation is different for other level lines, since there will be a much larger range of angles that will be attained by saddles lying close to this level. In particular, the angle between the principal axes of the saddle point will be close to 0 whenever one of the eigenvalues of the Hessian ∇ 2 Ψ at the saddle point, rather than both of them, is small. This is a much more likely event, and as a result no improvement can be made over the discretisation scheme in Theorem 3.

Remark 5.
The exact property of regular hexagonal lattices that is crucial to our result is stated in Lemma 13. To give a brief description, what is important is that a cone at a near-right-angle must induce an error pattern somewhere in the lattice (i.e. defining the percolation process via membership of the interior/exterior of the cone). Note that this property does not hold for, say, the square lattice, since if the cone was centred at the exact centre of a face, and was not exactly right-angled, then it may not induce an error pattern anywhere in the lattice; see Fig. 6. Fig. 6. A cone at near-right-angle induces an error pattern (in this case Type 2) on a regular hexagonal lattice (left), whereas it may not on a square or triangular lattice (centre), unless we widen the definition of Type 2 errors as described in Remark 6 (right) Remark 6. Although we do not establish it rigorously, it is easy to see from our proofs that a similar result to Theorem 5 should hold for any periodic lattice L, only with a different, perhaps more complicated, definition of the Type 2 error patterns used to identify the set of visible ambiguities in collections of nearby cells (in all cases, Type 1 error patterns are defined identically; this only has content if L is not a triangulation). To illustrate, let us give suitable definitions of Type 2 error patterns for the square lattice and the regular triangular lattice-i.e. the lattice in which each face is an equilateral triangle.
For the square lattice, a set of four horizontally or vertically adjacent faces f 1 , f 2 , f 3 , f 4 ∈ εF display a Type 2 error patten if: 1. The edge that is common to the faces f 2 and f 3 has endpoints that have the opposite sign in the percolation process P ε ; and 2. Both f 1 ∪ f 2 and f 3 ∪ f 4 have vertices of the opposite sign in P ε .
For the regular triangular lattice, the set of Type 2 error patterns can be defined on sets of 12 faces by grouping the triangular faces into pairs of hexagons (in each of the three ways of doing this), and applying the same definition as for the hexagonal lattice.
In general, for each lattice L there exists a visibility parameter μ > 0 such that we may define a certain set of Type 2 error patterns on collections of faces that lie inside balls of radius μ, for which the discretisation up to ambiguities holds by analogy to Theorem 5.

Discretisation for estimating the Nazarov-Sodin constant.
The final discretisation scheme gives a way to estimate the constant c N S = c N S (κ) defined in (1). As discussed above, under some general conditions on κ this constant is known to exist, but it has not yet been computed explicitly for any non-degenerate Gaussian field. Hence, there is much interest in developing ways of numerically approximating this constant.
Our discretisation scheme gives a way to compare the number of nodal domains of planar Gaussian fields with the number of excursion domains induced by the discretised nodal set N ε . To this end, for each ε > 0 and bounded domain D, recall the definition of the εL-compatible set P ε (D), and let N ε (D) denote the number of components of the set P ε (D) \ N ε . Theorem 6 (Discretisation for estimating the constant c N S ). Suppose that κ satisfies Assumption 1, and further assume that c N S = c N S (κ) exists and is finite. Fix a smooth bounded domain D ⊆ R 2 . Let ε = ε s > 0 be a sequence such that ε = o(1) as s → ∞. Then, as s → ∞, Remark 7. For applications, it may be useful to have quantitative bounds on the rate of convergence in (5). The proof of Theorem 6 yields the following bound. For each δ ∈ (0, 1), there exists a constant c = c(κ, L, D, δ) > 0 such that, for each s > 0 and ε ∈ (0, 1),

Overview of Proofs and Outline of Paper
The validity of our discretisation schemes are established by controlling certain 'bad' events involving the level sets, each requiring the mesh-size to be set at a maximum level in order to rule them out with overwhelming probability. In this section we give a brief description of these events, and explain how we control them to prove the main results; this information is summarised in Table 1.
The events that we control are the following (see Fig. 7): is called small if its boundary ∂ S intersects at most one edge in εE; 3. Four-crossings-A face f ∈ εF is said to have a four-crossing if four or more edges e ∈ ∂ f are crossed by N (this only has content if L is not a triangulation); 4. Tubular-crossings-An edge e ∈ εE is said to have a tubular-crossing if, denoting by f 1 , f 2 ∈ εF the cells adjacent to e, there is a component of ( f 1 ∪ f 2 ) \ N that separates the endpoints of e (this implies that e has a double-crossing); and 5. Invisible errors (only defined for the regular hexagonal lattice)-An edge e ∈ εE is said to give rise to an invisible error if it has a tubular-crossing and, denoting by f 1 , f 2 ∈ εF the cells adjacent to e, the faces f 1 and f 2 do not display a Type 2 error pattern. In fact, for technical reasons we actually use a slightly modified version of 'tubular-crossing' in our definition; we explain the necessary modification in Sect. 6.
To emphasise the novelty of our approach, we remark that in the discretisation schemes developed in [4] and [11] only the first of these error events were controlled.
For the local discretisation scheme (Theorem 2), it is sufficient to control the first three events, since it is easy to see that N and N ε being ε-homeomorphic in f is equivalent to the events (1)-(3) not occurring in f or its boundary edges. To control these events inside s D we must use a mesh-size ε = o(s −2 ). This is essentially due to double-crossings; to control small excursion domains and four-crossings the mesh-size For the global discretisation scheme (Theorem 3), we no longer seek to control double-crossings but instead use the following additional insight: a double-crossing that results in a discrepancy in the global topology of the level set must either (i) result in a tubular-crossing, or (ii) be located on the boundary ∂ P ε (s D). Hence, in addition to small excursion domains and four-crossings, it is sufficient to control tubular-crossings, which require the mesh-size ε = o(s −1−δ ), and double-crossings on the boundary, which require the mesh-size ε = o(s −1/2 ).
For the global discretisation scheme up to visible ambiguities (Theorem 5), we use the fact that discrepancies in the global topology of the nodal lines must either be correctable by resolving a visible ambiguity (e.g. four-crossings that lead to Type 1 errors, or tubularcrossings that lead to Type 2 errors), or must be the result of invisible errors (note that small nodal domains are deterministically ruled out). Hence we only need to control double-crossings on the boundary and invisible errors, both of which requires the mesh- size ε = o(s −1/2 ). Again, we note that we actually use a slight modification of the definition of tubular-crossings for this task; see Sect. 6.
For the discretisation scheme for estimating the constant c N S (Theorem 6), we control the exact same events as for the global discretisation scheme (Theorem 3). This is natural, since the number of excursion domains is a function of the global topology. On the other hand, since our approximation is in expectation, we need to control the multiplicity of these events, rather than simply their occurrence.
The rest of the paper is structured as follows. In Sect. 3, we collect preliminary properties of Gaussian vectors and fields; most of these results are standard, but they will serve as a significant input into our proofs. In Sect. 4, we undertake an analysis of events involving crossings of the level set and critical points near the level set, based around several Kac-Rice arguments; these arguments rely on certain matrix computations, and in order not to disrupt the flow of the paper we defer these computations to Appendix B. In Sect. 5, which is completely deterministic, we consider the effect of small perturbations on the nodal set of a generic function; the results in this section are only relevant to our analysis of the random plane wave, but are crucial in exploiting the extra control we have in that case, allowing us to reduce the study of the global topology of the nodal set to the study of invisible errors. Finally, in Sect. 6, we combine the above ingredients to establish the main results.

Gaussian Estimates
In this section we collect some preliminary Gaussian estimates.

Basic properties of Gaussian vectors.
In the first part, we state some basic properties of Gaussian vectors. Proposition 1 is standard, whereas Lemmas 1 and 2 follow easily from Hölder's inequality and the compactness of the set of normalised covariances matrices.
Suppose that Σ Y is invertible, and fix a vector y ∈ R n . Then, conditionally on Y = y, the vector X is Gaussian with Lemma 2 (Bound on the maximum of a Gaussian vector under conditioning). Let n ∈ N and ∈ R. Then there exists a c := c(n, ) > 0 such that, for each centred Gaussian vector (X, Y ) of dimension n + 1 with covariance Σ := (s i j ) 1≤i, j≤n+1 , it holds that, for each δ > 0,

Gaussian fields and their derivatives.
We now state some properties of Gaussian fields and their derivatives. Let Ψ : R 2 → R be a stationary, normalised Gaussian field with correlation kernel κ : R 2 → [−1, 1]. The following lemmas are standard; see [1,3].

Lemma 3 (Smoothness and non-degeneracy). The following hold:
1. Suppose that κ is C 2k at the origin for k ∈ N. Then almost surely Ψ is k-times differentiable and (k − 1)-times continuously differentiable. 2. Suppose that κ is C 2k at the origin and, for unit vectors v 1 , . . . v k ∈ (S 1 ) k , Lemma 5 (Covariances between the field and its derivatives). Let s ∈ R be a point and v ∈ S 1 a unit vector. Suppose that κ is C 2 at the origin and C 1 at s. Then, Lemma 6 (Derivative bounds). Suppose that κ is C 2k at the origin for k ∈ N. Fix r > 0.
Then there exists a c = c(κ, k, r ) > 0 such that The following proposition is less standard; we give a full proof in Appendix A. Note that in many important cases-e.g. isotropic κ or κ ∈ L 1 (R 2 )-the proposition follows immediately from the resulting strict positive-definiteness of κ (see [14,Theorem 3.8] and [16,Theorem 6.11] respectively, where for the latter also note that, under Assumption 1, κ is continuous).
3.3. Small ball estimate for derivatives of the random plane wave. To close this section, we give a small ball estimate for the first and second derivatives of the random plane wave; this will be necessary in completing the proof of Theorem 5.
Let Ψ : R 2 → R be the random plane wave, i.e. the stationary, normalised Gaussian field whose correlation kernel satisfies (2) for frequency parameter k > 0. Let λ 1 and λ 2 denote the eigenvalues of the Hessian matrix ∇ 2 Ψ (0). The following proposition essentially expresses the order-one repulsion between the eigenvalues λ 1 and λ 2 ; the existence of this repulsion is not restricted to the plane wave (it holds for any sufficiently smooth isotropic Gaussian field) but we state it only in this case. We give a full proof in Appendix A.

Kac-Rice Arguments
In this section we develop several Kac-Rice arguments which underpin the validity of our discretisation schemes. Throughout this section let Ψ : R 2 → R be a stationary, normalised Gaussian field with correlation kernel κ : R 2 → [−1, 1] satisfying Assumption 1, and let ∈ R be a fixed level. By Lemma 3, this implies that the components of the level set N are almost surely simple closed curves, and also that the distribution of ∇Ψ (s) is non-degenerate for any s ∈ R 2 .
For k ∈ N and a collection of points s = (s 1 , . . . , s k ) ∈ (R 2 ) k , let θ − (s) and θ + (s) denote, respectively, the size of the smallest and largest interior angle among all possible triangles formed by points in s; if k ≤ 2 or the points in s are not distinct, set θ − (s) := 0 and θ + (s) := π . The variable θ − (s) gives a quantitative measure of the degeneracy of s, and in addition θ + (s) identifies degeneracy that is due to co-linearity.

Kac-Rice formulae.
Kac-Rice formulae are a general way to express expectations involving the zeros of a random process in terms of certain integrals. The exact form of a given Kac-Rice formula will depend on the precise quantity that is desired to be computed. We state two Kac-Rice formulae, one for controlling the intersection of the level set with line segments, and another for controlling the set of critical points; since these are standard, we omit the proof.
where ϕ s ( ) denotes the density at of the k-dimensional Gaussian vector Ψ (s). Then where ϕ s (0) denotes the density at zero of the two-dimensional Gaussian vector ∇Ψ (s).

Level set crossings.
We first use the one-dimensional Kac-Rice formula to give bounds for events involving the intersection of level sets with line segments.
Proof. Let v ∈ S 1 be a unit vector in the direction of the line segment L. By the one-dimensional Kac-Rice formula in Proposition 4 (setting k := 1, μ − := 0, μ + := π, L 1 := L), Applying Lemma 4, Ψ v (εs) is independent of Ψ (εs). Moreover, by stationarity, the distribution of Ψ v (x) is identical for each x. Hence which proves the result.
Moreover, if L 1 and L 2 lie on a common line, the bound can be replaced with cm 2 ε 3 .
Proof. Remark that, by a linear rescaling of the plane, we may reduce to the case that κ vv (0) = −2 for all unit vectors v ∈ S 1 ; this simplification will be useful when we appeal to the computation in Appendix B. We may also choose ε sufficiently small so that, for each s = (s 1 , s 2 ) ∈ εL with s i distinct, the distribution of Ψ (s) is non-degenerate; this is possible since κ vv (0) = 0.
Let v 1 and v 2 be unit vectors in the direction of the line segments L 1 and L 2 respectively. By the one-dimensional Kac-Rice formula in Proposition 4 (setting k := 2, μ − := 0, μ + := π ) and using the fact that the density of a two-dimensional Gaussian vector with covariance matrix Σ is at most (2π) −1 |Σ| −1/2 , Fix for a moment s = (s 1 , s 2 ) ∈ L 2 such that s 1 = s 2 . Applying Proposition 1 and Lemma 5, conditionally on Ψ (εs 1 ) = Ψ (εs 2 ) = the vector has a Gaussian distribution with mean μ ε,s , where Applying Lemma 1, we deduce that there exists a c 0 = c 0 ( ) > 0 such that We first treat the case that L 1 and L 2 do not lie on a common line. Observe that, since A ε,s and B ε,s are positive definite, d ε,s ii ≤ 2. Moreover, by Taylor's theorem, there exists a constant c 1 = c 1 (κ, d) > 0 such that where the last inequality holds by Taylor's theorem. Integrating over s ∈ L 1 × L 2 , we have the result. Turning now to the case that L 1 and L 2 lie on a common line, we may choose v 1 = v 2 =: v. Remark that, by Taylor's theorem, there exists a constant c 1 = c 1 (κ, d) > 0 (i.e. independent of s) and a C 1 function f s : R + → R satisfying, for each x ∈ (0, 1), and Applying the computation in Proposition 14, there exists a c 2 = c 2 (κ, d) > 0, such that, for each ε ∈ (0, 1), Integrating over s ∈ L, we have the result.

Triple-crossings.
Fix d > 0 and an angle μ ∈ [0, π). Let L := (L i ) 1≤i≤3 be a triple of line segments in R 2 contained within the ball B(d). For each ε > 0, define the random variable Remark that N ε counts triple-crossings that are not too degenerate; such triple-crossings pose problems for our asymptotic analysis, and we prefer to deal with these in other ways.  (κ, , d, μ, L , δ) > 0 such that, for each ε ∈ (0, 1), Remark 8. Unlike the condition θ + (s) ≤ μ, which is crucial, we do not believe that the constraint θ − (s) ≥ ε 3/2 is necessary for Proposition 8 to hold; we impose it because it greatly simplifies the proof (see in particular the computation in Proposition 15). On the other hand, the result would be too weak for our purposes if we chose any stronger constraint; see Proposition 11.
Proof. As in the proof of Proposition 7, by a linear rescaling of the plane we may reduce to the case that κ vv (0) = −2 for all unit vectors v ∈ S 1 . We may also choose ε sufficiently small so that, for each s = (s 1 , s 2 , s 3 ) ∈ εL with {s i } 1≤i≤k distinct and not co-linear, the distribution of Ψ (s) is non-degenerate; this is possible by Proposition 2.

Critical points near a level.
We now use the two-dimensional Kac-Rice formula to give a bound on the number of critical points for which the Ψ has a value that is near the level . Let D be a bounded domain. For each s, δ > 0, define the random variable N s,δ := |{x ∈ s D : ∇Ψ (x) = 0, |Ψ (x) − | < δ}| .

Perturbation Analysis
In this section we develop a perturbation analysis that will be crucial in our study of the nodal set of the random plane wave. This section is completely deterministic, and may be of independent interest, although a reader that is not concerned with the special case of the random plane wave may prefer to skip ahead. Define a regular conic section C to be any one of the following objects: (i) the empty set; (ii) a line; (iii) a circle; or (iv) a hyperbola whose principal axes intersect at rightangles. The aim of this section is to give sufficient criteria under which a solution to the Helmholtz equation (3) has a nodal set that is locally well-approximated by a regular conic section.
First we define the sense in which we consider a set to be 'well-approximated' by a regular conic section (the results that follow actually hold in a stronger sense, but this is sufficient for our purposes). For each δ > 0 and set S ⊆ R 2 , define the δ-thickening of S to be the closed set of all points within a distance δ from S. For each δ > 0 and domain D ⊆ R 2 , we say that a set S ⊆ R 2 satisfies the δ-conic property in D if there exists a regular conic section C and a homeomorphism h : D → D mapping S ∩ D onto C ∩ D such that: (i) S ⊆ C δ ; (ii) |S ∩ ∂ D| = |C ∩ ∂ D|; and (iii) |h(s, 1) − s| < δ for each s ∈ S ∩ ∂ D.
See Fig. 8 for an illustration of the δ-conic property. Remark that it requires both that S lies inside a δ-thickening of a regular conic section, and also that the boundary points of S are close to the boundary points of the conic section. Proof. Without loss of generality we may assume that f (x, y) = y +d for some constant d ≥ 0. Remark that the nodal set N is contained in the δ-thickened line If d > 1 + δ then |C δ ∩ B(1)| = 0, and so it suffices to consider the case d ≤ 1 + δ and to set D := B(1 + 4δ). Remark that, for sufficiently small δ, the set C δ ∩ ∂ D consists of two disjoint arcs, and moreover Hence h is not stationary along either of the arcs in C δ ∩ ∂ D, and so N ∩ ∂ D contains exactly one point in each arc. Since these arcs have length at most cδ for some c > 0, the cδ-conic property is satisfied.

Lemma 8 (Nodal set of a perturbed quadratic function: The elliptic case)
. Suppose that f is quadratic with Hessian matrix ∇ 2 that has eigenvalues λ 1 and λ 2 of the same sign such that |λ 1 | = |λ 2 | = 2. Let δ 2 := g C 1 (B(2)) . Then there exists a c > 0 such that, for sufficiently small δ > 0, there is a simply-connected domain D with B(1) ⊆ D ⊆ B(1 + cδ) such that N satisfies the cδ-conic property in D.
Proof. Without loss of generality we may assume that for some x 0 , y 0 ∈ R 2 and d ∈ R. Remark that the nodal set N is contained in the set which is contained in a δ-thickened circle for sufficiently small δ. If d < 2δ 2 , then C δ is contained in a ball of radius 2δ, and so at least one of D = B(1) or D = B(1 + 4δ) is such that |C δ ∩ ∂ D| = 0. Suppose instead that d ≥ 2δ 2 . Then, on the set C δ , we have for sufficiently small δ > 0, |∇ f | ≥ 2 2δ 2 − δ 2 = 2δ > δ 2 =: g C 1 (B(2)) .
Hence we may find a D such that |N ∩ ∂ D| ≤ 2, by either taking D = B(1) or the smallest set that contains B(1) such that ∂ D passes through C δ along gradient flow lines of f . Since these gradient flow lines have length at most cδ for some c > 0, such a set is in B(1 + cδ). Moreover, since N ∩ ∂ D contains at most one point in each gradient flow line that is traversed in this procedure, the cδ-conic property is satisfied. Proof. Without loss of generality we may assume that for some x 0 , y 0 ∈ R 2 and d ≥ 0. Remark that the nodal set N is contained in the set The upshot of Lemmas 7-9 is the following general criteria for a C 3 function Ψ : R 2 → R to have a nodal set that locally satisfies the δ-conic property for some small δ.

Corollary 1.
There exists a c > 0 such that, for sufficiently small δ > 0, the following holds. Let Ψ : R 2 → R be a C 3 function, with N its nodal set. Let λ 1 and λ 2 denote the eigenvalues of the Hessian matrix ∇ 2 Ψ (0). Suppose that, for some ε > 0, either or λ 1 and λ 2 are both non-zero and  (x, y). Remark that the nodal set of h is ε −1 N . Remark also that f is linear with |∇ f (0)| = 1, and further, by Taylor's theorem, By Lemma 7, and after rescaling by ε, we deduce the result. Assume instead that the second condition holds. Without loss of generality we may assume the eigenvectors corresponding to λ 1 and λ 2 are in the direction of the x and y axes respectively, and that |λ 1 | ≥ |λ 2 |. Let μ := √ |λ 1 /λ 2 |, and note that the second condition together with the inequality and let f (x, y) and g(x, y) be respectively the second-order Taylor polynomial of h and its remainder. Then, f is quadratic whose Hessian matrix has eigenvalues μ 1 , μ 2 satisfying |μ 1 | = |μ 2 | = 2. Moreover, by Taylor's theorem, Writing Proof. This is clear from the definition of the δ-conic property.

Application of the perturbation analysis to plane waves.
We now show how to apply the above perturbation analysis to plane waves, i.e. solutions to the Helmholtz equation (3). In particular, we deduce sufficient conditions under which the nodal set locally satisfies the δ-conic property for some small δ > 0.
We begin by explaining the heuristics of the approach. Let Ψ : R 2 → R satisfy Eq. (3) for frequency parameter k > 0, with nodal set N . Let λ 1 and λ 2 be the eigenvalues of the Hessian matrix ∇ 2 Ψ (0). Corollary 1 suggests that, outside an error event on which both ∇Ψ (0) and the eigenvalues λ 1 and λ 2 are very small (relative to the C 3 norm of Ψ near the origin), the nodal set of Ψ satisfies the δ-conic property as long as the ratio |λ 1 /λ 2 | is close to one. To make the link with plane waves, observe in addition to the above that unless Ψ (0) is also small, the nodal set will be empty. On the other hand, if Ψ (0) is small, then since (3) implies that Ψ (0) is proportional to λ 1 + λ 2 , the ratio |λ 1 /λ 2 | must be close to one (outside the aforementioned error event that λ 1 and λ 2 are both very small), and so the δ-conic property holds. Formalising the above, we deduce the following result.
Corollary 2 (Sufficient conditions for nodal set of a plane wave to be locally wellapproximated by a regular conic section). There exists a c = c(k) > 0 such that, for sufficiently small δ > 0, the following holds. For each ε > 0, if the conditions and

do not all hold, then there exists a simply-connected domain D with B(ε) ⊆ D ⊆ B((1 + cδ)ε) such that N satisfies the cδε-conic property in D.
Proof. Let c 0 > 0 denote the constant appearing in Corollary 1 and fix δ > 0 sufficiently small that the conclusion of that corollary is available.
Choosing c > c 2 , this implies that λ 1 and λ 2 are non-zero and of opposite sign. Hence using the inequality 1 − |x| |y| ≤ |x + y| max{|x|, |y|} − |x + y| , valid for non-zero x and y of opposite sign, the second condition of Corollary 1 is then satisfied if c is chosen large enough.

Proof of the Main Results
In this section we complete the proofs of the main results. The variables κ, , L, D and δ are fixed throughout this section. Note that κ is assumed to satisfy the conditions in Assumption 1, and the domain D is always considered to be bounded, but it is not necessarily smooth unless explicitly stated. Throughout this section, whenever an edge e ∈ εE is considered, f 1 , f 2 ∈ εF are defined to be the faces adjacent to e.

Control of the 'bad' events.
We first show how to control each of the 'bad' events introduced in Sect. 2 using the results developed in Sects. 3-5. As explained in that section, for the proofs of Theorems 2, 3 and 5 we need only control the existence of these events. On the other hand, for the proof of Theorem 6 we shall also need to control their multiplicity. Proof. This is a direct application of Proposition 7 together with Markov's inequality.
Four-crossings and tubular-crossings. We control four-crossings and tubular-crossings by linking these to simpler crossing events which we treat with Kac-Rice formulae. In particular, we show that four-crossings and tubular-crossings necessarily imply either (i) non-linear crossings of three line segments, or (ii) crossings that are extremely close to lattice vertices or to each other; we control the former using Proposition 8, and the latter using Proposition 7.
For each ε > 0 and edge e ∈ εE, let T (e) be the event that e has a tubular crossing (as defined in Sect. 2), and define the multiplicity of the tubular-crossing of e, denoted T (e), to be the number of components of f 1 ∪ f 2 \ N that lie between the endpoints of e (i.e. two less than the smallest number of components crossed by any curve joining the endpoints).
Moreover, for each ε > 0 and face f ∈ εF, let F( f ) be the event that f has a four crossing (as defined in Sect. 2), and define the multiplicity of the four-crossing to bẽ Lemma 11 (Four-crossings and tubular-crossings imply non-degenerate triplecrossings). Suppose that L is strictly-convex, i.e. the interior angles of the faces F are strictly less than π . Then, there exists an angle μ = μ(L) ∈ [0, π) and a constant c = c(L) > 0 such that, for any δ ∈ (0, 1) and any e ∈ E and f ∈ F: Proof. This can be proved by combining Lemma 11 (setting δ := ε 3/2 ) with the bounds in Proposition 8 (to control cases (1a) and (2a)) and Proposition 7 (to control cases (1b), (2b) and (2c)), and applying Markov's inequality. Notice that the choice of δ := ε 3/2 is the maximum possible that will give the correct order O(ε 4 ) in (1b) and (2b)). It is for this reason that we are able to impose the constraint θ − (s) ≥ ε 3/2 in our statement of Proposition 8 (which eases the calculations), but no stronger constraint is possible.

Small excursion domains.
We control small excursion domains by using the derivative bounds in Lemma 6 and Proposition 9. For each ε > 0 and bounded domainD, define a small excursion domain inD to be a component ofD \ N whose boundary intersects at most one edge in εE and does not intersect ∂D. For each s > 0, letS(s, ε) to be the number of small excursion domains in s D.
Proof. Let S ⊆ s D be a small excursion domain in s D. Observe that there must be a local maxima or minima contained within S; suppose that x ∈ S is such a point. Let f ∈ εF denote the face that contains x, and F denote the union of f with all its adjacent faces in εF. By Taylor's theorem, each y ∈ F satisfies Since S is a small excursion domain, the level set N intersects F, and so we may select a y ∈ F such that Ψ (y) = . We conclude that x ∈ S satisfies which implies the result.
where in the final step we conditioned on the σ -algebra generated by Ψ C 2 (1+2d(L)) and used Proposition 9. Applying Lemma 6 together with Markov's inequality gives the result. Invisible errors. Finally, we control invisible errors by combining the deterministic analysis in Corollary 2 with the bounds in Lemma 6 and Proposition 3. As for Sect. 5, this subsection treats only the special case of the random plane wave, and may be skipped by a reader not concerned with this case.
In this subsection we assume that κ satisfies (2) for a frequency parameter k > 0, L is a regular hexagonal lattice, and the level is set at = 0. Without loss of generality we may rescale L so that d(L) = √ 13/4 ≈ 0.9 (this ensures that a circle of unit radius exactly circumscribes two adjacent hexagons).
First we give a formal definition of the 'invisible errors' that we consider. Recall the signed regions S + and S − defined in (4). Recall also the definition of a Type 2 (see in particular Fig. 4). Let ε > 0 and consider an edge e ∈ εE. LetD be a domain such that f 1 ∪ f 2 ⊆D. We say that e gives rise to an invisible error inD if the faces f 1 and f 2 do not display a Type 2 error pattern, and either: (i) there is a component ofD \ N that separates the endpoints of e inD; or (ii) there is a component of ( f 1 ∪ f 2 ) \ N that separates the endpoints of e in f 1 ∪ f 2 , and moreover this component intersects ∂D and at least one vertex in f 1 ∪ f 2 . See Fig. 10 for examples of these two cases.
Before discussing how to control invisible errors, we state a result that links them to the geometry of the regular hexagonal lattice, by way of the δ-conic property; this is the only time the specific geometry of regular hexagons is needed (with one small caveat mentioned below). Proof. This is clear from inspection; see Fig. 11.
Corollary 3 (The δ-conic property precludes invisible errors). There exists a sufficiently small δ > 0 such that, for each ε > 0 and each e ∈ εE, if N satisfies the δε-conic property in a simply-connected domainD with f 1 ∪ f 2 ⊆D, then e does not give rise to an invisible error. Proof. This can be deduced from Lemma 13; see Fig. 12. The most delicate case is when the regular conic section is a pair of lines intersecting at right-angles, in which case Lemma 13 applies directly. Note that we may assume without loss of generality that ε = 1.
For each ε > 0 and edge e ∈ E, let m e denote the mid-point of e. Let a > 1 be a constant such that (i) the ball B(m e , aε) does not contain any other vertices in V \ ( f 1 ∪ f 2 ); and (ii) 2a < 3d(L). It is easy to check that the first condition is possible (note that this uses a specific property of the regular hexagonal lattice, but one that is also shared by the square and regular triangular lattices), and the second condition is the origin of the bound 3d(L) in the definition of ε-homeomorphic in Sect. 1). Let I (e) be the event that e gives rise to an invisible error in every domainD such that B(ε) ⊆D ⊆ B(aε).
Proof. Assume without loss of generality that m e is the origin. Let c 1 = c 1 (k) > 0 be a constant satisfying Corollary 2, and select δ > 0 sufficiently small such that, if N satisfies the c 1 δ-conic property in a simply-connected domainD satisfying f 1 ∪ f 2 ⊆D, then e does not give rise to an invisible error. This is possible by Corollary 3. Also choose δ small enough to ensure that 1 + c 1 δ < a, where a > 1 is the constant used in the definition of the event I (e).

Proof of the main results.
We complete the proof of the main results by combining our control on the 'bad' events in Propositions 10, 11, 12 and 13 with simple topological arguments.
Proof of Theorem 2. Observe first that we may assume that L is strictly-convex, since if it is not we may simply apply the result to a strictly-convex graph L * formed by adding additional edges to L, and deduce the result. Note that the number of edges in s D∩εE is of order ε −2 s 2 . Hence, by Propositions 10, 11 and 12 and the union bound, there exists a c = c(κ, , L, D) such that, for each ε ∈ (0, 1) and s > 0, there is an event X of probability 1 − cs 2 ε on which: (i) no edge in s D ∩ εE has a double-crossing; (ii) no face in s D ∩ εF has a four-crossing; and (iii) there are no small excursion domains inside s D. We henceforth assume that the event X holds.
Consider a face f ∈ P ε (s D)∩εF. Since we are working on X , none of the boundary edges e ∈ ∂ f has a double-crossing, nor does f have a four-crossing or a small excursion domain. Hence the level set N | f is either empty or consists of exactly one curve which joins distinct edges of f . In either case, the level sets N ε and N are ε-homeomorphic in f , as required.
Remark that, by this argument, we extract a bound of cεs 2 on the rate of convergence in the statement of Theorem 2, for some c = c(κ, , L, D). Proof of Theorem 3. As in the proof of Theorem 2, we may assume that L is strictlyconvex. Note that, although the number of edges in s D ∩ εE is of order ε −2 s 2 , since the domain D is smooth the number of boundary edges in ∂ P ε (s D) is only of order ε −1 s. Hence by Propositions 10, 11 and 12 and the union bound, there exists a c = c(κ, , L, D, δ) such that, for each ε ∈ (0, 1) and s > 0, there is an event X of probability 1 − ε 2−δ s 2 on which: (i) no boundary edge in ∂ P ε (s D) ∩ εE has a double-crossing; (ii) no face in P ε (s D) ∩ εF has either a tubular-crossing or a four-crossing; and (iii) there are no small excursion domains inside P ε (s D). We henceforth assume that the event X holds.
As in the proof of Theorem 2, if a face f ∈ P ε (s D) ∩ εF does not have a boundary edge e ∈ ∂ f with a double-crossing, then the level sets N ε and N are ε-homeomorphic in f . So consider an edge e ∈ P ε (s D) ∩ εE that has a double-crossing, and let f 1 , f 2 ∈ εF denote the adjacent faces.
Since we are working on X , any such edge is: (i) in the interior of P ε (s D); and (ii) does not have a tubular-crossing. This implies the existence of a pair of adjacent points Since there are no small excursion domains, there therefore exists a homeomorphism h : f 1 ∪ f 2 → f 1 ∪ f 2 mapping the level line c onto a second curveĉ such that (i) c andĉ have the same endpoints on ∂( f 1 ∪ f 2 ), and (ii)ĉ crosses the edge e two fewer times than c. Moreover, this homeomorphism is local, in the sense that it moves points at most 2d(L)ε.
Sequentially applying such homeomorphisms to eliminate all double-crossing from e (see Fig. 13), and then repeating the procedure for each e ∈ P ε (s D) ∩ εF that has a double-crossing, we see that the level sets N ε and N are ε-homeomorphic in s D, as required.
Remark that, by this argument, we extract a bound of cε 2−δ s 2 on the rate of convergence in the statement of Theorem 3, for some c = c(κ, , L, D, δ).
Proof of Theorem 5. Without loss of generality we may rescale L so that d(L) = √ 13/4 (since L is regular, this is equivalent to rescaling the mesh-size ε > 0).
Assume first that ε > 0 is sufficiently small that nodal sets intersect at least two adjacent faces (i.e. there are no small nodal sets). Similarly to in the proofs of Theorems 2 and 3, by Propositions 10 and 13 and the union bound, there exists a c = c(k, D) such that, for each s > 0 and ε > 0, there is an event X of probability 1 − cε 2 s on which: (i) no boundary edge in ∂ P ε (s D) ∩ εE has a double-crossing; (ii) no edge e ∈ s D ∩ εE satisfies the invisible error event I (e). Note that the error of order ε 2 s arise from the control on double-crossings on the boundary; the control on invisible errors gives a term of order ε 4 s 2 . Henceforth assume that the event X holds.
As in the proof of Theorem 2, the level sets N ε and N are ε-homeomorphic in f unless a face f ∈ P ε (s D) ∩ εF either does not have a four-crossing, or does not have an edge e ∈ ∂ f with a double-crossing. We consider each of these cases in turn.
Consider an edge e ∈ P ε (s D) ∩ εE that has a double-crossing, and let f 1 , f 2 ∈ εF denote the adjacent faces. Since we are working on X , any such edge is (i) in the interior of P ε (s D), and (ii) does not satisfy the invisible error event I (e). Using the definition of the event I (e), and since there are no nodal domains, there are now two possibilities.
Either f 1 and f 2 display a Type 2 error, in which case they are contained in a component of the set of visible ambiguities π ε . Or, there must exist a domainD ⊆ P ε (s D) satisfying 2d(D) < 3ε in which a homeomorphism can be constructed that eliminates all doublecrossing of e (in the same sense as in the proof of Theorem 3) and that fixes the vertices Consider next a face f ∈ P ε (s D) ∩ εF with a four-crossing, but whose boundary edges do not have double-crossings. Then f must display a Type 1 error pattern, and hence is also contained in a component of π ε .
Combining the above, we have shown that the level sets N ε and N are εhomeomorphic up to resolutions of the set of visible ambiguities, as required. Remark that, by this argument, we extract a bound of cε 2 s on the rate of convergence in the statement of Theorem 5, for some c = c(k, L, D).
Proof of Theorem 6. As in the proof of Theorem 2, we may assume that L is strictlyconvex. Recall that P ε (s D) is the largest L-compatible set that is contained in s D, and let P ε (s D) be the smallest L-compatible set that contains s D. Recall also that N (s D) and N ε (s D) denote, respectively, the number of components of s D ∩ N and P ε (s D) ∩ N ε .
Observe that discrepancies in N (s D) and N ε (s D) come from three sources: (i) components of s D \ N that do not intersect any v ∈ s D ∩ εV; (ii) vertices v ∈ s D ∩ εV that belong to the same component of s D \ N but different components of s D \ N ε ; and (iii) vertices v ∈ s D ∩ εV that belong to the same component of s D \ N ε but different components of s D \ N .
The first of these is bound above by the number of small excursion domains in P ε (s D), the total multiplicity of double-crossings on the boundary ∂s D, and the total multiplicity of tubular-crossings in s D; the second is bound by the total multiplicity of tubular-crossings in s D; the third is bound by the total multiplicity of four-crossings in s D.
All in all, defining the random variables we see that As in the proof of Theorem 3, by Propositions 10, 11 and 12 and the union bound, there exists a c = c(κ, , L, D) > 0 such that, for ε ∈ (0, 1) and s > 0, Putting this together, we have

A. Appendix A: Proof of the Gaussian Estimates
In this appendix we give the proof of Propositions 2 and 3.
Proof of Proposition 2. Let ρ denote the spectral measure of Ψ , defined by the relation By Assumption 1 (as well as the remarks that follow), ρ is not supported on a line. Hence by linear rescaling we may assume, without loss of generality, that We next argue that it is sufficient to prove the result for the fieldΨ whose spectral measureρ is the sum of unit point-masses at each point in P. To see why, assume the claimed result fails for a given field Ψ , i.e. for each δ > 0 there exist distinct non-colinear points (s i ) i=1,2,3 inside a ball of radius δ such that the matrix K := (κ(s j − s k )) 1≤ j,k≤3 is singular. For each such set of points, the singularity of K is equivalent to the existence of non-zero ( which in turn is equivalent to supp(ρ) being contained in the zero set of g s (μ) := 3 j=1 c j e 2πi s i ,μ . Hence, since supp(ρ) ⊆ supp(ρ), the same is true for supp(ρ), and so the claimed result must also fail forΨ .
Finally we show that the result is true forΨ by direct computation. By stationarity, it is enough to show the existence of a δ > 0 such that the distribution of is non-degenerate for all distinct non-zero s 1 := (x 1 , y 1 ) and s 2 := (x 2 , y 2 ) in [0, 1] 2 ∩ B(δ) such that s 1 and s 2 are not co-linear with the origin; we actually prove the stronger statement that (8)  the result follows once we show that, for λ > 1, the function is strictly convex on the domain 0 ≤x < 1/λ 2 . This may be checked by direct computation, since Proof of Proposition 3. By Lemma 4, the two components of ∇Ψ (0) and the random vector (λ 1 , λ 2 ) are independent. Moreover, the distribution of ∇Ψ (0) is non-degenerate and so there exists a c > 0 such that Hence it is sufficient to show that there exists a c such that We use the fact that ∇ 2 Ψ (0) is distributed as where M is a 2×2 symmetric random matrix with all elements centred Gaussian random variables with covariance with δ i j the Kronecker delta function (see [9,Eq. (5)], for which it is useful to note that, in the notation of that paper (see, especially, Eq. (3)), J 2 = 2k 4 J (4) 0 (0)/3 = k 4 /4 in the case of the random plane wave, where J (4) 0 is the fourth derivative of the zeroth Bessel function). The density of the (ordered) eigenvalues μ 1 ≤ μ 2 of the matrix M is known [8,Theorem 2.2], and satisfies for a certain normalising constant c 1 > 0. Hence there exists a c 2 > 0 such that the density of λ 1 and λ 2 satisfies Integrating over the region {|λ 1 | < δ 2 , |λ 2 | < δ 2 , |λ 1 + λ 2 | < δ 3 } yields the result.

B. Appendix B: Kac-Rice Computations
In this appendix we give details of the matrix computations that we used in the Kac-Rice arguments in Sect. 4 above. These are used to control the effect of conditioning on certain Gaussian vectors involving the random field Ψ and its derivatives Ψ v .
For each x > 0, define the matrices a matrix C = (c i j ) 1≤i, j≤1 satisfying c 11 = c 22 = 2, the vector and the matrix Then there exists a constant c = c(δ) > 0 such that, for each x ∈ (0, 1), Proof. By explicit calculation, and Putting the bounds on f and g into the above formulae yields that, as x → 0, and Combining these gives the result.
Proposition 15. Fix a constant δ > 0, an angle μ ∈ (0, π), and a set of unit vectors For each 1 ≤ i = j ≤ 3, denote by η i j the angle between the line segment s i s j and the vector v j , let e s i j ∈ R be such that and let f s i j , g s i j : R + → R be such that, for each x ∈ (0, 1), f s i j (x) − 1 + x 2 + e s i j x 4 < δx 6 and g s i j (x) + cos(η i j ) 2x + 4e s i j x 3 < δx 5 .
For i = 1, 2, 3, let α i denote the angle in triangle s at vertex s i , and let the opposing side-length be εa i , for some a i ∈ (0, 1]. Without loss of generality we may assume that a 3 = 1, and notice that max{a 1 , a 2 } ≥ 1/2 by the triangle inequality. We may also assume, without loss of generality, that α 1 ≤ α 2 and α 1 ≤ α 3 .
Notice that, depending on the orientation of s 2 s 3 and v 1 , Hence because of the identity sin 2 (θ ) + 2 cos(α) cos(β) cos(θ ) − cos 2 (α) − cos 2 (β) = 0, valid for each pair of angles α and β, and any third angle which yields the result in the case that the angle α 1 is uniformly bounded away from zero.
So it remains to deal with the case in which angle α 1 is not uniformly bounded away from zero, but may instead be as small as ε 3/2 (this constraint is crucial in restricting the order to which we need to control errors in the Taylor expansions below). Suppose then that α 1 < μ/2. Then, since θ + (s) ≤ μ, and as the angles in the triangle sum to π , it must be the case that α i ≥ (π − μ)(1 − δ) ≥ (π − μ)/2 for each i = 2, 3. Using the sine rule and the bound sin(x) < x for x > 0, it follows that there exist a c 1 = c 1 (μ) such that a 1 < c 1 α 1 .
By the triangle inequality, it must also be the case that there exists a c 2 = c 2 (μ) such that |a 2 − a 3 | < c 2 α 1 .
Combining these gives the result.

C. Appendix C: Proof of the Russo-Seymour Welsh Estimates
In this appendix we give the proof of Theorem 4; for this, it is helpful to have read both [4,15]. We follow the proof in [4], but make three essentially distinct improvements to the arguments. The main improvement is to replace the discretisation scheme in [4,Theorem] with our enhanced scheme in Theorem 3; this is already enough to reduce the required decay exponent from α ≈ 325 to α ≈ 55. The other two improvements are relatively minor. Before explaining them, we state a preliminary proposition on the quasi-independence of signs of Gaussian vectors, which is a strengthened form of [4,Theorem 4.3]. This bound is completely independent of other arguments, and may be of independent interest. where (a i ) m i=1 and (b i ) n i=1 are vectors consisting of independent standard normal random variables, and (λ 1,i ) m i=1 and (λ 2,i ) n i=1 are the eigenvalues, in descending order, of Σ X and Σ Y respectively, with (u 1,i ) m i=1 and (u 2,i ) n i=1 the corresponding orthonormal eigenvectors. Fix a threshold λ > 0 to be set later, and truncate X and Y as X =X +X = where m 0 and n 0 are chosen so that λ 1,m 0 ≥ λ > λ 1,m 0 +1 and λ 2,n 0 ≥ λ > λ 2,n 0 +1 . The first step is to control the probability that X andX (respectively Y andỸ ) have the same signs. Set ε 1 := √ c 0 λm and ε 2 := √ c 0 λn for a constant c 0 > 1 to be determined later. Notice that by the union bound

C.1. Quasi-independence of the signs of Gaussian vectors.
On the other hand, applying Parseval's theorem is an upper bound on the probability that Y andỸ do not have the same signs. The second step is to control the total variation distance between (X ,Ỹ ) and the vector (X ,Ȳ ) formed from independent copies ofX andỸ ; on the event that X and X have the same signs, this gives an upper bound on the quantity in (12). By linear transformation, this total variation distance is the same as the total variation distance between the vector Z = (a 1 , . . . , a m 0 , b 1 , . . . , b n 0 ) and a vector of m 0 + n 0 independent standard normal random variables. By Pinsker's inequality, this distance is bound above by C.2. Proof of Theorem 4. Before we begin, we remark that the overall decay exponent of α = 16 arises essentially from three distinct sources: it can be understood as α = d(1 + a)b, where d is the dimension of R 2 , ε = s −a is the required mesh-size in the discretisation scheme (in our case a = 1 by Theorem 3), and η = m −b is the maximum size of the cross-correlation terms between two Gaussian vectors of size m in order to guarantee quasi-independence of the law of the signs (in our case b = 4 by Proposition 16). As such, any improvement in the exponents in either Theorem 3 or Proposition 16 would lower this decay exponent.
Assume that κ ≥ 0 and that κ is invariant under reflection through the horizontal axis and under rotation by π/2. Fix a periodic lattice L that is a triangulation, is invariant under reflection through the horizontal axis and under rotation by π/2, and is integral in the sense that V ⊆ ( 1 N Z) 2 for some N ∈ N. Fix a sequence ε = ε(s), a smooth bounded domain D and disjoint boundary arcs γ and γ on ∂ D.
Consider the discretised excursion domain s D ∩ N ε , and observe that there is a natural way to partition this set into a To interpret the constituent parts of (14), the term ε = ε 1+log(3/2)/ log(4/3) is an adjustment to the lattice size to ensure that Tassion's argument in [15] applies, the term m = sε −1 2 is the order of the number of vertices of εL inside s D, the term η = s −α gives the order of correlations on the scale s, and finally the relationship is a sufficient condition, by [4,Theorem 4.3], to ensure the asymptotic quasi-independent of the signs of εL ∩ s D 1  Our first improvement to the decay exponent is achieved by a direct substitution of the discretisation scheme in [4,Theorem 1.5] with the scheme in Theorem 3, which is valid if ε = o(s −1−δ ) and under Assumption 1. Then it is immediate that (17) holds as long as κ(x) = o(|x| −α−δ ) for α > (1 + 1 × (1 + log(3/2)/ log(4/3))) × 2 × 8 ≈ 55.
We then observe that [15, Lemma 2.2] still holds, with an identical proof, if we replace the constant 2/3 with any ρ ∈ (0, 1), and so the exponent c can be replaced with 1 + δ 1 for any δ 1 > 0 by taking ρ close enough to one.
Second, we substitute the strengthened quasi-independence result in Proposition 16 in place of [4,Theorem 4.3]. The consequence is that the exponent in (16) may be reduced from 8 to 4 + δ 2 , for any δ 2 > 0.