text:
# Self-stabilization of high frequency oscillations in semiconductor superlattices by time-delay autosynchronization

## Abstract

We present a novel scheme to stabilize high-frequency domain oscillations in semiconductor superlattices by a time–delayed feedback loop. Applying concepts from chaos control theory we propose to control the spatio-temporal dynamics of fronts of accumulation and depletion layers which are generated at the emitter and may collide and annihilate during their transit, and thereby suppress chaos. The proposed method only requires the feedback of internal global electrical variables, viz current and voltage, which makes the practical implementation very easy.

superlattice, autosynchronization, chaos control

Semiconductor superlattices ESA70  have been demonstrated to give rise to self-sustained current oscillations ranging from several hundred MHz KAW86 ; KAS95 ; HOF96  to 150 GHz at room temperature SCH99h . Various mechanisms with BUE77 ; PAT98 ; SCH02e ; BON02  or without KRO00  the involvement of propagating field domains have been discusssed. In any case, a superlattice constitutes a highly nonlinear system SCH00 , and instabilities are likely to occur. Indeed, chaotic scenarios have been found experimentally ZHA96 ; LUO98b ; BUL99  and described theoretically in periodically driven BUL95  as well as in undriven systems AMA02a . For a reliable operation of a superlattice as an ultra-high frequency oscillator such unpredictable and irregular conditions should be avoided. In principle, synchronization of oscillations in a superlattice by an external signal SCH02b  could be exploited to achieve a desired periodic behavior. However, in reality, the control of the forcing frequency in the ultrahigh range presents substantial technical problems.

Here, we propose a simple self-stabilizing scheme that is especially suitable for semiconductor devices like superlattices. It uses a profound concept of chaos control from nonlinear dynamics and chaos theory. Within this approach, an intrinsically unstable time-periodic motion is stabilized using a simple feedback loop with a time delay PYR92 . This type of control needs only small control forces initially, and they vanish once control has been achieved. A sound advantage is that the oscillation mode to be stabilized need not be known beforehand, in contrast to other chaos control schemes. Rather, a simple delay line leads to autosynchronization of the system. Methods of nonlinear control theory OTT90 ; SCH99c  have been usefully applied to real world problems in various areas of physics, chemistry and biology BIE94 ; PIE96 ; HAL97 ; SUK97 ; PAR99 ; BEN00 ; LUE01 ; BEN02 ; BET03 , but no use has been made of this in the field of semiconductor self-oscillators.

Control methods can be either local or global BEC02 . Local methods require our ability to measure, and apply forcing directly to, the spatially resolved state variables of the system under study. However, in nanotechnology such variables, being e.g. electron densities in some quantum wells, are not easily accessible, and thus local methods cannot be applied. Unlike those, global methods require access only to some macroscopic variable(s) characterizing some integral output of the system. Such output can generally be reliably measured, and thus global methods seem to be the only option for the control of devices like superlattices. However, as we will show below, they are not straightforwardly applicable to nano-systems whose structure is spatially discrete. In this Letter we present a general approach to self-stabilization of irregular oscillations in semiconductor devices based upon essentially discrete quantum structures.

We consider a model for nonlinear electronic transport in semiconductor superlattices that yields complex and chaotic dynamic behavior under fixed time-independent external voltage in a regime where self-sustained dipole waves WAC02  are spontaneously generated at the emitter. Those dipole waves are associated with traveling field domains, and consist of electron accumulation and depletion fronts that in general travel at different velocities and may merge and annihilate. Such moving fronts are widespread in nonlinear, spatially extended systems, and similar chaotic front patterns occur in many other systems, e.g., spatially continuous models describing bulk impurity impact ionization breakdown in semiconductors CAN01  or globally coupled heterogeneous catalytic reactions GRA93a . Thus the time-delay autosynchronization method proposed in this work could be readily applied to stabilize similar space-time patterns in a variety of systems.

Our model of a superlattice is based on sequential tunneling of electrons WAC02 . In the framework of this model the quantum wells are assumed to be only weakly coupled, and electrons are localized at these wells. The tunneling rate to the next well is lower than the typical relaxation rate between the different energy levels within one well. The electrons within one well are then in quasi–equilibrium and transport through the barrier is incoherent. The resulting tunneling current density $`J_{mm+1}(F_m,n_m,n_{m+1})`$ from well $`m`$ to well $`m+1`$ depends only on the electric field $`F_m`$ between both wells and the electron densities $`n_m`$ and $`n_{m+1}`$ in the respective wells (in units of $`cm^2`$). A typical dependence of $`J_{mm+1}`$ on the electric field between two consecutive wells is $`N`$-shaped and exhibits a pronounced regime of negative differential conductivity.

The rate of variation of electron density in well $`m`$ is governed by the continuity equation

$$e\frac{\text{d}n_m}{\text{d}t}=J_{m1m}J_{mm+1}\text{for }m=1,\mathrm{}N$$
(1)
and Gauss’s law

$$ϵ_rϵ_0(F_mF_{m1})=e(n_mN_D)\text{for }m=1,\mathrm{}N,$$
(2)
where $`N`$ is the number of wells in the superlattice, $`ϵ_r`$ and $`ϵ_0`$ are the relative and absolute permittivities, $`e<0`$ is the electron charge, $`N_D`$ is the donor density, and $`F_0`$ and $`F_N`$ are the fields at the emitter and collector barrier, respectively. The total applied voltage $`U`$ between emitter and collector imposes a global constraint $`U=_{m=0}^NF_md`$, where $`d`$ is the superlattice period. This, together with (2), allows us to eliminate the field variables $`F_m(n_1,\mathrm{},n_N,U)`$ from the dynamic equations.

At the contacts Ohmic boundary conditions $`J_{01}=\sigma F_0`$, $`J_{NN+1}=\sigma F_Nn_N/N_D`$ are chosen, where $`\sigma `$ is the Ohmic contact conductivity, and the factor $`n_N/N_D`$ is introduced in order to avoid negative electron densities at the collector. The value of $`\sigma `$ essentially determines the oscillation mode AMA02a .

If the contact conductivity $`\sigma `$ is chosen appropriately, electron accumulation and depletion fronts are generated at the emitter. Those fronts form a traveling high field domain, with leading electron depletion front and trailing accumulation front. This leads to self-generated current oscillations. A fixed voltage $`U`$ imposes a constraint on the lengths of the high-field domains and thus on the front velocities. If $`N_a`$ accumulations fronts and $`N_d`$ depletion fronts are present, the respective front velocities $`v_a`$ and $`v_d`$ must obey $`v_d/v_a=N_a/N_d`$. If the accumulation and depletion fronts have different velocities, they may collide in pairs and annihilate. At certain combinations of contact conductivity $`\sigma `$ and voltage $`U`$, chaotic motion arises, when the annihilation of fronts of opposite polarity occurs at irregular positions within the superlattice AMA02a . The inset of Fig. 1 shows the plane of $`\sigma `$ and $`U`$, where regions with distinct regimes are marked by different shading. Black regions are those where chaotic behavior has been found. As a computationally convenient criterion for chaos we have used the rapid decay of the autocorrelation function estimated from $`n_{20}(t)`$. Chaotic regimes are found at low contact conductivity and low voltages, where dipole oscillations with leading accumulation and trailing depletion fronts occur, and at higher contact conductivity and higher voltage, where the role of accumulation and depletion fronts is interchanged. In Fig. 1, a one-parameter bifurcation diagram is given, obtained by plotting the time differences $`\mathrm{\Delta }t`$ between two consecutive maxima of the electron density in a specified well. The value of $`\sigma `$ is $`0.5\mathrm{\Omega }^1\text{m}^1`$, and $`U`$ is changed. Chaotic bands and periodic windows can be clearly seen.

The transition from periodic to chaotic oscillations is enlightened by considering the space-time plot for the evolution of the electron densities (Fig. 2a). At $`U=1.15\text{V}`$ chaotic front patterns with irregular sequences of annihilation of front pairs at varying positions within the superlattice occur. We have calculated the largest Lyapounov exponent as $`1.1\times 10^9s^1`$ which is a clear indication of chaos.

We shall now introduce a time delayed feedback loop to control the chaotic front motion and stabilize a periodic oscillation mode which is inherent in the chaotic attractor. In the extended time–delay autosynchronization scheme first suggested by Socolar et al SOC94 , multiple time delays are used to improve the control performance. Analytical insight into those schemes has been gained only recently BLE96 ; JUS97 ; NAK97 , and various ways of coupling of the control force, including local and global schemes, have been compared FRA99 ; BEC02 ; BAB02 ; JUS03 . Whereas local coupling schemes usually lead to efficient control in a large control domain, they are not easily implemented in real systems since local, spatially resolved measurements are necessary. Therefore, here we propose a much simpler global scheme. In our problem, as a global output signal that is coupled back in the feedback loop, it is natural to use the total current density $`J`$ defined as follows: $`J=\frac{1}{N+1}_{m=0}^NJ_{mm+1}`$ WAC02 . For the uncontrolled chaotic oscillations, $`J`$ is given in Fig. 2a by grey, showing irregular spikes at those times when two fronts annihilate. Note that the grey current time trace is modulated by fast small-amplitude oscillations (due to well-to-well hopping of depletion and accumulation fronts in our discrete model) which are not resolved in the plot. However, as the variable $`J`$ is fed back to the system for the purposes of control, these high–frequency oscillations render the control loop unstable. They need to be filtered out by using e.g. the following low–pass filter:

$`\overline{J}(t)`$ $`=`$ $`\alpha {\displaystyle _0^t}J(t^{})e^{\alpha (tt^{})}dt^{}\text{,}`$ (3)
with a cut-off frequency $`\alpha `$.

The multiple time delays of the resulting signal $`\overline{J}`$ (Fig. 2, black curve) are then used to modulate the voltage $`U`$ across the superlattice:

$`U`$ $`=`$ $`U_0+U_c(t)`$ (4)
$`U_c(t)`$ $`=`$ $`K\left(\overline{J}(t)\overline{J}(t\tau )\right)+RU_c(t\tau )`$
$`=`$ $`K{\displaystyle \underset{\nu =0}{\overset{\mathrm{}}{}}}R^\nu \left(\overline{J}(t\nu \tau )\overline{J}(t(\nu +1)\tau )\right)`$
where $`U_0`$ is a time–independent external bias, and $`U_c`$ is the control voltage. $`K`$ is the amplitude of the control force, $`\tau `$ is the delay time, and $`R`$ is a memory parameter. Such a global control scheme is easy to implement experimentally. It is non-invasive in the sense that the control force vanishes when the target state of period $`\tau `$ has been reached. This target state is an unstable periodic orbit of the uncontrolled system. The period $`\tau `$ can be determined by observing the resonance-like behavior of the mean control force versus $`\tau `$. The result of the control is shown in Fig. 2b. The front dynamics exhibits annihilation of front pairs at fixed positions within the superlattice, and stable periodic oscillations of the current are obtained. In Fig. 3 the control domain is depicted in the parameter plane of $`R`$ and $`K`$. A typical horn-like control domain similar to the ones known from other coupling schemes BEC02  is found. Typically, the left-hand control boundary corresponds to a period-doubling bifurcation, while the right-hand boundary is associated with a Hopf bifurcation. These findings show that our control scheme is robust.

