C Barrabès
email: barrabes@celfi.phys.univ-tours.fr

P A Hogan

Generating a Cosmological Constant with Gravitational Waves

A technique is given to derive the well known Bell-Szekeres solution of the Einstein-Maxwell vacuum field equations describing the space-time and the Maxwell field following the head-on collision of two homogeneous, plane, electromagnetic shock waves. The analogue of this technique is then utilized to construct the space-time model of the gravitational field following the head-on collision of two homogeneous, plane, gravitational shock waves. The latter collision, which is followed by a pair of impulsive gravitational waves and a pair of light like shells traveling away from each other, provides a mechanism for generating a cosmological constant which may be important in the theoretical description of dark energy.

Introduction

This paper is concerned with exploring the strong analogy between the collision of electromagnetic shock waves and the collision of gravitational shock waves and with demonstrating that the latter is a mechanism for generating a cosmological constant, which may be important in the theoretical description of dark energy [1]. It is important in this context to note that the appearance of a cosmological constant term on the left-hand side of the Einstein field equations is equivalent to the appearance of an energy-momentum-stress tensor for a perfect fluid with an equation of state expressing the vanishing of the sum of the matter proper density and the isotropic pressure. Thus for the case of colliding gravitational shock waves the space-time consists of an anti-collision region which is a vacuum and a post-collision region which is a non-vacuum in this sense. Vacuum and non-vacuum regions of space-time are familiar from solving the field equations for so-called interior and exterior solutions. Immediately following the collisions considered here a pair of impulsive gravitational waves is formed, moving away from each other, for the case of colliding electromagnetic shock waves (which is a well-known phenomenon [2]) while following the collision of two gravitational shock waves a pair of impulsive gravitational waves and a pair of light-like shells, moving away from each other, are formed.

The collision of plane electromagnetic shock waves with a Heaviside step function profile has been solved in Einstein-Maxwell theory many years ago by Bell and Szekeres [2]. The collision is head-on and the waves are homogeneous. Such waves are described by a solution of the vacuum Einstein-Maxwell field equations starting with the line element ds 2 = -cos 2 au + (dx 2 + dy 2 ) + 2 du dv ,

where a is a real constant and u + = uϑ(u) with ϑ(u) the Heaviside step function which is equal to zero if u < 0 and equal to unity if u > 0. If we write the line element in terms of the basis 1-forms ϑ 1 = cos au + dx, ϑ 2 = cos au + dy, ϑ 3 = dv and ϑ 4 = du we have

ds 2 = -(ϑ 1 ) 2 -(ϑ 2 ) 2 + 2ϑ 3 ϑ 4 = g ab ϑ a ϑ b , (1.2)
where the constants g ab = g ba are the components of the metric tensor on the half null tetrad defined by the basis 1-forms. Tetrad indices (such as a, b here) will be lowered and raised using g ab and g ab respectively, with g ab the components of the inverse of the matrix with entries g ab and thus g ab g bc = δ a c . The components of the Ricci tensor on this tetrad are given by

R ab = -2a 2 ϑ(u)δ 4 a δ 4 b . (1.3)
The Maxwell field is described by the 2-form

F = 1 2 F ab ϑ a ∧ ϑ b = a ϑ(u)ϑ 1 ∧ ϑ 4 , (1.4)
(or equivalently by the Newman-Penrose [3] component Φ 2 = a ϑ(u)) and the corresponding electromagnetic energy-momentum tensor is

E ab = F ac F c b - 1 4 g ab F dc F dc = -a 2 ϑ(u)δ
F = aϑ(u)ϑ 2 ∧ ϑ 4 , (1.7)
which clearly satisfies the vacuum Maxwell field equations

d * F = 0 , (1.8)
where d denotes exterior differentiation. The Maxwell field (1.4) is type N (the radiative type) in the Petrov classification with the vector field ∂/∂v as degenerate principal null direction and thus represents pure electromagnetic radiation. The profile of the wave is the Heaviside step function and so the wave is a shock wave. The coefficient of the step function in (1.4) is a constant (in particular independent of the coordinates x, y) and thus the waves are homogeneous. The histories of the wave fronts in space-time are the null hyperplanes u = constant which are generated by the null, geodesic, shear-free (and, trivially, twist-free) integral curves of the vector field ∂/∂v.

In the space-time with line-element (1.1) the vector field ∂/∂u is null and generates the null hyperplanes v = constant. The latter can be the histories of the wave fronts of electromagnetic shock waves traveling in the opposite direction to those with histories u = constant and are given by the solution of the Einstein-Maxwell vacuum field equations consisting of the line element

ds 2 = -cos 2 bv + (dx 2 + dy 2 ) + 2 du dv , (1.9)
where b is a real constant and v + = vϑ(v), and the Maxwell field is given by the 2-form

F = 1 2 F ab ϑ a ∧ ϑ b = b ϑ(v)ϑ 3 ∧ ϑ 1 , (1.10)
or equivalently by the Newman-Penrose component Φ 0 = b ϑ(v). Bell and Szekeres considered the collision of two such families of electromagnetic shock waves, in the sense that if the space-time and Maxwell field for v < 0 is given by (1.1) and (1.4) and if the space-time and Maxwell field for u < 0 is given by (1.9) and (1.10) (the space-time for v < 0 and u < 0 is trivially flat Minkowskian space-time with vanishing Maxwell field) then what is the post-collision space-time for u > 0, v > 0 and post-collision Maxwell field for u > 0, v > 0? The post-collision space-time and Maxwell field constitute the Bell-Szekeres solution of the Einstein-Maxwell vacuum field equations. We give a derivation of this solution in section 2 which acts as a guide for solving the analogous problem for colliding gravitational shock waves in section 3.

The space-time model of the gravitational field of gravitational shock waves is given by the line element

ds 2 = -cos 2 au + dx 2 -cosh 2 au + dy 2 + 2 du dv , (1.11)
where a is a real constant. This takes the form (1.2) with now the basis 1-forms given by ϑ 1 = cos au + dx , ϑ 2 = cosh au + dy , ϑ 3 = dv , ϑ 4 = du .

(1.12)

The components R ab of the Ricci tensor on the tetrad defined by the basis 1-forms vanish, so that the metric tensor given via the line element (1.11) is a solution of Einstein's vacuum field equations. The corresponding Riemann curvature tensor components in Newman-Penrose notation are given by

Ψ 0 = Ψ 1 = Ψ 2 = Ψ 3 = 0 and Ψ 4 = a 2 ϑ(u) . (1.13)
This is a Petrov type N curvature tensor and thus represents pure gravitational radiation with propagation direction in space-time the degenerate principal null direction of the Riemann tensor ∂/∂v. We have again here plane, shock waves which are homogeneous and have a step function profile, but in this case the waves are gravitational waves. For similar waves traveling in the opposite direction

ds 2 = -cos 2 bv + dx 2 -cosh 2 bv + dy 2 + 2 du dv , (1.14)
where b is a real constant. The analogue of the Bell and Szekeres problem is to ask: if the space-time for v < 0 is given by (1.11) and if the space-time for u < 0 is given by (1.14) (the space-time for v < 0 and u < 0 is trivially flat Minkowskian space-time) then what is the post-collision space-time for u > 0, v > 0? Following the pattern of the Bell-Szekeres derivation in section 2 below we give the answer to this question in section 3. This is the main result of this paper. The most significant feature of this answer is that the appropriate field equations in the post-collision region of the space-time are Einstein's vacuum field equations with a cosmological constant. The cosmological constant, which could conceivably represent dark energy, is proportional to the product of the constants a and b associated with the incoming gravitational shock waves indicated above. Some important properties of the derivation are discussed in section 4 including a summary of the principal physical attributes of the solution derived in section 3.

Derivation of Bell-Szekeres Solution

For the head-on collision of homogeneous plane waves, the post-collision region u > 0, v > 0 of space-time is well known to take the Rosen-Szekeres form [4], [5], [6]

ds 2 = -e -U (e V dx 2 + e -V dy 2 ) + 2e -M du dv , (2.1)
where U, V, M are each functions of u, v. For the case involving electromagnetic waves, which interests us particularly in this section, the Maxwell field has in general only two Newman-Penrose components Φ 0 (u, v) and Φ 2 (u, v), or equivalently, the Maxwell field is given by a 2-form

F = 1 2 F ab ϑ a ∧ ϑ b = Φ 0 ϑ 3 ∧ ϑ 1 + Φ 2 ϑ 1 ∧ ϑ 4 , (2.2)
which also serves to define Φ 0 and Φ 2 in terms of F ab and the basis 1-forms.

It is convenient to write

Φ 0 = e M/2 φ 0 and Φ 2 = e M/2 φ 2 , (2.3)
and then Maxwell's equations in terms of φ 0 (u, v) and φ 2 (u, v) read:

∂φ 0 ∂u = 1 2 U u φ 0 - 1 2 V v φ 2 , (2.4
)

∂φ 2 ∂v = 1 2 U v φ 2 - 1 2 V u φ 0 , (2.5)
with the subscripts denoting partial differentiation, and the Einstein-Maxwell field equations read:

U uv = U u U v , (2.6) 2U uu = U 2 u + V 2 u -2U u M u + 4φ 2 2 , (2.7) 2U vv = U 2 v + V 2 v -2M v U v + 4φ 2 0 , (2.8) 2V uv = U u V v + U v V u + 4φ 0 φ 2 , (2.9) 2M uv = V u V v -U u U v .
(2.10)

We wish to solve these equations in the region u > 0, v > 0 with the following boundary conditions (which are obtained from (1.1), (1.4) and from (1.9), (1.10) respectively) :

When v = 0 : e -U = cos 2 au , V = 0 , M = 0 , φ 2 = a ; (2.11) When u = 0 : e -U = cos 2 bv , V = 0 , M = 0 , φ 0 = b , (2.12)
where a, b are constants. We note the well known fact that under a change of coordinates u → ū = ū(u) and v → v = v(v) the line element and field equations remain invariant in form with M → M , φ 0 → φ0 , φ 2 → φ2 with

e M = dū du dv dv e M , φ0 = dv dv -1 φ 0 , φ2 = dū du -1 φ 2 . (2.13)
Eq.(2.6) is equivalent to (e -U ) uv = 0 and solving this with

e -U = -1 + cos 2 au + cos 2 bv = cos(au -bv) cos(au + bv) , (2.14)
satisfies the boundary conditions. We see from (2.12) that when u = 0 we have V v = 0 but we shall require V u and φ 2 when u = 0 which we write as (V u ) u=0 and (φ 2 ) u=0 . Also when v = 0 we have V u = 0 and we shall require (V v ) v=0 and (φ 0 ) v=0 . To derive (V u ) u=0 and (φ 2 ) u=0 we use (2.5) and (2.9) evaluated at u = 0. First we note from (2.14) that U u = a tan(au -bv) + a tan(au + bv) ,

(2.15)

U v = -b tan(au -bv) + b tan(au + bv) , (2.16)
and so (U u ) u=0 = 0 and (U v ) u=0 = 2b tan bv. Hence (2.5) and (2.9) evaluated at u = 0 yield

d dv {(φ 2 ) u=0 cos bv} = - b 2 (V u ) u=0 cos bv , (2.17)
d dv {(V u ) u=0 cos bv} = 2b (φ 2 ) u=0 cos bv , (2.18)
respectively. Solving these we arrive at

(φ 2 ) u=0 = P tan bv + Q and (V u ) u=0 = -2P + 2Q tan bv , (2.19)
where P, Q are constants. When u = 0 and v = 0 we have φ 2 = a and V u = 0 and this implies that P = 0 and Q = a. Therefore (φ 2 ) u=0 = a and (V u ) u=0 = 2a tan bv .

(2.20)

Similarly using (2.4) and (2.9) evaluated at v = 0 leads to

(φ 0 ) v=0 = b and (V v ) v=0 = 2b tan au .
(2.21) Now Maxwell's equations can be written in the form

∂ ∂u (log φ 0 ) = 1 2 U u - 1 2 V v φ 2 φ 0 and ∂ ∂v (log φ 2 ) = 1 2 U v - 1 2 V u φ 0 φ 2 , (2.22)
from which we easily deduce that

2 ∂ 2 ∂u∂v log φ 0 φ 2 = ∂ ∂u V u φ 0 φ 2 - ∂ ∂v V v φ 2 φ 0 . ( 2

.23)

We now make the key assumption that

φ 0 φ 2 = A(u) B(v) , (2.24)
for some functions A(u), B(v). On account of (2.13) this means that there exists a frame of reference ū, v such that φ0 = φ2 . This implies that φ2 0 = φ2 2 and so if φ0 , φ2 are thought of as describing two families of backscattered electromagnetic radiation following the collision, then we are making the simplifying assumption that in the barred frame the energy densities of the two families of backscattered waves are equal. From (2.20) and (2.21) it follows that (2.24) implies

φ 0 φ 2 = b a , (2.25)
for u ≥ 0, v ≥ 0. Using this in (2.23) we arrive at the wave equation

1 a ∂ ∂u 1 a V u = 1 b ∂ ∂v 1 b V v , (2.26) or, with ū = au, v = bv, V ūū = V vv . (2.27)
Solving this for V (ū, v) using the d'Alembert formula, with

V (0, v) = 0 and V ū(0, v) = 2 tan v , (2.28)
by (2.12) and (2.20), we find that

V (ū, v) = 1 2 v+ū v-ū 2 tan ξ dξ = log cos(ū -v) cos(ū + v) , (2.29)
and thus

V (u, v) = log cos(au -bv) cos(au + bv) , (2.30) for u ≥ 0, v ≥ 0. Now V u = -a tan(au -bv) + a tan(au + bv) , (2.31) V v = b tan(au -bv) + b tan(au + bv) , (2.32)
and so with (2.15) and (2.16) we find that We can write the line element of the space-time in a unified form incorporating each of the four regions of the space-time I(u < 0, v < 0), II(u > 0, v < 0), III(u < 0, v > 0) and IV(u > 0, v > 0) as

U u U v = V u V v .
ds 2 = -cos 2 (au + -bv + ) dx 2 -cos 2 (au + + bv + ) dy 2 + 2 du dv .
(2.35)

On the tetrad defined by the basis 1-forms ϑ 1 = cos(au + -bv + ) dx, ϑ 2 = cos(au + + bv + ) dy, ϑ 3 = dv and ϑ 4 = du the Ricci tensor components are

R ab = -2 ab ϑ(u) ϑ(v) δ 1 a δ 1 b + 2 ab ϑ(u) ϑ(v) δ 2 a δ 2 b -2 b 2 ϑ(v)δ 3 a δ 3 b -2 a 2 ϑ(u)δ 4 a δ 4 b = 2 E ab , (2.36)
where E ab is the electromagnetic energy-momentum tensor calculated with the Maxwell field

