On the characterisation of the underlying event

The measurement of the underlying event (UE) and its separation from hard interactions in hadron-collider events is a conceptually and practically challenging task. We develop a simple, mostly analytical toy model for the UE in order to understand how different UE measurement approaches fare on the practical aspects of this problem, comparing the traditional approach used so far at Tevatron with a recently proposed"jet-area/median"approach. Both are found to perform comparably well in measuring average properties of the UE, such as the mean transverse momentum flow, but the jet-area/median approach has distinct advantages in determining its fluctuations. We then use the jet-area/median method to investigate a range of UE properties in existing Monte Carlo event-generator tunes, validating the main results of the toy-model and highlighting so-far unmeasured characteristics of the UE such as its rapidity dependence, as well as its intra- and inter-event fluctuations and correlations.


Introduction
The "underlying event" (UE) in high-energy hadron-hadron collisions can be thought of as the low transverse momentum (p t ) part of the event activity that is not naturally associated with the hard interaction. Despite being a low-p t phenomenon, the underlying event has a large impact on high-p t physics at hadron colliders. For example in measurements of the inclusive jet spectrum at p t ∼ 50 GeV it can affect the result by up to 50% [1]. It also biases kinematic reconstructions (for example in top mass measurements), degrades their resolution and affects the efficiency of isolation criteria that enter into the experimental identification of particles like photons and electrons. It has even been suggested that certain new physics scenarios might show up in "anomalous" characteristics of the underlying event [2]. For these reasons it is important to have a good understanding of its properties.  Figure 1: Views of a hard scattering plus UE; (a) shows multiple 2 → 2 interactions, as incorporated in the most successful models; (b) illustrates how a collinear splitting in the initial-state can lead to correlation between the partons involved in a double 2 → 2 scattering; from (c) one sees the similarity with a perturbative 1-loop 2 → 4 diagram; and (d) represents a BFKL-inspired picture for the UE.
The purpose of this article is to investigate some different ways in which the underlying event can be measured and/or constrained experimentally. Such measurements enter into tunes of Monte Carlo models of the UE [3][4][5][6][7][8][9][10]. They can also serve as an input to analytical methods of accounting for the average UE correction to a jet's transverse momentum [11,12] and to approaches that correct for the UE on a jet-by-jet and event-by-event basis [13], as well as for related work that seeks to optimise jet definitions.
A difficulty in discussing the measurement of the UE is that there exists no good definition of what the UE actually is, or how to distinguish it, in a conceptually unambiguous way, from the hard interaction. For instance, the most successful phenomenological models of the UE involve multiple 2 → 2 scattering as in fig. 1a ("multiple parton interactions" -MPI). Very simply, they supplement the one hard interaction in the event with multiple other lower-p t interactions, whose multiplicity is determined by the 2 → 2 cross section, regulated with a low-p t cutoff of the order of a couple of GeV. Fig. 1a can only be part of the picture because some of the partons entering multiple 2 → 2 scatterings are necessarily correlated, e.g. due to energy conservation (cf. ref. [14]), or because they can have a common origin from an initial-state collinear splitting, as in fig. 1b. However the contribution of fig. 1b is itself also part of the 2 → 4 1-loop scattering diagram fig. 1c (as discussed for example in [15]), which is relevant as a N 4 LO correction to the dijet cross section or at N 3 LO in its interference with tree-level 2 → 4 scattering. This means there are non-trivial questions of double-counting between multiple parton interaction and perturbative higher orders. In addition, the radiation that fills the event is not bound to come just from 2 → 2 scatterings, but may also arise from BFKL type configurations which can involve (multiple) chains with low-p t emissions spread in rapidity, as in fig. 1d (some work towards modelling this while retaining consistency with the total cross section is to be found in [16]). Though these ways of viewing the UE represent just a subset of the diversity of physics considerations that are of potential relevance for its modelling (more detailed reviews are to be found in Refs. [10,17]), they do illustrate the difficulties that arise in ascribing unambiguous physical meaning to it.
Given this complexity in discussing what the UE might be, how are we to go about measuring it? A feature present in most models is that, on some low p t scale, UE activity fills the whole event. One way then of characterising the UE is to say that it is whatever physical effect fills most of the event with radiation. To help understand what implications this picture has for UE measurements, we shall develop (section 3) a semi-analytical two-component toy model: one component will be purely soft and dispersed across the event, corresponding to the UE, while the second component will involve the hard scattering and a simple approximation for the perturbative radiation with which it is associated. Though our toy model is undoubtedly too simple to fully reflect reality, the fact that we know exactly what goes into it will make it quite powerful: it will, for example, allow us to take different UE-measurement approaches and examine to what extent their results are affected both by the radiation associated with the hard scattering, and by the techniques used to limit that hard contamination. This can be investigated both for averaged quantities, and for event-by-event extractions of information about the UE. A number of the results will be given in analytical form, in terms of the characteristics of the UE and the hard scattering and of the parameters of the measurement methods. This will give insight into the compromises that arise when measuring the UE, especially when extending existing methods to the greater phase space that is available at LHC relative to Tevatron.
The two UE measurement methods that we shall investigate are both reviewed in section 2. One, which we call the "traditional" approach (see for example [18,19]), is currently the default approach for most UE studies. Another, the jet-area/median based approach of [12,13], was originally developed for evaluating pp pileup (and backgrounds in heavy-ion collisions, as used for example in [20]), but may also have benefits for UE studies. The basic results for the two approaches will be derived in section 3. Some readers may prefer to skip most of this section and read just the final summary of these results, as presented in section 3.6.
For the purpose of determining the quantity that we call ρ, the UE's mean transverse momentum (p t ) flow per unit rapidity-azimuth area, both methods will turn out to have systematics that are under control and of similar magnitude, at about the 20% level (except at high p t for the traditional method). However it is also important to have knowledge of fluctuations of the underlying event (both intra and inter-event). Since it is the jetarea/median method that will prove to be the more robust tool for measuring them, it is this method that we will use when, in section 4, we look at a range of possibly interesting measurements of the UE. They will be carried out on Monte Carlo events simulated with Pythia 6.4 [21] and with Herwig 6.5 [22] with Jimmy 4.3 [4]. This part of the study will help validate the understanding developed with the toy-model, and illustrate determinations of the UE average p t flow as a function of rapidity, its event-to-event fluctuations, its intra-event fluctuations, and the degree of intra-event correlation. These observables go beyond the kinds of measurements that are commonly discussed for Tevatron or envisaged so far for the LHC and, as we shall see, there will be substantial differences between Jimmy and Pythia UE tunes on a number of them.
2 Overview of measurement approaches

Traditional approach
The main current approach [18] to measuring the properties of the underlying event involves considering the central part of the detector, say |η| < 1, where the pseudorapidity η is defined as η = − ln tan θ 2 . One tags events based on the presence of a jet (whose direction defines an azimuthal angle φ = 0), and divides the central part of the event into four blocks in azimuth: the "towards" region, typically |φ| < π/3, an away region 2π/3 < |φ| < π, and two transverse regions, covering π/3 < |φ| < 2π/3. This is illustrated in fig. 2 (left). Since the trigger and recoil jets will usually occupy the towards and away regions, one then restricts one's attention  Right: representative distribution of p tj /A j for the set of (real and ghost) jets in a single event, as used to determine ρ and σ in the area/median approach.
to the two transverse regions. There one measures the multiplicity of charged tracks above some transverse-momentum threshold as well as the total transverse momentum contained in the charged tracks (sometimes normalised per unit area, dηdφ). The results for the charged track multiplicity and charged p t flow are usually presented as averages across many events, as a function of the p t of the leading jet. One also sees measurements of the charged momentum flow as a function of the multiplicity.
Since there is a probability of order α s that at least one of the transverse regions is contaminated by perturbative radiation from the dijet event, which substantially affect the extracted information about the UE's p t flow. To work around this, it is usual to label the two transverse regions as TransMin and TransMax, respectively the less and more active of the two. The largest component of perturbative contamination should be restricted to the TransMax region, while TransMin should be less affected.
In the earliest variants of this "two-region" method [23,24], the two regions used for sampling the UE were actually placed at non-central rapidities, rather than central rapidities and transverse azimuth, and it was the total transverse energy flow that was considered rather than just its charged component. Another variant measured charged momentum flow in cones [25]. These differences reflect the freedom inherent to this method: the question of where to place the "transverse" regions, and the choice of their shape and size. In the rest of this article we will always assume that the transverse regions are well separated from the dijet system in an event, as shown in fig. 2 (left), and we will quote our results as a function of the area A Trans of each of the transverse regions.

