text:
# Homoclinic crossing in open systems: Chaos in periodically perturbed monopole plus quadrupolelike potentials

## I Introduction

Since the pioneering work of Poincaré poincare1  in celestial mechanics in which the mathematical basis of deterministic chaos in compact phase space systems was laid down, the study of homoclinic phenomena in closed systems with hyperbolic unstable periodic orbits has allowed the understanding of a rich variety of nonlinear effects in physics, chemistry, and biology bai . Due to its universality, models in which unstable periodic orbits are subjected to small periodic perturbations has become one of the main paradigms of deterministic chaos arn . An analytical tool to study such models is the Melnikov method poincare2 ; melnikov ; arnold  in connection with Smale-Birkhoff homoclinic theorem birkhoff ; smale , and Kolmogorov-Arnold-Moser (KAM) theory in the Hamiltonian case moser .

The Melnikov function describes the transversal distance between the stable and unstable manifolds associated to an unstable periodic orbit. Its isolated odd zeros indicate transversal intersections between these manifolds, and hence the onset of chaos guckenheimer . Examples of applications of the Melnikov method in gravitation are the motion of particles in perturbed two- and three-dimensional Stäckel potentials gerh1 ; gerh2 , the chaotic evolution of cosmological models koiller , the study of orbits around a black hole perturbed by either gravitational radiation vieira  or an external quadrupolar shell moeckel , and the bounded motion of particles in a periodically perturbed attractive center described by a monopole plus a quadrupolelike potential were considered in Ref. letelier .

The Melnikov method has also been used in many other branches of physics. We find examples of applications of this method to the study of Josephson junctions yao ; ken , planar periodic vortical flows tae , solitons abd , liquid crystals ste , and transfer dynamics of quasiparticles hen .

Even for Hamiltonian systems fundamental questions about chaos in non-compact phase space systems remain to be answered. Among the more important unsolved questions are the notion of chaos itself and the lack of an adequate theory to deal with it. Partial results obtained in this area are the fractal techniques in scattering processes ott , these are numerical techniques that present some difficulties due to the existence of different time scales for nearly bounded scattering. They are inadequate to the study of chaotic behavior arising from separatrices between bounded and unbounded orbits. Furthermore, they are unable to present a complete description of the chaotic motion as the one provided by analytical methods in the closed system case.

The aim of this paper is to study the homoclinic phenomenon for a class of open systems that by a suitable change of coordinates can be approached in terms of an adequate formulation of the Melnikov method and Smale-Birkhoff homoclinic theorem. The change of coordinates regularizes the unstable periodic orbit at infinity and it compactifies the region of interest of the phase space; however, the phase space as a whole remains noncompact. Alas the resulting unstable periodic orbit is typically nonhyperbolic and the standard stable manifold theorem, needed to state the Melnikov method guckenheimer , does not apply. McGehee gehee  extended this theorem to degenerated cases in the context of the Newtonian three body problem. Xia xia , and Dankowicz and Holmes dankowicz , among others, used McGehee’s result in connection with Melnikov method and Smale-Birkhoff homoclinic theorem to study the nonintegrability of the three body problem.

Here we consider the equatorial motion of a particle moving in a potential described by a monopolar term plus a quadrupolelike contribution. This potential models the gravitational attraction of a galaxy bulge or any nonspherical celestial body; it also arises in general relativity in the study of the motion of a test particle around a Schwarzschild black hole, the quadrupole term being a general relativistic effect associated to the angular momentum of the particle in the reduced two-dimensional phase space, see for instance Ref. chandra .

In Sec. II the fixed saddle points associated to the monopole plus quadrupole system, as well as the coordinate transformation that regularizes these points at infinity, are studied. In the next section we present some mathematical preliminaries and the Melnikov method. The equations of motion are used to reduce this method to the analysis of simple graphics. We find that the perturbation induces transverse homoclinic orbits in some ranges of the parameters, and we apply the modified Smale-Birkhoff homoclinic theorem to verify the presence of a symbolic dynamics equivalent to a Smale horseshoe map; see Ref. guckenheimer . In Sec. IV the study of the motion is completed with a presentation of Poincaré sections that reveal different levels of chaotic behavior as a function of the parameters. Finally, in the last section, we make some remarks about the class of system in which the same kind of analysis can be performed.

## II The homoclinic orbit and perturbations

We shall consider the orbit of a particle in a plane under the influence of a force modeled by a potential with inverse square law plus a quadrupolelike term. It is convenient to work with dimensionless quantities. The motion of the particle is described by letelier

$$H_0=\frac{p^2}{2}+\frac{1}{2r^2}\frac{1}{r}\frac{\beta }{r^3}\left(p=\frac{dr}{dt}\right),$$
(1)
where $`r`$, $`p`$, $`t`$, $`H_0`$, and $`\beta `$ are dimensionless quantities proportional to, respectively, the radius, the radial momentum, the time, the Hamiltonian function, and the quadrupole moment of the attraction center. The effective potential

$$V_{eff}=\frac{1}{2r^2}\frac{1}{r}\frac{\beta }{r^3},$$
(2)
is presented in Fig. 1 for different values of the parameter $`\beta 0`$.

The natural space to study a periodically perturbed planar system is $`IR^2\times S^1`$, where the unstable periodic orbits have a proper meaning guckenheimer . In the corresponding unperturbed autonomous case the phase space is in $`IR^2`$ and the unstable periodic orbits reduce to fixed saddle points. The unstable periodic orbits are also reduced to fixed saddle points for the maps defined on Poincaré sections $`IR^2\times \{t_0\}IR^2\times S^1`$.

The above system presents a homoclinic loop associated to the hyperbolic fixed saddle point at $`(r,p)=(r_M,0)`$, where $`1/r_M=1/6\beta +\sqrt{\left(1/6\beta \right)^21/3\beta }`$, for $`\beta `$ limited by $`1/16<\beta <1/12`$. This case is important in the study of bounded orbits and was explored in Ref. letelier .

In the present work we study the instabilities of unbounded orbits, the relevant values of the parameter are $`\beta =1/16`$ and $`0\beta <1/16`$. Let $`H_0=0`$, if $`\beta =1/16`$ the points $`r=r_M=1/4`$ and $`r=\mathrm{}`$ on the $`r`$ axis represent fixed points such that the particle takes infinity time to reach to or to depart from each of these points. For zero energy still, if $`0\beta <1/16`$ the motion of the particle is restricted to the region between $`r=r_{}`$ and $`r=\mathrm{}`$, where $`1/r_{}=1/4\beta \sqrt{(1/4\beta )^21/\beta }`$, and only the last point represents a fixed point in this range of values of $`\beta `$.

Since the orbits of interest are in a semi-infinity region bounded away from the origin, we can compactify this part of the phase space with a change of the position coordinate like $`r=u^\alpha `$ with $`\alpha <0`$. We find that the transformation $`r=1/u^2`$ allows us to model the problem in a way similar to Refs. gehee ; xia ; dankowicz . This new coordinate regularizes the fixed point at infinity that now is at the point $`(u,p)=(0,0)`$. The zero energy orbits generate a heteroclinic loop for $`\beta =1/16`$ associated to the hyperbolic fixed saddle point at $`(u,p)=(2,0)`$ and the degenerated fixed saddle point at $`(u,p)=(0,0)`$, and a homoclinic loop for $`0\beta <1/16`$ associated to the degenerated fixed saddle point at $`(u,p)=(0,0)`$. Degenerated in the sense that both eigenvalue of the linearized vector field are zero, as can be seen from the Hamiltonian equations

$`{\displaystyle \frac{du}{dt}}={\displaystyle \frac{1}{2}}u^3p,`$ (3)
$`{\displaystyle \frac{dp}{dt}}=u^4+u^63\beta u^8.`$ (4)
The points $`(u,p)=(2,0)`$ and $`(u,p)=(0,0)`$ correspond to, respectively, the hyperbolic fixed saddle point at $`(r,p)=(1/4,0)`$ for $`\beta =1/16`$ and the degenerated fixed saddle point at $`(r,p)=(\mathrm{},0)`$ for $`0\beta 1/16`$. The homoclinic and heteroclinic loops are defined by the intersection between the stable and unstable manifolds on the $`(u,p)`$ plane. The homoclinic loop for $`\beta =1/18`$ is shown in Fig. 2 and the heteroclinic one is presented in Fig. 3.

The explicit integration of the homoclinic and heteroclinic loops will be necessary to apply the Melnikov method and can be obtained from the first integral of motion ($`H_0=0`$). For $`\beta =0`$ we find

$$t(v)=\pm \frac{(1+v)(2v)^{1/2}}{3v^{2/3}},$$
(5)
where $`v=1/r`$, the time origin is take in the symmetry point of loop, and the sign refers to the upper ($`+`$) and lower ($``$) parts of the loop. Analogously, for $`0<\beta <1/16`$ it reads