To conclude, we have demonstrated that time-delay autosynchronization represents a convenient and simple scheme for the self-stabilization of high-frequency current oscillations due to moving domains in superlattices. This approach lacks the drawback of synchronization by an external ultrahigh-frequency forcing, since it requires nothing but delaying of the global electrical system output by the specified time lag. The proposed low-pass filtering of the output signal presents a solution of the problem one necessarily encounters when trying to control a nano-system with a crucially discrete quantum structure leading to superimposed fast well-to-well hopping oscillations in our case.

This work was supported by DFG in the framework of Sfb 555 and through grant no. JU261/3-1. NJ acknowledges partial support by EPSRC. We gratefully acknowledge discussion with A. Wacker.


format: text

text:
# On the nature of the ring-like structure in the outer Galactic disk

## 1 Introduction

Large-scale surveys like the Sloan Digital Sky Survey (SDSS) are starting to probe the limits of our understanding of the structure and dynamics of the Galaxy (Newberg et al. 2002). The most recent example is the discovery of a coherent structure at large Galactocentric distance and low galactic latitude spanning about 100 degrees on the sky (Yanny et al. 2003, hereafter Y03; Ibata et al. 2003, hereafter I03). The nature of this “ring” of low metallicity stars remains a puzzle. Y03 suggest that it could be the debris of a tidally disrupted satellite moving on roughly circular orbits, while I03 prefer the hypothesis that the ring results from perturbations in the disk such as the remnants of ancient warps. Neither scenario is clearly preferable. Could the debris of a disrupted galaxy remain a sharp feature at 20 kpc from the Galactic center while being stretched over one hundred degrees on the sky? What kind of disk instability would put metal-poor stars on such a distant ring around the Galaxy?

It is possible to shed some light on the viability of such scenarios by turning to cosmological simulations of galaxy formation. Over the past few years, these simulations have improved in resolution to the point of allowing a detailed analysis of morphological and structural features of galaxies and their relation to individual merger and accretion histories (Bekki & Chiba 2001; Springel & Hernquist 2003; Governato et al. 2002; Sommer-Larsen et al. 2002; Abadi et al. 2002a,b). Abadi et al. (2002b), in particular, point out that the bulk of old disk stars may originate from tidal debris of disrupted satellites whose orbital planes were roughly coincident with that of the disk. This result prompted us to analyze the evolution of satellite debris close to the galactic plane to explore whether it could explain the “ring-like” structure discovered by the SDSS.

The formation of ripples, arcs, and rings are common phenomena associated to mergers of galaxies. They arise from the folding and wrapping of tidally-stripped stars in the potential well of the host galaxy. In this Letter we show that the structure observed by Y03 and I03 could indeed be a tidal feature, and as such it may be considered a natural consequence of the hierarchical build-up of galaxies.

## 2 Tidal arcs in a cosmological simulation

Figure 1 illustrates the evolution of the debris of a tidally disrupted satellite in the cosmological simulation of the formation of a disk galaxy described by Abadi et al. (2002a,b). As discussed by these authors, this galaxy resembles more an early-type spiral than the Milky Way; still, the luminosity of the disk component ($`L_{\mathrm{disk}}2\times 10^{10}\mathrm{L}_{}`$) and the circular speed at $`R=10`$ kpc ($`250`$ km s-1 ) are comparable to those of our Galaxy. Hence, physical quantities quoted here may be taken as indicative of similar accretion events affecting the Galaxy without substantial rescaling. We emphasize, however, that the simulated galaxy is not a model of the Milky Way, and that our main goal is to understand qualitatively rather than quantitatively the nature of the structure observed by Y03 and I03. Indeed, the simulation precedes the announcement of the ring by Y03, and its results were thus in no way prejudiced by its discovery.

The satellite has an orbit roughly coplanar with that of the main disk of the galaxy, and its “face on” projection is shown in the left panels of Figure 1. The accretion of this satellite adds $`3.5\times 10^9M_{\mathrm{}}`$ (or about $`6\%`$) to the stellar mass of the preexisting galaxy. The first pericentric passage that results in significant stripping occurs at $`z0.67`$.

As is clear from Figure 1, sharp “ring-like” structures are readily seen at various times during the disruption of the satellite. In the top panels of Fig. 1, at $`z=0.67`$, most stars with $`4.6`$ kpc $`<R<7.2`$ kpc and $`100\mathrm{°}<\varphi <210\mathrm{°}`$ were stripped from the satellite during its first pericentric passage. The top panels (b) and (c) demonstrate that this incomplete ring (hereafter “tidal arc”, for short) is made out of particles with a range of energies that decrease systematically with increasing angle $`\varphi `$: particles at $`\varphi 175\mathrm{°}`$ are at the apocenter of their orbit, whilst those at $`\varphi 125\mathrm{°}`$ have yet to reach it, and those at $`\varphi 200\mathrm{°}`$ are already past it. Their coincidence on an “arc” of roughly constant radius is therefore a transitory event. Similar arcs are produced after every pericentric passage that strips a significant number of stars from the satellite; for example, a second one may be seen at $`z=0.60`$ between $`14.1`$ kpc $`<R<18.6`$ kpc, and $`45\mathrm{°}\varphi 120\mathrm{°}`$, as shown in panel (a’) of Fig. 1. These “tidal arcs” become progressively less well defined with time, as the orbits of particles in different energy levels precess around the main galaxy at different rates, weakening the radial coincidences shown in Fig. 1. After $`z=0.5`$, tidal arcs become very difficult to discern in the satellite debris in our simulation.

The tidal arc formation mechanism described above implies a number of well defined properties that may be verified observationally. Perhaps the most salient is illustrated in panels (b) and (b’) of Figure 1: a well defined radial velocity gradient is expected across the tidal arc. For the two tidal arcs discussed thus far, the (galactocentric) radial velocity varies by almost $`200`$ km s-1 across $`50`$ degrees in azimuthal angle. At any given location, the radial velocity distribution is fairly narrow, as shown in panels (d) and (d’) of Fig. 1. Typical velocity dispersions are of order $`20`$ to $`50`$ km s-1 .

Another notable feature is that tidal arcs span a limited azimuthal range. The actual range spanned depends both on the location of the observer as well as on the eccentricity of the satellite’s orbit at the time of stripping; more circular orbits lead typically to wider rings but somewhat less sharply defined in radius. Interestingly, the angular extent of the tidal arcs shown in Figure 1 ($`\mathrm{\Delta }\varphi 100\mathrm{°}`$) is comparable to that of the structure observed by Y03 and I03. We also note that particles on the tidal arcs are on rather eccentric orbits: the specific angular momentum of particles in the tidal arc seen at $`z=0.67`$ is almost a factor of $`2`$ lower than that of a circular orbit at the same radial location, and their $`R_{\mathrm{per}}/R_{\mathrm{apo}}`$ ratio averages $`0.18`$.

The fraction of the mass of the satellite contained within each arc varies according to the fraction of mass stripped from the satellite at each pericentric passage, which in turn depends on the orbit. Since it is possible to strip a substantial fraction of the satellite’s mass in one single passage (see, e.g., Hayashi et al. 2003), the rings might contain a significant fraction of the original mass of the satellite: in particular, the arcs identified at $`z=0.67`$ and at $`z=0.60`$ contain, respectively, $`8\%`$ and $`7\%`$ of the total mass of the satellite.

Finally, the orbit of the disrupting satellite is not exactly coplanar, leading to asymmetries above and below the plane of the disk: the rings identified in Figure 1 at $`z=0.67`$ and $`0.60`$ have, respectively, $`64\%`$ and $`60\%`$ of their mass above the plane of the disk, with a typical scaleheight of $`1`$ kpc. In summary, tidal arcs are structures characterized by being transitory, asymmetric in velocity and height, and of limited (but possibly substantial) angular extent.

## 3 Radial shells as relicts of minor mergers

One potential difficulty for the tidal-arc interpretation of ring-like features in the outer Galactic disk is that these arcs are transient configurations that phase-mix away within a few rotation periods. The actual timescale depends both on the initial properties as well as on the orbital period $`T`$ of the satellite. Typically the surface density will decrease in time as $`(T/t)^2`$ (Helmi & White 1999). The tidal arcs in Fig. 1 have $`T0.1`$ Gyr, hence they become a factor of ten fainter just one Gyr after the first pericentric passage. This explains why they cannot be recognized unequivocally in the cosmological simulation, particularly given our limited numerical resolution. We recall, however, that sharp features and ripples in physical space are expected to be long-lived under most circumstances, as shown, for example, by earlier work on the origin of shells around elliptical galaxies (e.g. Quinn 1984; Dupraz & Combes 1986; Hernquist & Quinn 1988, 1989; Heisler & White 1990).

To test the hypothesis that other ring-like structures might be long-lived, we have followed the disruption of a $`10^6`$-particle satellite (modeled after the one shown in Fig. 1) in a fixed galactic potential. We use a quadrupole expansion code (Zaritsky & White 1988), and represent the satellite with a King profile of concentration $`\mathrm{log}_{10}(r_t/r_c)=1.34`$, total mass $`3.44\times 10^9M_{\mathrm{}}`$, total one-dimensional velocity dispersion $`59.4`$ km s-1, and core radius $`r_c=0.82`$ kpc. The main galaxy is chosen to match that of the cosmological simulation, and is approximated by the superposition of a Miyamoto-Nagai disk (Miyamoto & Nagai 1975), a Hernquist bulge (Hernquist 1990) and an NFW halo (Navarro, Frenk & White 1997). The initial position and velocity of the center of mass of the satellite are those of the satellite in the cosmological simulation at $`z=0.74`$, before any significant disruption takes place.

Since the goal of this (re)simulation is to study the long-term survival of ring-like features, we focus our analysis on the structure of the debris at $`z=0`$ (7 Gyr after infall). Panel (a) in Fig. 2 shows a “face-on” view of the satellite debris at this time: several sharply-defined ring-like features are easily identified in this panel, especially in the outer regions. To illustrate the properties of these features, we choose the annulus defined by $`34.7`$ kpc $`<R<37.1`$ kpc and $`70\mathrm{°}<\varphi <150\mathrm{°}`$. The radial enhancement is dominated by particles at the apocenter of their orbits (i.e., those with $`v_R0`$ in panel (c) of Fig. 2). This density structure is a caustic surface (2-D singularity or fold) in the mapping of the phase-space density of a cold/small system into the 3-D physical space (Nulsen 1989; Hernquist & Quinn 1988). This structure is analogous the “shell” structures found around some elliptical and a few spiral galaxies (Malin & Carter 1983; Schweizer & Seitzer 1988; Wilkinson et al. 2000).

To assess the visibility of shells when a spread of magnitudes is present, we have convolved the distance distribution of particles with $`140\mathrm{°}\varphi 150\mathrm{°}`$ and $`R25`$ kpc with the magnitude distribution near the turnoff of the globular cluster M92, i.e. for $`0.16(gr)0.3`$ (Stetson & Harris 1988)111The $`B`$ and $`V`$ magnitudes have been transformed to the $`g`$ and $`r`$ passbands following Cohen (1985).. Here we assumed the observer is located at $`R25`$ kpc, which implies that its distance to the shell is the same as that of the structure observed by SDSS. The apparent magnitude distribution $`N(m)`$ of these particles, plotted in panel (2e), is similar to that shown in Newberg et al. (2002).