F = a ϑ(u)ϑ 1 ∧ ϑ 4 + b ϑ(v)ϑ 3 ∧ ϑ 1 .
(2.37)

Finally the Newman-Penrose components of the Weyl conformal curvature tensor for the Bell-Szekeres space-time are given by

Ψ 0 = -b δ(v) tan au + , (2.38) Ψ 1 = 0 , (2.39) Ψ 2 = 0 , (2.40) Ψ 3 = 0 , (2.41) Ψ 4 = -a δ(u) tan bv + .
(2.42)

For u > 0, v > 0 this space-time is conformally flat. The appearance of the Dirac delta functions in two of these components means that there is an impulsive gravitational wave with history v = 0, u > 0 described by Ψ 0 and an impulsive gravitational wave with history u = 0, v > 0 described by Ψ 4 . Thus the energy in the incoming electromagnetic shock waves is re-distributed following the collision into two impulsive gravitational waves followed by the superposition of two systems of electromagnetic shock waves described by (2.37).

Collision of Gravitational Shock Waves

Taking the Bell-Szekeres example as our guide we consider now the head-on collision of the homogeneous, plane gravitational waves described by (1.11) and (1.14). Again for the post collision space-time we work with a line element of the form (2.1) but now we are interested in Einstein's field equations with a cosmological constant Λ given by R ab = Λ g ab . Written out explicitly these equations read:

U uv = U u U v -Λ e -M , (3.1
) 2V uv = U u V v + U v V u , (3.2
) 2U uu = U 2 u + V 2 u -2M u U u , (3.3) 2U vv = U 2 v + V 2 v -2M v U v , (3.4
) 2M uv = V u V v -U u U v . (3.5)
The Newman-Penrose components of the Weyl conformal curvature tensor are:

Ψ 0 = - 1 2 e M {V vv + (M v -U v )V v } , (3.6)
Ψ 1 = 0 , (3.7)
Ψ 2 = 1 4 e M (V u V v -U u U v ) + 1 6 Λ , (3.8)
Ψ 3 = 0 , (3.9)
Ψ 4 = - 1 2 e M {V uu + (M u -U u )V u } . (3.10)
The Bianchi identities are:

∂Ψ 0 ∂u = M u + 1 2 U u Ψ 0 - 3 2 V v Ψ 2 , (3.11
)

∂Ψ 4 ∂v = M v + 1 2 U v Ψ 4 - 3 2 V u Ψ 2 , (3.12
)

∂Ψ 2 ∂u = 3 2 U u Ψ 2 - 1 2 V v Ψ 4 , (3.13
)

∂Ψ 2 ∂v = 3 2 U v Ψ 2 - 1 2 V u Ψ 0 . (3.14)
These can be calculated using the formulas given by Chandrasekhar [START_REF] Chandrasekhar | The Mathematical Theory of Black Holes, Oxford Classic Texts in the Physical Sciences[END_REF] or, alternatively, they can be calculated directly by taking the appropriate partial derivatives of (3. To begin with we shall require (U v ) v=0 , (V v ) v=0 , (M v ) v=0 and (Ψ 0 ) v=0 . It is useful to note from the boundary conditions that (U u ) v=0 = a tan au -a tanh au , (3.17) (V u ) v=0 = -a tan au -a tanh au .

(3.18) To determine (V v ) v=0 we evaluate (3.2) at v = 0 and arrive at

Evaluating (3.1) at v = 0 results in d du (U v ) v=0 + a(-tan au + tanh au)(U v ) v=0 = -Λ , (3.19)
d du (V v ) v=0 = a 2 (tan au -tanh au)(V v ) v=0 + Λ 4 (tan au + tanh au) 2 . (3.24)
To solve this we make use of the identity

(V v ) v=0 = Λ 2a (tan au -tanh au) + C 1 (cos au cosh au) -1/2 , (3.29)
where C 1 is a constant of integration. Since V v vanishes when u = 0 and v = 0 we have C 1 = 0 and thus Using this, along with (3.17) and (3.30), in the Bianchi identity (3.11) evaluated at v = 0 we find

(V v ) v=0 = Λ 2a (tan au -tanh au) . ( 3
d du (Ψ 0 ) v=0 - 1 2 a(tan au-tanh au)(Ψ 0 ) v=0 = - Λ 2 8a
(tan au-tanh au) , (3.36) which can be rewritten in the form

d du {(cos au cosh au) 1/2 (Ψ 0 ) v=0 } = Λ 2 4a 2 d du (cos au cosh au) 1/2 , (3.37)
which leads to

(Ψ 0 ) v=0 = Λ 2 4a 2 + C 3 (cos au cosh au) 1/2 , (3.38)
where C 3 is a constant of integration. Now we can summarize the current situation with regard to the Weyl tensor evaluated at v = 0 as follows:

(Ψ 0 ) v=0 = Λ 2 4a 2 + C 3 (cos au cosh au) 1/2 , (Ψ 2 ) v=0 = Λ 6 , (Ψ 4 ) v=0 = a 2 . (3.39)
We will return to these equations after we have discussed, for the present case, the analogue of (2.23) with (2.24). From (3.11) we have

∂ ∂u log Ψ 0 = M u + 1 2 U u - 3 2 V v Ψ 2 Ψ 0 , (3.40)
and from (3.12)

∂ ∂v log Ψ 4 = M v + 1 2 U v - 3 2 V u Ψ 2 Ψ 4 , (3.41)
from which we deduce that 2 3

∂ 2 ∂u∂v log Ψ 4 Ψ 0 = ∂ ∂v V v Ψ 2 Ψ 0 - ∂ ∂u V u Ψ 2 Ψ 4 . ( 3
∂ ∂u log Ψ 2 = 3 2 U u - 1 2 V v Ψ 4 Ψ 2 , (3.43) ∂ ∂v log Ψ 2 = 3 2 U v - 1 2 V u Ψ 0 Ψ 2 , (3.44)
and thus we have

∂ ∂v V v Ψ 4 Ψ 2 - ∂ ∂u V u Ψ 0 Ψ 2 = 0 . (3.45)
We have here two equations, (3.42) and (3.45), where we had one equation (2.23) in the electromagnetic case. If we make the assumption that

Ψ 2 = kΨ 1/2 0 Ψ 1/2 4 , (3.46)
for some constant k to be determined later, we can arrive at a situation exactly analogous to the electromagnetic case. With this assumption (3.42) becomes 2 3k

∂ 2 ∂u∂v log Ψ 4 Ψ 0 = ∂ ∂v   V v Ψ 1/2 4 Ψ 1/2 0   - ∂ ∂u   V u Ψ 1/2 0 Ψ 1/2 4   , (3.47)
while (3.45) now reads

∂ ∂v   V v Ψ 1/2 4 Ψ 1/2 0   - ∂ ∂u   V u Ψ 1/2 0 Ψ 1/2 4   = 0 . (3.48)
Thus we now have the analogue of the electromagnetic case in that

Ψ 4 Ψ 0 = A(u) B(v) , (3.49)
for some functions A(u) and B(v) and then V (u, v) is determined by (3.48) (the analogue of (2.26)).

Before proceeding it is useful to now compute the constant of integration C 3 in (3.38). We see that when (3.46) is evaluated at v = 0 we must have (Ψ 0 ) v=0 constant, since (Ψ 2 ) v=0 and (Ψ 4 ) v=0 are constants given in (3.39). Therefore C 3 = 0 and so

(Ψ 0 ) v=0 = Λ 2 4a 2 . (3.50)
But from (3.16) we have Ψ 0 = b 2 when u = 0 and v = 0 and so Λ 2 = 4a 2 b 2 and we can take

Λ = 2ab , (3.51)
from now on since, in particular, the signs of a and b are free to specify. This demonstrates the importance of including the cosmological constant in the field equations at the beginning of this derivation. Thus we have

(Ψ 0 ) v=0 = b 2 , (Ψ 2 ) v=0 = 1 3 ab , (Ψ 4 ) v=0 = a 2 , (3.52)
and so (3.46) evaluated at v = 0 yields k = 1 3 . We now have in particular (Ψ 0 ) v=0 = b 2 and (Ψ 4 ) v=0 = a 2 and similarly we find that (Ψ 4 ) u=0 = a 2 and (Ψ 0

) u=0 = b 2 . Hence it follows from (3.49) that for u ≥ 0, v ≥ 0, Ψ 4 Ψ 0 = a 2 b 2 . (3.53)
Consequently with ū = au, v = bv we see, following (3.48), that V (ū, v) satisfies the wave equation

V vv = V ūū . (3.54)
Using the d'Alembert formula with (by (3.16) and (3.31))

V (0, v) = log cos v cosh v , and V ū(0, v) = tan v -tanh v , (3.55)
we have

V (ū, v) = 1 2 log cos(ū + v) cos(v -ū) cosh(ū + v) cosh(v -ū) + 1 2 v+ū v-ū tan ξ -tanh ξ dξ , (3.56
) which simplifies to

V = log cos(au -bv) cosh(au + bv) , (3.57)
for u ≥ 0, v ≥ 0. We do not yet know Ψ 0 , Ψ 2 , Ψ 4 individually for u > 0, v > 0. However we do know the following relations between them which hold for u ≥ 0, v ≥ 0:

Ψ 4 Ψ 0 = a 2 b 2 , 3Ψ 2 = Ψ 1/2 0 Ψ 1/2 4 , ( 3
∂Ψ 2 ∂u = M u + 1 2 U u 3 b a Ψ 2 - 3 2 V v Ψ 2 , (3.60) 3 a b ∂Ψ 2 ∂v = M v + 1 2 U v 3 a b Ψ 2 - 3 2 V u Ψ 2 , (3.61) ∂Ψ 2 ∂u = 3 2 U u Ψ 2 - 3 a 2 b V v Ψ 2 , (3.62
)

∂Ψ 2 ∂v = 3 2 U v Ψ 2 - 3 b 2 a V u Ψ 2 . ( 3
M u -U u = - a b V v , (3.64)
(3.61) and (3.63) with Ψ 2 = 0 yield respectively. Since M u = 0 when u = 0 and when v = 0, and in addition M v = 0 when u = 0 and when v = 0, the only solutions of (3.68) and (3.69) which satisfy the boundary conditions are

M v -U v = - b a V u . ( 3
M u = 0 = M v , (3.70

Discussion

We have seen that an important part of the argument in section 3 is the assumption (3.46)

∂ ∂v V v Ψ 2 Ψ 0 - ∂ ∂u V u Ψ 2 Ψ 4 = 0 . (4.1) Putting ξ = Ψ 2 Ψ 1/2 0 Ψ 1/2 4 , (4.2)
we can rewrite (4.1) as

∂ ∂v   V v Ψ 1/2 4 Ψ 1/2 0 ξ   = ∂ ∂u   V u Ψ 1/2 0 Ψ 1/2 4 ξ   . (4.3)
Using (4.2) again we can rewrite (3.45) as

∂ ∂v   V v Ψ 1/2 4 Ψ 1/2 0 ξ -1   = ∂ ∂u   V u Ψ 1/2 0 Ψ 1/2 4 ξ -1   . (4.4)
Writing these equations out more explicitly we have

  V v Ψ 1/2 4 Ψ 1/2 0   ξ v -   V u Ψ 1/2 0 Ψ 1/2 4   ξ u = -    ∂ ∂v   V v Ψ 1/2 4 Ψ 1/2 0   - ∂ ∂u   V u Ψ 1/2 0 Ψ 1/2 4      ξ , (4.5)
-   V v Ψ 1/2 4 Ψ 1/2 0   ξ v +   V u Ψ 1/2 0 Ψ 1/2 4   ξ u = -    ∂ ∂v   V v Ψ 1/2 4 Ψ 1/2 0   - ∂ ∂u   V u Ψ 1/2 0 Ψ 1/2 4      ξ . ( 4
 V v Ψ 1/2 4 Ψ 1/2 0   ξ v -   V u Ψ 1/2 0 Ψ 1/2 4   ξ u = 0 . ( 4
V v ξ v = V ū ξ ū , (4.9)
and

V vv = V ūū , (4.10)
respectively. The general solution of the wave equation (4.10) is

V (ū, v) = f (ū + v) + g(ū -v) , (4.11)
where the functions f, g are arbitrary functions of their arguments. Substituting into (4.9) we have

(ξ ū -ξ v) f ′ + (ξ ū + ξ v) g ′ = 0 , (4.12)
where the primes on f, g denote differentiation of these functions with respect to their arguments. If (4.12) is to hold for f, g arbitrary then we must have ξ ū = 0 = ξ v which implies that ξ = constant and thus (4.2) agrees with (3.46).

As in section 2 we can write the line element of the space-time in section 3 in a unified form incorporating each of the four regions of the space-time I(u < 0, v < 0), II(u > 0, v < 0), III(u < 0, v > 0) and IV(u > 0, v > 0) as ds 2 = -cos 2 (au + -bv + ) dx 2 -cosh 2 (au + + bv + ) dy 2 + 2 du dv . For u > 0, v > 0 this Weyl tensor is Type D in the Petrov classification. The situation following the collision of the two gravitational shock waves is more complicated than in the electromagnetic case. Now the energy in the incoming waves is re-distributed into (1) two light-like shells (see [START_REF] Barrabès | Singular Null Hypersurfaces in General Relativity[END_REF]), which could correspond to bursts of neutrinos, described by the Dirac delta function terms in the Ricci tensor components It is well known that the Bell-Szekeres space-time with line element (2.34) can be written as the sum of the line elements of two 2-dimensional manifolds of equal constant curvature [START_REF] Stephani | Exact Solutions of Einstein's Equations[END_REF]. This is the so-called Bertotti-Robinson space-time [START_REF] Bertotti | [END_REF], [11]. Motivated by this we have derived the line element (3.73) as the sum of the line elements of two 2-dimensional manifolds of constant curvature differing only in sign (see [12], [START_REF] Barrabès | Advanced General Relativity: Gravity Waves, Spinning Particles and Black Holes[END_REF]). This is the so-called Nariai-Bertotti space-time [START_REF] Bertotti | [END_REF], [14]. A particularly simple example of a collision of plane, homogeneous, light-like signals combining gravitational and electromagnetic shock waves in each in-coming signal has recently been given in [15].

  6)-(3.10) and simplifying the results using the field equations (3.1)-(3.5) and using again (3.6)-(3.10). The boundary conditions in this case are: When v = 0 : e -U = cos au cosh au , e V = cos au cosh au , M = 0 , Ψ 4 = a 2 ; (3.15) When u = 0 : e -U = cos bv cosh bv , e V = cos bv cosh bv , M = 0 , Ψ 0 = b 2 , (3.16) where a, b are constants.

which we can rewrite as d du (cos au

  cosh au(U v ) v=0 ) = -Λ cos au cosh au , (3.20) and solve with cos au cosh au(U v ) v=0 = -Λ 2a (sin au cosh au + cos au sinh au) + C , (3.21)where C is a constant of integration. Since U v at u = 0 and v = 0 vanishes we have C = 0 and thus(U v ) v=0 = -Λ 2a (tan au + tanh au) .(3.22)Similarly by evaluating (3.1) at u = 0 we obtain(U u ) u=0 = -Λ 2b(tan bv + tanh bv) .(3.23)

Λ 4 (- a 2 (

 42 tan au + tanh au) 2 = Λ 2a d du (tan au -tanh au) -Λ 4 (tan au -tanh au) 2 . (3.25) This enables us to write (3.24) in the form dW du tan au -tanh au)W = 0 , (3.26) with W = (V v ) v=0 -Λ 2a (tan au -tanh au) . (3.27) Now (3.26) reads d du {(cos au cosh au) 1/2 W } = 0 , (3.28) and thus we arrive at

  .30) Proceeding similarly with (3.2) evaluated at u = 0 we find that (V u ) u=0 = Λ 2b (tan bv -tanh bv) . (3.31) Next we require (M v ) v=0 and this is obtained by evaluating (3.5) at v = 0. From (3.17), (3.18), (3.22) and (3.30) we see that the right hand side of (3.5) evaluated at v = 0 vanishes and so we have d du (M v ) v=0 = 0 , (3.32) and thus (M v ) v=0 = C 2 = constant. But M v vanishes when u = 0 and v = 0 and so C 2 = 0. Hence (M v ) v=0 = 0 . (3.33) Evaluating (3.5) at u = 0 similarly results in (M u ) u=0 = 0 . (3.34) Finally we shall require (Ψ 0 ) v=0 . The expression for Ψ 2 given by (3.8) evaluated at v = 0, remembering (3.15), yields (Ψ 2 )

  .42) Next (3.13) and (3.14) can be written as

  .63) Thus (3.60) and (3.62) with Ψ 2 = 0 yield

  .65) With V given by (3.57) the first of these two equations reads M u -U u = -a tan(au -bv) + a tanh(au + bv) , = ∂ ∂u (log cos(au -bv) + log cosh(au + bv)) . (3.66) Integrating and using the boundary conditions (3.16) results in M -U = log {cos(au -bv) cosh(au + bv)} . (3.67) We easily see that this satisfies (3.65) with again V given by (3.57). At this point we do not know M and U separately for u > 0, v > 0. Substituting for U in terms of M from (3.67) and for V from (3.57) into the field equations (3.3) and (3.4) results in the equations 2 M uu = -M 2 u ,

  ) for u ≥ 0, v ≥ 0 and thus it follows from the boundary conditions (3.15) and (3.16) that M = 0 , (3.71) for u ≥ 0, v ≥ 0. Hence we have from (3.67): U = -log{cos(au -bv) cosh(au + bv)} , (3.72) for u ≥ 0, v ≥ 0. With V, M, U given by (3.57), (3.71) and (3.72) respectively the line element (2.1) reads ds 2 = -cos 2 (au -bv) dx 2 -cosh 2 (au + bv) dy 2 + 2 du dv , (3.73) for u ≥ 0, v ≥ 0. Calculating Ψ 0 , Ψ 2 , Ψ 4 from (3.6), (3.8) and (3.10), using V, M, U given by (3.57), (3.71) and (3.72), results in Ψ 0 = b 2 , Ψ 2 = 1 3 ab , Ψ 4 = a 2 , (3.74) for u ≥ 0, v ≥ 0. It is now straightforward to check that (3.11)-(3.14) are satisfied since b U u = a V v and a U v = b V u with V given by (3.57) and U by (3.72).

. 6 )

 6 Adding (4.5) and (4.6) with ξ = 0 results in (3.48) and then (4.5) and (4.6) each become 

. 7 )- 1 Ψ 4

 714 We now have before us the wave equation(3.48) and this equation (4.7). Under the coordinate transformation u → ū = ū(u) and v → v = v(v), mentioned following eq.(2.12) above, Ψ 0 → Ψ0 = dū du dv dv-1 Ψ 0 , Ψ 2 → Ψ2 = Ψ 2 and Ψ 4 → Ψ4 = dvdv dū du and thus ξ → ξ = ξ. Hence if we choosedū du = (A(u)) 1/2 and dv dv = (B(v)) 1/2 ,(4.8)with A(u), B(v) appearing in the presently assumed separation of variables (3.49) then Ψ4 = Ψ0 and (4.7) and (3.48) simplify to read

(4. 13 )a δ 4 b

 134 On the tetrad defined by the basis 1-forms ϑ 1 = cos(au + -bv + ) dx, ϑ 2 = cosh(au + + bv + ) dy, ϑ 3 = dv and ϑ 4 = du the Ricci tensor components areR ab = b δ(v) (tan au + + tanh au + ) δ 3 a δ 3 b + a δ(u) (tan bv + + tanh bv + ) δ 4 +2 ab ϑ(u) ϑ(v) g ab . (4.14)The Weyl conformal curvature tensor components in Newman-Penrose notation areΨ 0 = b 2 ϑ(v) + 1 2 b (tanh au + -tan au + ) δ(v) ,bv + -tan bv + ) δ(u) . (4.19)

  (4.14),(2) two impulsive gravitational waves described by the Dirac delta function terms in the Weyl tensor components (4.15) and (4.19) and (3) dark energy in the region u > 0, v > 0 described by the Ricci tensor in this region of space-time having the form R ab = Λ g ab with Λ = 2 a b .(4.20)

  Now from (2.10) we have M uv = 0 from which it follows, using the boundary conditions (2.11) and (2.12), that M = 0 for u ≥ 0, v ≥ 0. Substituting U from (2.14) and V from (2.30) into (2.7) and (2.8) with M = 0 results in