$`t(v)=\pm \sqrt{{\displaystyle \frac{2}{\beta }}}\{{\displaystyle \frac{(2v_++v_{})F(\delta ,q)2(v_++v_{})E(\delta ,q)}{3v_{}^2v_+^{3/2}}}`$
$`+{\displaystyle \frac{v_+v_{}+(2v_++v_{})v}{3v_+v_{}^2}}\sqrt{{\displaystyle \frac{v_{}v}{(v_+v)v^3}}}\}`$ (6)
with

$$\delta =\mathrm{arcsin}\sqrt{\frac{v_+(v_{}v)}{v_{}(v_+v)}},q=\sqrt{\frac{v_{}}{v_+}},$$
(7)
where $`F(\delta ,q)`$ and $`E(\delta ,q)`$ are elliptic integrals of first and second type in the Legendre normal form (Ref. gradshteyn , p. 224), and $`v_\pm =1/4\beta \pm \sqrt{(1/4\beta )^21/\beta }`$ are the roots of $`V_{eff}=0`$. Moreover, for $`\beta =1/16`$ we get

$`t(v)=\pm \{{\displaystyle \frac{\sqrt{2}}{3}}({\displaystyle \frac{1}{v^{3/2}}}{\displaystyle \frac{1}{(4/3)^{3/2}}})+{\displaystyle \frac{1}{2\sqrt{2}}}({\displaystyle \frac{1}{\sqrt{v}}}{\displaystyle \frac{1}{\sqrt{4/3}}})`$
$`+{\displaystyle \frac{1}{8\sqrt{2}}}\mathrm{ln}\left[\left({\displaystyle \frac{2\sqrt{v}}{2+\sqrt{v}}}\right)\left({\displaystyle \frac{2+\sqrt{4/3}}{2\sqrt{4/3}}}\right)\right]\}`$ (8)
with the choice $`t(4/3)=0`$, where $`v=4/3`$ is the local minimum of $`V_{eff}`$.

Now, let us consider the Hamiltonian (1) perturbed by a periodic multipolar term of the form

$$H=H_0+\epsilon H_1,$$
(9)
$$H_1=r^n\mathrm{cos}(\mathrm{\Omega }t)(n2),$$
(10)
where $`n=2`$ is dipolar, $`n=3`$ is quadrupolar, etc. These perturbations can model the attraction due to a distribution of masses with periodic motions that are placed inside the planet orbit.

We shall consider our attraction center with a fix total mass. In other words, we excluded the monopolar case ($`n=1`$) that represents a periodic variation of the mass.

In the next sections we study how these perturbations can affect the dynamics of the system.

## III Melnikov method

Powerful tools to study near integrable systems are Melnikov type of techniques that detect transversal intersections between the stable and unstable manifolds associated to a unstable periodic orbit. The presence of such transversal intersections is a guarantee of complicated dynamics and in some cases leads to a symbolic dynamics equivalent to the Smale horseshoe guckenheimer .

In order to simplify the analysis it is convenient to abstract a little from the particular problem presented above. In what follows the loops are on the $`X`$ plane and the manifolds are in the $`(X,\theta )`$ space, where $`\theta t`$ mod $`2\pi /\mathrm{\Omega }`$ so that $`(X,\theta )IR^2\times S^1`$, and $`\mathrm{\Sigma }_{\theta _0}`$ denotes the section $`\theta =\theta _0`$. We consider a Hamiltonian of the form

$$\stackrel{~}{H}(X,t)=\stackrel{~}{H}_0(X)+\epsilon \stackrel{~}{H}_1(X,t)(\stackrel{~}{H}_1\mathrm{\hspace{0.33em}2}\pi /\mathrm{\Omega }\text{periodic in }t),$$
(11)
where $`\stackrel{~}{H}_0`$ is integrable with homoclinic (heteroclinic) loop $`\mathrm{\Gamma }`$ associated to some hyperbolic fixed saddle point(s). Under hyperbolicity hypothesis it can be shown that for sufficiently small $`\epsilon `$ the invariant manifolds are only deformed, and possibly their intersections become transversal, see for instance Ref. guckenheimer .

Let $`X_0`$ be a point on $`\mathrm{\Gamma }`$ and $`X^{s/u}(\theta _0,\epsilon )`$ be points on the stable/unstable manifolds such that they are on $`\mathrm{\Sigma }_{\theta _0}`$, in the line perpendicular to $`\mathrm{\Gamma }\times \{\theta _0\}`$ at $`X_0`$ and whose trajectories take the least amount of time to reach/depart any small neighborhood of the unstable periodic orbit. A computable measure of the transversal distance between the stable and unstable manifolds on $`\mathrm{\Sigma }_{\theta _0}`$, which defines the Melnikov function, is given by the zero order term of $`1/\epsilon [\stackrel{~}{H}_0(X^u(\theta _0,\epsilon ))\stackrel{~}{H}_0(X^s(\theta _0,\epsilon ))]`$ ozorio . In fact, if $`X^{s/u}(t;\theta _0,\epsilon )`$ denotes the time evolution under $`\stackrel{~}{H}`$ such that $`X^{s/u}(\theta _0;\theta _0,\epsilon )=X^{s/u}(\theta _0,\epsilon )`$, and $`X_0(t)`$ denotes the time evolution under $`\stackrel{~}{H}_0`$ such that $`X_0(0)=X_0`$,

$`\stackrel{~}{H}_0(X^{s/u}(\theta _0,\epsilon ))\stackrel{~}{H}_0(X^{s/u}(\pm \mathrm{};\theta _0,\epsilon ))={\displaystyle _\pm \mathrm{}^{\theta _0}}{\displaystyle \frac{d\stackrel{~}{H}_0}{dt}}\left[X^{s/u}(t;\theta _0,\epsilon )\right]𝑑t`$
$`={\displaystyle _\pm \mathrm{}^{\theta _0}}{\displaystyle \frac{d\stackrel{~}{H}_0}{dt}}[X_0(t\theta _0),t]𝑑t+O(\epsilon ^2).`$ (12)
Thus in the homoclinic case the Melnikov function can be written as

$$M(\theta _0)=\frac{1}{\epsilon }_{\mathrm{}}^{\mathrm{}}\frac{d\stackrel{~}{H}_0}{dt}[X_0(t),t+\theta _0]𝑑t,\left(\frac{1}{\epsilon }\frac{d\stackrel{~}{H}_0}{dt}=\{\stackrel{~}{H}_0,\stackrel{~}{H}_1\}\right),$$
(13)
where $`\{.,.\}`$ are the usual Poisson brackets.

The implicit function theorem allows us to conclude that if $`M(\theta _0)`$ has simple zeros, then, for sufficiently small $`\epsilon `$ the invariant manifolds intersect transversely for some $`\theta _0`$. On the other hand, if $`M(\theta _0)`$ is bounded away from zero, then the invariant manifolds do not intersect for all $`\theta _0`$.

Now let us take the map defined by system (3)-(4) on a arbitrary section. After a scale change in the $`p`$ coordinate it reads

$`u_{k+1}=u_kCu_k^3\left[p_k+O(4)\right],`$ (14)
$`p_{k+1}=p_kCu_k^3\left[u_k+O(3)\right],`$ (15)
where $`C=\sqrt{2}\pi /\mathrm{\Omega }`$. Here the standard Melnikov method breaks down because of the degeneracy of the saddle point. McGehee gehee , Xia xia , Dankowicz and Holmes dankowicz  studied systems of this class in the context of the three body problem, where they established the fundamental results needed to support the Melnikov method: structural stability of the unstable periodic orbits (trivial in our case since the degenerated unstable periodic orbit remains fixed); existence of local stable and unstable analytic manifolds $`C^{\mathrm{}}`$ close to those of the unperturbed case gehee ; dankowicz ; solutions on the perturbed and the unperturbed manifolds approach to the unstable periodic orbit at a similar rate xia ; dankowicz . Following the proofs step-by-step we can observe that all the essential hypothesis involved to achieve their results are also satisfied by the above system. Thus these statements apply to Eqs. (14)-(15) allowing the expansion in Eq. (12) that justifies the use of Melnikov method in the present problem.

Another important result whose standard form assumes hyperbolicity is the Smale-Birkhoff homoclinic theorem. It was given a formulation of this theorem that is valid for the degenerated problem of Sitnikov dankowicz , which is grounded on a suitable approximation of the linearization of the map in the neighborhood of the saddle point. Since such approximation results in the same expressions for the case of Eqs. (14)-(15), we conclude that the Smale-Birkhoff homoclinic theorem applies to the above system. Thus transversal homoclinic intersections in our problem lead to the Smale horseshoe.

Now we shall apply Melnikov method to Eq. (1) subjected to the perturbations (10). In the homoclinic case $`(0\beta <1/16)`$ we find for the Melnikov function

