Phase Singularities in Complex Arithmetic Random Waves

Complex arithmetic random waves are stationary Gaussian complex-valued solutions of the Helmholtz equation on the two-dimensional flat torus. We use Wiener-It\^o chaotic expansions in order to derive a complete characterization of the second order high-energy behaviour of the total number of phase singularities of these functions. Our main result is that, while such random quantities verify a universal law of large numbers, they also exhibit non-universal and non-central second order fluctuations that are dictated by the arithmetic nature of the underlying spectral measures. Such fluctuations are qualitatively consistent with the cancellation phenomena predicted by Berry (2002) in the case of complex random waves on compact planar domains. Our results extend to the complex setting recent pathbreaking findings by Rudnick and Wigman (2008), Krishnapur, Kurlberg and Wigman (2013) and Marinucci, Peccati, Rossi and Wigman (2016). The exact asymptotic characterization of the variance is based on a fine analysis of the Kac-Rice kernel around the origin, as well as on a novel use of combinatorial moment formulae for controlling long-range weak correlations.


Overview and main results
Let T := R 2 /Z 2 be the two-dimensional flat torus, and define ∆ = ∂ 2 /∂x 2 1 + ∂ 2 /∂x 2 2 to be the associated Laplace-Beltrami operator. Our aim in this paper is to characterize the high-energy behaviour of the zero set of complex-valued random eigenfunctions of ∆, that is, of solutions f of the Helmholtz equation ∆f + Ef = 0, (1.1) for some adequate E > 0. In order to understand such a setting, recall that the eigenvalues of −∆ are the positive reals of the form E n := 4π 2 n, where n = a 2 + b 2 for some a, b ∈ Z (that is, n is an integer that can be represented as the sum of two squares). Here, and throughout the paper, we set S := {n ∈ N : a 2 + b 2 = n, for some a, b ∈ Z}, and for n ∈ S we define Λ n := {λ = (λ 1 , λ 2 ) ∈ Z 2 : λ 2 := λ 2 1 + λ 2 2 = n} to be the set of energy levels associated with n, while N n := |Λ n | denotes its cardinality. An orthonormal basis (in L 2 (T)) for the eigenspace associated with E n is given by the set of complex exponentials {e λ : λ ∈ Λ n }, defined as e λ (x) := e i2π λ,x , x ∈ T, with i = √ −1. For every n ∈ S, the integer N n =: r 2 (n) counts the number of distinct ways of representing n as the sum of two squares: it is a standard fact (proved e.g. by using Landau's theorem) that N n grows on average as √ log n, and also that there exists an infinite sequence of prime numbers p ∈ S, p ≡ 1 mod 4, such that N p = 8. A classical discussion of the properties of S and N n can be found e.g. in Section 16.9 and 16.10]. In the present paper, we will systematically consider sequences {n j } ⊂ S such that N n j → ∞ (this is what we refer to as the high-energy limit).
The complex waves considered in this paper are natural generalizations of the real-valued arithmetic waves introduced by Rudnick and Wigman in [R-W], and further studied in [K-K-W, M-P-R-W, O-R-W, R-W2]; as such, they are close relatives of the complex fields considered in the physical literature -see e.g. Be3,N,, as well as the discussion provided below. For every n ∈ S, we define the complex arithmetic random wave of order n to be the random field Θ n (x) := 1 √ N n λ∈Λn v λ e λ (x), x ∈ T, (1.2) where the v λ , λ ∈ Λ n , are independent and identically distributed (i.i.d.) complex-valued Gaussian random variables such that, for every λ ∈ Λ n , Re(v λ ) and Im(v λ ) are two independent centered Gaussian random variables with mean zero and variance one 1 . The family {v λ : λ ∈ Λ n , n ∈ S} is tacitly assumed to be defined on a common probability space (Ω, F , P), with E indicating expectation with respect to P. It is immediately verified that Θ n satisfies the equation (1.1), that is, ∆Θ n + E n Θ = 0, and also that Θ n is stationary, in the sense that, for every y ∈ T, the translated process x → Θ n (y + x) has the same distribution as Θ n (this follows from the fact that the distribution of {v λ : λ ∈ Λ n } is invariant with respect to unitary transformations; see Section 1.2 for further details on this straightforward but fundamental point).
We will show below (Part 1 of Theorem 1.2) that, with probability one, I n is a finite collection of isolated points for every n ∈ S; throughout the paper, we will write I n := |I n | = Card(I n ), n ∈ S. (1.4) In accordance with the title of this work, the points of I n are called phase singularities for the field Θ n , in the sense that, for every x ∈ I n , the phase of Θ n (x) (as a complex-valued random quantity) is not defined.
As for nodal lines of real arithmetic waves [K-K-W, M-P-R-W], our main results crucially involve the following collection of probability measures on the unit circle S 1 ⊂ R 2 : µ n (dz) := 1 N n λ∈Λn δ λ/ √ n (dz), n ∈ S, (1.5) as well as the associated Fourier coefficients µ n (k) := S 1 z −k µ n (dz), k ∈ Z. (1.6) In view of the definition of Λ n , the probability measure µ n defined in (1.5) is trivially invariant with respect to the transformations z → z and z → i · z. The somewhat erratic behaviours of such objects in the high-energy limit are studied in detail in [K-K-W, K-W].
Here, we only record the following statement, implying in particular that the sequences {µ n : n ∈ S} and { µ n (4) : n ∈ S} do not admit limits as N n diverges to infinity within the set S.
Recall from [K-K-W, K-W] that a measure µ on (S 1 , B) (where B is the Borel σ-field) is said to be attainable if there exists a sequence {n j } ⊂ S such that N n j → ∞ and µ n j converges to µ in the sense of the weak-⋆ topology.
Proposition 1.1 (See [K-W, K-K-W]). The class of attainable measures is an infinite strict subset of the collection of all probability measures on S 1 that are invariant with respect to the transformations z → z and z → i · z. Also, for every η ∈ [0, 1] there exists a sequence {n j } ⊂ S such that N n j → ∞ and | µ n j (4)| → η.
Note that, if µ n j converges to µ ∞ in the weak-⋆ topology, then µ n j (4) → µ ∞ (4). For instance, one knows from [E-H, K-K-W] that there exists a density one sequence {n j } ⊂ S such that N n j → ∞ and µ n j converges to the uniform measure on S 1 , in which case µ n j (4) → 0. (Ω, F , P) will be denoted by law =⇒, whereas equality in distribution will be indicated by the symbol law = .
The main result of the present work is the following exact characterization of the first and second order behaviours of I n , as defined by (1.4), in the high-energy limit. As discussed below, it is a highly non-trivial extension of the results proved in [K-K-W, M-P-R-W], as well as the first rigorous description of the Berry's cancellation phenomenon [Be3] in the context of phase singularities of complex random waves.
3. Relations (1.8)-(1.9) are completely new and are among the main findings of the present paper; in particular, they show that the asymptotic behaviour of the variance of I n is non-universal. Indeed, when N n j → ∞, the fluctuations of the sequence d n j depend on the chosen subsequence {n j } ⊂ S, via the squared Fourier coefficients µ n j (4) 2 : in particular, the possible limit values of the sequence {d n j } correspond to the whole interval 5 128π 2 , 1 16π 2 . As discussed in the sections to follow, such a non-universal behaviour echoes the findings of [K-K-W], in the framework of the length of nodal lines associated with arithmetic random waves. We will see that our derivation of the exact asymptotic relation (1.8) follows a route that is different from the one exploited in [K-K-W] -as our techniques are based on chaos expansions, combinatorial cumulant formulae, as well as on a novel local Taylor expansion of the second order Kac-Rice kernel around the origin. 4. The support of the distribution of each variable J η is the whole real line, but the distribution of J η is not Gaussian (this follows e.g. by observing that law of J η has exponential tails). As discussed in the forthcoming Section 1.2, similar non-central and non-universal second order fluctuations have been proved in [M-P-R-W] for the total nodal length of real arithmetic random waves. We will show below that this striking common feature originates from the same chaotic cancellation phenomenon exploited in [M-P-R-W], that is: in the Wiener chaos expansion of the quantity I n , the projection on the second chaos vanishes, and the limiting fluctuations of such a random variable are completely determined by its projection on the fourth Wiener chaos. It will be clear from our analysis that, should the second chaotic projection of I n not disappear in the limit, then the order of Var(I n ) would be proportional to E 2 n /N n , as N n → ∞.
5. Choosing ǫ n j ≡ ǫ (constant sequence), one deduces from Point 3 of Theorem 1.2 that the ratio I n /n converges in probability to π, whenever N n → ∞. See e.g. [D, p. 261] for definitions.
6. It is easily checked (for instance, by computing the third moment) that the law of J η 0 differs from that of J η 1 for every 0 ≤ η 0 < η 1 ≤ 1. This fact implies that the sequence {|µ n (4)|} dictates not only the asymptotic behaviour of the variance of I n , but also the precise nature of its second order fluctuations.
7. Reasoning as in [M-P-R-W, Theorem 1.2], it is possible to suitably apply Skorohod representation theorem (see e.g. [D,Chapter 11]), in order to express relation (1.11) in a way that does not involve the choice of a subsequence {n j }. We leave this routine exercise to the interested reader.
In the next section, we will discuss several explicit connections with the model of real arithmetic random waves studied in [K-K-

Complex zeros as nodal intersections
For simplicity, from now on we will write T n (x) := Re(Θ n (x)), T n (x) := Im(Θ n (x)), ( 1.13) for every x ∈ T and n ∈ S; in this way, one has that We will also adopt the shorthand notation Our next statement yields a complete characterization of the distribution of the vectorvalued process T n , as a two-dimensional field whose components are independent and identically distributed real arithmetic random waves, in the sense of [K-K- Proposition 1.4. Fix n ∈ S. Then, T n and T n are two real-valued independent centered Gaussian fields such that E T n (x)T n (y) = E T n (x) T n (y) = 1 N n λ∈Λn cos(2π λ, x − y ) =: r n (x − y). (1.14) Also, there exist two collections of complex random variables (1.15) with the following properties: (i) A(n) and A(n) are stochastically independent and identically distributed as random vectors indexed by Λ n ; (ii) For every λ ∈ Λ n , a λ is a complex-valued Gaussian random variable whose real and imaginary parts are independent Gaussian random variables with mean zero and variance 1/2; (iii) If λ / ∈ {σ, −σ}, then a λ and a σ are stochastically independent; (iv) a λ = a −λ ; (v) For every x ∈ T, Proof. One need only show that T n and T n are two centered independent Gaussian fields such that (1.14) holds; the existence of the two families A(n) and A(n) can consequently be derived by first expanding T n and T n in the basis {e λ : λ ∈ Λ n }, and then by explicitly computing the covariance matrix of the resulting Fourier coefficients. Relation (1.14) follows from a direct computation based on the symmetric structure of the set Λ n , once it is observed that T n and T n can be written in terms of the complex Gaussian random variables {v λ } appearing in (1.2), as follows: The fact that r n only depends on the difference x − y confirms in particular that T n is a two-dimensional Gaussian stationary process.
Assumption 1.5. Without loss of generality, for the rest of the paper we will assume that, for n = m, the two Gaussian families are stochastically independent; this is the same as assuming that the two vector-valued fields T n and T m are stochastically independent.
As anticipated, relation (1.16) implies that T n and T n are two independent and indentically distributed real arithmetic random waves of order n, such as the ones introduced in [R -W], and then further studied in [K-K-W, M-P-R-W, O-R-W, R-W2]. We recall that, according to [C], with probability one both T −1 n (0) and T −1 n (0) are unions of rectifiable curves, called nodal lines, containing a finite set of isolated singular points. The following statement yields a further geometric characterisation of I n and I n : its proof is a direct by-product of the arguments involved in the proof of Part 1 of Theorem 1.2. Proposition 1.6. Fix n ∈ S. Then, with probability one the nodal lines of T n and T n have a finite number of isolated intersection points, whose collection coincides with the set I n ; moreover, I n does not contain any singular point for T n .
In view of Proposition 1.6, it is eventually instructive to focus on the random nodal lengths L n := length (T −1 n (0)), n ∈ S, for which we will present a statement collecting some of the most relevant findings from [R-W] (Point 1), [K-K-W] (Point 2) and [M-P-R-W] (Point 3).
As discussed e.g. in [K-K-W, R-W], relations (1.17)-(1.18) yield immediately a law of large numbers analogous to (1.10). We stress that Theorem 1.2 and Theorem 1.7 share three common striking features (explained below in terms of a common chaotic cancellation phenomenon), namely: (a) an inverse quadratic dependence on N n , as displayed in formulae (1.8) and (1.18), (b) non-universal variance fluctuations, determined by the quantities d n and c n defined in (1.9) and (1.19), respectively, and (c) non-universal and non-central second order fluctuations (see (1.11) and (1.20)).
The estimate (1.18) largely improves upon a conjecture formulated by Rudnick and Wigman in [R-W], according to which one should have Var(L n ) = O(E n /N n ). The fact that the natural leading term E n /N n actually disappears in the high-energy limit, thus yielding (1.18), is connected to some striking discoveries by Berry [Be3], discussed in the forthcoming section.

More about relevant previous work
Random waves and cancellation phenomena. To the best of our knowledge, the first systematic analysis of phase singularities in wave physics appears in the seminal contribution by Nye and Berry [N-B]. Since then, zeros of complex waves have been the object of an intense study in a variety of branches of modern physics, often under different names, such as nodal points, wavefront dislocations, screw dislocations, optical vortices and topological charges. The reader is referred e.g. to [D-O-P, N, U-R], and the references therein, for detailed surveys on the topic, focussing in particular on optical physics, quantum chaos and quantum statistical physics.
One crucial reference for our analysis is Berry [Be2], where the author studies several statistical quantities involving singularities of random waves on the plane. Such an object, usually called the (complex) Berry's random wave model (RWM), is defined as a complex centered Gaussian field, whose real and imaginary parts are independent Gaussian functions on the plane, with covariance where E > 0 is an energy parameter, and J 0 is the standard Bessel function (see also [Be1]). Formula (1.21) implies in particular that Berry's RWM is stationary and isotropic, that is: its distribution is invariant both with respect to translations and rotations. As discussed e.g. in [K-K-W, Section 1.6.1], if {n j } ⊂ S is a sequence such that N n j → ∞ and µ n j converges weakly to the uniform measure on the circle, then, for every x ∈ T and using the notation (1.14), showing that Berry's RWM is indeed the local scaling limit of the arithmetic random waves considered in the present paper. Reference [Be3], building upon previous findings of Berry and Dennis [B-D], contains the following remarkable results: , respectively. Following [Be3], the estimates at Points (b) and (d) are due to an 'obscure' cancellation phenomenon, according to which the natural leading term in variance (that should be of the order of √ E and E 3/2 , respectively) cancels out in the high-energy limit. The content of Point (b) has been rigorously confirmed by Wigman [W] in the related model of real random spherical harmonics, whose scaling limit is again the real RWM. See also [A-L-W].
As explained in [K-K-W], albeit improving conjectures from [R-W], the order of the variance established in (1.18) differs from that predicted in (b): this discrepancy is likely due to the fact that, differently from random spherical harmonics, the convergence in (1.22) does not take place uniformly over suitable regions. As already discussed, in [M-P-R-W] it was shown that the asymptotic relation (1.18) is generated by a remarkable chaotic cancellation phenomenon, which also explains the non-central limit theorem stated in (1.20).
The main result of the present paper (see Theorem 1.2) confirms that such a chaotic cancellation continues to hold for phase singularities of complex arithmetic waves, and that it generates non-universal and non-central second order fluctuations for such a random quantity. This fact lends further evidence to the natural conjecture that cancellation phenomena analogous to those described in [Be3,W,Ro] should hold for global quantities associated with the zero set of Laplace eigenfunctions on general mani-folds, as long as such quantities can be expressed in terms of some area/co-area integral formula.
We stress that the fact that the order of the variance stated in (1.8) differs from the one predicted at Point (d) above, can once again be explained by the non-uniform nature of the scaling relation (1.22).
Leray measures and occupation densities. While the present paper can be seen as a natural continuation of the analysis developed in [K-K-W, M-P-R-W], the methods implemented below will substantially diverge from those in the existing literature. One fundamental difference stems from the following point: in order to deal with strong correlations between vectors of the type (T n (x), ∂/∂ 1 T n (x), ∂/∂ 2 T n (x)) and (T n (y), ∂/∂ 1 T n (y), ∂/∂ 2 T n (y)), x = y, the authors of [K-K-W] extensively use results from [O-R-W] (see in particular Section 4.1]) about the fluctuations of the Leray measure which is defined as the limit in L 2 (P) of the sequence k → T ϕ k (T n (x)) dx, with {ϕ k } a suitable approximation of the identity; on the other hand, following such a route in the framework of random phase singularities is impossible, since the formal quantity cannot be defined as an element of L 2 (P). The technical analysis of singular points developed in Section 6 is indeed how we manage to circumvent this difficulty. We observe that, in the parlance of stochastic calculus, the quantity A n (resp. B n ) is the occupation density at zero of the random field T n (resp. T n ) -in particular, the fact that A n is well-defined in L 2 (P) and B n is not -follows from the classical criterion stated in [G-H, Theorem 22.1], as well as from the relations where we have used the fact that, according e.g. to [O-R-W, Lemma 5.3], the mapping x → (1 − r 2 n (x)) −1 behaves like a multiple of 1/ x − x 0 2 around any point x 0 such that r n (x 0 ) = ±1.
Nodal intersections of arithmetic random waves with a fixed curve. A natural problem related to the subject of our paper is that of studying the number of nodal intersections with a fixed deterministic curve C ⊂ T whose length equals L, i.e. number of zeroes of T n that lie on C: Z n := T −1 n (0) ∩ C. In [R-W2], the case where C is a smooth curve with nowhere zero-curvature has been investigated. The expected number of nodal intersections is E[|Z n |] = (π √ 2) −1 × E n × L, hence proportional to the length L of the curve times the wave number, independent of the geometry. The asymptotic behaviour of the nodal intersections variance in the high energy limit is a subtler matter: it depends on both the angular distribution of lattice points lying on the circle with radius corresponding to the given wavenumber, in particular on the sequence of measures {µ n }, and on the geometry of C. The asymptotic distribution of |Z n | isanalyzed in [Ro-W]. See [Ma] for the case where C is a segment.
Zeros of random analytic functions/Systems of polynomials. To the best of our expertise, our limit result (1.11) is the first non-central limit theorem for the number of zeros of random complex analytic functions defined on some manifold M. As such, our findings should be contrasted with the works by Sodin and Tsirelson [S-T, N-S], where one can find central limit results for local statistics of zeros of analytic functions corresponding to three different models (elliptic, flat and hyperbolic). As argued in [W,Section 1.6.4], these results are roughly comparable to those one would obtain by studying zeros of complex random spherical harmonics, for which a central high-energy behaviour is therefore likely to be expected. References [S-Z1, S-Z2], by Shiffman and Zelditch, contain central limit result for the volume of the intersection of the zero sets of independent Gaussian sections of high powers of holomorphic line bundles on a Kähler manifold of a fixed dimension.
In view of Proposition 1.6, our results have to be compared with works dealing with the number of roots of random system of polynomials. The first important result about the number of roots of random systems is due to Shub and Smale [S-S], where the authors compute the expectation of the number of roots of a square system with independent centered Gaussian coefficients with a particular choice of the variances that makes the distribution of the polynomials invariant under the action of the orthogonal group of the parameter space. This model is called the Shub-Smale model. Later, Edelman and Kostlan, see [K] and references therein, and Azaïs and Wschebor [A-W1] extend these results to more general Gaussian distributions. Wschebor [Ws1] studies the asymptotic for the variance of the number of roots of a Shub-Smale system in the case where the number of equations and variables tends to infinity. Armentano and Wschebor [Ar-W] consider the expectation of non-centered (perturbed) systems. McLennan [McL] studies the expected number of roots of multihomogeneous systems. Rojas [Ro] consider the expected number of roots of certain sparse systems.

Short plan of the paper
Section 2 gathers some preliminary results that will be needed in the sequel. In Section 3 we explain the main ideas and steps of the proof of our main result (Theorem 1.2). Finally, the remaining sections are devoted to the detailed proofs. In particular, we collect in Section 8 some technical computations and proofs of intermediate results.

Wiener chaos
Let {H k : k = 0, 1, ...} be the sequence of Hermite polynomials on R, recursively defined as follows: H 0 ≡ 1, and, for k ≥ 1, It is a standard fact that the collection dt is the standard Gaussian measure on the real line. By construction, for every k ≥ 0, one has that (2.24) In view of Proposition 1.4 (recall also Assumption 1.5), every random object considered in the present paper is a measurable functional of the family of complex-valued Gaussian random variables where A(n) and A(n) are defined in (1.15). Now define the space A to be the closure in L 2 (P) of all real finite linear combinations of random variables ξ of the form where λ, τ ∈ Z 2 , z, u ∈ C and c 1 , c 2 ∈ R. The space A is a real centered Gaussian Hilbert subspace of L 2 (P).
Definition 2.1. For a given integer q ≥ 0, the q-th Wiener chaos associated with A, denoted by C q , is the closure in L 2 (P) of all real finite linear combinations of random variables of the type k j=1 H p j (ξ j ), with k ≥ 1, where the integers p 1 , ..., p k ≥ 0 verify p 1 + · · · + p k = q, and (ξ 1 , ..., ξ k ) is a centered standard real Gaussian vector contained in A (so that C 0 = R).
In view of the orthonormality and completeness of H in L 2 (γ), it is not difficult to show that C q ⊥ C q ′ (where the orthogonality holds in L 2 (P)) for every q = q ′ , and moreover the previous relation simply indicates that every real-valued functional F of A can be uniquely represented in the form where F [q] := proj(F | C q ) stands for the the projection of F onto C q , and the series converges in L 2 (P). By definition, one has

About gradients
Differentiating both terms in (1.16) yields that, for j = 1, 2, (where we used the shorthand notation ∂ j = ∂ ∂x j ). It follows that, for every n ∈ S and every x ∈ T, Another important fact (that one can check by a direct computation) is that, for fixed x ∈ T, the six random variables appearing in (2.27) are stochastically independent.
Routine computations (see also [R-W, Lemma 2.3]) yield that for any j = 1, 2, any n and any x ∈ T. Accordingly, we will denote by ∂ j the normalised derivative and adopt the following (standard) notation for the gradient and its normalised version:

Leonov-Shiryaev formulae
In the proof of our variance estimates, we will crucially use the following special case of the so-called Leonov-Shiryaev combinatorial formulae for computing joint moments. It follows immediately e.g. from [P-T, Proposition 3.2.1], by taking into account the independence of T n and T n , the independence of the six random variables appearing in (2.27), as well as the specific covariance structure of Hermite polynomials (see e.g. [N-P, Proposition 2.2.1]).

Arithmetic facts
We will now present three results having an arithmetic flavour, that will be extensively used in the proofs of our main findings. To this end, for every n ∈ S we introduce the 4-and 6-correlation set of frequencies The first statement exploited in our proofs yields an exact value for |S 4 (n)|; the proof is based on an elegant geometric argument due to Zygmund [Zy], and is included for the sake of completeness.
Lemma 2.3. For every n ∈ S, every element of S 4 (n) has necessarily one of the following four (exclusive) structures: In particular, |S 4 (n)| = 3N n (N n − 1).
We will also need the following elementary fact about the behaviour of the correlation function r n , as defined in (1.14), in a small square containing the origin.
Finally, we will make use of the following result, corresponding to a special case of [Ko,Corollary 9,p. 80]: it yields a local ersatz of Bézout theorem for systems of equations involving trigonometric polynomials. We recall that, given a smooth mapping U : R 2 → R 2 and a point x ∈ R 2 such that U (x) = (0, 0), one says that x is nondegenerate if the Jacobian matrix of U at x is invertible.
The next section contains a precise description of the strategy we will adopt in order to attack the proof of Theorem 1.2 3 Structure of the proof of Theorem 1.2

Chaotic projections and cancellation phenomena
We will start by showing in Lemma 4.1 that I n can be formally obtained in L 2 (P) as where δ 0 denotes the Dirac mass in 0 = (0, 0), J Tn is the Jacobian matrix and |J Tn | is shorthand for the absolute value of its determinant. Since I n is a squareintegrable functional of a Gaussian field, according to the general decomposition (2.25) one has that where I n [q] = proj(I n |C q ) denotes the orthogonal projection of I n onto the q-th Wiener chaos C q . Since I n [0] = E[I n ], the computation of the 0-order chaos projection will allow us to conclude the proof of Part 1 of Theorem 1.2 in Section 4.2. One crucial point in our analysis is that, as proved in Lemma 4.4, the projections of I n onto odd-order Wiener chaoses vanish and, more subtly, also the second chaotic component disappears. Namely, we will show that, for every n ∈ S, it holds Our proof of (3.31) is based on Green's identity and the properties of Laplacian eigenfunctions (see also [Ro,Section 7.3 and p.134]).

Leading term: fourth chaotic projections
The first non-trivial chaotic projection of I n to investigate is therefore I n [4]. One of the main achievements of our paper is an explicit computation of its asymptotic variance, as well as a proof that it gives the dominant term in the asymptotic behaviour of the total variance Var(I n ) = q≥2 Var(I n [2q]). The forthcoming Propositions 3.1, 3.2 and 3.3, that we will prove in Section 5, are the key steps in order to achieve our goals.
In view of Remark 1.3(2), Proposition 3.1 coincides with Part 2 of Theorem 1.2, once we replace I n j [4] with I n j . Let us now set, for n ∈ S, where S n (4), S n (6) are the sets of 4-and 6-correlation coefficients defined in Section 2.4, and we have used Lemma 2.3 in (3.32). The following result (Proposition 3.2), combined with Proposition 3.1 and Lemma 2.4 allows us to conclude that, as N n → ∞, where J η is defined in (1.11).

Controlling the variance of higher-order chaoses
In order to prove Proposition 3.2, we need to carefully control the remainder given by We partition the torus into a union of disjoint squares Q of side length 1/M , where M is proportional to √ E n . Of course where proj (·|C ≥6 ) denotes the orthogonal projection onto q≥6 C q , that is, the orthogonal sum of Wiener chaoses of order larger or equal than six. We now split the double sum on the RHS of (3.36) into two parts: namely, one over singular pairs of cubes and the other one over non-singular pairs of cubes. Loosely speaking, for a pair of non-singular cubes (Q, Q ′ ), we have that for every (z, w) ∈ Q × Q ′ , the covariance function r n of the field T n and all its normalized derivatives up to the order two ∂ i r n , ∂ ij r n := (E n /2) −1 ∂/∂ x i x j r n for i, j = 1, 2 are bounded away from 1 and −1, once evaluated in z − w (see Definition 6.1 and Lemma 6.2).
Lemma 3.4 (Contribution of the singular part). As N n → +∞, In order to show Lemma 3.4 (see Section 6), we use the Cauchy-Schwarz inequality and the stationarity of T n , in order to reduce the problem to the investigation of nodal intersections in a small square Q 0 around the origin: for the LHS of (3.37) we have (Q,Q ′ ) sing.

Cov proj
Thus, we need to (i) count the number of singular pairs of cubes, (ii) compute the expected number of nodal intersections in Q 0 and finally (iii) calculate the second factorial moment of I n | Q 0 . Issue (i) will be dealt with by exploiting the definition of singular pairs of cubes and the behavior of the moments of the derivatives of r n on the torus (see Lemma 6.3), thus obtaining that |{(Q, Q ′ ) sing.}| ≪ E 2 n R n (6).
Relations (1.12) and (3.35) yield immediately that E I n | Q 0 is bounded by a constant independent of n.
To deal with (iii) is much subtler matter. Indeed, we need first to check the assumptions for Kac-Rice formula (see [A-W2, Theorem 6.3] ) to hold in Proposition 2.5. The latter allows us to write the second factorial moment E I n | Q 0 I n | Q 0 − 1 as an integral on Q 0 × Q 0 of the so-called two-point correlation function K 2 , given by K 2 (x, y) := p (Tn(x),Tn(y)) (0, 0)E |J Tn (x)| |J Tn (y)| T n (x) = T n (y) = 0 , where x, y ∈ T and p (Tn(x),Tn(y)) is the density of (T n (x), T n (y)).
The stationarity of T n then reduces the problem to investigating K 2 (x) := K 2 (x, 0) around the origin. Hypercontractivity properties of Wiener chaoses and formulas for the volume of ellipsoids (see [K-Z]) yield the following estimation where |Ω n (x)| stands for the absolute value of the determinant of the matrix Ω n (x), defined as the covariance matrix of the vector ∇T n (0), conditionally on T n (x) = T n (0) = 0. An explicit Taylor expansion at 0 for Ψ n (made particularly arduous by the diverging integral in (1.23) -see Lemma 8.1) will allow us to prove that E I n | Q 0 I n | Q 0 − 1 is also bounded by a constant independent of n. This concludes the proof of Lemma 3.4.
To achieve the proof of Theorem 1.2, we will eventually show the following result, whose proof relies on Proposition 2.2, on the definition of non-singular cubes, as well as on the behavior of even moments of derivatives of the covariance function r n .
Lemma 3.5 (Contribution of the non-singular part). As N n → +∞, we have The rest of the paper contains the formal proofs of all the statements discussed in the present section.

Phase singularities and Wiener chaos 4.1 Chaotic expansions for phase singularities
In this part we find the chaotic expansion (3.30) for I n . The first achievement in this direction is the following approximation result.

An integral expression for the number of zeros
For ε > 0 and n ∈ S, we consider the ε-approximating random variable where 1 [−ε,ε] 2 denotes the indicator function of the square [−ε, ε] 2 . The following result makes the formal equality in (3.29) rigorous.
Lemma 4.1. For n ∈ S, with probability one, I n is composed of a finite number of isolated points and, as ε → 0, I n (ε) → I n , (4.40) both a.s. and in the L p (P)-sense, for every p ≥ 1.
Proof. Fix n ∈ S. In order to directly apply some statements taken from [A-W2], we will canonically identify the random field (x 1 , x 2 ) → T n (x 1 , x 2 ) with a random mapping from R 2 to R 2 that is 1-periodic in each of the coordinates x 1 , x 2 . In what follows, for x ∈ R 2 we will write T n (x, ω) to emphasize the dependence of T n (x) on ω ∈ Ω. We subdivide the proof into several steps, numbered from (i) to (vi).
(i) First of all, since T n is an infinitely differentiable stationary Gaussian field such that, for every x ∈ R 2 , the vector T n (x) has a standard Gaussian distribution, one can directly apply [A-W2, Proposition 6.5] to infer that there exists a measurable set Ω 0 ⊂ Ω with the following properties: P(Ω 0 ) = 1 and, for every ω ∈ Ω 0 and every x ∈ R 2 such that T n (x, ω) = 0, one has necessarily that the Jacobian matrix J Tn (x, ω) is invertible.
(ii) A standard application of the inverse function theorem (see e.g. [A-T, p. 136]) implies that, for every ω ∈ Ω 0 , any bounded set B ⊂ R 2 only contains a finite number of points x such that T n (x, ω) = 0. This implies in particular that, with probability one, I n (as defined in (1.3)) is composed of a finite number of isolated points and I n < +∞.
(iii) Sard's Lemma yields that, for every ω ∈ Ω 0 , there exists a set U ω ⊂ R 2 such that U c ω has Lebesgue measure 0 and, for every u ∈ U ω there is no x ∈ R 2 such that T n (x, ω) = u and J Tn (x, ω) is not invertible. Note that, by definition, one has that 0 ∈ U ω for every ω ∈ Ω 0 .
(iv) Define B := {x = (x 1 , x 2 ) ∈ R 2 : 0 ≤ x i < 1/L}, i = 1, 2}, where L is any positive integer such that L > √ n. For every u ∈ R 2 , we set I n,u (B) to be the cardinality of the set composed of those x ∈ B such that T n (x) = u; the quantity I n,u (T) is similarly defined, in such a way that I n,0 (T) = I n . Two facts will be crucial in order to conclude the proof: (a) for every ω ∈ Ω 0 and every u = (u 1 , u 2 ) ∈ U ω , by virtue of Lemma 2.6 as applied to the pair (P, Q) given by as well as of the fact that B ⊂ W , one has that I n,u (B)(ω) ≤ α(n), and (b) as a consequence of the inverse function theorem, for every ω ∈ Ω 0 there exists η ω ∈ (0, ∞) such that the equality I n (ω) = I n,u (T)(ω) holds for every u such that u ≤ η ω . Indeed, reasoning as in [A-T, Proof of Theorem 11.2.3] if this was not the case, then there would exist a sequence u k → 0, u k = 0, and a point x ∈ T such that: (1) T n (x, ω) = 0, and (2) for every neighborhood V of x (in the topology of T) there exist k ≥ 1 and x 0 , x 1 ∈ V such that x 0 = x 1 and T n (x 0 ) = T n (x 1 ) = u k -which is in contradiction with the inverse function theorem.
(v) By the area formula (see e.g. [A-W2, Proposition 6.1 and formula (6.2)]), one has that, for every ω ∈ Ω 0 , where we used the property that the complement of U ω has Lebesgue measure 0. Since the integral on the right-hand side of (4.41) equals I n whenever ε ≤ η ω / √ 2, we conclude that (4.40) holds P-a.s. .
The fact that (4.40) holds also in L p (P) now follows from Point (v) and dominated convergence.

Chaotic expansions
Let us consider the collections of coefficients {β l : l ≥ 0} and {α a,b,c,d : a, b, c, d ≥ 0} defined as follows. For l ≥ 0 where (as before) H 2l is the 2l-th Hermite polynomial. For instance, with (X, Y, V, W ) a standard real four-dimensional Gaussian vector. Note that on the right-hand side of (4.44), |XW − Y V | is indeed the absolute value of the determinant of the matrix X Y V W . Proof. Let us assume without loss of generality (by symmetry) that a is odd and that at least one integer among b, c and d is even. We will exploit (2.24). If b is even, then, using leading to α a,b,c,d = 0. If c (resp. d) is even, the same reasoning based on (X, Y, V, W ) ) leads to the same conclusion.
Proof. The main idea is to deduce the chaotic expansion for I n from the chaotic expansion for (4.39) and Lemma 4.1. Let us first rewrite (4.39) as where, for l ≥ 0 (4.50) and φ is still denoting the standard Gaussian density. For the indicator function of [−ε, ε] 2 appearing in (4.49), we thus have The chaotic expansion for the absolute value of the Jacobian determinant appearing in (4.49) is, thanks to Lemma 4.2, where α a,b,c,d are given in (4.44). In particular, observe that Lemma 4.2 ensures that the odd chaoses vanish in the chaotic expansion for the Jacobian. It hence follows from (4.51) and (4.52) that the chaotic expansion for I n (ε) in (4.49) is (taking sums over even i 1 , j 1 and i 2 , i 3 , j 2 , j 3 with the same parity) I n (ε) = E n 2 q≥0 i 1 +i 2 +i 3 +j 1 +j 2 +j 3 =2q (4.53) Noting that, as ε → 0, β ε l → β l , where β l are given in (4.42) and using Lemma 4.1, we prove both (4.45) and (4.46).