Panel (b) of Figure 2 shows that there are no gradients in the radial velocity across the “shell”. The radial velocity distribution at any given location on the shell has several components, as can be readily seen from panels (2c) and (2d): it is dominated by a central core with $`|v_R|<100`$ km s-1 and $`\sigma _R40`$ km s-1 comprising $`70\%`$ of the stars, together with a number of small peaks corresponding to particles of higher energy on their way to or from their (more distant) apocenters. The presence of these peaks at high $`|v_R|`$ increases the radial velocity dispersion of the ring to $`120`$ km s-1 .

Shells are relaxed features that are expected to be symmetric above/below the plane of the disk: roughly $`50\%`$ of the particles in the shell identified in Fig. 2 lie above the disk plane. Shells are long-lived, but it is unlikely that they may carry a significant fraction of the original mass of the satellite: the amount of mass they contain decreases rapidly as $`(T/t)^3`$. The radial shell shown in Fig. 2 contains only $`0.7\%`$ of the satellite’s mass, an order of magnitude smaller than the fraction of mass carried by the tidal arcs shown in Fig. 1. Other long-lived shell structures that are easily identified in the same figure carry a similarly low fraction ($`<1\%`$) of the original mass of the satellite.

## 4 Discussion

The main purpose of the simulations presented here is not to match quantitatively the properties of the structure discovered by SDSS, but rather to show that qualitatively similar features arise naturally as a result of the accretion of satellites on orbits roughly coplanar with the Galactic disk. Tidal arcs (§ 2) are prevalent during the disruption process (which lasts a few orbital times, or $`1`$ Gyr) whereas shells (§ 3) are long-lasting features that may be identified long after the debris has reached equilibrium in the main galactic potential.

A cursory comparison with the properties of the structure observed by Y03 and I03 shows that both interpretations are consistent with observations. Indeed, the extent in Galactic longitude is comparable to both arcs and shells shown in Figs. 1 and 2. The systematic variation in the mean radial velocity with longitude claimed by Y03, from $`22`$ km s-1 at $`(l,b)=(182\mathrm{°},28\mathrm{°})`$, to $`103`$ km s-1 at $`(l,b)=(225\mathrm{°},27\mathrm{°})`$, could be, but is not necessarily, a signature of a “tidal arc”. This “gradient” can also be entirely attributed to perspective effects resulting from the projection of the solar motion and the tangential velocity of stars in the ring along these two different lines of sight. This degeneracy can be broken with radial velocity measurements along more widely separated lines of sight, and a gradient, if detected, would then favor the tidal arc interpretation. On the other hand, larger spectroscopic samples would allow to test whether the radial velocity distribution is smooth or if it consists of multiple peaks as in the shell scenario. Measurements of the extent of the structure in Galactocentric longitude, and above/below the plane would provide other possible ways of distinguishing amongst the two interpretations.

Discriminating between these two scenarios is important in order to estimate robustly the time of accretion as well as the original mass of the satellite. The tidal arc interpretation would suggest a very recent accretion event of a fairly small satellite; if the arc comprises $`10\%`$ of the original mass of the satellite, the mass estimated for the ring by Y03 would imply a progenitor mass of roughly $`2.5\times 10^8\mathrm{M}_{}`$. Thus, the Sagittarius dwarf galaxy (Ibata, Gilmore & Irwin 1994) might not be the only satellite being accreted today by the Milky Way. The shell scenario, on the other hand, would imply the past accretion of a satellite roughly ten times more massive ($`2.5\times 10^9\mathrm{M}_{}`$), and the presence of considerable debris in the inner (thick) disk. Such stars are expected to have quite distinct kinematics and to be on rather eccentric orbits, perhaps consistent with the peculiar kinematics of the thick disk in the solar neighborhood reported by Gilmore, Wyse & Norris (2002).

We end by noting that tidal ripples of the kind discussed here are not a new idea, and were discussed by Hernquist & Quinn (1988) and others well before our study. The cosmological simulation analyzed here suggests that satellites will often have orbits approximately coplanar with the disk of the primary, and thus a significant fraction of their debris should be concentrated in the Galactic disks (Abadi et al. 2002b). Identifying and cataloging such structures will be possible with large surveys like RAVE (Steinmetz 2002) and GAIA (Perryman et al. 2001), and should help us unravel the puzzle of the assembly of the Galaxy.

This work has been partially supported by grants from the U.S. National Aeronautics and Space Administration (NAG 5-10827) and the David and Lucile Packard Foundation.


text:
# Inclusive Jet Production, Parton Distributions, and the Search for New Physics

## I Introduction: Jet cross sections at the Tevatron

Jet production at the Tevatron probes the highest momentum transfer region currently accessible. As this region is potentially sensitive to a wide variety of new physics, there was great interest when the inclusive jet cross section measured by the CDF collaboration in Run 1b exhibited an excess in the highest $`E_T`$ range, when compared to NLO predictions using then-current parton distribution functions cdf1A . In an attempt to determine if conventional physics could explain the deviation, the CTEQ PDF fitting group carried out a global analysis using information from deep-inelastic scattering and Drell-Yan data, as well as from jet data from the Tevatron, but giving a higher statistical emphasis to the high $`E_T`$ jet data from CDF cteq4hj . NLO predictions using the resulting fit (CTEQ4HJ) reduced the size of the excess observed by CDF. The jet data from both CDF and DØ were also used in a more conventional fit (CTEQ4M) where no special emphasis was given to the high $`E_T`$ data cteq4 .

The dominant subprocess that contributes to jet production at high $`E_T`$ is quark-(anti)quark scattering. However, the quark distributions in the corresponding $`x`$ range are very well constrained by the precise DIS and DY data used in the global fits. Only the gluon distribution has the flexibility to change significantly in the high $`x`$ region, and indeed the gluon distribution increases by roughly a factor of 2 for $`x`$ values of about 0.5 in the CTEQ4HJ fit. The gluon-quark scattering subprocesses increase from approximately 20% of the total jet cross section at high $`E_T`$ using CTEQ4M to 40% using CTEQ4HJ cteq4hj .

The next group of PDF’s from CTEQ (CTEQ5) cteq5  also contained two sets: CTEQ5M, the standard lowest $`\chi ^2`$ solution, and CTEQ5HJ, defined by a similar statistical enhancement applied to the high $`E_T`$ data from CDF cdf1b  and DØ d01b .

More recently, DØ has measured the inclusive jet cross section as a function of rapidity $`y`$ over the range $`0|y|3`$ d0rapidity . This data set, comprised of 90 data points, has a greater statistical power in the global fits than the CDF and DØ jet cross section measurements in the central region alone. The CTEQ6M fit cteq6  utilizes these DØ jet cross section measurements, along with the CDF measurements in the central rapidity region, as well as the most recent DIS data from HERA and existing fixed target DIS and DY data. An important step was taken to quantify the uncertainties of the parton distributions and their physical predictions with the inclusion of 40 eigenvector basis PDF sets, along with the best fit, to characterize the uncertainties in the neighborhood of the global minimum, using the Hessian method.

This paper is devoted to a detailed and more focused study of the uncertainties of the inclusive jet cross section due to global constraints, both at current energies, and more importantly for further QCD studies and the search for new physics at future collider programs. During the course of this work, we have produced an improved version of the CTEQ6 PDF’s. The new minimum set, called CTEQ6.1M henceforth, provides a global fit that is almost equivalent in every respect to the published CTEQ6M cteq6 , although some parton distributions (e.g., the gluon) may deviate from CTEQ6M in some kinematic ranges by amounts that are well within the specified uncertainties. The more significant improvements are associated with some of the 40 eigenvector sets, which are made more symmetrical and reliable in CTEQ6.1M111The CTEQ6.1M parton distributions will be available in the CTEQ table format at www.phys.psu.edu/ cteq/ and in the LHAPDF format at vircol.fnal.gov.

The comparison between the DØ data for the 5 different rapidity intervals, and the predictions using the CTEQ6.1M PDF’s, is shown in Fig. 1. The theory predictions for CTEQ6M, CTEQ5M and CTEQ5HJ parton distributions are also shown. The greater statistical power that results from including the new DØ jet data, and the fact that all of the rapidity intervals prefer a larger gluon at high $`x`$, results in the increases in the cross sections found using NLO predictions based on CTEQ6M or CTEQ6.1M rather than CTEQ5M.

This point deserves further emphasis. It is crucial to measure the jet cross section over a wide rapidity range: any new physics should contribute mostly to the central rapidity region, while a PDF explanation will contribute to all rapidity ranges. This will be addressed in more detail below.

As discussed in Ref. cteq6 , it is important to have a parameterization for the PDF’s that is flexible enough to allow an increase in the high $`x`$ gluon, without modifying too much the shape of the gluon distribution at lower values of $`x`$. A comparison of the CTEQ5M, CTEQ5HJ and CTEQ6M gluon distributions is shown in Fig. 2. As can be observed, the CTEQ6M gluon at high $`x`$ lies between the gluons in CTEQ5M and in CTEQ5HJ, just as the CTEQ6M predictions for the inclusive jet cross sections at high $`E_T`$ lie between the predictions using the two other PDF’s. The CTEQ6.1M gluon distribution is shown by dashed curves in Fig. 2. At the highest $`x`$ values the function resembles CTEQ5HJ, and in the intermediate $`x`$ range it resembles CTEQ6M.

The subprocess contributions to jet production are shown in the next figures, for the central rapidity bin (Fig. 3) and the forward rapidity bin (Fig. 4), for the parton distributions CTEQ5M and CTEQ6M.

The dominant difference between the CTEQ5M and CTEQ6M sets is the increased gluon distribution in CTEQ6M at large values of $`x`$; the quark distributions are nearly unchanged. This shows up as a decrease in the fractional contributions of the quark-quark subprocesses and as increases in the fractional contributions of the quark-gluon and gluon-gluon subprocesses. It is interesting to note that the $`gq`$ scattering subprocess is even more important at high $`E_T`$ for the highest rapidity bin than for the central rapidity bin.

The cross section for production of high-$`E_T`$ jets at the Tevatron depends on the parton distributions at large $`x`$. The relevant range of $`x`$ may be estimated by considering the leading-order $`22`$ kinematics. Let $`x_1`$ and $`x_2`$ be the momentum fractions of the two incoming partons. The outgoing parton 4-momenta are

$`p_1^\mu `$ $`=`$ $`(p_T\mathrm{cosh}y_1,\stackrel{}{p}_T,p_T\mathrm{sinh}y_1),`$ (1)
$`p_2^\mu `$ $`=`$ $`(p_T\mathrm{cosh}y_2,\stackrel{}{p}_T,p_T\mathrm{sinh}y_2);`$ (2)
the outgoing partons have rapidities $`y_1`$ and $`y_2`$, and transverse momenta $`\stackrel{}{p}_T`$ and $`\stackrel{}{p}_T`$, respectively. If the first parton is identified with the observed jet then the jet rapidity is $`y_j=y_1`$ and its transverse momentum is $`\stackrel{}{p}_T`$. The $`22`$ kinematics determines the relation between the momentum fractions and outgoing parton variables,