$`M(\theta _0)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}nr^{n1}\mathrm{cos}[\mathrm{\Omega }(t+\theta _0)]{\displaystyle \frac{dr}{dt}}𝑑t`$
$`=`$ $`2n\mathrm{sin}(\mathrm{\Omega }\theta _0)K(\mathrm{\Omega }),`$
$`K(\mathrm{\Omega })`$ $``$ $`{\displaystyle _0^v_{}}v^{n1}\mathrm{sin}[\mathrm{\Omega }t(v)]𝑑v,`$ (17)
where we have fixed $`X_0=[v(0),p(0)]`$ as the symmetry point on the loop and the integrand $`t(v)`$ means the positive branch of Eqs. (5)-(6). Thus the Melnikov function has simple zeros as long as $`K(\mathrm{\Omega })0`$. With the change $`tv`$ we pass from an infinite interval in Eq. (III) to a finite one in Eq. (17), and it allows us to study $`K(\mathrm{\Omega })`$ using graphics. Although $`u`$ ($`r=1/u^2`$) is important to justify the Melnikov method, we have great freedom in the choice of a coordinate to study the results. The most simple one that is adequate to this end is the coordinate $`v=1/r`$.

The integrand of $`K(\mathrm{\Omega })`$ is formed by the product of an oscillating function and a polynomial. Near the origin this oscillation is a rapid one, since $`t\mathrm{}`$ as the particle goes to the unstable periodic orbit. In Fig. 4 we show a graph of $`\mathrm{sin}[\mathrm{\Omega }t(v)]`$ for $`\beta =1/25`$, and several values of $`\mathrm{\Omega }`$. For $`\mathrm{\Omega }>4`$ we will have more zeros in the interval shown in the figure. For $`\mathrm{\Omega }<1`$ the curve will look like the one for $`\mathrm{\Omega }=1`$. Since the area under the curves are clearly not null the integral of $`\mathrm{sin}[\mathrm{\Omega }t(v)]`$ is nonzero for $`\beta =1/25`$ and $`0<\mathrm{\Omega }4`$. The cases of interest are $`n2`$, where we will have a more favorable situation. We have transverse homoclinic orbits in all these cases. To better understand this behavior we show in Fig. 5 a graph of the integrand of $`K(\mathrm{\Omega })`$ for $`\beta =1/25`$, $`\mathrm{\Omega }=3`$ and different values of $`n`$. For $`0\beta <1/16`$ the graphics of $`\mathrm{sin}[\mathrm{\Omega }t(v)]`$ will look like the one for $`\beta =1/25`$, with almost the same upper bound for $`\mathrm{\Omega }`$. See for instance Fig. 6, where we plot the graphics of $`\mathrm{sin}[\mathrm{\Omega }t(v)]`$ for the monopolar attraction center ($`\beta =0`$) for the same values of $`\mathrm{\Omega }`$ employed in Fig. 4. Therefore we have transverse homoclinic orbits for all $`\beta `$ limited by $`0\beta <1/16`$.

For these range of parameters the Smale-Birkhoff homoclinic theorem implies the existence of a hyperbolic invariant set for which the action of an $`N`$th iterate of the map has a symbolic dynamic equivalent to that of the Smale horseshoe map. Some important consequences of this result are sensitive dependence on initial condition (a characteristic of chaos); nonexistence of real analytic integral of motion (nonintegrable system); existence of infinitely many periodic orbits with arbitrary large periods (whose number increases exponentially with the period); capture of orbits by the system (in both directions of time).

It is illustrative to see how the chaotic orbits look in the original noncompactified coordinates $`r,p`$. Due to the lack of an extra integral of motion the particle can have a highly erratic motion and have access to a two dimensional region of the phase space. The sensitive dependence on initial condition implies that the evolution of two infinitesimally near points in the space $`r,p`$ can result in two completely different bounded orbits, in two completely different unbounded orbits, or even in one bounded and one unbounded orbit. Orbits that are bounded for all $`t<0`$ can go to infinite for $`t+\mathrm{}`$, and orbits of particles coming from infinite can remain bounded for all $`t>0`$. However, regular orbits are also present and in particular there is a family of periodic orbits.

Let us consider the heteroclinic case $`\beta =1/16`$. The Melnikov method obtained from Eq. (12) applies to each branch of the loop. The distance between $`X^{s/u}(\pm \mathrm{};\theta _0,\epsilon )`$ and $`X^{s/u}(\pm \mathrm{};\theta _0,0)`$ is of order $`O(\epsilon )`$. Therefore $`H_0(X^{s/u}(\pm \mathrm{};\theta _0,\epsilon ))`$ is of order $`O(\epsilon ^2)`$. Then, this term can be neglected and the Melnikov function reads

$`M(\theta _0)=2n\mathrm{sin}(\mathrm{\Omega }\theta _0)K_s(\mathrm{\Omega })2n\mathrm{cos}(\mathrm{\Omega }\theta _0)K_c(\mathrm{\Omega }),`$ (18)
$`K_s(\mathrm{\Omega }){\displaystyle _0^v_{}}v^{n1}\mathrm{sin}[\mathrm{\Omega }t(v)]𝑑v,K_c(\mathrm{\Omega }){\displaystyle _0^v_{}}v^{n1}\mathrm{cos}[\mathrm{\Omega }t(v)]𝑑v.`$ (19)
The orbits in the heteroclinic loop are less symmetrical than in the homoclinic one and, consequently, the Melnikov function has simple zeros when at least one of the integrals of the previous formula is different from zero. The integrands of $`K_s`$ and $`K_c`$ oscillate very rapidly near the unstable periodic orbits. To better understand this behavior we show in Fig. 7 a graphic of the positive branch of Eq. (8). The function $`t(v)`$ has small values for a large range of values of $`v`$, and hence $`K_c`$ will result nonzero values for small $`\mathrm{\Omega }`$. Indeed, we show in Fig. 8 a graph of $`\mathrm{cos}[\mathrm{\Omega }t(v)]`$ for different values of $`\mathrm{\Omega }`$. Since cosine is an even function, the same figure is valid for the negative branch. For $`\mathrm{\Omega }>3`$ we will have more zeros in the interval $`1<v<4`$. For $`\mathrm{\Omega }<0.5`$ the curves will look like the one for $`\mathrm{\Omega }=0.5`$. For small values of $`n`$ the integral $`K_c`$ will be nonzero for $`0<\mathrm{\Omega }3`$. Due to the change of sign of $`\mathrm{cos}[\mathrm{\Omega }t(v)]`$ near $`v=4`$, the upper bound for $`\mathrm{\Omega }`$ decrease with $`n`$. But it is clear from Fig. 8 that for each $`n`$ will exist a upper bound $`\mathrm{\Omega }_o(n)>0`$ such that $`K_c`$ will be nonzero, leading to the presence of transversal heteroclinic orbits, for $`0<\mathrm{\Omega }\mathrm{\Omega }_o(n)`$.

## IV Poincaré sections method

The system (9) has four parameters, $`\beta `$, $`n`$, $`\mathrm{\Omega }`$, and $`\epsilon `$. In opposition to the Melnikov method, Poincaré sections method is able to predict results only for fixed values of these parameters. However, the Poincaré method can locate the regular and irregular regions and provide a qualitative idea of chaotic behavior.

The perturbed phase space is in $`IR^2\times S^1`$ and the maps defined on different values of the angular variable $`\theta =t`$ mod $`2\pi /\mathrm{\Omega }`$ are topologically conjugated. So we have restricted the study to sections built on $`\theta =0`$, but for a large number of different values of the parameters $`\beta `$, $`n`$, $`\mathrm{\Omega }`$, and $`\epsilon `$. We included $`\mathrm{\Omega }`$ values which would require a more elaborate numerical computation of the Melnikov function. In Fig. 9 we show the Poincaré sections for some select values of $`\beta `$, $`n`$, $`\mathrm{\Omega }`$, and $`\epsilon `$. The general aspect of the sections is represented in this figure, with the predominance of irregular behavior near the destroyed invariant loop and the presence of a regular region near the center.

A careful analysis of the Poincaré sections reveals that chaos increases with $`n`$, which is natural since the major contribution of the perturbation comes from $`v>1`$. In Figs. 9(a) and 9(b) we show the Poincaré sections for $`n=2`$ and $`n=5`$, respectively, where $`\beta =1/25`$, $`\mathrm{\Omega }=1`$, and $`\epsilon =0.003`$. As a function of $`\mathrm{\Omega }`$, the most chaotic behavior occurs for frequencies of the same order of the angular frequencies of the orbits in the local minimum of $`V_{eff}`$, from $`(4/3)^2`$ for $`\beta =1/16`$ to $`1`$ for $`\beta =0`$. It is a reasonable result since for small $`\mathrm{\Omega }`$ the system is almost autonomous, and for large frequencies the particle feels only an average of the multipolar motion that goes to zero for $`\mathrm{\Omega }\mathrm{}`$. In Figs. 9(a), 9(c), and 9(d) we show Poincaré sections for $`\mathrm{\Omega }=1`$, $`\mathrm{\Omega }=0.1`$, and $`\mathrm{\Omega }=5`$, respectively, where $`\beta =1/25`$, $`n=2`$, and $`\epsilon =0.003`$. Also, the chaos increases strongly with the “size” of the perturbation $`\epsilon `$. In Figs. 9(a) and 9(e) we show the Poincaré sections for $`\epsilon =0.003`$ and $`\epsilon =0.01`$, respectively, where $`\beta =1/25`$, $`n=2`$, and $`\mathrm{\Omega }=1`$.