φ 0 = b and φ 2 = a ,

(2.33)

for u ≥ 0, v ≥ 0, and now (2.9) is automatically satisfied, as are the Maxwell equations (2.4) and (2.5) since bU u -aV v = 0 = aU v -bV u on account of (2.15), (2.16), (2.31) and (2.32). The line element for u > 0, v > 0 is now

ds 2 = -cos 2 (au -bv) dx 2 -cos 2 (au + bv) dy 2 + 2 du dv .

(2.34)

Equations (2.33) and (2.34) constitute the Bell-Szekeres solution of the Einstein-Maxwell vacuum field equations.

  relating Ψ 2 to Ψ 0 and Ψ 4 . It is clear that, as a consequence of the field equations (3.42) and (3.45), (3.46) implies the separation of variables (3.49) in the function Ψ 4 /Ψ 0 . The converse is also true.

	If the
	separation of variables (3.49) is assumed then, as a consequence of the field
	equations (3.42) and (3.45), the equation (3.46) holds. To see this we note
	that if the separation of variables (3.49) is assumed then (3.42) reads


Charles Knessl
email: knessl@uic.edu

Wojciech Szpankowski

Quicksort algorithm again revisited

Keywords: Algorithms, Analysis of algorithms, Asymptotic analysis, Binary search tree, Quicksort, Sorting

We consider the standard Quicksort algorithm that sorts n distinct keys with all possible n! orderings of keys being equally likely. Equivalently, we analyze the total path length ✂ n in a randomly built binary search tree. Obtaining the limiting distribution of ✂ n is still an outstanding open problem. In this paper, we establish an integral equation for the probability density of the number of comparisons ✂ n . Then, we investigate the large deviations of ✂ n . We shall show that the left tail of the limiting distribution is much "thinner" (i.e., double exponential) than the right tail (which is only exponential).

Our results contain some constants that must be determined numerically. We use formal asymptotic methods of applied mathematics such as the WKB method and matched asymptotics.

Introduction

Hoare's Quicksort algorithm [START_REF] Hoare | Quicksort[END_REF] is the most popular sorting algorithm due to its good performance in practise. The basic algorithm can be briefly described as follows [START_REF] Hoare | Quicksort[END_REF][START_REF] Knuth | The Art of Computer Programming. Sorting and Searching[END_REF][START_REF] Mahmoud | Evolution of Random Search Trees[END_REF]: A partitioning key is selected at random from the unsorted list of keys, and used to partition the keys into two sublists to which the same algorithm is called recursively until the sublists have size one or zero.

To justify the algorithm's good performance in practise, a body of theory was built. First of all, every undergraduate learns in a data structures course that the algorithm sorts "on average" n keys in Θ✄ n logn☎ steps. To be more precise, one assumes that all n! possible orderings of keys are equally likely. It is, however, also known that in the worst case the algorithm needs O✄ n 2 ☎ steps (e.g., think of an input that is given in a decreasing order when the output is printed in an increasing order). Thus, one needs a more detailed probabilistic analysis to understand better the Quicksort behavior. In particular, one would like to know how likely (or rather unlikely) it is for such pathological behavior to occur.

A large body of literature is devoted to analyzing the Quicksort algorithm [START_REF] Cramer | A Note Concerning the Limit Distribution of the Quicksort Algorithm[END_REF][START_REF] Devroye | A Note on the Height of binary Search Trees[END_REF][START_REF] Eddy | How Many Comparisons Does Quicksort Use[END_REF][START_REF] Hennequin | Combinatorial Analysis of Quicksort Algorithm[END_REF][START_REF] Knuth | The Art of Computer Programming. Sorting and Searching[END_REF][START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF][START_REF] Mahmoud | Evolution of Random Search Trees[END_REF][START_REF] Rachev | Probability Metrics and Recursive Algorithms[END_REF][START_REF] Régnier | A Limiting Distribution for Quicksort[END_REF][START_REF] Rösler | A Limit Theorem for Quicksort[END_REF][START_REF] Tan | Some Properties of a Limiting Distribution in Quicksort[END_REF]. However, many aspects of this problem are still largely unsolved. To review what is known and what is still unsolved, we introduce some notation. Let ✝ n denote the number of comparisons needed to sort a random list of length n. It is known that after selecting randomly a key, the two sublists are still "random" (cf. [START_REF] Knuth | The Art of Computer Programming. Sorting and Searching[END_REF]). Clearly, the sorting time depends only on the keys' ranking, so we assume that the input consists of the first n integers ✞ 1✟ 2✟ ✡✠ ☛✠ ☞✠ ✌✟ n✍ , and key k is chosen with probability 1✎ n. Then, the following recurrence holds

✝ n ✏ n ✑ 1 ✒ ✓✝ k ✒ ✔✝ n✕ 1✕ k ✠ (1)
Now, let L n ✄ u☎ ✏ Eu✖ n ✏ ∑ k✗ 0 ✘ ✚✙ ✜✛ ✝ n ✏ k✢ u k be the probability generating function of ✝ n . The above recurrence implies that

L n ✄ u☎ ✏ u n✕ 1 n n✕ 1 ∑ i✣ 0 L i ✄ u☎ L n✕ 1✕ i ✄ u☎ (2)
with L 0 ✄ u☎ ✏ 1. Observe that the same recurrences are obtained when analyzing the total path length ✝ n of a binary search tree built over a random set of n keys (cf. [START_REF] Knuth | The Art of Computer Programming. Sorting and Searching[END_REF][START_REF] Mahmoud | Evolution of Random Search Trees[END_REF]). Finally, let us define a bivariate generating function L✄ z✟ u☎ ✏ ∑ n✗ 0 L n ✄ u☎ z n . Then ( 2) leads to the following partial-differential functional equation

∂L✄ z✟ u☎ ∂z ✏ L 2 ✄ zu✟ u☎ ✤✟ ∂L✄ 0✟ u☎ ∂z ✏ 1✠ (3)
Observe also that L✄ z✟

1☎ ✏ ✄ 1 ✑ z☎ ✕ 1 .
The moments of ✝ n are relatively easy to compute since they are related to derivatives of L n ✄ u☎ at u ✏ 1. Hennequin [START_REF] Hennequin | Combinatorial Analysis of Quicksort Algorithm[END_REF] analyzed these carefully and computed the first five cumulants. He also conjectured an asymptotic formula for the cumulants as n ✥ ∞ which he later proved in [START_REF] Hennequin | Tri Rapide at Arbres de Recherche[END_REF].

The main open problem is to find the limiting distribution of ✝ n . Régnier [START_REF] Régnier | A Limiting Distribution for Quicksort[END_REF] proved that the limiting distribution of ✄ ✌✝ n ✑ E ✛ ✝ n ✢ ✦☎ ✧✎ n exists, while Rösler [START_REF] Rösler | A Limit Theorem for Quicksort[END_REF][START_REF] Rösler | A Fixed Point Theorem for Distributions[END_REF] characterized this limiting distribution as a fixed point of a contraction satisfying a recurrence equation. A partial-differential functional equation seemingly similar to ( 3) was studied recently by Jacquet and Szpankowski [START_REF] Jacquet | Asymptotic Behavior of the Lempel-Ziv Parsing Scheme and Digital Search Trees[END_REF]. They analyzed a digital search tree for which the bivariate generating function L✄ z✟ u☎ (in the so-called symmetric case) satisfies

∂L✄ z✟ u☎ ∂z ✏ L 2 ✄ 1 2 zu✟ u☎ (4)
with L✄ z✟ 0☎ ✏ 1. The above equation was solved asymptotically in [START_REF] Jacquet | Asymptotic Behavior of the Lempel-Ziv Parsing Scheme and Digital Search Trees[END_REF], and this led to a limiting normal distribution of the path length ✝ n in digital search trees. While the above equation and ( 3) look similar, there are crucial differences. Among them, the most important is the contracting factor 1 2 in the right-hand side of the above. Needless to say, we know that ( 3) does not lead to a normal distribution since the third central moment is not asymptotically equal to zero (cf. [START_REF] Mahmoud | Evolution of Random Search Trees[END_REF]). More precisely, the third (and all higher odd) moments of ✄ ✌✝ n ✑ E ✛ ✝ n ✢ ★☎ ☞✎ n does not tend to zero as n ✥ ∞.

In view of the above discussion, a less ambitious goal was set, namely that of computing the large deviations of ✝ n , i.e., ✘ ✩✙ ✧✛ ✌✪ ✝ n ✑ E ✛ ✝ n ✢ ✪ ✬✫ εE ✛ ✝ n ✢ ✭✢ for ε ✮ 0. Hennequin [START_REF] Hennequin | Combinatorial Analysis of Quicksort Algorithm[END_REF] used Chebyshev's inequality to show that the above probability is O✄ 1✎ ✯✄ ε log 2 n☎ ✧☎ . Recently, Rösler [START_REF] Rösler | A Limit Theorem for Quicksort[END_REF] showed that this probability is in fact O✄ n ✕ k ☎ for any fixed k, and soon after McDiarmid and Hayward [START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF] used the powerful method of bounded differences to obtain an even better estimate, namely that the tail is approximately equal to n ✕ 2εlog log n (see the comment after Theorem 1 of Section 2).