Proof of Part 1 of Theorem 1.2
According to Lemma 4.3 and Equation (4.43), for every n ∈ S one has that thus yielding the desired conclusion.

Investigation of the fourth chaotic components
In this section we shall investigate fourth chaotic components. In particular, we shall prove Proposition 3.1 and Proposition 3.3.
Putting everything together, we arrive at the following expression for B n . Now, let us prove that each component of W n j is asymptotically Gaussian as N n j → +∞. Since all components of W n j belong to the same Wiener chaos (the second one) and have a converging variance (see indeed the diagonal part of B n just above), according to the Fourth Moment Theorem (see, e.g., Theorem 5.2.7]) it suffices to show that the fourth cumulant of each component of W n j goes to zero as N n j → +∞. Since we are dealing with sum of independent random variables, checking such a property is straightforward. For sake of illustration, let us only consider the case of W 2 (n j ) which is representative of the difficulty. We recall that, given a real-valued random variable Z with mean zero, the fourth cumulant of Z is defined by κ 4 (Z) : Since the a λ are independent except for the relation a λ = a −λ , we can write, setting Λ + n = {λ ∈ Λ n : λ 2 > 0}, to obtain the last inequality, we have used that λ 2 2 ≤ λ 2 1 +λ 2 2 = n. As a result, κ 4 (W 2 (n j )) → 0 as N n j → +∞ and it follows from the Fourth Moment Theorem that W 2 (n j ) is asymptotically Gaussian. It is not difficult to apply a similar strategy in order to prove that, actually, each component of W n j is asymptotically Gaussian as well; details are left to the reader.
Finally, we make use of [N-P, Theorem 6.2.3] to conclude the proof of Lemma 5.3. Indeed, (i) all components of W n belong to the same Wiener chaos (the second one), (ii) each component of W n j is asymptotically Gaussian (as N n j → +∞), and finally (iii) Σ k,l (n j ) → M k,l (η) for each pair of indices (k, l).
is a sequence of random variables belonging to a fixed Wiener chaos and converging in distribution, by standard arguments based on uniform integrability, we also have the proof of Proposition 3.1 is then concluded, once computing 10 − 4G 2 11 + 4G 2 9 − 2G 2 12 − 2G 2 13 + 8G 2 14 − 12G 12 G 13 = 8(3η 2 + 5), and noting that the latter variance is the same in both cases (i) and (ii).

