Maciej Dunajski
email: m.dunajski@damtp.cam.ac.uk

Prim Plansangkate

Scalar-flat Kähler metrics with conformal Bianchi V symmetry

, and construct all four-dimensional Kähler metrics with vanishing scalar curvature which are invariant under the conformal action of Bianchi V group. The construction is based on the combination of twistor theory and the isomonodromic problem with two double poles. The resulting metrics are non-diagonal in the left-invariant basis and are explicitly given in terms of Bessel functions and their integrals. We also make a connection with the LeBrun ansatz, and characterise the associated solutions of the SU (∞) Toda equation by the existence a non-abelian two-dimensional group of point symmetries.

Introduction

Let (M, g) be a Riemannian four-manifold with anti-self-dual (ASD) Weyl curvature. The metric g is said to have cohomogeneity-one if (M, g) admits an isometry group acting transitively on codimension-one surfaces in M. We shall say that a conformal structure [g] on M has cohomogeneity-one, if there exists a cohomogeneity-one metric g ∈ [g].

The four-dimensional cohomogeneity-one metrics can be classified according to the Bianchi type of the three-dimensional real Lie algebra [START_REF] Bianchi | Sugli spazii a tre dimensioni che ammettono un gruppo continuo di movimenti[END_REF] of the isometry group G. Locally M = R × G, and the problem of finding ASD cohomogeneity-one metrics reduces to solving a system of ODEs, as the group coordinates do not appear in the Weyl tensor. Moreover, the reduction-integrability dogma suggests that the resulting ODEs will be in some sense solvable as the underlying anti-self-duality equations are integrable by twistor transform [START_REF] Penrose | Nonlinear gravitons and curved twistor theory[END_REF][START_REF] Mason | Integrability, selfduality, and twistor theory[END_REF][START_REF] Dunajski | Solitons, Instantons & Twistors[END_REF]. This is indeed the case. If G = SU(2) then the ODEs reduce to Painleve VI, or (if additional assumptions are made about the metric) Painleve III ODEs [START_REF] Tod | Self-dual Einstein metrics from the Painlevé VI[END_REF][START_REF] Hitchin | Twistor spaces, Einstein metrics and isomonodromic deformations[END_REF][START_REF] Tod | Scalar-flat Kähler and hyper-Kähler metrics from Painlevé III[END_REF][START_REF] Tod | Cohomogeneity-one metrics with self-dual Weyl tensor[END_REF][START_REF] Maszczyk | Selfdual Bianchi metrics and the Painleve transcendents[END_REF][START_REF] Dancer | Scalar-flat Kähler metrics with SU(2) symmetry[END_REF]. In [START_REF] Tod | Cohomogeneity-one metrics with self-dual Weyl tensor[END_REF] Tod has analysed the general case of conformal ASD cohomogeneity-one metrics with arbitrary group, diagonal in the basis of left-invariant one forms on the group. He has shown that cohomogeneity-one diagonal Bianchi V metrics are conformally flat. The diagonalisabilty assumption can only be made without loss of generality if the underlying metric is Einstein, so the existence of non-diagonal Bianchi V conformal structures has not been ruled out in Tod's work. In particular he has raised a question whether such structures (if they exist) can admit a Kähler metric in a conformal class. In this paper we shall use the method of isomonodromic deformations to construct all such conformal structures.

Let g be a Lie algebra of Bianchi type V, and let G be a corresponding simply connected Lie group (see Appendix). The left-invariant vector fields on G satisfy

[L 0 , L 1 ] = L 1 , [L 0 , L 2 ] = L 2 , [L 1 , L 2 ] = 0.
(1.1)

A general cohomogeneity-one metric on M = R × G will be have the form

g = dt 2 + h jk (t)λ j λ k , j,k = 0, 1, 2
where λ j are the Maurer-Cartan one-forms on G such that L j λ k = δ j k . We will say that a group action on a complex manifold is holomorphic if Lie derivatives of the complex stucture along any of the generators vanish.

Theorem 1.1 Any cohomogeneity-one Bianchi V ASD conformal structure which admits a Kähler metric such that the group acts holomorphically can be locally represented by a cohomogeneity-one metric g = e 1 e 1 + e 2 e 2 , (1.2)

where the complex one-forms e 1 , e 2 are dual to the vector fields

E 1 = ∂ t - i t (r 0 L 0 + r 1 L 1 + r 2 L 2 ), E 2 = p 0 L 0 + p 1 L 1 + p 2 L 2
and the real functions r(t) = (r 0 (t), r 1 (t), r 2 (t)) and complex functions p(t) = (p 0 (t), p 1 (t), p 2 (t)) satisfy the linear system of ODEs

dr dt = tIm (p 0 p), t dp dt = i(r 0 p -p 0 r), dp 0 dt = 0, dr 0 dt = 0 (1.3) such that det [r, p, p] = 0.
We shall prove this theorem in Section 2, where we shall also show that equations (1.3) reduce to the Bessel equation and quadratures.

All ASD Kähler metrics have vanishing scalar curvature -we shall follow LeBrun [START_REF] Lebrun | Explicit self-dual metrics on CP 2 # • • • #CP 2[END_REF] in calling them scalar-flat Kähler -and conversely all scalarflat Kähler metrics have anti-self-dual Weyl tensor (we choose the natural orientation, where the Kähler two-form is self-dual). In Section 3 we shall construct the Kähler structure in the conformal class of Theorem 1.1. We shall demonstrate that although a G-invariant conformal factor turning the metric (1.2) into a Kähler metric does not exist, there is a simple function on G which does the job. We shall prove Proposition 1.2 Any conformally Bianchi-V scalar-flat Kähler metric can locally be put into a form

g K = Ω 2 g, (1.4)
where g is given by (1.2) and Ω : G → R is a function on the group such that

L 0 (Ω) = Ω, L 1 (Ω) = 0, L 2 (Ω) = 0. (1.5)
The holomorphic [START_REF] Bianchi | Sugli spazii a tre dimensioni che ammettono un gruppo continuo di movimenti[END_REF]0) vector fields are given by E 1 and E 2 . Moreover g K is Ricci flat if and only if it is flat.

The right-invariant vector fields R j on G generate the conformal transformations of the Kähler metric (1.4)

L R 0 g K = 2g K , L R 1 g K = 0, L R 2 g K = 0,
where L denotes the Lie derivative. If we choose coordinates (ρ, x 1 , x 2 ) on G (see Appendix) such that

L 0 = ρ ∂ ∂ρ , L 1 = ρ ∂ ∂x 1 , L 2 = ρ ∂ ∂x 2 , ρ ∈ R + , (x 1 , x 2 ) ∈ R 2
then Ω = ρ. Equations (1.3) imply that r 0 and p 0 are constants. The remaining equations reduce to the Bessel equation possibly with imaginary order. A particularly simple class of solutions characterised by r 0 = 0 is

r = (0, tJ 1 (t), tY 1 (t)), p = (1, iJ 0 (t), iY 0 (t)),
where J α (t) and Y α (t) are Bessel functions of first and second type of order α. This leads to the following example of a Bianchi-V scalar flat Kähler metric

g K = 1 4 (dρ 2 + ρ 2 dt 2 ) + G AB (t)dx A dx B , A,B = 1, 2 (1.6)
where

G(t) = π 2 t 2 16
Y 2 0 + Y 2 1 -J 0 Y 0 -J 1 Y 1 -J 0 Y 0 -J 1 Y 1 J 2 0 + J 2 1 .
All scalar-flat Kähler metrics admitting a Killing vector which also preserves the Kähler form arise from the LeBrun's ansatz [START_REF] Lebrun | Explicit self-dual metrics on CP 2 # • • • #CP 2[END_REF], and are determined by a solution to the SU(∞)-Toda field equation

u xx + u yy + (e u ) zz = 0, where u = u(x, y, z), (1.7)
together with a solution to its linearisation. Kähler metrics arising from Proposition 1.2 admit a two-dimensional group of symmetries (generated by the right translations R 1 and R 2 ) preserving the Kähler from. Thus any linear combination of R 1 and R 2 should lead to a solution of (1.7). In Section 4 we shall characterise the solutions corresponding to the metric (1.6). They admit a two-dimensional non-abelian group of Lie point symmetries and can be found by making an ansatz u = u(y/z).

ASD structures and isomonodromy

A Kähler structure on a four-dimensional real manifold M consists of a pair (g, I) where g is a Riemannian metric and I : T M → T M is a complex structure such that, for any vector fields X, Y , g(X, Y ) = g(IX, IY ) and the two-form ω defined by ω(X, Y ) = g(IX, Y ) is closed. Given an orientation on M the Hodge operator * : Λ 2 → Λ 2 satisfies * 2 = Id and gives a decomposition Λ 2 = Λ 2 + ⊕ Λ 2 -of two-forms into self-dual (SD) and anti-self-dual (ASD) components. The two-form ω induces a natural orientation on M given by the volume form ω ∧ ω. With respect to this orientation ω is self-dual. It is well known [START_REF] Pontecorvo | On twistor spaces of anti-self-dual Hermitian surfaces[END_REF][START_REF] Dunajski | Solitons, Instantons & Twistors[END_REF] that if the scalar curvature of a Kähler metric vanishes, then the Weyl tensor of the underlying conformal structure is ASD. Conversely ASD Kähler metrics are scalar-flat.

A convenient way to express the ASD condition on a conformal structure is summarised in the following theorem. The theorem below is originally due to Penrose [START_REF] Penrose | Nonlinear gravitons and curved twistor theory[END_REF], but taken in this form from [START_REF] Mason | Integrability, selfduality, and twistor theory[END_REF][START_REF] Dunajski | Solitons, Instantons & Twistors[END_REF].

