Laurent Bétermin

Ladislav Šamaj

Igor Travěnec

Interplay between critical and off-critical zeros of two-dimensional Epstein zeta functions

Keywords: Epstein zeta function on the rectangular lattice, analytic continuation, zeros off the critical line, Jacobi theta functions AMS CLASSIFICATION 11E45, 11M41, 11R42

The two-dimensional Epstein zeta function associated to a rectangular lattice with spacings ax = 1 and ay = ∆ ∈ R * + , defined by

, where the sum goes over all integers except of the origin (j, k) = (0, 0), is studied. It can be analytically continued to the whole complex s-plane except for the point s = 1. The nontrivial zeros {ρ = ρx + iρy} of the Epstein zeta function, defined by ζ (2) (ρ, ∆) = 0, split into "critical" zeros (on the critical line ρx = 1 2 ) and "off-critical" zeros (ρx = 1 2 ). This work presents rigorous asymptotic and analytic results as well as numerical investigations. According to the present numerical calculation, the critical zeros form open or closed curves ρy(∆) in the plane (∆, ρy). Nearest critical zeros merge at special points (∆ * , ρ * y ), referred to as left/right edge zeros, which are defined by a divergent tangent dρy/d∆. Each of these critical edge zeros gives rise to a continuous curve of off-critical zeros which can thus be generated systematically. As a rule, each curve of off-critical zeros joins a pair of left and right edge zeros. It is shown that in the regions of small/large values of the anisotropy parameter ∆ the Epstein zeta function can be approximated adequately by a function which reveals an equidistant distribution of critical zeros along the imaginary axis in the limits ∆ → 0 and ∆ → ∞. It is also numerically found that for each ∆ ∈ (0, ∆

141733 there exists a pair of real off-critical zeros, their ρx components going to the borders 0 and 1 of the critical region in the limits ∆ → 0, ∞.

We consider the two-dimensional (2D) Epstein zeta function [START_REF] Epstein | Zur Theorie allgemeiner Zetafunctionen[END_REF][START_REF] Epstein | Zur Theorie allgemeiner Zetafunctionen II[END_REF] associated to the rectangular lattice a x Z ⊕ a y Z with spacings (a x , a y ) = (1, ∆), ∆ ∈ R * + , defined by

ζ (2) (s, ∆) := 1 2 ′ (j,k)∈Z 2 1 (j 2 + ∆ 2 k 2 ) s ℜ(s) > 1, (1.1)
where ′ means that the term (j, k) = (0, 0) is excluded from the summation. This represents a natural 2D extension of the one-dimensional (1D) Riemann zeta function [START_REF] Riemann | Über die Anzahl der Primzahlen unter einer gegebenen Grösse[END_REF] defined by

ζ(s) = 1 2 ′ j∈Z 1 |j| s = j∈N 1 j s .
(1.

2)

The sum (1.1) naturally arises in Physics as the energy per point of a system of identical particles located on the lattice sites of Z ⊕ ∆Z and interacting pairwisely through the Riesz potential r → 1/r 2s [START_REF] Brauchart | Optimal discrete Riesz energy and discrepancy[END_REF]. Notice that the prefactor 1 2 is present because each interaction energy is shared by a pair of particles.

The function ζ (2) (s, ∆) possesses the obvious symmetry ζ (2) (s, ∆) = 1 ∆ 2s ζ (2) (s, 1/∆) (1.3) which means that the values of ∆ can be constrained to either of the intervals (0, 1] or [1, ∞).

Similarly as the Riemann zeta function, the Epstein zeta function can be analytically continued to the whole complex s-plane [START_REF] Elizalde | Regularization of general multidimensional Epstein Zeta-functions[END_REF][START_REF] Ennola | On a problem about the Epstein zeta function[END_REF][START_REF] Epstein | Zur Theorie allgemeiner Zetafunctionen[END_REF][START_REF] Epstein | Zur Theorie allgemeiner Zetafunctionen II[END_REF]. The region 0 < ℜ(s) < 1 is referred to as the critical strip, the critical line is defined by ℜ(s) = 1 2 . For ∆ 2 ∈ {1, 2, 3, 4, 7}, the 2D lattice sum (1.1) can be expressed as a product of simpler 1D sums, namely Dirichlet L-series, whereas for other special integer values of ∆ it is expressible as a sum of products of Dirichlet L-functions [START_REF] Borwein | Lattice Sums Then and Now[END_REF]. For instance, in the isotropic square lattice case ∆ 2 = 1, it holds ζ (2) (s, 1) = 2ζ(s)β(s), where β(s) :

= ∞ j=0 (-1) j (2j + 1) s = 1 2 2s ζ s, 1 4 -ζ s, 3 4 (1.4)
is the Dirichlet beta function and ζ(s, a) := ∞ j=0 (j + a) -s denotes the Hurwitz zeta function. The zeros {ρ = ρ x +iρ y } of the Epstein zeta function s → ζ (2) (s, ∆) associated to the parameter ∆ are defined by the equality ζ (2) (ρ, ∆) = 0. Besides the trivial zeros at ρ ∈ Z -there exist two kinds of nontrivial zeros:

• the "critical" zeros which lie on the critical line ρ x = 1 2 ; • the "off-critical" zeros which have ρ x = 1 2 . The Riemann hypothesis states that all nontrivial zeros of the 1D Riemann zeta function ζ are constrained to the corresponding critical line ρ x = 1 2 [START_REF] Riemann | Über die Anzahl der Primzahlen unter einer gegebenen Grösse[END_REF]. Provided that the Riemann hypothesis holds, all nontrivial zeros of the Dirichlet beta function (1.4) are constrained to the critical line ρ x = 1 2 as well [START_REF] Lander | The Zeros of the Dirichlet beta function encode the odd primes and have real part 1/2[END_REF]. Consequently, all nontrivial zeros of the Epstein zeta function ζ (2) (s, 1) associated to the square lattice, given by (1.4), lie on the critical line. Similar phenomenon is expected also for ∆ 2 ∈ {2, 3, 4, 7} when the Epstein zeta function factorizes itself into the product of Dirichlet L-functions; for analytic and numerical studies of the statistics of gaps between critical zeros, see [START_REF] Baier | A note on the gaps between zeros of Epstein's zetafunctions on the critical line[END_REF][START_REF] Bogomolny | Statistical properties of the zeros of zeta functions -beyond the Riemann case[END_REF][START_REF] Bombieri | Sur des zeros des fonctions zêta d'Epstein[END_REF][START_REF] Jutila | Gaps between the zeros of Epstein's zeta-functions on the critical line[END_REF]. This is no longer true for anisotropic (i.e. non-square) rectangular lattices (1.1) with other integer values of ∆ 2 . The first off-critical zero of ζ (2) (s, ∆) was detected for ∆ 2 = 5 [START_REF] Potter | The zeros of Epstein's zeta functions[END_REF], since then many other off-critical zeros were identified [START_REF] Bateman | On Epstein's zeta function[END_REF][START_REF] Davenport | On the zeros of certain Dirichlet series I[END_REF][START_REF] Mcphedran | Zeros of lattice sums: 1. Zeros of the critical line[END_REF][START_REF] Stark | On the zeros of Epstein's zeta function[END_REF].

Main results and related works

In this paper, we perform a numerical investigation of the critical zeros of s → ζ (2) (s, ∆) for ∆ ∈ (0, 1] in a compact set and explain how the off-critical zeros of this Epstein zeta function emerge from certain critical "edge" zeros. Most of our numerical findings are based on rigorous asymptotic and analytic results. Furthermore, we systematically check the degree of accuracy between our approximations and the direct numerical computations from more complicated equations. More precisely, our observation and results are the following:

(1) For ζ (2) (s, ∆) one obtains numerically in the plane (∆, ρ y ) continuous curves of critical zeros ρ y (∆) (remember that ρ x = 1/2) and identifies critical edge zeros on them that correspond to the points of the curve with vertical tangents. These curves and critical edge zeros are depicted in Figure 1 and their construction is explained in Section 3. (2) The critical edge zeros are the starting points to generate systematically (here numerically) continuous curves of off-critical zeros. These continuous curves, connecting a pair of left and right edge zeros, are depicted in Figure 2 and their construction is explained in Section 4.

(3) The analysis of the limits ∆ → 0 and ∆ → ∞ reveals an equidistant distribution of critical zeros along the imaginary axis, with spacing between the nearest zeros going to zero as O(1/| log ∆|). This analysis is done in Section 3.3. (4) A pair of real off-critical zeros is numerically found for each ∆ ∈ (0, ∆

* c ] ∪ [1/∆ * c , ∞) with ∆ * c ≈ 0.
141733. This is explained in Section 4.2. Furthermore, our asymptotic study for small values of ∆ combined with a result by Montgomery [START_REF] Montgomery | Minimal Theta Functions[END_REF] on theta functions (see also [START_REF] Faulhuber | Optimal Gabor frame bounds for separable lattices and estimates for Jacobi theta functions[END_REF]) yields the following conjecture that we have also numerically checked.

Conjecture 1.1. The unique solution of ζ (2) (1/2, ∆) = 0 for ∆ ∈ (0, 1) is ∆ = ∆ * c = e γ 4π
where γ is the Euler-Mascheroni constant.

Our work is motivated by another recent paper on d-dimensional Epstein zeta function where surprising similarities appear. Indeed, the Epstein zeta function can be defined for any lattice structure and the two last authors have recently studied in [START_REF] Travěnec | Generation of off-critical zeros for hypercubic Epstein zeta functions[END_REF] the properties of the Epstein zeta function on an isotropic hypercubic lattice

ζ (d) (s) = 1 2 ′ (n1,n2,...,nd)∈Z d 1 (n 2 1 + n 2 2 + • • • + n 2 d ) s/2 ℜ(s) > d, (1.5)
where d = 1, 2, . . . is the spatial dimension and the Riesz interaction energy between two particles at distance r was chosen as 1/r s to reproduce in d = 1 the Riemann zeta function, ζ (1) (s) = ζ(s).

+ + + + + + + + + + + + + + +++++++ + + # # # # # # # # # # # # # # # # # # # # # # # # ## # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # @ @ @
@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ 0 0 0 0 000000000000000000000 00 0 0 0 0 0.2 0.4 0.6 0. An analytic continuation of the Epstein zeta function ζ (d) (s) to the whole complex s-plane was constructed; the corresponding formula is applicable for the spatial dimension d being a continuous variable ranging from 0 to ∞. Numerical calculations of critical zeros (with ρ x = d 2 ) indicate that they form closed or semi-open curves in the plane ρ x = d 2 , ρ y . Each curve involves a number of left/right "edge" points at which a couple of nearest critical zeros merge. The critical edge zeros give rise to two tails of off-critical zeros, coupled via a symmetry, with continuously varying dimension d and the ρ x -component along each tail. This fact permits one to avoid a "blind" search for off-critical zeros, the off-critical zeros are generated systematically starting from the critical edge zeros. Another benefit of the method is an exact treatment of the limits d → 0 and d → ∞. An exact formula for lim d→0 ζ (d) (s)/d was derived. In the limit d → ∞, an equidistant distribution of critical zeros along the imaginary axis was obtained, with spacing between the nearest zeros going to zero as 2π/ log(d). As a by-product of the formalism, a conjugate pair of real off-critical zeros was found for each dimension d > 9.24555 . . .. The problem of determining zeros for the d-dimensional isotropic Epstein zeta function (1.5) seems at first sight to be unrelated to that for the 2D Epstein zeta function (1.1). In analogy with [START_REF] Travěnec | Generation of off-critical zeros for hypercubic Epstein zeta functions[END_REF], each of the curves of critical zeros in the (∆, ρ y ) plane pictured in Figure 1 contains a number of left/right edge zeros. The mechanism of generation of off-critical zeros from these edge points, depicted in Figure 2 and explained in Section 4.1, is the same as in the previous case of the d-dimensional isotropic Epstein zeta function.

However, there are small differences in the form of curves of off-critical zeros. In the present 2D case, there is only one curve going out of the classical edge zero, at each point of the curve there are two conjugate off-critical zeros with ρ x and 1 -ρ x coordinates and each curve of off-critical zeros connects just one left and one right edge points. In the previous d-dimensional case, there are two tails of off-critical zeros starting at each edge point, coupled via a symmetry, and a tail + + + + + + + + ++ + 0 0 0 0 0 0 00 000000000000000 000 00 0 0 0 0 @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ 0. either interconnects two edge zeros or starts from one edge zero and ends at some unspecial point in d = 0 or d → ∞ dimensions.

Another analogy with the previous paper [START_REF] Travěnec | Generation of off-critical zeros for hypercubic Epstein zeta functions[END_REF] is an equidistant distribution of critical zeros along the imaginary axis in special regions of model's parameters. In the case of the d-dimensional isotropic Epstein function (1.5), the critical zeros are distributed equidistantly in the limit d → ∞ with the spacing 2π/ log(d) between the nearest-neighbor zeros. In the present case of the 2D anisotropic Epstein function (1.1), the critical zeros are distributed equidistantly in the limits ∆ → 0 and ∆ → ∞, with the spacing of order π/| log(∆)| between the nearest-neighbor zeros (see Section 3.3).

The next similarity with [START_REF] Travěnec | Generation of off-critical zeros for hypercubic Epstein zeta functions[END_REF] consists in the appearance of a couple of real off-critical zeros with the component ρ y = 0, discussed in Section 4.2.

We would like to emphasize that the presented mechanism of generation of continuous curves of off-critical zeros from the edge critical ones might be not the only possible one. Our crucial assumption was that the deviation of ρ x from its critical value 1 2 changes continuously when the off-critical zero goes out the edge zero. A discontinuous change of ρ x was excluded from our analysis.

Plan of the paper. Section 2 concerns technicalities, like recalling the analytic continuation of the 2D Epstein zeta function ζ (2) (s, ∆) to the complex s-plane (Section 2.1) and basic equations for determining critical and off-critical zeros (Section 2.2). Section 3 deals with critical zeros of ζ (2) (s, ∆). Based on numerical calculation of open and closed curves of critical zeros in the (∆, ρ y ) plane, the critical edge zeros are introduced in Section 3.1, together with an explicit form of two coupled integral equations determining their position in the (∆, ρ y ) plane. The singular expansion of ρ(y) around the edge points is discussed in Section 3.2. An accurate approximation of ζ (2) (s, ∆) for small (and large) values of ∆, indicating an equidistant distribution of critical zeros along the imaginary axis, is presented in Section 3.3. Section 4 is about off-critical zeros of ζ (2) (s, ∆). The most important result of this work, the generation mechanism of curves of offcritical zeros starting from critical edge zeros, is explained in Section 4.1. Pairs of real off-critical zero are found in Section 4.2.

Technicalities

2.1. Regularization of ζ (2) (s, ∆)
This part is devoted to a short review of the analytic continuation of the Epstein zeta function (1.1), defined when the real part ℜ(s) > 1, to the whole complex s-plane. We have the following theorem.

Theorem 2.1 (Analytic continuation of the Epstein zeta function [START_REF] Blanc | The crystallization conjecture: a review[END_REF][START_REF] Epstein | Zur Theorie allgemeiner Zetafunctionen[END_REF]). Let ∆ > 0 and s ∈ C. Let us define

Z(s, ∆) := ∆ π s Γ(s)ζ (2) (s, ∆). (2.1)
We have

Z(s, ∆) = - 1 2 1 1 -s + 1 s + 1 2 1 0 t s-1 + t -s θ 3 e -πt∆ θ 3 e -πt/∆ - 1 t dt, (2.2)
where the Jacobi elliptic theta function (with zero argument, see [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF]) is given by

θ 3 (q) = j∈Z q j 2 . (2.3)
In particular, (2.2) represents an analytic continuation of ζ (2) (s, ∆) to the whole complex plane, except for the simple pole at s = 1. Z(s, ∆) fulfills the following duality relation

Z(s, ∆) = Z(1 -s, ∆). (2.4)
Remark 2.2. Note that the symmetry

Z(s, ∆) = Z(s, 1/∆) (2.5)
is automatically ensured by formula (2.2).

Integral equations determining zeros

We start by defining what is a zero of s → ζ (2) (s, ∆) associated to the parameter ∆.

Definition 2.1. We say that ρ ∈ C is:

• a critical zero (of s → ζ (2) (s, ∆)) associated to the parameter ∆ if ρ = 1 2 + iρ y , ρ y ∈ R and ζ (2) (ρ, ∆) = 0; • an off-critical zero (of s → ζ (2) (s, ∆)) associated to the parameter ∆ if ρ = ρ x + iρ y , (ρ x , ρ y ) ∈ R × R, ρ x = 1 2 and ζ (2) (ρ, ∆) = 0.
The nontrivial zeros of ζ (2) (s, ∆) are related to the nullity condition of the right-hand side of (2.2). In the case of critical zeros ρ = 1 2 + iρ y , the right-hand side of (2.2) is real and we obtain the following result.

Theorem 2.3 (Equation for critical zeros). Let ρ = 1 2 + iρ y be a critical zero associated to the parameter ∆, then ρ y satisfies the following equation:

- 2 1 + 4ρ 2 y + 1 0 cos (ρ y log t) θ 3 e -πt∆ θ 3 e -πt/∆ - 1 t dt √ t = 0. (2.6)
Remark 2.4 (Symmetries of the equation and consequences). The symmetry of (2.6) with respect to the transformation ∆ → 1/∆ tells us that the set of critical zeros is the same for the couple of values ∆ and 1/∆. The symmetry ρ y → -ρ y means that to each critical zero 1 2 , ρ y there exists a complex conjugate critical zero 1 2 , -ρ y .

In the case of off-critical zeros with ρ x = 1 2 , it is useful to introduce the deviation of ρ x from its critical value

δρ x := ρ x - 1 2 . (2.7)
The right-hand side of (2.2) becomes complex in this case and the off-critical zeros are determined by a pair of coupled integral equations (corresponding to the real and imaginary parts).

Theorem 2.5 (Equation for the off-critical zeros). Let ρ = 1 2 + δρ x + iρ y , δρ x = 0, be an off-critical zero associated to the parameter ∆, then ρ satisfies the following equations:

- 1 + 2δρ x (1 + 2δρ x ) 2 + 4ρ 2 y + 1 -2δρ x (1 -2δρ x ) 2 + 4ρ 2 y + 1 0 cos (ρ y log t) cosh (δρ x log t) θ 3 e -πt∆ θ 3 e -πt/∆ - 1 t dt √ t = 0, (2.8)
2ρ y 1 (1 + 2δρ x ) 2 + 4ρ 2 y - 1 (1 -2δρ x ) 2 + 4ρ 2 y + 1 0 sin (ρ y log t) sinh (δρ x log t) θ 3 e -πt∆ θ 3 e -πt/∆ - 1 t dt √ t = 0. (2.9)
Furthermore, these equations are invariant with respect to the transformations ∆ → 1/∆, ρ y → -ρ y and δρ x → -δρ x . In particular, to each off-critical zero 1 2 + δρ x , ρ y there exist the related off-critical zeros 1 2 + δρ x , -ρ y , 1 2 -δρ x , ρ y and 1 2 -δρ x , -ρ y . Each curve of critical zeros in Figure 1 is defined by the function ρ y (∆). Varying the parameter ∆ in the positive or negative direction, the distance between a couple of nearest zeros can go to zero and the zeros merge at points referred to as the critical "edge" zeros. They originate at specific values of ∆ = ∆ * and have imaginary part ρ * y = ρ y (∆ * ). More precisely, the edge zeros are defined by a divergent tangent as follows.

Definition 3.1 (Critical edge zeros). We call (critical) edge zero any critical zero ρ = 1 2 + iρ y such that there exists ∆ * > 0 satisfying

d d∆ ρ y (∆) ∆=∆ * = ±∞. (3.1)
Furthermore, assuming that the curve of critical zeros is defined by the inverse function ∆(ρ y ), the condition for the edge zeros (3.1) is equivalent to

d dρ y ∆(ρ y ) ρy=ρ * y = 0 (3.2)
and ∆ * = ∆(ρ * y ).

The edge zeros split into two groups: the left/right edge zeros are situated on the left/right with respect to the curve ρ y (∆). More precisely: Definition 3.2 (Left/right critical edge zeros). Let ρ * = 1 2 + iρ * y be a (critical) edge zero associated to the parameter ∆ * and the curve of critical zeros is defined by the inverse function ∆(ρ y ). We say that ρ is a right (resp. left) edge zero if

d 2 dρ 2 y ∆(ρ y )| ρy=ρ * y < 0 resp. d 2 dρ 2 y ∆(ρ y )| ρy =ρ * y > 0 .
The right and left edge zeros are denoted in Figure 1 as 2a, 3a, 4a,. . . and 1, 2b, 3b, 4b,. . ., respectively; the close connection between the right and left edge zeros 2a-2b, 3a-3b, etc. will become clear later (see Section 4.1). Notice that the critical zero numbered 1 with coordinates (∆ * , ρ * y ) = (∆ * c ≈ 0.141733, 0) is a left edge zero because the corresponding curve continues reflection-symmetrically across the ∆-axis into the lower quadrant. We can easily deduce closed-form equations for specifying edge zeros.

Theorem 3.1 (Equation satisfied by a critical edge zero). If ρ * = 1 2 + iρ * y is an edge zero associated to the parameter ∆ * , then

f (ρ * y , ∆ * ) = 0, (3.3)
where the function f is given by

f (ρ y , ∆) := 16ρ y (1 + 4ρ y 2 ) 2 - 1 0 log t sin (ρ y log t) θ 3 e -πt∆ θ 3 e -πt/∆ - 1 t dt √ t . (3.4)
Proof. We first recall that an edge zero has to satisfy the general equation (2.6) for critical zeros, i.e.,

- 2 1 + 4ρ * y 2 + 1 0 cos ρ * y log t θ 3 e -πt∆ * θ 3 e -πt/∆ * - 1 t dt √ t = 0. (3.5)
Taking into account the edge-zero condition (3.2), the derivative of (2.6) with respect to ρ y leads to the result.

It has to be noticed that the coupled equations (3.3) and (3.5) have an infinite number of real solutions for ∆ * and ρ * y . The characteristics of edge critical zeros from Figure 1, constrained to the intervals 0 < ∆ ≤ 1 and 0 ≤ ρ y ≤ 21, are summarized in Table 1.

Singular expansion around critical edge zeros

The map ∆ → ρ y (∆) exhibits locally an analytic expansion in ∆, except for the edge zeros where it shows a singular expansion in ∆. The singular expansion around edge zeros can be documented by performing the Taylor series expansion of the inverse function ∆(ρ y ) on the corresponding curve of critical zeros. We obtain the following result. Lemma 3.2 (First order asymptotics around a critical edge zero). Let ρ * = 1 2 + iρ * y be an edge zero associates to the parameter ∆ * . Then, as ρ y → ρ * y where ρ = 1 2 + iρ y is a critical zero associated to the parameter ∆, we have

ρ y -ρ * y = O |∆ -∆ * | .
Proof. The result directly follows from the order two Taylor expansion of the inverse function ∆(ρ y ) given by

∆(ρ y ) = ∆(ρ * y ) + 1 2! d 2 dρ 2 y ∆(ρ y ) ρy=ρ * y ρ y -ρ * y 2 + o ρ y -ρ * y 2 , (3.6)
where the condition (3.2) was taken into account. We now perform a general analysis for critical zeros, deriving an asymptotic expansion of Equation (2.6) around a critical zero ρ = 1 2 + iρ y associated to the parameter ∆.

Lemma 3.3 (Asymptotic expansion around a general critical zero). Let ρ = 1 2 + iρ y and ρ = 1 2 + i(ρ y + δρ y ) be two critical zeros respectively associated to the parameters ∆ and ∆ + δ∆. Then we have, as δ∆ → 0 (and therefore δρ y → 0), (1 + 4ρ y 2 ) 3 -

aδ∆ + f δρ y + c (δρ y ) 2 -b(δ∆)(δρ y ) -d (δρ y ) 3 + O (δ∆) 2 + O δ∆ (δρ y ) 2 = 0, ( 3
1 2 1 0 cos (ρ y log t) (log t) 2 θ 3 e -πt∆ θ 3 e -πt/∆ - 1 t dt √ t , (3.10)
d = 128ρ y 1 -4ρ y 2 (1 + 4ρ y 2 ) 4 - 1 6
1 0 sin (ρ y log t) (log t) 3 θ 3 e -πt∆ θ 3 e -πt/∆ - 1 t dt √ t (3.11)
with θ 3 is defined by (2.3) and ϑ(q) := q d dq θ 3 (q) = j∈Z j 2 q j 2 .

(3.12)

Proof. One substitutes ∆ → ∆ + δ∆, and ρ y → ρ y + δρ y (3.13) into equation (2.6) and expands in Taylor series in powers of small deviations δ∆ and δρ y . Furthermore, the θ 3 -functions appearing in (2.6) are expanded as follows

θ 3 e -πt∆ → θ 3 e -πt∆ -πt(δ∆)ϑ e -πt∆ , θ 3 e -πt/∆ → θ 3 e -πt/∆ + πtδ∆ ∆ 2 ϑ e -πt/∆ . (3.
14) It is now straightforward to obtain the desired result.

It is easy to derive the following asymptotics for δρ y in the case of a critical zero which does not belong to the set of edge zeros. Lemma 3.4 (Asymptotic expansion around a critical non-edge zero). Let ρ = 1 2 + i(ρ y + δρ y ) be a critical zero associated to the parameter ∆ + δ∆ such that the critical zero 1 2 + iρ y associated to the parameter ∆ is not an edge zero. Then we have, as δ∆ → 0 (and consequently δρ y → 0),

δρ y = - a f δ∆ + O (δρ y δ∆) .
Proof. Since ρ is not an edge zero, we have f (ρ y , ∆) = 0 and therefore the leading term of the expansion of the deviation δρ y in small δ∆ is determined by the equation aδ∆ + f δρ y + O (δρ y δ∆) = 0, which leads to our result.

We can also derive from (3.7) the following asymptotics around an edge zero.

Lemma 3.5 (Asymptotic expansion around a critical edge zero). Let ρ = 1 2 + i(ρ * y + δρ y ) associated to the parameter ∆ be a critical zero where ρ * = 1 2 + iρ * y , associated to the parameter ∆ * with ∆ > ∆ * , is a left edge zero. Then we have, as ∆ → ∆ * (and therefore ρ y → ρ * y ),

ρ y (∆) = ρ y (∆ * ) ± - a c √ ∆ -∆ * + 1 2c b - ad c (∆ -∆ * ) + O (∆ -∆ * ) 3/2 , (3.15)
where the sign ± determines the up/down branches of ρ y (∆).

Remark 3.6. In particular, the leading order term of our asymptotics is given, as ∆ → ∆ * , by

ρ y (∆) = ρ y (∆ * ) ± - a c √ ∆ -∆ * + o √ ∆ -∆ * . (3.16)
Here, the ratio -a/c must be positive for the real component ρ y to exist and it was checked numerically for all left edge zeros presented in Table 1 that it is so. Notice that singular expansion of type (3.16) with the critical exponent 1 2 occurs in a mean-field description of classical statistical systems at the second-order phase transition [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF][START_REF] Šamaj | Introduction to the Statistical Physics of Integrable Many-body Systems[END_REF]. An analogous analysis can be done for right edge zeros.

Proof. For a critical edge zero ρ * = 1 2 + iρ * y , the linear term of order δρ y is absent in (3.7) since f (ρ * y , ∆ * ) = 0. Terms on the left-hand side of (3.7) can be classified according to their power in the smallness parameter δ∆ and the leading order can be easily derived. To go to the next order in δ∆, one adds to δρ y in (3.16) a higher-order term αδ∆ with the as-yet-undetermined constant α. Expanding all functions in (3.7) up to the order (δ∆) 3/2 one ends up with the relation 2cα -b + ad/c = 0 which fixes α. It is now straightforward to get our final expression of the asymptotics.

Application and comparison with the direct computations. Let us choose one of the left edge critical zeros in Figure 1, say the edge point 3b with characteristics ∆ * and ρ * y listed in Table 1. After the numerical evaluation of the coefficients a, b, c and d, the expansion (3.15) takes the form

δρ y (∆) = ±4.87411 √ ∆ -∆ * + 22.493 (∆ -∆ * ) + O (∆ -∆ * ) 3/2 .
(3.17)

The fitting of numerical data for the up branch (the + sign) implies the prefactor 4.87412 to √ ∆ -∆ * and the one 22.686 to (∆ -∆ * ), while for the down branch (the -sign) the corresponding prefactors -4.87415 and 22.230 deviate a bit more from the "exact" ones in (3.17), but they are still suitable.

Approximation of Epstein zeta function for small values of ∆

In this part, for small values of the anisotropy parameter ∆, the Epstein function is approximated well by a function which reveals an equidistant distribution of critical zeros along the imaginary axis in the limit ∆ → 0. We start by showing the following asymptotic result as ∆ → 0. Theorem 3.7 (Asymptotics for small values of ∆). As ∆ → 0, we have

ζ (2) (s, ∆) = 1 ∆ 2s 1 π 1 2 -2s Γ 1 2 -s Γ(s) ζ(1 -2s) + 1 ∆ √ πΓ s -1 2 Γ(s) ζ(2s -1) + O(∆). (3.18)
Proof. The sum in (1.1) can be straightforwardly converted to

1 2 ∞ ′ j,k=-∞ 1 (j 2 + ∆ 2 k 2 ) s = ζ(2s) + 1 ∆ 2s ζ(2s) + 2 ∞ j,k=1 1 (j 2 + ∆ 2 k 2 ) s . (3.19) Writing ∞ k=1 1 (j 2 + ∆ 2 k 2 ) s = ∞ k=0 1 (j 2 + ∆ 2 k 2 ) s - 1 j 2s (3.20)
the Epstein zeta function is expressed as

ζ (2) (s, ∆) = 1 ∆ 2s ζ(2s) -ζ(2s) + ∞ j=1 ∞ k=0 f j (k), (3.21)
where

f j (x) = 2 (j 2 + ∆ 2 x 2 ) s .
(3.22)

In the limit ∆ → 0, the difference between successive terms in the sum ∞ k=0 f j (k) is also negligibly small. The sum can be thus treated as an integral according to the Euler-Maclaurin formula [START_REF] Apostol | An elementary view of Euler's summation formula[END_REF]

∞ k=0 f j (k) = ∞ 0 f j (x)dx+ f j (0) + lim M →∞ f j (M ) 2 + [ p 2 ] l=1 B 2l (2l)! lim M →∞ f (2l-1) j (M ) -f (2l-1) j (0) +R j (p), ( 3
|R j (p)| ≤ 2ζ(p) (2π) p ∞ 0 f (p) j (x) dx. (3.24)
f j (x) is an even function of x and therefore its odd derivatives with respect to x vanish at x = 0. Since the derivatives of f j (x) vanish also at x → ∞, it holds that

∞ k=0 f j (k) = 2 ∞ 0 dx (j 2 + ∆ 2 x 2 ) s + 1 j 2s + R j (p). (3.25)
Evaluating the integral

∞ 0 dx (j 2 + ∆ 2 x 2 ) s = 1 ∆ 1 j 2s-1 √ π 2 Γ s -1 2 Γ(s) , (3.26)
we end up with 3.33) (dashed curves) with the ones obtained by using the exact equation (2.6) (the symbol notation is taken from Figure 1).

ζ (2) (s, ∆) = 1 ∆ 2s ζ(2s) + 1 ∆ √ πΓ s -1 2 Γ(s) ζ(2s -1) + ∞ j=1 R j (p). ( 3
By applying the dual relation (see e.g. [START_REF] Blanc | The crystallization conjecture: a review[END_REF])

1 π s Γ(s)ζ(2s) = 1 π 1 2 -s Γ 1 2 -s ζ(1 -2s) (3.28)
we end up with the following asymptotic expansion as ∆ → 0:

ζ (2) (s, ∆) = 1 ∆ 2s 1 π 1 2 -2s Γ 1 2 -s Γ(s) ζ(1 -2s) + 1 ∆ √ πΓ s -1 2 Γ(s) ζ(2s -1) + ∞ j=1 R j (p). (3.29)
To estimate the error term, writing s = s x + is y , we insert f j (x) (3.22) into the bound (3.24) for p = 2 and obtain by a straightforward computation that

∞ j=1 R j (2) ≤ ∞ j=1 2ζ(2) (2π) 2 ∞ 0 f (2) j (x) dx = 8∆ζ(2)ζ(2s x + 1) (2π) 2 s 2
x + s 2 y I(s x , s y ), (3.30) where the integral

I(s x , s y ) = 1 0 [1 -(1 + 2s x )t 2 ] 2 + 4s 2 y t 4 + t 2sx [t 2 -(1 + 2s x )] 2 + 4s 2 y dt (1 + t 2 ) 2+sx (3.31) converges if s x > -1
2 . Note that ζ(2s x + 1) diverges for s x = 0. This completes the proof. Approximate equation for zeros in the small ∆ regime and comparison with direct computations. Let us neglect the error term of order ∆ and consider the approximate equation for zeros {ρ} of the Epstein zeta function in the region of small ∆: In Figure 3, the results for the first critical zeros as the functions of the anisotropy parameter ∆ calculated from this approximative equation (dashed curves) are compared with the ones obtained by using the exact equation (2.6) (the symbol notation is taken from Figure 1). It is seen that the critical zeros obtained from the approximate equation (3.33) agree with the exact ones unexpectedly far away, up to ∆ ≈ 0.5. Roughly speaking, the approximate formula (3.33) works well until approaching an edge zero at which two critical zeros merge, the phenomenon which is out of reach of this formula. This can be seen in the upper right corner of Fig. 3 where the two "exact" curves tend to the right edge point 3a. High accuracy of the approximative values of the first eight critical zeros for the Epstein zeta function is documented for ∆ = 1/ √ 7 ≈ 0.378 in Table 2. As is intuitively expected, the accuracy of the approximative results deteriorates as the value of ρ y increases.

∆ π 2ρ-1 = - Γ 1 2 -ρ ζ(1 -2ρ) Γ ρ -1 2 ζ(2ρ -1) . ( 3
The following result illustrates what we observe on Figure 3 and shows that, in the limit ∆ → 0, the set of critical zeros is equidistant along the imaginary axis, with the spacing π | log(∆/π)| between the nearest neighbors. Lemma 3.9 (Asymptotic behavior of ρ y as ∆ → 0). Let ρ = 1 2 + iρ y be a critical zero associated to the parameter ∆. Then we have lim ∆→0 ρ y (∆) = 0. Furthermore, in the limit ∆ → 0, the asymptotic critical zeros are given by {ρ(n) := 1 2 + ρ y (n)} where

ρ y (n) := πn | log(∆/π)| , n ∈ Z * . (3.34)
Proof. Let us consider the two first orders of the asymptotics in (3.29) when s = 1 2 + iρ y :

ζ (2) (s, ∆) = √ πΓ(-iρ y )ζ(-2iρ y ) ∆Γ 1 2 + iρ y ∆ π -2iρy + Γ(iρ y )ζ(2iρ y ) Γ(-iρ y )ζ(-2iρ y ) . (3.35)
Let us assume that in the limit ∆ → 0 also the component ρ y of critical zeros goes to 0, as is seen in Figure 3. Since

lim ρy→0 Γ (iρ y ) ζ(2iρ y ) Γ (-iρ y ) ζ(-2iρ y ) = -1, (3.36)
the critical zeros are given by (∆/π) 2iρy = 1, in agreement with (3.34). As

√ πΓ(-iρ y )ζ(-2iρ y ) ∆Γ 1 2 + iρ y ∼ ρy→0 - i 2ρ y , (3.37)
the zeta function ζ (2) 1 2 + ρ y , ∆ diverges as | log ∆|/∆ on the curves of critical zeros as ∆ → 0, but this has no impact on the location of its zeros. The proof is complete.

The distance between the nearest neighbors in the asymptotic sequence of critical zeros (3.34) is predicted to be π/| log(∆/π)|; the dependence on the inverse of log(∆) indicates that one has to take extremely small values of ∆ to obtain reliable results. We have performed numerical evaluation of the first few critical zeros with ρ y > 0 for ∆ = 0.0001 by using the exact formula (2.6). The distance between the first and second zeros is 0.375, between the second and third zeros is 0.357, between the third and fourth zeros is 0.347, between the fourth and fifth zeros is 0.340, between the fifth and sixth zeros is 0.337, etc., which means that the spectrum of zeros is almost equidistant as was anticipated. Our asymptotic result (3.34) suggests that the distance should be 0.303 which is a reasonable estimate for the considered (not small enough) value of ∆ = 0.0001. Remark 3.10 (The large ∆ case). With regard to the symmetry ∆ → 1/∆ of basic equations (2.6) for critical and (2.8), (2.9) for off-critical zeros, one can accomplish an analogous analysis in the opposite limit ∆ → ∞, with the result

(π∆) 2ρ-1 = - Γ ρ -1 2 ζ(2ρ -1) Γ 1 2 -ρ ζ(1 -2ρ) . ( 3

.38)

Similarly as in the limit ∆ → 0, the spectrum of critical zeros is equidistant along the imaginary axis in the limit ∆ → ∞, with the spacing π log(π∆) between the nearest neighbors. We observe the following in Figure 1: given an edge zero associated to the parameter ∆ * , there exists δ 0 > 0 such that for all 0 < δ < δ 0 , for a left (resp. right) edge zero, there is no critical zero associated to the parameter ∆ * -δ (resp. ∆ * + δ). In the specific case of left edge zeros, this is caused by the fact that Eq. (3.7) with the numerically verified inequality -a/c > 0 has no real solution for δρ y if δ∆ = ∆ -∆ * < 0, see also Eq. (3.16). Therefore, the only way to have a zero of ζ (2) (s, ∆) corresponding to these values of the parameter ∆ close to an edge zero is to allow the ρ x -component to deviate from its critical value 1 2 . We therefore obtain the following result. Lemma 4.1 (Asymptotic expansion of an off-critical zero around a critical zero). Let ρ = 1 2 + δρ x + i (ρ y + δρ y ) be an off-critical zero associated to the parameter ∆ + δ∆ such that ρ = 1 2 + iρ y is a critical zero associated to the parameter ∆. Then we have, as δ∆ → 0 (and therefore δρ y → 0 and δρ x → 0), the following two equations:

aδ∆ + f δρ y -c (δρ x ) 2 + c (δρ y ) 2 -bδ∆δρ y + 3d (δρ x ) 2 δρ y -d (δρ y ) 3 + o(δρ 3 y ) = 0, (4.1)
δρ x -f + bδ∆ -2cδρ y -d (δρ x ) 2 + 3d (δρ y ) 2 + o(δρ 2 y ) + o(δρ 2 x ) = 0, (4.2)
where the function f is defined by (3.4) and the functions a, b, c and d by equations (3.8)- (3.11). Furthermore if ρ is not an edge zero, then there is no such off-critical zero ρ in its neighborhood.

Proof. We simply use (3.13) in (2.8) and (2.9) and Taylor expanding in powers of small variables δ∆, δρ x and δρ y . Furthermore, if the critical zero ρ is not an edge zero, it holds that f (ρ y , ∆) = 0. The second equation (4.2) with sufficiently small δ∆ and δρ y has the only solution δρ x = 0. In other words, there are no off-critical zeros in the neighborhood of the critical zero which is not an edge zero.

Lemma 4.2 (Asymptotic expansion of an off critical zero around a left edge zero). Let ρ = 1 2 + δρ x + i ρ * y + δρ y be an off-critical zero associated to the parameter ∆ such that ρ * = 1 2 + iρ * y is a left edge zero associated to the parameter ∆ * . Then, we obtain, for ∆ < ∆ * , ∆ → ∆ * (and then δρ y → 0 and δρ x → 0),

δρ x = ± - a c √ ∆ * -∆ + o √ ∆ * -∆ , (4.3)
and

δρ y = - 1 2c b - ad c (∆ * -∆) + o (∆ * -∆) , (4.4)
Remark 4.3. The relations (4.3) and (4.4) are the asymptotic formulas for a curve of off-critical zeros starting from the considered left edge point which are valid for ∆ close to ∆ * and ρ y close to ρ * y .

The ± sign for δρ x in (4.3) means that at each point along a curve of off-critical zeros ρ y (∆) there exist a conjugate pair of solutions ρ x = 1 2 + δρ x and 1 -ρ x = 1 2 -δρ x . A similar analysis can be made for right edge zeros.

Proof. Since ρ * is an edge zero, we have f (ρ * y , ∆ * ) = 0 and it follows that (4.2) is satisfied also for δρ x = 0. Indeed, note that in the first equation (4.1) the term c (δρ y ) 2 , which was dominant in the previous analysis of critical zeros in Section 3.2, has a counterpart with the opposite sign -c (δρ x ) 2 . This latter term becomes dominant when the difference δ∆ = ∆ -∆ * changes its positive sign to the negative one, implying that

a(∆ -∆ * ) -c(δρ x ) 2 + c(δρ y ) 2 + o(δρ 2 y ) + o(δρ 2 x ) = 0, Since δρ x = 0, the second equation (4.2) implies that b(∆ -∆ * ) -2cδρ y -d (δρ x ) 2 + o (∆ * -∆) 2 + o(δρ 2 x ) = 0.
These two equations exhibit the solutions with expansions of type (4.3) and (4.4).

Numerical method to generate off-critical zeros curves. The fact that each curve of off-critical zeros starts/ends at edge points simplifies very much the numerical evaluation of offcritical zeros by using Mathematica. As a function to deal with we take the sum of the squared left-hand side of Eqs. (2.8) and (2.9). Applying the command FindMinimum to this function, the zero is taken as sure if the function value is less than 10 -23 . To avoid escape from a local minimum, one starts from a (say right) edge point and increases ∆ by a tiny amount 0.0001, after few steps the shift can be augmented to 0.001 -0.01. The search for a minimum takes around 60 sec of CPU time on a conventional PC. For integer values of ∆ 2 = 5, 6, 8, . . ., when off-critical zeros can be calculated with a high precision from exact sums of products of Dirichlet L-functions [START_REF] Borwein | Lattice Sums Then and Now[END_REF], our numerical results agree with these analytic predictions by at least 20 decimal digits.

Numerical observations. As is seen in Figure 2, each curve of off-critical zeros joins a pair of critical right (notation "a") and left ("b") edge zeros. As a rule, the ρ y -coordinate of the right edge point is smaller than that of the corresponding left edge point. In the large majority of cases the curves of off-critical zeros go up monotonously when increasing ∆; the only exception from the curves presented in Figure 2 is the curve starting at the right edge point 11a which first goes down in a short interval of ∆-values and then goes up to the left edge point 11b (not in the figure). The intersection of a solid curve of off-critical zeros with a curve of critical zeros (symbols) is not contradictory: the component ρ x varies along the solid curves (not indicated in the figure) while it is constant 1 2 along the critical curves. The vertical dashed lines, pictured at the values 1/ √ 2, 1/ √ 3, 1/ √ 4 and 1/ √ 7 of ∆, correspond, together with ∆ = 1, to the special cases when ζ (2) (s, ∆) factorizes itself into product of a zeta function, a Dirichlet L function and a prefactor function whose zeros (lying on the critical line only) can be determined analytically [START_REF] Borwein | Lattice Sums Then and Now[END_REF][START_REF] Mcphedran | Zeros of lattice sums: 1. Zeros of the critical line[END_REF]. According to the generalized Riemann hypothesis [START_REF] Davenport | On the zeros of certain Dirichlet series I[END_REF][START_REF] Mcphedran | Zeros of lattice sums: 1. Zeros of the critical line[END_REF], ζ (2) (s, ∆) exhibits only critical zeros for these values of ∆. This fact is clearly seen in Figure 2 where no solid curve of off-critical zeros intersects dashed and ∆ = 1 lines, although some of the edge zeros are localized very close to dashed lines. of critical zeros, ρ x = 1 2 at the edge points 2a and 2b. Note that at each point of the curve there exists another solution with the component 1 -ρ x . The intersection of the curve of off-critical zeros (full squares) with ρ x = 1 2 and the curve of critical zeros (open triangles) with ρ x = 1 2 is artificial because the corresponding components ρ x do not coincide.

Comparison of our asymptotic formulas and analytic data obtained directly. A check of the asymptotic formulas (4.3) and (4.4) for the curve of off-critical zeros close to a left edge zero was made for the edge point denoted as 3b in Figure 1, similarly as in the previous case of the expansion formula (3.17) for the curve of critical zeros. The expected dependences

δρ x (∆) = ±4.87411 √ ∆ * -∆, -δρ y (∆) = 22.493 (∆ * -∆) (4.5)
are reproduced very well by fitting numerical data, namely the prefactor obtained for δρ x equals to 4.87412 and the prefactor for -δρ y is estimated to 22.498. The agreement of the asymptotic relations (4.5) (dashed lines) with the numerical data (open symbols) is pictured in Figure 5.

Real off-critical zeros

In this section, we are interested in real off-critical zeros. It has already been mentioned that the critical zero numbered by 1 in Figure 1, lying on the ∆-axis, is a left edge zero because the curve of critical zeros passes across the ∆-axis into the lower quadrant in a reflection-symmetrical way.

Numerical observations. The imaginary part of this point ρ * y = 0 as well as its associated Proof. It has been shown in [START_REF] Faulhuber | Optimal Gabor frame bounds for separable lattices and estimates for Jacobi theta functions[END_REF][START_REF] Montgomery | Minimal Theta Functions[END_REF] that, for all t > 0, the function ∆ → θ 3 e -πt∆ θ 3 e -πt/∆ is strictly decreasing on (0, 1). Furthermore, for ∆ = 1 we have

-2 + 1 0 θ 3 e -πt 2 - 1 t dt √ t ≈ -1.9501325 < 0.
To study the opposite ∆ → 0 limit, we recall that Z(s, ∆) is related to the Epstein zeta function where γ = 0.5772156 . . . is the Euler-Mascheroni constant. Thus,

Z 1 2 , ∆ ∼ ∆→0 - log ∆ √ ∆ → +∞.
It follows that ∆ * c is unique and the proof is complete.

As is seen in Table 1, the numerical solution of this equation is ∆ * c ≈ 0.141733239663887. According to our numerical observations:

(1) it turns out that the curve of off-critical zeros going out of the critical edge zero 1 stays on the ∆-axis, i.e. ρ y = 0. (2) As concerns the coupled equations for off-critical zeros (2.8) and (2.9), the second one is automatically fulfilled for ρ y = 0 while the first one implies

- 1 1 + 2δρ x + 1 1 -2δρ x + 1 0
cosh (δρ x log t) θ 3 e -πt∆ θ 3 e -πt/∆ -1 t dt √ t = 0.

(4.7) This equation has two conjugate real solutions for δρ x and -δρ x (or, equivalently, ρ x and 1 -ρ x ), only if 0 ≤ ∆ < ∆ * c as observed in Figure 6. (3) The value of ρ x is 1 2 at the critical edge point 1 corresponding to ∆ = ∆ * c . As observed again in Figure 6, decreasing the value of ∆ < ∆ * c , the two values of ρ x split and tend to the borders 0 (down branch) and 1 (up branch) of the critical strip in the limit ∆ → 0.

Heuristic derivation of ∆ * c from small ∆ approximation. It is useful to test how Eq. (3.32), which is accurate but certainly only approximate for complex zeros (see Table 2 We checked that ∆ * c evaluated by using this analytic relation coincides with the previous numerical estimate ∆ * c ≈ 0.141733239663887 obtained from the exact Eq. (4.6) by at least 22 decimal digits; to go further a computer facility more powerful than the one at our disposal is needed. This indicates that the result (4.11) might be exact which is difficult to prove directly by using (4.6). We have stated the corresponding open problem as Conjecture 1.1.

Approximation in the ρ x , ∆ → 0 regime and comparison with numerical data. As concerns the accuracy of the real off-critical zeros implied by Eq. (4.8), for various values of ∆ from the interval (0, ∆ * c ) they coincide with the ones obtained from the exact Eq. (4.7) up to 27 decimal digits, which supports the hypothesis that the real off-critical zeros generated from Eq. (4.8) are exact.

For the down branch in Figure 6, to obtain the asymptotic tendency of ρ x to 0 as ∆ → 0 by using the relation (4.8), one writes δρ x = ρ x -1 2 and expands the right-hand side in small ρ x :

- Γ 1 2 -ρ x ζ (1 -2ρ x ) Γ ρ x -

1 .

 1 The zeta functions and their zeros.

Figure 1 .

 1 Figure 1. Numerical results for the critical zeros in the region of 0 < ∆ ≤ 1 and 0 ≤ ρy ≤ 21. Zeros lying on the same (closed or open) curve are denoted by a common open symbol (circle, square, triangle,. . .). The right and left edge points are denoted as 2a, 3a, 4a,. . . and 1, 2b, 3b, 4b,. . ., respectively.

Figure 2 .

 2 Figure 2. A fragment of Figure 1. Numerical results for solid curves of off-critical zeros interconnecting the critical right (notation "a") and left ("b") edge zeros. The evolution of the component ρx along the solid curves of off-critical zeros is not indicated. The vertical dashed lines correspond to the values of ∆ 2 = 1, 2, 3, 4, 7 with no off-critical zeros.

3. Zeros on the critical line 3 . 1 .

 31 Critical edge zerosLet us first comment on Figure1. The critical zeros of the Epstein zeta function ζ(2) (s, ∆), calculated numerically by using Eq. (2.6), are represented by open symbols in Figure1. The set of zeros lying on the same (closed or open) curve are denoted by a common open symbol (circle, square, triangle,. . .). There is an infinite sequence of loop circuits in the left-down corner of the figure with the common point at the origin (0, 0) which are not drawn. With regard to the symmetry ∆ → 1/∆, each critical zero in the considered interval ∆ ∈ (0, 1) has a counterpart in the complementary interval (1, ∞). To maintain high accuracy of the results (8-20 decimal digits), only critical zeros with the component ρ y smaller than 21 are presented. The numerical evaluation of one critical zero by using Mathematica takes around 5 seconds of CPU time on a conventional PC.

. 7 )= π 1 0 √ 0 √

 700 where the function f = f (ρ y , ∆) is defined by (3.4) and the other prefactor functions a, b, c and d, depending on (ρ y , ∆), are given by a t cos (ρ y log t) θ 3 e -πt∆ 1 ∆ 2 ϑ e -πt/∆ -θ 3 e -πt/∆ ϑ e -πt∆ dt, (3.8)b = π 1 t sin (ρ y log t) (log t) θ 3 e -πt∆ 1 ∆ 2 ϑ e -πt/∆ -θ 3 e -πt/∆ ϑ e -πt∆ dt, (3.9)

  .23) where p ≥ 2 is an integer , [• • • ] denotes the integer part, {B 2l } are Bernoulli numbers and the absolute value of the error term R j (p) is bounded by

Figure 3 .

 3 Figure 3. The comparison of the results for the first critical zeros as the functions of ∆ calculated from the approximative equation (3.33) (dashed curves) with the ones obtained by using the exact equation (2.6) (the symbol notation is taken from Figure1).

For the critical zeros ρ = 1 2 +

 2 iρ y , this equation takes the form ∆ π 2iρy = -Γ (-iρ y ) ζ(-2iρ y ) Γ (iρ y ) ζ(2iρ y ) . (3.33)

4. Zeros off the critical line 4 . 1 .

 41 Generation of off-critical zeros from critical edge zeros

Figure 4 Figure 4 .

 44 Figure 4. Numerical results for the curve off-critical zeros (full squares) between the right edge point 2a (lying on the curve of critical zeros represented by open squares) and the left edge point 2b (lying on the curve of critical zeros represented by open triangles). The evolution of the component ρx along the curve is indicated by short lines. At each point of the curve there exists another solution with the component 1 -ρx.

Figure 5 . 6 )

 56 Figure 5. The log-log plot of numerical and analytic data for the curve of off-critical zeros going from the left edge denoted as 3b in Figure 1. Numerical dependences of δρx (open circles) and -δρy (open squares) on the small deviations from the edge point ∆ * -∆ are compared with the analytic predictions (4.5) represented by dashed lines.

Lemma 4 . 4 .

 44 Equation (4.6) admits a unique solution ∆ * c on (0, 1].

ζ ( 2 ) 1 ∆ s 1 π 1 2 -s Γ 1 2 - 1 2 1 2Z

 2111211 (s, ∆) via equation (2.1) and the small-∆ behavior of ζ(2) (s, ∆) is given by the asymptotic relation(3.29) where the p = 2 error term of the order O(∆) can be neglected in the limit ∆ → 0. Thus one arrives at the asymptotic relationZ(s, ∆) ∼ ∆→0 s ζ(1 -2s) + 1 ∆ 1-s π -s Γ s -1 2 ζ(2s -1)which exhibits the required duality symmetry (2.4). Consequently, lim s→ (s, ∆) = γ -log(4π∆) √ ∆ ,

4 . 8 )Figure 6 .

 486 Figure 6. The dependence of the ρx-component of real off-critical zeros (ρy = 0) on the anisotropy parameter ∆, calculated numerically by using equation (4.7). As explained in the text, ρx = 1 2 at ∆ * c ≈ 0.141733239663887 and there are two conjugate solutions ρx and 1 -ρx for 0 < ∆ < ∆ * c . As ∆ → 0, the two values of ρx tend to the borders 0 (the down branch) and 1 (the up branch) of the critical strip. The inset concerns the down branch and the region of small ∆ where numerical data for ρx versus ∆ (open circles) satisfy the asymptotic relation ρx ∼ (3/π)∆ (dashed line).

Table 1 .

 1 The coordinates of edge critical zeros appearing in the region of 0 < ∆ ≤ 1 and 0 ≤ ρy ≤ 21 in Figure1.

	edge point	∆ *	ρ * y
	1	0.141733239663887	0
	2a	0.798382429865856 9.17479405815734
	2b	0.950672823506692 13.5092488680816
	3a	0.578095740200051 11.2961629757333
	3b	0.688797339793161 12.7134082666419
	4a	0.708261915413478 14.3461052173020
	4b	0.810471985748564 17.3035168808027
	5a	0.577833206956181 15.7904269230734
	5b	0.625830051933379 16.7721421891791
	6a	0.397042034784957 14.9386821841068
	6b	0.499955572107973 16.3629327845743
	7a	0.634086781531453 17.8588321271621
	7b	0.690295752437308 19.4462462865857
	8a	0.356573014664413 16.4816098051657
	8b	0.375454386384881 16.8495675287149
	9a	0.337272867689201 18.5674591768417
	9b	0.374296061779980 19.3629136770424
	10a	0.560652822542094 19.8540419510498
	10b	0.577320038404815 20.4237238736290
	11a	0.896821462590355 19.0008766867965
	12a	0.778481639573212 20.0611304186419
	12b	0.789270563104711 20.2667094854061
	13a	0.578437965650995 20.8178435639014
	14a	0.309679721075915 20.1102459521285

Table 2 .

 2 The comparison of the exact and approximate values of the first 8 zeros for ∆ = 1/ √ 7.

	zero #	ρ y exact	ρ y approx
	1	1.133090035457 1.133090358285
	2	4.475738283729 4.475726461185
	3	6.845491712491 6.845712742060
	4	7.931630248198 7.930996972746
	5	10.19781031911 10.20336832640
	6	11.16018454312 11.14537554655
	7	12.48960334303 12.51829228147
	8	14.13472514173 14.05004856679
	Remark 3.8. The expansion (3.29) then has the meaning of a systematic Laurent series expan-
	sion in ∆ provided that ℜ(s) > 0.

& & & & & & &

Acknowledgment

Ladislav Šamaj and Igor Travěnec acknowledge the support received from VEGA Grant No. 2/0092/21.

Consequently, again as ∆ → 0

As is seen in the inset of Figure 6, numerical data (open circles) agree well with this analytic prediction (dashed line).

Remark 4.5. With regard to the symmetry ∆ → 1/∆ of basic equations for zeros of ζ (2) (s, ∆), there exists a pair of continuous curves of real off-critical zeros also for each ∆ > 1/∆ * c ≈ 7.055507955448192.


Refet A Yalcin

Cagatay Haratoka

David Babonneau

Sophie Camelio

Karl Joulain

Jeremie Drevillon

The effect of size and shape-dependent optical properties on colored radiative cooling in metal island films

We study visible and infrared optical absorption of metal island films that are commonly used for applications such as color filters and sensors. Calculations were performed using the T-matrix method and compared to experimental transmittance measurements at visible spectrum on a real sample. While the radiation transfer solution through these films often assumes (i) no plasmon coupling, (ii) spherical particles, (iii) monodispersed size particles and/or (iv) bulk refractive index, we carefully study the influence of these assumptions on spectral transmittance. After a validation of the method, the resulting performance in terms of solar energy absorption, color lightness and hue are evaluated on a radiative cooling structure. Alternative designs that provide stronger coloration and thermal performance are proposed and discussed.

Introduction

Radiative cooling is the radiative heat exchange between an object and the space which is at 3K. As demonstrated earlier [START_REF] Raman | Passive radiative cooling below ambient air temperature under direct sunlight[END_REF], an object can passively cool below ambient temperature even under direct sunlight. In order to perform efficiently, radiative coolers are expected to be spectrally selective. For radiative cooling applications, objects should be reflective in the solar spectrum and emissive in the infrared, especially at wavelengths between 8 and 13 µm where the atmosphere is transparent. As a result of high solar reflectance, effective radiative coolers have whitish colors.

However, due to aesthetic purposes, different works have considered the design of colored radiative coolers [START_REF] Yalcin | Colored Radiative Cooling Coatings with Nanoparticles[END_REF][START_REF] Yalçın | Colored radiative cooling coatings with fluorescence[END_REF]. The challenge to introduce color into radiative coolers lies in performing a trade-off between thermal performance and color intensity (lightness) [START_REF] Yalcin | Colored Radiative Cooling Coatings with Nanoparticles[END_REF]. Color is introduced by absorption of light in the visible spectrum which also increases the absorbed solar heat flux and decreases the thermal efficiency of radiative coolers. Having a narrow absorption peak is therefore an efficient way in order to obtain low solar absorption coloration. Metal nanoparticles (NPs) are good candidates for this application [START_REF] Yalcin | Colored Radiative Cooling Coatings with Nanoparticles[END_REF]. Nevertheless, effects of localized surface plasmons and plasmon coupling, that can negatively affect thermal and visual performance of the radiative coolers, should be further investigated.

Metal island films consisting of metal particles dispersed on a dielectric surface are known to exhibit absorbance peaks in the visible range due to the excitation of localized surface plasmons when interacting with light. They have a wide application range such as spectrally selective solar light absorbing [START_REF] Fabijanić | Plasmonic Nanoparticles and Island Films for Solar Energy Harvesting: A Comparative Study of Cu, Al, Ag and Au Performance[END_REF], filtering [START_REF] Freund | Clusters and islands on oxides: from catalysis via electronics and magnetism to optics[END_REF], sensing applications [START_REF] Hoa | Towards integrated and sensitive surface plasmon resonance biosensors: A review of recent progress[END_REF] and catalyzing the growth of unidimensional objects such as carbon nanotubes [START_REF] Takagi | Single-Walled Carbon Nanotube Growth from Highly Activated Metal Nanoparticles[END_REF].

While simulating the light transfer through metal NPs in a plane, because of the complexity and computational intensiveness of the problem, several assumptions are usually performed such as;

(i) neglecting the coupling effects between neighboring particles, (ii) assuming spherical particles, (iii) using monodispersed particles, (iv) using bulk refractive index. The present study aims to identify the accuracy of these assumptions while simulating light transfer in a plane-parallel slab of layered silver NPs that are small (typically below 10 nm in radius). Finally, an alternate approach to achieve colorization introducing thin film filters to radiative coolers is proposed for a better thermal and visual performance.

Background

2.1.Solution of Maxwell's equations

Solution of Maxwell's equations for a single sphere under plane wave excitation was developed by Gustav Mie [START_REF] Mie | Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen[END_REF] and known as Lorenz-Mie theory. The solution was later extended [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF] for cylinders and spheroids. In a medium where the NPs are well-separated from each other the socalled "independent scattering assumption" can be applied [START_REF] Mishchenko | Radiative transfer theory verified by controlled laboratory experiments[END_REF]. Under this assumption, superposition of the absorption cross-sections of separate particles that are calculated by Lorenz-Mie theory can be superposed to find the total absorption cross-section of the domain. This approach does not account for the coupling effects among the neighboring particles and particle, surface interaction.

Lorenz-Mie theory applies to a single particle in vacuum. Electromagnetic wave transport in heterogeneous structures, such as NP layers is also governed by Maxwell's equations. Solution of these equations in complex structures can be obtained by various methods including the finite difference time domain (FDTD) [START_REF] Yee | The finite-difference time-domain (FDTD) and the finite-volume timedomain (FVTD) methods in solving Maxwell's equations[END_REF], the discrete dipole approximation (DDA) [START_REF] Draine | Discrete-Dipole Approximation For Scattering Calculations[END_REF], the finite element (FEM) [START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF], and the superposition T-matrix [START_REF] Mackowski | A multiple sphere T-matrix Fortran code for use on parallel computer clusters[END_REF] methods. In all these methods the problem gets computationally more expensive when metallic NPs are involved, distance between them gets smaller and the number of particles as well as size of the domain with respect to the incident wavelength gets larger. As a result, a proper method in terms of accuracy, computational resource usage and computational time should be selected. The T-matrix method is known to be an effective method compared to the others [START_REF] Baran | Calculation of the single-scattering properties of randomly oriented hexagonal ice columns: a comparison of the T-matrix and the finitedifference time-domain methods[END_REF]. Some software based on the T-matrix method like Smuthi [START_REF] Egel | SMUTHI: A python package for the simulation of light scattering by multiple particles near or between planar interfaces[END_REF] and CELES [START_REF] Egel | CELES: CUDA-accelerated simulation of electromagnetic scattering by large ensembles of spheres[END_REF] are also profited from GPU acceleration that bring down computation times up to 100-fold for systems with large number of particles.

2.2.Size dependent refractive index

A common approach while setting the refractive index of NPs is to use bulk refractive index. For NPs under 20 nm radii, the dielectric function of metals becomes size dependent and as the particle size gets smaller, the refractive index of the NP deviates more from the bulk material's optical properties [START_REF] Herrera | Determination of plasma frequency, damping constant, and size distribution from the complex dielectric function of noble metal nanoparticles[END_REF]. Size dependency accounts for both intraband transitions (free electrons) and interband transitions (bound electrons). The modification that arises from bound electrons is significant only for NPs smaller than 2 nm as pointed out by Herrera et al. [START_REF] Herrera | Determination of plasma frequency, damping constant, and size distribution from the complex dielectric function of noble metal nanoparticles[END_REF]. Since NPs smaller than 2 nm cannot be determined accurately with the applied methods in this study, interband transitions were not considered. Only free electron contribution was considered in this study for the size dependent dielectric function which is given by [START_REF] Herrera | Determination of plasma frequency, damping constant, and size distribution from the complex dielectric function of noble metal nanoparticles[END_REF]:

𝜀(𝜔, 𝑟 𝑒𝑓𝑓 ) = 𝜀 𝑏𝑢𝑙𝑘 (𝜔) -𝜔 𝑝 2 ∑(-1) 𝑛 ∞ 𝑛=1 ( 𝑖𝜔𝐴𝜈 𝐹 𝑟 𝑒𝑓𝑓 ) 𝑛 (𝜔 2 + 𝑖𝜔𝛾 𝑓𝑟𝑒𝑒 ) 𝑛+1 (1)
where 𝑟 𝑒𝑓𝑓 and 𝜔, 𝜔p, γfree, A, υf are effective radius of the particle, electromagnetic wave frequency, plasma frequency, damping constant (for noble metals A  0.8), and Fermi velocity respectively.

2.3.Radiative cooling

Main

𝑃 𝑠𝑜𝑙𝑎𝑟 ′′ = ∫ 𝐼 𝑆𝑜𝑙𝑎𝑟 (𝜆) ∞ 0 𝛼(𝜆, 𝜃) 𝜃=0 𝑑𝜆 (4)
where,

𝜖 𝑎𝑡𝑚 (𝜆, 𝜃) = 1 -𝜏 𝑎𝑡𝑚 (𝜆) 1/ cos(𝜃) (6)
Here, τatm() is the atmospheric transmittance, Isolar is normal direct irradiance, Ibb,Tsur is the blackbody intensity at the surface temperature. Table 1 lists the investigated cases of the study; (i) assuming oblate spheroid and (ii) spherical particles (iii) using monodispersed particles, and (iv) assuming bulk refractive index. The measured spectral transmittance of the real sample is compared with the calculated transmittance that are accounting for the defined cases. Finally, an additional case study, that weakens the coupling effects between neighboring particles by changing the vertical positions of particles and increasing the minimum distance, cmin, between them, was carried out.

L LAg SiO2 L Al2O3

3.2.HAADF-STEM image digitization

The location and size of the particles were digitized by using circular Hough transform [START_REF] Yuen | Comparative study of Hough Transform methods for circle finding[END_REF]. and in-plane size of the particles were identified by digitizing the HAADF-STEM image with surface area L x L using circular Hough transform [START_REF] Yuen | Comparative study of Hough Transform methods for circle finding[END_REF]. Moreover, in addition to the in-plane size distribution, quantitative information about the height of the particles can be obtained from the intensity distribution in the HAADF-STEM image [START_REF] Babonneau | Gold and silver nanoparticles embedded in dielectric-capping layers studied by HAADF-STEM[END_REF]. A correlation exists between a and c as 𝑐 = 𝑓 𝑐 (a) = 0.36 𝑎 + 1.27 nm if a > 2.5 nm (7) 𝑐 = 𝑓 𝑐 (a) = 𝑎 otherwise

3.3.Monodisperse assumption

A single volume averaged particle with periodic boundary conditions is considered to represent all of the metal islands. This approach can be preferred since modelling a high number of spheroids is computationally expensive. The volume equivalent average radius rvol, can be calculated as;

𝑟 𝑣𝑜𝑙 = √ ∑ (𝑎 𝑘 • 𝑎 𝑘 • 𝑐 𝑘 ) 𝑁 𝑝 𝑘=1 𝑁 𝑝 3 (8)
Here, the in-plane area (Lvol x Lvol) of the periodic domain was found by keeping the silver volume per unit area constant which yielded;

𝑟 𝑣𝑜𝑙 • 𝑟 𝑣𝑜𝑙 • 𝑓 𝑐 (𝑟 𝑣𝑜𝑙 ) 𝐿 𝑣𝑜𝑙 𝑥 𝐿 𝑣𝑜𝑙 = ∑ (𝑎 𝑘 • 𝑎 𝑘 • 𝑐 𝑘 ) 𝑁 𝑝 𝑘=1 𝐿 𝑥 𝐿 (9)
Following Equations ( 7) and (8), the volume equivalent average radius, 𝑟 𝑣𝑜𝑙 = 5.92 nm and corresponding short axis c = fc(rvol) = 3.40 nm was found. Side of the periodic square domain was also found from Equation (9) as 𝐿 𝑣𝑜𝑙 = 7.66 nm for the monodisperse assumption case calculations.

3.4.Methods of solution

In the current study, we implemented open-source Python code Smuthi [START_REF] Egel | SMUTHI: A python package for the simulation of light scattering by multiple particles near or between planar interfaces[END_REF] for utilizing the Tmatrix method. The script creates the metal islands by reading the spatial coordinates, [xk, yk, zk] and radius, ak. Then short axis of spheroids, ck, (Equation 7) and size dependent refractive indices of particles (Equation 1) are calculated. Finally, the script solves Maxwell's equations using the software Smuthi and calculates transmittance and reflectance. A recently developed approach [START_REF] Theobald | Simulation of light scattering in large, disordered nanostructures using a periodic T-matrix method[END_REF] to implement periodic boundary condition capability to the T-matrix method is used to investigate the effect of monodisperse assumption.

The T-matrix method cannot solve particles with intersecting surfaces. Figure S2 illustrates the geometries which T-matrix can and cannot be used to solve. Intersecting spheroids are removed from the digitized dataset before utilizing the T-matrix method. Particles that are smaller than 2 nm cannot be determined accurately with the applied methods in this study, therefore they are removed too. After the circular Hough transform digitizes the HAADF-STEM image, the particle data list is ordered from largest sized circle to smallest. Then starting from the largest one, each circle is controlled with other circles to see if there is an overlap. In case of an overlap, the smaller one of the overlapped particles is omitted. With this algorithm, larger circles are prioritized while removing the overlapping circle. Figure 3(c) illustrates the spatial particle distribution, where the intersected and small ones are omitted, the omitted 8 particles constitute less than 0.5% of the particles in the disk in terms of mass, which is insignificant. The finalized digitized particulate disk (Figure 3(c)) is considered while calculating the spectral transmittance using the T-matrix method. Particles are considered sitting on the silica layer as illustrated in Figure 1, therefore particle-substrate interaction is accounted for. We further calculate extinction cross sections of different sized particles on silica and alumina. Since the surrounding (capping) medium is also alumina, the simulations considering particles sitting over alumina disregards particle-substrate interaction. Figure S4(c) illustrates the low difference between two approaches and we attribute weak particle-substrate interaction to low refractive index mismatch between silica and alumina at the spectrum of interest.

We performed a preliminary study to determine which maximum spherical harmonics order, lmax, to use in our study. Light transfer through the original monolayer design was solved at different lmax, from 1 to 5. Their root mean square errors are found for 1 to 4 respectively; 6.8%, 5.5%, 2.2%, 1.2% with respect to lmax = 5. lmax > 5 cases could not be investigated due to excess RAM requirement of 824 particles. lmax = 4 is selected and applied for each case in this study since the error it produced was found to be low enough.

Transmittance and reflectance of the metal island layer in Al2O3 that is presented in Figure 3(c

𝑃(𝑟, 𝑧) = 𝑃 0 [1 -𝑒 -2𝑟 2 𝜔 2 (𝑧) ⁄ ] (10)
where

𝑃 0 = 1 2 𝜋𝐼 0 𝜔 0 2 (11)
In that case, r = ω, the power passed through the circle becomes

𝑃(𝑧) 𝑃 0 = 1 -𝑒 -2 = 0.8647
It means that 13.53% of the gaussian beam's intensity pass through from outside of the circle that contains metal islands (Figure 3(c)). Therefore, a correction is applied to the spectral transmittance, Tmatrix,, and reflectance Rmatrix,, that are calculated by the T-matrix method.

Corrected spectral transmittance of particulate layer, T, and reflectance, R, are calculated as;

𝑇  = 𝑇 𝑚𝑎𝑡𝑟𝑖𝑥, 1 -𝑒 -2 + (1 - 1 1 -𝑒 -2 ) ( 12
)
𝑅  = 𝑅 𝑚𝑎𝑡𝑟𝑖𝑥, 1 -𝑒 -2 (13)
For thickness with orders of magnitude larger than the incident wavelength, there are instabilities at the T-matrix method due to interference. As a result, the ray tracing approach detailed in [START_REF] Siegel | Thermal radiation heat transfer[END_REF] was adopted to account for the boundary reflection effects. In this approach, we consider 1D ray tracing since (i) scattering is small, at least with one order of magnitude, compared to absorption as presented in Figure S4(d) and (ii) the appearance of the sample is not translucent but clear.

Sample reflects the light specularly as presented in Figure S3. At the solar spectrum, we consider a 5-layer model; air, alumina, pigmented layer, silica, and silver. Boundary reflections between air-alumina and silicasilver found by Fresnel formalism. Spectral absorption in these nonscattering layers is calculated by Beer-Lambert law, Aλ = 1 -exp(-4••kλ•L/λ). Here kλ is the spectral imaginary refractive index of the layer. On the other hand, absorption in the non-homogenous pigmented layer and reflection on its boundaries are calculated from the T-matrix calculations and then implemented to the ray tracing calculations. Between 2.5 m -20 m spectra, by running the T-matrix simulations we found pigmented layer neither absorbs nor reflects the light; therefore, we disregard particles and consider 4-layer structure for emissivity calculations; air, alumina, silica, and silver. Bulk components of the refractive index of silver and silica were obtained from Wu et al. [START_REF] Wu | Intrinsic Optical Properties and Enhanced Plasmonic Response of Epitaxial Silver[END_REF] and Kitamura et al. [START_REF] Kitamura | Optical constants of silica glass from extreme ultraviolet to far infrared at near room temperature[END_REF] Refractive index of alumina was retrieved from ellipsometric measurements on an Al2O3 thin film grown on a silicon substrate.

Material and Methods

4.1.Sample Preparation, Materials, and Characterization

In the present study, silver island film deposited onto amorphous substrates (carbon-coated copper grids for transmission electron microscopy and fused silica for optical measurements as sketched in Figure 1(a)) were fabricated. Silver NPs capped with Al2O3 layers were prepared at 200 °C in an ultrahigh vacuum chamber (2×10 -8 mbar before deposition) by electron-beam evaporation [START_REF] Babonneau | Gold and silver nanoparticles embedded in dielectric-capping layers studied by HAADF-STEM[END_REF].

In the first step, a thin metal layer was grown (nominal thickness t = 3 nm, deposition rate 0.03 nm/s) with the aim of forming the metal island by a Volmer Weber growth mechanism. Then, after a delay time of 30 min, the nano islands were covered with a 20-nm-thick Al2O3 capping layer (deposition rate 0.05 nm/s).

The morphology (size and shape) and organization (interparticle distance) of the resulting capped NPs were investigated from in-plane views by high-angle annular dark field scanning transmission electron microscopy (HAADF-STEM) performed with a JEOL 2200FS microscope operated at 200 kV (probe size of 1.5 nm and inner collection angle of 50 mrad). Structural information retrieved from the analysis of the plane-view HAADF-STEM images were used as input parameters to simulate optical transmittance spectra and radiative cooling properties.

4.2.Transmittance Measurements.

The far-field optical response of the capped NPs grown on fused silica substrates (Figure 1(a))

was investigated by transmittance measurements at normal incidence (specular beam with diameter ∼3 mm, forward direction) in the spectral range between 350 nm and 800 nm with a 5 nm monochromator step. Integrating sphere is not used since the material is not translucent but transparent. After coating it with silver, we observed that it also reflects the light specularly as we presented in Figure S3.

Results and Discussion

5.1.HAADF-STEM image digitization and post process

Providing the accurate spatial coordinates and particle size is crucial while simulating the spectral transmittance of the layer composed of silver NPs. From the initial digitization to the final state as the data is provided to the T-matrix method, the following steps are taken. First, the HAADF-STEM image is digitized (Figure 3(a)). Then, the particles that lie outside of the inner tangent circle are removed (Figure 3

5.2.Radiative cooling calculations

The calculated spectral emittance of the radiative cooler (Figure 1 Figure 5. Spectral atmospheric transmittance, emittance and blackbody intensity at T = 300 K.

5.3.Effect of shape, size distribution and using size dependent refractive index

Here, we investigated the effects of aforementioned assumptions listed in Table 1. For this, we compared the calculated transmittances of different cases with experimental measurements of the real sample, which exhibits a dip in transmittance located around  = 515 nm (Figure 6). The optical model of the approximated cases was constructed based on the real sample design that is presented in Figure 1(a). First, we investigated the effect of using spherical particles instead of spheroids on spectral transmittance. Figure 6 shows that a deeper dip in transmittance occurs for the "sphere assumption" compared to the "real sample", which can be attributed to an overestimation of the total amount of silver when spheres are considered. Furthermore, we can see that using spheres instead of oblate spheroids causes a significant blueshift of the surface plasmon resonance, as expected [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF].

Using the bulk refractive index instead of the size dependent one is the common practice since defining each sphere in the design with a different refractive index is a tedious process. The result of the "bulk refractive index assumption" case study is also presented in Figure 6 and Table 2, which show that use of size dependent refractive index is required for better accuracy, especially at the spectrum where transmittance dip was observed. The wiggles around the dip disappeared and absorption at the peak damped after the utilization of size dependent refractive index as expected [START_REF] Kreibig | Optical Properties of Metal Clusters[END_REF]. index compared to size dependent ones, which results in wiggles at transmittance calculations. We attribute these wiggles to the limited number of spheres considered in simulations. Especially absorption is expected to be dominated by the particles at the center of the disk where gaussian beam's intensity is stronger. It is worth to note that, for the spectrum between 430 and 460 nm, the T-matrix method did not converge for the bulk refractive index assumption case and the data among this portion was linearly interpolated.

Figure 6.

Experimental and calculated transmittance spectra at normal incidence. Numerical results were obtained by considering the cases listed in Table 1.

We observed a considerable difference of spectral transmittance between the monodisperse case and the real case. Monodisperse distribution assumption predicted a narrower and a deeper dip at transmittance. Here, we should note that the monodisperse case was not only considering a single sized spheroid that disregards particle size distribution, but also assumes that the location of each particle is well aligned like a photonic structure.

Spheroid particle assumption resulted in a relatively accurate result with a root mean square error less than 3%. The location and width of the transmittance dip compared to the real case is also promising. The outcomes show that spheroid assumption is acceptable and enables the experimental data to be fairly well-reproduced by the T-matrix method. It is worth remembering that the spheroid assumption should be applied with proper particle size distribution and size dependent refractive index.

5.4.Effect of plasmonic coupling

To demonstrate the significance of the coupling effect, we performed a case study, in which we first calculated the absorption cross-section of the digitized particle disk (Figure 3

∑

As described briefly in the problem statement, an additional case study is performed that gradually weakens the coupling effect to observe its effect on spectral characteristics like reflectance and color. 6 different "center to center inter-particle distance with the nearest neighbor", cmin, are chosen, as 10, 15, 20, 25, 30, 35 nm. Then, an algorithm that loops through each particle and shifts them upwards through the light source if they are close to another particle more than cmin. Their new location becomes [xk, yk, zk+ cmin]. Algorithm starts with the digitized particle disk (Figure 3(c)) and spheres can shift to upper levels many times if they are still close to a particle more than cmin. This process is finished until it is ensured that no particle is close to the other more than cmin. In our previous study [START_REF] Simonot | Evolution of plasmonic nanostructures under ultra-low-energy ion bombardment[END_REF], we showed that a post-deposition treatment using an ultra-low energy ion bombardment (<100eV) can be applied to obtain a dilute layer. This method can be applied after each silver deposition, then make a capping layer of dielectric layer then again silver deposition followed by plasma treatment. This procedure can be repeated until each dilute layer is obtained.

Conclusion

This study demonstrates the significance of neighboring particles' coupling effects to the spectral absorption of metal islands. As an application, we investigated the metal islands as color introducing filters to radiative coolers. We found that considering multilayers of dilute silver NPs instead of a single dense layer was more efficient for introducing color to radiative cooling structures in terms of thermal efficacy and color implementation. Furthermore, we investigated the commonly considered assumptions while modeling electromagnetic radiation transfer through metal islands. The studied assumptions were using (i) oblate spheroid and (ii) spherical particles (iii) considering monodisperse size distribution, and finally (iv) using bulk refractive index for silver NPs. The effects of these assumptions were identified visually from spectral transmittance graphs and numerically by calculation of root mean square errors. We found that using these assumptions resulted in errors that can exceed 10% and therefore should be avoided for applications requiring a high accuracy. Finally, the contact angle of the NPs was set to 180° (nonwetting) in this study. In a future study, a realistic contact angle for spheroids can be investigated for a more accurate model.

Figure 1 Figure 1 .

 11 Figure 1 illustrates the silver metallic islands exposed to visible light investigated in the present study and the different parameters considered. The structure was characterized by (i) shape of the particles, (ii) size distribution of the particles, (iii) the particle spatial arrangement, (iv) the particle refractive index, and (v) the continuous medium refractive index at wavelength λ.

Figure 2 3 .Figure 2 .

 232 Figure 2. HAADF-STEM image of the sample.

)

  are calculated by exposing the metal island with a gaussian beam having a beam waist, ω0, equal to the particulate disk radius, L/2 = 247 nm. The particles outside of the inscribed circle are not considered in order to reduce computational load. The resulting change in particles' spatial distribution can be observed by comparing Figure3(a) and Figure3(b). Considering power of a gaussian beam passing through a circle of radius r in the transverse plane at position z is:

4 .Figure 3 .

 43 Figure 3. (a) Digitized HAADF-STEM image, (b) the NPs inside the inner tangent circle of the HAADF-STEM image. (c) NPs that are accounted for in the T-matrix method. Smaller of the intersecting particles are omitted.

Figure 4 .

 4 Figure 4. Particle size distribution of silver NPs.

Figure 5 .

 5 (b)) between 2 -25 µm is presented in Corresponding daytime radiative cooling without solar absorption is 𝑃 𝑟𝑎𝑑 ′′ -𝑃 𝑎𝑡𝑚 ′′ = 137.5 W m -2 . Solar heat flux absorption, 𝑃 𝑠𝑜𝑙𝑎𝑟 ′′ , due to the 500 µm silica layer and silver reflector is around 47 W m -2 . Therefore, the baseline, the power flux of the whitish radiative cooler in this study, has a 𝑃 𝑐𝑜𝑜𝑙 ′′ = 90.6 W m -2 .

  Figure S4(a) illustrates size dependent and bulk refractive indices whereas Figure S4(b) illustrates extinction cross sections of single particles at different radii that consider size dependent and bulk refractive indices. Results show that peaks are narrower with bulk refractive

  (c)) which has 824 spheroids. Then we calculated the absorption cross-sections of each 824 spheroids separately (without considering neighboring particles) and added the absorption cross-sections up. Error! Reference source not found. indicated that interparticle coupling effect resulted in a red shift as well as widening and damping in the absorption peak as expected[START_REF] Jain | On the Universal Scaling Behavior of the Distance Decay of Plasmon Coupling in Metal Nanoparticle Pairs: A Plasmon Ruler Equation[END_REF].

Figure 7 .

 7 Figure 7. Absorption cross-section of the digitized particle disk (shown in Figure 3(c)) and cumulative individual absorption cross-sections of the spheroids inside the disk.

Figure 8 Figure 8 .

 88 Figure 8 illustrates 6 different geometries are constructed with 6 different cmin. It is worth noting that all of the 6 different geometries have the same top view which is the same as Figure 3(c).

Figure 6 Figure 9 .

 69 Figure 6 illustrates the spectral transmittance for the multilayered case at cmin = 35 nm. The transmittance dip is much narrower and deeper compared to the single layered densely located particles. This behavior is promising for colored radiative cooling applications where high lightness, L, and low solar absorptance Psolar is favored. The case study for surface plasmons was also performed for the radiative cooling structure (Figure 1(b)) to observe the coupling effect on the structure's solar heat absorption and color. D65 illumination and a CIELUV standard representing the outside illumination accounting for the environment is considered while

Table 1 .

 1 The investigated cases in this study.

	Case name	Particle shape	Particle size distribution	Refractive index
	Real sample		Polydisperse	Size dependent
	Spheroid assumption	Oblate spheroid Polydisperse	Size dependent
	Sphere assumption	Sphere	Polydisperse	Size dependent
	Monodisperse assumption	Oblate spheroid Monodisperse	Size dependent
	Bulk ref. index assumption	Oblate spheroid Polydisperse	Bulk data

Table 2 .

 2 The investigated cases in this study and their corresponding root mean square error of transmittance with respect to experimental result.

	Case name	Root mean square error
	Real sample	-
	Spheroid assumption	2.9%
	Sphere assumption	11.3%
	Monodisperse assumption	9.4%
	Bulk ref. index assumption	6.8%

Acknowledgement

We are grateful to Dr. David Lantiat for performing the HAADF-STEM characterizations. R.A.Y is also grateful for financial support from the Marie Skłodowska-Curie Actions Postdoctoral Individual Fellowship Program (project no: 892456)

Supplementary Information

S.1. Calculating the color

Estimation of the color is performed according to the position of eye. For radiative coolers (Figure 1(b)) normal reflected light, R(0,𝜆). For the source, D65 CIE standard illuminant representing a phase of daylight with a correlated color temperature of approximately 6504 K. Spectrum to color conversion is then performed by using the CIELUV color space.

First, using the spectral reflected (or transmitted) intensity R(0,𝜆) ( and scattering (dashed line) cross sections of alumina capped single particle sitting on silica.