The variance of higher order chaoses
In this section we shall prove Proposition 3.2. Let us decompose the torus T as a disjoint union of squares Q k of side length 1/M (where M ≈ √ E n is a large integer 2 ), obtained by translating along directions k/M , k ∈ Z 2 , the square Q 0 := [0, 1/M ) × [0, 1/M ) containing the origin. By construction, the south-west corner of each square is therefore situated at the point k/M .

Singular points and cubes
Let us first give some definitions, inspired by [O-R-W, §6.1] and [R-W2, §4.3]. Let us denote by 0 < ε 1 < 1 10 10 a very small number 3 that will be fixed until the end. From now on, we shall use the simpler notation r j := ∂ j r n , and r ij := ∂ ij r n for i, j = 1, 2.
ii) A pair of cubes (Q, Q ′ ) is called singular if the product Q × Q ′ contains a singular pair of points.
For instance, (0, 0) is a singular pair of points and hence (Q 0 , Q 0 ) is a singular pair of cubes. In what follows we will often drop the dependence of k from Q k .
Proof. First note that the function T ∋ s → r(s/ √ n) and its derivatives up to the order two are Lipschitz with a universal Lipschitz constant c = 8π 3 (in particular, independent of n). Let us denote by (x, y) the singular pair of points contained in Q × Q ′ and suppose that r(x − y) > ε 1 . For every (z, w) ∈ Q × Q ′ , The case r(x − y) < −ε 1 in indeed analogous. The rest of the proof for derivatives follows the same argument.
Let us now denote by B Q the union of all squares Q ′ such that (Q, Q ′ ) is a singular pair. The number of such cubes Q ′ is M 2 Leb(B Q ), the area of each cube being 1/M 2 . Lemma 6.3. It holds that Leb(B Q ) ≪ T r(x) 6 dx.
Proof. Let us first note that where B 0 Q is the union of all cubes Q ′ such that there exists (x, y) ∈ Q × Q ′ enjoying |r(x − y)| > ε 1 and for i, j = 1, 2, B i Q is the union of all cubes Q ′ such that there exists (x, y) ∈ Q × Q ′ enjoying |r i (x − y)| > ε 1 √ n and finally B ij Q is the union of all cubes Q ′ such that there exists (x, y) ∈ Q × Q ′ enjoying |r ij (x − y)| > ε 1 n. We can hence write Let us now fix z ∈ Q; then Lemma 6.2 yields Moreover, for i = 1, 2, where r i := r i / √ E n are the normalized derivatives. Since An analogous argument applied to B ij Q for i, j = 1, 2 allows to conclude the proof.
The number of cubes Q ′ such that the pair (Q, Q ′ ) is singular is hence negligible with respect to E n R n (6).