Theorem 2.1 Let E 1 , E 2 be two complex vector fields in T M ⊗ C and let e 1 , e 2 be the corresponding dual one-forms. The conformal structure defined by g = e 1 e 1 + e 2 e 2 is ASD if and only if there exists functions f 0 , f 1 on M × CP 1 holomorphic in λ ∈ CP 1 such that the distribution

l = E 1 -λE 2 + f 0 ∂ ∂λ , m = -E 2 -λE 1 + f 1 ∂ ∂λ (2.8)
is Frobenius integrable, that is, [l, m] = 0 modulo l and m.

A general ASD conformal structure [g] does not admit a Kähler metric in its conformal class. The existence of such metric is characterised by vanishing of higher order conformal invariants of [g] [START_REF] Dunajski | Four Dimensional Metrics Conformal to Kähler[END_REF]. In what follows we shall use a simpler twistor characterisation of ASD metrics which are conformal to Kähler.

Twistors and divisors

Let us complexify (M, g) and regard M as a holomorphic four-manifold with a holomorphic metric g

g = V 00 V 11 -V 01 V 10 , (2.9)
where V AA , A, A = 0, 1 is the null tetrad of one forms written in the twocomponent spinor notation. The reality conditions can be imposed to recover the Riemannian metric by setting

V 00 = e 1 , V 10 = -e 2 , V 01 = e 2 , V 11 = e 1 .
Let V AA be the corresponding tetrad of vector fields so that the integrable distribution of Theorem 2.1 is

l = V 00 -λV 01 + f 0 ∂ ∂λ , m = V 10 -λV 11 + f 1 ∂ ∂λ .
(2.10)

A twistor space of (M, g) is the space of two-dimensional totally null surfaces spanned by V 00 -λV 01 , V 10 -λV 11 in M. It is a three dimensional complex manifold Z which arises as a quotient of M × CP 1 by l, m. The points in M correspond to rational curves (called twistor curves) in Z with normal bundle O(1) ⊕ O(1), where O(n) is a line bundle over CP 1 with Chern class n. The holomorphic canonical line bundle κ of Z restricted to any of these curves is isomorphic to O(-4). To reconstruct a real four-manifold M, the twistor curves must be invariant under an anti-holomorphic involution τ on Z which restricts to an antipodal map on each curve. A theorem of Pontecorvo [START_REF] Pontecorvo | On twistor spaces of anti-self-dual Hermitian surfaces[END_REF] states that an ASD conformal structure (M, [g]) admits a Kähler metric if and only if there exists a section D of the bundle O(2) = κ -1/2 over the twistor space Z of (M, [g]), with exactly two distinct zeros on each twistor line. The section must be τ invariant, so the its zeroes lie on the antipodal points on each twistor curve. Now consider the group G acting on an ASD conformally Kähler manifold M by holomorphic conformal isometries with generically three-dimensional orbits. This gives rise to a complexified group action of G C on the twistor space Z. Let R j , j = 0, 1, 2 be holomorphic vector fields on Z corresponding to the right-invariant conformal Killing vectors R j on M. The subset of Z where R j are linearly dependent is given by the zero set of s = vol Z ( R 0 , R 1 , R 2 ). As the canonical bundle κ = O(-4), the divisor s = 0 defines a quartic and vanishes at four points on each twistor line. Hitchin ([8], Proposition 3) showed 1 that s is not identically zero if g is not Ricci-flat and that the divisor when s vanishes is equal to the Pontecorvo's divisor D in the case when [g] is cohomogeneity one and contains a Kähler class. In the double fibration picture

M ←-M × CP 1 -→ Z the section s pulls back to s = (dλ ∧ ν)(l, m, R 0 , R 1 , R 2 ), (2.11)
where ν = V 01 ∧ V 10 ∧ V 11 ∧ V 00 is the volume form on M and R j are the lifts of the three generators of

G to M × CP 1 such that [l, R j ] = 0, [m, R j ] = 0 modulo l, m,.
Thus in the proof of Theorem 1.1 we will require that this quartic has two distinct zeros of order two. This will guarantee the existence of a Kähler metric in the cohomogeneity-one conformal class.

Proof of Theorem 1.1. We will work in the complexified category and impose the reality conditions at the end. Let G C be a three-complex dimensional Lie group (we will eventually take G C to be a complexification of the Bianchi V group, but the first part of the proof does not depend on the choice of G C ). We assume that the orbits of G C are three-dimensional, and the metric g on C × G C is invariant under the left translations of G C on itself. One can write the null tetrad V AA in terms of the vector field ∂ t and three linearly independent vector fields P, Q, R tangent to G C (in the complexified setting the fields Q and P are independent. Once the reality conditions are imposed at the end of the proof, we shall set Q = -P ) which are t-dependent and invariant under the left translations, as

V 00 = ∂ t + i R t , V 11 = ∂ t -i R t , V 01 = P, V 10 = Q. (2.12)
Now, let R 0 , R 1 , R 2 be the right-invariant vector fields on G corresponding to three generators of the left translations. Since R j are independent of t, one has

[l, R j ] = -R j (f 0 )∂ λ , [m, R j ] = -R j (f 1 )∂ λ .
A direct calculation shows that there is no lift of R j of the form

R j + Q j ∂ λ for some function Q j such that [l, R j + Q j ∂ λ ] = 0, [m, R j + Q j ∂ λ ] = 0 modulo l, m. Hence, we conclude that [l, R j ], [m, R j ]
, are identically zero. This implies that f 0 and f 1 are constant on G, and hence they are functions of λ and t only.

We claim that the quartic s as defined in (2.11) is proportional to λf 0 +f 1 , with the proportionality factor given by a function h on M. Indeed, using the fact that V 01 , V 00 -V 11 , V 10 are linearly independent and invariant under the left translations of G C , the right-invariant vector fields R j can be written in the basis of

P = V 01 , Q = V 10 , R = it 2 (V 11 -V 00 ).
Thus, performing all the contractions we are left with s ∝ (det H) (λf 0 + f 1 ), where H is the matrix of coefficients of R 0 , R 1 , R 2 written in the basis of P, Q, R.

As det H does not depend on λ, the quartic s has two distinct zeros of order two if and only if λf 0 + f 1 has two distinct zeros of order two (for the moment we rule out the case where f 0 and f 1 are both zero which corresponds to a hyper-Kähler metric [START_REF] Dunajski | Solitons, Instantons & Twistors[END_REF]). We shall assume that this is the case so that the ASD conformal structure admits a Kähler metric. It is now possible to use Möbius transformation to put the two zeros at 0 and ∞. The Möbius transformation in λ corresponds to a change of null tetrad by a right rotation V AA → V AA r A B , where r is an SL(2, C)-valued function. Since the coefficients of the quartic λf 0 + f 1 are functions of t only, the required r A B will only depend on t. Thus, the rotated tetrad is still G C -invariant.

With the zeros at 0 and ∞, λf 0 + f 1 is of the form a(t)λ 2 . Let us first assume that a does not vanish identically. One still has a Möbius degree of freedom that preserves (0, ∞), that is, the multiplication of λ by a function of t. Let us use this freedom to set a(t) = 2/t. The current tetrad is some right rotation of the original one. It is possible to use another freedom: a left rotation

V AA → l A B V AA , l ∈ SL(2, C) to keep V 00 -V 11 , V 01 , V 10 tangent to G C .
The right rotation does not change the quartic λf 0 + f 1 , and we now have f 0 , f 1 of the form

f 0 = b(t)λ 2 + c(t)λ + d(t), f 1 = -b(t)λ 3 + 2 t -c(t) λ 2 -d(t)λ (2.13)
for some functions b(t), c(t), d(t).

Now, consider a pair of linear combinations of l and m (2.10)

L = λl + m λf 0 + f 1 = ∂ ∂λ + 2λiRt -1 -λ 2 P + Q λf 0 + f 1 (2.14) M = f 1 l -f 0 m λf 0 + f 1 = ∂ ∂t + (f 1 -λf 0 )iRt -1 -λf 1 P -f 0 Q λf 0 + f 1 .
Since the conformal class is ASD, Theorem The compatibility conditions [L, M ] = 0 are then given by

tP t -i[R, P ] + (tc(t) -1)P = 0, 2iR t -t[P, Q] = 0, (2.15) tQ t + i[R, Q] -(tc(t) -1)Q = 0.
This shows that a cohomogeneity-one metric (1.2), in the tetrad (2.12) is ASD if the vector fields P, Q, R satisfy the system (2.15), where c(t) is defined in (2.13). Now, let

R = R, P = h(t)P, Q = h -1 (t)Q, where h(t) = e R (c(t)-1 t )dt . (2.16)
The system (2.15) implies that the vector fields P , Q, R satisfy

tP t -i[R, P ] = 0, 2iR t -t[P, Q] = 0, (2.17) tQ t + i[R, Q] = 0,
where we have dropped the hat from the rescaled vector fields. Moreover, the tetrad (2.12) constructed from a solution (P, Q, R) of (2.17) gives the same metric (1.2) as the tetrad determined from (h -1 (t)P, h(t)Q, R) which satisfy (2.15) with c(t) = ht h + 1 t . Neither the metric nor the equations (2.17) depend on h, so we can set it equal to 1. The resulting Lax pair (2.14) is

L = ∂ ∂λ + (tQ + 2iλR -λ 2 tP ) 2λ 2 , M = ∂ ∂t - (λQ + λ 3 P ) 2λ 2 , ( 2.18)
where we again dropped the hat from the rescaled vector fields.

We conclude that any cohomogeneity-one metric (1.2), which belongs to an ASD conformal structure admitting a quartic s defined in (2.11) with two distinct zeros of order two, can be written in terms of a null tetrad (2.12), where the vector fields P, Q, R satisfy the system (2.17).