The form of the orbits change with $`\beta `$, since there is a change in the destroyed loop, but the relative chaotic area of the section is almost independent of the value of $`\beta `$. In Figs. 9(a), 9(f), and 9(g) we show the Poincaré sections for $`\beta =1/25`$, $`\beta =0`$ and $`\beta =1/16`$, respectively, where $`n=2`$, $`\mathrm{\Omega }=1`$ and $`\epsilon =0.003`$. In particular, the case $`\beta =1/16`$ presents two unstable periodic orbits and no higher chaotic behavior seems to be associated to it.

## V Final remarks

Making use of a compactification procedure and extensions of the standard Melnikov method and Smale-Birkhoff homoclinic theorem, we have detected transverse homoclinic and heteroclinic orbits and Smale horseshoe in periodically perturbed monopole $`+`$ quadrupolelike potential. Moreover the compactified coordinates have proved useful to work out the Poincaré section method.

Although we have studied only one particular system, we stress that the approach is general and can be applied to a large class of open systems. The mathematical results depend essentially on the presence of the monopolar term in the unperturbed equations. Thus the same approach and with the same compactification applies to every inner-multipolar expansion model with nonzero monopolar contribution. In such cases we can allow non-Hamiltonian perturbations of the form $`f(u,t,\epsilon )`$ in the equation corresponding to Eq. (4), where $`f`$ is an analytic function of its variables, periodic in $`t`$, satisfying $`f(u,t,0)=0`$ and of order $`O(u^6)`$.

These systems represent very general situations. In the gravitational case, for instance, they can model a potential due to a mass distribution moving periodically with reflection symmetry. This includes all planar potential written as Fourier expansion in time together with an inner-multipolar expansion in space variable. A specific example is provided by a Newtonian binary system perturbed by gravitational radiation. This system models the long-term dynamical evolution of binary systems of stars due to the emission and absorption of gravitational radiation chicone .

In electromagnetic systems we can have a different situation. The possible absence of monopolar contribution creates a difficulty in transforming the map equations into the analytic form (14)-(15). The solution may be found using the Cassasayas-Fontich-Nunes’s result, which establishes the Melnikov method for systems with parabolic degenerated saddle points fontich . A suitable compactification to treat the electromagnetic case is under study.

ACKNOWLEDGMENTS

The authors are thankful to CNPq and Fapesp for financial support.


format: text

text:
# 1 Introduction

## 1 Introduction

A major motivation for the CERN Large Hadron Collider (LHC) experiments is to shed light into the mechanism of electroweak symmetry breaking. As the LEP precision measurements favour a relatively light Higgs boson, its decay into a pair of photons is a prominent search channel which deserves intense consideration from the theoretical side. The inclusion of QCD radiative corrections is mandatory not only for the signal but also for the background if one is aiming towards reliable predictions. In this respect it is important to note that especially at LHC, but also at the Tevatron, prompt photon signals — photons which do not stem from mesonic decays — are in general contaminated by $`\pi ^0`$ events. The photon pairs from the mesonic decay which are very collimated, if the pions have a high momentum compared to the pion mass, may appear as single (fake) photons inside the detector. This background is especially worrisome in collider experiments, such as the Collider Detector at Fermilab (CDF)  and D0  experiments at the Fermilab Tevatron, where the pions cannot all be reconstructed event by event. At the LHC $`\gamma \pi ^0`$ and $`\pi ^0\pi ^0`$ production represent an important fraction of the reducible background to the search of neutral Higgs bosons in the channel $`H\gamma \gamma `$ in the intermediate mass range 80–140 GeV. Since the jet-jet cross section is about eight orders of magnitude larger than the expected Higgs boson signal for a standard model Higgs boson, this reducible background is overwhelming before any selection cut is imposed.

Regarding both issues of prompt diphoton production and Higgs boson search, various cuts, especially isolation cuts in transverse energy have to be imposed to the candidate events, in order to reject most of this background. Yet, the actual fraction of this background which passes these cuts is not negligible, and it is not well known. Up to now its estimations rely on Monte Carlo simulations encoding the partonic subprocesses at the leading logarithmic (LL) accuracy only. We may thus worry about large higher order corrections which may be not accurately reproduced by parton showers. Furthermore, the hadronisation models matched with the partonic ingredients in these simulations are tuned to describe correctly the bulk of pion production, but they may not be accurate in the kinematical region relevant for the pions which pass the cut. A tentatively quantitative study of isolated $`\gamma \pi ^0`$ and $`\pi ^0\pi ^0`$ production, confronted to a careful analysis of data collected at the Tevatron could improve this situation. The motivation of this article is to provide a study of full next-to-leading order (NLO) accuracy, based on the computer code DIPHOX as an attempt in this direction.

The computer program DIPHOX on which our study relies is a Monte Carlo code of partonic event generator type. This code has been designed to describe the production of pairs of particles in hadronic collisions, accounting for all contributing partonic processes at full NLO accuracy. These particles can be prompt photons or hadrons, in particular neutral pions. This code is flexible enough to accommodate various kinematic or calorimetric cuts. Especially, it allows to compute cross sections for isolated $`\gamma \pi ^0`$ and $`\pi ^0\pi ^0`$ pairs, for any infrared and collinear safe isolation criterion which can be implemented at the partonic level. The physical content and schematical description of DIPHOX has been given in Ref. , where we have concentrated on the case of diphoton production . In Ref.  we recalled briefly how to use this code also in the case of $`\gamma \pi ^0`$ and $`\pi ^0\pi ^0`$ production, and we focused on the phenomenology of dipion production at fixed target energies, with a comparison to the data of the E706 experiment  at the Tevatron. A similar study was also presented in . The aim of the comparison was to confront the theoretical framework on which DIPHOX relies to existing dipion data in order to validate this framework before transporting this knowledge to collider energies  using the factorisation property of perturbative QCD.

The main goal of the present article is to provide a NLO prediction for the reducible background to Higgs boson searches in the mass range 80–140 GeV at the LHC. In section 2 we present estimates for distributions which are of interest regarding the Higgs boson search, such as the distributions of invariant mass and transverse momentum of a pair. In section 3 we will extensively discuss the uncertainties of our next–to–leading order approach induced by severe experimental cuts. Finally section 4 gathers our conclusions and some outlook.

## 2 Predictions for the $`\gamma \pi ^0`$, $`\pi ^0\pi ^0`$ background at the LHC

In this section we discuss NLO predictions for the $`\gamma \pi ^0`$ and $`\pi ^0\pi ^0`$ reducible background to the neutral Higgs boson search channel $`H\gamma \gamma `$ in the intermediate mass range 80–140 GeV at the LHC. All NLO curves where computed with the DIPHOX code, version (1.2) 111In this version a bug was removed compared to the first version described in . In that paper the main issue was the prediction of photon pairs. Only fragmentation contributions in the case of severe isolation cuts were affected by the bug. These are numerically irrelevant for the diphoton rate predictions.. For all figures we are using the MRST2  parton distribution functions, and, if not stated differently, the KKP fragmentation functions . For all plots given below the following cuts on the transverse momenta, pseudo rapidities of the photons/pions and their invariant mass were applied:

$`p_T(\gamma _1,\pi _1)`$ $`>`$ $`40\text{GeV},`$
$`p_T(\gamma _2,\pi _2)`$ $`>`$ $`25\text{GeV},`$
$`|y(\gamma ,\pi )|`$ $`<`$ $`2.5,`$
$`80\text{GeV}`$ $`<`$ $`m_{\pi \pi ,\gamma \pi ,\gamma \gamma }<140\text{GeV}.`$ (1)
In addition, isolation cuts are imposed to reject events which are accompanied by a certain amount of hadronic energy, as this is typically not the case for photonic Higgs decay. A photon is said to be isolated if, in a cone in rapidity and azimuthal angle about the photon direction, the amount of deposited hadronic transverse energy $`E_T^{had}`$ is smaller than some value $`E_{Tmax}`$ fixed by the experiment:

$$E_T^{had}E_{Tmax}\text{inside}\left(yy_\gamma \right)^2+\left(\varphi \varphi _\gamma \right)^2R^2$$
(2)
The actual experimental criteria, which are defined at the detector level, are much more complicated than the simple one given by Eq. (2), see . They cannot be implemented in a partonic calculation. The schematical criterion of Eq. (2) is a modelization of their effects at the parton level. The topic of isolation of photons based on the above criterion is extensively discussed in the literature . We stress that there is a subtle interplay between higher order corrections and the experimental cuts. These subtleties will be discussed in the next section.