Variance and cubes
We write the total number I n of nodal intersections as the sum of the number I n | Q of nodal intersections restricted to each square Q, i.e.
We are going to separately investigate the contribution of the singular pairs and the nonsingular pairs of cubes: Var(proj(I n |C ≥6 )) = (Q,Q ′ ) sing.

The contribution of singular pairs of cubes
Proof of Lemma 3.4. By Cauchy-Schwarz inequality and the stationarity of T n , recalling moreover Lemma 6.3, we have where, from now on, Q 0 denotes the square containing the origin. Now, It is immediate to check that in particular E I n | Q 0 = O(1). Note that A is the 2-th factorial moment of I n | Q 0 : Applying [A-W2, Theorem 6.3], we can write where K 2 (x, y) := p (Tn(x),Tn(y)) (0, 0) E |J Tn (x)| · |J Tn (y)| T n (x) = T n (y) = 0 is the so-called 2-point correlation function. Indeed, Proposition 2.5 ensures that for (x, y) ∈ Q 0 × Q 0 , the vector (T n (x), T n (y)) is non-degenerate except on the diagonal x = y. Note that, by stationarity of the model, we can write (6.59) as where K 2 (x) := K 2 (x, 0) and Q 0 is 2Q 0 . Let us first check that the function x → K 2 (x) is integrable around the origin. Note that, by Cauchy-Schwarz inequality, Hypercontractivity on Wiener chaoses [N-P] ensures that there exists c > 0 such that where Ω n (x) denotes the covariance matrix of ∇T n (0) conditioned to T n (x) = T n (0) = 0 (see (8.66) for a precise expression). The Taylor expansion in Lemma 8.1 gives that, as x → 0, for some other constant c > 0, where the constants involving in the 'O' notation do not depend on n, so that which is the result we looked for.
However, we need the speed of convergence to 0 of det Ω. Hence, we will use a Taylor expansion argument around 0. We denote x = (x, y). Lemma 8.1. As x → 0, we have det Ω n (x) = cE 3 n x 2 + E 4 n O( x 4 ), and hence where both c > 0 and the constants in the 'O' notation do not depend on n.
Proof. Let us start with r.
More precisely, the remainder R r n is of the form R r n = O sup ∂ 6 r n x 6 .
Also here, the remainders R 1 n and R 2 n are both of the form