The three vector fields P, Q, R can be written in the basis of left-invariant vector fields L 0 , L 1 , L 2 satisfying (1.1) as

P = p 0 (t) L 0 + p 1 (t) L 1 + p 2 (t) L 2 Q = q 0 (t) L 0 + q 1 (t) L 1 + q 2 (t) L 2 , (2.19) R = r 0 (t) L 0 + r 1 (t) L 1 + r 2 (t) L 2 ,
for some functions p j (t), q j (t), r j (t). Then using the commutation relation (1.1) the system (2.17) implies that p 0 , q 0 , r 0 are constant and that t(p j ) t = i(r 0 p j -p 0 r j ), t(q j ) t = i(-r 0 q j + q 0 r j ), 2i(r j ) t = t(p 0 q j -q 0 p j ).

(2.20) The reality conditions corresponding to Riemannian metrics come down to choosing R real (so that (r 0 , r 1 , r 2 ) are all real functions) and Q = -P so that q j = -p j ). Equations (2.20) then give the linear system (1.3).

Let us now return to the case a = 0 which corresponds to f 0 and f 1 vanishing in the distribution (2.10). The resulting conformal structure must then be hyper-Hermitian [START_REF] Dunajski | Solitons, Instantons & Twistors[END_REF] and (as the vector fields in (2.10) are volume preserving) it is actually conformal to hyper-Kähler. The integrability [l, m] = 0 modulo l, m implies the system of Nahm equations

P t -i[R/t, P ] = 0, 2i(R/t) t -[P, Q] = 0, Q t + i[R/t, Q] = 0. (2.21)
Imposing the reality condition that R is real and Q = -P , (2.21) become

d(rt -1 ) dt = Im (p 0 p), t dp dt = i(r 0 p -p 0 r), dp 0 dt = 0, d(r 0 t -1 ) dt = 0.
The general solutions for r and p are given in terms of trigonometric and exponential functions depending on the constants p 0 and r 0 /t, and the resulting metric is conformally flat (in the proof of Proposition 1.2 we shall find the conformal factor which makes it flat).

2

Our derivation of (2.17) and (2.18) did not depend on the choice of the isometry group. The system (2.17) describes the general isomonodromic deformation equations with two double poles. The corresponding Lax pair (2.18) with P, Q, R given by 2 × 2 matrices was shown by Jimbo and Miwa [START_REF] Miwa | Monodromy preserving deformations of linear ordinary differential equations with rational coefficients[END_REF] to give rise to the Painlevé III equation. The same Lax pair also arises as the reduced Lax pair of the ASDYM equation, by the Painlevé III group. It is shown to be the isomonodromic Lax pair for the Painlevé III equation when the gauge group of the ASDYM connection is SL(2, C) [START_REF] Mason | Self-duality and the Painlevé transcendents[END_REF], or (with certain algebraic constraints on normal forms) SL(3, C) [START_REF] Dunajski | Strominger-Yau-Zaslow geometry, Affine Spheres and Painleve III[END_REF].

Bessel equation

For generic values of the constants (r 0 , p 0 ), the solutions to (1.3) are determined by two linearly independent solutions of the Bessel equation2 .

If p 0 = 0 then r is a constant vector and equations for p 1 , p 2 can be easily integrated. The resulting metric is conformally flat.

If p 0 = 0 we rescale the metric by a constant |p 0 | 2 and redefine (r, p 1 , p 2 , t) to set p 0 = 1. Differentiating the first set of equations in (1.3) twice, and using the second set of equations shows that the real functions f 1 = 2t -1 (r 1 ) t and f 2 = 2t -1 (r 2 ) t satisfy a pair of Bessel equations

t 2 d 2 f k dt 2 + t df k dt + (t 2 + r 2 0 )f k = 0, k = 1, 2. ( 2

.22)

If r 0 = 0 the general solution of the Bessel equation (2.22) is given in terms of Bessel functions of pure imaginary order ir 0

f 1 = c 1 J ir 0 (t) + c 2 Y ir 0 (t), f 2 = c 3 J ir 0 (t) + c 4 Y ir 0 (t).
The constant complex coefficients c 1 , c 2 , c 3 , c 4 can be chosen so that the functions f 1 and f 2 are real, see for example [START_REF] Dunster | Bessel functions of purely imaginary order, with an application to second-order linear differential equations having a large parameter[END_REF]. Given functionally independent f 1 and f 2 we find r 1 , r 2 , p 1 , p 2 by integrations and algebraic manipulations. The case r 0 = 0 is special. The linear system (1.3) now reduces to a pair of ODEs

t 2 d 2 r k dt 2 -t dr k dt + t 2 r k = 0, k = 1, 2
whose solutions are given by

r 1 = c 1 tJ 1 (t) + c 2 tY 1 (t), r 2 = c 3 tJ 1 (t) + c 4 tY 1 (t),
where J α are Y α are Bessel functions of the first and second kind respectively of order α, and c 1 , . . . , c 4 are real constants of integrations. In this case p 1 , p 2 must be purely imaginary and the recursion relations

∂ t J 0 = -J 1 , ∂ t Y 0 = -Y 1 , ∂ t (tJ 1 ) = tJ 0 , ∂ t (tY 1 ) = tY 0 (2.23) imply that p 1 = i(c 1 J 0 (t) + c 2 Y 0 (t)), p 2 = i(c 3 J 0 (t) + c 4 Y 0 (t)). (2.24)
Performing a linear transformation of the vector fields L 1 and L 2 we can always set

c 1 = 1, c 2 = 0, c 3 = 0, c 4 = 1 which yields r = (0, tJ 1 (t), tY 1 (t)), p = (1, iJ 0 (t), iY 0 (t)).
To write down the metric we invert the vector fields E 1 , E 2 from Theorem 1.1 and use the fact that 2t(Y 0 J 1 -J 0 Y 1 ) = 4/π. This yields

g = 1 4 (dt 2 + (λ 0 ) 2 ) + G AB (t)λ A λ B , where A, B = 1, 2
G(t) = π 2 t 2 16
Y 2 0 + Y 2 1 -J 0 Y 0 -J 1 Y 1 -J 0 Y 0 -J 1 Y 1 J 2 0 + J 2 1
, where λ 0 , λ 1 , λ 2 are the left-invariant one forms on G which satisfy (A1).

Kähler structure

Proof of Proposition 1.2. The Kähler structure on M can be read off from the divisor (2.11). In the proof of Theorem (1.1) we have moved the double zeros of s to 0 and ∞, which in spinor notation means that

ω A B = o (A ι B ) and s ≈ (ω A B π A π B ) 2
where π A are the homogeneous coordinates on CP 1 such that λ = -π 0 /π 1 , and the spinor basis is o A = (0, 1), ι A = (-1, 0). Thus, in the null tetrad (2.12), the Kähler form is proportional to

ω = i 2 ε AB ω A B V AA ∧ V BB = i 2 (e 1 ∧ e 1 + e 2 ∧ e 2 ),
and the space T 1,0 (M) of holomorphic vector fields on M is spanned by V 11 , V 01 (equivalently by the vectors E 1 and E 2 in Theorem 1.1). The Frobenius integrability conditions

[T 1,0 , T 1,0 ] ⊂ T 1,0
guaranteeing the vanishing of Nijenhuis torsion follows from the construction, but we can also verify it directly as

[E 1 , E 2 ] = dp 1 dt - i t r 0 p 1 + i t r 1 p 0 L 1 + dp 2 dt - i t r 0 p 2 + i t r 2 p 0 L 2 = 0,
where we have used the constancy of (r 0 , p 0 ) and equations (1.3).

To determine the conformal factor we look for a function Ω : M → R such that d(Ω 2 ω) = 0. Once this has been found the Kähler metric g K and the associated two-form ω will be given by

g K = Ω 2 g, ω = iΩ 2 2 (e 1 ∧ e 1 + e 2 ∧ e 2 ).
It can be verified by explicit calculation that there is no G-invariant Ω (i.e. there is no conformal factor which depends only on t). To demonstrate that Ω : G → R such that (1.5) holds gives the correct conformal factor, consider 2ΩdΩ ∧ ω + Ω 2 dω = 0.

(3.25)

Since ∂ t Ω = 0, one can write dΩ in the basis of left-invariant one-forms

dΩ = L 0 (Ω)λ 0 + L 1 (Ω)λ 1 + L 2 (Ω)λ 2 .
The three-form dω can be calculated using the expressions for the dual oneforms of E 1 and E 2 in Theorem 1.1 and the Maurer-Cartan's structure equation (A1). Finally, the system (1.3) is used to simplify the LHS of (3.25), and one finds that (3.25) is satisfied if and only if (1.5) holds. This conformal factor is in fact unique -it could have also been read off from the divisor as it is proportional to a power of ω A B ω A B [START_REF] Dunajski | Four Dimensional Metrics Conformal to Kähler[END_REF].

In the conformally flat case where s in (2.11) is identically zero, which corresponds to the system (2.21), the same conformal factor makes g flat. Moreover we verify by explicit calculation that g K given by (1.4) where s = 0 is Ricci flat if and only if it is flat, which proves the last part the Proposition.

2

The conformal Killing vectors generating the group action on g K are given by the right-invariant vector fields on G. If the coordinates are chosen for the group, these vectors are given by (A3). Example. The simplest explicit example of the scalar-flat Kähler metric corresponds to r 0 = 0 in Theorem 1.1 and is given by (1.6). The determinant of the metric (1.6) given by

det g K = π 2 ρ 2 t 2 1024 .
Since by definition ρ = 0, g K may be degenerate only at t = 0 or t = ∞.

The Ricci scalar R is identically zero because the Kähler metric is ASD. The remaining curvature invariants are

