text:
# The Underlying Event in Hard Scattering Processes

## I Introduction

FIG. 1 illustrates the way QCD Monte-Carlo models simulate a proton-antiproton collision in which a ”hard” $`2`$-to-$`2`$ parton scattering with transverse momentum, $`p_T(\mathrm{hard})`$, has occurred. The resulting event contains particles that originate from the two outgoing partons (plus initial and final-state radiation) and particles that come from the breakup of the proton and antiproton (i.e., “beam-beam remnants”). The “underlying event” is everything except the two outgoing hard scattered “jets” and receives contributions from the “beam-beam remnants” plus initial and final-state radiation. The “hard scattering” component consists of the outgoing two “jets” plus initial and final-state radiation.

The “beam-beam remnants” are what is left over after a parton is knocked out of each of the initial two beam hadrons. It is the reason hadron-hadron collisions are more “messy” than electron-positron annihilations and no one really knows how it should be modeled. For the QCD Monte-Carlo models the “beam-beam remnants” are an important component of the “underlying event”. Also, it is possible that multiple parton scattering contributes to the “underlying event”. FIG. 2 shows the way PYTHIA pythia  models the “underlying event” in proton-antiproton collision by including multiple parton interactions. In addition to the hard $`2`$-to-$`2`$ parton-parton scattering and the “beam-beam remnants”, sometimes there is a second “semi-hard” $`2`$-to-$`2`$ parton-parton scattering that contributes particles to the “underlying event”.

Of course, from a certain point of view there is no such thing as an “underlying event” in a proton-antiproton collision. There is only an “event” and one cannot say where a given particle in the event originated. On the other hand, hard scattering collider “jet” events have a distinct topology. On the average, the outgoing hadrons “remember” the underlying the $`2`$-to-$`2`$ hard scattering subprocess. An average hard scattering event consists of a collection (or burst) of hadrons traveling roughly in the direction of the initial beam particles and two collections of hadrons (i.e., “jets”) with large transverse momentum. The two large transverse momentum “jets” are roughly back to back in azimuthal angle. One can use the topological structure of hadron-hadron collisions to study the “underlying event” field\_prd ; field\_dpf ; tano . The ultimate goal is to understand the physics of the “underlying event”, but since it is very complicated and involves both non-perturbative as well as perturbative QCD it seems unlikely that this will happen soon. In the mean time, we would like to tune the QCD Monte-Carlo models to do a better job fitting the “underlying event”. The “underlying event” is an unavoidable background to most collider observables. To find “new” physics at a collider it is crucial to have Monte-Carlo models that simulate accurately “ordinary” hard-scattering collider events. In this talk I will compare collider observables that are sensitive to the “underlying event” with the QCD Monte-Carlo model predictions of PYTHIA 6.115 pythia , HERWIG 5.9 herwig , and ISAJET 7.32 isajet  and discuss the tuning of PYTHIA.

## II The “Transverse” Region

In a proton-antiproton collision large transverse momentum outgoing partons manifest themselves, in the laboratory, as a clusters of particles (both charged and neutral) traveling in roughly the same direction. These clusters are referred to as “jets”. In this analysis we examine only the charged particle component of “jets”. Our philosophy in comparing the QCD Monte-Carlo models with data is to select a region where the data is very “clean” so that “what you see is what you get” (almost). Hence, we consider only charged particles measured by the CDF central tracking chamber (CTC) in the region $`p_T>0.5\mathrm{GeV}/\mathrm{c}`$ and $`|\eta |<1`$, where the track finding efficiency is high and uniform (estimated to be $`92\%`$ efficient) and we restrict ourselves to charged particle jets with transverse momentum less than $`50\mathrm{GeV}/\mathrm{c}`$. The data presented here are uncorrected. Instead the theoretical Monte-Carlo models are corrected for the track finding efficiency by removing, on the average, $`8\%`$ of the charged particles. The theory curves have an error (statistical plus systematic) of about $`5\%`$. Thus, to within $`10\%`$ “what you see is what you get”.

Charged particle “jets” are defined as clusters of charged particles ($`p_T>0.5\mathrm{GeV}/\mathrm{c}`$, $`|\eta |<1`$) in “circular regions” of $`\eta `$-$`\varphi `$ space with radius $`R=0.7`$. Every charged particle in the event is assigned to a “jet”, with the possibility that some jets might consist of just one charged particle. The transverse momentum of a charged jet, $`P_T(\mathrm{chgjet})`$, is the scalar $`p_T`$ sum of the particles in the jet. We use the direction of the leading charged particle jet to define correlations in azimuthal angle, $`\mathrm{\Delta }\varphi `$. The angle $`\mathrm{\Delta }\varphi =\varphi \varphi _{\mathrm{chgjet}\mathrm{\#}1}`$ is the relative azimuthal angle between a charged particle and the direction of chgjet#1. The“toward” region is defined by $`|\mathrm{\Delta }\varphi |<60^{}`$ and $`|\eta |<1`$, while the “away” region is $`|\mathrm{\Delta }\varphi |>120^{}`$ and $`|\eta |<1`$. The “transverse” region is defined by $`60^{}<|\mathrm{\Delta }\varphi |<120^{}`$ and $`|\eta |<1`$. The three regions “toward”, “transverse”, and “away” are shown in FIG. 3. Each region has an area in $`\eta `$-$`\varphi `$ space of $`4\pi /3`$. As illustrated in FIG. 4, the “toward” region contains the leading charged particle jet, while the “away” region, on the average, contains the “away-side” jet. The “transverse” region is perpendicular to the plane of the hard $`2`$-to-$`2`$ scattering and is therefore very sensitive to the “underlying event”.

FIG. 5 and FIG. 6 compare the “transverse” $`N_{\mathrm{chg}}`$ and the “transverse” $`P_T\mathrm{sum}`$, respectively, with the QCD Monte-Carlo predictions of HERWIG, ISAJET, and PYTHIA 6.115 with their default parameters and $`p_T(\mathrm{hard})>3\mathrm{GeV}/\mathrm{c}`$. The solid points are Min-Bias data and the open points are the JET20 data. The JET20 data connect smoothly to the Min-Bias data and allow us to study observables over the range $`0.5<P_T(\mathrm{chgjet}\mathrm{\#}1)<50\mathrm{GeV}/\mathrm{c}`$. The average number of charged particles in the “transverse” region doubles in going from $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`=1.5\mathrm{GeV}/\mathrm{c}`$ to $`2.5\mathrm{GeV}/\mathrm{c}`$ and then forms an approximately constant “plateau” for $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`>5\mathrm{GeV}/\mathrm{c}`$. If we suppose that the “underlying event” is uniform in azimuthal angle $`\varphi `$ and pseudo-rapidity $`\eta `$, the observed $`2.3`$ charged particles at $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`=20\mathrm{GeV}/\mathrm{c}`$ translates to $`3.8`$ charged particles per unit pseudo-rapidity with $`p_T>0.5\mathrm{GeV}/\mathrm{c}`$ (multiply by $`3`$ to get $`360^{}`$, divide by $`2`$ for the two units of pseudo-rapidity, multiply by $`1.09`$ to correct for the track finding efficiency). We know that if we include all $`p_T>50`$ MeV/c there are, on the average, about four charged particles per unit rapidity in a “soft” proton-antiproton collision at 1.8 TeV CDF3 . The data in FIG. 5 imply that in the “underlying event” of a hard scattering there are, on the average, about $`3.8`$ charged particles per unit rapidity with $`p_T>0.5\mathrm{GeV}/\mathrm{c}`$! Assuming a charged particle $`p_T`$ distribution of $`e^{2p_T}`$ (see FIG. 17) implies that there are roughly $`10`$ charged particles per unit pseudo-rapidity with $`p_T>0`$ in the “underlying event” (factor of e). Since we examine only those charge particles with $`p_T>0.5\mathrm{GeV}/\mathrm{c}`$, we cannot accurately extrapolate to low $`p_T`$, however, it is clear that the “underlying event” has a charge particle density that is at least a factor of two larger than the four charged particles per unit rapidity seen in “soft” proton-antiproton collisions at this energy.