In Fig. 1 we show the invariant mass spectra of $`\gamma \gamma `$, $`\gamma \pi ^0`$ and $`\pi ^0\pi ^0`$ pairs for a very loose isolation criterion of $`E_{Tmax}=100`$ GeV in a cone of $`R=0.4`$ at NLO. Note that the fully inclusive case would be plagued by uncertainties stemming from the low–$`z`$ regime of the pion fragmentation functions (see discussion below). To avoid this source of uncertainty we do not present figures showing the inclusive rate. The plot shows the importance of the reducible $`\pi ^0\gamma `$, $`\pi ^0\pi ^0`$ backgrounds before isolation cuts. The $`\pi ^0\gamma `$ component is more than an order of magnitude above the $`\gamma \gamma `$ background. The $`\pi ^0\pi ^0`$ background is of the order of several hundreds of pb/GeV. Tightening the isolation requirement by imposing $`E_{Tmax}=10`$ GeV changes the production rate of these components. This is shown in Fig. 2. We remind that the very narrow Higgs signal will be a sharp peak — its width is determined by the detector resolution of around 1 GeV — ranging from 50 to 100 fb/GeV in the given mass window . The $`\pi ^0\pi ^0`$ background is reduced considerably whereas the $`\gamma \pi ^0`$ background is still of the order of the irreducible diphoton background. Concerning the prompt photon pairs, the used isolation criterion reduces its rate by about a factor of two, because the fragmentation component of the prompt photon pair is vetoed by isolation cuts. Note that for the given isolation criteria not only direct photon pairs remain. For $`E_{Tmax}=10`$ GeV almost 20 per cent of the total rate still come from fragmentation in the given plot for the scale choice used. Fig. 3 shows the suppression of the irreducible diphoton background when isolation is imposed. To suppress photons from fragmentation more efficiently harder isolation cuts than the ones shown in the plot have to be imposed. A remaining uncertainty for the prompt photon pairs comes from higher order corrections to the box contribution. These corrections have been computed recently  but are not yet implemented in DIPHOX. To give an idea of the rejection factors due to isolation we have plotted the invariant mass distributions of the photon–pion pair for maximally allowed hadronic energies of $`E_{Tmax}=100,20`$, and 10 GeV in a cone of $`R=0.4`$ in Fig. 4. Due to the collinear fragmentation model used the hadronic energy is mainly located along the fragmenting particle. This leads to the fact that there is no significant $`R`$ dependence present. Only genuine higher order $`23`$ matrix elements can induce such a dependence and we only state that they are on the per cent level. In the case of photon–pion production two mechanisms for the production of the photon are possible. The photon can be produced directly in a hard interaction or through fragmentation. In the case of isolation the latter component is reduced by an order of magnitude as can be seen in Fig. 5. For completeness we provide with the plot for the invariant mass distribution of pion pairs, Fig. 6. In a realistic experimental situation the pion rates can be suppressed further. In the relevant $`p_T`$ range, an extra reduction factor of 2 to 3 per pion can be expected in fine-grained detectors with recognition of isolated $`\pi ^0`$ as two overlapping showers .

Another spectrum of relevance is the transverse momentum distribution of the pair produced particles, $`q_T=|\stackrel{}{p}_{T\mathrm{\hspace{0.33em}1}}+\stackrel{}{p}_{T\mathrm{\hspace{0.33em}2}}|`$. In Fig. 7 this distribution is shown for two isolation criteria. Recall that the full pion photon rate is built up out of the direct component and the fragmentation component. At colliders which operate at lower energies the fragmentation component is typically only a fraction of the full rate. Through isolation the fragmentation component dies out quickly. Fig. 8 gives an idea of the suppression factors due to isolation. Note that this theoretical curve is not reliable in the first few bins, as multiple soft gluon emission alters the divergent behaviour of the fixed order calculation. Near $`q_T=0`$ the partonic cross section diverges logarithmically. This divergence is diluted by convoluting the partonic cross section with the pion fragmentation function. Also at $`q_T=E_{Tmax}`$ there is a critical point inside the physical spectrum (see  for a discussion). The large bin width prevents from observing this singular behaviour. To deal with these issues in a rigorous way is beyond the scope of this study222The corresponding resummations to be done are technically very involved and no satisfying procedure is defined yet..

Hard isolation cuts test fragmentation functions at high $`z`$ values, where current parametrisations are not at all constrained by experimental data. This problem is visualised in Fig. 9 where we show the photon–pion invariant mass distribution calculated with both, the KKP  and the older BKK  fragmentation functions. The KKP parametrisation leads to a higher prediction than the BKK one. The ratio of the two curves (KKP/BKK) is about 3/2. For dipion observables the situation is much worse. In Fig. 10 the discrepancy between the BKK and KKP parametrisations in the critical region is shown. One observes that at high $`z`$ values the fragmentation gluon–to–pion is strongly suppressed. Focusing on the quark–to–pion fragmentation, one sees that for high $`z`$ values the parametrisations differ considerably. This will be discussed in more detail below.

We have to comment now on the remaining theoretical uncertainty due to higher order corrections. As will be pointed out below the large $`z`$ region is dangerous for the predictive power of NLO calculations because terms are produced which are divergent as $`z1`$. In Fig. 11 , we plot LO and NLO predictions333LO refers here to the Born matrix elements whereas NLO parton densities and fragmentation functions are used. for the invariant mass distribution of the photon–pion pair for several isolation criteria. Note that the ratio of both curves is pretty stable with decreasing $`E_{Tmax}`$. The ratio of both curves lies around 1.7, a typical value for a NLO correction. We just mention that the situation worsens as $`E_{Tmax}`$ is reduced further which indicates the presence of large NLO contributions induced by pressing $`z`$ to one. As will be pointed out below there are two reasons for that. Phase space regions which are excluded at LO are populated at NLO and become dominant, and large logarithms of $`(1z)`$ are introduced. To have a proper control on the latter fixed order effect, a resummation of the relevant contributions is needed. For the isolation criteria used in the given plots we conclude from the NLO to LO ratios that these effects are not yet dominant. Still, for harder isolation criteria the predictivity of our fixed order approach has to be critically investigated.

As the isolation cuts spoil compensations between scale dependent LO and NLO contributions the scale dependence remains important. To give an idea of the uncertainties we have plotted the variation of the invariant mass spectrum of $`\pi ^0\gamma `$ pairs under the change of renormalization ($`\mu `$), factorisation ($`M`$) and fragmentation scale ($`M_f`$) for a strict isolation criterion in Fig. 12. By multiplying our standard scale choice by 1/2 and 2 one finds variation of about +40 per cent to –30 per cent. This variation is reduced if one uses less restrictive experimental cuts.

The importance of NLO calculations is best visible in distributions which are restricted by LO kinematics. No LO parton shower can provide for a correct description of events far away from the LO kinematics. A relevant observable to exemplify this point is the $`q_T`$ distribution of a photon–pion pair. Increasing the isolation criterion, suppresses the LO phase space more and more, such that finally the whole tail of the distribution is a pure NLO effect, see Fig. 13. For completeness we show the effect of isolation on this observable in Fig. 8. For a comparison of the shape of our DIPHOX prediction to PYTHIA see .

## 3 Experimental cuts and higher order corrections

To understand the reliability of the presented predictions it is necessary to discuss the interplay between experimental cuts and theoretical NLO results. QCD cross sections are stabilised by the inclusion of higher order corrections, as scale dependencies are reduced. This works generally the better the more inclusive the observables are. Yet, the collider experiments at the Tevatron and the forthcoming LHC do not perform inclusive photon measurements. The experimental selection of prompt photons requires isolation cuts. We always use the cone isolation criterion introduced in Eq. (2). For the hadronic energy propagating into the isolation cone two sources can be distinguished in an experimental setup. First there is a component which stems from the actual partonic, hard reaction which may be mapped to certain Feynman diagrams – this component is calculable inside a QCD calculation to a given order. And second there is a component coming from multiple-interactions and other events in the same particle bunch (pile–up). One has to be aware of the fact that an experimentally used $`E_{Tmax}`$ is thus always bigger than a partonic value used in our study. It is even so that the experimental values $`E_{Tmax}`$ are so small that a large fraction of the allowed hadronic transverse energy accompanying the candidate photon may come from such events. This leaves not much phase space for partonic fragmentation reactions. In the following we mean with $`E_{Tmax}`$ only the hadronic energy stemming from the hard interaction. The sensitivity on $`E_{Tmax}`$ defined in Eq. (2) is best seen by looking at the fragmentation variable $`z`$, which is the momentum fraction of a hadron, $`p_{Th}`$ ($`h=\pi ,\gamma `$), and its parent parton $`p_T`$. At LO the kinematical constraints on this variable are induced by the experimental cuts through the collinearity assumption in the fragmentation model, $`p_T=zp_T+(1z)p_T=p_{Th}+E_T^{had}`$. One finds