In this paper, we obtain some new results for the tail probabilities. First of all, we establish an integral equation for the probability density of ✝ n , and using this we derive a left tail and a right tail of the large deviations of ✝ n . We demonstrate that the left tail is much thinner (i.e., double exponential) than the right tail, which is roughly exponential.

We establish these results using formal asymptotic methods of applied mathematics such as the WKB method and matched asymptotics. By "formal" we mean that we do not rigorously establish error bounds on the various asymptotic expansions. The main merit of these methods is that they can be used to obtain asymptotic information directly from the underlying equations. Similar asymptotic approaches were used for enumeration problems by Knessl and Keller [START_REF] Knessl | Partition Asymptotics from Recursions[END_REF] and Canfield [START_REF] Canfield | From Recursions to Asymptotics: On Szekeres" Formula for the Number of Partitions[END_REF].

The paper is organized as follows. In the next section we describe our main findings and compare them with other known results. In Section 3 we derive the integral equation for the asymptotic probability density of ✝ n . In Section 4 we obtain our large deviations results.

Formulation and Summary of Results

As before, we let ✝ n be the number of key comparisons made when Quicksort sorts n keys. The probability generating function of ✝ n becomes

L n ✄ u☎ ✏ ∞ ∑ k✣ 0 ✘ ✚✙ ✜✛ ✝ n ✏ k✢ u k ✏ E ✛ u✖ n ✢ ✧✠ (5)
The upper limit in this sum may be truncated at k ✏ ✱✰ n 2 ✲ , since this is clearly an upper bound on the number of comparisons needed to sort n keys.

The generating function in ( 5) satisfies ( 2 on the page before) which we repeat below (cf. also [START_REF] Eddy | How Many Comparisons Does Quicksort Use[END_REF][START_REF] Mahmoud | Evolution of Random Search Trees[END_REF][START_REF] Régnier | A Limiting Distribution for Quicksort[END_REF][START_REF] Rösler | A Limit Theorem for Quicksort[END_REF])

L n✳ 1 ✄ u☎ ✏ u n n ✒ 1 n ∑ i✣ 0 L i ✄ u☎ L n✕ i ✄ u☎ ✴✟ L 0 ✄ u☎ ✏ 1✠ (6)
Note that L n ✄ 1☎ ✏ 1 for all n ✫ 0, and that the probability ✘ ✚✙ ✜✛ ✝ n ✏ k✢ may be recovered from the Cauchy integral

✘ ✩✙ ✵✛ ✝ n ✏ k✢ ✏ 1 2πi ✶ C u ✕ k✕ 1 L n ✄ u☎ du✠ ( 7
)
Here C is any closed loop about the origin.