The Min-Bias data were collected with a very “loose” trigger. The CDF Min-Bias trigger requirement removes elastic scattering and most of the single and double diffraction events, but keeps essentially all the “hard-scattering” events. In comparing with the QCD Monte-Carlo models we do require that the models satisfy the CDF Min-Bias trigger, however, for $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`>5\mathrm{GeV}/\mathrm{c}`$ essentially all the generated events satisfy the Min-Bias trigger (i.e., the Min-Bias trigger is unbiased for large $`p_T`$ “jets”). If we had enough Min-Bias events we would not need the JET20 data, but because of the fast fall-off of the cross section we run out of statistics at around $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`=20\mathrm{GeV}/\mathrm{c}`$ (that is why the Min-Bias data errors get large at around $`20\mathrm{GeV}/\mathrm{c}`$). The JET20 data were collected by requiring at least $`20\mathrm{GeV}`$ of energy (charged plus neutral) in a cluster of calorimeter cells. We do not use the calorimeter information, but instead look only at the charged particles measured in the CTC in the same way we do for the Min-Bias data. The JET20 data is, of course, biased for low $`p_T`$ jets and we do not show the JET20 data below $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$ around $`20\mathrm{GeV}/\mathrm{c}`$. At large $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$ values the JET20 data becomes unbiased and, in fact, we know this occurs at around $`20\mathrm{GeV}/\mathrm{c}`$ because it is here that it agrees with the (unbiased) Min-Bias data.

We expect the “transverse” region to be composed predominately of particles that arise from the break-up of the beam and target and from initial and final-state radiation. This is clearly the case for the QCD Monte-Carlo models as can be seen in FIGS. 7-9, where the predictions for the “transverse” region are divided into two categories: charged particles that arise from the break-up of the beam and target (beam-beam remnants), and charged particles that result from the outgoing jets plus initial and final-state radiation (hard scattering component). For PYTHIA the “beam-beam remnant” contribution includes contributions from multiple parton scattering. It is interesting to see that in the QCD Monte-Carlo models it is the “beam-beam remnants” that are producing the approximately constant “plateau”. The contributions from initial-state and final-state radiation increase as $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$ increases. In fact, for ISAJET it is the sharp rise in the initial-state radiation component that is causing the disagreement with the data for $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`>20\mathrm{GeV}/\mathrm{c}`$. The hard scattering component of HERWIG and PYTHIA does not rise nearly as fast as the hard scattering component of ISAJET.

There are two reasons why the hard scattering component of ISAJET is different from HERWIG and PYTHIA. The first is due to different fragmentation schemes. ISAJET uses independent fragmentation, which produces too many soft hadrons when partons begin to overlap. The second difference arises from the way the QCD Monte-Carlo produce “parton showers”. ISAJET uses a leading-log picture in which the partons within the shower are ordered according to their invariant mass. Kinematics requires that the invariant mass of daughter partons be less than the invariant mass of the parent. HERWIG and PYTHIA modify the leading-log picture to include “color coherence effects” which leads to “angle ordering” within the parton shower. Angle ordering produces less high $`p_T`$ radiation within a parton shower which is what is seen in FIG. 10.

Of course, the origin of an outgoing particle (“beam-beam remnant” or “hard-scattering”) is not an experimental observable. Experimentally one cannot say where a given particle comes from. However, we do know the origins of particles generated by the QCD Monte-Carlo models and FIGS. 7-9 show the composition of the “transverse” region as predicted by ISAJET, HERWIG, and PYTHIA.

## III Maximum and Minimum “Transverse” Regions

We now break up the “transverse” region into two pieces. As illustrated in FIG. 11, on an event by event basis, we define “transMAX” (“transMIN”) to be the maximum (minimum) of the two “transverse” pieces, $`60^{}<\mathrm{\Delta }\varphi <120^{}`$, $`|\eta |<1`$, and $`60^{}<\mathrm{\Delta }\varphi <120^{}`$, $`|\eta |<1`$. Each has an area in $`\eta `$-$`\varphi `$ space of $`2\pi /3`$ and what we previously referred to as the “transverse” region is the sum of “transMAX” and “transMIN”. One expects that “transMAX” will pick up more of the initial and final state radiation while “transMIN” should be more sensitive to the “beam-beam remnant” component of the “underlying event”. Furthermore, one expects that the “beam-beam remnant” component will nearly cancel in the difference, “transMAX” minus “transMIN”. If this is true then the difference, “transMAX” minus “transMIN”, would be more sensitive to the “hard scattering” component (i.e., initial and final-state radiation). I believe that this idea was first proposed by Bryan Webber and then implemented in a paper by Jon Pumplin pumplin .

FIG. 12 and FIG. 13 show the data on the $`N_{\mathrm{chg}}`$ and $`P_T\mathrm{sum}`$, respectively, for the“transMAX” and “transMIN” region as a function of the $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$ compared with QCD Monte-Carlo predictions of HERWIG, ISAJET, and PYTHIA with their default parameters and $`p_T(\mathrm{hard})>3\mathrm{GeV}/\mathrm{c}`$. FIG. 14, FIG. 15, and FIG. 16 show the data on $`N_{\mathrm{chg}}`$ for “transMAX”, “transMIN”, and the difference “transMAX” minus “transMIN”, respectively, compared with QCD Monte-Carlo predictions of HERWIG. The predictions of HERWIG are divided into two categories: charged particles that arise from the break-up of the beam and target (beam-beam remnants), and charged particles that result from the outgoing jets plus initial and final-state radiation (hard scattering component). It is clear from these plots that in the QCD Monte-Carlo models the “transMAX” is more sensitive to the “hard scattering component” of the “underlying event” while “transMIN” is more sensitive to the “beam-beam remnants”, especially at large values of $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$. For example, for HERWIG at $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`=40\mathrm{GeV}/\mathrm{c}`$ the hard scattering component makes up $`62\%`$ of the “transMAX” $`N_{\mathrm{chg}}`$ with $`38\%`$ coming from the “beam-beam remnants”. On the other hand, the hard scattering component makes up only $`42\%`$ of the “transMIN” $`N_{\mathrm{chg}}`$ with $`58\%`$ coming from the “beam-beam remnants” at $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`=40\mathrm{GeV}/\mathrm{c}`$. Taking difference between “tansMAX” and “transMIN” does not completely remove the “beam-beam remnant” component, but reduces it to only about $`32\%`$ at $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`=40\mathrm{GeV}/\mathrm{c}`$.

## IV The Transverse Momentum Distribution in the “Transverse” Region

FIG. 17 shows the data on the transverse momentum distribution of charged particles in the “transverse” region for $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`>2\mathrm{GeV}/\mathrm{c}`$, $`5\mathrm{GeV}/\mathrm{c}`$, and $`30\mathrm{GeV}/\mathrm{c}`$. Each point corresponds to $`dN_{\mathrm{chg}}/dp_T`$ and the integral of the distribution gives the average number of charged particles in the “transverse” region, $`N_{\mathrm{chg}}`$, shown in FIG. 5. FIG. 5 shows only mean values, while FIG. 17 shows the distribution from which the mean is computed. The data are compared with the QCD hard scattering Monte-Carlo models predictions HERWIG. Since these distributions fall off sharply as $`p_T`$ increases, it is essentially only the first few points at low $`p_T`$ that determine the mean. The approximately constant plateau seen in FIG. 5 is a result of the low $`p_T`$ points in FIG. 17 not changing much as $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$ changes. However, the high $`p_T`$ points do increase considerably as $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$ increases. This effect cannot be seen by simply examining the average number of “transverse” particles. FIG. 17 shows the growth of the “hard scattering component” at large $`p_T`$ in the ”transverse region” (i.e., three or more hard scattering jets). For the QCD Monte-Carlo models the $`p_T`$ distribution in the “transverse” region, at low values of $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$, is dominated by the “beam-beam remnant” contribution with very little “hard scattering” component. This can be seen in FIG. 18 which shows both the “beam-beam remnant” component together with the total overall prediction of HERWIG for $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`>5\mathrm{GeV}/\mathrm{c}`$. For the QCD Monte-Carlo models the $`p_T`$ distribution in the “transverse” region, at low values of $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$, measures directly the $`p_T`$ distribution of the “beam-beam remnants”. Both ISAJET and HERWIG have the wrong $`p_T`$ dependence in the “transverse” region due to a “beam-beam remnant” component that falls off too rapidly as $`p_T`$ increases. It is, of course, understandable that the Monte-Carlo models might be slightly off on the parameterization of the “beam-beam remnants”. This component cannot be calculated from perturbation theory and must be determined from data.

FIG. 19 shows both the “beam-beam remnant” component together with the overall prediction of HERWIG for $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`>30\mathrm{GeV}/\mathrm{c}`$. Here the QCD Monte-Carlo models predict a large “hard scattering” component corresponding to the production of more than two large $`p_T`$ jets. HERWIG, ISAJET, and PYTHIA all do well at describing the high $`p_T`$ tail of this distribution.