R abcd R abcd = 256 ρ 4 t 2 , W abcd W abcd = 128 ρ 4 t 2
which indicates that t = 0 is a singularity. Rescaling the metric by a conformal factor ρ -2 (tf (t)) -1 gives W abcd W abcd = 128f 2 , which needs to be regular if the conformal class contains a complete metric, but the regularity of the conformal factor requires that tf (t) is also regular and non-zero. Thus the norm of the Weyl tensor blows up at 0. The asymptotic behaviour of g K for large t is

g K = 1 4 dρ 2 + ρ 2 dt 2 + πt 8 (d(x 1 ) 2 + d(x 2 ) 2 ).

SU (∞) Toda equation

LeBrun [START_REF] Lebrun | Explicit self-dual metrics on CP 2 # • • • #CP 2[END_REF] has shown that any Kähler metric g K with symmetry preserving the Kähler form admits a local coordinate system {τ, x, y, z} such that

g K = W h + 1 W (dτ + θ) 2 , ω = W e u dx ∧ dy + dz ∧ (dτ + θ), (4.26)
where h = e u (dx 2 + dy 2 ) + dz 2 . (4.27)

Here τ is a coordinate along the orbits of the Killing vector K = ∂ τ , {x, y, z} are coordinates on the space of orbits, and (u, W ) and θ are functions and a one-form on the space of orbits such that u satisfies the SU(∞) Toda equation (1.7), W satisfies the so-called monopole equation

W xx + W yy + (W e u ) zz = 0 (4.28)
and θ is determined by W together with the condition dω = 0. The ansatz (4.26) can be understood as follows. Given that K = ∂ τ is a Killing vector, the metric necessarily takes the form (4.26), where 1 W = g K (K, K) and h is a metric on the three-dimensional space of orbits. Now, since the Kähler form ω Lie derives along K, we have

K ω = dz (4.29)
for some function z on the space of orbits of K. The isothermal coordinates x, y on the orthogonal complement of the space spanned by K and I(K) (where I is the complex structure) can be used together with z to parametrise the space of orbits. The metric then takes the form (4.26) with h given by (4.27) for some u = u(x, y, z). The integrability of the complex structure and the closure of the ω imply (4.28). The scalar-flat condition gives (1.7).

Bessel solutions to SU(∞) Toda equation

The scalar-flat Kähler metric (1.6) has two Killing symmetries ∂/∂x 1 and ∂/∂x 2 preserving the Kähler form. We can follow the algorithm described above and find the solution u of the SU(∞) Toda equation and the associated monopole W corresponding to a linear combination

K = c 1 ∂/∂x 1 + c 2 ∂/∂x 2 .
We shall set c 2 = 0 for simplicity. Set τ = x 1 so that the vector field K = ∂ τ and dx 1 = dτ . Then the metric (1.6) takes the form (4.26), where

h = 1 4W ρ 2 dt 2 + dρ 2 + t 2 π 2 64 d(x 2 ) 2 , W = 16 t 2 π 2 (Y 2 0 + Y 2 1 )
, and θ = -