$`x_1`$ $`=`$ $`{\displaystyle \frac{p_T}{\sqrt{s}}}\left(e^{y_j}+e^{y_2}\right),`$ (3)
$`x_2`$ $`=`$ $`{\displaystyle \frac{p_T}{\sqrt{s}}}\left(e^{y_j}+e^{y_2}\right).`$ (4)
These parametric equations yield a curve in the ($`x_1,x_2`$) plane, parametrized by the rapidity $`y_2`$ of the second parton, for the specified jet variables.

Figure 5 shows the $`(x_1,x_2)`$ values for central rapidity ($`y_j=0`$) on the left and forward rapidity ($`y_j=2`$) on the right. In the central case, the solid curve is the locus of $`(x_1,x_2)`$ points for $`p_T=400`$ GeV and the dashed curve is for $`p_T=200`$ GeV. Both $`x_1`$ and $`x_2`$ must be $`0.25`$ to produce the high-$`p_T`$ jet. In the forward case the solid curve corresponds to $`p_T=200`$ GeV and the dashed curve to $`p_T=100`$ GeV. Here one parton must have large momentum fraction to produce the jet at large $`p_T`$ and forward rapidity. The second parton can have a low value of $`x`$, leading to the second jet being close in rapidity to the first jet, or a larger value of $`x`$, leading to configurations where the two jets are on opposite sides of the detector. Consider the contribution of the quark-gluon scattering subprocesses. For the configuration where one $`x`$ is small and the other is large, the smaller $`x`$ value most likely corresponds to a gluon and the larger to a quark. In this case the difference between the CTEQ5M and CTEQ6M distributions is small. On the other hand, for the configuration where both $`x`$ values are large the increased gluon in CTEQ6M will enhance the cross section. The conclusion is that the increase in the cross section when going from CTEQ5M to CTEQ6M occurs in configurations where the forward jet is balanced in $`p_T`$ by a jet with roughly the opposite rapidity.

Next, consider the question of where in the jet phase space one is most likely to observe signs of compositeness or other new physics. The angular distributions of the dominant QCD subprocesses are sharply peaked at small scattering angles, whereas the models usually used to estimate compositeness have much flatter distributions. This means that the signal/QCD ratio will be largest for 90 scattering in the parton-parton rest frame. In the overall hadron-hadron frame, this corresponds both jets having the same rapidity. Using the lowest order kinematics described previously, the squared dijet mass is given by

$$M_{jj}^2=2p_T^2[1+\mathrm{cosh}(y_1y_2)].$$

Setting $`y_1=y_2`$ yields $`M_{jj}^2=4p_T^2.`$ This shows that maximizing the dijet mass means maximizing $`p_T`$. Since the $`p_T`$ reach is maximized in the central region, this is the region where the highest dijet masses can be reached with the lowest background for new physics signals (at least those with a flat angular distribution). The model discussed for compositeness later in this paper gives an explicit example of this behavior. On the other hand, as shown by the kinematic examples in Figs. 5, modifications of the gluon distribution at high values of $`x`$ will affect the high-$`p_T`$ jet cross section at all rapidities.

## II Uncertainty analysis of the cross section for inclusive jet production at the Tevatron

As discussed in the previous section, the cross section for inclusive jet production has been measured by the CDF and DØ collaborations at the Tevatron. These results have been used in the latest CTEQ global analysis of PDF’s, leading to the CTEQ6 (and now 6.1M) parton distributions. In this section we analyze the uncertainty of the theoretical QCD cross section due to PDF uncertainty.

### II.1 The Hessian method of uncertainty analysis

The Hessian method of uncertainty analysis has been described in detail previously PSTHessian . For completeness, we summarize the method briefly. The parton distributions are constructed by a method of chi-square minimization with fitting of systematic errors. A chi-square function $`\chi ^{\mathrm{\hspace{0.17em}2}}`$ is defined by