$`z>z_{min}={\displaystyle \frac{p_{Tmin}}{p_{Tmin}+E_{Tmax}}}`$ (3)
where $`p_{Tmin}`$ is the experimental cut on $`p_{Th}`$. Evidently, the variable $`z`$ is pressed towards one, as the ratio $`E_{Tmax}/p_{Tmin}`$ decreases, means for hard isolation. At NLO, where extra radiated partons are present, this bound is even more severe as the veto is more easily fulfilled. The bulk of the pions which are produced in collider events correspond to values of $`z`$ substantially smaller than 1. Hard isolation tests only the very end of the $`z`$ spectrum. The question is, wether we can reliably estimate the production rate of these high–$`z`$ pions, and especially of hard $`\gamma \pi ^0`$ and $`\pi ^0\pi ^0`$ pairs. We find three aspects important in this respect, which are discussed in more detail in the following.

### 3.1 Reliability of fragmentation functions

The usual fixed order perturbative QCD calculations of hadro production – as e.g. the one coded in DIPHOX – rely on a model of independent collinear fragmentation through fragmentation functions of a parton of species $`j`$ (quark of any flavour, or gluon) into a pion, $`D_{j\pi ^0}(z,M_f^2)`$, $`z`$ is defined above. Let us first comment on the uncertainties at small $`z`$, where this treatment is questionable. In particular the fixed order calculations do not resum the $`\mathrm{ln}(1/z)`$ arising order by order in perturbation theory, which become large in the small $`z`$ region. The latter are associated with coherence effects in the multiple emission of soft gluons, which control the inclusive production of small $`z`$ pions . Therefore this usual framework is notoriously inadapted to account for the production of these small $`z`$ pions. Moreover the dominant $`z`$ values at colliders are smaller than for fixed targets, as a consequence of scaling violations. Therefore the erroneous treatment of the fragmentation process at small $`z`$ will plague the estimates of inclusive dipion observables at colliders. This is the reason why the fits of independent fragmentation functions $`D_{i\pi ^0}(z,M_f^2)`$ to experimental data presented in the literature , and which we use in this work, start above some minimum value for $`z`$, e.g. $`z_{min}0.1`$. The independent fragmentation framework is expected to account correctly for the production of pions with larger $`z`$ values. The requirement of a loose isolation criterion such as, for example $`E_{Tmax}`$ 100 GeV for pions with $`p_{T\pi ^0}`$ 25 GeV imposes that $`z0.2`$ over the whole $`p_{T\pi ^0}`$ range of interest. To avoid uncertainties due to small $`z`$ values, as was pointed out already, we do not provide estimates for inclusive observables, but only for observables corresponding to pions submitted to isolation cuts.

The behaviour near the end point, $`z1`$, is responsible for most of the uncertainties that plague the present study. Indeed, the use of very tight isolation cuts such as those used by the collider experiments select the vicinity of the elastic boundary $`z1`$. For example, in the Run Ib of the Tevatron, for prompt photon measurements, CDF required $`E_{Tmax}=2`$ GeV in a cone $`R=0.7`$. For $`p_T15`$ GeV, this means $`z_{min}(p_T)0.9`$. The value $`E_{Tmax}=2`$ GeV required by the CDF experiment includes also the deposition of hadronic transverse energy coming from underlying events, while the latter are ignored in our partonic calculation as noted above. Consequently, the actual effective value of $`E_{Tmax}`$ which should be used in a partonic calculation should be even smaller, the effective $`z_{min}`$ would be even larger. However the fragmentation functions $`D_{j\pi ^0}(z,M_f^2)`$ are fitted to data corresponding essentially to $`z0.7`$, i.e. much below the relevant region for the present study. Moreover, the data used in  concern hadronic final states in $`e^+e^{}`$ annihilation. The gluon–to–pion fragmentation function is poorly constrained at moderate and large $`z`$ by this process because the subprocesses involving gluons in $`e^+e^{}`$ annihilation appear at higher order in the perturbative expansion in powers of $`\alpha _s`$, while subprocesses involving outgoing gluons contribute at lowest order in hadronic collisions. Although the contributions from gluon fragmentation are not expected to dominate, large uncertainties on their magnitude may still affect quantitative predictions. In order to illustrate how poor the present-day control over the large $`z`$ region is, it is instructive to compare the extrapolations of the KKP and BKK sets of parton-to-$`\pi `$ fragmentation functions. The KKP set is the most recent one, and it is constrained by more data, especially the LEP1 data at the $`Z`$ peak, in the region $`0.1z0.7`$. These two sets use different analytic ansätze. These different ansätze do not differ very much in the $`z`$ range in which they are most constrained by the $`e^+e^{}`$ data of course, but their respective extrapolations turn out to differ substantially from each other above $`z0.7`$, especially between $`z=0.9`$ and 1, as shown in Fig. 10. This is true for the gluon, but also, though to a lesser extend, for the $`u`$ and $`d`$ fragmentation functions. This translates into a large uncertainty for the predictions of the reducible Higgs boson background in the channel $`H\gamma \gamma `$ in the mass range 80–140 GeV. Clearly a determination of fragmentation above $`z0.7`$ from data sensitive to this range is required for any further accurate prediction. Such an experimental input could in principle be provided by data on inclusive $`\pi ^0`$ production at fixed target energy 444In Ref. , the comparison between E706 data and BKK fragmentation functions reveals an overall normalisation discrepancy. However, as noticed in Ref. , this problem seems to be largely reduced by using the KKP set., as well as measurements of isolated pions by the collider experiments CDF and D0, in association with their prompt photon measurements.

It is worth mentioning that the present inaccuracy of the description of fragmentation at large $`z`$ is not only a problem confined to the type of perturbative QCD calculations which are discussed in this article. They also potentially affect the fragmentation stage of the “full hadronic event generators” such as PYTHIA  and HERWIG  as well. Indeed, the respective fragmentation ingredients of these event generators have been tuned to describe accurately the bulk of the data, especially at LEP, whereas the tails of the distributions near the end point have, to our knowledge, never been analysed in detail. In addition, the parton showers implemented in these event generators provide an effective summation of the soft gluon effects near $`z1`$ at a leading logarithmic accuracy; on the other hand they miss a part of the higher order corrections to the partonic hard subprocesses which may be large. The respective approximations in the “partonic calculations” on one hand, the full event generators on the other hand, being different, it would be instructive to compare the results given by both approaches.

### 3.2 Large $`z`$ values and the need for end-point summation

As the production of isolated pions at colliders and the one of inclusive pions at fixed target experiment, deal with $`z`$ values near the end point $`z=1`$, the theoretical calculation of the partonic subprocesses is confronted with the appearance of powers of large logarithmic terms $`\mathrm{ln}(1z)`$ to any fixed order in perturbation theory. These large logarithms have in principle to be summed to all orders in order to improve – if not restore – the predictive power of the partonic calculation. The general framework to resum infrared logarithms in hadronic collisions has been settled in , however further work is needed to treat the endpoint summation in the case of one particle production, even more for pair production, where the kinematics is more constrained. In the present study, we investigate what can be done without such summations, relying on the extrapolations of the fragmentation functions available in the literature, and examining critically the limitations of our approach.

To be explicit let us consider the process $`PP\gamma \pi ^0+X`$ in detail. The photon can be either direct or from fragmentation. The second case is considerably suppressed if one imposes severe isolation; thus we deal with the direct process only. At leading order there are two types of partonic processes: $`qgq\gamma `$555Here $`q`$ stands for quark or antiquark., $`q\overline{q}\gamma g`$. At small $`x`$ values typical for the LHC the first is dominant. Radiative corrections give rise to initial state and final state singularities. Only the latter are of relevance for the point to be made here. They arise from the splitting of the outgoing parent quark into a quark and a gluon. Both may fragment into a pion but, again, for severe isolation one may focus on the dominant $`q\pi ^0`$ fragmentation only. Thus we are looking at the case where a collinear gluon is radiated of a hard quark, indicated by the subscript $`g//q`$, which subsequently fragments into a pion. Following Ref.  the finite remnant of the collinear approximation of this correction is of the form