Jet-area/median approach
In the jet-area/median approach, one first clusters the event with a Cambridge/Aachen (C/A) [26,27] or inclusive k t [28,29] type jet algorithm. To each jet, j, one attributes an "active" jet-area, A j , as described in more detail in [12]. This is calculated by adding a large number of "ghost" particles to the event (each with negligible p t ∼ 10 −100 GeV) and including them in the clustering. The area of a jet is then proportional to number of ghosts it contains. Some jets contain just ghost particles ("pure ghost jets") and are considered to have p t = 0.
A proposal [13] for a way to measure the level ρ of uniform background noise in an event is then to take it to be the median of the distribution of the p tj /A j for the ensemble of jets in that event as shown schematically in fig. 2 (right). The logic of the use of the median is that it is much less susceptible to contamination from outliers (i.e. hard perturbative jets) than the mean. In addition to measuring ρ one can also determine the intra-event fluctuations of the UE. We introduce a quantity σ, defined such that a fraction X/2 of jets satisfy where X = Erf(1/ √ 2) ≃ 0.6827 is the fraction of a Gaussian distribution within one standard deviation of the mean. 1 The approach to measuring σ is analogous in spirit to the use of the median for determining ρ. As shown in fig. 2 (right), it is one sided (i.e. just considering jets with p tj /A j < ρ). This choice has been made so as to limit contamination of σ from hard jets when the total number of jets is small.
This method was originally suggested in [13] as a way of measuring average pileup noise across an entire event (e.g. for |y| < 5). One difficulty when using it for the UE will come from the fact that the UE is significantly softer, and has relative fluctuations that are larger and less Gaussian than those in (say) 20 pileup minimum-bias events, both consequences of the lower density of particles.
Another point is that in measuring the UE it is important to obtain differential information on the UE's rapidity (y = ln E+pz E−pz ) dependence. This leads us to consider ρ(y) = median j∈jets, |y j −y|<δy where 2δy is the width of a rapidity window in which one carries out the measurement. The fact that one has a relatively limited number of jets in a given rapidity window will be one further challenge that we will face in using this method for studying the UE, because the relative impact of the presence of a hard jet in the region of interest is amplified by the small total number of jets. A modification of the method that helps address these difficulties is "hard jet removal," first employed in the STAR collaboration's [30] use of the techniques of [13] for estimating the (very large) UE for jet measurements in heavy-ion collisions. We will investigate the impact of this modification here for events with dijet topologies, in which case we will simply remove the two hardest jets in the event from the overall list of jets. A choice that is present in both pileup and UE measurements with the area/median method is that of the jet algorithm and jet radius R. Both the k t [28,29] and Cambridge/Aachen [26,27] algorithms are suitable options, because they produce jets whose area distribution is quite regular. In contrast, algorithms that give mostly conical jets (like antik t [31] and, to a lesser extent, SISCone [32]) tend not to be, because they fill in the "holes" between the cones with jets with very small areas, which can have unrepresentative p t /A values. The question of what R value to use is one of the freedoms of the method and will be discussed in the coming section.

A toy model
To understand the strengths and weaknesses of different UE-measurement approaches, it is helpful to consider events as consisting of two components: a low-p t ("soft") noise component, defined to be the UE, supplemented with hard jets from a perturbative scattering and associated higher-order corrections. We will introduce models for each of these two components and investigate how the methods behave when either of the two components is present alone and when both are present together. Our guiding principle in designing these models has been to keep them sufficiently simple as to be treatable analytically, while also maintaining a reasonable degree of realism. While our combined hard and soft models will not quite have the continuous transition between hard and soft components that is present in Monte Carlo MPI models, we will see that they nevertheless lead to certain signature behaviours in UE measurement methods that correspond nicely to what is observed in Monte Carlo simulations.

Low-p t component
As a simple model for the underlying event let us imagine that on average in a patch of unit of area there are ν particles, that the probability distribution for the number of particles n in a specific patch of area A follows a Poisson distribution, P n = (νA) n e −νA /n!, and that the single-particle transverse-momentum probability distribution is given by where µ is the mean transverse momentum per particle. This particular form has been chosen mainly because it will allow us to carry out analytical calculations. 2 If a patch contains n particles then the probability distribution for its p t is given by For a patch of area A, summing over the Poisson distribution for the number of particles in the patch ( n = νA) gives us the overall probability distribution for the transverse momentum in the patch as where I 1 is the (first order) modified Bessel function of the first kind. The mean and standard deviation of the distribution are given by νAµ and √ 2νAµ. It is convenient to express this as saying that the transverse momentum in a patch of area A is where, in the model discussed here, ρ and σ are given by In the limit in which νA ≫ 1, the distribution in eq. (6) tends to a Gaussian with mean and standard deviation as given in eqs. (7) and (8). This is a consequence of the central-limit theorem.
Traditional approach. Taking the area of each of the transverse regions to be A Trans , the traditional approach will extract the following results for ρ in the transverse Average, Min and Max regions where the Min and Max results have been derived using the Gaussian limit of eq. (6). The only one of the above results that correctly estimates ρ is the average. The Min and Max results tend slowly to the correct answer in the limit in which A Trans is large. In the literature (e.g. [18]), A Trans has usually been taken equal to 2π/3. As we shall see in section 4, a typical value for σ at the LHC ( √ s = 10 TeV) is σ ≃ 0.5ρ − 0.75ρ, which implies that the σ/ √ πA Trans term is about 20 − 30% of ρ.
Area/median-based approach. To help understand the behaviour of the area/medianbased approach, let us replace the jets (which have a range of areas) with uniform rectangular tiles, each of which has a fixed area A tile . It is important to use the full distribution dP/dp t (A tile ) as given by eq. (6) rather than the Gaussian distribution, because a physically interesting domain is that in which νA tile is of order 1. The extracted value ρ ext of the UE p t density in the tiled approximation is given by the median value of p t,tile /A tile across the many tiles in a single event. It can be determined from the solution of the equation This result has been obtained in the limit of there being a large number of tiles, i.e. large A tot , which allows us to approximate the distribution of tile transverse momenta in a specific event with the average probability distribution dP dpt (A tile ) (see also appendix B.1). The integral in eq. (10) is non-trivial to evaluate analytically, however an approximation to the solution for ρ ext that is accurate to a couple of percent and has the correct asymptotic behaviours is given by The result is non-zero only for νA tile > ln 2, which stems from the requirement that tiles with no particles, i.e. the δ(p t ) contribution in eq. (6), should not account for more than half of the total number of tiles. This property of a sudden turn-on, as well as the fact that at large νA tile the offset from the correct ρ goes as 1/(νA tile ), are features that we have found to hold for certain other analytic forms of dP 1 /dp t , notably all those with a structure p m t e −(m+1)pt/µ (for arbitrary positive m). Other characteristics, such as the particular coefficient of the 1/(νA tile ) offset in eq. (12), or the analytic structure close to the turn-on, do depend on the form taken for dP 1 /dp t .
The determination of ρ from the median of jets' p t densities differs from the above "tiled" model in that jets do not have a fixed area. There is no simple way of extending the analytical model so as to account for this, however one can study the impact of the question numerically, by examining toy events in which many soft particles have been generated according to eq. (4), with random positions in the y − φ plane. On each event one runs the jet-based procedure to determine ρ ext and compares it to ρ ext for the tile-based procedure and to the analytic approximation eq. (11). We assume that tiles of A tile are comparable to jets of the same average area, A tile = A jet . One might think that the average jet area should be given by the result for ghost jets in Ref. [12], A ghost−jet ≃ 0.55πR 2 . However, this only holds in the limit of very dense UE; for the typical kinds of configuration that are of interest to us, it will be more appropriate to use where we have defined a constant c J that will reappear in several places below. Given this relation between the typical jet area and radius R, we can deduce the critical radius, R crit , below which ρ is zero, We can also rewrite eq. (12) in terms of R and σ/ρ or R crit : The above results depend on the specific form of toy model that one chooses. To estimate the importance of this model dependence, one can replace dP 1 /dp t as given in eq. (4) with the alternative form 1 P 1 as discussed in more detail in appendix A. The essential relations for this model are σ/ρ = 3/(2ν), R crit = σ ρ (2 ln 2)/(3c J ) ≃ 0.48 σ ρ and ν = 5, exponential p t distr., |y| < π dP 1 / dp t = 4p t e -2p t analytic approx. median tile ρ median jet ρ dP 1 / dp t = e -p t analytic approx. median tile ρ median jet ρ Figure 3: Comparison of jet and tile median-based determinations of ρ for two toy soft-event models (as described in the text), together with the corresponding analytical approximations. The mean number of particles per unit area, ν, was 5 (commensurate with typical expectations for the UE at LHC) and the jets that were used had |y| < π (this choice stems from our use of square tiles). To achieve sub-percent-level agreement between the analytical and jet-based results at larger R values, it was necessary to carry out the determination of ρ ext using p tscheme recombination of 4-momenta in the jets.
These results involve the same analytic structures as for the original form of the toy model, with numerical coefficients that imply slightly smaller corrections for finite R. In fig. 3, the approximate analytical toy model results as a function of R (or the equivalent A tile ) are compared to the average results obtained by applying the "tiled" median approach as well as the jet-area median approach to toy-model configurations. This is done for both variants of the toy model. There is near perfect agreement between the analytical approximations and the average median tile-based results. This is indicative both of the quality of the analytical approximation and of the limited impact of the practical use of finite A tot (here 4π 2 ), as opposed to the large A tot limit that went into the analytical results. The median jet-based results are rather close to the tile-based results for larger R values, 3 though the precise shape in the turn-on region differs a little with respect to the tile-based expectation. Moderate differences exist between the results with the two choices for dP 1 /dp t , and these can be taken as indicative of the magnitude of the model-dependence in the above analysis.
One message to take from fig. 3 is that for an R value that is twice that where the turn-on occurs, ρ ext underestimates ρ by about 10 − 20%. This can be kept in mind as a ballpark value for the accuracies that we will be able to achieve and can be compared to the 20 − 30% effect discussed above for the TransMin ρ determination in the traditional approach.