(J 0 Y 0 + J 1 Y 1 ) Y 2 0 + Y 2 1 d(x 2
τ = x 1 , x = - πx 2 8 , y = - πρY 0 8 , z = πρtY 1 8 . ( 4

.31)

The solution to SU(∞) Toda equation is now implicitly given by e u = t 2 . (4.32) Formulae (4.31) imply that u = u(v), where v = z/y. Thus u is constant on the plane y v(u) -z = 0 (compare [START_REF] Tod | Scalar-flat Kähler and hyper-Kähler metrics from Painlevé III[END_REF] where solutions constant on quadrics were constructed) and is invariant under a two-dimensional group of Lie point symmetries generated by vector fields ∂/∂x, x∂ x + y∂ y + z∂ z (4.33)

We shall now show that the existence of these symmetries uniquely characterises (4.32). Any solution u of (1.7) which is invariant under symmetries generated by the vector fields (4.33) is a function u = u(v). Then (1.7) becomes an ODE (v 2 + e u )u vv + 2vu v + u 2 v e u = 0. Equivalently, interchanging the dependent and independent variables we have

(v 2 + e u ) 2 ∂ u v u v 2 + e u = 0, which integrates to v u = c(v 2 + e u ) ( 4 . 3 4 )
for some constant c. Now we shall argue that this constant c can always be set to -1/2, provided that we consider solutions to the SU(∞) Toda equation as equivalent if they determine the conformally equivalent metrics (4.27). First, note that a transformation

u -→ u + 2β, z -→ ±e β z, (4.35)
where β is a constant is a symmetry of (1.7) which rescalles the three-metric by a constant factor. Now, under (4.35), the variable v transforms as v -→ ±e β v. This can be used to set c = -1/2. It was shown in [START_REF] Tod | The SU(∞)-Toda field equation and special fourdimensional metrics[END_REF] that any ASD Einstein metric with symmetry and non-zero scalar curvature can be written as

R 0 = τ ∂ τ + x∂ x + y∂ y + z∂ z , R 1 = ∂ τ , R 2 = ∂ x . ( 4
g E = W z 2 e u (dx 2 + dy 2 ) + dz 2 + 1 W z 2 (dτ + θ) 2 , (4.40)
where W = const (zu z -2), and u is a solution to the SU(∞) Toda equation. Take u given by (4.32) and the Einstein monopole W with const= 1/2. Then

W = Y 0 Y 1 t(Y 2 0 + Y 2 1 ) -1, θ = 1 2 Y 2 0 -Y 2 1 Y 2 0 + Y 2 1 + ln πρ 8 dx.
The resulting Einstein metric (4.40) has negative scalar curvature and is non-conformally flat. Note that the metric is not conformal to a Bianchi V metric. It only admits a two-dimensional group of symmetries.

The right-invariant one forms and vector fields can be found analogously from dgg -1 . The algebra of isometries of metrics from Theorem 1.1 is spanned by the right-invariant vector fields, which in our coordinate basis are given by

R 0 = ρ ∂ ∂ρ + x 1 ∂ ∂x 1 + x 2 ∂ ∂x 2 , R 1 = ∂ ∂x 1 , R 2 = ∂ ∂x 2 . ( A3
)
The left and right invariant vector fields satisfy the commutation relations

[L 0 , L 1 ] = L 1 , [L 0 , L 2 ] = L 2 , [L 1 , L 2 ] = 0, [R 0 , R 1 ] = -R 1 , [R 0 , R 2 ] = -R 2 , [R 1 , R 2 ] = 0, [R j , L k ] = 0.

  2.1 means [l, m] = (. . . )l+(. . . )m. This in turn implies that [L, M ] = 0, modulo L and M. However, one sees that [L, M ] does not contain ∂ λ or ∂ t , thus [L, M ] must be identically zero. It turns out that [L, M ] = 0 implies that b(t) = 0 = d(t).

  The equation (4.34) with c = -1/2 is equivalent to the Bessel equation of order 0t 2 Y tt + tY t + t 2 Y = 0. (4.36)To see it use t 2 = e u and v = 2∂ u ln Y . The four-dimensional metric (4.26) resulting from the solution u = u(z/y) together with the monopole (4.30) admits a conformal action of the Bianchi V group generated by

Example: ASD Einstein metric with symmetry.

  .37)We shall now establish that given a solution u = u(z/y) to the SU (∞) Toda equation, the function W given by (4.30) is (up to gauge transformation) the only solution to the linearised SU(∞) Toda equation (4.28) such that the resulting metric (4.26) in four dimensions defines a cohomogeneity-one Bianchi V conformal class. To show it, observe that the metric h in (4.27) is invariant under R 1 , R 2 and conformally invariant h → c 2 h under the one parameter group of transformations generated by R 0 . This implies that the function W must be invariant under the Bianchi V group G and the one-form θ is invariant under the translations generated by R 1 , R 2 , but transforms as θ → c θ under the action generated by R 0 . Thus the metric (4.26) is conformally invariant under G if and only if the function W and the components of the one-form θ are functions of t only. Therefore the PDE (4.28) for W becomes a second order ODEW tt v t (v 2 + t 2 ) + 2W t (v t (2t + vv t ) -v tt (v 2 + t 2 )) + 2W (v t -tv tt ) = 0.The invariance properties of θ lead to one further constraintW t (t 2 + v 2 ) + 2tW = 0.(4.38)The relation (4.31) yields v = -t Y 1 Y 0 and using (2.23) we find the general solution W of (4.38) Changing the proportionality constant k in W amounts to changing the coordinate τ -→ τ /k in (4.26) and rescaling the metric by the constant k. Hence we can always set k = 64/π 2 , which gives the monopole W in (4.31). Here we shall present an example of a metric (4.26) which is obtained from a monopole W different from (4.39).

	W =	k t 2 (Y 2 0 + Y 2 1 )	, k = const.	(4.39)

His result was derived for G C = SL(2, C) but it remains valid for any three-dimensional group acting transitively on M , as the corresponding action of G C preserves the points where Pontecorvo's divisor vanishes as long as the group action is holomorphic. Moreover we are allowed to work in a conformal class of g K , as the twistor equation underlying the Pontecorvo's divisor is conformally invariant[START_REF] Pontecorvo | On twistor spaces of anti-self-dual Hermitian surfaces[END_REF][START_REF] Dunajski | Four Dimensional Metrics Conformal to Kähler[END_REF].

The relation with the Bessel equation is already expected from the result of[START_REF] Maszczyk | The classification of self-dual Bianchi metrics[END_REF], who classified all reductions of anti-self-dual Yang Mills equations leading (by switch map) to cohomogeneity-one ASD conformal structures without however determining a Kähler class.

For a Killing vector given by a general linear combination of ∂/∂x 1 and ∂/∂x 2 , a linear combination of Bessel functions appears in the final formula, and in particularz y = -t c 2 J 1 -c 1 Y 1 c 2 J 0 -c 1 Y 0 .

Acknowledgements. We are grateful to Philip Boalch, Andrew Dancer and Paul Tod for helpful discussions.

Appendix

The real three-dimensional Lie algebra of Bianchi type V is defined by commutation relations

We can choose its representation by 3 × 3 matrices

The corresponding Lie group G is the multiplicative group of real matrices of the form

The left-invariant one-forms {λ j , j = 0, 1, 2} corresponding to a basis {X j } of a Lie algebra are given by g -1 dg = λ j X j , where g ∈ G.

Hence

and

The left invariant vector fields defined by L j λ k = δ j k are found to be


P M Hansen

J Alvarez-Muñiz

R A Vázquez

A comprehensive study of shower to shower fluctuations

Keywords: Cosmic rays, Extensive air showers, Ground detector, Simulation, Muon component, Electromagnetic component

By means of Monte Carlo simulations of extensive air showers (EAS), we have performed a comprehensive study of the shower to shower fluctuations affecting the longitudinal and lateral development of EAS. We split the fluctuations into physical fluctuations and those induced by the thinning procedure customarily applied to simulate showers at EeV energies and above. We study the influence of thinning on the calculation of the shower to shower fluctuations in the simulations. For thinning levels larger than R thin = 10 -5 -10 -6 , the determination of the shower to shower fluctuations is hampered by the artificial fluctuations induced by the thinning procedure. However, we show that shower to shower fluctuations can still be approximately estimated, and we provide expressions to calculate them. The influence of fluctuations of the depth of first interaction on the determination of shower to shower fluctuations is also addressed.

Introduction

Extensive air showers (EAS) have been studied over the last 70 years [1]. They result from the interaction in the atmosphere of high-energy protons and nuclei arriving from space. The product of these collisions are a set of secondary particles carrying a fraction of the primary energy. These secondaries move through the atmosphere and interact again generating new secondaries. The process continues, increasing the number of secondary particles, until their energies are too low to contribute to the generation of new particles. Particles reaching ground are sampled with arrays of detectors, and their properties are used to infer the properties of the primary initiating the shower. Measurements of the electron and muon density, of the arrival time of the particles at ground, and of the depth at which the shower has the maximum number of particles (X max ), give information on the arrival direction, primary energy, and on the mass of the primaries [1].

The complexity of the cascade phenomena, and the poor knowledge of the hadronic interactions at very high energy [2], make the experimental determination of the properties of the primaries very difficult. Moreover, primary particles with the same energy, mass and direction produce secondary particles with parameters that vary from shower to shower. This feature is called "shower to shower fluctuations". An understanding of the shower to shower fluctuations will help to improve the interpretation of cosmic-ray data.

The calculation of shower to shower fluctuations can in principle be addressed with Monte Carlo simulations of extensive air showers. However, the number of particles that are produced in an air shower at ultra high energy (above ∼ 10 18 eV) is so large (∼ 10 10 ), that it is almost impossible to follow the propagation to ground level of all the secondaries in the Monte Carlo in a reasonable amount of time, or even to store the large amount of information produced. For this reason, a statistical sampling procedure called "thinning" [START_REF] Hillas | Proc of the Paris Workshop on Cascade simulations[END_REF] is used in the simulations. Thinning algorithms typically consist on propagating only a small, representative fraction of the total number of particles in the shower, assigning statistical weights to the sampled particles to compensate for the rejected ones. However, thinning algorithms introduce artificial fluctuations in the simulated showers, hampering the determination of the intrinsic, physical shower to shower fluctuations with Monte Carlo simulations. For this reason the study of fluctuations using Monte Carlo simulations is quite difficult and uncertain. This is of utmost importance in cosmic-ray physics, since an incorrect assumption on the shower to shower fluctuations can lead to systematic errors on the determination of the parameters of the primary particles.

In this work we address the problem of determining the true, physical shower to shower fluctuations in Monte Carlo simulations, and quantify the effect of thinning on their determination. We give expressions that allow the estimation of physical fluctuations from Monte Carlo simulations, even in the case of relatively strongly thinned showers. Other recent approaches which study the effect of artifical fluctuations due to the thinning procedure are given in [4,5] The paper is organized as follows: In Section 2 we describe the simulations performed in this work, and the thinning algorithm adopted. In Section 3 we identify the different sources of fluctuations in shower simulations. In Section 4 we perform a comprehensive study of the fluctuations in the longitudinal and lateral shower development, and give expressions that allow to separate physical shower to shower fluctuations from the artificial fluctuations induced by the thinning procedure, and we quantify for which thinning levels these expressions can be applied. In Section 5 we apply our results to study the behaviour of the physical fluctuations at ground. In Section 6, we quantify the influence on the shower to shower fluctuations of the fluctuations of the depth of first interaction of the primary initiating the shower. Finally, we summarize our conclusions in section 7. In the Appendix we give an explicit mathematical derivation of the expressions presented in Section 4.

The simulations

In this work we have used the air shower simulation program, AIRES [START_REF] Sciutto | Proc. 27th ICRC[END_REF][START_REF] Sciutto | AIRES User's Manual and Reference Guide[END_REF], along with the hadronic model QGSJET01 [START_REF] Kalmykov | [END_REF] to simulate proton and ironinduced showers with primary energy 10 19 eV. As explained above, due to the large number of particles that are created in the simulation, AIRES includes a statistical sampling algorithm, that consists on propagating a small, representative fraction of the total number of particles, assigning a statistical weight w to the sampled particles to compensate for the rejected ones. The weight is adjusted in such a way that both the total energy and the average number of particles is guaranteed to be conserved.

Before the simulation starts, the user indicates, as an input to AIRES, the relative thinning level R thin . The thinning energy E thin -the energy below which the thinning process starts -is defined as E thin =R thin × E p where, E p is the primary energy. For ultra high energy cosmic ray shower simulations convenient values for the relative thinning are R thin = 10 -5 -10 -9 , but the actual choice depends on the purpose of the simulation. The thinning level affects both the simulation CPU time and the size of the output produced in the simulation, both typically behaving linearly with R -1 thin . If we increase R thin by a factor of 10 the simulation speeds up by a similar factor, and the output is reduced accordingly, but the price to pay is an enhancement of the artifical fluctuations in the simulated showers as discussed below.

We describe here the thinning algorithm implemented in the AIRES code [START_REF] Sciutto | AIRES User's Manual and Reference Guide[END_REF], originally due to Hillas [START_REF] Hillas | Proc of the Paris Workshop on Cascade simulations[END_REF]. At the beginning of the simulation, the primary particle is assigned a weight w = 1. Then the primary is propagated and interacts in the atmosphere producing n secondary particles. Before incorporating any secondary particle in the simulation, the energy of the primary E p which has generated that secondary is compared to E thin . If E p > E thin , then all the secondaries with energy greater or equal than E thin are kept, and their weight is equal to the weight of the primary particle. Secondaries with energy less than E thin are kept with a probability p i = E i /E thin (E i is the energy of i th secondary), and their weight is adjusted so that w i = (1/p i ) × w, with w being the weight of the mother particle producing that secondary. On the other hand if E p < E thin , it means that the particle came from a previous thinning operation. Then, one and only one of all the produced secondaries -say the j th -is kept, with probability p j = E j / n i=1 E i . Again, the weight of this particle is increased by a factor w j = (1/p j ) × w.

To avoid confusion between particles and weights, we identify an entry with a particle explicitely followed in the simulation which is associated a weight w. Hence, an entry represents w particles. It is important to stress that once the thinning energy is reached, the number of entries N e is no longer increased in the shower processes (only one secondary particle is followed in each interaction), while the number of particles N does however increase, since the weight of each entry typically increases in the showering process. When evaluating a physical observable, each entry must be weighted with its corresponding statistical weight.

In AIRES, the thinning algorithm is complemented with an "extended thinning algorithm" [START_REF] Sciutto | AIRES User's Manual and Reference Guide[END_REF], designed to ensure that all the statistical weights are always smaller than a certain positive number (other algorithms based on this same idea are possible, see for instance [9]). To ensure this, an external parameter called statistical weight factor W is available in the simulation, which limits the maximum weight w max of the particles. To further optimize the procedure of sampling, different weight factors for electromagnetic and heavy particles are defined. These translate into different maximum weights w max for different particle species, which also depend on the primary energy and the thinning energy. For instance, for a primary energy of 10 19 eV, and a thinning level R thin = 10 -6 the maximum weight for electrons is approximately w max ∼ 10 6 while for muons w max ∼ 10 4 .

We have simulated proton and iron-induced showers with primary energy E p = 10 19 eV, zenith angle θ = 0 • , 30 • , 45 • and 60 • and relative thinning R thin = 10 -5 , 10 -6 , 10 -7 and 10 -8 (the two latter R thin for proton only). We have also simulated showers with relative thinning of R thin = 10 -7 but with weight factor for electrons W (EM) = 0.1 which translates into a smaller maximum weight than the default value given above. Finally, we have also simulated two sets of proton and iron-induced showers at θ = 0 • , with fixed depth of first interaction, starting at the corresponding mean interaction depth for protons and iron at 10 19 eV.

Fluctuations in EAS

In a real shower or in a simulation of an EAS, there are a number of different fluctuations that can occur. Rather generally, we can make a simple classification as shown in Table 1.

"Physical fluctuations" are those due to physical processes in the shower. Here we split them into those due to the first interaction, and those occuring in the secondary interactions, as is customary, and because it has recently been suggested that "universal" shower properties may emerge when considering only the fluctuations in the first interaction point [10]. Physical fluctuations occuring in the first interaction are further divided into those affecting the depth of the first interaction, and those that arise from fluctuations of multiplicity or inelasticity also in the first interaction.

Physical fluctuations -Depth of first interaction.

-Multiplicity, inelasticity, etc, in 1 st interaction.

-Secondary interactions.

Experimental fluctuations -Detector response.

-Sampling fluctuations.

Artificial fluctuations -Thinning.

-Un-thinning.

Table 1 Classification of the fluctuations in a shower, arising from the physical processes in the shower and the measurement process, and those that appear only in shower simulations.

In the case of real data, fluctuations are enlarged due to the detector response, and to the fact that the detector usually only samples a small fraction of the shower front. This "sampling fluctuation" is a statistically well known problem, and sampling fluctuations are rather well studied [11,12]. We will not consider them in this work. Also the detector response introduces an additional source of fluctuations, which are detector dependent, and will not be considered here.

On the other hand, Monte Carlo simulated data is affected by artificial fluctuations due to the thinning and un-thinning (re-sampling) procedures. For the purposes of this work we do not need to consider the effect of fluctuations induced by the unthinning procedure [11].

4 Fluctuations of the longitudinal and lateral shower development

Fluctuations of the longitudinal profile

In Fig. 1 we show the average longitudinal profile of the number of electrons (left panels) and muons (right panels) N, obtained in simulations of 100 proton-induced showers with E p = 10 19 eV, for thinning levels R thin = 10 -6 and 10 -7 and θ = 0 • , 60 • . Also shown are the relative shower to shower fluctuations σ/ N for electrons and muons. As it is well known [START_REF] Gaisser | Cosmic Rays and Particle Physics[END_REF], the relative fluctuation has a minimum close to the depth of shower maximum. Also and as it is apparent from the figures, the dependence of the relative fluctuation on the thinning is small, at least for depths close to the depth of maximum.

In Fig. 2 we show for the same showers in Fig. 1, the skewness and the kurtosis of the distribution of the number of particles N at different depths. It is worth recalling that the skewness of the distribution of a variable x is defined as

γ 3 = (x -x) 3 σ 3 x , (1)
where x (σ x ) is the average (standard deviation) of x. The kurtosis is defined as

γ 4 = (x -x) 4 σ 4 x -3. (2)
where the "-3" in the definition is a convention to make γ 4 = 0 for a Gaussian distribution. Both the kurtosis and the skewness can be positive or negative.

The skewness is a measure of the asymmetry of the distribution with respect to the mean value. A negative sign implies that the distribution is "deformed" towards values of x smaller than the mean. The contrary applies for a positive sign. The kurtosis is a measure of the length of the tails of the distribution.

Positive values imply that the distribution has tails longer than those of a Gaussian, while negative values imply that the tails are shorter (for instance a flat distribution, a box, has kurtosis -1.2). A Gaussian distribution has γ 3 = γ 4 = 0.

Several remarks can be made from Fig. 2. Firstly, it is apparent that both, the skewness and the kurtosis of the distribution of the number of particles depend strongly on the thinning level, contrary to what happens to the mean N and to the relative fluctuations σ/ N. Close to the depth of shower maximum both the skewness and the kurtosis have local extrema very different from zero, implying that the distribution is strongly non-Gaussian. The skewness is negative and this implies that the distribution is asymmetric towards smaller values of N than average. The positive values of the kurtosis imply that the distribution of N has tails longer than those of a Gaussian, at least close to shower maximum. Remarkably, the log-Gaussian distribution, widely used to parameterize fluctuations in the number of electrons, has both γ 3 > 0 and γ 4 > 0, while the fluctuations predicted by Monte Carlo simulations near the maximum of the shower, have negative skewness. For muons and at large depths the skewness is close to zero, so that a Gaussian or a log-Gaussian distribution is a good approximation. For electrons, this is never the case.

Fluctuations of the lateral profile at ground

Of special importance for cosmic-ray physics performed with arrays of detectors is the study of fluctuations in the number of particles at ground.

In Fig. 3 we show the relative fluctuations (σ/ N) of the total number of electrons (left panels) and muons (right panels) at ground (upper panels), and in a ring of width Δr at a distance r = 1000 m from the shower axis (lower panels). In both cases the fluctuations are shown as a function of the number of showers simulated. The ring was taken from r min = 912 m to r max = 1092 m, i.e. Δr = 180 m corresponding to a symmetric interval in the logarithm of r around r = 1000 m, chosen so that it compensates the decreasing density of particles with a larger area as r increases. As expected, fluctuations in the ring Δr are larger than the fluctuations in the whole ground. Also, the fluctuations in the ring have a stronger dependence with the thinning level used than those in the whole ground. This is easy to understand. An entry of weight w falling in the ring represents w particles, so that by losing or gaining just a single entry, one would lose or gain w particles and the fluctuations are enlarged. This effect is not so strong when accounting for all the particles falling anywhere on the ground. In Fig. 3 it can also be clearly seen that no reliable evaluation of the fluctuations can be done with less than about 20 simulated showers, especially in the case of fluctuations in the ring. It can also be seen that a thinning level of R thin = 10 -6 or larger, introduces large artificial fluctuations, so that the shower to shower physical fluctuations cannot be evaluated reliably. This is however no obstacle to approximately estimate the physical shower to shower fluctuations as will be shown in the following.

Physical shower to shower fluctuations at ground

In the Appendix we prove that the distribution of the number of particles N as obtained in thinned Monte Carlo simulations of extensive air showers has a mean N and a standard deviation σ given by:

N = Ne w, (3)
and,

σ 2 = Ne Ω 2 + w2 s 2 . ( 4
)
where:

• Ne and s are respectively the mean and the standard deviation of the distribution of the number of entries N e falling in a given ring around shower axis in each shower, i.e., the distribution of the number of non-thinned (explicitely sampled) particles in the simulation. The standard deviation of N e is simply the RMS of the distribution of N e over many showers. • w and Ω are respectively the mean and the standard deviation of the distribution of weights assigned to the entries. The standard deviation of w can be calculated from the RMS of the distribution of weights in a single shower (provided there are enough entries in that particular shower falling in the region of interest at ground), or from the RMS of the distribution of weights over many showers.

Of course Eq. ( 3) is exact since the thinning algorithm is designed to reproduce it. For Eq. ( 4), the proof only assumes that the probability for an entry to have a given weight w is independent of the probability of a shower to have a given number of entries N e . This is only approximate since the total number of entries and their weights are constrained by energy conservation.

One can interpret Eq. ( 4) as follows. If all entries (sampled particles) had the same weight equal to w in all the simulated showers, then the distribution of weights would not fluctuate from shower to shower and we would get Ω = 0.

In this limit we would obviously have σ = ws. This special case of Eq. ( 4) was also found in [START_REF] Risse | Proc. of 27th ICRC 2001[END_REF]. In the particular case in which all weights are equal to w = 1, i.e. the shower is fully simulated and the thinning procedure is not applied, then clearly N = N e , σ = s and the fluctuations would be obviously dominated by the true, physical shower to shower fluctuations. In the opposite limit, if we imagine that showers always have the same number of entries equal to Ne , i.e. there are no physical shower to shower fluctuations, then we would get s = 0, and the fluctuation in the number of particles would be solely due to the fluctuations of the weight of the entries, i.e., the fluctuations would be dominated by the thinning procedure and σ 2 = Ne Ω 2 . Clearly in these limits the first term of Eq. ( 4) can be identified with the artificial fluctuations introduced by the thinning procedure (σ 2 thin = Ne Ω 2 ), and the second term with the true (physical) shower to shower fluctuations (σ 2 phys = w2 s 2 ).

In the following, and by means of our Monte Carlo simulations, we numerically study to which extent the two terms in Eq. ( 4) allow us to split the shower-toshower fluctuations into artificial and true ones (as the limiting cases suggest), and hence allow us to estimate the effect of thinning in a particular set of simulations. This will clearly depend on the thinning level as shown below.

Before that, we numerically demonstrate that Eq. ( 4) accounts for all the fluctuations (artificial and physical) appearing in simulations of EAS with thinning.

To see this in detail, we have calculated in our Monte Carlo simulations the average weight w and the sigma of the distribution of weights Ω, the average number of entries Ne and the corresponding sigma of its distribution s, as well as the average number of particles N and the sigma of its distribution σ, both for electrons and muons, and compared them to what is predicted by Eq. ( 4).

Firstly, we have verified that the average number N of electrons and muons at different distances to the shower axis r are rather independent of the thinning level, as expected since the thinning algorithm is designed to reproduce the average values of the numbers of particles in the shower. This is however not the case for the relative fluctuations σ/ N shown in Fig. 4, which depend strongly on the thinning level used in the simulations. For both electron and muons σ/ N has a minimum at the distance at which the number of particles is largest. Also, as expected, the fluctuations can be seen to converge to a common value at each distance as the thinning level decreases, because the effect of thinning is increasingly less important, fluctuations are increasingly dominated by the physical shower to shower fluctuations. The artificial fluctuations introduced by thinning in the case of electrons, do not contribute equally to the total fluctuation at all distances from the core as expected.

It can be seen for instance that the relative fluctuation rises with r but the increase is smaller the smaller the thinning level.

From our simulations we also obtain, as expected, that the average weight w assigned in the process of thinning to electrons and muons is simply proportional to the thinning level w ∝ R thin . In Fig. 5 we plot the relative fluctuations of the distribution of weights Ω/ w for electrons, and muons as a function of distance to the shower core. As expected the fluctuation of the weight decreases as R thin decreases and the showers are less thinned. Also, for values of R thin < 10 -5 the relative fluctuation Ω/ w is roughly independent of the thinning level, and as a consequence we have that approximately Ω ∝ R thin .

Regarding the average number of entries Ne , we have verified in the simulations that Ne ∝ R -1 thin , as imposed by the constraint in Eq. ( 3) that the average number of particles N should be independent of R thin , together with the fact that w ∝ R thin . In Fig. 6, we show the relative fluctuation in the number of entries s/ Ne . For small thinning levels (R thin < 10 -6 ), we find that the relative fluctuation is approximately independent of the thinning level, implying that s ∝ R -1 thin .

Finally, in Fig. 7, we compare the relative fluctuation of the number of particles σ/ N obtained directly in Monte Carlo simulations, with that predicted by Eq. ( 4), using the values of w, Ω, Ne and s obtained in the same simulations.

The comparison is shown for thinning levels R thin = 10 -5 , 10 -6 and 10 -7 .

The agreement between the σ/ N obtained in Monte Carlo simulations, and that predicted by Eq. ( 4) is at the level of < 20% for electrons and < 5% for muons, confirming that Eq. ( 4) accounts for all the fluctuations (artificial and physical) appearing in the simulations of EAS with thinning.

In Fig. 8 we plot the two terms in Eq. 4 as a function of the distance to the shower core for different thinning levels. For thinning level R thin ∼ 10 -7 (and smaller) the artificial fluctuations introduced by thinning are very small, and we approach the limiting case explained above in which the identification of the second term of Eq. 4 with the physical fluctuations is exact. For intermediate thinning levels R thin ∼ 10 -6 the artificial fluctuations are a factor of 3 larger than those at R thin = 10 -7 , however the second term of Eq. ( 4) predicts that the fluctuations σ phys are very similar to those obtained with smaller thinning levels, despite the fact that the artificial fluctuations are not negligible -as a matter of fact they are of the same order of the physical fluctuations. This gives support to our identification of σ phys with the true physical fluctuations even for a relatively large thinning level such as 10 -6 which is commonly used in shower simulations. Finally, for thinning levels R thin ∼ 10 -5 or larger, artificial fluctuations dominate, the physical fluctuations are almost completely erased, and the second term in Eq. ( 4) cannot reproduce them.

From the scalings with R thin of the different magnitudes involved in Eq. ( 4) obtained before, it is straightforward to deduce that σ thin ∝ R thin , while σ phys should be approximately independent of R thin for thinning levels 10 -6 and smaller. This is also seen in Fig. 8.

To further check the conclusions above, we have run full simulations (i.e. with no thinning) of air showers and we have compared the shower-to-shower fluctuations obtained in this case -which are solely due to the physics fluctuations since the thinning procedure is not applied -with the fluctuations predicted by the second term of Eq. ( 4) in simulations with different thinning levels.

Given that it is impractical to run full simulations at an energy of 10 19 eV, we have instead run 100 proton showers at lower energies, namely 10 15 eV. Our results are shown in Fig. 9 where it can be seen that the physics fluctuations predicted by the second term of Eq. ( 4) reproduce the fluctuations in showers with no thinning -at the 25% level or better -for thinning levels ∼ 10 -6 and smaller, in agreement with our previous conclusions.

After having shown that the second term in Eq. ( 4) reproduces the physical shower-to-shower fluctuations for thinning levels R thin ∼ 10 -6 and smaller, in the next two sections we study the behaviour of physical fluctuations.

Dependence of fluctuations at ground on the number of particles

Let us now consider the dependence of fluctuations on the size of the ring around a distance to the shower core r, where particles are collected in the simulation. Let us first consider how the density of particles is evaluated.

If N(r, Δr) is the average number of particles at distance r in a small bin Δr (here we assume cylindrical symmetry around the shower axis, but the argument does not depend on this simplification), then the density of particles ρ(r) can be defined as

ρ(r) = lim Δr→0 N (r, Δr) 2πrΔr . ( 5
)
In the limit Δr → 0, ρ(r) is finite, at least for r = 0. However, the same is not true for the fluctuations, so that in general one can not define a "density of fluctuations" ρ σ (r). We can see this in a simple example. Assume that σ(r, Δr) is the standard deviation of the distribution of the number of particles in a bin of size Δr and at a distance r. One could try to define the density of fluctuations as

ρ σ (r) = lim Δr→0 σ(r, Δr) 2πrΔr . ( 6
)
If the fluctuations in the number of particles at ground were purely Poissonian, we would have σ(r, Δr) = N(r, Δr), [START_REF] Sciutto | AIRES User's Manual and Reference Guide[END_REF] and therefore,

ρ σ (r) = lim Δr→0 N(r, Δr) 2πrΔr = lim Δr→0 2πrΔrρ(r) 2πrΔr = lim Δr→0 ρ(r) √ 2πrΔr → ∞, (8)
i.e., we can not define a density of fluctuations for Poissonian shower to shower fluctuations in the number of particles. The dependence on bin-size has been identified with the fractal structure of showers [START_REF] Kempa | [END_REF]. In the case of Poissonian fluctuations in the number of particles, the deduced behaviour is ρ σ (r) ∝ Δr α with α = -1/2, and we would be tempted to identify the coefficient α with a fractal exponent. However notice that for a Poissonian process no fractal structure is implied at all and it would be erroneous to call it a fractal exponent.

On the other hand, in the case in which the fluctuations behave as:

σ(r, Δr) = f (r)Δr + O(Δr 2 ), (9)
where f (r) is a function that does not depend on Δr, one can define a density of fluctuations as can be shown trivially applying Eq. ( 6). Remarkably, this is precisely the case of Furry's fluctuations [16], in which σ ∼ N = 2πrρ(r)Δr and then ρ σ would be independent of Δr, i.e. α = 0. Recall that Furry statistics appears as a extremely simplified model of shower fluctuations [16], but it does take into account the branching structure of the shower (and therefore has an implicit fractal structure included).

To our knowledge, the actual behaviour of the fluctuations of showers and its dependence with the bin size is an open theoretical problem, with the theoretical prejudice ranging between those two extremes: purely Poissonian fluctuations σ ∼ √ N , and Furry fluctuations σ ∼ N. For instance, for the longitudinal development of showers one can show that in fact both types of behaviour occur [17], i.e.,

σ 2 = a N 2 + b N, ( 10
)
where a and b vary slowly with primary energy [17]. Near the maximum of the shower, the first term dominates and fluctuations are not Poissonian. For the lateral distribution no such result exists but one would expect a similar conclusion.

In Fig. 10 we show the relative fluctuations in the number of particles as predicted by the two terms in Eq. ( 4) for thinning level R thin = 10 -6 . for which we have shown before that the two terms in Eq. ( 4) predict the physical and artificial fluctuations. It can be seen that σ phys / N is consistent with being flat with Δr, while σ thin / N behaves as a power law (a fit gives σ thin / N ∝ Δr -1/2 ). These results suggest that the artificial fluctuations are Poissonian, while physical fluctuations behave as σ phys ∝ N.

Dependence of shower to shower fluctuations on composition and depth of first interaction

In this Section we study the influence of the fluctuations in the depth of first interaction on the overall shower to shower fluctuations of the number of particles. In Fig. 11 we plot the relative fluctuations σ/ N in 10 19 eV proton and iron-induced showers. In all panels we show the results of our regular simulations, together with the results of a special set of simulations performed by fixing the depth of first interaction of the primary particle (proton or iron) at the value of its mean interaction depth predicted by the QGSJET model (namely 44.9 g/cm 2 for proton at 10 19 eV and 10.7 g/cm 2 for iron at the same energy). In all cases we use Eq. ( 4) to split the fluctuations into artificial and physical fluctuations, and we also show them in the figures.

Firstly, it is interesting to see that the artificial fluctuations in the number of electrons in iron showers are approximately equal to those in proton showers, while the physical fluctuations are smaller in iron than in proton-induced showers. The latter observation is a well-known effect which is attributed to the fact that showers initiated by a nuclei can be considered, in a first approximation, as a superposition of A (atomic mass) nucleons, each with an energy E/A with E the energy of the primary nucleus.

It is rather remarkable that the relative fluctuations in the number of particles σ/ N in the two different sets of simulations (fixing or varying the depth of the first interaction point), are essentially the same. One could think that this is due to the fluctuations induced by thinning which mask the effect of the fluctuations of the depth of first interaction, however this does not seem to be the case, since as can be seen in Fig. 11 neither the first term of Eq. ( 4) (the thinning fluctuations), nor the second term (the physical fluctuations) change much when varying or fixing the depth of first interaction. We conclude that the relative shower to shower fluctuations of the number of particles at ground σ/ N are rather insensitive to the physical fluctuations of the depth of first interaction, and that these are of the same order as the fluctuations occuring in the subsequent secondary interactions in the shower. This is related to the fact that the maximum of a shower at θ = 0 • , where the fluctuations are minimum [START_REF] Gaisser | Cosmic Rays and Particle Physics[END_REF], occurs near the ground. For other zenith angles a small difference appears in the physical fluctuations of the simulations performed with fixed and fluctuated first interaction point.

Conclusions

In this work we have performed a comprehensive study of shower to shower fluctuations by means of Monte Carlo simulations of extensive air showers. An understanding of the shower to shower fluctuations will help to improve the interpretation of cosmic-ray data.

The determination of the true, physical shower to shower fluctuations is hampered by the thinning procedure necessary to simulate in a practical manner air showers at EeV energies and above. However, we have shown that the artificial fluctuations induced by thinning (σ thin ) can be identified and splitted from the physical fluctuations (σ phys ) with the aid of Eq. ( 4), provided thinning levels R thin ∼ 10 -6 or smaller are used in the simulations. It is important to remark that for R thin = 10 -6 the artificial fluctuations are not negligible, in fact they are of the same order as the physical fluctuations, and hence Eq. ( 4) is able to reproduce them for a thinning level commonly used in shower simulations. We have also shown that Eq. ( 4) accounts for all, true and artificial fluctuations appearing in the simulations regardless of the thinning level. Moreover, Eq. ( 4) reproduces the expectation that as the thinning level decreases R thin → 0, and showers are less thinned, then the artificial fluctuations decrease, the physical ones become dominant, and they do not depend on R thin .

Eq. ( 4) indicates that the physical shower to shower fluctuations of the number of particles at ground behave proportionally to the number of particles N, while the artificial fluctuations are Poissonian, i.e., behave as √ N .

We have also shown that the size of the relative fluctuations due to the depth at which the first interaction initiating the shower occurs, is smaller or of the same order as the fluctuations occuring in the subsequent secondary interactions in the shower.

Appendix

In this Appendix, we calculate the probability distribution of the number of particles N in a given bin of r distance to shower core, or of energy. making some simplifying assumptions. We want to calculate the probability distribution of particles, possibly in a given bin of r or of energy. We will assume that the probability for an entry (a non-thinned particle) to have a weight w i is given by P w (w i ). In addition the number of entries, N e in a given shower is a random variable with probability distribution P e (N e ). Our main simplifying assumption is the following: we will assume that P w and P e are independent of each other. This assumption is only approximate, because in a shower simulated with thinning both the entries and the weight assigned to each particle are controlled by the branching of the shower and therefore they must be related. However, as we will see a posteriori, the approximation is good enough for our purposes here, and it serves to clarify the role of thinning.

Under this approximation we can write the probability P (N) of having N particles as,

P (N) = Ne P e (N e ) dw 1 • • • dw Ne P w (w 1 ) • • • P w (w Ne ) δ(w 1 + • • • + w Ne -N) .
(11) where the δ-function expresses the constraint that the sum of weights is equal to the total number of particles.

In what follows we evaluate this expression first by making further assumptions about the shape of P e and P w , and afterwards in the general case using the characteristic function, related to the probability distribution. The definitions of cumulants and of the characteristic function can be found in any text book on statistics, for instance [START_REF] Eadie | Statistical Methods in experimental physics[END_REF].

We start with the integral (12) and introduce the Fourier representation for the delta function.

I(N e , N) = dw 1 • • • dw Ne P w (w 1 ) • • • P w (w Ne ) δ(w 1 + • • • + w Ne -N),
I(N e , N) = dw 1 • • • dw Ne P w (w 1 ) • • • P w (w Ne ) 1 2π dk e ik(w 1 +•••+w Ne -N ) .
(13) Changing the order of integration gives

I(N e , N) = 1 2π
dk e -ik N dw 1 P w (w 1 )e ik w

) 14
To further continue with the evaluation of P (N) we need to make additional approximations. We assume that P w is a Gaussian distribution with average w and rms Ω P w (w

) = A e -(w-w) 2 /(2Ω 2 ) , ( 15
)
where A = 1/ √ 2πΩ 2 is the probability normalization. Its Fourier transformation is given by dw e ikw P w (w) = e ik w e -k 2 Ω 2 /2 . (

Therefore,

I(N e , N) = 1 2π dk e -ik(N -Ne w) e -Nek 2 Ω 2 /2 , ( 17
)
For large values of Ne , the sum in Eq. ( 11) can be approximated by an integral

) 18
If we assume that P e (N e ) is also a Gaussian with average Ne and standard deviation s and inserting Eq. ( 17) in Eq. [START_REF] Eadie | Statistical Methods in experimental physics[END_REF] gives