$`\sigma (PPqg\gamma )|_{g//q}`$ $`=`$ $`{\displaystyle 𝑑y_\pi 𝑑y_\gamma 𝑑p_{T\pi }\underset{x_{min}}{\overset{1}{}}\frac{dx}{x}\underset{z_{min}}{\overset{1}{}}\frac{dz}{z}\frac{\alpha _s(\mu )}{2\pi }}`$ (4)
$`\times D_{\pi ^0/q}(x,M_f^2){\displaystyle \frac{F_{q/P}(x_1,M^2)}{x_1}}{\displaystyle \frac{F_{g/P}(x_2,M^2)}{x_2}}|_{qgq\gamma }^{LO}|`$
$`\times \left[2\left({\displaystyle \frac{\mathrm{ln}(1z)}{1z}}\right)_+a_{qq}^{(4)}(z)+\mathrm{ln}\left({\displaystyle \frac{p_{T\pi }^2}{x^2z^2M_f^2}}\right)P_{qq}^{(4)}(z)+\mathrm{}\right]`$
$`P_{qq}^{(4)}(z)`$ $`=`$ $`{\displaystyle \frac{a_{qq}^{(4)}(z)}{(1z)_+}}+{\displaystyle \frac{3}{2}}C_F\delta (1z)`$
$`a_{qq}^{(4)}(z)`$ $`=`$ $`C_F(1+z^2)`$
Here $`y_\pi ,p_{T\pi }`$ are the rapidity and transverse momentum of the pion, $`y_\gamma `$ the rapidity of the photon, $`xz`$ the momentum fraction of the pion with respect to the parent quark, $`p_{T\pi }=xzp_{Tq}`$, and $`z`$ the momentum fraction of the gluon with respect to the parent quark, $`p_{Tg}=zp_{Tq}`$. In the given collinear approximation the momentum fractions of the initial partons are

$`x_1`$ $`=`$ $`{\displaystyle \frac{p_{Tq}+p_{Tg}}{\sqrt{S}}}[\mathrm{exp}(y_\pi )+\mathrm{exp}(y_\gamma )]`$
$`x_2`$ $`=`$ $`{\displaystyle \frac{p_{Tq}+p_{Tg}}{\sqrt{S}}}[\mathrm{exp}(y_\pi )+\mathrm{exp}(y_\gamma )]`$
$`z_{min}`$ $`=`$ $`x_{min}/x`$
The value of $`x_{min}`$ will be discussed below, as it is the sensitive variable to isolation cuts. The dots stand for terms which are less dominant as $`z1`$. They are not necessarily small. By writing $`a_{qq}^{(4)}(z)/z=C_F[(1z)^2/z+2]`$, one sees that the first term in the bracket only leads to non logarithmic corrections. Thus, approximating the bracket by 2 and using that for $`z1`$, $`\mathrm{ln}(z)=𝒪(1z)`$, the $`z`$ integration can be done explicitly and one finds

$`\sigma (PPqg\gamma )|_{g//q}`$ $`=`$ $`{\displaystyle 𝑑y_\pi 𝑑y_\gamma 𝑑p_{T\pi }\underset{x_{min}}{\overset{1}{}}\frac{dx}{x}\frac{\alpha _s(\mu )C_F}{\pi }}`$ (5)
$`\times D_{\pi ^0/q}(x,M_f^2){\displaystyle \frac{F_{q/P}(x_1,M^2)}{x_1}}{\displaystyle \frac{F_{g/P}(x_2,M^2)}{x_2}}|_{qgq\gamma }^{LO}|`$
$`\times [\mathrm{ln}(1x_{min}/x)^2+\mathrm{ln}\left({\displaystyle \frac{p_{T\pi }^2}{x^2M_f^2}}\right)\mathrm{ln}(1x_{min}/x)+\mathrm{}]`$
The important point is that the behaviour of these terms changes drastically as one goes from the inclusive case, where $`x>x_{min}=2p_{T\pi }/\sqrt{S}`$ is much smaller than one, to the case of severe isolation which press $`x`$ to values near one. Using the fact that the fragmentation function behaves as $`D_{\pi ^0/q}(x)N(1x)^\beta `$ with $`\beta 1`$, one can write

$`{\displaystyle \frac{D_{\pi ^0/q}(x)}{x}}={\displaystyle \frac{D_{\pi ^0/q}(x_{min})}{x_{min}}}\left(1{\displaystyle \frac{xx_{min}}{1x_{min}}}\right)+𝒪(1x_{min},1\beta )`$ (6)
and perform the integral over $`x`$. Focusing on the double logarithmic terms, one finds

$`\sigma (PPqg\gamma )|_{g//q}`$ $`=`$ $`{\displaystyle 𝑑y_\pi 𝑑y_\gamma 𝑑p_{T\pi }\frac{\alpha _s(\mu )C_F}{\pi }}`$ (7)
$`\times (1x_{min}){\displaystyle \frac{D_{\pi ^0/q}(x_{min},M_f^2)}{x_{min}}}{\displaystyle \frac{F_{q/P}(x_1,M^2)}{x_1}}{\displaystyle \frac{F_{g/P}(x_2,M^2)}{x_2}}|_{qgq\gamma }^{LO}|`$
$`\times {\displaystyle \frac{1}{2}}[\mathrm{ln}(1x_{min})^2+\mathrm{ln}\left({\displaystyle \frac{p_{T\pi }^2}{M_f^2}}\right)\mathrm{ln}(1x_{min})+\mathrm{}]`$
the dots stand for single logarithmic terms and terms which are suppressed by an extra factor of $`(1x_{min})`$ and $`(1\beta )`$. Note that outside the bracket, up to the factor $`\alpha _s(\mu )C_F/(2\pi )`$, one has just the approximation for the leading order contribution. This means that the bracket itself times this factor could be taken as a good approximation for the higher order part of a $`K`$ factor666Note the way this $`K`$ factor depends on the experimental cuts applied and on the scale choice. This is a nice example of the fact that for non inclusive quantities the naive use of $`K`$ factors is not at all adequate to take into account higher order corrections..

$`K=1+{\displaystyle \frac{\alpha _s(\mu )}{\pi }}[\mathrm{ln}(1x_{min})^2+\mathrm{ln}\left({\displaystyle \frac{p_{T\pi }^2}{M_f^2}}\right)\mathrm{ln}(1x_{min})+\mathrm{}]`$ (8)
In here a factor of two is included from the second process with initial partons $`\overline{q}g`$. To get a feeling for the typical scales in the case of isolation one may try to choose $`M_f`$ such that the large logarithms compensate for a $`p_{T\pi }p_{Tmin}`$. This can be achieved by a (optimal) scale choice

$`M_f^{opt}=p_{T\pi min}\sqrt{{\displaystyle \frac{E_{Tmax}}{p_{Tmin}+E_{Tmax}}}}`$ (9)
Plugging in $`p_{Tmin}=25`$ GeV 777Due to the NLO kinematics $`p_{T\gamma }p_{T\pi }`$. For asymmetric $`p_T`$ cuts the photon will always be the particle restricted by the larger $`p_T`$ cut. and $`E_{Tmax}=10`$ GeV one finds $`M_f^{opt}10`$ GeV. This motivated us to choose for $`M_f`$ in the figures presented not the canonical choice $`M_{\pi \gamma }/2`$ but rather $`M_{\pi \gamma }/8`$. Note that for a much higher scale choice the logarithms become constructive. For the upper choice of parameters and using $`M_f=2p_{Tmin}`$ arbitrarily large K factors are induced in the limit $`x_{min}1`$. The fact that the optimal scale in Eq. (9) goes to zero, as $`E_{Tmax}`$ goes to zero indicates that the cone isolation criterion is not infrared safe in this limit. In other words, if the experimental value of $`E_{Tmax}`$ is already saturated by events from the hadronic environment such that there is no phase space left for the hard parton, fixed order calculations are not applicable.

It should be stressed that the whole discussion was neglecting constant terms and sub-leading logarithms. Thus it should only give a rough idea of the much more involved, complete situation. Using the same approximations the exercise can be worked through keeping also non–leading terms. The result is that one encounters in realistic situations additive contributions which are of the same size as the considered ones. All this deserves further investigations especially in the light of a possible resummation of the large logarithms and the role of these terms.

### 3.3 Hard isolation and asymmetric $`p_T`$ cuts

In LHC studies concerning Higgs boson search the used $`p_T`$ cuts are commonly asymmetric. We want to make the point here that this jeopardizes the predictive power of NLO calculations. The interplay between such asymmetric cuts and hard isolation criteria leads at LO to vetoed phase space regions which are filled by NLO processes only. The effect is visualised in Fig. 14, where the LO and NLO DIPHOX prediction for the dipion $`p_T`$ spectrum is plotted for symmetric, both $`p_T>25`$ GeV, and asymmetric cuts, $`p_{T1}>40`$ GeV, $`p_{T2}>25`$ GeV. In the case of asymmetric cuts one observes that the LO component is suppressed for low $`p_T`$ values. This can be understood as follows:
Due to transverse momentum conservation, one has

$$p_{T\mathrm{\hspace{0.17em}1}}+E_{T\mathrm{\hspace{0.17em}1}}^{had}=p_{T\mathrm{\hspace{0.17em}2}}+E_{T\mathrm{\hspace{0.17em}2}}^{had}$$
(10)
where the pion 1 is assumed to have higher $`p_T`$ than pion 2. This means for the transverse momentum of the pair, $`q_T`$ that