In Section 3, we analyze ( 6) asymptotically for n ✥ ∞ and for various ranges of u. We use asymptotic methods of applied mathematics, such as the WKB method and matched asymptotics [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers[END_REF][START_REF] Froman | JWKB Approximation[END_REF]. The most important scale is where

n ✥ ∞ with u ✑ 1 ✏ O✄ n ✕ 1 ☎ , which corresponds to k ✏ E ✛ ✝ n ✢ ✷✒ O✄ n☎ ✏ 2n logn ✒
O✄ n☎ . Most of the probability mass is concentrated in this range of k. As mentioned before, the existence of a limiting distribution of ✄ ✌✝ n ✑ E ✛ ✝ n ✢ ✦☎ ✧✎ n as n ✥ ∞ was established in [START_REF] Régnier | A Limiting Distribution for Quicksort[END_REF][START_REF] Rösler | A Limit Theorem for Quicksort[END_REF], though there seems to be little known about this distribution (cf. [START_REF] Cramer | A Note Concerning the Limit Distribution of the Quicksort Algorithm[END_REF][START_REF] Eddy | How Many Comparisons Does Quicksort Use[END_REF][START_REF] Hennequin | Combinatorial Analysis of Quicksort Algorithm[END_REF][START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF][START_REF] Tan | Some Properties of a Limiting Distribution in Quicksort[END_REF]). Numerical and simulation results in [START_REF] Cramer | A Note Concerning the Limit Distribution of the Quicksort Algorithm[END_REF][START_REF] Eddy | How Many Comparisons Does Quicksort Use[END_REF] show that the distribution is highly asymmetric and that the right tail seems much thicker than the left tail. It is also of interest to estimate these tails (cf. [START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF][START_REF] Rösler | A Limit Theorem for Quicksort[END_REF]), as they give the probabilities that the number of key comparisons will deviate significantly from E ✛ ✝ n ✢ , which is well known to be asymptotically equal to 2n logn as n ✥ ∞ (cf. [START_REF] Hennequin | Combinatorial Analysis of Quicksort Algorithm[END_REF][START_REF] Mahmoud | Evolution of Random Search Trees[END_REF]).

For u ✑ 1 ✏ w✎ n ✏ O✄ n ✕ 1 ☎ and n ✥ ∞, we derive in Section 3 the asymptotic expansion

L n ✄ u☎ ✏ exp ✄ A n w✎ n☎ ✤✸ G 0 ✄ w☎ ✹✒ log n n G 1 ✄ w☎ ✹✒ 1 n G 2 ✄ w☎ ✺✒ o✄ n ✕ 1 ☎ ✼✻ (8)
where A n ✏ E ✛ ✝ n ✢ . The leading term G 0 ✄ w☎ satisfies a non-linear integral equation. Indeed, in Section 3 we find that (cf. ( 49 on page 52))

e ✕ w G 0 ✄ w☎ ✏ ✶ 1 0 e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx (9) G 0 ✄ 0☎ ✏ 1; G✿ 0 ✄ 0☎ ✏ 0 ( 10
)
where

φ✄ x☎ ✏ x log x ✒ ❀✄ 1 ✑ x☎ log✄ 1 ✑ x☎ (11)
is the entropy of the Bernoulli(x) distribution. Furthermore, the correction terms G 1 ✄ ✧❁ ❂☎ and G 2 ✄ ☞❁ ☎ satisfy linear integral equations (cf. ( 50 on page 52)-( 51 on page 52)).

By using ( 8) in ( 7 on the page before) and asymptotically approximating the Cauchy integral we obtain

✘ ✩✙ ✧✛ ✝ n ✑ E✄ ✌✝ n ☎ ✏ ny✢ ❄❃ 1 n P✄ y☎ (12)
where

P✄ y☎ ✏ 1 2πi
✶ c✳ i∞ c✕ i∞ e ✕ yw G 0 ✄ w☎ dw✟ [START_REF] Knessl | Partition Asymptotics from Recursions[END_REF] or equivalently

G 0 ✄ w☎ ✏ ✶ ∞ ✕ ∞
e yw P✄ y☎ dy [START_REF] Knuth | The Art of Computer Programming. Sorting and Searching[END_REF] and c is a constant. Hence, G 0 ✄ w☎ is the moment generating function of the density P✄ y☎ . Now, we can summarize our main findings. We establish the results below under the assumptions that P✄ y☎ has certain forms as y ✥ ❆❅ ∞ (cf. Section 4).

Theorem 1 Consider the Quicksort algorithm that sorts n randomly selected keys.

(i) The limiting density P✄ y☎ satisfies

P✄ y ✒ 1☎ ✏ ✶ 1 0 ✶ ∞ ✕ ∞ P ✸ xt ✒ y ✑ 2φ✄ x☎ 2✄ 1 ✑ x☎ ✻ P ✸ ✑ ❇✄ 1 ✑ x☎ t ✒ y ✑ 2φ✄ x☎ 2x ✻ dtdx ( 15
)
and

✶ ∞ ✕ ∞ P✄ y☎ dy ✏ 1✟ ✶ ∞ ✕ ∞ yP✄ y☎ dy ✏ 0✠ (16)
(ii) The left tail of the distribution satisfies

✘ ✚✙ ✜✛ ✝ n ✑ E✄ ❈✝ n ☎ ✤❉ nz✢ ✯❃ 2 π 1 ❊ 2 log2 ✑ 1 exp ✸ ✑ α exp ✸ β ✑ z 2 ✑ log ✕ 1 2 ✻ ❋✻ (17)
for n ✥ ∞, and z ✏ z✄ n☎

that asymptotically becomes

w■ ✏ log ◆ y 2 ❖ ✒ loglog ◆ y 2 ❖ ✒ log log✄ y✎ 2☎ log✄ y✎ 2☎ ✄ 1 ✒ o✄ 1☎ ☞☎ (20)
for y ✥ ∞ (cf. ( 89 on page 60) in Section 4).

Finally, we relate our results for the tails to those of McDiarmid and Hayward [START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF]. These authors showed that

✘ ✚✙ ✜✛ ✌✪ ✝ n ✑ E✄ ✌✝ n ☎ ✪ ✮ εE✄ ❈✝ n ☎ ✜✢ ✏ exp ✛ ✑ 2ε logn✄ log log n ✑ log✄ 1✎ ε☎ ✹✒ O✄ log loglog n☎ ☞☎ ✢ ✧✟ (21)
which holds for n ✥ ∞ and ε in the range

1 log n ◗ ε ❉ 1✠ (22)
As pointed out in [START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF], this estimate is not very precise if, say, ε ✏ O✄ log log n✎ logn☎ .

¿From Theorem 1 we conclude that (since the right tail is much thicker than the left tail)

✘ ✚✙ ✧✛ ✌✪ ✝ n ✑ E✄ ❈✝ n ☎ ✪ ✬✫ ny✢ ✯❃ Ck ✄ y☎ e φ✽ y✾ ✟ y ✥ ∞ ( 23
)
where C is a constant and

φ✄
✑ yw■ ✯✒ ✶ w❏ 1 2e u u du ✏ y ❙ ✑ log ◆ y 2 ❖ ✑ log log ◆ y 2 ❖ ✒ 1 ✒ o✄ 1☎ ✵❚ (25)
✏ ✑ 2ε logn ✛ log log n ✑ log✄ 1✎ ε☎ ✹✒ log log✄ ε logn☎ ✯✑ 1✢ ❯✒ o✄ logn☎ which agrees precisely with the estimate [START_REF] Tan | Some Properties of a Limiting Distribution in Quicksort[END_REF], and also explicitly identifies the O✄ log loglog n☎ error term.

This suggests that ( 23) applies for y as large as 2 logn, though it cannot hold for y as large as n✎ 2 in view of the fact that

✘ ✩✙ ✵✛ ✝ n ✏ k✢ ✏ 0 for k ✮ ✰ n 2 ✲ .
An important open problem is obtaining an accurate numerical approximation to the constant C. This would likely involve the numerical solution of the integral equation for G 0 ✄ w☎ .

3 Analysis of the Generating Function for n ❱ ∞

We study ( 6 on page 45) asymptotically, for various ranges of n and u, namely: (i) u ✥ 1 with n fixed; (ii)

w ❲ n✄ u ✑ 1☎ fixed when n ✥ ∞ and u ✥ 1; (iii) w ✥ ❍❅ ∞; and (iv) u ◗ 1 or u ✮ 1.
We study these cases below.

A. CASE n FIXED AND u ✥ 1 First we consider the limit u ✥ 1 with n fixed. Then using the Taylor expansion

L n ✄ u☎ ✏ 1 ✒ A n ✄ u ✑ 1☎ ✹✒ B n ✄ u ✑ 1☎ 2 ✒ O✄ ✧✄ u ✑ 1☎ 3 ☎ (26) ✏ e A n ✽ u✕ 1✾ ✛ 1 ✒ ❀✄ B n ✑ 1 2 A 2 n ☎ ✴✄ u ✑ 1☎ 2 ✒ O✄ ☞✄ u ✑ 1☎ 3 ☎ ✢
we find from ( 6 on page 45) that A n ✏ L✿ n ✄ 1☎ and B n ✏ L✿ ✿ n ✄ 1☎ ☞✎ 2 satisfy the linear recurrence equations

A n✳ 1 ✏ n ✒ 1 n ✒ 1 n ∑ i✣ 0 ✛ A i ✒ A n✕ i ✢ ✏ n ✒ 2 n ✒ 1 n ∑ i✣ 0 A i ; A 0 ✏ 0✟ (27) B n✳ 1 ✏ ◆ n 2 ❖ ✒ 2n n ✒ 1 n ∑ i✣ 0 A i 1 n ✒ 1 ✒ n ∑ i✣ 0 ✛ 2B i ✒ A i A n✕ i ✢ ; B 0 ✏ 0✠ (28)
These are easily solved using either generating functions, or by multiplying ( 27) and ( 28) by n ✒ 1 and then differencing with respect to n. The final result is (cf. [START_REF] Knuth | The Art of Computer Programming. Sorting and Searching[END_REF][START_REF] Mahmoud | Evolution of Random Search Trees[END_REF])

A n ✏ 2✄ n ✒ 1☎ H n ✑ 4n (29) B n ✏ 2✄ n ✒ 1☎ 2 H 2 n ✑ ❳✄ 8n ✒ 2☎ (✄ n ✒ 1☎ H n ✒ n 2 ✄ 23n ✒ 17☎ ✯✑ 2✄ n ✒ 1☎ 2 H ✽ 2✾ n ✠ (30)
Here

H n ✏ 1 ✒ 1 2 ✒ 1 3 ✒ )❁ ☛❁ ✡❁ ✵✒ 1 n
is the harmonic number, and

H ✽ 2✾ n ✏ ∑ n k✣ 1 k ✕ 2
is the harmonic number of second order. In terms of A n and B n , the mean and variance of ✝ n are given by

E ✛ ✝ n ✢ ✏ A n ✏ 2✄ n ✒ 1☎ H n ✑ 4n (31) Var ✛ ✝ n ✢ ✏ L✿ ✿ n ✄ 1☎ ✹✒ L✿ n ✄ 1☎ ❬✑ ✛ L✿ n ✄ 1☎ ✢ 2 ✏ 2B n ✒ A n ✑ A 2 n (32) ✏ 7n 2 ✑ 2✄ n ✒ 1☎ H n ✒ 13n ✑ 4✄ n ✒ 1☎ 2 H ✽ 2✾ n ✠
Asymptotically, for n ✥ ∞, we obtain

A n ✏ 2n logn ✒ ❀✄ 2γ ✑ 4☎ n ✒ 2 logn ✒ 2γ ✒ 1 ✒ 5 6 n ✕ 1 ✒ O✄ n ✕ 2 ☎ (33) B n ✑ 1 2 A 2 n ❲ C n ✏ ✸ 7 2 ✑ π 2 3 ✻ n 2 ✑ 2n logn ✒ n ✸ 21 2 ✑ 2γ ✑ 2 3 π 2 ✻ ✒ o✄ n☎ ✴✠ (34)
These expressions will be used in order to asymptotically match the expansion for u ✥ 1 and n fixed, to those that will apply for other ranges of n and u. Since it is well known that the lth moment of

✝ n ✑ E ✛ ✝ n ✢ is of order
O✄ n l ☎ as n ✥ ∞, all terms in the series ( 26) will be of comparable magnitude when (roughly)

u ✏ 1 ✏ O✄ n ✕ 1 ☎
and n ✥ ∞. If we view ( 26) as an asymptotic series for n fixed and u ✥ 1, then it will cease to be valid when

n ✏ O✄ ☞✄ u ✑ 1☎ ✕ 1 ☎
, which motivates the analysis that follows.

B. CASE w ❲ n✄ u ✑ 1☎ FIXED WHEN n ✥ ∞ AND u ✥ 1
Next we consider the limit u ✥ 1, n ✥ ∞ with w ❲ n✄ u ✑ 1☎ held fixed. This scaling is necessary to obtain a non-trivial limiting problem. We define G✄ ✧❁ ❂☎ by

L n ✄ u☎ ✏ exp ✄ A n w✎ n☎ G✄ w; n☎ ✏ e A n ✽ u✕ 1✾ G✄ n✄ u ✑ 1☎ ; n☎ ✴✠ (35)
With this change of variables, we rewrite ( 6 on page 45) as

exp ✄ A n✳ 1 w✎ n☎ (❭ G✄ w ✒ w n ; n ✒ 1☎ ✏ ✄ 1 ✒ w✎ n☎ n n ✒ 1 (36) ❭ ) n✕ m ∑ i✣ m exp ◆ A i w n ✒ A n✕ i w n ❖ G ✸ w i n ; i✻ G ✸ w ✸ 1 ✑ i n ✻ ; n ✑ i✻ ✒ 2 m✕ 1 ∑ i✣ 0 exp ◆ A i w n ✒ A n✕ i w n ❖ ✸ 1 ✒ C i w 2 n 2 ✒ ❀❁ ☛❁ ☛❁ ✻ G ✸ w ✸ 1 ✑ i n ✻ ; n ✑ i✻ ❵❴ where C i ✏ B i ✑ 1 2 A 2 i .
Here we have broken up the sum in ( 6 on page 45) into the three ranges 0

❉ i ❉ m ✑ 1, m ❉ i ❉ n ✑ m, and n ✑ m ✒ 1 ❉ i ❉ n,
and used the symmetry ✄ i ✥ n ✑ i☎ of the summand. We expect ( 26 on the preceding page) to also be valid for large values of n, as long as n✄ u ✑ 1☎ ✩✥ 0 as n ✥ ∞. Thus, for 0 ❉ i ❉ m ✑ 1 we replaced L i ✄ u☎ in the sum by the approximation ( 26 on the page before). The integer m may be chosen arbitrarily, since the right side of ( 36) must ultimately be independent of m. For now we assume that m ✥ ∞ but m✎ n ✥ 0 as n ✥ ∞. For n large we have (cf. ( 33 on the page before))

A i n ✒ A n✕ i n ✑ A n✳ 1 n ✏ 2 ✸ i n log ✸ i n ✻ |✒ ❛✸ 1 ✑ i n ✻ log ✸ 1 ✑ i n ✻ ❜✻ (37) ✒ 2 n log ✸ i n ✸ 1 ✑ i n ✻ ❋✻ |✒ 3 n ✒ o✄ n ✕ 1 ☎
with which we rewrite ( 36) as

✄ n ✒ 1☎ e ✕ w ❭ ✸ 1 ✒ w 2 2n ✒ O✄ n ✕ 2 ☎ ✻ G ◆ w ✒ w n ; n ✒ 1 ❖ (38) ✏ n✕ m ∑ i✣ m e 2φ✽ i❝ n✾ w ✸ 1 ✒ w ✸ 2 n ψ ✸ i n ✻ ✓✒ 3 n ✻ ✓✒ O✄ n ✕ 2 ✟ n ✕ 1 m ✕ 1 ☎ ✼✻ ❭ G ✸ w i n ; i✻ G ✸ w ✸ 1 ✑ i n ✻ ; n ✑ i✻ ✒ 2 m✕ 1 ∑ i✣ 0 ◆ 1 ✒ A i w n ✒ O ✰ B i n ✕ 2 ✲ ❖ exp ◆ ❞✄ A n✕ i ✑ A n✳ 1 ☎ w n ❖ ❭ G ✸ w ✸ 1 ✑ i n ✻ ; n ✑ i✻ where φ✄ x☎ ✏ x log x ✒ ❀✄ 1 ✑ x☎ log✄ 1 ✑ x☎ (39) ψ✄ x☎ ✏ log ✛ x✄ 1 ✑ x☎ ✢ ✧✠
We now evaluate the two sums in ( 38) asymptotically and show that when the two results are added, the dependence on m disappears.

From ( 35 on the preceding page) and the identity L✿ n ✄ 1☎ ✏ A n ✏ E ✛ ✝ n ✢ we find that for all n G✄ 0; n☎ ✏ 1 and G✿ ✄ 0; n☎ ✏ 0✠ (40)

We assume that for n ✥ ∞, G✄ w; n☎ has an asymptotic expansion of the form

G✄ w; n☎ ✏ G 0 ✄ w☎ ✹✒ a 1 ✄ n☎ G 1 ✄ w☎ ✺✒ a 2 ✄ n☎ G 2 ✄ w☎ ✹✒ ❀❁ ☛❁ ☛❁ (41)
where a j ✄ n☎ is an asymptotic sequence as n ✥ ∞, i.e., a j✳ 1 ✄ n☎ ✧✎ a j ✄ n☎ ❬✥ 0 as n ✥ ∞. The appropriate sequence is determined by balancing terms in ( 38 on the page before). This will eventually yield

a 1 ✄ n☎ ✏ log n n ; a 2 ✄ n☎ ✏ 1 n (42)
so we use this form from the beginning. Note that G 0 ✄ w☎ is the moment generating function of the limiting density of ✄ ❈✝ n ✑ E ✛ ✝ n ✢ ✦☎ ✧✎ n, which is discussed in [START_REF] Rösler | A Limit Theorem for Quicksort[END_REF]. The conditions ( 40) imply that

G 0 ✄ 0☎ ✏ 1; G 1 ✄ 0☎ ✏ G 2 ✄ 0☎ ✏ ❁ ☛❁ ✡❁ ✏ 0 (43) G✿ 0 ✄ 0☎ ✏ G✿ 1 ✄ 0☎ ✏ G✿ 2 ✄ 0☎ ✏ ❁ ✡❁ ☛❁ ✏ 0✠ (44)
We consider the first sum in ( 38 on the page before), which we denote by S 1 ✏ S 1 ✄ n; m☎ . Using ( 41), ( 42)

and the Euler-MacLaurin formula we obtain

S 1 ✏ n✕ m ∑ i✣ m e 2φ✽ i❝ n✾ w ❡ 1 ✒ w ✸ 2 n ψ ✸ i n ✻ ✒ 3 n ✻ ✒ o✄ n ✕ 1 ☎ ❢ (45) ❭ ❣ G 0 ✸ w i n ✻ G 0 ✸ w ✸ 1 ✑ i n ✻ ❜✻ ✒ G 0 ✸ w i n ✻ log✄ n ✑ i☎ n ✑ i G 1 ✸ w ✸ 1 ✑ i n ✻ ❜✻ ✒ G 0 ✸ w ✸ 1 ✑ i n ✻ ❋✻ log i i G 1 ✸ w i n ✻ ✔✒ G 0 ✸ w i n ✻ 1 n ✑ i G 2 ✸ w ✸ 1 ✑ i n ✻ ❋✻ ✒ G 0 ✸ w ✸ 1 ✑ i n ✻ ❋✻ 1 i G 2 ✸ w i n ✻ ✒ ❀❁ ☛❁ ✡❁ ✐❤ ✏ n ✶ 1✕ m❝ n m❝ n e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✒ e 2φ✽ m❝ n✾ w G 0 ◆ w m n ❖ G 0 ◆ w ◆ 1 ✑ m n ❖ •❖ ✒ w ✶ 1✕ m❝ n m❝ n ✛ 2ψ✄ x☎ ✹✒ 3✢ e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✒ ✶ 1✕ m❝ n m❝ n log n ✒ log✄ 1 ✑ x☎ 1 ✑ x e 2φ✽ x✾ w G 0 ✄ wx☎ G 1 ✄ w ✑ wx☎ dx ✒ ✶ 1✕ m❝ n m❝ n log n ✒ logx x e 2φ✽ x✾ w G 0 ✄ w ✑ wx☎ G 1 ✄ wx☎ dx ✒ ✶ 1✕ m❝ n m❝ n 1 1 ✑ x e 2φ✽ x✾ w G 0 ✄ wx☎ G 2 ✄ w ✑ wx☎ dx ✒ ✶ 1✕ m❝ n m❝ n 1 x e 2φ✽ x✾ w G 2 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✒ o✄ 1☎ (✠
We note that all the integrals remain finite as m✎ n ✥ 0, in view of ( 43) and ( 44). However, if we were to consider higher order terms in the expansion ( 41), which would involve terms of order n ✕ 2 and n ✕ 3 , then the corresponding higher order terms in ( 45 on the page before) would involve integrands not integrable over [0,1]. It then becomes essential that we integrate only over the range ✛ m✎ n✟ 1 ✑ m✎ n✢ and consider the contribution from the second sum S 2 ✏ S 2 ✄ n; m☎ in ( 38 on page 49). We can further simplify S 1 by evaluating each term in the limit m✎ n ✥ 0, which we assumed to be true. We have

T 11 ❲ ✶ 1✕ m❝ n m❝ n e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✏ ✶ 1 0 e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✑ 2 ✶ m❝ n 0 e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✏ ✶ 1 0 e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✑ 2 ✶ m❝ n 0 ✛ 1 ✒ 2w✄ x logx ✑ x☎ ✹✒ O✄ x 2 ☎ ✢ ❭ ✛ 1 ✒ O✄ x 2 ☎ ✢ ✛ G 0 ✄ w☎ ❬✑ wxG✿ 0 ✄ w☎ ✹✒ O✄ x 2 ☎ ✢ dx ❃ ✶ 1 0 e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✑ 2G 0 ✄ w☎ ❡ m n ✒ 2w ✸ m 2 2n 2 log ◆ m n ❖ ✑ 3 4 m 2 n 2 ✻ ❢ ✒ m 2 n 2 wG✿ 0 ✄ w☎ ; T 12 ❲ 2 ✶ 1✕ m❝ n m❝ n logn ✒ log x x e 2φ✽ x✾ w G 0 ✄ w ✑ wx☎ G 1 ✄ wx☎ dx ❃ 2 ✶ 1 0 log n ✒ logx x e 2φ✽ x✾ w G 0 ✄ w ✑ wx☎ G 1 ✄ wx☎ dx ✑ 2 m n ✄ log n☎ G 1 ✄ w☎ ;
and

T 13 ❲ 2 ✶ 1✕ m❝ n m❝ n e 2φ✽ x✾ w 1 x G 2 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ❃ 2 ✶ 1 0 e 2φ✽ x✾ w 1 x G 2 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✑ 2 m n G 2 ✄ w☎ (✠
Thus S 1 simplifies to

S 1 ✏ n ✶ 1 0 e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✒ G 0 ✄ w☎ (46) ✒ w ✶ 1 0 ✛ 2ψ✄ x☎ ✺✒ 3✢ e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✒ 2 ✶ 1 0 log n ✒ logx x e 2φ✽ x✾ w G 0 ✄ w ✑ wx☎ G 1 ✄ wx☎ dx ✒ 2 ✶ 1 0 1 x e 2φ✽ x✾ w G 0 ✄ w ✑ wx☎ G 2 ✄ wx☎ dx ✒ ❣ ✑ 2mG 0 ✄ w☎ ✹✒ w m 2 n G✿ 0 ✄ w☎ ❬✑ 2m 2 n wG 0 ✄ w☎ ❡ log ◆ m n ❖ ✑ 3 2 ❢ ✑ 2 m n ✄ log n☎ G 1 ✄ w☎ ✯✑ 2 m n G 2 ✄ w☎ ✴• ❦✒ ❀❁ ☛❁ ☛❁
where we have grouped the terms involving m inside the "✞ •✍ ". The error term in ( 46) approaches zero as

n ✥ ∞.
Now we consider the second sum S 2 in ( 38 on page 49).

Using A✄ n ✑ i☎ ❑✑ A✄ n ✒ 1☎ ✚❃ ♠✑ ❇✄ i ✒ 1☎ A✿ ✄ n☎ ✚❃ ✑ ❇✄ i ✒ 1☎ (✄ 2 logn ✒ 2γ ✑ 2☎ , we obtain S 2 ❃ 2 m✕ 1 ∑ i✣ 0 ❙ 1 ✒ A i w n ❚ ❇❙ 1 ✑ w n ✄ i ✒ 1☎ ✴✄ 2 logn ✒ 2γ ✑ 2☎ ❚ (47) ❭ ❡ G 0 ✄ w☎ ✯✑ i n wG✿ 0 ✄ w☎ ✹✒ logn n G 1 ✄ w☎ ✹✒ 1 n G 2 ✄ w☎ ❢ ❃ 2G 0 ✄ w☎ ❋♥ m ✒ w n m✕ 1 ∑ i✣ 0 ✄ A i ✑ 2✄ log n☎ ✴✄ i ✒ 1☎ ❬✑ ❳✄ 2γ ✑ 2☎ (✄ i ✒ 1☎ ☞☎ ♣♦ ✒ 2 m n ✄ log n☎ G 1 ✄ w☎ ✹✒ G 2 ✄ w☎ ✢ ❞✑ m 2 n wG✿ 0 ✄ w☎ ❃ 2G 0 ✄ w☎ m ✒ 2 m 2 n wG 0 ✄ w☎ ❡ log ◆ m n ❖ ✑ 3 2 ❢ ✑ m 2 n wG✿ 0 ✄ w☎ ✒ 2 m n ✛ ✄ logn☎ G 1 ✄ w☎ ✺✒ G 2 ✄ w☎ ✜✢ ✧✠
Here we have used 2

∑ n✕ 1 i✣ 0 A i ✏ n ✛ A n ✑ |✄ n ✑ 1☎ ✜✢ .
Upon adding ( 46 on the page before) to ( 47 on the preceding page), we see that all the terms involving m cancel, and that the leading three terms in the expansion of S 1 ✒ S 2 (i.e., the right side of ( 38 on page 49)) are of order O✄ n☎ , O✄ log n☎ and O✄ 1☎ , respectively. Using ( 41 on page 50) and ( 42 on page 50), the left side of ( 38 on page 49) becomes

✄ n ✒ 1☎ e ✕ w ✸ 1 ✒ w 2 2n ✻ q✸ G 0 ✄ w☎ ✹✒ log n n G 1 ✄ w☎ ✹✒ 1 n ✄ wG✿ 0 ✄ w☎ ✹✒ G 2 ✄ w☎ ☞☎ ✺✒ o✄ n ✕ 1 ☎ ✼✻ |✠ (48)
Thus, comparing the above to S 1 ✒ S 2 , we find that

e ✕ w G 0 ✄ w☎ ✏ ✶ 1 0 e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx✟ (49) e ✕ w G 1 ✄ w☎ ✏ 2 ✶ 1 0 1 x e 2φ✽ x✾ w G 0 ✄ w ✑ wx☎ G 1 ✄ wx☎ dx✟ (50) e ✕ w ✛ G 0 ✄ w☎ ✹✒ 1 2 w 2 G 0 ✄ w☎ r✒ wG✿ 0 ✄ w☎ ✹✒ G 2 ✄ w☎ ✢ ✏ G 0 ✄ w☎ (51) ✒ w ✶ 1 0 ✛ 2ψ✄ x☎ ✹✒ 3✢ e 2φ✽ x✾ w G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ dx ✒ 2 ✶ 1 0 log x x e 2φ✽ x✾ w G 0 ✄ w ✑ wx☎ G 1 ✄ wx☎ dx ✒ 2 ✶ 1 0 1 x e 2φ✽ x✾ w G 0 ✄ w ✑ wx☎ G 2 ✄ wx☎ dx✠
Equations ( 49)-( 51), along with ( 39 on page 49), ( 43 on page 50) and ( 44 on page 50) are integral equations for the first three terms in the series ( 41 on page 50). Below we discuss some aspects of the solutions to these problems. The leading order equation ( 49) was previously obtained in [START_REF] Eddy | How Many Comparisons Does Quicksort Use[END_REF], using more probabilistic arguments.

We observe that the solution to ( 49) is not unique: G 0 ✄ w☎ ✏ 0 is one solution and if G 0 ✄ w☎ is any solution, then so is e cw G 0 ✄ w☎ for any constant c. We can construct the solution as a Taylor series:

G 0 ✄ w☎ ✏ 1 ✒ ∞ ∑ j✣ 1 g j w j ✠ (52)
This eliminates the trivial solution G 0 ✄ w☎ ✏ 0 and satisfies the normalization G 0 ✄ 0☎ ✏ 1. Using ( 52) in ( 49) and noting that ✶ 1 0 ✛ 2φ ✄ x☎ ✯✒ 1✢ dx ✏ 0, we see that g 1 remains arbitrary, and then we can easily calculate g j for j ✫ 2 in terms of g 1 . But, ( 44 on page 50) forces g 1 ✏ 0 and then all the Taylor coefficients in ( 52) are uniquely determined. They may be evaluated from the recursion

✸ 1 ✑ 2 n ✒ 1 ✻ g n ✏ n✕ 1 ∑ i✣ 1 B✄ i ✟ n ✑ i ✟ 0☎ g i g n✕ i (53) ✒ n✕ 1 ∑ s ✣ 0 s ∑ i✣ 0 B✄ i✟ ★t ✩✑ 1✟ n ✑ ✉t ✴☎ g i g s ✕ i
for n ✫ 2 where

B✄ i✟ j ✟ k☎ ✏ ✶ 1 0 x i ✄ 1 ✑ x☎ j 1 k! ✛ 2φ✄ x☎ ✹✒ 1✢ k dx✠ (54)
In particular,

g 2 ✏ 7 2 ✑ π 2
3 . Next we consider the equations ( 50 on the page before) and ( 51 on the preceding page) for the correction terms G 1 and G 2 . These are linear, Fredholm integral equations of the second kind. Their solutions may also be constructed as Taylor series in w. In view of ( 43 on page 50) and ( 44 on page 50) we must have

G 0 ✄ w☎ ✈❃ 1 ✒ α 0 w 2 ; G 1 ✄ w☎ ✈❃ β 0 w 2 ; G 2 ✄ w☎ ✈❃ γ 0 w 2 ✟
as w ✥ 0, where we have already computed α 0 ✏ g 2 . Given β 0 , we can easily compute the higher order

Taylor coefficients of G 1 ✄ w☎ from ( 50 on the page before), in terms of the (now uniquely determined) Taylor coefficients of G 0 ✄ w☎ . However, the constant β 0 cannot be determined solely from ( 50 on the preceding page), ( 43 on page 50) and ( 44 on page 50). To fix β 0 we use the principle of asymptotic matching.

We require that expansions ( 26 on page 48) and ( 35 on page 49) (with ( 41 on page 50)) agree in some intermediate limit, where u ✥ 1, n ✥ ∞ and n✄ u ✑ 1☎ ✩✥ 0. Then the behavior of ( 26 on page 48) as n ✥ ∞ must agree with the behavior of ( 35 on page 49) (with ( 41 on page 50)) as w ✥ 0. Writing ( 26 on page 48)

as

L n ✄ u☎ ✏ exp ✛ A n ✄ u ✑ 1☎ ✜✢ ✛ 1 ✒ C n ✄ u ✑ 1☎ 2 ✒ O✄ ✧✄ u ✑ 1☎ 3 ☎ ✜✢ the matching condition becomes 1 ✒ C n ✄ u ✑ 1☎ 2 ✒ ❀❁ ☛❁ ☛❁ ☛✇ ✇ ✇ ✇ n① ∞ ❃ G 0 ✄ w☎ ✹✒ log n n G 1 ✄ w☎ ✺✒ 1 n G 2 ✄ w☎ ✹✒ ❀❁ ☛❁ ☛❁ ♣✇ ✇ ✇ ✇ w① 0 ✠ (55)
Setting u ✑ 1 ✏ w✎ n, using ( 34 on page 48), and noting that the right side of ( 55) is 1 ✒ n ✕ 2 ✛ α 0 n 2 ✒ β 0 ✄ log n☎ n ✒ γ 0 n✢ w 2 , we obtain

β 0 ✏ ✑ 2; γ 0 ✏ 21 2 ✑ 2γ ✑ 2 3 π 2 ✠ ( 56
)
Now consider equation ( 51 on the preceding page) in the limit w ✥ 0. We write ( 50 on the page before) abstractly as T G 1 ✏ 0 where T is the linear integral operator in ( 50 on the preceding page). Then ( 51 on the page before) may be written as

T G 2 ✏ f ✄ G 0 ✟ G 1 ☎
where f is a known function of the first two terms.

Since T G 1 ✏ 0 has a non-zero solution (made unique by the condition β 0 ✏ ✑ 2☎ , we expect that T G 2 ✏ f will have a solution only if a solvability condition is met. To obtain this solvability condition we expand ( 51 on the preceding page) as w ✥ 0. Obviously ( 51 on the page before) is satisfied as w ✥ 0, and, since

✶ 1 0 ✛ 2ψ✄ x☎ ✯✒ 3✢ dx ✏ ✑ 1, ( 51
on the preceding page) also holds to order O✄ w☎ . Comparing O✄ w 2 ☎ terms in ( 51 on page 52) using ( 55 on the preceding page) gives

1 ✒ 2α 0 ✏ ✑ 1 2 β 0 ✑ 3 ✒ 4 ✶ 1 0 φ✄ x☎ ψ✄ x☎ dx (57) ✏ ✑ 1 2 β 0 ✑ 3 ✒ 4 ✸ 5 2 ✑ π 2 6 ✻
and thus β 0 ✏ ✑ 2, which regains the result we obtained by using asymptotic matching. Note, however, that this argument required that we derive the equation satisfied by G 2 ✄ w☎ in order to uniquely specify the previous term G 1 ✄ w☎ . In contrast, asymptotic matching made no use of ( 51 on page 52). Presumably, by deriving the equation for the next term G 3 ✄ w☎ in ( 41 on page 50) and examining its solvability condition, we would have

an alternate way of computing γ 0 . Given β 0 and γ 0 , we can easily obtain the Taylor expansions of G 1 and G 2 from ( 50 on page 52) and ( 51 on page 52).

To summarize, we have obtained the expansion ( 35 on page 49), with ( 41 on page 50)-( 44 on page 50), ( 49 on page 52)-( 51 on page 52) and ( 56 on the preceding page), for the scaling n ✥ ∞, u ✥ 1 with n✄ u ✑ 1☎ ✏ w fixed. We have not been able to explicitly solve these integral equations. However, we can derive some approximate formulas in the limits w ✥ ❆❅ ∞, and these may be used to obtain approximations to the tail probabilities of the Quicksort distribution.

C. CASE w ❲ n✄ u ✑ 1☎ ❑✥ ❆❅ ∞
We shall only examine the leading term G 0 , and we first consider the limit w ✥ ②✑ ∞. As w ✥ ②✑ ∞, the "kernel" exp ✛ 2wφ✄ x☎ ✢ in ( 49 on page 52) is sharply concentrated near x ✏ 1✎ 2, and behaves as a multiple of

✪ w ✪ ✕ 1❝ 2 δ✄ x ✑ 1✎ 2☎ .
Thus we treat ( 49 on page 52) as a Laplace type integral (cf. [START_REF] Henrici | Applied and Computational Complex Analysis[END_REF]).

Assuming that G 0 ✄ w☎ has a weak (e.g., algebraic) dependence on w, we approximate the right side of ( 49 on page 52) by Laplace's method, which (to leading order) yields the functional equation

e ✕ w G 0 ✄ w☎ •③ ❙ G 0 ◆ w 2 ❖ ❚ 2 ④ π ✑ 4w e ✕ 2wlog 2 ✟ w ✥ ⑤✑ ∞✠ (58)
But, if G 0 varies weakly with w, then the exponential orders of magnitude of the right and left sides of ( 58) do not agree. In order to get agreement, we would need G 0 ✄ w☎ to vary much more rapidly as w ✥ ⑤✑ ∞, of the order exp ✛ O✄ w log✄ ✼✑ w☎ ✧☎ ✜✢ . But this then contradicts the assumption used to obtain ( 58). Therefore, we return to ( 49 on page 52) and allow more rapid variation of G 0 . Specifically, we assume that as

w ✥ ❍✑ ∞ G 0 ✄ w☎ •❃ e α 1 wlog✽ ⑥✕ w✾ ⑥✳ βw ✄ ✵✑ w☎ γ 1 δ 1 ✠ (59)
Using ( 59) in ( 49 on page 52) yields

e ✕ w ✄ ✼✑ w☎ γ 1 δ 1 ❃ δ 2 1 ✄ ✵✑ w☎ 2γ 1 ✶ 1 0 x γ 1 ✄ 1 ✑ x☎ γ 1 e φ✽ x✾ ⑥✽ 2✳ α 1 ✾ w dx✠ (60)
If 2 ✒ α 1 ✮ 0 we again use Laplace's method to approximate the integral, thus obtaining

e ✕ w ❃ δ 1 ✸ 1 4 ✻ γ 1 ✄ ✼✑ w☎ γ 1 ✕ 1 2 ④ π 2✄ 2 ✒ α 1 ☎ e ✕ ▲✽ α 1 ✳ 2✾ ⑥✽ log 2✾ w ✠ (61)
Hence, we must have

α 1 ✏ ✑ 2 ✒ 1 log 2 ; γ 1 ✏ 1 2 ; δ 1 ✏ 2 ❊ 2 ❊ π log2 (62)
and β remains arbitrary. Since the solution of ( 49 on page 52) is not unique, we can never determine β using only this integral equation. We thus have

G 0 ✄ w☎ •❃ 2 ❊ 2 ❊ π log 2 ❊ ✑ w exp ❡ βw ✒ ❛✸ 1 log 2 ✑ 2✻ w log✄ ✼✑ w☎ ❢ ✟ w ✥ ❍✑ ∞✠ (63)
By computing higher order terms in the expansion of the integral in ( 49 on page 52), we can extend ( 59 on the page before) and ( 63) to an asymptotic series in powers of w ✕ 1 . To fix β, we must use the condition in ( 43 on page 50) and ( 44 on page 50), and this probably requires a numerical solution of ( 49 on page 52).

By comparing this numerical solution to ( 49 on page 52) for large, negative w to ( 63), we can obtain an approximation to β.

Now consider ( 49 on page 52) in the limit w ✥ ⑤✒ ∞. Then exp ✛ 2wφ✄ x☎ ✜✢ is concentrated near the endpoints

x ✏ 0 and x ✏ 1. We assume that G 0 ✄ w☎ has an asymptotic expansion in the WKB form [START_REF] Froman | JWKB Approximation[END_REF]:

G 0 ✄ w☎ •❃ K ✄ w☎ exp ✛ Ψ✄ w☎ ✜✢ ✼✠ (64)
The major contribution to the integral will come from where 1 ✑ x ✏ o✄ w ✕ 1 ☎ . Thus we use the Taylor series

G 0 ✄ wx☎ ✏ 1 ✒ O✄ ☞✄ wx☎ 2 ☎
, and ( 64) to approximate G 0 ✄ w ✑ wx☎ . This yields Note that ( 66) is consistent with wΨ✿ ✄ w☎ ❬✥ ∞ and shows that G 0 grows very rapidly (as a double exponential) as w ✥ ❆✒ ∞. We note that the asymptotic expansion of the integrand in ( 65) is only valid for x ✏ o✄ w ✕ 1 ☎ , but the major contribution to the integral turns out to be from x ✏ O✄ e ✕ w ☎ , which is certainly o✄ w ✕ 1 ☎ . To obtain the next term K ✄ w☎ , we refine the expansion ( 65) to

K ✄ w☎ e Ψ✽ w✾ ⑥✕ w ❃ 2 ✶ 1❝ 2 0 ✛ 1 ✒ 2wφ✄ x☎ ✹✒ O✄ w 2 φ 2 ✄ x☎ ✧☎ ✜✢ ✛ 1 ✒ O✄ ☞✄ wx☎ 2 ☎ ✢ K ✄ w ✑ wx☎
K ✄ w☎ e Ψ✽ w✾ ✭✕ w ✏ 2 ✶ 1❝ 2 0 ✛ 1 ✒ 2wφ ✄ x☎ ✹✒ O✄ w 2 φ 2 ✄ x☎ ☞☎ ✢ ✛ 1 ✒ O✄ ✧✄ wx☎ 2 ☎ ✜✢ (67) ❭ ✛ K ✄ w☎ ❬✑ wxK✿ ✄ w☎ ✹✒ O✄ w 2 x 2 K✿ ✿ ✄ w☎ ✧☎ ✜✢ ❭ exp ✸ Ψ✄ w☎ ❬✑ wΨ✿ ✄ w☎ x ✒ 1 2 w 2 Ψ✿ ✿ ✄ w☎ x 2 ✒ O✄ w 3 Ψ✿ ✿ ✿ ✄ w☎ x 3 ☎ ✻ dx✠
For x small we have 2wφ✄ x☎ ✏ 2w✄ x logx ✑ x ✒ O✄ x 2 ☎ ✜✢ and if x ✏ O✄ e ✕ w ☎ then w 2 x 2 ✏ O✄ w 2 e ✕ 2w ☎ and

w 3 x 3 Ψ✿ ✿ ✿ ✄ w☎ ✏ O✄ w 2 e ✕ 2w ☎ . Setting x ✏ ηe ✕ w ✎ 2 we have exp ❡ 1 2 x 2 w 2 Ψ✿ ✿ ✄ w☎ ❢ ✏ 1 ✒ 1 4 η 2 ✄ w ✑ 1☎ e ✕ w ✛ 1 ✒ o✄ 1☎ ✢ ✧✠
Thus ( 67 on the preceding page) becomes, after cancelling the common factor exp ✛ Ψ✄ w☎ ✜✢ ,

e ✕ w K ✄ w☎ ✏ e ✕ w ❣ K ✄ w☎ ✶ ∞ 0 e ✕ η dη (68) ✒ K ✄ w☎ ✶ ∞ 0 e ✕ η ❡ 1 4 e ✕ w η 2 ✄ w ✑ 1☎ ✹✒ we ✕ w η ◆ ✷✑ w ✒ log ◆ η 2 ❖ ✑ 1 ❖ ❢ dη ✑ K✿ ✄ w☎ we ✕ w ✶ ∞ 0 1 2 ηe ✕ η dη ✒ O✄ e ✕ 2w K ✄ w☎ w α ☎ ❯❤
for some constant α. To leading order ( 68) is obviously satisfied and then collecting terms that are of order O✄ e ✕ w ☎ as w ✥ ∞ (modulo some algebraic factors) we obtain the following differential equation

w 2 K✿ ✄ w☎ K ✄ w☎ ✏ ✶ ∞ 0 ❡ 1 4 η 2 e ✕ η ✄ w ✑ 1☎ ✹✒ wηe ✕ η ◆ ✷✑ w ✒ log ◆ η 2 ❖ ✑ 1 ❖ ❢ dη (69)
and thus

K✿ ✄ w☎ K ✄ w☎ ✏ ✑ 2w ✒ 1 ✑ 2γ ✑ 2 log2 ✑ 1 w ✠ (70)
Solving ( 70) and using the result along with ( 66 on the page before), we have

G 0 ✄ w☎ •❃ C■ exp ✸ ✶ w 1 2e u u due ✕ w 2 ✳ ❑✽ 1✕ 2γ✕ 2log 2✾ w ✻ 1 w ✠ ( 71
)
Here C■ is an undetermined constant and we have chosen the lower limit on the integral in ( 66 on the preceding page) as one. An alternate choice would only change C■ , which we cannot determine using only the integral equation ( 49 on page 52).

Our analysis shows that as w ✥ ⑧✒ ∞, the nonlinear equation ( 49 on page 52) may be approximated by a linear one. To fix C■ it would seem that we will again need an accurate numerical solution to ( 49 on page 52).

We have thus obtained formal asymptotic results as w ✥ ❍❅ ∞ for the solution to ( 49 on page 52). Using our procedure we can derive full asymptotic series in these limits, but the constants β and C■ will remain undetermined.

In Section 4, we will use our results for G 0 ✄ w☎ to obtain asymptotic expansions for the limiting density P✄ y☎ as y ✥ ⑨❅ ∞.

D. CASE u ◗ 1 OR u ✮ 1
We next study ( 6 on page 45) for n ✥ ∞ but for fixed u ✮ 1 or u ◗ 1. First we assume that 0 ◗ u ◗ 1. Note that u appears in ( 6 on page 45) only as a parameter. We assume an expansion of the form

L n ✄ u☎ ✚❃ e A✽ u✾ nlog n✳ B✽ u✾ n n C✽ u✾ D✄ u☎ (72)
for n ✥ ∞ and u ◗

1. The major contribution to the sum in ( 6 on page 45) will come from the midpoint i ③ n✎ 2. Using ( 72) in ( 6 on page 45) and noting that

✄ n ✒ 1☎ log✄ n ✒ 1☎ = n logn ✒ log n ✒ 1 ✒ O✄ n ✕ 1 ☎ and i logi ✒ ❀✄ n ✑ i☎ log✄ n ✑ i☎ = n logn ✑ n log2 ✒ 2n ✕ 1 ✄ i ✑ n✎ 2☎ 2 + O✄ n ✕ 2 ✄ i ✑ n✎ 2☎ 3 ☎ , we obtain e Anlog n e A log n e A✳ B e Bn n C D ❃ D 2 ◆ n 2 ❖ 2C e Bn u n n e Anlog n e ✕ Anlog 2 ❭ n ∑ i✣ 0 exp ✰ ✑ ❇✄ 2✎ n☎ ✪ A ✪ ✄ i ✑ n✎ 2☎ 2 ✲ for A ◗ 0.
The sum is asymptotically equal to ⑩ nπ✎ 2 ✪ A ✪ , and thus

A ✏ log u log 2 ✟ C ✏ 1 2 ✒ log u log 2 ✟ (73) D ✏ 2e B u 2 ❶ ✑ 2 logu π log2 exp ✸ log u log 2 ✻
and B ✏ B✄ u☎ remains arbitrary.

Some further information may be obtained by asymptotically matching ( 72 on the preceding page) to the expansion valid for n ✥ ∞, u ✑ 1 ✏ O✄ n ✕ 1 ☎ . In the intermediate limit where u ✥ 1, n ✥ ∞ and n✄ u ✑ 1☎ ✹✥ ❷✑ ∞, we use ( 63 on page 55) with w ✏ n✄ u ✑ 1☎ to get

e A n ✽ u✕ 1✾ G 0 ✄ w☎ •❃ e 2✽ u✕ 1✾ nlog n e ✽ 2γ✕ 4✾ n✽ u✕ 1✾ 2 ❊ 2 ❊ π log2 e βn✽ u✕ 1✾ ❭ ⑩ n✄ 1 ✑ u☎ exp ❣ ✸ 1 log 2 ✑ 2✻ n✄ u ✑ 1☎ ✛ logn ✒ log✄ 1 ✑ u☎ ✜✢ ☞❤ •✠
As u ❸ 1, ( 72 on the preceding page) and ( 73) yields

exp ❣ ❡ u ✑ 1 log 2 ✒ O✄ ✧✄ u ✑ 1☎ 2 ☎ ❢ n logn ✒ Bn❤ n 1❝ 2 2e B ❶ 2✄ 1 ✑ u☎ π log2
and these two expressions agree provided that as u ❸ 1

B✄ u☎ ✚❃ ✸ 1 log 2 ✑ 2✻ ✄ u ✑ 1☎ log✄ 1 ✑ u☎ ✺✒ ❹✄ 2γ ✑ 4 ✒ β☎ (✄ u ✑ 1☎ ✴✟ u ✥ 1✠ (74)
This relates the behavior of B✄ u☎ as u ❸ 1 to the constant β in ( 63 on page 55). Now consider n ✥ ∞ with u ✮ 1. The dominant contribution in the sum ( 6 on page 45) now comes from the terms with i ✏ 0 and i ✏ n. Thus

L n✳ 1 ✄ u☎ ✏ u n n ✒ 1 ✛ 2L 0 ✄ u☎ L n ✄ u☎ ✹✒ 2L 1 ✄ u☎ L n✕ 1 ✄ u☎ ✹✒ ❀❁ ☛❁ ☛❁ ✢ ❄❃ 2u n n ✒ 1 L n ✄ u☎ (75)
and

L n ✄ u☎ •❃ u ✄ n 2 ☎ 2 n n! k 1 ✄ u☎ ; n ✥ ∞✟ u ✮ 1 (76)
where k 1 ✄ u☎ is an undetermined function. Since for u ✥ ∞ and all n ✫ 2,

L n ✄ u☎ ✤❃ u ✄ n 2 ☎ 2 n ✎ ✛ 4n!✢ , we have k 1 ✄ u☎ ✩✥ 1✎ 4 as u ✥ ∞.
We examine asymptotic matching between the expansions for u ✮ 1 and u ✑ 1 ✏ O✄ n ✕ 1 ☎ . If ( 76) matches to e A n ✽ u✕ 1✾ G 0 ✄ w☎ then ( 71 on the page before) would agree with the expansion of ( 76) as u ❺ 1. However this cannot be true as the dominant exponential term in ( 71 on the page before) is

O ❡ exp ✸ 2e w w ✻ ❢ ✏ O ♥ exp ❻ 2e n✽ u✕ 1✾ n✄ u ✑ 1☎ ✬❼ ♦ ( 77
)
while the dominant term in ( 76) as u ✥ 1 is

O ❡ exp ✸ n 2 2 ✄ u ✑ 1☎ ✻ ❢ ✠ (78)
This suggests that yet another expansion is needed on some scale where n ✥ ∞ and u ❺ 1 with n✄ u ✑ 1☎ ❑✥ ▼✒ ∞.

By comparing ( 77 on the page before) to ( 78 on the preceding page), this new scale is likely to be

w ✏ n✄ u ✑ 1☎ ✏ log n ✒ 2 loglog n ✒ O✄ 1☎ (✠ (79)
With this scaling both ( 77 on the page before) and ( 78 on the preceding page) are exp ✛ O✄ n logn☎ ✜✢ . We have not examined this intermediate scale in any detail. However, we would guess that with ( 79), the approximation of ( 6 on page 45) will involve retaining an infinite number of terms in this sum (rather than just the 2 terms in ( 75 on the page before)), but not approximating the sum by an integral, as was possible

when u ✑ 1 ✏ O✄ n ✕ 1 ☎ .
To summarize this section, we have analyzed ( 6 on page 45) in various asymptotic limits. These include

(i) u ✥ 1, n fixed; (ii) u ✥ 1, n ✥ ∞, n✄ u ✑ 1☎ ✏ w fixed; (iii) w ✥ ❆❅ ∞; (iv) 0 ◗ u ◗ 1, n ✥ ∞; and (v) u ✮ 1, n ✥ ∞.
In the next section, we use these results to obtain information about the distribution ✘ ✚✙ ✵✛ ✝ n ✏ k✢ .

Tails of the Limiting Distribution

Using the approximation ( 35 on page 49) for u ✑ 1 ✏ O✄ n ✕ 1 ☎ , we obtain

✘ ✚✙ ✜✛ ✝ n ✑ E✄ ✌✝ n ☎ ✏ ny✢ ✏ 1 2πi ✶ C z ✕ ❑❽ ny✳ A n ✳ 1❾ L n ✄ z☎ dz✠ (80) ❃ 1 n 1 2πi ✶ C ◆ 1 ✒ w n ❖ ✕ ny✕ 1 e ✕ A n log✽ 1✳ w❝ n✾ ⑥✳ A n w❝ n G 0 ✄ w☎ dw ❃ 1 n 1 2πi ✶ Br e ✕ yw G 0 ✄ w☎ dw ✏ 1 n P✄ y☎
where Br ✏ ✄ c ✑ i∞✟ c ✒ i∞☎ for some constant c, is any vertical contour in the w-plane. Here we have set z ✏ 1 ✒ w✎ n in the integral. It follows that

G 0 ✄ w☎ ✏ ✶ ∞ ✕ ∞
e wy P✄ y☎ dy (81) so that G 0 ✄ w☎ is the moment generating function of the density P✄ y☎ . In view of ( 43 on page 50) and ( 44 on page 50) we have ❿ ∞ ✕ ∞ P✄ y☎ dy ✏ 1 and ❿ ∞ ✕ ∞ yP✄ y☎ dy ✏ 0. Observe that, using ( 35 on page 49), ( 41 on page 50), and ( 42 on page 50), we can refine the approximation ( 80) to

✘ ✚✙ ✧✛ ✝ n ✑ E✄ ❈✝ n ☎ ✏ ny✢ ✏ 1 n ✸ P✄ y☎ ✹✒ log n n ✛ P 1 ✄ y☎ ✹✒ P✿ ✿ ✦✄ y☎ ✢ ✒ 1 n ✛ P 2 ✄ y☎ ✹✒ P✿ ✄ y☎ ✺✒ 1 2 yP✿ ✿ ✄ y☎ ✹✒ ❀✄ γ ✑ 2☎ P✿ ✿ ✄ y☎ ✢ ➀✒ o✄ n ✕ 1 ☎ ✻ where P k ✄ y☎ ✏ 1 2πi ✶ Br e ✕ yw G k ✄ w☎ dw (82) for k ✏ 1✟ 2.
An integral equation for P✄ y☎ can easily be derived from ( 49 on page 52). We multiply ( 49 on page 52) by e ✕ wy ✎ ✯✄ 2πi☎ and integrate over a contour Br in the w-plane:

P✄ y ✒ 1☎ ✏ 1 2πi ✶ Br e ✕ w✽ y✳ 1✾ G 0 ✄ w☎ dw (83) ✏ 1 2πi ✶ Br ✶ 1 0 e 2wφ✽ x✾ G 0 ✄ wx☎ G 0 ✄ w ✑ wx☎ e ✕ wy dxdw ✏ ✶ 1 0 1 2πi ✶ Br e ✕ wy ✶ ∞ ✕ ∞ P✄ ξ☎ e wxξ dξ ✶ ∞ ✕ ∞
P✄ η☎ e ✽ w✕ wx✾ η dη e 2wφ✽ x✾ dwdx

✏ ✶ 1 0 ✶ ∞ ✕ ∞ ✶ ∞ ✕ ∞ P✄ ξ☎ P✄ η☎ δ ✄ y ✑ xξ ✑ ❳✄ 1 ✑ x☎ η ✑ 2φ ✄ x☎ ☞☎ dηdξdx ✏ ✶ 1 0 1 x ✶ ∞ ✕ ∞ P✄ η☎ P ✸ y x ✑ 1 ✑ x x η ✑ 2φ ✄ x☎ x ✻ dηdx ✏ ✶ 1 0 ✶ ∞ ✕ ∞ P ✸ xt ✒ y✎ 2 ✑ φ✄ x☎ 1 ✑ x ✻ P ✸ ✑ ❇✄ 1 ✑ x☎ t ✒ y✎ 2 ✑ φ✄ x☎ x ✻ dtdx✟
where we used the well known identity (cf. [START_REF]Handbook of Mathematical Functions[END_REF])

1 2πi ✶ Br e ✕ yw dw ✏ δ✄ y☎ (84)
where δ✄ y☎ is Dirac's delta function. The last expression is precisely ( 15 on page 46). The solution to this integral equation is not unique: if P✄ y☎ is a solution, so is P✄ y ✒ c☎ for any c.

We study P✄ y☎ ✏ ✄ 2πi☎ ✕ 1 ❿ Br e ✕ yw G 0 ✄ w☎ dw as y ✥ ➁❅ ∞. We argue that the asymptotic expansion of the integral will be determined by a saddle point w ✏ s✄ y☎ , which satisfies s✄ y☎ ➂✥ ➁❅ ∞ as y ✥ ②❅ ∞. Thus for y ✥ ❍✑ ∞, we can use the approximation ( 63 on page 55) for G 0 ✄ w☎ , which yields

P✄ y☎ ✩❃ 1 2πi ✶ Br 2 ❊ 2 ❊ π log2 ❊ ✑ w exp ❡ ✄ β ✑ y☎ w ✒ ❛✸ 1 log 2 ✑ 2✻ w log✄ ✵✑ w☎ ❢ dw✠ (85)
This integrand has a saddle point where

d dw ❡ ✑ ❇✄ y ✑ β☎ w ✒ ✸ 1 log 2 ✑ 2✻ w log✄ ✵✑ w☎ ❢ ✏ 0 so that w ✏ ✑ 1 e exp ❡ ✑ y ✑ β 2 ✑ 1✎ log2 ❢ ❲
w✄ y☎ which satisfies w✄ y☎ ✩✥ ⑤✑ ∞ as y ✥ ➃✑ ∞. Then the standard saddle point approximation to ( 85) yields

P✄ y☎ ✚❃ 2 ❊ 2 ❊ π log2 ❊ ✑ w exp ❡ ✄ β ✑ y☎ w ✒ ✸ 1 log 2 ✑ 2✻ w log✄ ✵✑ w☎ ✜❢ ❭ 1 2πi ✶ Br exp ❡ 1 2 w ✸ 1 log 2 ✑ 2✻ |✄ w ✑ w☎ 2 ❢ dw (86) ✏ 2 πe 1 ❊ 2 log2 ✑ 1 exp ❡ β ✑ y 2 ✑ 1✎ log2 ✑ 2 ✑ 1✎ log2 e exp ✸ β ✑ y 2 ✑ 1✎ log2 ✻ ❢
for y ✥ ➄✑ ∞. Thus, the left tail is very small and the behavior of P✄ y☎ as y ✥ ➄✑ ∞ is similar to that of an extreme value distribution (i.e., double-exponential distribution).

w■ ✏ log ◆ y 2 ❖ ✒ loglog ◆ y 2 ❖ ✒ log log✄ y✎ 2☎ log✄ y✎ 2☎ ✄ 1 ✒ o✄ 1☎ ✧☎ ✹✟ y ✥ ∞ ✠ (89)
For fixed z and y we have, as n ✥ ∞,

✘ ✚✙ ✵✛ ✝ n ✑ E✄ ❈✝ n ☎ ✚❉ nz✢ ➆❃ ✶ z ✕ ∞ P✄ y☎ dy (90) ✘ ✚✙ ✧✛ ✝ n ✑ E✄ ❈✝ n ☎ ✫ ny✢ ➇❃ ✶ ∞ y P✄ u☎ du✠ (91)
If z ✥ ➈✑ ∞ or y ✥ ➁✒ ∞, then these integrals may be evaluated asymptotically using ( 86 on the preceding page) and ( 88), and we obtain the results ( 17 on page 46) and ( 18 on page 47), respectively.

This derivation of the expansions of P✄ y☎ as y ✥ ⑨❅ ∞ has the disadvantage in that it assumes the existence of certain saddle points. However, we can obtain the identical results simply by using the integral equation for P✄ y☎ , which we now show.

Let us write ( 83 on the preceding page) in the form After some calculation, we obtain the following linear differential equation for G: where we have used ( 100 on the preceding page) to change the variable of integration. But then the dominant exponential terms in ( 88 on page 60) and ( 98 on the preceding page) agree precisely. To solve ( 101 on the page before) we change variables y ➌ w■ and write ( 101 on the preceding page) as

P✄ y ✒ 1☎ ✏ ✶ 1 0 1 x✄ 1 ✑ x☎ ✶ ∞ ✕ ∞ P ✸ ξ ✒ y✎ 2 ✑ φ ✄ x☎ 1 ✑ x ✻ P ✸ ✑ ξ ✒ y✎
G✿ ✒ 1 2 F ✿ ✿ G ✏ G yF ✿ ✑ 2G y ✛ γ ✒ log✄ ✼✑ yF ✿ ☎ ✢ ❞✑ G✿ F ✿ ✒ F ✿ ✿ ✄ F ✿ ☎
G✿ ✄ y☎ G✄ y☎ ✏ ✑ F ✿ ✿ ✄ y☎ G dG dw■ ✏ ✸ 1 ✒ 1 F ✿ ✻ ✕ 1 ❡ F ✿ ✿ ✄ F ✿ ☎ 2 ✑ 1 2 F ✿ ✿ ✒ 1 yF ✿ ✒ 2F✿ y ✑ 2✄ γ ✒ log2☎ y ❢ (103)
where we have used log✄ ✵✑ yF ✿ ☎ ✏ log✄ 2e

  •✥ ❍✑ ∞ sufficiently slowly, where α ✏

								2log 2✕ 1 e log 2 ✏ 0✠ 205021 ✠ ✡✠ ☞✠ and β is a constant.
	(iii) The right tail becomes
	✘ ✚✙ ✜✛ ✝ n ✑ E✄ ❈✝ n ☎ ✫ ny✢ ❄❃	C■ 2 ❊ π	w■	❊	1 w■ ✑ 1	e ✕ w❏ ✧✽ w❏ ✧✕ 1 2 ✳ 2γ✳ log2✾ exp ✸ ✑ yw■ ❑✒	✶ w❏ 1	2e u u	du✻	(18)
	for n ✥ ∞, and y ✏ y✄ n☎ ▲✥ ▼✒ ∞ sufficiently slowly. Here C■ is a constant, γ is Euler's constant and w■ ✏ w■ ✄ y☎
	is the solution to				y ✏	2 w■	exp✄ w■ ☛☎

  ✏ εE✄ ✌✝ n ☎ ✧✎ n ✏ 2ε logn ✒ |✒ ε✄ 2γ ✑ 4☎ ✹✒ O✄ ε✄ log n☎ ✧✎ n☎ and use[START_REF] Rösler | A Fixed Point Theorem for Distributions[END_REF] to approximate w■ as y ✥ ∞. This yields

	y☎ ✏	✑ yw■ ✯✒	✶ w❏ 1	2e u u	du✟	y ✏	2e w❏ w■	✟	(24)
	k ✄ y☎ ✏	2 ✕ 2γ✕ log2☎ e ✕ w 2 ❏ e w❏ ✄ 1 w■ ❊ w■ ❑✑ 1	✠
	We have not been able to determine the upper limit on y for the validity of ( 23). However, it is easy to see
	that ( 23) reduces to ( 21) if we set y

  e Ψ✽ w✕ wx✾ dx

		wΨ✿ ✄ w☎ •❃ 2e w
	and hence we define	Ψ✄ w☎ ✏ 2	✶ w e u u	du✠	(66)

❃ 2K ✄ w☎ e Ψ✽ w✾ ✶ 1❝ 2 0 e ✕ wΨ⑦ ✐✽ w✾ x dx✠ (65)

If furthermore Ψ✄ ☞❁ ❂☎ is such that wΨ✿ ✄ w☎ •✥ ∞ as w ✥ ∞, then ( 65) yields (after a slight rearrangement)

  Now take y ✥ ❆✒ ∞ and use ( 71 on page 56) to get ✏ 2e w ✎ w. Let w■ ✏ w■ ➅✄ y☎ be the solution to ( 19 on page 47) that satisfies w■ ➂✥ ∞ as y ✥ ∞. Then expanding the integrand in ( 87) about w ✏ w■ ✷✄ y☎ and using the standard saddle point approximation yields ✥ ∞, from which ( 18 on page 47) easily follows. Thus for y ✥ ∞ we have P✄ y☎ ✏ exp ✛ O✄ ✵✑ y logy☎ ✜✢ and hence the right tail is thinner than the right tail of an extreme value distribution. ¿From ( 19 on page 47) it is easy to show that

		P✄ y☎ ✩❃	1 2πi	✶ Br C■ exp ✸ ✶ w 1	2e u u	du✻	1 w	e ✕ yw✕ w 2 ✳ ❑✽ 1✕ 2γ✕ 2log 2✾ w dw✠	(87)
	The saddle point now satisfies
						d dw	❡ ✑ yw ✒ ✶ w 1	2e u u	du❢ ✏ 0
	or y P✄ y☎ ✚❃	2	C■ ❊	2π	❊ y ⑩ 1 ✑ 1✎ w■	exp ❡ ✑ yw■ ✒ ✶ w❏ 1	2e u u	du ✑ w 2 ■ ✑ ❳✄ 2γ ✒ 2 log2☎ w■ ✼❢	(88)
	as y

  2 ✑ φ✄ x☎ x ✻ dξdx (92) and assume that for y ✥ ➃✑ ∞, P✄ y☎ has the form P✄ y☎ ✩❃ de cy exp ✛ ✑ ae ✕ by ✢ (93) for some constants a✟ b✟ c✟ d. Using ( 93) in ( 92) we find that the major contribution to the double integral will come from ✄ x✟ ξ☎ ✏ ✰ 1 2 ✟ 0 ✲ . After scaling x ✏ 1 2 ✄ 1 ✒ e by❝ 2 u☎ and ξ ✏ e by❝ 2 η we obtainde c e cy exp ✛ ✑ ae ✕ b e ✕ by ✢ ✺❃ 2e by d 2 exp ✛ ✑ 2ae ✕ by e ✕ 2blog 2 ✢ ✒ 2 log2☎ ✺✒ 2η☎ 2 b 2 ✑ 2νb✄ u 2 ✄ y ✒ 2 log 2 ✑ 1☎ ✹✒ 2ηu☎ ✵➊ dudη where ν ✏ ae ✕ 2blog 2 . The double integral is equal to π✄ 2ν☎ ✕ 1 b ✕ 3❝ 2 ✄ 2 ✑ 1✎ b☎ ✕ 1❝2 and thus the above yields 2e ✕ 2blog 2 ✏ e ✕ b ✟ e cy ✏ e ✽ b✳ 2c✾ y (94) de c ✏ 2d 2 π 2a e 4c log 2 b ✕ 3❝ 2 e 2blog 2 ✸ 2 ✑ 1 b remains arbitrary. But then ( 93 on the page before) is equivalent to ( 86 on page 59) if we identify a as Next we consider y ✥ ⑤✒ ∞ and write the equation for P✄ y☎ in the formP✄ y ✒ 1☎ ✏ 2 P✄ η☎ P ✸ y ✒ x✄ y ✑ η☎ ✯✑ 2φ✄ x☎ 1 ✑ x ✻ dηdx✠ (97)We seek solutions in the WKB formP✄ y☎ ✚❃ G✄ y☎ e F ✽ y✾ ✟ y ✥ ⑨✒ ∞(98)where log G varies less rapidly than F. Then ( 97) may be approximated bye F ✽ y✾ ⑥✳ F⑦ ✽ y✾ ✛ G✄ y☎ ✺✒ G✿ ✄ y☎ ✺✒ ❹❁ ✡❁ ☛❁ ✢ ❡ 1 ✒ P✄ η☎ e F ✽ y✾ exp ❣ ❡ x✄ y ✑ η☎ 1 ✑ xAfter scaling x ✏ ω✎ ✛ ✑ yF ✿ ✄ y☎ ✢ , the leading term in the right side of (99) is e F ✽ y✾ G✄ y☎ (✄ ✵✑ 2☎ ✧✎ ✯✄ yF ✿ ✄ y☎ ✧☎ and

												✕ 1❝ 2
												✻	✠
	We thus have				b ✏	1 2 ✑ 1✎ log2	✟ c ✏ ✑ b✟	d ✏	2a π	b 3❝ 2	④	2 ✑	1 b	(95)
	and a a ✏	2 ✑ 1✎ log2 e	exp ❡ 2 ✑ 1✎ log2 β	❢ ✠	(96)
						✶ 1❝ 2	1		✶ ∞
						0	1 ✑ x	✕ ∞
												1 2	F ✿ ✿ ✄ y☎ ✹✒ ❀❁ ☛❁ ☛❁ ❢ ✏	(99)
	✏	2	✶ 1❝ 2 0	1 1 ✑ x	✶ ∞ ✕ ∞							✑	2φ✄ x☎ 1 ✑ x	❢ F ✿ ✄ y☎ ❤
	❭	❡ G✄ y☎ ✹✒	x✄ y ✑ η☎ ✯✑ 2φ✄ x☎ 1 ✑ x	G✿ ✄ y☎ ✹✒ ❀❁ ☛❁ ☛❁ ❢ ♥ 1 ✒	1 2
	hence								e F ⑦ ❂✽ y✾ ✏ ✑	2 yF ✿ ✄ y☎

❭

e 2cy e 4c log 2

✶ ∞ ✕ ∞ ✶ ∞ ✕ ∞ exp ➉ ✌✑ ν✄ u✄ y ✸ x✄ y ✑ η☎ ✯✑ 2φ✄ x☎ 1 ✑ x ✻ 2 F✿ ✿ ✄ y☎ ✹✒ ❀❁ ☛❁ ✡❁ ♦ dηdx✠ ✠

(100)

  2 G✠ (101) Equation (100) is a first order non-linear ordinary differential equation for F, which is really just a transcendental equation for F ✿ . Setting F ✿ ✄ y☎ ✩❲ ➋✑ w■ ✬✄ y☎ we have

	F ✄ y☎ ✏ ✑	✶ y	w■ ✬✄ t ☎ dt ✏	✑	✶ w❏ ✧✽ y✾	2ζ ❻ e ζ ζ	✑	e ζ ζ 2 ❼	dζ	(102)
			✏ ✑ yw■ ✯✒ 2	✶ w❏ e ζ ζ	dζ

  ✕ F⑦ ☎ ✏ log 2 ✑ F✿ . From ( 100 on the page before) we get 2w■ ❑✑ 2✄ γ ✒ log2☎ (✠ (105) Solving ( 105) we see that ( 98 on the page before) is equivalent to ( 88 on page 60), noting that ❊ y ✏ ⑩ 2✎ w■ exp ✄ w■ ✎ 2☎ . Of course, G is determined only up to a multiplicative constant, which corresponds to C■ ☛✎ ✯✄ 2 ❊ π☎ in ( 88 on page 60).

			e w❏ ✸	1 w■	✑	1 w 2 ■	✻ ✄ ✼✑ F ✿ ✿ ✄ y☎ ✧☎ ✏	1 2	✠	(104)
	Using ( 104) in ( 103) we obtain, after some calculation,
	1 G	dG dw■	✏ 1 2	✑	1 2	1 w■ ✑ 1

✑

† The work of this author was supported by NSF Grant DMS-93-00136 and DOE Grant DE-FG02-93ER25168. ‡ The work of this author was supported by NSF Grants NCR-9415491 and CCR-9804760, and NATO Collaborative Grant CRG.950060. 1365-8050 c ✆ 1999