## V “Transverse Regions” Versus “Transverse Cones”

In a complementary CDF analysis Valeria Tano tano  has studied the “underlying event” in hard scattering processes by defining “transverse cones” instead of “transverse regions”. As illustrated in FIG. 20, the “transverse cones” (with radius in $`\eta `$-$`\varphi `$ space of $`R=0.7`$) are located at the same pseudo-rapidity as the leading jet but with azimuthal angle $`\mathrm{\Delta }\varphi =+90^{}`$ and $`\mathrm{\Delta }\varphi =90^{}`$ relative to the leading “jet”. In the cone analysis the “jet” is a “calorimeter jet” (charged plus neutrals) defined using the standard CDF cluster algorithm. Maximum (MAX) and minimum (MIN) “transverse” cones are determined, on an event-by-event basis, similar to the “transMAX” and “transMIN” regions described in Section III. Each “transverse cone” has an area in $`\eta `$-$`\varphi `$ space of $`\pi R^2=0.49\pi `$ (compared with $`0.67\pi `$). FIG. 21 shows data at $`1.8\mathrm{TeV}`$ on the average scalar $`p_T`$ sum of charged particles ($`p_T>0.4\mathrm{GeV}/\mathrm{c}`$, $`|\eta |<1`$) within the MAX and MIN “transverse cones” versus the transverse energy of the leading (highest $`E_T`$) “calorimeter jet” compared with the QCD hard scattering Monte-Carlo models predictions HERWIG. The “transverse cone” analysis covers the range $`50<E_T(\mathrm{calorimeter}\mathrm{jet}\mathrm{\#}1)<300\mathrm{GeV}`$, while the “transverse region” analysis examines only charged particles and covers the range $`0<`$$`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`<50\mathrm{GeV}/\mathrm{c}`$. One cannot directly compare the two analysis, but if one scales the low $`E_T(\mathrm{jet}\mathrm{\#}1)`$ points in FIG. 21 by the ratio of areas $`0.67\pi /0.49\pi =1.36`$, one gets approximate agreement with the high $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$ points in FIG. 13. Both FIG. 13 and FIG. 21 indicate that the $`P_T\mathrm{sum}`$ of charged particles generated by HERWIG is slightly too small. Together the two CDF analyses give us a good handle on the “underlying event” in hard scattering processes.

## VI Tuning PYTHIA to Fit the “Underlying Event”

Now that we have constructed collider observables that are sensitive to the “underlying event” we would like to tune the multiple parton interaction parameters of PYTHIA to fit the data. There are many tunable parameters. Here we consider only the parameters given in Table 1. The default values of the parameters are given in Table 2. Note that the PYTHIA default values sometimes change as the version changes new\_pythia .

FIG. 22 and FIG. 23 show data on the average number of charged particles in the “transverse” region compared with the QCD Monte-Carlo predictions of PYTHIA 6.115 with different structure functions and different multiple parton interaction parameters and with $`p_T(\mathrm{hard})>0\mathrm{GeV}/\mathrm{c}`$. For PYTHIA the amount of multiple parton scattering depends on the parton distribution functions (i.e., the structure functions) and hence the number of particles produced in the “transverse” region (i.e., the “underlying event”) changes if one changes the structure functions. HERWIG and ISAJET do not include multiple parton scattering and for them the number of particles in the “transverse” is essentially independent of the choice of structure functions.

FIGS. 2429 show the results of a “tuned” version of PYTHIA 6.115 with MSTP(82) = $`4`$ and PARP(82) = $`2.4\mathrm{GeV}/\mathrm{c}`$ using the CTEQ4L structure functions. One must first choose a structure function and then tune the multiple parton scattering parameters for that structure function. In generating the PYTHIA curves in FIGS. 2429 we have taken $`p_T(\mathrm{hard})>0\mathrm{GeV}/\mathrm{c}`$. In general the perturbative $`2`$-to-$`2`$ parton scattering subprocesses diverge as $`p_T(\mathrm{hard})`$ goes to zero. PYTHIA regulates these divergences using the same cut-off parameters that are used to regulate the multiple parton scattering cross section (see Table 1). This allows for the possibility of using PYTHIA to simultaneously describe both “soft” and “hard” collisions. Most of the CDF Min-Bias events are “soft”, with less than $`3\%`$ of the events having $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`>5\mathrm{GeV}/\mathrm{c}`$. There is no clear separation between “soft” and “hard” collisions, but roughly speaking $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`<2\mathrm{GeV}/\mathrm{c}`$ corresponds to “soft” Min-Bias collisions and demanding $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`>5\mathrm{GeV}/\mathrm{c}`$ assures a “hard” collision. FIGS. 2427 show that the “tuned” version of PYTHIA with $`p_T(\mathrm{hard})>0\mathrm{GeV}/\mathrm{c}`$ describes fairly well the transition between “soft” and “hard” collisions. The QCD Monte-Carlo models with $`p_T(\mathrm{hard})>3\mathrm{GeV}/\mathrm{c}`$ cannot describe the data for $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`<3\mathrm{GeV}/\mathrm{c}`$ (see FIG. 5 and FIG. 6), whereas PYTHIA with $`p_T(\mathrm{hard})>0\mathrm{GeV}/\mathrm{c}`$ seems to do a good job on the “transverse” observables as $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$ goes to zero.

FIG. 28 shows the data on the transverse momentum distribution of charged particles in the “transverse” region compared with the “tuned” version of PYTHIA 6.115 (CTEQ4L, MSTP(82) = $`4`$, PARP(82) = $`2.4\mathrm{GeV}/\mathrm{c}`$). The fit is not perfect, but it is much better than the HERWIG prediction shown in FIG. 17. Multiple parton scattering produces more large $`p_T`$ particles in the “transverse” region, which is what is needed to fit the data. As seen in FIG. 29, the $`p_T`$ distribution in the “transverse” region, at low values of $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$, for the “tuned” version of PYTHIA is also dominated by the “beam-beam remnant” contribution as is the case for HERWIG (see FIG. 18). However, for PYTHIA the “beam-beam remnant” component includes contributions from multiple parton scattering, which results in a less steep $`p_T`$ distribution.

## VII Summary and Conclusions

The “underlying event” in a hard scattering process is a complicated and interesting object which involves aspects of both non-perturbative and perturbative QCD. Studying the “transMAX” and “transMIN” pieces of the “transverse” region provides additional information not contained in the sum. In the QCD Monte-Carlo models the various components that make up the “underlying event” are weighted differently in “transMAX” and “transMIN” terms. The “transMAX” term preferentially selects the “hard component” of the “underlying event” (outgoing jets plus initial and final-state radiation) while the “transMIN” term preferentially selects the “beam-beam remnant” component. Unfortunately one cannot cleanly isolate a single component of the “underlying event” since all components contribute to both “transMAX”, “transMIN”, and to the difference. However, requiring the Monte-Carlo models to fit both “transMAX” and “transMIN” (or the sum and difference) puts additional constraints on the way the generators model the “underlying event”.

ISAJET (with independent fragmentation) produces too many (soft) particles in the “underlying event” with the wrong dependence on $`P_T(\mathrm{chgjet}\mathrm{\#}1)`$. HERWIG and PYTHIA modify the leading-log picture to include “color coherence effects” which leads to “angle ordering” within the parton shower and they do a better job describing the “underlying event”. Both ISAJET and HERWIG have the too steep of a $`p_T`$ dependence of the “beam-beam remnant” component of the “underlying event” and hence do not have enough “beam-beam remnants” with $`p_T>0.5\mathrm{GeV}/\mathrm{c}`$. PYTHIA with multiple parton scattering does the best job at fitting the data.

The increased activity in the “underlying event” of a hard scattering over that observed in “soft” collisions cannot be explained solely by initial-state radiation. Multiple parton interactions provide a natural way of explaining the increased activity in the “underlying event” in a hard scattering. A hard scattering is more likely to occur when the “hard cores” of the beam hadrons overlap and this is also when the probability of a multiple parton interaction is greatest. For a soft grazing collision the probability of a multiple parton interaction is small. However, multiple parton interactions are very sensitive to the parton structure functions (PDF). You must first decide on a particular PDF and then tune the multiple parton interactions to fit the data.

One should not take the “tuned” version of PYTHIA 6.115 (CTEQ4L, MSTP(82) = $`4`$, PARP(82) = $`2.4\mathrm{GeV}/\mathrm{c}`$) presented here too seriously. It is encouraging that it describes fairly well the ”transverse” region over the range $`0<`$$`P_T(\mathrm{chgjet}\mathrm{\#}1)`$$`<50\mathrm{GeV}/\mathrm{c}`$ including the transition from “soft” to “hard” collisions. However, it is still not quite right. For example, it does not reproduce very well the multiplicity distribution of “soft” collisions. More work needs to be done in tuning the Monte-Carlo models butterworth . In addition, more work needs to be done before one can say that the multiple parton interaction approach is correct. HERWIG without multiple parton scattering is not that far off the data. Maybe we simply need to change and improve the way the Monte-Carlo models handle the “beam-beam remnant” component.


format: text

text:
# Long-range fractal correlations in literary corpora

## I Introduction

The human language faculty reveals as a phenomenon of remarkable complexity, intimately interwoven with all our superior mental functions . Its role in the modelling and categorisation of factual experience by the mind cannot be downplayed. Moreover, it has been suggested that the outburst of symbolic and abstract imagery pervading the archeological records signalling the transition between the middle and upper Paleolithic periods, might be interpreted as the hallmark of an associated, and rather sudden, transition in the linguistic performance of modern humans . The evolutionary advantage of complex syntactic communication stems from the capacity that it confers to language of coding complex information using moderate symbolic resources . In the traditional areas of linguistics much success has been attained in dissecting language structure from sentence level down the morphological rules of word formation. However, effective linguistic communication is a complex cognitive process  where sequences of sentences cohere in the assemblage of larger meaningful structures. At this point the study of language becomes susceptible of being complemented by using techniques that have been successfully applied on other complex systems in nature that also allow a variety of interdisciplinary approaches.

In particular, one question that still remains open to discussion is the ultimate origin of long-range correlations in complex information systems like the genetic code or human language. The main goal of this work is to examine the existence of long-range correlations in the use of words in written language drawing on methods from statistical time series analysis. To that end, we first describe the particular mapping scheme used to translate the given sequence of words present in a text sample into a time series. Then, we devote a section to explain the rescaled range analysis and how it can be used to infer the presence of long-range correlations in time series. Finally, the results of our systematic analysis are presented followed by concluding remarks.

## II The mapping of texts onto time series

The statistical analysis of long symbolic structures had an accelerated development right after the first complete DNA sequences started to become available, around a decade ago . Basically, the original techniques tried to create a random walk out of the particular symbolic sequence by mapping each distinct symbol onto a step in an independent direction in space. Thus, in a genetic sequence, made up of the concatenation of four independent nucleotides, represented by four letters, each instance of these symbols was translated onto an elementary jump along independent directions in a four dimensional space. Although different alternatives for the particular details of the mappings were put forward, all were devised to keep the original sequence’s structure at symbol level translated into the time series. Then, by applying a set of standard statistical techniques, such as power spectrum analysis, detrended fluctuation analysis or rescaled range analysis among others, it was possible to unravel the existence of long-range correlations in biosequences .

Written samples of natural languages, are also complex information carriers to be read sequentially. Thus, similar techniques as those employed in genetic sequence analysis were adapted with appropriate modifications . The different procedures shared one common point, namely, that the mappings were invariably performed at letter level. However, it has never been discussed in the literature so far whether the mapping at letter level is indeed adequate in the case of human language. It is crucial to note that the specific coding of words in some particular spelling or phonetic system is alien to the linguistic structure of communication, as is clear by noting that a given text can be readily coded in sign languages, for instance. Therefore, a better mapping, founded on linguistics basis, should recognise the symbols at the minimum level intrinsic to the communication process. In what follows we propose the simplest step in that direction, which consists in taking word tokens as the fundamental units of communication. By doing this we can assert that the mapped sequences do not carry any structure below word level.

In order to make clear the details of the actual translation of a given text onto a time series, a brief detour around Zipf’s analysis is required . It basically consists in counting the number of occurrences of each distinct word in a given text sample, and then producing a list of all these words ordered by decreasing frequency. Each word in this list can be identified by an index equal to its rank, $`r`$, in the list, that is, the most frequent word has index $`r=1`$, the second in the list is given $`r=2`$, and so on. Words that occur the same number of times are ordered arbitrarily in their corresponding rank interval.

By means of Zipf’s analysis each word in the original text can be replaced by its corresponding index $`r`$. Then, at position $`t`$ starting from the beginning of the text, we have the corresponding index $`r(t)`$. According to this mapping, the whole text becomes a time series of ranks, namely $`\{r(t)\}_{t=1}^T`$, where $`T`$ stands as the total length of the text. This particular numerical assignment may seem arbitrary at first, and in fact a different choice would certainly work as well. However, the selection of the rank as an indexing key to the words may be rendered more natural if we think of it as the assignment that minimises the effort in lexical access in the rank-ordered list of words in writing the whole text. When a word is required at position $`t`$ in the text, then it must be picked from position $`r(t)`$ in the list of words ordered by decreasing frequency.

In order to compare different time series generated by this mapping, it is useful to rescale the series to zero mean and unit variance. Thence, we define the standard deviation of the rank sample values as follows $`\widehat{\sigma }=\sqrt{r(t)^2r(t)^2}`$, where the symbol $`\mathrm{}`$ denotes an arithmetic average over the whole series of ranks. Then, the quantity $`\xi (t)`$ is defined as

$$\xi (t)=\frac{r(t)r(t)}{\widehat{\sigma }}.$$
(1)
In this way, the sequence of normalised increments $`\xi (t)`$ can now be regarded as the sequence of elementary jumps in a random walk each taken at a discrete time $`t`$.

## III Rescaled Range Analysis

Hydrology is the oldest discipline in which the presence of noncyclic very long-term dependence has been reported. Particularly, the rescaled range analysis was introduced by Harold E. Hurst  when he was studying the Nile in order to describe the long-term dependence of water levels in the river and reservoirs. Later, this statistical technique was further developed and applied by Mandelbrot and Wallis .

Let $`\xi (t)`$, $`1tT`$, be the normalised sequence of increments of a process in discrete time. In our case, $`\xi (t)`$ is the normalised ordered sequence of ranks from a text corpus of $`T`$ words, as described in the previuos section. From this sequence, the record

$$X(t)=\underset{u=1}{\overset{t}{}}\xi (u)$$
(2)
is constructed. Viewing the $`\xi (t)`$ as spatial increments in a one-dimensional discrete random walk, $`X(t)`$ is the position of the walker from the starting point at time $`t`$. For any given integer span $`s>1`$ and any initial time $`t`$, a detrended subrecord $`D(u,t,s)`$, for $`0us`$, can be defined as

$$D(u,t,s)=X(t+u)X(t)\frac{u}{s}\left(X(t+s)X(t)\right).$$
(3)
In this quantity, the mean $`_{w=1}^s\xi (t+w)/s`$ was substracted to remove the trend in the subrecord. The cumulated range $`R(t,s)`$ of the subrecord is defined by

$$R(t,s)=\underset{0us}{\mathrm{max}}D(u,t,s)\underset{0us}{\mathrm{min}}D(u,t,s)$$
(4)
and the variance $`S^2(t,s)`$ of the subrecord is defined by

$$S^2(t,s)=\frac{1}{s}\underset{w=1}{\overset{s}{}}\xi ^2(t+w)\left(\frac{1}{s}\underset{w=1}{\overset{s}{}}\xi (t+w)\right)^2.$$
(5)
For many time series of natural phenomena the average of the sample values of $`R(t,s)/S(t,s)`$, carried over all admissible starting points $`t`$ within the sample, follows the Hurst’s law: $`[R(t,s)/S(t,s)]s^H`$ with $`H>1/2`$. Hurst’s observation is remarkable considering the fact that in the absence of long-run statistical dependence one should find $`H=1/2`$, for processes with finite variance. For example, for a stationary Gaussian process $`\xi (t)`$ with $`\xi (t)=0`$ and $`\xi ^2(t)=1`$, Feller  has analytically proved that

$$\underset{s\mathrm{}}{lim}s^{1/2}[R(t,s)/S(t,s)]=\sqrt{\frac{\pi }{2}}.$$
(6)
Additionally, Mandelbrot and Wallis  showed that the $`s^{1/2}`$ law also applies to processes of independent increments having a variety of distributions: truncated Gaussian, hyperbolic, and (highly skewed) log-normal. Moreover, they also showed that when the increments of the process are statistical dependent but the dependence is limited to the short run, the $`s^{1/2}`$ law holds asymptotically. The effect of strong cyclic components was also studied. When a white Gaussian noise (of zero mean and unit variance) is superimposed with a purely periodic process of amplitude $`A`$, it can be seen that

$$\underset{s\mathrm{}}{lim}s^{1/2}[R(t,s)/S(t,s)]=\sqrt{\frac{\pi }{2}}\left(1+\frac{A}{2}\right)^{1/2}.$$
(7)
The value of the limit as well as the speed with which this limit is attained is dependent on $`A`$. When $`A`$ increases, $`[R(t,s)/S(t,s)]`$ takes an even longer time to reach the $`s^{1/2}`$ law. Moreover, the transition to the asymptote is non monotonic, but typically exhibits a series of oscillations.

The above comments support the key idea that $`s^{1/2}`$ law holds for every process for which long-term dependence is unquestionably absent and does not hold for processes exhibiting noncyclic long-term statistical dependence. Thus, a Hurst exponent of $`H=1/2`$ corresponds to the vanishing of correlations between past and future spatial increments in the record. For $`H>1/2`$ one has persistent behaviour, which means a positive increment in the past will on the average lead to a positive increment in the future. Conversely a decreasing trend in the past implies on the average a sustained decrease in the future. Correspondingly, the case $`H<1/2`$ denotes antipersistent behaviour.

In all the above discussion the focus was on the scaling properties of the ratio $`R(t,s)/S(t,s)`$. It is important to note that the introduction of the denominator $`S(t,s)`$ becomes ever more relevant for processes that deviate from the Gaussian and/or have long-term dependence . Futhermore, the ratio $`R(t,s)/S(t,s)`$ has a better sampling stability in the sense that the relative deviation, defined by $`\sqrt{Var[R(t,s)/S(t,s)]}/[R(t,s)/S(t,s)]`$, is smaller than any alternative expression used to study long-term dependence.

## IV Estimation of $`H`$ in literary corpora

In this section we shall explain the implementation of the rescaled range analysis, as well as present the results obtained after performing the experiments on texts coded as a normalised sequences of ranks.

As it was stated in the previous section, the method is based on the estimation of sample averages of the ratio $`R(t,s)/S(t,s)`$, for proper choice of different values of the integer span $`s`$. We shall follow the nomenclature and prescriptions introduced by Mandelbrot and Wallis  in the construction of the R/S diagrams from the experimental data. We first compute the exponent $`m`$ such that $`2^mT<2^{m+1}`$, and then we select the value of the span $`s`$ from the decreasing sequence of integers: $`\{T/2^p,p=0,1,\mathrm{},m2\}`$. For each $`s`$, we choose the starting points: $`t=qs+1,q=0,1,\mathrm{},2^p1`$ in order to construct $`2^p`$ nonoverlapping detrended subrecords $`D(u,t,s)`$ from the Eq. (3). Thus, for a given value of the index $`p`$, we have specified a value of $`\mathrm{log}s`$ which is marked on the axis of abscissas. In this way, we can reckon $`2^p`$ values of $`\mathrm{log}[R(t,s)/S(t,s)]`$ corresponding to the different starting points $`t`$, which are plotted as ordinates (marked by + signs). For each $`s`$, the logarithm of the sample average of the quantities $`R(t,s)/S(t,s)`$ over the different starting points is also marked in the R/S diagram. Then, the value of $`H`$ is estimated by linear regression of $`\mathrm{log}[R(t,s)/S(t,s)]`$ vs. $`\mathrm{log}s`$. The slope is calculated using the linear least-square method and the error is evaluated taking the uncertainty in the ordinates, associated with each value of the abscissa ($`\mathrm{log}s`$), equal to one quarter of the amplitude between the corresponding extreme points.

There are two pitfalls related to the R/S diagrams. First, for small values of $`\mathrm{log}s`$ (short subrecords) there is large scattering in the values of $`\mathrm{log}[R(t,s)/S(t,s)]`$. Second, for large values of $`\mathrm{log}s`$ we have very few nonoverlapping subrecords, thus narrowing the R/S diagram. This tightening involves a deceiving evaluation of uncertainty in the ordinate. Hence, when fitting a straight trendline, small and large values of $`\mathrm{log}s`$ should be neglected. Given a value of $`p`$ we have $`2^p`$ subsamples each of which have $`T/2^p2^{mp}`$ ranks. Therefore, we will retain the values of the span $`s`$ given by $`4pm4`$ as our criterion for fitting. In this manner, we will only consider average over 16 or more samples and each subrecord will have at least 16 ranks.

### A Results from literary corpora

In Fig. 1 we show the R/S diagram corresponding to the coded sequence consisting in 885534 ranks from 36 plays by William Shakespeare. The total number of different ranks, associated with different words in the corpus, is in this case equal to 23150. The mean and standard deviation ($`\widehat{\sigma }`$) of the original sequence and the third (M3) and fourth (M4) moments of the normalised sequence are shown in an inset. For simplicity’s sake, we only mark in the diagrams the minimum and maximum points for a given $`\mathrm{log}s`$, and the corresponding sample average is plotted as a small circle. The trendline is drawn as a solid line along the points effectively used in the fitting. The measured value of $`H`$ and its error are displayed in Table I, and confirm the presence of long-range correlation in the series. In the same figure, we can also see as small squares the diagram for the sample average from the sequence after deleting all ranks outside the interval $`(100,2000)`$. Strikingly, the corresponding value of $`H`$ is statistically indistinguishable from that of the original sequence. This fact tells us that the core of long-range correlation is neither supported by the most frequent words nor the least used.

There are two additional experiments which can provide information in tracking down the source of the correlated behaviour. The first one is a simple random shuffling of all the ranks in the sequence which has the effect of recasting Shakespeare’s plays into a nonsensical realisation, keeping the same original words without discernible order at any level. As is clear, after this experiment the analysis must yield a Hurst exponent indicating the uncorrelated character of the sequence. This is corroborated also in Fig. 1 where the sample average is plotted as small triangles. Yet more interesting is the analysis of a shuffling of Shakespeare’s plays that preserves sentence structure, and therefore English grammar. That is, by defining a sentence as the sequence of words between two periods, we can reorder them in a random fashion and thus produce a grammatically correct, though hardly meaningful, version of the corpus. The sample averages are plotted in this case as small diamonds. The resulting value of $`H`$ shows that grammar is not sufficient to induce long-range correlations as we can see again from Table I. Let us note that for maintaining the readability of the graph we only included the extreme R/S points of the original sequence.

At this point it may be illustrative to compare a plot of the record $`X(t)`$, the position of the walker as a function of time, both for the Shakespeare’s plays and for a stochastic sequence generated with $`H=0.7`$ by the successive random addition method . As can be seen in Figs. 2 and 3 both records show strong persistence that manifests itself in the long spans of average monotonous behaviour in the records.

In Figs. 4 and 5 we reproduce an identical analysis from two text corpora gathering a collection of works by Dickens and Darwin respectively, whose results are summarised in Table I. The Dickens sequence was obtained from 56 books by the author, and has 5616403 ranks (words) in length (44700 different ranks), whereas the Darwin sequence was obtained from 11 books and has a length of 1508483 ranks (30120 words in the vocabulary). It is worth noticing that the values of $`H`$ from the original sequences are indistinguishable for these two texts, written in prose and with different styles, although they are slightly greater than the value obtained for the Shakespeare sequence. However, the three values of $`H`$ corresponding to the truncation ($`100<\text{rank}<2000`$) are statistically equivalent. This fact suggests that the long-range correlations associated to words in the interval considered are a robust phenomenon over different styles and authors. Finally, the shuffling experiments show the same behaviour for the three authors, and are consistent with uncorrelated sequences.

### B Tests of consistence

In this subsection we present three tests of consistency of our analysis. First, we performed the same constructions as developed in the previous subsection over portions extracted from the original text corpora. Thus, we split the Shakespeare corpus into two parts of equal length and measured $`H`$ for each part. On the other hand, from the Dickens corpus we extracted a succession of 861038 words from an arbitrary origin and then measured the corresponding value of $`H`$. In Table II the values from both the original corpora and the parts are confronted. Values indistinguishable of $`H`$ result from the original corpus and its portions. In interpreting this result we should have in mind that Zipf’s analysis of a portion of the text generates a quite different table of ranks from the one corresponding to the entire corpus. This serves as an indication that we are quantifying a robust phenomenon inherent to the fractal structure of texts.

In the insets of Figs. 1 and 4 we present the values of mean and standard deviation ($`\widehat{\sigma }`$) corresponding to the set of ranks in the entire corpus of Shakespeare and Dickens respectively. For the Shakespeare corpus we can read that the mean of ranks is equal to 1003, whereas $`\widehat{\sigma }=2632`$. On the other hand, Fig. 4 shows for the Dickens corpus that the mean of ranks is equal to 1327 and $`\widehat{\sigma }=3793`$. These quantities are calculated from sample values of a probability distribution $`P(r)`$, for which Zipf’s law represents a crude approximation. However, by using more accurate analytical descriptions of the probability density $`P(r)`$, it is possible to evaluate the statistical moments for the whole distribution, that is, in the limit of infinite vocabulary. In particular, let us mention that the standard deviation calculated for both the Shakespeare and Dickens corpora, using the analytical expressions for $`P(r)`$ obtained in Ref. , result finite in the infinite vocabulary limit. The reason for this is that for large values of $`r`$, the probability density develops a fast exponential decay in the case of single-author corpora .

In Table I we also report the values of $`H`$ corresponding to corpora generated by means of stochastic processes. Particularly, we did the analysis explained above over a sequence of ranks generated by the Simon’s model  and other corresponding to a Markovian text , with a memory of seven letters. This short-range memory is enough to string out groups of few words in grammatical order. However, as we can see from the values on the table, we do not obtain correlations from neither type of sequence, which is consistent with the expected behaviour for processes with short-range correlations. Consequently, our implementation of the rescaled range analysis allows to clearly distinguish between real and stochastic version of texts .

### C The Dictionary

So far, we still do not have enough evidence as to assert a precise source for the phenomenon being analysed in this paper. We have characterised with a robust quantifier the presence of long-range correlations in literary texts beyond sentence level, and in fact over ranges spanning more than one individual work. The shuffling experiments attest where the long-range correlation does not originate, but say nothing on where it does. In previous attempts to study this phenomenon, a variety of possible origins for the long-range dependence have been put forward, though the conclusions in all cases have been inferred from the observed statistical behaviour of single literary works, and from the mapping of texts at letter level. In the present work we use a robust mapping with a strong footing on linguistics, and focus on the correlations that arise in large corpora comprising many individual books with no thematic linkage. Therefore, the existence of correlations overarching sets of entire works should be attributed to something else than either the relation between ideas expressed by the author as it was proposed in Ref. , or nonuniformities in the distribution of word’s lengths and the associated densities of blank spaces as it was suggested in Ref. . In particular, the latter alternative is ruled out from the outset since our mapping does not carry information on word structure.

In order to gain more insight into the origin of the attested long-range correlations, we decided to perform another experiment with a special type of text. A dictionary is a large collection of entries, one for each different word, arranged in alphabetical order. Clearly, in general, sequences of thematic affinity span at most a few entries, for small groups of words associated by common roots. Nevertheless, as is confirmed in Fig. 6 where we display the R/S diagram generated from the Webster’s abridged dictionary (1913 edition), the value of $`H`$ certifies the presence of long-range correlations. Table III compiles the exponent $`H`$ from the entire dictionary and from the text without the abbreviations (which indicate the grammatical category of entry words), obtaining in both cases totally consistent values for $`H`$. This result is rather striking and hints at the presence of a type of long-range organisation in the layout of lexical entries that may be connected with the phenomenon observed in literary works. This last assertion is further supported by performing a shuffle over the whole dictionary keeping integrity at entry level. That is, we do not disrupt the structure of entries made up of each word together with the associated definition, but we do mix them in a random manner. The resulting collection of definitions now lacks the alphabetical ordering; however, the explicit information contained in the dictionary is still intact. Notwithstanding, the long-range order is completely obliterated by the shuffle as accounted for by the fall in the $`H`$ value down to that of an uncorrelated sequence. This experiment corroborates that beyond the local structures derived from semantic affinity of small groups of related words, the ordered lexicon also possesses an overall macrostructure that emerges, as a hidden layer, out of the alphabetical order of entries.

A close examination at the structure of the dictionary reveals the source of the long-range order. Let us take for example the following sequence of related entries from the Webster’s abridged dictionary (1913 edition): advice, advisability, advisable, advisableness, advisably, advise, advised, advising, advisedly, advisedness, advisement, adviser, advisership, advisor, and advisory. These entries form a small cohesive block of information, and all their definitions share many words. In turn, these shared terms point to, possibly, remote locations in the dictionary where, again, words with common roots form small clusters. This process carries on to many levels of depth, thus building a complex network of relations at word level, which may be closely related to the presence of long-range correlation in the succession of words in the dictionary. The entry–level shuffling destroys that emerging order, and thereby the correlations.

## V Conclusions

In this work we have addressed the important issue of the emergence of long-range correlations in human written communication. To that end we proposed a simple mapping of texts onto random walks that keeps the word as the basic unit of communication. Therefore, the time series generated by this mapping retain the structure relevant to the linguistic phenomenon being analysed.

We addressed in some detail the rescaled range analysis, which has been successfully applied to the analysis of a vast variety of time series exhibiting long-term correlations. By applying this analysis to the coded texts we found conclusive evidence for long-range correlation in the use of words over spans as long as the whole corpora under study. It is worth noticing that the corpora used in this work are made up of the concatenation of independent literary works by individual authors, and, therefore the long-range effects must emerge as a phenomenon independent of the particular bounds of single literary works. This observation led us to analyse the case of a dictionary as a special kind of text, which provided insight into possible sources for the long-range order. The alphabetical order of entries in a dictionary corresponds only to the first visible layer of structural organisation. Yet, by performing the shuffling at entry level we realised that the whole of the long-range order is dependent upon a more complex structural layer related to the network of associations among word clusters.

As it has been convincingly shown in Ref.  sets of literary works also possess higher level structures associated with systematic patterns in word usage, which might also give rise to a complex network topology of relations among groups of words. The detailed mechanisms whereby the onset of these structures takes place in language requires further interdisciplinary research. In the light of the evidence supplied in this work, it is plausible that the ultimate source of these correlations is deeply related to the structural patterns of word distribution in written communication. In this view, groups of words associated by affine semantic hierarchies form a complex arrangement of cohesive clusters even over spans of entire corpora, thereby giving rise to long-range order in human written records.

###### Acknowledgements.

The authors are grateful to D. H. Zanette for his critical reading of the manuscript and to F. A. Tamarit for useful suggestions. This work has been partially supported by “Secretaría de Ciencia y Tecnología de la Universidad Nacional de Córdoba” (fellowship and grant 194/00). The digital texts analysed in this paper were obtained from Project Gutenberg Etext


text:
# Variation of Spot-induced Anomalies in Caustic-crossing Binary Microlensing Event Light Curves

## 1 Introduction

Microlensing experiments (MACHO: Alcock et al. 1993; EROS: Aubourg 1993) were originally proposed to search for Galactic dark matter in the form of massive compact halo objects by monitoring lensing-induced light variations of stars located in the Large Magellanic Cloud (Paczyński 1986). Besides this original goal, it was demonstrated that microlensing can also be applied to other fields of astronomy, especially stellar astrophysics (see the recent review of Gould 2001). Over the last decade, this aspect of microlensing achieved important progress thanks partially to the theoretical studies on various methods to extract additional information about lensed source stars and more importantly to the large number of event detections from additional experiments directed towards the Galactic bulge (OGLE: Udalski et al. 1993; DUO: Alard & Guibert 1997) and detailed light curves of events obtained from intensive followup observations (PLANET: Albrow et al. 1998; GMAN: Alcock et al. 1997; MPS: Rhie et al. 1999; MOA: Bond et al. 2001).

One of the lensing applications to stellar astrophysics is the detection and characterization of stellar spots. Spot detection via microlensing is possible for high magnification events produced by the source’s crossing of the lens caustic, which refers to the source position on which the magnification of a point source event becomes infinity. For these events, one can resolve the source star surface because different parts of the source is magnified by different amounts due to the large gradient of magnification over the source during the caustic crossing (Gould 1994; Nemiroff & Wickramasinghe 1994; Witt & Mao 1994). For a single lens, the location of the caustic is that of the lens itself (point caustic). For a binary lens, the set of caustics forms closed curves in which each curve is composed of three or more concave line segments (fold caustic) that meet at cusps. Heyrovský & Sasselov (2000) first investigated the possibility of lensing spot detections and showed that, for point-caustic-crossing single lens events, spots can cause fractional deviations in magnification larger than 2 %, which are detectable from followup observations. Han et al. (2000) further investigated the possibility of detecting stellar spots from the observations of fold-caustic-crossing binary lens events and showed that the fractional deviations are comparable to those of point-caustic-crossing events. However, these works were concentrated only on the possibility of spot detections, and thus no detailed investigation about the various forms of spot-induced anomalies in lensing light curves has been done. Bryce & Hendry (2000) mentioned some of the variations in their unpublished paper, but these analysis was only for point-caustic-crossing events, which are much less common than fold-caustic-crossing events. In addition, the previous studies are based on very simplified assumptions of spots and their host stars, e.g. a circular spot with a uniform surface brightness on also a uniform background stellar surface. Therefore, extensive study about the microlensing signature of stellar spots for caustic-crossing binary lens events based on realistic models of spots and host stars is required.

Another reason for the necessity of studying spot-induced light curve deviations was recently raised by Lewis (2001). He pointed out that if a source star possessing a planet is microlensed when the planet is transiting the source star surface, the resulting light curve will have similar deformation to those induced by a spot. Then, unless one can distinguish the two types of deviations, analyses based on one naively adopted assumption for the cause of the deviation might result in false information about the source star.

In this paper, by performing realistic simulations of microlensing light curves of events occurred on spotted source stars, we investigate how the pattern of the spot-induced anomalies in lensing light curves varies for different states of spots. Based on the results of this investigation, we also examine the feasibility of distinguishing the two types of anomalies caused by a spot and a transiting planet.

The paper is organized as follows. In § 2, we describe the basics of microlensing, which are required to describe the lensing behavior of caustic-crossing binary lens events occurred on spotted source stars. In § 3, by simulating lensing events occurred on source stars with various states of spots, we investigate the dependencies of the spot-induced anomaly pattern on various spot parameters, which characterize the physical state of spots. In § 4, we discuss the features of the spot-induced anomalies that can be used to distinguish the spot-induced anomalies from those caused by a transiting planet. We conclude in § 5.

## 2 Basics of Microlensing

If a point source located at $`𝒓_\mathrm{S}`$ on the projected plane of the sky is lensed by a binary lens system, where the individual lens components have masses of $`m_i`$ and projected locations of $`𝒓_{\mathrm{L},i}`$ ($`i=1`$ and 2), the positions of the resulting images, $`𝒓_\mathrm{I}`$, are obtained by solving the lens equation, which is expressed by

$$𝒓_\mathrm{S}=𝒓_\mathrm{I}\theta _\mathrm{E}^2\underset{i=1}{\overset{2}{}}\frac{m_i}{M}\frac{𝒓_\mathrm{I}𝒓_{\mathrm{L},i}}{\left|𝒓_\mathrm{I}𝒓_{\mathrm{L},i}\right|^2},$$
(1)
where $`M=m_1+m_2`$ are the total mass of the lens system and $`\theta _\mathrm{E}`$ represents the combined angular Einstein ring radius. The Einstein ring radius is related to the total mass and the geometry of the lens system by

$$\theta _\mathrm{E}=\left[\frac{4GM}{c^2}\left(\frac{1}{D_{\mathrm{OL}}}\frac{1}{D_{\mathrm{OS}}}\right)\right]^{1/2},$$
(2)
where $`D_{\mathrm{OL}}`$ and $`D_{\mathrm{OS}}`$ are the distances to the lens and the source, respectively. The magnifications of the individual images are given by the inverse of the Jacobian of the lens equation evaluated at each image position $`𝒓_{\mathrm{I},j}`$;

$$A_j=\left(\frac{1}{\left|detJ\right|}\right)_{𝒓_{\mathrm{I},j}};detJ=\left|\frac{𝒓_\mathrm{S}}{𝒓_\mathrm{I}}\right|,$$
(3)
and the total magnification is obtained by the sum of the magnifications of the individual images, i.e. $`A=_j^{N_\mathrm{I}}A_j`$, where $`N_\mathrm{I}`$ is the total number of images. The number of images is three when the source is outside the caustic curve. As the source crosses the caustic, two additional images appear and thus the number of images becomes five when the source is inside the caustic curve. Due to the infinite magnification during the source’s crossing of the lens caustic, the light curves of a caustic-crossing binary lens event is characterized by sharp spikes.

In reality, the lensed source star is not a perfect point source, and thus the lensing magnification does not become infinity even for a caustic-crossing event. For an event involved with an extended source, the magnification is given by the intensity-weighted magnification averaged over the source star surface, i.e.

$$A_0=\frac{I(𝒓)A(𝒓;𝒓_{\mathrm{L},i})𝑑\mathrm{\Sigma }_{}}{I(𝒓)𝑑\mathrm{\Sigma }_{}},$$
(4)
where $`𝒓`$ is the displacement vector of a point on the source surface with respect to the source center, $`I(𝒓)`$ is the surface brightness over the source star surface111Due to limb darkening, the surface brightness profile of a source star is not uniform but decreases with the increasing distance from the source center., $`A(𝒓;𝒓_{\mathrm{L},i})`$ is the point-source magnification at the point $`𝒓`$ caused by a binary lens with masses located at $`𝒓_{\mathrm{L},i}(i=1\mathrm{and}2)`$, and the notations $`\mathrm{}𝑑\mathrm{\Sigma }_{}`$ represents the surface integral over the source star surface. The existence of a spot (spots) causes further deviations in the lensing light curve. For an event involved with a spotted source star, the lensing light curve is represented by

$$A_{}=\frac{I(𝒓)A(𝒓;𝒓_{\mathrm{L},i})𝑑\mathrm{\Sigma }_{}f(𝒓^{})I(𝒓^{})A(𝒓^{};𝒓_{\mathrm{L},i})𝑑\mathrm{\Sigma }_{}}{I(𝒓)\left[1f(𝒓)\right]𝑑\mathrm{\Sigma }_{}},$$
(5)
where $`𝒓^{}`$ is the displacement vectors of a point on the spot with respect to the source center, $`f(𝒓^{})`$ represents the fractional decrement in the surface intensity due to the spot, and $`\mathrm{}𝑑\mathrm{\Sigma }_{}`$ represents the surface integral over the spot region of the source star.

## 3 Simulations

In this section, we investigate the pattern of spot-induced anomalies in lensing light curves resulting from various physical states of spots. For this, we perform simulations of events by using the formalism described in the previous section. The details of the simulations are described below.

### 3.1 Common Factors

Since we are interested mainly in the dependency of the anomaly pattern on the spot parameters characterizing the physical states of spots, the tested events in the simulations are produced by a common binary lens system with a common source trajectory involved with also a common source star (unspotted background part). The lens system tested in the simulations has equal mass components separated by $`\theta _\mathrm{E}`$. The source trajectory of the tested events is presented in Figure 1 along with the locations of the lens caustics. The source star has an angular radius $`r_{}=0.1\theta _\mathrm{E}`$ and its surface brightness profile is modeled by a linear form of

$$I_\nu (r)=1𝒟_\nu \left(1\sqrt{1(r/r_{})^2}\right),$$
(6)
where $`𝒟_\nu `$ is the limb-darkening coefficient whose value depends on the observed band. Assuming that events are observed both in $`B`$ and $`I`$ bands, the adopted values of the coefficients are $`𝒟_B=0.912`$ and $`𝒟_I=0.503`$222We consider multi-band observations because the color information might be useful to break the possible degeneracy between the anomalies caused by a spot and a transiting planet, as mentioned by Lewis (2001), but we note that the light curves to be presented in the following subsection are those measured in $`I`$ band., which corresponds to those of a K giant star with $`T_{\mathrm{eff}}=4750\mathrm{K}`$, $`\mathrm{log}g2.0`$, and a metallicity similar to the sun, i.e. $`[\mathrm{Fe}/\mathrm{H}]0`$ (Van Hamme 1993). Since the spin rotation period of a giant star is much larger than the caustic sweeping time, we do not consider the effect of changing spot position during caustic crossings in our simulations.

### 3.2 Spot Models and Resulting Anomaly Patterns

The physical state of a spot is characterized by various factors. These factors include the size, the shape, and the surface brightness contrast with respect to the unspotted stellar region. In addition, spots are likely to appear in groups and may have umbra/penumbra structures. To investigate the individual effects of these factors on the pattern of spot-induced anomalies in lensing light curves, in each set of simulations, we produce a series of light curves by varying one spot parameter but fixing other parameters.

First, we investigate the dependency of the spot-induced anomaly pattern on the surface brightness contrasts of spots. In Figure 2, we present the light curves produced with spot models having different contrast parameters. The contrast parameter is defined by $`C=I_{}/I_{}`$, where $`I_{}`$ and $`I_{}`$ are the surface brightnesses of the source star at the spot position with and without the presence of the spot and it is related to the decrement factor $`f`$ in eq. (5) by $`C=1/(1f)`$. With this definition of the contrast parameter, the anomaly pattern induced by a transiting planet is produced by setting $`C=0.0`$, because the planet completely blocks the light from the source star. In the simulations, we assume that the spots are located at the center of the source star and have a common radius of $`r_{}=0.2r_{}`$ (and thus covers 4% of the source star surface area). The time is expressed in units of the time scale that is required for the source to transit $`\theta _\mathrm{E}`$ (Einstein time scale, $`t_\mathrm{E}`$). To better show the deviation from the unspotted event light curve, we define a fractional magnification deviation by

$$ϵ=\frac{A_{}A_0}{A_0},$$
(7)
and its change as a function of time (deviation curve) is presented in the lower panel. We test four different surface brightness contrast parameters of $`C=1.0`$ (no spot), 0.4, 0.2, and 0.0 (a completely dark spot or a transiting planet).333For comparison, we note that the contrast parameter of a typical solar spot (umbra part) is $`C0.21`$ in $`I`$ band (Allen 1973). We note that the deviation curve is asymmetric with respect to the time of source center’s crossing of the caustic because not only the magnifications but also the gradients of the magnification on the left and right sides of the caustic are different each other. We also note that the value of $`ϵ`$ is larger than zero for some parts of the light curve due to the decrease in the source baseline flux caused by the spot, and thus smaller value of the denominator in eq. (5). From the comparison of the light and deviation curves, one finds that darker spots, as expected, induce deeper hollows in the deviation curves.

Second, we investigate the dependency of the anomalies on the spot size. Figure 3 shows the light curves produced with three different spot sizes of $`r_{}=0.1r_{}`$, $`0.2r_{}`$, and $`0.3r_{}`$. For all these events, we assume that the spots are located at the center of the source and have a common contrast parameter of $`C=0.2`$. From the figure, one finds that as the spot size becomes bigger, the deviation curve becomes not only deeper but also broader.

Third, we probe the dependency of the anomaly pattern on the umbra/penumbra structure of spots. We test three different spot structure models where the area of the umbra part, $`\mathrm{\Sigma }_\mathrm{u}`$, is fixed while the area of the penumbra region has three different values of $`\mathrm{\Sigma }_\mathrm{p}=1.67\mathrm{\Sigma }_\mathrm{u}`$, $`2.5\mathrm{\Sigma }_\mathrm{u}`$, and $`5.0\mathrm{\Sigma }_\mathrm{u}`$. In each model, we assume that the umbra region covers $`4\%`$ of the total source surface area. The boundaries of the umbra and penumbra regions have a concentric circular shape centered at the center of the source star. The adopted values of the contrast parameters are $`C_\mathrm{u}=0.2`$ and $`C_\mathrm{p}=0.7`$ for the umbra and the penumbra parts, respectively. One finds that the existence of the penumbra causes the deviation curve to become broader, but the increase in depth is less important due to the small surface brightness contrast of the penumbra region. Although the curvature of the deviation curve changes when the caustic crosses the umbra/penumbra boundary (at around $`t=0.12`$ and 0.17 in the curve), the change is too small to be noticed.

Fourth, we investigate the anomaly pattern dependency on the shape of spots. The shape of a spot can differ from a circle either due to the projection effect or because of its intrinsic shape. Moreover, spots tend to exist in a group and if small spots are closely located, the group as a whole will appear to be elongated. We test two elliptical spots with eccentricities of $`\epsilon =\sqrt{a^2b^2}/a=0.9`$ and 0.5, where $`a`$ and $`b`$ are the semi-major and semi-minor axes of the ellipse, respectively. We also investigate the effect of the spot orientation with respect to the sweeping caustic. For this, we test two cases for each spot shape where the major axis of the spot is parallel with and normal to the fold caustic, respectively (see the spots on the icons of Figure 5). The area of each elongated spot is normalized so that it is same as that of a circular spot with $`r_{}=0.2r_{}`$. We assume that the surface brightness contrast is $`C=0.2`$. The resulting light curves are presented in Figure 5. From the figure, one finds that the deviation curve becomes deeper and narrower when the elongated direction of the spot is parallel with the sweeping caustic. This is because as the elongated direction of the spot aligns more with the caustic, larger part of the spot crosses the caustic, but during shorter period of time.

Finally, we investigate how the multiplicity of spots affects the anomaly pattern. We test dual spot model, where the two spots have the same size, surface brightness contrast, and shape. In Figure 6, we present the light curves resulting from three different sets of spot locations, which are marked on the corresponding icons. We assume that the spots are intrinsically circular, but they are off-centered and thus appear to be elliptical. The size of each spot is normalized such that the area is same as that of a circular spot with $`r_{}=0.2r_{}`$. We assume a surface contrast of $`C=0.2`$. One finds that except for the case where the two spots are aligned with the caustic, the individual spots cause separate hollows in the deviation curve.

## 4 Transiting Planets Versus Spots

As mentioned, the lensing light curve anomalies induced by a transiting planet can mimic those induced by a spot. However, there are several characteristics of a spot that make the spot differentiated from a transiting planet. Then, for some of the spot-induced lensing light curve anomalies, it may be possible to isolate them from the anomalies induced by transiting planets. One may think of several such distinctive characteristics of spots. First, unlike transiting planets, spots have in general non-zero surface brightness contrasts. Second, spots normally have non-circular shapes, while the silhouettes of transiting planets are circular. Third, spots are likely to have complex umbra/penumbra structures unlike the uniform contrast of transiting planets. Fourth, spots would appear, by nature, in groups of two or more, while the chance for two planets simultaneously transit their host star is extremely low. Based on the patterns of spot-induced anomalies presented in the previous section, however, one finds that the spot characteristics of the non-zero surface brightness, the non-circular shape, and the umbra/penumbra structure do not produce distinctive features in the anomaly pattern. The major effects of these spot characteristics are making the width and the depth of the resulting deviation curve change, but these variations can be reproduced by modulating the planet size. On the contrary, the multiple deviation feature in the anomaly pattern caused by the multiplicity of spots is unique only for the spot-induced anomalies. Since the multiplicity is a general tendency for spots, we think that the spot-induced anomalies, in many cases, can be identified from this feature.

Besides the anomaly pattern in lensing light curves, we additionally check the possibility to distinguish the two types of anomalies by using the color information obtained from multi-band photometry. For this check, we simulate color changes of events occurred on source stars with a spot and a transiting planet. For each case, we test two models where the spot (planet) is located at the center and the limb of the source star during the caustic crossing (see the icons in Figure 7). We assume that both the spot and the planet covers 4% of the source star surface area. We also assume that observations are carried out in $`B`$ and $`I`$ bands. We note that while the contrast parameter for the transiting planet is zero regardless of the observed band, the spot contrast parameter is wavelength dependent. The spot contrast parameter have different values depending on the observed wavelength band not only because the spot region of the stellar surface has a different temperature from that of the background stellar surface but also because in different wavelengths one sees different layers of the stellar atmosphere (e.g. in longer wavelength, one sees cooler upper layer of the atmosphere). We assume that the spot contrast parameters are $`C=0.05`$ and 0.21 in $`B`$ and $`I`$ bands, respectively, by adopting those of a typical solar spot (Allen 1973). In Figure 7, we present the color change curves obtained from the simulations. In th lower panel, we also present the difference (in units of magnitude) from the color change curve of the event produced without the spot or the planet. We note that due to the limb darkening of the source star, color change occurs even without the presence of the spot or the planet on the source star.444For the description of the general color change patterns of both the point- and fold-caustic-crossing microlensing events occurred on limb-darkened source star, see Han & Park (2001). Therefore, the difference presented in the lower panel represents the additional color changes induced by the spot and the planet. From the figure, one finds that the pattern of the additionally induced color changes takes various forms. However, the amount of these additional color changes are expected to be very small ($`<0.01`$ mag), and thus we think that color information will be useful only when the the uncertainty in the measured color is very small.

## 5 Conclusion

We have investigated the patterns of spot-induced anomalies in caustic-crossing binary lens events. For this purpose, we performed simulations of events with various models of the physical state of spots which is characterized by the surface brightness contrast, the size, the number, the umbra/penumbra structure, the shape, and the orientation with respect to the sweeping caustic. From these simulations, we learned that the spot-induced anomalies take various forms depending on these factors. We also examined the feasibility of distinguishing the two possibly degenerate types of anomalies caused by a spot and a transiting planet and found that the degeneracy in many cases can be broken from the characteristic multiple deviation feature in the spot-induced anomaly pattern caused by the multiplicity of spots.

This work was supported by a grant (2001-DS0074) from Korea Research Foundation (KRF).