$$0q_T=p_{T\mathrm{\hspace{0.17em}1}}p_{T\mathrm{\hspace{0.17em}2}}=E_{T\mathrm{\hspace{0.17em}2}}^{had}E_{T\mathrm{\hspace{0.17em}1}}^{had}E_T^{max}$$
(11)
and because $`p_{T\mathrm{\hspace{0.17em}1}}>p_{T\mathrm{\hspace{0.17em}1}min}`$ (= 40 GeV in our case)

$$p_{T\mathrm{\hspace{0.17em}2}}p_{T\mathrm{\hspace{0.17em}1}}E_T^{max}p_{T\mathrm{\hspace{0.17em}1}min}E_T^{max}$$
(12)
This means that in the case of $`p_{T\mathrm{\hspace{0.17em}1}min}E_T^{max}>p_{T\mathrm{\hspace{0.17em}2}min}`$ one has a phase space restriction coming from the isolation criterion. Note that a saturation of the bound, $`p_{T\mathrm{\hspace{0.17em}2}}=p_{T\mathrm{\hspace{0.17em}1}min}E_T^{max}`$ is only possible if $`E_{T\mathrm{\hspace{0.17em}1}}^{had}0`$, i.e. $`z_11`$. But in this limit the fragmentation function goes to zero and suppresses the cross section additionally. This will not happen in the case of photon–pion production where the same formulas apply with $`E_{T\mathrm{\hspace{0.17em}1}}^{had}=0`$. Eq. (12) is also true and leads to a restriction of the phase space but there is no additional suppression. The problem aggravates for lower $`E_{Tmax}`$ values but becomes less relevant when $`E_{Tmax}>p_{Tmin1}p_{Tmin2}`$ and vanishes in the case of symmetric cuts $`p_{Tmin1}=p_{Tmin2}`$. The hole in the $`p_T`$ spectrum will be filled with events from higher order corrections. As the low $`p_T`$ region is energetically preferred, one gets numerically important NLO contribution to the cross section in the dipion case. The ratio between LO and NLO can become very big due to this effect. As pointed out, in photon–pion production the LO order contribution is less suppressed and the effect is milder. In the $`q_T`$ distribution shown in Fig. 8, one observes that only the first bin contains the LO contribution due to the same kinematical restriction mentioned above. On the other hand the NLO contribution is basically zero there (even below the plotted range). This is because pions with $`p_T`$’s as large as the photon $`p_T`$ are disfavoured. The steep fall of the quark–to–pion fragmentation functions makes it very unlikely to find a photon–pion pair with $`p_{T\gamma }p_{T\pi }>40`$ GeV which is needed to have $`q_T0`$.

From the discussion one sees that symmetric cuts are preferable from theoretical side for processes which have to be suppressed by isolation. Whereas this effect is for sure important for dipion production the effect is numerically sizable in the photon–pion case only for isolation criteria which are harder than the ones given in the presented figures. As can be inferred from Fig. 11, the ratio of the NLO to LO curves does not explode because the discussed phase space effect is not yet numerically sizable. Still we want to make the point that our NLO analysis suggest to choose symmetric transverse momentum cuts as theoretical higher order predictions tend to be more reliable in that case.

## 4 Conclusions and outlook

We have presented a full next to leading order study of the production of $`\pi ^0\gamma `$ as well as $`\pi ^0`$ pairs with a large invariant mass relying on the computer code DIPHOX. We have presented estimates for invariant mass and transverse momentum of the pair ($`q_T`$) distributions taking into account realistic experimental cuts. The NLO corrections are large. Theoretical calculation improves the scale dependences in the case of inclusive observables. On the other hand, the results remain plagued by large uncertainties even in the NLO approximation, if realistic isolation criteria are applied. One source of large uncertainties is our poor present knowledge of the parton–to–pion fragmentation functions at large $`z`$. The estimates which we have presented here rely on extrapolations of analytical ansätze for these fragmentation functions for values of $`z`$ above 0.8, i.e. outside the region where they were fitted to the data. The spreading of the extrapolations of different sets of parametrisations available in the literature for $`z`$ above 0.8 illustrates this fact. Experimental information may be provided by inclusive measurements of single pions and pions pairs at fixed targets, e.g. as performed by the E706 experiment , as well as by measurements of isolated pions at colliders, as could be performed by the CDF  and D0  experiments.

Apart from this more experimental issue we pointed out that our next-to-leading order approach is still plagued by relatively large uncertainties due to phase space effects related to the different LO and NLO kinematics and also the presence of potentially large logarithms $`\mathrm{ln}(1x_{min})`$, if severe isolation cuts are imposed. Concerning the phase space effect, we discussed that symmetric $`p_T`$ cuts for the pions/photons are preferable from theoretical side. Otherwise NLO corrections are magnified in phase space regions vetoed at lowest order. Our predictions for $`\gamma \pi ^0`$ production seems still reliable, up to scale variations of around 40 per cent, for the used isolation cuts. The same cannot be said for $`\pi ^0\pi ^0`$ production. The region $`z1`$ is probed in a more critical way in this case. Further, predictions are very sensitive to the unconstrained part of fragmentation functions. We have to conclude that the uncertainties for this component of the reducible background for Higgs searches are not under control from theoretical side even at the next–to–leading order level. We only quote that we find large K factors and large scale variations, factors of 3 to 4 respectively, if we use the same isolation cuts as was done for the isolated $`\gamma \pi ^0`$ production case. Fortunately this component to the diphoton background can be suppressed experimentally by separating pions from photons by distinguishing the shapes of the corresponding electromagnetic showers. Reduction factors of 2 to 3 per pion are possible and compensate somewhat the poor knowledge of the pion rates. Still, a more satisfactory theoretical knowledge of the dipion and pion/photon rates at LHC requires a better understanding of fragmentation functions near $`z=1`$, the resummation of large endpoint logarithms of $`(1z)`$, together with the use of experimental cuts which do not destabilise compensations between LO and higher order contributions.

From the discussion above it should be clear that it is a delicate issue to make definite predictions for the $`\gamma \pi ^0`$ and $`\pi ^0\pi ^0`$ component of the $`H\gamma \gamma `$ background. All statements depend strongly on the experimental cuts. Especially severe isolation cuts restrict the predictivity of our result. To make some explicit statements we just want to remind how the $`\gamma \pi ^0`$, $`\pi ^0\pi ^0`$ rates compare to the irreducible $`\gamma \gamma `$ background for a definite isolation criterion, see Fig. 2. As the Higgs boson of the Standard Model is already restricted by the LEP experiments to have a mass of $`M_H115`$ GeV, the $`\gamma \pi ^0`$ rate is about a factor of two below the irreducible $`\gamma \gamma `$ background. Uncertainties due to fragmentation functions and scale variations may be factors of around two. On the other hand additional experimental capabilities based e.g. on calorimetric studies indicate that the rejection factor of pions w.r.t isolated photons can be improved by at least a factor two . These numbers suggest a conservative upper bound on the ratio of the $`\gamma \pi ^0`$, $`\gamma \gamma `$ rates of $`0.5`$. If one tightens the isolation criterion this will be further reduced, though one has to keep in mind that the theoretical uncertainty is getting more important. In this respect we note that a comparison with PYTHIA distributions shows that our rates are lower if compared to default PYTHIA settings . Though we want to stress the relevance of a next–to–leading order analysis, from our discussion one may not expect that the $`\gamma \pi ^0`$ process becomes a dangerous background. As already noted above the uncertainties from fragmentation functions and scale variation for the $`\pi ^0\pi ^0`$ rate are substantial. To be conservative factors as big as ten should be considered. Assuming again a reduction factor of two per pion from pion–photon separation would lead in Fig. 2 to a $`\pi ^0\pi ^0`$ to $`\gamma \gamma `$ ratio of around 1/3 for $`M_{\pi \pi }>115`$ GeV. More severe isolation will suppress this ratio considerably again with the price of even increased uncertainties. Note that the given numbers should only give an idea of the situation, as they depend strongly on the experimental settings. We can tentatively conclude that next–to–leading order corrections should not change the qualitative statement that the $`\pi ^0\pi ^0`$, $`\pi ^0\gamma `$ backgrounds are sufficiently under control at the LHC, although a couple of interesting QCD issues plaguing them still deserve to be addressed.

## Acknowledgements

We thank K. Lassila–Perini, V. Tisserand and E. Tournefier for discussions and comments. This work was supported in part by the EU Fourth Training Programme “Training and Mobility of Researchers”, Network “Quantum Chromodynamics and the Deep Structure of Elementary Particles”, contract FMRX-CT98-0194 (DG 12 - MIHT). LAPTH is a “Unité Mixte de Recherche (UMR 5108) du CNRS associée à l’Université de Savoie”.