$`\chi ^{\mathrm{\hspace{0.17em}2}}`$ $`=`$ $`{\displaystyle \underset{e}{}}\chi _e^2(a,r)`$ (5)
$`\chi _e^2(a,r)`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{\left[D_i_kr_k\beta _{ki}T_i(a)\right]^2}{\alpha _i^2}}+{\displaystyle \underset{k}{}}r_k^2`$ (6)
where $`e`$ labels an experimental data set and $`i`$ labels a data point in the data set. $`D_i`$ is the data value, $`\alpha _i`$ is the uncorrelated error, and $`\beta _{ki}`$ is the $`k`$th correlated systematic error; these numbers are published by the experimental collaboration. Then $`T_i(a)`$ is the theoretical value, a function of a set of $`n`$ PDF parameters, $`\{a_1,\mathrm{},a_n\}`$. Also, $`\{r_k\}`$ is a set of Gaussian random variables and $`r_k\beta _{ki}`$ is a (correlated) shift applied to $`D_i`$ to represent the $`k`$th systematic error. We minimize the function $`\chi ^{\mathrm{\hspace{0.17em}2}}(a,r)`$ with respect to both the PDF parameters $`\{a\}`$ and the systematic shift variables $`\{r_k\}`$. The result yields both the standard PDF model with parameters $`\{a_0\}`$, and the optimal shifts $`\{\widehat{r}_k\}`$ to bring theory and data into agreement. This minimum of $`\chi ^{\mathrm{\hspace{0.17em}2}}`$ represents the best fit to the data cteq6 .

The Hessian matrix is

$$H_{ij}=\frac{1}{2}\frac{^2\widehat{\chi }^2}{a_ia_j}\text{where}\widehat{\chi }^2(a)\chi ^{\mathrm{\hspace{0.17em}2}}(a,\widehat{r}(a)).$$
(7)
This matrix determines the behavior of $`\widehat{\chi }^2(a)`$ in the neighborhood of the minimum. We have developed an iterative method for calculating the Hessian accurately, leading to an improvement in the $`\chi ^2`$ minimization PSTIterative . The point $`\{a_0\}`$ in the $`n`$-dimensional parameter space—where $`\widehat{\chi }^2(a)`$ is minimum—is the best fit to the global data set. However, points in some small neighborhood of $`\{a_0\}`$ are also acceptable fits. The variation within the acceptable neighborhood, of a PDF or a physical prediction, is our measure of the uncertainty of the quantity.

In the Hessian method of uncertainty analysis, we first define $`n`$ special directions in the parameter space, namely the $`n`$ eigenvectors of the Hessian. Then for each eigenvector we have two displacements from $`\{a_0\}`$ (in the $`+`$ and $``$ directions along the vector) denoted $`\{a_i^+\}`$ and $`\{a_i^{}\}`$ for the $`i`$th eigenvector. At these points, $`\widehat{\chi }^2=\widehat{\chi }_0^2+T^2`$ where $`\widehat{\chi }_0^2=\widehat{\chi }^2(a_0)`$ $`=`$ the minimum, and $`T`$ is a parameter called the tolerance. We consider any PDF set with $`\widehat{\chi }^2\widehat{\chi }_0^2<T^2`$ to be an acceptable fit to the global data set. In particular, the $`2n`$ PDF sets $`\{a_i^\pm \}`$, which are called the eigenvector basis sets, span the parameter space in the neighborhood of the minimum.

The appropriate choice of tolerance $`T`$ cannot be decided without a further, more detailed, analysis of the quality of the global fits. After studying a number of examples PSTHessian ; PSTLagrange  we conclude that a rather large tolerance, $`T10`$, represents a realistic estimate of the PDF uncertainty.222Different choices of the tolerance have been made by other groups, e.g., $`T\sqrt{40}`$ by MRST mrst2002 .

Any quantity $`X`$ that depends on PDF’s has a predicted value $`X_0=X(a_0)`$ and a range of uncertainty $`\delta X`$. A simple measure of $`\delta X`$ is the spread of the values of $`X(a_i^\pm )`$ for the $`2n`$ eigenvector basis sets. However, a more complete uncertainty range is between the minimum and maximum values of $`X`$ for all points with $`\widehat{\chi }^2\widehat{\chi }_0^2<T^2`$. It can be shown that in a linear approximation these bounds are $`X_0\pm \delta X`$ where

$$\left(\delta X\right)^2=T^2\underset{i,j}{}\left(H^1\right)_{ij}\frac{X}{a_i}\frac{X}{a_j};$$
(8)
or, in terms of the eigenvector basis sets,

$$\left(\delta X\right)^2=\frac{1}{4}\underset{k=1}{\overset{n}{}}\left[X(a_i^+)X(a_i^{})\right]^2.$$
(9)
We refer to (9) as the Master Equation for the uncertainty of $`X`$ in the Hessian method PSTHessian . Equation (9) is based on a linear approximation: $`\widehat{\chi }^2(a)`$ is assumed to be a quadratic function of the parameters $`\{a\}`$, and $`X(a)`$ is assumed to be linear. We find that these approximations are not strictly valid, so instead we calculate asymmetric bounds by

$$\delta X_\pm =\underset{k=1}{\overset{n}{}}\left[X(a_i^\pm )X(a_0)\right]^2;$$
(10)
the range of uncertainty of $`X`$ is $`(X_0\delta X_{},X_0+\delta X_+)`$. (In (10) $`\{a_i^+\}`$ and $`\{a_i^{}\}`$ are the displaced points where $`X>X_0`$ and $`X<X_0`$, respectively.)

The uncertainties on the gluon and $`u`$-quark distributions are shown in Fig. 6, for $`Q^2=10`$ GeV2. The $`u`$-quark distribution is tightly constrained for $`x0.8`$, whereas the uncertainty on the gluon distribution grows to a factor of 2 or larger for $`x`$ values greater than $`0.4`$. For comparison, the CTEQ5M, CTEQ5HJ, MRST2001 mrst2001  and MRST2002 mrst2002  PDF’s are shown, compared to CTEQ6.1M, in Fig. 7, for $`Q^2=10`$ GeV2.

### II.2 The Tevatron jet cross section

To illustrate our method, and the importance of the systematic errors, Fig. 8 shows the fractional difference between the data obtained by CDF in Run 1b, and the theory with CTEQ6.1M partons. The error bars are the statistical errors only. The left panel is $`(D_iT_i)/T_i`$. Note that there are some systematic differences between the CDF jet data and the theory, both in terms of normalization and in shape.333It is important to note all Run 1 CDF cross sections are intrinsically larger than comparable DØ cross sections by approximately 6% due to different assumptions regarding the total inelastic cross section at 1.8 TeV cdf\_d0\_sigma . The right panel shows $`(\overline{D}_iT_i)/T_i`$ where $`\overline{D}_i`$ is the shifted data,

$$\overline{D}_i=D_i\underset{k=1}{\overset{n}{}}\widehat{r}_k\beta _{ki};$$
(11)
i.e., the optimal systematic shifts determined by the global fit have been subtracted from the data values. There is no systematic difference between $`\overline{D}_i`$ and $`T_i`$, as expected. The sizes of the systematic shifts are, as expected, comparable to the standard deviations published by the CDF collaboration; that is, $`\widehat{r}_k`$ is of order $`1`$ for each of the systematic errors.

Figure 9 shows the cross section for inclusive jet production at the Tevatron, integrated over the rapidity range $`0.1<|\eta |<0.7`$, as a function of jet transverse energy $`E_T`$, for CTEQ6.1M partons and for the 40 eigenvector sets. (The rapidity interval is that for which the cross section has been published by the CDF collaboration.) The spread of the 41 curves in Fig. 9 is a simple estimate of the PDF uncertainty. The CDF data points are superimposed on the calculation.

Figure 10 shows the fractional differences $`\left[T(a_i^\pm )T(a_0)\right]/T(a_0)`$ between the 40 sets and the standard prediction (CTEQ6.1M). The points are the CDF data compared to CTEQ6.1M; the error bars are statistical errors and the data points are plotted without the systematic shifts. Also shown is the prediction using CTEQ5HJ partons. The range of uncertainty encompasses both the CDF data and the CTEQ5HJ predictions.

Figure 11 shows the 40 eigenvector basis sets separately. For each eigenvector the two curves are the positive and negative displacements. The Hessian method assumes that the variations are approximately linear. If the linear approximation were strictly valid the curves would be mirror images about zero. Figure 11 implies that the linear approximation is reasonably accurate for most eigenvectors, especially those corresponding to the ten largest eigenvalues (1 – 10).

Figure 12 shows the overall PDF uncertainty calculated from the “Master Equation,” asymmetrically for positive and negative differences. The ordinate is the fractional difference between the extreme cross section and the central (CTEQ6.1M) prediction. The outer curves are the upper and lower bounds of the uncertainty band. The CDF data is superimposed.

The overall uncertainty in the calculation of $`(d\sigma /dE_T)_{\mathrm{CDF}}`$ is dominated by one of the eigenvector directions. The dominant eigenvector is “eigenvector 15,” i.e., that for which the eigenvalue of the Hessian is 15th in order of magnitude. The $`+`$ and $``$ displacements in parameter space along this direction have a much larger effect on the value of the cross section than any of the other eigenvector directions. In other words, $`(d\sigma /dE_T)_{\mathrm{CDF}}`$ is most sensitive to the parameter variations along this eigenvector. Figure 13 shows the gluon distribution for CTEQ6.1M and for the eigenvector basis sets in the $`+`$ and $``$ directions along eigenvector 15. Eigenvector 15, like the 19 others, corresponds to some linear combination of all of the 20 fitting parameters. The exceptional feature of eigenvector 15 is that it represents the largest change in the high-$`x`$ gluon behavior.

In the new CTEQ6.1M analysis, the eigenvector basis sets have been carefully constructed to have $`\widehat{\chi }^2=\widehat{\chi }_0^2+100`$.444The earlier CTEQ6 analysis cteq6  made an approximation of the NLO cross sections—holding the $`K`$-factor tables for DIS and DY processes constant as the PDF parameters varied in the neighborhood of the minimum—such that $`\mathrm{\Delta }\widehat{\chi }^2`$ was not strictly constant at $`T^2`$ for the Hessian PDF sets. The CTEQ6.1M analysis is preferred although the changes from CTEQ6 are small. But $`\widehat{\chi }^2`$ is the global chi-square function. The criterion $`\mathrm{\Delta }\widehat{\chi }^2=100`$ is only a reasonable rule of thumb, and a more thorough analysis—inspecting the fit for each data set—is needed to justify the extreme sets $`a_{15}^+`$ and $`a_{15}^{}`$. We must ask whether these variations could be ruled out by any specific experiment.

The theoretical cross section most affected by displacements along eigenvector 15 is the inclusive jet cross section. Table 1 shows the $`\widehat{\chi }^2`$ values for the CDF and DØ data for three PDF sets ($`a_0`$, $`a_{15}^+`$ and $`a_{15}^{}`$.) Because of the significant systematic uncertainties in the experimental data, represented by the error correlation matrices published by the collaborations, we might conclude that the extreme sets $`a_{15}^+`$ and $`a_{15}^{}`$ are still acceptable fits to the data within the systematic errors.

As implied by the table above, much of the increase in $`\chi ^2`$ (50 out of 100) for eigenvector +15 comes from the DØ jet cross section. The Run 1b DØ cross section can tolerate some increase compared to the CTEQ6.1M predictions, but the increase corresponding to eigenvector +15 is in fact too extreme. It is not surprising that such a large portion of the increase in $`\chi ^2`$ comes from the DØ jet data: as shown in Fig. 13, this eigenvector probes primarily the high-$`x`$ gluon distribution and the other data sets (with the exception of the CDF jet data) are fairly insensitive to this direction. The opposite direction, eigenvector $`15`$, gives the lowest $`\chi ^2`$ for the DØ jet data. (Of course the $`\chi ^2`$’s of other experiments are larger so that the net change is $`100`$.) The cross section in this case is less than that for CTEQ6.1M, by as much as $`10`$$`20`$%. The fact that $`\widehat{\chi }^2/N`$ is much less than 1 for the DØ data suggests that the size and correlation of systematic errors are not accurately understood and thus CTEQ6.1M gives a satisfactory fit although not the best fit.

Figures for the DØ jet data analogous to Figs. 9 and 10 are shown in Figs. 14 and 15. Figure 14 shows the inclusive jet cross section as a function of $`E_T`$ for the five rapidity bins of the DØ measurements. Calculations for all 40 eigenvector basis sets are superimposed. The points are the DØ data and the error bars equal the single-point statistical and systematic errors added in quadrature.555The systematic errors for the DØ data are supplied by the collaboration only in the form of correlation matrices. Unlike the CDF data, for which the standard deviations of individual systematic errors have been published, we are unable to separate the optimal systematic shifts for subtraction from the DØ data. Figure 15 shows the inclusive jet cross section as a function of $`E_T`$ for the five rapidity bins, plotted as fractional differences compared to the central (CTEQ6.1M) model.

## III Reliability of next-to-leading-order predictions

Thus far, the discussion of uncertainties has dealt primarily with the issue of the propagation of experimental errors on the observables and their effect on the parameters of the parton distributions. There are, of course, issues related to theoretical uncertainties both in the global fit itself and in the theoretical calculations for jet production. However, it is far more difficult to quantify these. We have previously described studies of some theoretical uncertainties cteq6 . For example, we have chosen cuts on $`Q`$ and $`W`$ for the DIS data sets used for CTEQ6 and earlier analyses in order to reduce the theoretical uncertainties due to effects from higher twist and related terms. Also, we have varied the functional form of the parametrization to ensure that sufficient flexibility to cover the wide kinematic range of $`x`$ and $`Q`$ is available. In this section we concentrate on the theoretical uncertainties related to NLO calculations for jet production. The response of the predictions to variations in the renormalization and factorizations scales yields some information on the theoretical uncertainties. In this section the scale dependence of the jet cross sections and the effects on the fits of varying the scales will be discussed. Some estimates of the possible effects of two-loop corrections will also be presented.

### III.1 Scale dependence of the jet cross section

Leading-logarithm calculations of the jet cross section at high-$`E_T`$ generally decrease monotonically as the factorization and renormalization scales are increased. The scale variation of next-to-leading-logarithm calculations is generally reduced as compared to that of the leading-log calculations. A brief review of how this comes about is presented in the Appendix.

The increased rapidity coverage available with the recent DØ data d01b  has provided additional constraints on the gluon distribution. However, if the theoretical uncertainties become large at high rapidity, then these constraints may not be as useful as first believed. As one approaches the edge of phase space, it might be anticipated that there will be large logarithms due to the constraints on multiple gluon emission. The effects of these large logs will be investigated in Sec. (III.3). Their effects may show up already in the next-to-leading-log calculations used in the global fits, so it is important to look for any anomalous behavior in the predictions.

The first step is to examine the scale dependence of the jet cross section predictions in the highest rapidity bin of the DØ data. For ease of presentation, the factorization and renormalization scales have been chosen to be the same. In the jet cross section calculations used for the CTEQ PDF fits, and in this paper, a scale $`\mu _R=\mu _F=E_T/2`$ has been used. Fig. 16 shows the comparison of the scale dependences of the leading-log and next-to-leading-log calculations at $`E_T=70`$ GeV/$`c`$. The monotonically decreasing behavior of the leading-log results is evident, as is the roughly quadratic dependence of the next-to-leading-log curve. Both behaviors are consistent with the expectations discussed in the Appendix. Similar results are shown in Fig. 17 for $`E_T=170`$ GeV/$`c`$. Both figures show that the scale dependence of the next-to-leading-log calculation is greatly reduced with respect to that of the leading-log calculation, provided that one avoids scale choices much smaller that $`\mu =E_T`$.

It is interesting to note that the $`K`$-factor (NLO/LO) for $`\mu =E_T/2`$ is approximately unity for $`E_T=70`$ GeV/$`c`$, while it is significantly less for $`E_T=170`$ GeV/$`c`$. The reasons for this are discussed in more detail in Appendix A.

### III.2 Fits with different scales

The next point to be addressed is the effect on the fits of varying the choice of scale. To investigate this point, the CTEQ6M fits666This procedure has not been repeated with the CTEQ6.1M PDF’s. However, since the changes in the region of the DØ and CDF data were minimal, the scale dependence in the 6M and 6.1M sets should be essentially the same. were repeated with scale choices of $`\mu =E_T`$ and $`2E_T`$. As compared to the original fit with $`\mu =E_T/2`$ the latter two yielded increases in chi-square of 16 and 70, respectively. Thus, the choice of $`E_T/2`$ used in the CTEQ6M fits yields the best chi-square of the three choices, but this does not constitute a fit for the optimal scale. Note, however, that the rapid increase in the scale dependence for choices below about $`E_T/2`$ shown in Figs. 16 and 17 precludes fits with significantly smaller scales while the larger scale choices show an increase in chi-square.

The gluon distributions from these fits are shown in Fig. 18 relative to the standard CTEQ6M distribution corresponding to $`\mu =E_T/2`$. The effects at values of $`x`$ below about $`0.4`$ are minimal, but the effects at larger $`x`$ values are significant. In particular, the reduction in the cross section due to the larger scale choice of $`2E_T`$ is compensated by an increase of the gluon distribution. However, it should be noted that the chi-square increased by 70 units for this fit. This is still less than the chi-square increase of 100 that determines the 40 error PDF’s.

A comparison of Fig. 18 to Fig. 6 indicates that the uncertainty on the gluon distribution due to the scale choice for jet production in the global fits is everywhere less than the uncertainty from the treatment of the experimental errors.

It is important to bear in mind that the gluon distribution itself is not an observable quantity. Rather, the gluon and quark distributions contribute via convolutions with appropriate hard scattering cross sections. Accordingly, the impression provided by Fig. 18 may be somewhat misleading as much of the gluon variation is offset by changes in the hard scattering cross section with which it is convoluted, with the end result being that the variation in the physical observables is much less.

The jet cross section predictions using scales of $`E_T/2`$, $`E_T`$ and $`2E_T`$ are shown in Figs. 19 and 20 for the DØ rapidity regions of $`0.0`$$`0.5`$ and $`2.0`$$`3.0`$. These curves may be compared to those for bins 1 and 5 in Fig. 15. The scale dependence variations are generally within the bands formed by the 40 eigenvector basis sets, indicating that the uncertainty introduced by the scale dependence of the jet cross sections is less than or on the order of that coming from the experimental errors.

### III.3 Threshold resummation

As noted in the introduction to this section, there are potentially significant contributions to the jet cross resulting from phase space limitations on gluon emission as one goes to large $`E_T`$ and large rapidity. Such corrections can be treated by threshold resummation techniques resum . In Ref. twoloop  these corrections for jet production were estimated at the two-loop level for central rapidity values and found to be small. The calculation has been repeated for the largest rapidity bin covered in the DØ data d01b  and the results are shown in Figs. 21 and 22 by the dashed lines. As was the case in the central region, the two-loop estimates of the next-to-next-to-leading-logarithm contributions are not large for the scale choices considered in this analysis.

To summarize the results of this section, the next-to-leading-logarithm jet calculations show the expected reduction in scale dependence, even for large values of $`E_T`$ in the highest rapidity bin in the existing data set. Nevertheless, the remaining scale variation does show up in an increased uncertainty in the gluon distribution at large values of $`x`$. This uncertainty is, of course, partially offset by the compensating scale dependence of the hard scattering cross section. Even so, the scale dependence of the next-to-leading-log calculations remains as a contributing factor for the overall uncertainty of the predictions for new processes. Finally, the next-to-next-to-leading-log corrections estimated by expanding the threshold resummation predictions for jet production remain small, even at high $`E_T`$ in the highest rapidity bin considered, lending support to the use of the next-to-leading-log expressions in the global fit.

## IV Room for New Physics in Run 1b

Hadron collider jet data have traditionally been used to constrain models of new physics, usually models of quark substructure. The measurements include dijet angular distributions that are insensitive to theoretical uncertainties due to parton distributions angdist . Inclusive jet cross sections or dijet mass distributions typically provide better probes of quark compositeness, but are sensitive to parton distributions at large $`x`$ which have large uncertainties dijetmass . Currently the best limits are provided by a hybrid approach that uses the ratio of dijets in different rapidity bins d0dijet . This approach has the advantage that many theoretical and experimental systematic errors cancel out in the ratio. It has the disadvantage that it removes some useful information, such as the absolute values ot the cross sections in the two rapidity intervals.

There are several reasons to revisit this issue. First, the published limits are based on parton distributions several generations old, and the best-fit gluon distribution at large $`x`$ is significantly changed. Second, the recent jet measurements at large rapidity provide constraints on the large $`x`$ partons while being largely insensitive to new physics. Finally, and most important, the global analysis machinery has never been used in these analyses to refit the parton distributions in the presence of a new component to the QCD calculation. If the constraints on the parton distributions are not sufficient, it is possible that new physics could hide inside the PDF uncertainties, or current limits may be weakened.

Compositeness signals are typically parameterized by a mass scale $`\mathrm{\Lambda }`$ which characterizes the quark substructure coupling. The coupling is approximated by a four-Fermi contact interaction giving rise to an effective Lagrangian. Only the term describing left-handed coupling of quarks and anti-quarks has been calculated comp ,

$$_{qq}=A(2\pi /\mathrm{\Lambda }^2)(\overline{q}_L\gamma ^\mu q_L)(\overline{q}_L\gamma _\mu q_L)\text{where}A=\pm 1.$$
(12)
The simplest squared amplitude is the quark anti-quark s-channel annihilation process, with the first term the normal QCD interaction and the second due to substructure:

$$A(u\overline{u}d\overline{d})^2=\frac{4}{9}\alpha _S^2(Q^2)\frac{(\widehat{u}^2+\widehat{t}^2)}{\widehat{s}^2}+\frac{A\widehat{u}}{\mathrm{\Lambda }^2}.$$
(13)
Since $`\widehat{u}`$ is negative, constructive (destructive) interference in this process occurs when $`A=1(+1)`$. This is generally true for the other scattering terms as well, which will be included in the calculations. As in the experimental papers, we will use the quark compositeness implementation in PYTHIA pythia  to form the ratio of QCD+substructure/QCD, then multiply this ratio by the NLO QCD calculation for each iteration of parton distribution as the PDF parameters are being fit. As expected, we find larger cross sections with the constructive interference choice $`A=1`$, approximately a factor of 10 larger. However, the ratio of the first two eta bins used by DØ is slightly larger with the destructive interference term $`A=1`$, and this is what is used for the current best limit of $`\mathrm{\Lambda }>2.7`$ TeV d0dijet . Therefore we have used the destructive interference choice for illustration in the rest of this paper.

Figure 23 shows the cross section for jet production including a contact interaction, for three values of $`\mathrm{\Lambda }`$: $`1.6`$, $`2.0`$ and $`2.4`$ TeV. The cross section is plotted as a ratio to the pure QCD prediction with the CTEQ6.1M parton distributions. The range of kinematic parameters are the same as the DØ data, and the DØ data is superimposed on the plots with error bars of statistical and systematic errors combined in quadrature.

If the effect of the contact interaction is smaller than the PDF uncertainty of the QCD prediction then the data cannot rule out the compositeness model; any difference between theory and data could be explained either as a PDF effect or as a sign of new physics. However, if the compositeness model disagrees with the data by an amount that is larger than the PDF uncertainty, the model is ruled out. The PDF uncertainty of the QCD prediction is comparable to the experimental error bars (see Fig. 15), so we may use those error bars to judge whether the effect of the contact interaction is smaller than the PDF uncertainty. The contact interaction is seen to be ruled out for $`\mathrm{\Lambda }=1.6`$ TeV, primarily from the disagreement with the central rapidity interval, $`0<\eta <0.5`$. The cross section at large rapidity ($`\eta >1.5`$) is not sensitive to the contact interaction.

The largest two $`\mathrm{\Lambda }`$ values, $`\mathrm{\Lambda }=2.0,2.4`$ TeV, are not inconsistent with the DØ data on inclusive jet production because the effect of the contact interaction is within the PDF uncertainty. The global fit for $`\mathrm{\Lambda }=2.4`$ TeV is even slightly better than for the pure QCD model, although the difference in the overall quality of the fit ($`\mathrm{\Delta }\chi _{\mathrm{global}}^2=9`$) is much less than our standard tolerance for PDF uncertainty ($`\mathrm{\Delta }\chi _{\mathrm{global}}^2100`$).

If there is a contact interaction, then the CTEQ6.1M PDF’s would be inaccurate because the data were fit to a pure QCD model. Therefore a complete analysis of the contact interaction requires that the PDF’s be refit, comparing the data to the theory with compositeness. Figure 24 shows the final result of this analysis. The cross section for jet production with the DØ kinematic parameters is again shown for three values of the $`\mathrm{\Lambda }`$ parameter ($`1.6`$, $`2.0`$ and $`2.4`$ TeV) as a ratio to the pure QCD model with CTEQ6.1M. In Fig. 24 the PDF’s for each $`\mathrm{\Lambda }`$ value come from a global fit that includes the contact interaction in the theory of the jet cross section. The smallest value of $`\mathrm{\Lambda }`$ is now more consistent with the data, but still ruled out.

The limit derived from this analysis, $`\mathrm{\Lambda }>1.6`$ TeV, is not as strong as those from either the CDF or DØ analyses that use angular distribution information angdist . This is a result of the remaining uncertainty in the PDFs, and the fact that in this model of new physics the angular distributions are quite different than QCD. Other deviations from the Standard Model might change the absolute cross sections while maintaining a more QCD-like angular distribution. In this case this analysis is the most in-depth study so far of allowed deviations in jet cross sections due to new physics.

The limit derived from this analysis is also less than that obtained by the DØ collaboration when using the ratio of the dijet cross sections in the two rapidity bins ($`0.00.5`$ and $`0.51.0`$). However, it should be noted that at least part of the higher limit from DØ may be due to the dijet cross section for $`0.5<\eta <1.0`$ being somewhat larger than the standard model prediction, while the dijet cross section for $`0.0<\eta <0.5`$ is in better agreement with the SM prediction. The larger the ratio of the $`0.51.0`$ to $`0.00.5`$ rapidity intervals, the stronger the compositeness limit. The detailed information about the absolute values of the cross sections in each rapidity bin has not been used.

Table 2 shows how the $`\chi ^2`$’s are affected by the contact interaction.

## V Predictions for Run 2

The increase in the center-of-mass energy (from 1.8 to 1.96 TeV) and the increased luminosity expected in Run 2 (from approximately $`100`$ pb-1 to $`2`$ – $`15`$ fb-1 will result in a dramatically larger kinematic range for measuring jet production, and thus for searching for possible new physics as well. Also, in Run 2, the new CDF endplug calorimeter will allow for an extension of measurements of the inclusive jet cross section to the forward rapidity region, similar to the analysis that has been performed by DØ in Run 1b.

The predictions for the Run 2 jet cross section as a function of transverse momentum (using the DØ rapidity intervals) are shown in Fig. 25. With $`2`$ fb-1 (Run 2a) a measurement of the inclusive jet cross section out to a transverse energy of $`600`$ GeV is possible, while with $`15`$ fb-1 (Run 2b) the measurement can be extended to almost $`650`$ GeV. The predictions are displayed on the left graph in the form of the curves for the 40 eigenvector basis sets of CTEQ6.1M. The central prediction and overall uncertainty band from the Master Equation (10) are shown separately on the right graph. Figure 26 shows the uncertainty range of the Run 2 cross section by displaying the ratios of the eigenvector basis sets to the central prediction. Near the kinematic limit, the PDF uncertainties allow the cross section predictions to be as much as a factor of 2 larger than those from the central fit (CTEQ6.1M), again due to eigenvector 15.

For the Run 2 jet cross section predictions, we have used the Run 2 cone jet algorithm run2work ; ellis  which uses 4-vector kinematics both to define the jet and to specify its transverse momentum and rapidity. This algorithm has been adopted by both the CDF and DØ collaborations for their inclusive Run 2 jet results. Note that this algorithm specifies the transverse momentum $`p_T`$ rather than the transverse energy $`E_T`$.777The Run 2 jet algorithm results in cross section predictions for CDF that are $`7`$% lower than those using the Run 1 algorithm, over most of the kinematic range.

For completeness we also show the plots of the predictions of the inclusive jet cross section at the Tevatron Run 2 for the CDF choice of rapidity bins. Figure 27 shows the cross section as a function of $`p_T`$ on a log scale. Figure 28 shows the uncertainty band in the form of the ratios of $`d\sigma /dp_T`$ for the 40 eigenvector basis sets compared to the central prediction.

While the statistical errors in Runs 2a and 2b will be small, so that the jet cross section can be extended to high $`p_T`$, the experience from Run 1b shows that the systematic errors will have a dominant effect on the significance of the results. However, we can anticipate that the PDF’s, especially the gluon distribution, will be pinned down at large $`x`$ by these data. The experimental collaborations will make the data most useful for global analysis by providing detailed information on systematic uncertainties.

### V.1 The ratio of Run 2 to Run 1b cross sections

As mentioned previously, the Run 2 center-of-mass energy $`\sqrt{s}`$ is higher than that from Run 1b. The increase in $`\sqrt{s}`$ will have only a modest impact on most physics cross sections. However, it will lead to a rather large increase in the production of high-$`E_T`$ jets. It may be useful to examine the ratio of jet production (for the same $`E_T`$ values) at the two different energies. In such a ratio, many of the theoretical (and experimental) errors will cancel. The Run 2 to Run 1b ratio for the 40 sets of PDF’s is shown in Fig. 29 for the central rapidity range of CDF ($`0.1<|\eta |<0.7`$). As can be observed the ratio has a rather narrow theoretical error band from PDF uncertainties, i.e., smaller uncertainty than that of the absolute prediction of the inclusive cross section.888Here we have multiplied our Run 2 predictions by a factor of $`1.07`$ so that the comparison effectively uses the Run 1 algorithm for both the numerator and denominator. This is done for convenience since CDF will present its first results for the Run 2 jet cross section using the Run 1 algorithm.

## VI Inclusive jet production at the LHC

The increase of the center-of-mass energy to $`14`$ TeV at the LHC will result in a dramatically larger accessible kinematic range. Inclusive jet cross sections will be measured out to transverse momentum of $`5`$ TeV/$`c`$ in the central rapidity range and out to $`1.5`$ TeV/$`c`$ in the forward rapidity region. Figure 30 shows the predictions with uncertainty ranges for three rapidity intervals. Figure 31 shows the ratios of the 40 eigenvector basis set predictions to the central prediction. The cross section uncertainties near the kinematic limit at the LHC due to PDF uncertainties are similar in magnitude to those obtained for Run 2 at the Tevatron. Again the extremes of the predictions are provided by the eigenvector basis sets $`+15`$ and $`15`$, which correspond to extremes of the gluon distribution at large $`x`$.

The jet cross section at the LHC will be a discovery mode for new physics. Supersymmetry would produce an enhanced cross section from production of gluino and squark jets. Extra dimensions would show up from production of Kaluza-Klein modes of the graviton. So the pure QCD prediction will be compared to the data and any difference will be the first evidence for these or other new physics scenarios.

## VII Conclusions

Jet production serves as a probe of the highest $`Q^2`$ scales accessible at a hadron-hadron collider. NLO QCD provides a good description of the inclusive jet cross sections measured by CDF and D0 in Run 1 at the Tevatron, provided that the gluon distribution is enhanced at high $`x`$ relative to earlier estimates. Such an outcome is now the standard consequence of the global fitting procedure due to the statistical power of the Tevatron jet data, and the preference for a larger gluon in all rapidity regions. The uncertainties on the jet cross section predictions still remain large, however, primarily due to the remaining uncertainties on the gluon distribution. In most regions, they are comparable to the experimental systematic uncertainties.

The standard choice for the renormalization and factorization scales in evaluating inclusive jet cross sections has been $`E_T^{jet}/2`$. We have re-done the global fitting using alternate scales of $`E_T^{jet}`$ and $`2E_T^{jet}`$ and have found changes to the gluon distribution that are well within the gluon PDF uncertainty band. The effects of the resummation of threshold logarithms has been calculated for the first time in the forward rapidity regions. As for the central rapidity region, the effects of the threshold logarithms are relatively small.

The increase of the center-of-mass energy to 1.96 TeV (from 1.8 TeV) and the expected greater integrated luminosity will lead to an expanded kinematic range for Run 2 at the Tevatron. Predictions have been made for both the inclusive jet cross section and its uncertainty for Run 2 and for the LHC. The PDF uncertainty on the jet cross section near the limits of the data reach in both cases are similar to those obtained for Run 1.

We have carried out an exercise where we have examined the allowed range for any new physics such as compositeness, given the current scale of uncertainties on the jet cross section predictions. Given the Run 1b jet data, any new physics must have a scale larger than 1.6 TeV. When data is available from Run 2 at the Tevatron, and from the LHC, analogous analyses will be used to identify evidence for new physics or place even stricter limits on new physics.

Appendix

One component of the theoretical uncertainty of any perturbative QCD calculation is the dependence of the answer on the various renormalization and factorization scales. One virtue of using NLO calculations in global determinations of parton distributions is the reduction in scale dependence usually exhibited by such calculations. The following brief discussion outlines the origin of the reduction of the scale dependence which is often observed with NLO calculations.

Consider a large transverse momentum process such as the single jet inclusive cross section involving only massless partons. Furthermore, in order to simplify the notation, suppose that the transverse momentum is sufficiently large that only the quark distributions need be considered. In the following, a sum over quark flavors is implied. Schematically, one can write the lowest order cross section as

$$E\frac{d^3\sigma }{dp^3}\sigma =a^2(\mu )\widehat{\sigma }_Bq(M)q(M)$$
(14)
where $`a(\mu )=\alpha _s(\mu )/2\pi `$ and the lowest order parton-parton scattering cross section is denoted by $`\widehat{\sigma }_B`$. The renormalization and factorization scales are denoted by $`\mu andM`$, respectively. In addition, various overall factors have been absorbed into the definition of $`\widehat{\sigma }_B`$. The symbol $``$ denotes a convolution defined as

$$fg=_x^1\frac{dy}{y}f(\frac{x}{y})g(y).$$
(15)
When one calculates the $`𝒪(\alpha _s^3)`$ contributions to the inclusive cross section, the result can be written as

$`\sigma `$ $`=`$ $`a^2(\mu )\widehat{\sigma }_Bq(M)q(M)`$ (16)
$`+`$ $`2a^3(\mu )b\mathrm{ln}(\mu /p_T)\widehat{\sigma }_Bq(M)q(M)`$
$`+`$ $`2a^3(\mu )\mathrm{ln}(p_T/M)P_{qq}\widehat{\sigma }_Bq(M)q(M)`$
$`+`$ $`a^3(\mu )Kq(M)q(M).`$
In writing Eq. (16), specific logarithms associated with the running coupling and the scale dependence of the parton distributions have been explicitly displayed; the remaining higher order corrections have been collected in the function $`K`$ in the last line of Eq. (16). The $`\mu `$ dependence of the running coupling is obtained from

$$\mu \frac{a(\mu )}{\mu }=\beta (a(\mu ))$$
(17)
where $`\beta =ba^2(1+ca)`$ with

$$b=\frac{332f}{6}\mathrm{and}c=\frac{15319f}{2(332f)}$$
(18)
where $`f`$ denotes the number of flavors. The scale dependence of the quark distributions is given by the nonsinglet DGLAP equation DGLAP

$$M\frac{q(x,M)}{M}=a(M)P_{qq}q(M).$$
(19)
Now, the NLO expression in Eq. (16) should be independent of the choice of scale up to corrections of order $`a^4`$. That this is indeed the case is easy to see. For convenience, set $`M=\mu `$ and calculate the logarithmic derivative of $`\sigma `$ with respect to $`\mu `$,

$`\mu {\displaystyle \frac{\sigma }{\mu }}=`$ $`2a(\mu )[ba^2(\mu )(1+ca(\mu ))]\widehat{\sigma }q(\mu )q(\mu )`$
$`+2a^3(\mu )P_{qq}\widehat{\sigma }_Bq(\mu )q(\mu )`$
$`+2ba^3(\mu )\widehat{\sigma }_Bq(\mu )q(\mu )`$
$`2a^3(\mu )P_{qq}\widehat{\sigma }_Bq(\mu )q(\mu )+\mathrm{}`$
$`=0+𝒪(a^4).`$ (21)
The derivative of the lowest order contribution yields the first two lines of Eq. (21) while the derivative of the logarithms in the second and third lines of Eq. (16) give rise to the last two lines in Eq. (21) which cancel the $`\mu `$ dependent terms coming from the derivative of the lowest-order contribution. All other derivatives yield contributions which are of order $`a^4`$. Thus, the structure of the NLO calculation is such that the scale dependence is cancelled up to terms of order $`a^4`$.

The preceding discussion shows that one might expect a reduced scale dependence for the NLO calculation as opposed to that observed at lowest order. However, the uncancelled $`𝒪(a^4)`$ and higher terms will still result in some scale dependence. It is easy to understand the systematics of the remaining dependence for the case at hand. Consider, first, the lowest-order contribution shown in the first line of Eq. (16). At large values of $`p_T`$ as $`x_T=2p_T/\sqrt{s}`$ tends to one, the scaling violations of the parton distributions lead to a decrease in the cross section as the scale increases. In addition, the running coupling decreases as the scale is increased. Therefore, at fixed $`p_T`$ the lowest order contribution exhibits a monotonically decreasing behavior as the scale increases.

The second line of Eq. (16) contains the contribution which partially offsets the scale dependence of the running coupling. Apart from the explicit logarithm, this contribution has much the same scale dependence as the lowest order contribution. Thus, for $`\mu <p_T`$ this term is negative, becoming positive as the scale is increased beyond $`p_T`$.

The third line contains the contribution which partially offsets the scale dependence originating in the parton distributions. In the kinematic region under consideration, the convolution of the splitting function $`P_{qq}`$ with the lowest contribution gives a negative result, reflecting the fact that the parton distributions decrease with increasing scale at large values of $`x`$. Therefore, this contribution is negative when the factorization scale $`M`$ is chosen to be less than $`p_T`$.

Finally, the remaining term on line four of Eq. (16) has much the same scale dependence as the lowest order term, as the dominant scale dependence has already been exhibited in the terms on lines two and three. Thus, lines one and four give contributions which decrease monotonically with increasing scale while lines two and three start out negative, reach zero when the scales are equal to $`p_T`$, and are positive for larger scales. This behavior explains why the typical NLO scale dependence in the high-$`p_T`$ region yields a result which first increases with increasing scale, reaches a maximum, and then decreases. For many high-$`p_T`$ processes involving massless partons this maximum occurs near $`\mu =M=p_T`$.

Next, it is interesting to consider how the location of the maximum changes with $`p_T`$. Suppose that one goes to small values of $`x_T`$ such that the typical values of $`x`$ in the parton distributions is around .2 or smaller. In this region the scaling violations change sign, corresponding to a growth of the parton distributions with increasing scale, rather than a decrease. The contribution of the third line of Eq.(1) then starts to compete against that of the second line. Accordingly, one must go to smaller values of $`\mu /p_T`$ before seeing the turnover and decrease with decreasing scale. That is, the position of the maximum at fixed $`p_T`$ moves to lower values of $`\mu /p_T`$ as $`p_T`$ decreases. The converse is also clearly true; as $`p_T`$ increases the position of the maximum moves to higher values of $`\mu /p_T`$. Thus, if one were to try to find an optimized scale at each value of $`p_T`$, the pattern that would emerge would be that the optimal scale would be a fractional multiple of $`p_T`$ at lower $`p_T`$ values and would increase to a multiple of $`p_T`$ greater than one at higher $`p_T`$ values. This pattern is precisely what is observed when using the single scale version of the Principle of Minimal Sensitivity PMS .

Another manifestation of this pattern concerns the behavior of the “K factor”, defined as the ratio of the NLO to LO result. At a fixed value of $`p_T`$, the K factor will first rise and then fall with increasing scale. Suppose that at some $`p_T`$ value the maximum occurs for the choice $`\mu =M=p_T/2`$. Then, as one goes to higher $`p_T`$ values the previous discussion shows that the maximum will move to a higher scale choice such as $`p_T`$ or even $`2p_T`$. At the larger $`p_T`$ values it will then appear that the K factor may increase with increasing scale. This is just the behavior seen in the calculations done for the high rapidity D0 jet data.


text:
# Analysis of microwave propagation in arrays of dielectric cylinders

## Abstract

We rigorously compute the propagation and scattering of microwaves by regular and defected arrays of dielectric cylinders in a uniform medium. Comparison with the previous experimental results is made, yielding agreements in the transmission coefficients for a certain range of frequencies. The results also show that localized states are possible for defected arrays. However, the experimentally claimed localized states are not observed. Moreover, comparison is made with previous theoretical results. The agreements and differences are pointed out.

Waves surround us. Direct human communication is mainly conveyed by acoustic waves, and is enriched by gestures which are passed into our eyes through optical waves. Nowadays electronic waves are also everywhere in our daily experiences such as audio & video systems, computers, and Nintendo games. Among many interesting properties pertinent to waves, the phenomenon of wave localization is probably most intriguing and remains as one of the unsolved mysteries of last century.

The concept of wave localization was first proposed by Anderson nearly half a century ago, initially for electronic systems Anderson58 , then extended to classical waves including microwaves to be addressed in this paper. It refers to situations that purely due to multiple scattering of waves by disorders, transmitted waves are confined or frozen in space and the wave envelopes decay nearly exponentially along any directionLee . Since its inception, wave localization has attracted substantial efforts from all areas of expertise, signified by the great body of literatures ranging from small scales such as electronic systems, nano-structured materials, to large scales such as seismic and ocean waves, as reviewed and summarized in, for instance, Refs. Lee ; Ishimaru ; Meso ; John ; Ad ; Datta ; Nano ; EA ; loc ; Seis ; ocean . In spite of the efforts, however, some fundamental issues concerning wave localization are still open Rama .

A few prominent problems in the investigation lie in the localization in two dimensional (2D) random media. First, although it has been conjectured, thought to be genuinely valid, by the scaling analysisgang4  that all waves are localized in 2D random media for any given amount of disorders or impurities, the experimental observation of 2D localization is scarce. Even for the limited experimental results that are available McCall , the disparity Meade ; Ho  still remains in the interpretation of the results. A main concern is whether the localization in 2D has indeed been observed or not Meade . Some authors further pointed out that there still lacks a definite observation of 2D wave localization, particularly for classical waves Sigalas . In addition, how waves are localized in 2D and whether the localization features can be explained or fitted by the existing theory DV  are among a few important issues yet to be addressed.

It would be too ambitious to answer all the unsolved issues dwelling on 2D localization in a single endeavor. However, a further look at some of previous experimental observations of 2D localization, in the hope of clarifying some confusions, seems more realistic and imperative, and will therefore be the main theme of the present paper.

Here we would like to concentrate on the experimental observation of microwave propagation in 2D dielectric lattices, reported by McCall et al. McCall . The reason is not only because the experiment is the first ever reported, but it has stimulated so many later investigations, and its importance could be inferred from the constantly increasing citation number. In McCall , the authors measured transmission of microwaves through regular and defected arrays of dielectric cylinders. The transmission bands and photonic band gaps for the regular arrays are found to be in excellent agreement with theoretical predictions. For the defected arrays, the authors suggested the observation of a localized state. This observation, however, has been contrarily discussed and supported by two independent theoretical inspections Meade ; Ho . In this paper, we will use the standard multiple scattering approach first formulated by Twersky Twersky  to further examine the experimental results. We will show that localized states are indeed possible for the defected arrays. However, the experimentally claimed localized states are not observed. Comparison will also be made with previous theoretical results. Agreements and differences are found and discussed. Some contraries in the previous analysis are clarified.

The systems considered here are from Ref. McCall  and they are composed of almost ideal arrays of low-loss high dielectric constant cylinders. The cylinders are placed in parallel in an uniform medium to form a square lattice of lattice constant denoted as d. The ratio of dielectric constants between the medium and cylinders is 9. Two types of measurement are carried out in McCall . One is to measure the E-polarized microwave transmission across the samples which are shaped as a rectangular slab with thickness of 9d and width of 18d and the lattice constant is 1.27 cm. The other is to measure the intensity distribution at 11.2 GHz for a sample size of 7.6 cm x 7.6 cm with a cylinder being removed from the center of the sample, i. e. defected array. However, the lattice arrangement for the second measurement differs from the first in that the lattice constant is 1.59 cm; such an inconsistency in the lattice arrangements has been discussed in Meade . The authors in McCall  suggested that the localized states are observed for the defected array.

In our computation, we will take the physical parameters from McCall . However, we will consider larger lattice sizes when appropriate to ensure the stability of the results. Both lattice arrangements, i. e. d = 1.27 and 1.59 cm corresponding to filling factors 0.4487 and 0.2863 respectively, will be considered. The cylinder radius is 0.48 cm. In addition to computing the transmission across rectangular slabs, for comparison purpose we also compute the transmitted intensity when the source in placed inside a sample that takes a circular shape. The spatial distribution of energy will also be calculated to compare with the experiment. Moreover, the phase diagram method Emile  will be used to investigate a coherence phenomenon associated with localization. To be consistent with the experiment, we only consider the E-polarized microwaves.

Figure 1 shows the band structures for the two 2D square lattices aforementioned. The results on the left panel reproduce that in McCall , confirming our numerical codes. For later discussion, the band structures of lattice constant 1.59 cm are also plotted. Here we see that the two results are similar for frequencies below the second gap. There are three prominent complete bandgaps for d = 1.27, while there are only two complete bandgaps for d = 1.59.

The transmission intensity is plotted as a function of frequencies in Fig. 2. The situations in (a1), (a2) and (a3) are considered in McCall . From this figure, we observe the following. First for d = 1.27 cm: (1) The results for situations considered by the experiment agree both qualitatively and quantitatively well with the measurement for frequencies up to the beginning of the third bandgaps at about 14 GHz; (2) While we did not plot the case of propagation along the  direction, we plot the case of the transmission from a source located inside the sample, with a cylinder being deleted at the center (a4). The result is in accordance with the case that the transmission is measured across the sample. In both cases, there is a transmission peak inside the valley at 11.27 GHz, slightly above the measured value, which is outside the complete gap, but within the gap along the  direction. The reason for the difference in the location and magnitude of the peak from the experiment is due to the resolution of plotting. (3) There are some significant differences between the theory and experiment for frequencies starting from the third gap. For these frequencies, we observe clearly a well defined reduction in the transmission, well matching the band structures in Fig. 1, and the reduction is even deeper than the second gap. The experimental results, however, show some smeared reduction. Further calculations point to a few factors contributing to the disagreement, such as disorders, inhomogeneities in the radii of the cylinders, and the size of samples.

The results of the transmission through a slab with d = 1.59 cm are depicted on the right panel of Fig. 2. These diagrams show that the transmission is significantly prohibited in the second gap, consistent with the band structure calculation. There is also a peak in the reduction regime. The frequency at which the energy distribution is measured is outside this forbidden regime.

An important feature in Figs. (a3) and (a4) is the peak inside the second gap along the  direction but outside the complete bandgap, i. e. within the pseudo gap. McCall et al. McCall  attempted to explain this peak as the result of a localized state. To support this view, McCall et al. further measured the energy spatial-distribution at this peak and attempted, believed to be successful, to show that there is indeed a localized wave at this frequency. Unfortunately, this time they used the different lattice with d = 1.59 cm. The inconsistency has been pointed out in Meade . Meade et al. Meade  argued that the peak is due to the resonance of extended states, and the apparent localized wave shown by Fig. 4 in McCall  is superficial and is due to the resonance.

To resolve the disagreement between the measurement and the previous theoretical analysis, we show in Fig. 3 the 2D spatial distribution of energy at 11.27 GHz for the lattice of d = 1.27 cm and 1.59 cm with one cylinder being taken out at the center of the sample. We also plot the associated phase diagram. The source transmits at the center.

According to Refs. Emile ; Ye , there will appear a coherence for the phase of the localized wave. We brief this feature here. The energy flow of EM waves is $`\stackrel{}{J}\stackrel{}{E}\times \stackrel{}{H}`$. By invoking the Maxwell equations to relate the electrical and magnetic fields, we can derive that the time averaged energy flow is $`<\stackrel{}{J}>_t\frac{1}{T}_0^T𝑑t\stackrel{}{J}|\stackrel{}{E}|^2\theta ,`$ where the electrical field is written as $`\stackrel{}{E}=\stackrel{}{e}_E|\stackrel{}{E}|e^{i\theta }`$, with $`\stackrel{}{e}_E`$ denoting the direction, $`|\stackrel{}{E}|`$ and $`\theta `$ being the amplitude and the phase respectively. It is clear that when $`\theta `$ is constant, at least by spatial domains, while $`|\stackrel{}{E}|0`$, the flow would come to a stop and the energy will be localized or stored in the space. The phase can be represented by a vector defined as $`\stackrel{}{v}(\stackrel{}{r})=\stackrel{}{e}_x\mathrm{cos}\theta +\stackrel{}{e}_y\mathrm{sin}\theta `$, and the phase vectors at various space points can be drawn on the 2D coordinates. The coherence in the phases is a unique indication of localization.

Fig. 3(a) shows that the wave is indeed localized at 11.27 GHz, as suggested in McCall  for the case of d = 1.27 cm. Here is shown that the energy is mostly confined within the lattice and the prescribed phase coherence prevails. The localization due to scattering by the cylinders is so strong that the magnitude of the intensity surpasses the energy directly from the source. However, this localized wave is not what was observed in the experiment. With the experimental setup of d = 1.59 cm, the experimentally claimed localization is actually absent, as suggested by Fig. 3 (b) where the phase vectors point to different directions indicating the flow of energy. Nevertheless, the ripples in the energy distribution mimics that observed in the experiment. We note that for all cases, the sample size considered here is larger than that in the experiment.

To further confirm that localization is indeed possible at 11.27GHz, we also considered the case of a complete random array which is otherwise regular with d = 1.27 cm. The results are shown in Fig. 4. Here the wave localization indeed appears. The energy is localized at the center of sample with phase vectors aligning along one direction. Due to the randomness, the peak value of the localized energy is reduced significantly.

We further consider two other frequencies for the case of regular array d = 1.27 cm with one cylinder being removed: one is at 11.45 GHz, slightly above 11.27 GHz but still within the pseudo gap along the  direction, and the other one lies at 13 GHz between the second and third gaps. The results are illustrated in Fig. 5. They show that the localization is apparent for 11.45 GHz while absent for 13 GHz. It is interesting to notice that the pointing direction of the phase vectors is mostly reversed from the case with the lower frequency, comparing to Fig. 3 (a). It is also worth to noting that the localization shown by both Figs. 3 and 5 occurs outside the complete band gap. These results will eventually give insight to the relation between gap and localization. The absence of localization in defected arrays may be due to the finite size of samples, in the view of previous studies. But it might also hint that waves are not necessarily always localized in 2D, as suggested before Emile .

Finally, to examine the difference between the theoretical and experimental transmission results when comparing the left panel of Fig. 2 and Fig. 3 in McCall , we investigate the effect of disorders. For brevity, we only show the effect of the positional disorder which is introduced by moving the cylinders around their regular lattice position to a certain extent. Here the disorder is measured as 12% of random displacement of the cylinders within a circle around their original regular positions with regard to the lattice constant. The result is shown in Fig. 6. The result indeed shows that the disorder tends to degrade the transmission reduction in the third gap to a certain degree, while the disorder effect is not so significant for the second gap.

In summary, in this paper we have shown a numerical analysis of the previous experimental results on microwave propagation through arrays of dielectric cylinders embedded in a uniform medium Meade  by the standard multiple scattering method. Comparison has been also made with other theoretical results. The results reveal that localized states are possible for defected arrays. However, the experimentally claimed localized states are not observed. The disorder effects are also presented. The present work is hopped to stimulate further experimental observations.

Thanks to NSC and NCU for supports.