P (N) = 1 2π dk e -ikN dN e 1 √ 2πs 2 e -(Ne-Ne) 2 /(2s 2 ) e ikNe w e -k 2 NeΩ 2 /2 . ( 19
)
The integral over N e can be done analytically

P (N) = 1 2π dk e -ikN exp[ i Ne wk - 1 2 ( Ne Ω 2 + s 2 w2 )k 2 + O(k 3 ) ]. ( 20
)
where we neglect in the exponential powers of k larger than 2, after applying the saddle point approximation. We arrive at the final expression

P (N) = 1 √ 2πσ 2 e -(N -N ) 2 /(2σ 2 ) , (21)
where N = Ne w,

σ 2 = Ne Ω 2 + w2 s 2 . ( 22
)
Notice that the above expressions have the correct asymptotic behaviour. If the particles have no weight, w → 1 and Ω → 0, then the average number of particles is equal to the average number of entries (non-thinned particles) N = Ne and σ = s. On the extreme case of a strongly thinned shower in which all particles are grouped together in a single entry ( Ne = 1, s = 0) then N = w and σ = Ω, as expected.

The above result is general and valid for any probability distribution for P w and P e . The only requirement is the "factorization" property given in Eq. ( 11). From Eq. ( 14), we introduce the characteristic function for the probability distribution P w , Pw (k) = dw P w (w) e ikw , (

which, in general can be written as

Pw (k) = exp[ ika 1 - 1 2 a 2 k 2 + • • • ] = e ig(k) , (24)
where the coefficients of the expansion of g(k) are related to the cumulants of the distribution of P w . For instance a 1 = w, a 2 = Ω 2 , etc. Then Eq. ( 11) reads

P (N) = 1 2π
dk e -ikN dN e P e (N e ) e iNeg (k) .

(25)

We define Pe (q) = dN e P e (N e ) e iqNe = exp[

ib 1 q - 1 2 b 2 q 2 + • • • ], (26)
where as before b 1 = Ne and b 2 = s 2 . Then we get

P (N) = 1 2π dk e -ikN Pe (g(k)) = 1 2π dk e -ikN exp[ ib 1 g(k) - 1 2 b 2 g(k) 2 + • • • ].
Where the function g

(k) = k w + i/2 k 2 Ω 2 + • • • .
Then after some algebra

P (N) = 1 2π dk e -ikN exp[ ik w Ne - 1 2 k 2 (s 2 w2 + Ne Ω 2 ) + • • • ] (27)
The coefficients of the expansion around k = 0 are again the cumulants of the distribution P (N), therefore we simply read the result given above in Eq. 22.

But, as a bonus, we obtain also all the other cumulants. For instance for the skewness we obtain

γ 3 = M 3 σ 3 = 1 σ 3 × (m 3 Ne + 3 ws 2 Ω 2 + w3 M e ), (28)
where m 3 is the third central moment ( m 3 = (ww) 3 ) ) of the weight distribution and M e is the third central moment of the distribution of the number of entries. In the same way, one can easily obtain other cumulants from the above expressions.