Purely perturbative events
In our definition of purely perturbative events there is no underlying event and the only regions of the event that have non-zero momentum are those that contain a perturbative emission. Nevertheless UE-determination methods may still give a non-zero result for their estimate of the UE energy density, because of the way they are affected by those extra jets.
We work here with the assumption that we have selected events with at least two hard jets (or with a W or Z boson), and that extra jets may be present at higher perturbative orders. A crude, but illustrative approximation for those higher orders will be obtained as follows. We will take the dominant source of emissions to be radiation from the initial-state partons that enter the hard reaction. 4 Furthermore we will assume that those emissions are soft, distributed uniformly in rapidity and azimuth, and with a p t spectrum given by α s (p t )/p t : where C i is the colour factor associated with the incoming partons (C i = C A ≡ 3 for gluons, . For a dijet event whose two hard jets have transverse momenta p t,hard , we will take eq. (18) to be valid independently of rapidity for Q 0 < p t < Q, where Q ∼ 1 2 p t,hard and Q 0 ∼ 1 GeV. 5 We will also assume that emissions are independent. Thus, the probability distribution for the number of emissions will be a Poisson distribution whose mean is obtained by the integral of eq. (18) over the relevant phase space.
Traditional approach. The average p t densities in the Average and the Max regions will both receive contributions from the emission of one gluon (relative to the Born diagram for the process). In contrast, the Min region only receives a contribution when at least two gluons have been emitted. One can obtain the ρ ext,Av value just by integrating eq. (18) up to Q, which we do in a fixed coupling approximation, α s = α s (Q), since the integral is dominated by values of p t ∼ Q: Here we have neglected the (small) impact of Q 0 and, in the fixed-coupling approximation, the result is complete to all orders in α s . One feature to note about the result is that it scales with Q.
To determine ρ ext,Min to O α 2 s , we assume that the left and right transverse regions each contain one gluon, and that the left-hand gluon (L) is harder than the right-hand one (R); ρ ext,Min is then given by the average transverse momentum in the right-hand region, 4 Ignoring radiation from any outgoing (Born) partons is not too poor an approximation, because a significant part of that radiation will be contained within the corresponding jets. 5 Based on collinear factorisation, one would expect that the upper limit on the pt of emissions to have significant rapidity dependence. For example, if the hard process takes place at central rapidities, then one might write pt p t,hard e −|y| . The rapidity-independent approximation is instead inspired by a high-energy factorisation picture, relevant when pt ≪ √ se −|y| . Studies with Herwig at parton-level (based on collinear factorisation) give a distribution for the upper pt limit on extra jets that is intermediate between these two expressions.
with an additional factor of 2 to account for the case where the right-hand gluon is the harder one. The result is proportional to α 2 s , i.e. suppressed by an extra factor of α s compared to ρ ext,Av , however it is enhanced by a factor of A Trans .
Finally we can estimate ρ ext,Max using the relation ρ ext,Min + ρ ext,Max = 2ρ ext,Av , giving us To appreciate the impact of the various terms, let us take C i ≡ C A = 3, Q = 50 GeV, α s (Q) = 0.13, and A Trans = 2π/3. Then we obtain These numbers scale roughly linearly with Q. The crudeness of our approximations for the perturbative part of the event means that they are not be trusted to better than within a factor of two (worse in the case of ρ ext,Min ). However the rough orders of magnitude are still instructive and highlight the advantage of the Min region.
The above analytic estimates can be verified by using more realistic events from a Monte Carlo generator at parton level. For the case of dijets from Herwig pp collisions at √ s = 10 TeV, with p t of the partons in the hard 2→2 process required to be above 100 GeV (consistent with Q = 50 GeV) and the soft underlying event turned off, one obtains These numbers are very close to those from the simple model of purely perturbative underlying event described above for gluon jets (though Herwig has an admixture of quark jets here).
We have verified that if we double the area over which the measurement is carried out, the ρ MC ext,Min result roughly doubles, as expected from eq. (20). One comment concerning the above results is that in the pure soft UE case it was ρ ext,Av that was the least biased estimate of ρ. Here it is ρ ext,Min that is the least biased by hard perturbative radiation. If one restricts one's attention to ρ ext,Min , then a further property of interest is that in the soft UE case, the bias is reduced by increasing the transverse region's area A Trans , while for hard perturbative contamination increasing A Trans increases the bias. This trade-off between the two issues is characteristic of the difficulty of accurately estimating ρ.
Area/median-based approach Let us suppose that we extract ρ based on jets contained in a region of area A tot . Assuming the typical area for jets as introduced in eq. (13), then the typical number of jets N in the region (including ghost jets) should be given by The exact value of N in each given event will depend on that event's detailed structure (and the exact set of ghosts), but eq. (24) should be adequate for our illustrative discussion here.
Of the N jets, we will assume that n h are "hard" jets, of which n b correspond to the final-state born particles (n b = 2 for dijet events, n b = 0 for Drell-Yan events) and n p stem from perturbative radiation. It is convenient, albeit somewhat simplistic, to model n p as being given by the number of emitted gluons 6 where we consider the number of perturbative emissions between some non-perturbative scale Q 0 and an upper limit Q related to the hard scale of the process (e.g. half the p t of the hardest jet, as before) and we have used a 1-loop running approximation for α s (p t ).
For the median estimator of ρ to be non-zero, at least half the jets should contain perturbative radiation, i.e. n p + n b ≥ N/2. Since the number of primary emissions follows a Poisson distribution, we get the probability of non-zero ρ from the following sum where we have also made the approximation that the sum is dominated by its first term, on the grounds that n p /(N/2 − n b ) ≪ 1. Given P , one can estimate ρ ext by observing that the (N/2) th jet will be the softest of all the perturbative jets, and therefore have p t ∼ Q 0 , giving This is plotted with thick lines in fig. 4 (left) as a function of A tot for R = 0.6, using a 1-loop 5-flavour coupling with Λ = 0.1 GeV 7 and with C i = C A = 3, Q = 50 GeV, Q 0 = 1 GeV and two values for n b , 0 and 2.
To estimate the uncertainties of our analytic formula introduced by approximations (26) and (27) we also plot with thin lines in fig. 4 (left) the result of numerical studies of the same simple set of perturbative emissions (equivalent to the full sum over n). As before, we associate each parton (i.e. gluon in our case) with one jet and assume that this jet has area c J R 2 . In the case n b = 2, we make sure that the two Born particles are always present in the region where the underlying event is measured. One sees in fig. 4 that the exact numerical results (equivalent to taking the full sum in eq. (26)) are moderately higher than the corresponding analytic approximations, as would be expected. In both cases, however, the contribution to ρ ext is negligible except for small values of A tot with n b = 2.
The approximation that we have used for the distribution of the number of jets is rather simplistic. To estimate the order of magnitude of the uncertainties that are present in our toy model of the perturbatively induced "UE", notably those associated with secondary radiation from the partons, we have also applied the area/median method to realistic dijet events from Herwig with soft UE turned off, taken at parton level. To make the comparison, we required that the partons in the hard 2→2 process have p t > 100 GeV, as before (and consistent with our choice Q = 50 GeV). The result is shown in fig. 4 (right) with black and open circles for the cases with and without removal of the two hardest jets from the ensemble used for the median. The A tot -dependence of the Herwig results differs somewhat from the expectations in eqs. (26) and (27). We attribute this to the fact that each parton (be it a Born parton or a primary emission) can itself radiate extra gluons, some of which will lead to additional jets. A simple way of accounting for this is to replace n p in eq. (26) with (1 + X) n p + Xn b , where X is the number of extra jets obtained per parton. 8 Using a modest X = 0.3 (1.3 jets per parton) brings the analytical result into accord with the Herwig results.
In terms of the practical impact of perturbative emission on the extracted ρ values, one sees that for A tot ≥ 4π (the minimal value that we shall use in section 4) it remains a small effect, though the curves also highlight that at small total area A tot , the area/median method can start to become sensitive to perturbative radiation, especially when Born partons are present.
A final point to comment on is the relation between these results and the discussion in [13], where it was argued that ρ ext would be non-zero starting only at O (α N/2−n b s ). Examining eq. (26), one sees that this statement stems from the fact that n p ∼ α s (in a fixed-coupling approximation). However, each power of α s is compensated in part by a power of A ln Q/Q 0 (which is large), and ultimately the small value of P in eq. (26) (and hence ρ ext ) cannot be solely attributed to an α s power-counting argument, but rather involves a more subtle interplay of all the factors in eq. (27).

Two-component events
Realistic events are neither purely perturbative nor consist of pure soft noise. It is instructive to examine what happens if one considers events that have both components together.
Traditional approach. The transverse-momentum density extracted from the average of the two transverse regions is straightforward to calculate in the two component model: it is just given by the sum of the soft and perturbative components, The results for the Min and Max regions are more complex: in the pure soft-component case, it was the soft radiation that determined which of the two regions was Min/Max; analogously, in the perturbative case, it was the perturbative radiation that determined this. When both can be present one has to consider which of the two fixes the Min/Max regions. It is useful to define P to be the fraction of events in which the amount of perturbative radiation in each of the two transverse regions is smaller than the size of soft fluctuations of those regions, p t,L , p t,R ≪ σ √ A Trans . In this set of events, it is the soft component that defines which region is Min/Max, and the bias in the extraction of ρ is just the soft bias, eqs. (9), with no perturbative bias. In the remaining events, in which one or more perturbative emissions are much harder than σ √ A Trans , it is those perturbative emissions that will determine which of the two regions is Min/Max. For these events, there will be no bias in the contribution from the soft component. This implies that the average bias in the Min/Max regions for the soft component will be ∓P · σ/ √ πA Trans . As concerns the perturbative contamination, the average results in eqs. (20) and (21) are already dominated by the set of events in which there is at least one hard emission, so these contributions remain unchanged in the two-component case. The final result for the Min/Max regions is therefore where P is given by .
Note that the choice of lower scale in the logarithm, σ √ A Trans , or Q 0 if that is larger, is only controlled to within a factor of O (1), just as is the choice of Q for the upper limit.
Area/median-based approach. The combination of low-p t and perturbative components is non-trivial also in the case of the area/median-based approach. To treat it analytically, it will be convenient to work at R values that are sufficiently large that the distribution of p t,jet /A jet for the many jets can be considered approximately Gaussian -i.e. we will work away from the turn-on region νA jet = ln 2 that was discussed in section 3.1.
Of the N jets that are used in determining the median, some will contain hard perturbative radiation with transverse momentum significantly above the scale of the fluctuations of the UE. Assuming that there are on average n h hard partons, and that the probability distribution of hard partons in a jet (or tile) is given by a Poisson distribution with mean n h /N , then the average number of jets not contaminated by the hard partons will be given by N e − n h /N ≃ N − n h + n h 2 /(2N ). These uncontaminated jets will have a distribution of values of ρ jet = p t,jet /A jet that is governed just by the soft component and is roughly Gaussian Assuming n h < N ln 2, the median procedure implies finding ρ ext such that N/2 of the N e − n h /N Gaussian-distribution jets have ρ jet < ρ ext , i.e. one must determine the value of ρ ext such that (the unphysical negative lower limit of the integral, an artefact of the Gaussian approximation, doesn't perturb the argument). In the small n h /N limit, this is easily solved and gives where ρ (soft) ext is the result obtained in the pure soft case, eq. (11) of section 3.1. 9 One can then use eqs. (13) and (24) to express A jet and N in terms of R and A tot . Keeping only the first two terms in the R expansion gives Features to note here are that the discrepancy is proportional to σR, with the next correction going as σR 3 . Eq. (25) provides a result for the average number of hard emissions n h , and Q 0 there should be replaced with max(Q 0 , O σ A jet ), or equivalently max(Q 0 , √ c J σR), because perturbative emissions whose transverse momentum is much smaller than the scale of fluctuations of the underlying event will not bias the median. This then gives us We see that the Born particles contribute when A tot is not very large. Perturbative emissions always contribute, essentially because the number of emissions scales with the total area. 9 The additivity of soft and hard results is an approximation, justified only when ρ (soft) ext is close to ρ. An additional point is that when plotting soft+hard results for ρext , we will eliminate the Θ-function in eq. (11) and use the prescription that ρext is well-defined only when the sum of soft and hard contributions gives a positive result.
Substituting physically reasonable numbers into eq. (35), i.e. which corresponds to setting the logarithm L equal to 1, gives which then gives the following numerical result for the bias, Ignoring the n b /A tot terms, for R = 0.6 in gluon-initiated processes (C i = C A ), the bias introduced in ρ is about 0.29σ. For values of σ = 0.5ρ − 0.7ρ, as we will obtain from the Monte Carlo study in section 4, the positive bias due to these perturbative effects is in the same 15−20% ballpark as the negative bias due to finite-particle density that was discussed in section 3.1 for pure soft events. There is an expectation that these two sources of bias should combine linearly, at least when R is sufficiently far above the critical turn-on point that the Gaussian approximation used above is valid. Since they have different R-dependences, ∼ +R and ∼ −1/R 2 respectively, there exists only a limited range of R in which they compensate each other. In this respect the numbers given above tend to confirm the choice R ∼ 0.5 − 0.6 that had originally been recommended based on Monte Carlo studies in [13]. 10 (36) together with the numerical results from the median tile-based and median jet-based approaches. Two sets of results are presented, corresponding to the events with and without Born particles. We see the very good agreement between the analytic and the numerical tile-based approaches and the median jet-based result (except in the threshold region, as in section 3.1). To achieve this it was essential to include the term ∼ R 3 in eq. (34). Fig. 5 suggests that the region of R with the least bias for the determination of ρ is R = 0.4 − 0.6. If one requires that the biases in eqs. (15) and (34) cancel each other, then one finds that R should be chosen proportional to ν 1/6 or equivalently proportional to (σ/ρ) 1/3 ∼ R 1/3 crit . Ignoring the n b /A tot component and the R 3 term in the equations of this section, one finds where the numerical values have been obtained setting L = 1 in eq. (35). The result for R zero-bias can be seen to be consistent with fig. 6, which shows the analytical approximation for ρ ext as a function of R for a broad range of particle densities ν. For a variation in the 10 An alternative method of extraction of ρ comes to mind at this point, fitting a formula motivated by the results of this section and appendix A: with the fit parameters ρ, Rcrit, n, c. One might choose to forgo n, or try a finite number of choices, e.g. n = 1, 2. We have found that such a procedure eliminates a substantial part of the biases in the extraction of ρ for some events, but in others statistical fluctuations lead to poor fits, with results for ρ that are highly skewed towards low values. For this reason, we do not adopt this approach for the current study. without Born median jets median tiles analytic appr. pure soft Figure 5: Underlying event from the two component toy model analysed in the range of rapidity |y| < π and the full span of the azimuthal angle. The results for the average ρ ext (normalised to ρ), extracted in tile and jet median-based approach (see text) are shown for two sets of events. The first set contains soft and perturbative particles while in the second set, two Born particles with p t,hard = 100 GeV are present in addition. For reference, the analytic curve from fig. 3 for the pure soft case is also shown. The mean number of soft particles per unit area, ν, was 5 and their average transverse momentum µ = 0.4 GeV, which corresponds to ρ = 2 GeV. The perturbative emissions are distributed between scales Q 0 = 1 GeV and Q = 50 GeV.
particle density (and ρ) by a factor of 50 (σ/ρ by a factor of 7), the zero-bias R value changes only moderately and in close accord with the expectations of eq. (39). Fig. 6 also illustrates that a fixed R ∼ 0.6 leads to a ρ ext value to within 20% of the correct result for a whole range of ν, with the relative impact of the biases steadily decreasing as the particle density is increased.

Fluctuations in the estimation of ρ
In the simple model studies discussed here, the same ρ value has been used to generate all events. Nevertheless, extracted values of ρ vary from one event to another. This is because any method of measuring UE can use only a limited part of an event (restricted A tot , A Trans , according to y and φ cuts) and works with a finite number of objects (particles, jets). The magnitude of the observed event-to-event fluctuations in this case is an important characteristic of a method, because it sets a lower limit on the uncertainty of the event-toevent ρ measurement. These intrinsic fluctuations also affect the measurement of the true fluctuations and correlations of realistic UE models. Therefore, in order to measure properties of the underlying event that can then be used for efficient subtraction or tuning of simulation programs, one is interested in reducing the fluctuations that come with the method itself.
For the soft underlying event from section 3.1, ignoring the small systematic biases that we found there in the determination of ρ, one can show that the standard deviations S d of The result for the median case is derived in appendix B.1. The formula for the Min/Max regions was obtained in the Gaussian approximation, valid in the limit νA Trans ≫ 1. Substituting realistic values for the areas, i.e. A Trans = 2π/3, A tot = 4π, and the density ν = 5, one arrives at the following numerical estimates The lower expected fluctuations for the area/median based approach are a consequence of the larger area used in the UE determination. The use of a larger area is possible in the first place because the method's dynamical separation of UE limits the need to cut away regions from the y − φ plane to reduce contamination from hard jets. In fig. 7, we show histograms of ρ extracted in the traditional and area/median approaches for the case of purely soft underlying event from section 3.1 (left) and for the combined model described in section 3.3 with the soft and perturbative components and 2 Born particles (right). The corresponding standard deviations of the ρ ext values are given in table 1. For the case of the purely soft UE, we see that these standard deviations follow the pattern from   eq. (41), with the area/median approach performing best. Note however that finite-density (ν) effects do tend to slightly reduce the standard deviation results as compared to expectations, especially in the TransMin case. Adding perturbative and Born particles increases fluctuations in all cases. From the shapes of the curves in fig. 7 (right), it seems that the area/median and TransMin results are both degraded by similar amounts, while the TransAv result suffers significantly more, with a long tail to large ρ ext values. The strong degradation of the TransAv result is an expected consequence of its sensitivity to perturbative radiation, where one observes the dependence on α s /A Trans (to be compared to the bias on the mean which goes as α s ). However, if one examines the results for S d in table 1, one sees that the TransMin standard deviation is also significantly increased by perturbative radiation. The toy-model expectation is S and for the particular parameters used in table 1, the result for S d turns out to be as large as ρ itself. Compared to the O α 2 s A suppression for the bias to the ρ ext,Min result, here the perturbative radiation bias is much stronger, O (α s ). The physical explanation is simple: while it is relatively rare for perturbative radiation to affect the TransMin region (hence the acceptable peak-region of the ρ ext distribution), when it does, the effect on ρ ext is large, contribution significantly to the final result for S d . 11 The area/median approach is much more robust in this respect, because the hard emissions' contribution to the standard deviation does not have significant enhancements compared to the average bias on ρ: as derived in appendix B.2. In particular, the larger numerical coefficient compared to the O (R) term of eq. (37) is compensated by the factor of 1/ √ A tot . This good behaviour is visible in fig. 7, and also in the values of table 1, which are roughly consistent with the above analytical estimate (maybe 20% higher). They also highlight the further improvement to be had with the hard-jet removal procedure discussed at the end of section 2.2 ("all-2 jets" result), which benefits not just S d , but also the peak position and height in fig. 7.
Overall, the results of this section suggest that for any measurement of fluctuations of the UE, it will be preferable to use the area/median method, with hard-jet removal able to provide some extra benefit.

Extraction of σ
The measurement of intra-event fluctuations, σ, has only been discussed so far in the context of the area/median approach. 12 Though we shall not go into full analytical detail, it is easy to convince oneself that many of the considerations that arise in the extraction of ρ apply also when determining σ. In particular, for pure soft events one underestimates σ when R is too small and the presence of perturbative radiation will bias the extracted σ at larger R.
One point to be aware of is that our method for extracting it, cf. section 2.2, has the characteristic that σ/ A jet never exceeds ρ. This translates into an R-dependent upper bound on σ ext /ρ ext , Fig. 8 (left) shows the average σ ext as a function of R in the toy model (ν = 5, with and without the perturbative radiation and the Born jets), normalised to the correct result for the soft component. In the case of just the soft component, one sees a threshold region followed by a slow approach towards the correct value, much as for ρ ext even if the detailed shapes differ (in part owing to the one-sided determination of σ, which causes the residual bias to be proportional to 1/R, rather than 1/R 2 for the bias to ρ). With the inclusion of the perturbative component there is an additional bias, which grows towards large R, again much as happens for ρ ext .
To compare the biases on σ ext and ρ ext , it is convenient to examine the ratio σ ext / ρ ext , fig. 8 (right), normalised to the soft-component result for σ/ρ. The impact of the bound on σ/ρ, eq. (45), is clearly visible up to R ≃ 0.5 (consistent, roughly, with σ/ρ ≃ 0.63). Beyond this point, over a reasonable range of R, σ ext / ρ ext remains compatible with the true value to within roughly 10 − 20%, and, as expected, deviations are larger in the presence of hard radiation. Overall, in the region of R that is suitable for extracting ρ, figure 8 suggests that the extraction of σ should also be quite acceptable. Table 2 summarises the main results of this section for the biases, δρ, and the event-to-event fluctuations, S d , that occur within the toy model in extracting ρ, the transverse momentum flow per unit area, with each of three UE estimation methods: the traditional approach in its TransAv and TransMin variants, whose main parameter is the area A Trans of each of the two transverse regions; and the area/median approach whose parameters are the total area under consideration, A tot and the jet radius R. Results are given both for the biases and fluctuations intrinsic to the soft component and for the additional biases that arise due to the presence of hard radiation in the event, expressed as

Summary of main results
The analytical formulae help illustrate the dependence on the parameters of the measurement methods (the area of the regions used, the jet radius) and the physical scales present in the events (hard scale Q, and the values for ρ and for the level of intra-event UE fluctuations, σ). One sees, for example, how the hard biases and event-to-event fluctuations in the traditional approach are always proportional to Q and to some power of α s , whereas in the area/median approach they are proportional to σ (with modest coefficients and weak additional ln ln Q dependence contained in the parameter L). Table 2 also gives numerical results for the biases and S d values, based on the default set of measurement parameter choices and physical scales that have been used throughout this section. This helps illustrate the expected orders of magnitude of different effects under realistic conditions. The large results for S (hard) d in the traditional method (i.e. unreliable event-by-event extraction of ρ), together with the proportionality to Q of the traditional method's biases and fluctuations and the fact that it offers no easy way of determining σ, lead us to prefer the area/median method for the Monte Carlo UE measurement studies that we will perform below.
Other results of this section that are not summarised in table 2 include: R crit , eq. (14) the R-value below which the area/median approach gives ρ ext = 0; the upper bound on σ ext /ρ ext as a function of R in the median-area method, eq. (45); and R zero-bias , eq. (39), the R value for which the soft and hard biases cancel out in the area/median approach. Finally, fig. 5 shows the characteristic shape of the R-dependence for ρ ext in the median/area method, while fig. 6 helps illustrate how R = 0.6 is a reasonable default choice for a wide range of UE conditions.

Illustration with Monte Carlo events
Given the area/median method to determine ρ and σ on an event-by-event basis, let us now explore what kinds of observables we might construct from them. The choices that we shall make are motivated by considerations of how the UE affects jet measurements at hadron colliders. For example the UE leads to an average shift in jet energy; it's important that quantity method result numerical value equation(s) TransAv Q is the hard scale, Q 0 the IR cutoff on perturbative emissions, 2b 0 ln α s (max(Q 0 , σ √ A Trans ))/α s (Q) , c J = 2.04, C i is the colour factor of the incoming partons. Numerical values are given for ρ = 2 GeV, σ/ρ = 0.63, Q = 50 GeV, Q 0 = 1 GeV, C i = C A = 3, n b = 0, R = 0.6, A tot = 4π and A Trans = 2π/3 (the area of a single transverse region), corresponding to L ≃ 1 and P ≃ 0.35. The result marked with a * is specific to the form of the soft toy model discussed in section 3.1.
one knows that shift as a function of rapidity, and hence one should determine ρ(y) . The UE also affects jet energy resolution, through a term of the form σ A jet . Thus we will also want to look at σ as a function of rapidity. A second way in which the UE affects jet energy resolution is that ρ itself is different event by event, so one might therefore examine its event-by-event distribution and its standard deviation S d . The different ways in which UE affects jets was summarised in [12] with the following equation for the variance of the change in jet p t due to the underlying event (neglecting back-reaction): where A JA,R is the average jet area and Σ 2 JA,R ∼ R 4 is the variance of the jet area. One sees that each of the terms involves a different characteristic of the UE: ρ , σ and S 2 d . Measurements of UE characteristics, as well as being of direct relevance to jet measurements, also have the power to constrain UE models. This has been the motivation for most UE studies to date and we believe that the range of UE characteristics discussed above would complement existing types of measurement and so provide an additional powerful set of constraints. Furthermore, as well as examining "local" quantities, such as ρ(y) , one can also, for example, ask the question of whether there are long-range correlations between the magnitude of the UE in different parts of an event. In (prevalent) UE models that involve multiple independent 2 → 2 scatterings, one might expect these correlations to be modest, whereas in a BFKL-inspired model (as might derive from work like Ref. [16]), where one or more gluon ladders stretch across a whole event, one might expect them to be larger.
So far in this article we have treated the traditional and area/median based UE measurement approaches on a similar footing. The results of section 3 suggest that both provide reasonably adequate information about average characteristics of the underlying event: both methods introduced biases and different sources of biases cancelled partially, limiting their overall impact (though in the area/median approach the biases had much weaker dependence on the hard scale of the event). In contrast, as we saw in section 3.4 and its fig. 7, in the traditional approach the fluctuations due to hard perturbative emission were likely to dominate over fluctuations in the soft component of the UE, therefore for observables that are sensitive to fluctuations, we believe that the area/median approach is more robust. For this reason, in what follows, we shall concentrate on this second approach.
This section will be structured as follows: after outlining the set of Monte Carlo simulations that we use (section 4.1) and the event selection cuts (section 4.1), we examine to what extent the toy model is qualitatively similar to the realistic simulations (section 4.3). The motivation is that the toy model guides our intuition about the measurement procedure, and it is important to establish that this intuition is well founded. Having done so we then study (section 4.4) a selection of observables that are relevant for eq. (47), including their rapidity dependence and also their degree of correlation across different parts of a same event. This set of observables (together with some of those in section 4.3) would, we believe, be interesting to examine experimentally, both in terms of the direct information it provides for understanding the impact of UE on jets (and isolation, etc.) and in terms of its ability to constrain models.

Monte Carlo models used
We shall examine a series of Monte Carlo UE models: the UE that comes with the "old" virtuality-ordered shower in Pythia 6.4 [21], in the DW and DWT tunes by R. Field [18], identical at Tevatron energies, but with different energy-dependences; the UE that comes with the "new" transverse-momentum-ordered shower in Pythia 6.4, in the S0A tune [5,34,35] by P. Skands; Herwig 6.5's [36,37] default "soft" UE, which fails to reproduce various aspects of Tevatron data, but instructive for comparisons between different types of models; Herwig 6.5 with Jimmy 4.3 [4] in an ATLAS tune 13 by Moraes [18]. All models are based on multiple interactions except for Herwig's soft UE.
By default we use a pp centre of mass collision energy of √ s = 10 TeV, though we will also consider the energy dependence of some of the observables we study.
A comment is due on the fact that we will carry out our investigations at particle (hadron) level. Experiments may carry out measurements on tracks only (well measured, though with some low p t threshold), on calorimeter information (subject to noise and noise-suppression thresholds) or on some combination of the the two. It is beyond the scope of this article to estimate the potentially substantial impact of detector effects on the results presented here. Nevertheless, we believe that the differences that we will see between various event-generator tunes should persist even after detector effects.

Event selection
We consider dijet events, where the leading jet, reconstructed with the anti-k t algorithm [31] with R = 0.6 has p t > 100 GeV, the next hardest jet has p t > 80 GeV and both jets are required to lie in the rapidity window |y| < 4. Note that since the cross section for jets falls steeply, a cut on jet p t introduces a "trigger-bias", i.e. favours events where the UE is slightly larger than average. 14 The above cut is a default used for the study presented in this section. In some places we employ tighter cuts on the rapidities of the two hardest jets in order to study the impact of these jets on UE determination. Whenever this is the case, it will be indicated explicitly. The choice of the anti-k t algorithm for event-selection here is motivated by the fact that it will be the first to be supported by ATLAS and CMS [39] and therefore will have been used by the experiments for their initial event selection.
For the determination of the properties of the UE within the jet-area/median approach, the anti-k t algorithm is not suitable, as discussed in section 2.2. We will therefore use jets from the C/A algorithm with R = 0.6 for this task. All jet finding is performed with FastJet 2.4.1 [40,41], and we use the ActiveAreaExplicitGhosts option to calculate areas, because this ensures the safest treatment of pure ghost jets. The ghost area that is used is 0.01 and ghosts are placed up to |y| < 9 (we could, however, have used a smaller upper limit for most of the plots). The other parameters are left at their default values.

Comparisons of characteristics of MC and toy model
Given the importance of the toy model in guiding our understanding of the measurement of the underlying event, let us start by examining whether realistic Monte Carlo events bear any similarity to toy-model events. One way of doing this is to examine the R-dependence 13 The non-default parameter setting are: PRSOF=0, JMRAD(73)=1. 8 [18]). 14 To investigate the impact of the trigger bias, we used the following procedure: we first measured the underlying event density ρ in a rapidity window |δy| < 1 around each hard jet, and then placed our hard jet cuts on pt,jet − ρAjet rather than on pt,jet. The result for ρ in events that passed these cuts on "subtracted" jets came out about 10% lower than for events where we cut on unsubtracted jets. Since this is not too large an effect, and for reasons of simplicity, we have however chosen not to apply this procedure to our analysis as a whole. 1/n dn/(dp t /A) [GeV Pythia DWT pp, √ s = 10 TeV a-k t +C/A, |y| < 5 Figure 9: The left-hand plot shows ρ(R) extracted with C/A jets in the area/median approach for 3 representative Pythia (DWT) events, which passed the selection cuts of section 4.2.
The right-hand plot shows the corresponding histograms of p t /A for the same 3 events.
of the extracted ρ, which, in the toy model, had a characteristic shape, fig. 5. Fig. 9 (left) shows ρ(R) for three representative Pythia events (DWT tune). In each case one observes the turn-on at some threshold R value, followed by a roughly linear slope at larger R, precisely as expected. There is substantial variation in the curves from one event to the next, and one can trace this back to their distributions of p t /A shown in fig. 9 (right): the blue (dashed) line which has small ρ and little R dependence corresponds to an event in which there is no activity at high p t /A values. The green (dotted) curve, which has large ρ and strong R-dependence, has a correspondingly broad distribution of p t /A values, with much activity at intermediate p t values. In this event, the toy-model picture of a clean separation between soft and hard physics is somewhat challenged, though the general pattern of ρ(R) having a turn-on followed by linear R-dependence still holds. If we average ρ(R) over many events we get fig. 10 (left), which shows results for several generators/tunes. One observation is that the Pythia tunes have a larger slope: based on the toy-model calculations, eq. (34), this can be interpreted as meaning that they have a larger value for the quantity σn h , where n h is the number of hard jets. The Herwig default curve, with no MPI, has the smallest slope (and ρ value). The average over many events for σ(R) (normalised to ρ(R) ) is shown in fig. 10 (right). An interesting feature is the linear rise for R = 0.2 − 0.4 (and up to R = 0.6 for the Pythia tunes). In this region, one saturates the bound eq. (45), which implies that R is too small for a proper measurement of σ. For R 0.6 there is a shallower rise, which we interpret as being due to the presence of hard jets, as is the case in fig. 8 (though there is less curvature in the MC events than in the toy model).
To get an idea of the event-to-event variations of ρ(R) we use the following procedure: given N events, for a given R, we sort the events according to ρ(R). We then define the 10 th percentile result for ρ(R) to be the value of ρ(R) for event N/10, and similarly for other percentiles. Fig. 11 shows the 10 th , 20 th , etc. percentile results for ρ(R), as a function of R for our 4 generator/tune combinations (together with the average, for comparison). 15 One observes the sharp turn-on as a function of R (washed out in the ρ(R) ). The smaller the turn-on point, R crit , the larger than value of ρ (and the larger the slope). The spread of events is noticeably large, both in the values of ρ and for R crit , especially, for the latter, in the context of Herwig+Jimmy. Though it is not our intention to highlight figures 10 and 11 as main results of this paper, we do note that it would be possible to make corresponding experimental measurements, and use them as input to MC tunes.
In the remaining parts of this section we shall concentrate on results extracted with R = 0.6, which, as anticipated in the toy-model section, seems to offer a reasonable compromise between being sufficiently large as to be well beyond the turn-on in most events, while not being too severely affected by the rise at large R that is induced by semi-hard radiation in the event.

Mean energy flow
Let us start by examining ρ and its rapidity dependence, fig. 12, for Tevatron and LHC energies. 16 The results for ρ in the central rapidity bin for Tevatron (left) should be strongly constrained by the standard UE measurements at Tevatron, and this is reflected in the small difference between S0A and DW, though there is a somewhat larger difference with the Her-wig+Jimmy tune. The rapidity dependence is quite strong, with stronger suppression at 15 A subtle point in the production of fig. 11 is that the event that provides the 10 th percentile for (say) R = 0.6 is usually not the same event that provides it for R = 0.8. Thus the curves in fig. 11 are not the curves that would be obtained for individual events (these are far less smooth, cf. fig. 9), but can be thought of as some idealisation of these curves in a world free of fluctuations. 16 To simplify the comparisons, we use the cuts of section 4.2 in both cases, though they involve a rapidity range that extends beyond the Tevatron's coverage. forward rapidities for S0A and Herwig+Jimmy than for the DW tune. One should remember in examining the rapidity dependence that there are essentially no experimental constraints on the level of the underlying event at forward rapidities -it is therefore a model-dependent extrapolation. At LHC energies we see, fig. 12 (right), that differences appear between models also at central rapidities, reflecting an uncertainty in the extrapolation in energy. The DWT tune's energy-dependence is disfavoured based on RHIC [42] and lower-energy Tevatron data [25], but we include it to give an idea of the magnitude of possible differences. One question that arises in the measurement of ρ is the possible bias from hard jets and the relation of this bias with the rapidity bin size or equivalently the area A tot in which one measures ρ. In this context, recall that eq. (37) implies a bias from the Born jets that scales as n b /A tot . In fig. 12 the hard jets were free to lie anywhere within |y| < 4 and were included in the sample of jets used to obtain the median. In fig. 13 (left) the red (dark) dashed curve is the same as the DW result in fig. 12 (right). The red (dark) solid curve shows the impact of removing the 2 hardest jets from the median sample. One sees that this has a rather small effect. Next we examine an event sample in which the 2 hardest jets are in the central rapidity bin, so that that bin receives all the bias from the Born jets. If those jets are included in the set used to calculate the median (dotted green histogram), then the impact on ρ(y) in that bin becomes noticeable, O (30%). Removing them from the set for the median (solid green curve) brings us almost back to the basic "all-2" result for the full dijet sample. That the  result is not fully identical is a consequence of the fact that when the two hard jets are central there is an increased probability that the 3rd hardest jet will also be central, thus biasing very slightly the central-y bin. The right-hand plot of fig. 13 shows what happens if we reduce the rapidity bin size, causing ρ to be measured in regions of area A tot = 2π rather than A tot = 4π. Since the impact of the Born particles is inversely proportional to A tot , requiring the Born particles to be in the central bin has a noticeably larger effect for the smaller rapidity bin size. 17 Discarding the two hard jets brings us back to a result that is roughly in accord with that for the larger bin size.
The conclusion from fig. 13 is that if one's event selection does not constrain the hard jets to be in the same bin as that used for measuring ρ and if the bin area is sufficiently large, A tot 12, then biases from the hard jets are quite small. In what follows, we will normally use A tot = 4π and leave in the hard jets, in order to keep the analysis as simple as possible.

Fluctuations
In fig. 14 we examine fluctuations, now only for LHC. The left-hand plot shows the size of intra-event fluctuations, through the ratio of σ(y) to ρ(y) , while the right hand plot shows inter-event fluctuations, through S d (y)/ ρ(y) . Fluctuations have not been as directly tuned as energy flow. Despite this the intra-event fluctuations are very similar across all the Pythia tunes and almost independent of rapidity when normalised to ρ(y) . Herwig+Jimmy's intraevent fluctuations are somewhat smaller, but do have rapidity dependence.
Concerning inter-event fluctuations, fig. 14 (right), the two virtuality ordered (DW/DWT) Pythia models are again flat, whereas the p t -ordered shower has increasing fluctuations at forward rapidities. Herwig is intermediate between the two sets of Pythia results.
One observation is that whereas DW/DWT have almost identical intra and inter-event fluctuations, Herwig's intra-event fluctuations are nearly 40% smaller than the inter-event  Figure 15: Correlation of ρ(y 2 ) with ρ(y 1 ), shown as a function of y 2 for y 1 in a given rapidity bin: (a) −1 < y < 1, (b) 0 < y < 1 and (c) 3 < y 1 < 5. In plots (a) and (c) ρ has been determined in bins of size δy = 2, while in (b) it has been determined in bins of size δy = 1.
fluctuations. This is reflected also when we examine correlations between ρ in different parts of the event, as shown in fig. 15. The correlations are noticeably larger for Herwig+Jimmy than they are for all the Pythia tunes. In determining the correlations it was important that we used sufficiently large rapidity bins. Comparing the upper-left plot from fig. 15 (δy = 2) and the upper-right plot (δy = 1) one sees that the smaller rapidity bins lead to noticeably smaller measured correlations. We interpret this as follows: in small rapidity bins, the "statistics" of jets for measuring the ρ value are more limited. This increases the error on the determination of ρ, thus reducing the maximum amount of correlation that can be observed between different bins. The final quantity that we examine is the event-by-event distribution of ρ, fig. 16, for a central rapidity bin (left) and a forward rapidity bin (right). Perhaps the most striking characteristic of these plots is the very broad nature of these distributions, which are far   Figure 16: The event-by-event distribution of ρ for a selection of generators, in a central rapidity bin, |y| < 1 (left-hand plot) and a forward rapidity bin, 3 < y < 5 (right).
from being Gaussian distributions of width S d centred on ρ . The right-hand plot also has a significant bin at ρ = 0 -i.e. there is a substantial number of events for which at least half of the jets at forward rapidities are pure ghost jets.

Energy dependence of results
Fig . 17 summarises the energy-dependence of the average energy flow ρ (left) and of the intra-event fluctuations σ normalised to ρ (right) in the central rapidity bin |y| < 1. The features of note are that for ρ , the Herwig+Jimmy Atlas tune has a significantly steeper energy dependence than the DW and S0A Pythia tunes, somewhat more like DWT. For σ / ρ , the Pythia tunes have almost no energy-dependence, whereas Herwig has substantial energydependence. We have not shown the energy dependence of S d / ρ because it is essentially independent of √ s for all generators.

Conclusions
The distinction between underlying event and the "hard" part of hadron-collider events is illdefined in QCD. Nevertheless, physically, one may picture the UE as being the component of a hadron-hadron scattering that fills the event fairly uniformly with low-p t radiation. A goal of this article was to investigate how different UE-measurement strategies fare in separating such a low-p t component from the hard part of the event.
To do so, we developed a simple toy-model for events, with two distinct components, one soft, corresponding to the UE, the other perturbative and hard. Within this model it is quite straightforward to establish to what extent a given UE measurement strategy correctly extracts just the soft component.
The two UE measurement strategies that we investigated are the "traditional" approach, measuring radiation in regions transverse to a leading jet, and the "jet-area/median" approach. Both involve strategies to help separate out the soft and hard components: the use of   the Av/Min/Max transverse regions in the traditional approach, the use of median activity rather than average in the jet-area/median approach. One result from the toy model is a quantification of how those strategies fare in the extraction of its soft component. A second result is a determination of the nature of the residual effects due to perturbative (hard) radiation. These two results could be expressed analytically in terms of the parameters of the measurement procedures (transverse-region area, jet radius, providing useful guidance in choosing them) and of characteristics of the hard scattering (notably the value of the hard scale and the properties of the soft component). We also examined the question of event-to-event fluctuations in the extraction of the characteristics of the soft component. Our toy-model results are summarised in section 3.6. Practically, one conclusion from this work is that for determinations of averaged quantities, for example the mean transverse-momentum density per unit area, ρ , both the TransMin and the area/median measurement methods give a fair determination of the soft component, as long as the momentum transfer of the hard scattering is not too large ( 100 GeV for the LHC; for higher momentum transfers, the TransMin method is affected by a rapidly increasing hard contribution). In particular, for the parameter choices used or advocated in the literature, the two kinds of bias seen in the toy model, mismeasurement of the soft component and contamination from the hard component, tend to partially cancel each other, giving a limited overall bias, of order 20%. In contrast, for event-by-event measurements and determinations of fluctuations of the soft component, the traditional approach is significantly affected by rare "outliers", cf. table 1.
For this reason, in the full Monte Carlo studies of section 4 we concentrated on the area/median approach. The results included a validation of the main qualitative prediction from the toy model, namely the structure of the R-dependence of the extracted ρ, section 4.3. Section 4.4 showed a range of possible observables whose measure we advocate at LHC: rapidity dependence of the UE, nature of the event-to-event fluctuations, and intra-event fluctuations and correlations. Though existing measurements may indirectly constrain some of these features of the UE, we believe that they are of sufficient practical importance that they deserve dedicated measurements, especially as they differ noticeably between various Monte Carlo models.

A.2 Variant of toy model
We can also consider a model with 1 P 1 dP 1 dp t = 4p t µ 2 e −2pt/µ .
At large ν one finds, again numerically, The two above equations both imply that ρ ext turns on and approaches its asymptotic value somewhat faster than in the model with dP 1 dpt ∝ e −pt/µ /µ. A reasonable analytic approximation for ρ ext over the whole domain is ρ ext ≃ ρ νA − ln 2 νA − ln 2 + 2 3 Θ(νA − ln 2) , which has the correct large-ν behaviour, and very nearly the correct coefficient for the √ νA − ln 2 turn-on.
B Fluctuations in area/median extraction of ρ

B.1 Pure soft case
To determine the event-to-event fluctuations S d in the area/median extraction of ρ when the intrinsic event-to-event fluctuations are zero, it is convenient to work in the limit νA jet ≪ 1 such that the probability distribution of ρ jet is close to a Gaussian, 18 The corresponding cumulative probability distribution for the jets is given by If the number N of jets is odd, N = 2m + 1, then the probability distribution of the median δρ is obtained from the product of the probability of having one jet with ρ jet = ρ − δρ, m jets with ρ smaller than this, and m jets with ρ larger than this: Making use of the expansion of the error function working in the large m and small δρ/σ limit, and making use also of Stirling's formula m! ≃ √ 2πm(m/e) m , one can approximate eq. (63) as where in the last step we have replaced m ≃ N/2. Using the relation N A jet = A tot , we finally obtain the following result for the standard deviation of the extracted median ρ values, This is about 25% larger than the standard deviation that would be obtained for ρ extracted as an average of the ρ jet values over the same total area. This moderate enhancement of fluctuations in the pure soft case is part of the price that one pays in exchange for the median's greater resilience to hard contamination.

B.2 Hard contamination
The result eq. (34) for the average discrepancy in ρ ext due to hard contamination can be obtained in an alternative manner, which will be more useful for estimating fluctuations. In this approach we imagine some distribution of soft jets, and then add in the hard partons. Some of the hard partons will enter jets whose ρ jet is already above the median value for ρ. These hard partons will not affected ρ ext . The remaining hard partons (a number k) will enter jets that were below the median. These jets will acquire much larger transverse momenta, taking them well above the median. Thus it becomes necessary to recalculate the median, which will be shifted by k soft jets' worth. From eq. (61), and working in the large-N limit (throughout this section), this translates to an average shift in ρ ext of δρ = kδ 1 , where δ 1 is the 1-jet shift. Substituting k = n h /2 gives in accord with the first order term in n h /N in eq. (33). We will consider two main sources of fluctuations in this result. Let us first imagine that k = 1. The median will shift up by one jet, and the distribution of δρ will be simply be given by the distribution of the difference in ρ jet between two neighbouring jets in the sorted sequence of jets (at position in the sequence that is near the median). That distribution is an exponential distribution with mean δ 1 , The distribution of the shift for k jets is with the standard deviation √ kδ 1 . The second source of fluctuations comes from the fact that k is not constant but rather has a Poisson distribution with mean n h /2 and standard deviation n h /2. 19 Combining this with the fluctuation on ρ for fixed k, leads to an overall standard deviation of In terms of R and A tot , this becomes where L ≃ 1 was defined in eq. (35).