P

  (N) = Ne P e (N e ) I(N e , N) ≈ dN e P (N e ) I(N e , N).(

7 Figure 1 . 7 Figure 2 .

 7172 Figure 1. Upper panels: Longitudinal development of the average number of electrons (left) and muons (right) as a function of the slanted depth. Lower panels: Relative fluctuations σ/ N as a function of the slant depth for electrons (left) and muons (right). In all panels: 100 proton-induced showers with E p = 10 19 eV and θ = 0 • and 60 • , were simulated with relative thinning R thin = 10 -6 , 10 -7 and 10 -8 .

Figure 4 .Figure 5 .Figure 6 . 4 .

 4564 Figure 4. The relative fluctuation σ/ N in the number of electrons (upper panel) and muons (lower panel) at ground, as a function of the logarithm of the distance to the shower core. 100 proton-showers of energy E p = 10 19 eV and θ = 0 • were simulated with different thinning levels as indicated in the insets.

Figure 7 .

 7 Figure 7. Relative fluctuations σ/ N of the distribution of number of electron (upper panel) and muons (lower panel) at ground, as a function of the logarithm of the distance to the core, for the same sets of shower simulations as in Fig. 4. The σ/ N obtained in the simulations (MC) is compared to that predicted by Eq. (4).

Figure 10 .

 10 Figure 10. Relative fluctuations σ/ N of the distribution of the number of electrons (upper panel) and muons (lower panel) in a ring of width Δr, centered at r = 100 m at ground, as a function of the logarithm of size of the bin Δr. The simulations were performed with relative thinning level R thin = 10 -6 . The two terms in Eq. (4) corresponding to fluctuations induced by thinning σ thin and physical fluctutations σ phys are also shown, see insets.

  1

  • • • dw Ne P w (w Ne )e ik w Ne

	=	1 2π	dke -ik N	dwP w (w)e ikw	Ne	.	(

Acknowledgements

We thank V. Canoa, G. Rodriguez-Fernandez, T. Tarutina, I. Valiño and E. Zas for discussions and comments. P. M. H. was supported by Juan de la Cierva grant. We thank Centro de Supercomputación de Galicia (CESGA) for computer resources. This work was made possible with support from the Ministerio de Ciencia e Innovación, Spain under grant FPA 2007-65114 and Consolider CPAN; and of ALFA-EC funds in the framework of the HELEN (High Energy Physics Latin-American-European Network) project. J. A-M also thanks Xunta de Galicia (INCITE09 206 336 PR) for financial support.

 4)) of the distribution of number of electrons at ground, as a function of the logarithm of the distance to the core, for the same sets of shower simulations as in Fig. 4. Figure 11. Relative fluctuations σ/ N of the distribution of number of electrons at ground, as a function of the logarithm of the distance to the core, for 10 19 eV proton showers (upper panel) and iron showers (lower panel) with θ = 0 • , and relative thinning level R thin = 10 -6 . We show (squares) the result of fixing the first interaction depth at 44.9 g/cm 2 (mean interaction depth of 10 19 eV proton-air collisions predicted by the QGSJET01 model -upper panel) or 10.7 g/cm 2 (corresponding to the mean interaction depth of 10 19 eV iron-air collisions predicted by the QGSJET01 model -lower panel), and also the case in which the depth of first interaction fluctuates (triangles). The two terms in Eq. ( 4) corresponding to fluctuations induced by thinning σ thin and physical fluctutations σ phys are also shown in all cases, see insets.