text:
# Evidence for a Large Stellar Bar in the Low Surface Brightness Galaxy UGC 7321. Based on observations obtained at the German-Spanish Astronomical Center (DSAZ), Calar Alto, jointly operated by the Max-Planck-Institut für Astronomie Heidelberg and the Spanish National Commission for Astronomy.

## 1 Introduction

Nearly all disk galaxies have both a stellar disk and a bulge. Traditionally, bulges were seen as similar to small ellipticals: kinematically hot, spheroidal or triaxial, and old. Therefore, they probably formed during the first collapse of the proto-cloud or at least in an early stage of their evolution. Multicolor studies of highly inclined galaxies (eg. Balcells & Peletier, 1994; Peletier et al., 1999) support this picture for early-type disk galaxies: their bulges have ages similar to cluster ellipticals, though they also tend to be bluer than ellipticals of similar luminosity. However, there are galaxies which do not have a bulge component. These pure disk galaxies show clearly that the bulge formation is not a necessary outcome of the early galaxy formation process.

Other possible origins of bulges are merger scenarios where either a small satellite galaxy is accreted and spirals into the center of the pre-existing disk Aguerri et al. (2001) or where two equal mass galaxies Kauffmann et al. (1996) merge and form an elliptical structure. Later, this object again accretes gas, which builds the present day surrounding disk. Another currently favoured scenario is the bulge growth by bar-driven mass inflow. Starting with a pure disk galaxy, numerical simulations have shown that the disks are, almost always, unstable against bar formation. The existing bar is then able to move material into the center. A buckling vertical instability pumps the stars above the disk, resulting in a three-dimensional bulge-like object (eg. Friedli & Martinet, 1993). The growing central mass eventually weakens and dissolves the bar Norman et al. (1996); Sellwood & Moore (1999). Thus we are left with an unbarred disk galaxy possessing a spheroidal bulge as observed today. In this scenario, pre-existing spheroidal bulges would become contaminated by disk stars, which might explain why bulges and the inner parts of disks have similar colour Balcells & Peletier (1994).

Bars are associated with the appearance of the so called box- or peanut-shaped (b/p) “bulges” visible in edge-on disk galaxies. The b/p part of these bulges are almost certainly the resonant off-plane thickening of a bar Combes et al. (1990); Raha et al. (1991); Pfenniger & Friedli (1991); Patsis et al. (2002a). The connection between bars and b/p bulges has been observationally verified by gas kinematics Kuijken & Merrifield (1995); Bureau (1998), and Lütticke et al. (2000b) found a significant correlation between b/p bulges and bar signatures seen in cuts along, and parallel to, the major axis in NIR images.

Photometrically observed peanut-shaped contours therefore provide evidence for the thickening process caused by the buckling instability.

In this paper, we analyze the case of UGC 7321, a nearly edge-on Sd galaxy. Global properties of this galaxy are given in Table 1. UGC 7321 does not have a prominent bulge as one can infer from its contour map (Fig. 1). However, the outer isophotes in deep $`R`$-band images show peanut-type distortions. UGC 7321 provides a useful case for studying the origin of peanut-type isophotal distortions in a galaxy without a prominent bulge. It could provide clues on the very early phases of formation of a central dense component. A key question is whether the galaxy harbours a bar. UGC 7321 was studied in detail in a series of papers by L. Matthews Matthews et al. (1999); Matthews (2000); Matthews & Wood (2001). Matthews et al. (1999) found a red nuclear feature in the $`BR`$ colour map of UGC 7321, only a few arcseconds across, and raised the possibility that it could be a kinematically distinct disk subsystem analogous to a bulge.

The existence, definition, and properties of bulges in late-type spiral galaxies were recently described by Böker et al. (2003). Late-type bulges are far from uniform and care needs to be taken regarding what may rightfully be called a bulge. Our analysis (Sect. 3) shows that, if the peanut distortions are indeed the result of a bar, the inferred bar length is uncomfortably long as compared to other late-type galaxies. However, alternative origins for the peanut distortions, discussed in Sect. 4, fail to provide convincing explanations.

## 2 Data

Images of UGC 7321 were obtained during an observing run in May 2001 at the 2.2 m telescope of the Calar Alto observatory (Spain) equipped with CAFOS. The dithered observations of this edge-on galaxy were used as night sky flatfields for the original face-on imaging study without spending extra dark time and are now discussed here. The data reduction is described in detail in Pohlen et al. (2002a) and only briefly recalled here. The images were taken in a special R-band filter showing a nearly rectangular filter characteristic combined with a high peak efficiency. Landolt standard fields Landolt (1992) were observed to achieve a photometric calibration to the standard Johnson system (cf. Pohlen et al., 2002a). During the data reduction process, special emphasis was given to performing the crucial flatfielding and sky subtraction. The final image used for this study is a combination of six stacked 600 s exposures, which allows one to reliably trace the profiles and contours down to $`\mu 27.2`$ R-mag$`/`$$`\mathrm{}^{\prime \prime }`$.

## 3 Results

A contour map of the inner 120″of UGC 7321 ($`R`$-band) is shown in Fig. 2. The six outer isophotes clearly deviate from pure ellipticity, with a depression along the minor axis above and below the galaxy center. According to Lütticke et al. (2000a), this depression along the minor axis is the characteristic feature of a classical, type 1 peanut-shaped distortion. The presence of peanut-shaped distortions is made obvious in cuts parallel to the major axis. In Fig. 3 we show such cuts at 0″, $`\pm `$6″, $`\pm `$12″and $`\pm `$18″. The peanut-shape signature (the central ”dip” in the profile) is visible in the $`\pm 12^{\prime \prime }`$ and obvious in the $`\pm 18^{\prime \prime }`$ profiles. The maximum peanut distortion occurs around $`\pm 50`$″from the center.

The term “peanut-shaped distortion” is commonly used to describe a property of the galaxy’s bulge. Here, however, the galaxy clearly lacks a kinematically hot, spheroidal $`R^{1/4}`$ system, and appears to lack a “pseudobulge” Kormendy (1993) or “exponential bulge” Carollo et al. (1998); we show below that the galaxy has little central light in excess of the best-fit disk exponential model. Therefore, the peanut-shaped distortion in UGC 7321 seems to be a property of the galaxy’s disk.

Given the absence of a strong central light concentration, we model the light distribution of UGC 7321 as a pure disk with the following hypotheses: the face-on surface brightness profile is exponential with scalelength $`h_{\mathrm{in}}`$ , out to a break radius $`R_{\mathrm{br}}`$ , and again exponential outside $`R_{\mathrm{br}}`$ , with an outer scalelength $`h_{\mathrm{out}}`$$`<`$$`h_{\mathrm{in}}`$ . The vertical surface brightness profile is also taken to be exponential, with scaleheight $`z`$. We employ a slightly modified version (cf. Pohlen, 2001) of the three-dimensional fitting routine described in Pohlen et al. (2000) to derive a best-fit model of the galaxy’s isophotes.

Radial profiles of the resulting model overplotted on the data are shown in Fig. 3. The best-fit parameters are: inner scalelength $`h_{\mathrm{in}}`$$`=79.2`$″, mean break radius $`R_{\mathrm{br}}`$$`=136.5`$″($`+135.0`$″and $`138.0`$″), outer scalelength $`h_{\mathrm{out}}`$$`=15.5`$″, and vertical scaleheight $`z=19.7`$″. Assuming a distance of 10.0 Mpc, this yields $`R_{\mathrm{br}}`$$`=6.6`$ kpc, $`h_{\mathrm{in}}`$$`=3.8`$ kpc, and $`z=1.0`$ kpc. The three-dimensional modeling allows for a rather precise determination of the galaxy’s inclination in nearly edge-on views. We derive an inclination of $`i=88.5\mathrm{°}`$, a central surface brightness of $`\mu _0=21.3`$ R-mag$`/`$$`\mathrm{}^{\prime \prime }`$and $`\mu _{\mathrm{br}}=23.5`$ R-mag$`/`$$`\mathrm{}^{\prime \prime }`$at the break radius. The corresponding face-on values are: $`\mu _0=23.8`$ R-mag$`/`$$`\mathrm{}^{\prime \prime }`$and $`\mu _{\mathrm{br}}=25.7`$ R-mag$`/`$$`\mathrm{}^{\prime \prime }`$. A contour plot of the model light distribution is shown in Fig. 2. Comparison of the two panels in Fig. 2 shows that the standard horizontal- and vertical-exponential disk model lacks the outer peanut distortions seen in UGC 7321. Figure 3 shows model cuts parallel to the major axis overplotted on the corresponding cuts for UGC 7321. Here again, the peanut distortions in the data at $`R<50`$″ deviate from the model in the lowest, outermost two cuts ($`z=\pm 12`$″, $`z=\pm 18`$″). Two further aspects are worth emphasizing. First, the good match in the outer parts between data and model provides evidence that UGC 7321, like other LSB galaxies (Pohlen et al. (2002b)), has the truncated, two-slope disk structure found in high-surface brightness spirals (cf. Pohlen et al., 2002a). Second, the major-axis cut (top profile in Fig. 3) shows a central enhancement over the best-fit exponential model at $`\pm 40`$″. This structure indicates that an additional component exists besides the double-exponential light distribution in UGC 7321. This component is affected by the in-plane dust extinction on the order of 0.4 mag Matthews & Wood (2001). Using a simple-minded definition of bulges as any excess light over the inward extrapolation of the exponential disk profile (e.g. Carollo et al., 1999), we could call the inner excess “the bulge” of UGC 7321. Interestingly, the radial position of the peanut distortions ($`\pm 50`$″) is comparable to the radial extent of the inner brightening, suggesting a parallel with the isophotal structure of classical box-peanut bulges (cf. Lütticke et al., 2000b).

Alternatively, the inner brightening of the major-axis profile might trace a bar in UGC 7321. Its profile and radial extent are hard to infer from the major-axis profile given the extinction by dust. The properties of such a bar might be inferred from the peanut distortions if we assume that the peanut feature is produced by the resonant off-plane thickening of a bar (e.g. Combes et al. 1990). Following Lütticke et al. (2000b), we derive the length of the peanut structure (BPL$`=83`$″) as the radial distance between the maxima of the peanut distortion (in detail: $`39`$″ and $`+44`$″) measured in radial cuts.

For a sample of normal HSB galaxies, Lütticke et al. (2000b) find a mean value of $`2.7\pm 0.3`$ for the ratio of projected bar length (BAL) to BPL. Assuming a similar scaling, the projected length of the bar in UGC 7321 would be $`224\pm 25`$″ in diameter and $`112\pm 13`$″ in radius. The two extreme values for the same sample (BAL/BPL$`=2.2`$ or $`3.2`$) would allow projected bar radii in the range of $`91`$″ to $`133`$″.

The estimated bar length is quite large. It is comparable to the total galaxy diameter, and would indicate that most of the body of UGC 7321 comprises a bar seen side-on with a ratio of bar length to galaxy size of BAL/$`D_{25}=224/302=0.7`$ ($`D_{25}`$ observed 25mag$`/`$$`\mathrm{}^{\prime \prime }`$B-band diameter uncorrected for inclination, from Matthews et al., 1999). In comparison, for a sample of normal disk galaxies observed in the NIR, Lütticke et al. (2000b) derive a mean value of BAL/$`D_{25}=0.4\pm 0.2`$. $`D_{25}`$ is probably not a good indicator for comparing LSBs with normal disk galaxies, but any inclination correction for $`D_{25}`$ will increase the ratio BAL/$`D_{25}`$ even more.

The numbers given above depend somewhat on the precise definition used for bar length and peanut length. Lütticke et al. (2000b) define the beginning of the bar in radial profiles as the location where the bar extends above the extrapolated outer disk. If bars are defined to start at their shoulders or bumps, BAL will be a factor of 1.3 smaller. In addition, using the maximum of the peanut rather than the full structure, BPL will be also a factor of 1.2 larger. This results in BAL/BPL$`=2.1`$ or even BAL/BPL$`=1.7`$, which is closer to the value obtained in N-body simulations Pfenniger (1984). For any of the given definitions, the case remains that the inferred bar length for UGC 7321 is quite large.

Comparing BAL with the radial scalelength ($`h`$) instead yields a ratio of BAL/$`h=224/79=2.8`$. Even assuming the lowest measured BAL/BPL ratio of $`2.2`$ from Lütticke et al. (2000b), BAL/$`h`$ will still be $`182/79=2.3`$. This is comparable to a bar of an early type galaxy. Erwin (2003) derives a mean value of BAL/$`h=2.4\pm 1.0`$ for a sample of face-on S0 to Sab galaxies and obtained a ratio of BAL/$`h=1.2\pm 0.5`$ for Sd galaxies taken from a sample of Martin (1995). Another example is the edge-on Sb galaxy NGC 2424. It exhibits a ratio of BAL/$`h=2.6`$, taking BAL $`=76`$ from Lütticke et al. (2000b) and $`h=29`$ from Pohlen (2001)111 The scalelength $`h`$ is determined only one-dimensionally and is therefore just a lower limit. In the worst case it is off by a factor of $`+20`$% compared to the true 3D scalelength (cf. Pohlen, 2001), which allows a BAL/$`h`$ as high as 3.3 for NGC 2424.. Note that the apparently unknown orientation of a bar in an edge-on galaxy is restricted by the actual peanut appearance of the contours. Lütticke et al. (2000b) have shown that the bar-produced b/p structure appears peanut-shaped only for a small range of aspect angles between $`77\mathrm{°}90\mathrm{°}`$ for the associated bar, assuming a galaxy inclination of $`90\mathrm{°}`$.

UGC 7321 is classified as a non box/peanut-shaped (b/p) bulge (type 4) in the catalogue of Lütticke et al. (2000a) and is actually the first Sd galaxy with a peanut bulge. However, they classify “bulges” to be either elliptical (4), boxy (3+2), or peanut-shaped (1), and therefore concentrate only on the inner part of UGC 7321 (which is indeed type 4) visible in DSS (Digitized Sky Survey) images. The peanut structure here is much larger than for a typical galaxy with a b/p-bulge. In addition, the characteristic angle $`\theta `$ between the major axis and the radial ray passing through the maximum of the peanut distortion is much smaller than normal (cf. Lütticke, 1999; Shaw et al., 1990). For UGC 7321, $`\theta `$ is $`18\mathrm{°}`$ while the mean value for all galaxies with peanut-shaped bulges is $`39\mathrm{°}\pm 10\mathrm{°}`$ with a minimum of $`31\mathrm{°}`$ Lütticke (1999). The conclusion would be that we are seeing an edge-on bar without the additional (superimposed) spheroidal bulge component as is typical for earlier (S0-Sc) galaxy types. Therefore, this large, slightly peanut-shaped, outer structure is the bar and one has to consider that the term b/p-bulge is a combination of a b/p structure, which we have here, and a spheroidal bulge component.

## 4 Discussion

Inferring the presence of a bar from the weak isophotal distortion observed in UGC 7321 might be controversial, especially given the implied large bar size compared to the galaxy disk. Here, we examine what other evidence exists for or against the presence of the bar, and what other configurations may create the observed peanut distortion.

Strong support for our photometric bar detection comes from the recent analysis of HI-data for UGC 7321 by Uson & Matthews (2003). They find that the major axis position-velocity (P-V) profile shows a clearly visible “figure-of-eight” structure (cf. their Fig. 10). This is typical for a barred galaxy viewed edge-on Kuijken & Merrifield (1995). Uson & Matthews (2003) favour the explanation proposed by Brinks & Shane (1984) that this feature may also result from the warping and flaring of its HI disk. However, there is only a very mild — compared to other later type galaxies García-Ruiz et al. (2002) — HI warp on the west side222here: positive side. Uson & Matthews (2003) assign a decrease in the HI intensity at $`25\mathrm{}`$ to the possible presence of a small bar extending from about $`30\mathrm{}`$ to $`+40\mathrm{}`$ which nicely coincides333Note: We use a different orientation of the image. with our “bulge” feature in the major axis light profile. Nevertheless, if this is a bar it would be too small to account for the peanut-shaped outer structure. We do not find indications of a flaring in the stellar disk. By introducing a flaring in our model we confirm that it could not be responsible for the observed peanut distortion in the inner part. The effect of a radially increasing scaleheight is most prominent towards the edges of the disk.

Binney & Petrou (1985) and Whitmore & Bell (1988) have shown that the accretion of satellite galaxies is a possible formation process for b/p bulges. For normal b/p-bulges these scenarios are excluded due to the overwhelming evidence for bar-driven evolution (eg. Bureau, 2002; Lütticke et al., 2000b), where the role of interaction is only to trigger the bar itself. However, Lütticke et al. (2003) found a special class of thick boxy bulges (TBB) which suggest a merger origin. These galaxies are often disturbed, frequently showing prominent irregularities and asymmetries, and possess significantly more projected satellites compared to a control sample. All these characteristics are in contrast to UGC 7321’s smooth and undisturbed disk and its unusually isolated environment. There are no signs for a recent merger event, such as an increase in star-formation, asymmetric dust distribution, arcs, shells, filaments, or a strong stellar warp. Uson & Matthews (2003) estimate that a possible encounter $`10^9`$ years ago of UGC 7321 with its closest neighbour could be responsible for the observed mild HI warp.

Taken together, nearly all of the observational evidence points against a merger scenario. The peanut-shaped structures found in simulations (eg. Hernquist & Quinn, 1989) resemble typical early-type spirals and S0s such as NGC 128, IC 4767, or HCG 87A and are often located in dense groups; some of them are even “hybrid” scenarios where the bar is triggered by a merger (eg. Mihos et al., 1995).

Patsis et al. (2002a), analysing N-body models, describe a new vertical resonance producing an edge-on ‘boxy’ or ‘peanut’ structure for non-barred galaxies. However, the examples from the N-body simulation have only an oval distortion as their main morphological feature and could be characterized morphologically as S0 galaxies. In addition, their model allows boxes, but not strong peanuts or X-shaped structures.

Therefore, the presence of a large bar extending over the main body of the galaxy provides the best explanation for the peanut-shaped contours in UGC 7321 and the measured “figure-of-eight” pattern in the recent HI measurements.

A bar, especially one as big as we photometrically infer, in a low surface brightness galaxy, a system expected to be dark-matter dominated, could solve the problem that many rotation curves of LSB galaxies cannot (eg. de Blok & Bosma, 2002) be fitted by the cuspy halo profiles predicted by CDM cosmogonies. Specifically, Swaters et al. (2003) discuss non-circular, bar-induced kinematics as the source of deviations between rotation curve-based mass models and the NFW-type CDM halos. Weinberg & Katz (2002) propose that a dark matter cusp might be removed due to dynamical interaction of the halo with the bar. However, several works argue against such a mechanism and in favor of non-circular kinematics to explain the deviations. For example, Mayer & Wadsley (2003) used high resolution N-Body/SPH simulations to address the question of bar formation and evolution in LSB galaxies. They conclude that LSB disks may go bar-unstable in gas-rich disks embedded in low-concentration halos. However, the resulting bars are too short and probably of too low angular momentum content to affect the inner halo density profile. They also conclude that a bulge-like component must be present in any LSB galaxy that has got bar unstable. In the case of UGC 7321 the inner component could be tracing this process, so that we are observing just the beginning of the morphological evolution from its first bar into a bulge-like central structure, according to the secular evolution scenario.

###### Acknowledgements.

Part of this work was supported by the German *Deutsche Forschungsgemeinschaft, DFG*. We would like to thank L. Matthews for stimulating discussions about UGC 7321 and D. Pfenniger and P. Patsis for their views about a possible bar scenario.


format: text

text:
# Sub and Super-Luminal Propagation of Intense Pulses in Media with Saturated and Reverse Absorption

## Abstract

We develop models for the propagation of intense pulses in solid state media which can have either saturated absorption or exhibit reverse absorption . We show that the experiments of Bigelow et al.\[Phys. Rev. Lett. 90, 113903 (2003); Science 301, 200 (2003).\] on subluminal propagation in Ruby, and superluminal propagation in Alexandrite are well explained by modelling them as three level and four level systems respectively, coupled to Maxwell equations. We present results well beyond the traditional pump-probe approach.

Since the discovery of ultraslow light with a group velocity 17 m/sec in a Bose condensate by Hau et al. many experiments have reported slow light in a variety of media . Kash et al.  demonstrated light propagation with a group velocity of 90 m/sec at room temperature in Rb vapor. Using Zeeman coherences, Budker et al.  reported slow light with group velocity 8 m/sec in Rb vapor. Hemmer et al.  first reported slow light in solid state material namely Pr doped Y2SiO5, maintained at a cryogenic temperature of 5K. The slow light ideas have been successfully used in storage and retrieval of light pulses . The understanding of storage and retrieval of light pulses has been provided by Dey and Agarwal , using the adiabaton theory of Grobe, Hioe and Eberly .

Work on pulse propagation continues to produce interesting results. Recently, Bigelow et al.  showed the propagation of light pulses in Ruby at a group velocity of 57.5 m/sec. This experiment differs considerably from all earlier experiments which were based on electromagnetically induced transparency . Bigelow et al. recognize that a two level system driven by a strong field and a probe gives rise to a hole in the probe response function with a width of the order of $`1/T_1`$, where $`T_1`$ is the longitudinal relaxation time . Note that in a material like Ruby the transverse relaxation time $`T_2T_1`$. These authors also discovered that they need not use separate pump and probe fields. A field with peak power of the order of saturation intensity could be slowed down considerably to about 57.5 m/sec. Further Bigelow et al.  extended their work to a material like Alexandrite, where they reported superluminal propagation. There have been several earlier reports of superluminal propagation in solid state systems ; and in vapors .

The purpose of this letter is to study the propagation of intense pulses in a homogeneously broadened medium, like Ruby, which can exhibit saturated absorption or a medium like Alexandrite, which can exhibit reverse absorption. Note that the pulse propagation in a nonlinear transparent medium has been extensively studied. The systems studied here differ from the previous studies as our systems posses very strong transverse and longitudinal relaxation effects. In order to model the experiments, we model Ruby as a three level system and Alexandrite as a four level system. We solve the system of coupled equations numerically to delineate various aspects of pulse propagation. We do not make any approximation on the strength of the pulses so that we can model experimental observations on strong pulses. We calculate group velocity from the relative delay or advancement between the reference pulse and the output pulse. We present numerical results on the propagation of Gaussian and modulated pulses and show good agreement with the experimental data of Refs.. The experiments of Bigelow et al. fall in two categories. One consists of a weak probe and a strong coherent cw field. These were explained in terms of the response to a probe field of a two level medium pumped by a coherent field. The other category consists of the self delay of a strong pulse. The latter requires solutions of the coupled Maxwell-Bloch equations and this is the case we concentrate on.

Media with saturated absorption—-For pulse propagation in the three level model in Ruby, we represent the ground state as $`|g`$, the $`{}_{}{}^{4}F_{2}^{}`$ absorption band as $`|e_1`$ and the levels $`2\overline{A}`$ and $`\overline{E}`$ as $`|e_2`$. In Ruby one has very rapid decay of the level $`|e_1`$ to $`|e_2`$, and as a result some of the coherences become irrelevant on experimental time scale. The density matrix equations for the model of the Fig. 1 are

$`\dot{\rho }_{_{gg}}`$ $`=`$ $`2\mathrm{\Gamma }_2\rho _{_{_{22}}}+i\mathrm{\Omega }(\rho _{_{_{1g}}}\rho _{_{_{g1}}})`$ (1)
$`\dot{\rho }_{_{_{22}}}`$ $`=`$ $`2\mathrm{\Gamma }_1\rho _{_{_{11}}}2\mathrm{\Gamma }_2\rho _{_{_{22}}}`$ (2)
$`\dot{\rho }_{_{_{1g}}}`$ $`=`$ $`\mathrm{\Gamma }_1\rho _{_{_{1g}}}+i\mathrm{\Omega }(\rho _{_{gg}}\rho _{_{_{11}}})`$ (3)
$`\rho _{_{gg}}`$ $`+`$ $`\rho _{_{_{11}}}+\rho _{_{_{22}}}=1,`$ (4)
where $`\rho _{_{_{ij}}}=e__i|\rho |e__j`$, $`i,j=1,2`$. The Rabi frequency $`2\mathrm{\Omega }`$ is defined by $`2\mathrm{\Omega }(z,t)=2\stackrel{}{d}_{_{1g}}\stackrel{}{}(z,t)/\mathrm{}`$, where $`\stackrel{}{d}_{_{1g}}`$ is the dipole matrix element and $`\stackrel{}{}(z,t)`$ is the envelop of the pulse. We assume that the carrier frequency, $`\omega `$, is on resonance with the frequency of the $`|e_1|g`$ transition. Under the approximations, $`\mathrm{\Gamma }_1\mathrm{\Gamma }_2,\mathrm{\Omega }`$$`\dot{\rho }_{_{1g}}0`$, we derive the approximate equation for the evolution of the ground state population as

$$\dot{\rho }_{gg}=2\mathrm{\Gamma }_2(1\rho _{gg})2\frac{\mathrm{\Omega }^2}{\mathrm{\Gamma }_1}\rho _{gg}$$
(5)
Note that we can prove that $`\dot{\rho }_{_{11}}0`$, if $`\mathrm{\Gamma }_1\mathrm{\Gamma }_2,\mathrm{\Omega }`$. Under the same conditions and the slowly varying envelop approximation, the evolution equation for the Rabi frequency of the field is governed by

$$\frac{\stackrel{~}{\mathrm{\Omega }}}{z}=\frac{\alpha _0}{2}\stackrel{~}{\mathrm{\Omega }}\rho _{gg},\stackrel{~}{\mathrm{\Omega }}=\mathrm{\Omega }/\mathrm{\Omega }_{sat},$$
(6)
where $`\alpha _0=4\pi \omega |d_{_{1g}}|^2/c\mathrm{}\mathrm{\Gamma }_1`$ and $`\mathrm{\Omega }_{sat}=2\sqrt{\mathrm{\Gamma }_1\mathrm{\Gamma }_2}`$. In Eqs. (2) and (3) we have used the pulse coordinates i.e, $`tz/c,z`$. The time derivative in Eq. (2) is with respect to $`(tz/c)`$. The time $`t`$ can be expressed in units of $`1/2\mathrm{\Gamma }_2`$. For numerical computation, we consider two types of input pulses, viz, a Gaussian pulse with a temporal width $`1/\mathrm{\Gamma }_2`$

$$\stackrel{~}{\mathrm{\Omega }}_{in}=\stackrel{~}{\mathrm{\Omega }}^0e^{[t^2/2\sigma ^2]}$$
(7)
and amplitude modulated pulse

$$\stackrel{~}{\mathrm{\Omega }}_{in}^2(t)=I=I_0(1+m\mathrm{cos}[\mathrm{\Delta }t]).$$
(8)
The Equations (2) - (5) are our working equations. We use these for numerical computations. We calculate the evolution of the pulse for arbitrary values of $`\stackrel{~}{\mathrm{\Omega }}^0`$ or $`I_0`$. Some typical results for the Gaussian pulses are shown in the Fig. 2. We get group velocities in the range 50 m/sec for $`\mathrm{\Omega }/\mathrm{\Omega }_{sat}1`$ and the transmission is rather small. In Fig. 3, we exhibit the behavior of $`v_g`$ and transmission as a function of the input intensity. These results are in agreement with the experimental findings of transmission in the range $`0.1\%`$. In Fig. 3, we also show for comparison the results of the group velocity and the transmission for the propagation of an intense pulse through a two level system described by the traditional Bloch equations. The coupled Maxwell-Bloch equations under the approximation $`T_1T_2`$ are given by

$`\dot{\rho }_{_{11}}T__1`$ $`=`$ $`\rho _{_{11}}+2|\stackrel{~}{\mathrm{\Omega }}|^2(12\rho _{_{11}})`$ (9)
$`\rho _{_{gg}}`$ $`+`$ $`\rho _{_{11}}=1`$ (10)
$`{\displaystyle \frac{\stackrel{~}{\mathrm{\Omega }}}{z}}`$ $`=`$ $`{\displaystyle \frac{\alpha _0}{2}}\stackrel{~}{\mathrm{\Omega }}(12\rho _{_{11}}),`$ (11)
where $`\stackrel{~}{\mathrm{\Omega }}=\mathrm{\Omega }\sqrt{T_1T_2}`$ and dot denotes $`/(tz/c)`$. As seen from the Fig. 3, there are substantial differences in the propagation of pulses in two level and three level media. Note that the time $`T_1`$ is equal to $`1/2\mathrm{\Gamma }_2`$. We believe that, in the light of the energy level diagram of Ruby, it is more appropriate to model it as a three level system.

We next consider input pulse as a modulated pulse given in Eq. (5). The output pulse is modulated with a phase shift (time delay). We show this time delay as a function of modulation frequency for two different pump powers in the Fig. 4. The results in Figs. 2 - 4 are in excellent agreement with the experimental data(cf for example Fig. 3 of Ref. with our Fig. 4).

Media with reverse saturation —-For the superluminal propagation in Alexandrite, Bigelow et al. recognized how the reverse saturation mechanism  can be at work in a material like BeAl2O4 doped with Cr3+ ions and with some Cr3+ ions replaced by Al3+. The reverse saturation produces an antihole in the susceptibility for the probe in presence of a pump field. The antihole can result in the superluminal propagation. In what follows we show how the measurement can follow by modelling the system as a four level system to account for reverse absorption. The model is shown in the Fig. 5, where state $`{}_{}{}^{4}A_{2}^{}`$ as $`|g`$, the absorption bands $`{}_{}{}^{4}T_{2}^{}`$ and $`{}_{}{}^{4}T_{1}^{}`$ as $`|e_1`$ and the level $`{}_{}{}^{2}E`$ as $`|e_2`$. The density matrix equations are now given by

$`\dot{\rho }_{_{_{gg}}}`$ $`=`$ $`2\mathrm{\Gamma }_2\rho _{_{_{22}}}+i\mathrm{\Omega }(\rho _{_{_{1g}}}\rho _{_{_{g1}}})`$ (12)
$`\dot{\rho }_{_{_{22}}}`$ $`=`$ $`2\mathrm{\Gamma }_1\rho _{_{_{11}}}2\mathrm{\Gamma }_2\rho _{_{_{22}}}+2\mathrm{\Gamma }_3\rho _{_{_{33}}}`$ (13)
$`\dot{\rho }_{_{_{33}}}`$ $`=`$ $`2\mathrm{\Gamma }_3\rho _{_{_{33}}}+i\mathrm{\Omega }(\rho _{_{_{23}}}\rho _{_{_{32}}})`$ (14)
$`\dot{\rho }_{_{_{32}}}`$ $`=`$ $`\mathrm{\Gamma }_3\rho _{_{_{32}}}+i\mathrm{\Omega }(\rho _{_{_{22}}}\rho _{_{_{33}}})`$ (15)
$`\dot{\rho }_{_{_{1g}}}`$ $`=`$ $`\mathrm{\Gamma }_1\rho _{_{_{1g}}}+i\mathrm{\Omega }(\rho _{_{gg}}\rho _{_{_{11}}})`$ (16)
$`\rho _{_{gg}}`$ $`+`$ $`\rho _{_{_{11}}}+\rho _{_{_{22}}}+\rho _{_{_{33}}}=1.`$ (17)
Following the same procedure as in the case of Ruby, we have derived the working equations

$`{\displaystyle \frac{\dot{\rho }_{_{gg}}}{2\mathrm{\Gamma }_2}}`$ $`=`$ $`(1\rho _{_{_{gg}}})2\stackrel{~}{\mathrm{\Omega }}^2\rho _{_{_{gg}}}`$ (18)
$`{\displaystyle \frac{\stackrel{~}{\mathrm{\Omega }}}{z}}`$ $`=`$ $`{\displaystyle \frac{\alpha _0}{2}}\stackrel{~}{\mathrm{\Omega }}\rho _{_{_{gg}}}{\displaystyle \frac{\stackrel{~}{\alpha _0}}{2}}\stackrel{~}{\mathrm{\Omega }}(1\rho _{_{_{gg}}}),`$ (19)
where $`\stackrel{~}{\alpha }_0`$ gives the reverse saturation. Following the experimental data of Bigelow et al., we estimate $`(\stackrel{~}{\alpha }_0/\alpha _0)4`$ . The Eqs. (10) and (11) are numerically integrated for the input Gaussian pulse given by Eq. (4). A representative set of results is shown in the Fig. 6. This Figure also shows how the group velocity and net transmission depends on the peak intensity of the Gaussian pulses. It should be borne in mind that in the range of the intensities of Fig. 6, no perturbation theory can be used. One has to study the full nonlinear behavior. We also notice that the input pulses get distorted in shape. The distortion becomes more pronounced as the nonlinearity of the medium becomes more pronounced.

In conclusion, we have shown how to model the propagation of intense pulse in solid state media with very strong relaxation effects. The media can exhibit either saturated absorption or reverse absorption.Our modelling goes well beyond the traditional pump-probe approach. We specifically present results on the propagation of pulses in Ruby and Alexandrite. Our model would also be applicable to other systems where reverse absorption could be dominant.

GSA is grateful to R. Boyd for extensive discussions on experiments and the contents of this paper. GSA also thanks E. Wolf for partially supporting visit to Rochester through US Air Force Office of Scientific Research, grant no. F49620-03-1-0318 and the Engineering Research Program of the Office of Basic Energy Sciences at the US Department of Energy under grant no. DE-FG02-ER 45992.


text:
# Contents

## I Introduction

Parton distributions provide the essential link between the theoretically calculated partonic cross-sections, and the experimentally measured physical cross-sections involving hadrons and mesons. This link is crucial if we are to make incisive tests of the standard model, and search for subtle deviations which might signal new physics. Since perturbative calculations are renormalization scheme dependent,111We shall use the term *renormalization scheme* in the general sense to include both the usual renormalization (ultra-violet subtraction) scheme, and the factorization (infrared and collinear subtraction) scheme. it is important to use properly matched hard-scattering cross-sections and parton distribution functions (PDFs) in evaluating factorized cross-sections for physical applications. This issue is particularly relevant for applications involving heavy quarks, since the heavy quark introduces a new mass scale which leads to complications of the PQCD formalism.

For most physical applications, it is convenient (and a good approximation) to use hard-scattering cross-sections calculated in the zero parton-mass limit. In fact, next-to-leading-order hard-scattering cross-sections in the *non-zero quark-mass case* have only been calculated for some very basic processes, such as deep inelastic scattering (DIS); they are not yet available for most interesting “new physics” applications. Thus, the most useful parton distributions for general applications are those determined in a global QCD analysis using the zero quark-mass approximation. This was the choice made for the latest series of CTEQ6 parton distributions: CTEQ6M, CTEQ6D, and CTEQ6L. cteq6  For definiteness, this scheme will be referred to as the zero-mass variable-flavor-number scheme (ZM-VFNS).

This paper extends Ref. cteq6  by performing a similar global analysis using the generalized (non-zero quark-mass) $`\overline{MS}`$ perturbative QCD framework of Refs. ColTun ; AOT ; ACOT ; Collins98 ; ASTW , which we label the general-mass variable-flavor-number scheme (GM-VFNS). When matched to the corresponding hard-scattering cross-sections calculated in the same scheme, the combination should provide a more accurate description of the precision DIS structure function data, as well as other processes which are sensitive to charm and bottom mass effects. The main result we present here is a new set of parton distributions evaluated in this GM-VFNS which we identify as CTEQ6HQ.222A fit named CTEQ6F3 employing the three fixed-flavor-number scheme (3-FFNS) for DIS (cf., Refs. F3NLO ; grs ) which is suitable for special applications in that scheme, will be presented separately.

Perturbative QCD with non-zero quark-masses is well-established in the literature. Witten:bh ; CWZ ; ColTun ; ACOT ; Collins98 ; BuzaEtal ; ChuvakinEtal  Its implementation is, however, necessarily more complicated than the corresponding zero-mass case. In addition to the familiar renormalization-scheme dependence, more ambiguities arise from the implementation of the charm and bottom thresholds (in relation to their masses), and the detailed way the 3-flavor, 4-flavor, and 5-flavor renormalization schemes are matched at the transition scales. Several schemes proposed in the recent literature are formally equivalent to each other at high energies; but, they can differ in the threshold region, as well as in the intermediate energy range, where much of the experimental data lie. ACOT ; BuzaEtal ; ChuvakinEtal ; ASTW ; NasonEtal ; TKS  Some of these schemes are more natural and numerically robust, e.g. insensitive to variations of the remaining parameters (such as scale), than others. We follow the procedure proposed in Refs. TKS ; acotchi  which has been shown to be particularly stable as it naturally takes into account the dominant (logarithmic and non-logarithmic) mass-dependent kinematic effects with a physically motivated choice of the scaling variable.

Section 2 summarizes the method of analysis. Because this work represents an extension of the zero quark-mass CTEQ6 global QCD analysis, we shall focus mostly on the new elements related to non-zero quark mass; details that are in common with the massless case can be found in Ref. cteq6 . Section 3 describes the global analysis of data. Section 4 compares the new parton distributions with data. Section 5 makes comparisons between the CTEQ5M, CTEQ6M and CTEQ6HQ PDFs. At asymptotically high energies, mass effects become negligible; therefore, all generalized parton model prescriptions should approach the same zero quark-mass limit. We show that the differences between CTEQ6M and CTEQ6HQ PDFs indeed vanish with increasing $`Q`$, as expected. This section also contains a discussion of applications and issues related to uncertainties of physics predictions resulting from these parton distributions.

## II The Generalized $`\overline{MS}`$ Scheme with non-zero mass heavy flavor partons

The main feature of this analysis, compared to the standard zero-mass parton CTEQ6 PDF sets, is the adoption of the generalized perturbative QCD framework incorporating non-zero mass effects for the charm and bottom quarks.333The top quark can, for all practical purposes, be treated as a heavy particle, not a parton. Collins has shown that, to all orders of the perturbation expansion, factorization of the DIS structure functions holds for massive partons to the same degree of rigor as in the familiar zero-mass case. Collins98  The relatively simple physical ideas underlying the full factorization proof are described in Refs. ASTW ; TKS ; acotchi . To implement this formalism, which we shall refer to as the generalized $`\overline{MS}`$ scheme (or GM-VFNS, as mentioned earlier), we follow the explicit procedures outlined in the Appendix of Ref. ASTW , and supplement this with the threshold prescription of Refs. TKS ; acotchi . The latter reference also contains a review of alternate approaches.

The NLO DIS structure functions are given by the generic formula (suppressing the structure function label 1,2,3),

$$F(x,Q)=\underset{a}{}\frac{dz}{z}f_a(z,\mu )\widehat{\omega }^a(\frac{x}{z},\frac{Q}{\mu },\frac{m_H}{\mu },\alpha _s(\mu ))+𝒪(\alpha _s^2,\frac{\mathrm{\Lambda }^2}{Q^2},\frac{\mathrm{\Lambda }^2}{m_H^2})$$
(1)
where “$`a`$” is the initial state parton label, $`f_a`$ is the parton distribution function, $`\widehat{\omega }^a`$ is the hard-scattering cross-section calculated in PQCD to order $`\alpha _s^1`$, $`\mathrm{\Lambda }`$ is the QCD-lambda parameter, and $`m_H`$ (generically) represents heavy quark masses, if present. The prescription-dependence allowed by the PQCD formalism is associated with possible implementations of the first term on the right-hand side of Eq. 1, within the accuracy specified by the remainder term (which will be dropped from now on).444Note, according to the factorization proof of Ref. Collins98 , the remainder term contains corrections of order $`\mathrm{\Lambda }^2/Q^2`$ and $`\mathrm{\Lambda }^2/m_H^2`$, but no corrections of order $`m_H^2/Q^2`$. In the sum over partons in Eq. 1, contributions due to the light quark ($`u,d,s`$) are standard; contributions arising from massive charm and bottom quarks, along with that of the gluon, are treated differently, as explained below.

### II.1 Leading contributions and the ACOT($`\chi `$) implementation

To be specific, we will consider DIS neutral-current charm production process. The leading-order process is a virtual-photon scattering off a charm parton: $`\gamma ^{}cc`$. The NLO QCD corrections for this process are the boson-gluon fusion process ($`\gamma ^{}gc\overline{c}`$), and the gluon radiation process ($`\gamma ^{}ccg`$) as well as the charm-initated one-loop virtual process ($`\gamma ^{}cc`$).

First, let us examine the kinematic region near the charm production threshold: $`W=\sqrt{Q^2(1/x1)}2m_c`$. To make the notation precise, we will use the renormalization scale $`\mu `$ to generically label our characteristic energy scale; separately, we use $`\mu _m`$ to label the matching point, and $`\mu _t`$ to label the transition point (cf. ASTW  and below). To implement the generalized $`\overline{MS}`$ scheme, the first step is to consider a 3-flavor scheme (appropriate for energy scales $`\mu \stackrel{<}{}m_c`$), a 4-flavor scheme (for scales $`\mu \stackrel{>}{}m_c`$), and choose a matching point $`\mu _m`$ (on the order of $`m_c`$) where the two schemes are matched (i.e., where the discontinuities of $`\alpha _s(\mu )`$ and $`f_a(x,\mu )`$ are calculated). We choose the conventional value for the matching point, $`\mu _m=m_c`$, which yields the simplest matching conditions. ColTun

Having precisely defined the 3-flavor and 4-flavor schemes, which co-exist in the region of the charm threshold, we can still choose a *transition point*, $`\mu _t`$, where one makes the transition from the 3-flavor to the 4-flavor scheme in the calculation of physical quantities.555The distinction between the matching and the transition points is not addressed in most papers. This distinction is made in Ref. Collins98 ; it is discussed at length in the Appendix of Ref. ASTW . As mentioned above, arguably, there is a good case for choosing the transition point $`\mu _t`$ to be higher than the simple matching point $`\mu _m=m_c`$. The choice of the transition point is arbitrary—its choice is part of the definition of the composite renormalization scheme. While there is much flexibility in the choice of $`\mu _t`$, clearly it must lie within the overlapping region of applicability of the 3-flavor and 4-flavor schemes—therefore, we typically choose $`\mu _tm_c`$. Although a case can be made, in principle, for choosing $`\mu _t`$ far above $`m_c`$ (since the 4-flavor scheme is not a natural scheme just above $`m_c`$), or for considering an $`x`$-dependent transition scale, $`\mu _t(x,m_c)`$ (because the threshold $`W>2m_c`$ depends on $`x`$ as well as $`Q`$), we make the plain choice $`\mu _t=m_c`$ to simplify the calculation. This choice is reasonable only within a prescription for handling $`x`$-dependent threshold effects of Eq. 1 which naturally suppresses the 4-flavor contribution to the physical structure function when $`\mu m_c`$ and/or $`W\stackrel{>}{}\mathrm{\hspace{0.17em}2}m_c`$. The specific renormalization prescription of Ref. TKS  has precisely this feature, as we shall outline below.

Detailed descriptions of the various NLO contributions to the the generalized $`\overline{MS}`$ scheme were discussed in Ref. ASTW . While we refrain from reproducing them here, we do need to specify the precise way the threshold effects are implemented when heavy quarks are involved. Consider first the numerically significant NLO *boson-gluon fusion* contribution ($`\gamma ^{}gc\overline{c}`$):

$`F^{(1)}(\gamma ^{}gc\overline{c})`$ $`=`$ $`{\displaystyle \frac{dz}{z}g(z,\mu )\widehat{\omega }_g(\frac{x}{z},\frac{Q}{\mu },\frac{m_c}{\mu })}`$
$`=`$ $`\alpha _s(\mu )\left[{\displaystyle _\chi ^1}{\displaystyle \frac{dz}{z}}g(z,\mu )\omega _g^1({\displaystyle \frac{x}{z}},{\displaystyle \frac{m_c}{Q}})\mathrm{ln}\left({\displaystyle \frac{\mu }{m_c}}\right){\displaystyle _\zeta ^1}{\displaystyle \frac{dz}{z}}g(z,\mu )P_{gc}\left({\displaystyle \frac{\zeta }{z}}\right)\omega _c^0\left({\displaystyle \frac{m_c}{Q}}\right)\right]`$
where the first term on the right-hand side corresponds to the unsubtracted boson-gluon-fusion diagram contribution (with the $`\alpha _s`$ factor explicitly taken out),666We use the notation that, for each flavor, $`\omega _a`$ denotes the parton-level cross-section before infrared and collinear subtraction, while $`\widehat{\omega }_a`$ denotes the corresponding infra-red safe cross-section after subtraction. and the second term represents the subtraction term which renders $`\widehat{\omega }_g`$ infra-red safe. The variable $`\chi `$ is dictated by the kinematics of the boson-gluon-fusion partonic subprocess ($`\gamma ^{}gc\overline{c}`$) to be a generalized scaling (or “rescaling”) variable

$$\chi =x\left(1+\frac{4m_c^2}{Q^2}\right)$$
(3)
where $`x`$ is the conventional Bjorken variable. In contrast, the other scaling variable, $`\zeta `$, in Eq. II.1 can, in principle, be chosen arbitrarily provided the following two constraints are imposed.

1. $`\zeta x`$ in the high energy limit.
2. The same $`\zeta `$ variable is used here as in the order $`\alpha _s^0`$ (“leading-order”) process $`\gamma ^{}cc`$.

In the high energy limit ($`m_c^2/Q^20`$), the first constraint ensures that both $`\zeta x`$ and $`\chi x`$ such that the complete boson-gluon fusion contribution is infra-red safe. The leading-order process $`\gamma ^{}cc`$ is:

$$F^{(0)}(\gamma ^{}cc)=\frac{dz}{z}c(z,\mu )\widehat{\omega }_c^0(\frac{Q}{\mu },\frac{m_c}{\mu })\delta \left(1\frac{\zeta }{z}\right)=c(\zeta ,\mu )\widehat{\omega }_c^0(\frac{Q}{\mu },\frac{m_c}{\mu }).$$
(4)
For the second constraint, if the same $`\zeta `$ variable is used for both the boson-gluon fusion process of Eq. (II.1) ($`\gamma ^{}gc\overline{c}`$) and the leading-order ($`\gamma ^{}cc`$) contribution of Eq. (4), then the gluon subtraction term (last term in Eq. II.1) cancels the leading-order contribution in the limit $`\mu m_c`$; this ensures that the 4-flavor formula will reduce to the 3-flavor formula, and the dominant contribution in this region will be the unsubtracted gluon fusion term, as required by the correct physics to this order of PQCD.

A natural choice of $`\zeta `$ that ensures the best behavior at *both* the high and the low energy limits is proposed by Ref. TKS :

$$\zeta =\chi x(1+4m_c^2/Q^2).$$
(5)
This choice has the property that as the heavy quark production threshold is approached from above ($`W2m_c`$), then $`\zeta 1`$ such that the integration measure in Eq. II.1 vanishes.777This is more easily seen from the equivalent formula $`\chi =1(1x)(14m_c^2/W^2)`$. Therefore, with the choice $`\zeta =\chi `$, the subtraction term in Eq. II.1 and the LO quark term, Eq. 4, individually vanish at threshold by kinematics. In addition, the cancellation between these two terms due to the dynamics of PQCD, as explained in Ref. ACOT , still operates.

This results in unequaled stability of the predictions in the threshold region, compared to previous implementations of the GM-VNFS in the threshold region ACOT ; ThorneRob ; ChuvakinEtal ; NasonEtal . These usually give rise to rather prescription-sensitive behavior just above threshold because the perturbation-expansion-inspired prescriptions generally overlook heavy quark production kinematic requirements, while the physically important gluon fusion term (in the threshold region) is strongly kinematically-suppressed. The robust implementation we adopt will be referred to as the “ACOT($`\chi `$)” scheme, following Refs. TKS ; acotchi .888So named since it is a variant of the original ACOT approach ACOT , replacing the naive $`x`$ with the kinematically natural rescaling variable $`\chi `$ of Eq. 3.

### II.2 Other contributions to the full NLO calculation

For the complete order $`\alpha _s`$ calculation, we must also include the real and virtual NLO quark-initiated contribution to Eq. 1,

$$F^{(1)}(\gamma ^{}cc+X)=\frac{dz}{z}c(z,\mu )\widehat{\omega }_c^1(\frac{x}{z},\frac{Q}{\mu },\frac{m_c}{\mu }).$$
(6)
The hard-scattering cross-section, $`\widehat{\omega }_c^1(x,\frac{Q}{\mu },\frac{m_c}{\mu })`$, has been computed for general masses nloquark  in the ACOT scheme ACOT . For scales around the canonical DIS choice of $`\mu Q`$, this term is numerically insignificant compared to the terms considered above. We include it in our calculation to ensure (perturbative) consistency.

It is also worth mentioning that the calculation of this term, which is computation-intensive in spite of the smallness of its numerical contribution, can be considerably streamlined by using a simplifying feature of PQCD in the generalized $`\overline{MS}`$ scheme: one can set $`m_c=0`$ in $`\widehat{\omega }_c^1`$ and $`\omega _c^0`$ in Eq. (II.1) and Eq. (4) without sacrificing accuracy. This point is discussed in detail in Ref. SACOT  where the numerical insensitivity of the calculation to the choice of $`m_c`$ is explicitly demonstrated.999The simplified scheme SACOT  can also be recovered by investigation of the resummed perturbation series. acotchi

### II.3 Remaining scheme and scale choices

The only remaining arbitrariness of our prescription is now associated with the familiar choice of the factorization (and renormalization) scale $`\mu `$. The scale dependence of perturbative calculations in general arise from imperfect matching between the terms in the truncated perturbation series. Because of the near-perfect matching in the ACOT($`\chi `$) prescription discussed above (in contrast to alternate choices), the dependence of the physical structure functions on the choice of $`\mu `$ is completely negligible (cf., Ref. TKS ) throughout the energy range of our global analysis. For definiteness, we use $`\mu ^2=Q^2+m_c^2`$. In summary, the quark initiated terms in the simplified ACOT($`\chi `$) scheme adopted in this work are:

$$_\chi ^1\frac{dz}{z}c(z,\mu )\widehat{\omega }_c^{0,1}\left(\frac{\chi }{z},\frac{Q}{\mu },\frac{m_c}{\mu }=0\right)|_{\mu ^2=Q^2+m_c^2}=_\chi ^1\frac{dz}{z}c(z,\mu )C_{\overline{\mathrm{MS}}}^{0,1}(\frac{\chi }{z},\frac{Q}{\mu })|_{\mu ^2=Q^2+m_c^2}$$
(7)
with $`C_{\overline{\mathrm{MS}}}^{0,1}`$ representing the standard coefficient functions for massless quarks in LO and NLO, defined in the conventional $`\overline{MS}`$ scheme.

The same treatment is applied to the $`b`$-quark threshold region. Here, the 4-flavor scheme is matched onto the 5-flavor scheme for scales larger than $`m_b`$, and the production of open bottom particles in the final state is added to the inclusive DIS structure functions. In principle, a similar treatment needs to be applied for other processes included in the global analysis (currently, Drell-Yan and jet production) using appropriate hard-scattering cross-sections calculated in the massive scheme. In practice, this is not done for two reasons: (i) NLO calculations of the hard-scattering cross-section involving massive quarks do not yet exist for these processes; and (ii) unlike the case of DIS, the experimental errors in these processes are relatively large compared to the anticipated differences between the massless and massive scheme calculations.

## III Global Fitting and the CTEQ6HQ Parton Distributions

The new global fitting performed in the generalized $`\overline{MS}`$ formalism follows the same procedure as that of the earlier CTEQ6 analysis cteq6  using the zero-mass parton formalism. The data sets used before are supplemented by: (i) the H1 $`F_2^{e^+p}`$ set H1positron99  (which was inadvertently left out in Ref. cteq6 ); and (ii) the H1 H1f2c01  and ZEUS Zeusf2c99  data sets for the structure function $`F_2^c`$ with tagged charm particles in the final state. The additional H1 $`F_2^{e^+p}`$ data set does not have much influence on the new analysis, since one already obtains an excellent fit to these data just by comparing them with predictions of CTEQ6M. The $`F_2^c`$ data sets are quite relevant for this analysis since $`F_2^c`$ is sensitive to the charm and gluon distributions, which are tightly coupled in the generalized $`\overline{MS}`$ formalism.101010These data sets were not used in the CTEQ6M analysis because $`F_2^c`$ is not well-defined theoretically in the the zero-mass parton formalism.

The parametrization of the non-perturbative parton distribution functions at $`Q_0=m_c=1.3`$ GeV is the same as in Ref. cteq6 . For this study, we assume as usual that charm partons are entirely “radiatively generated” (i.e., through QCD evolution) from the starting scale $`Q_0`$ onward. This assumption is somewhat arbitrary, and it is obviously dependent on the choice of $`Q_0`$. The possibility for having a small component of non-perturbative charm incc  at low $`Q`$, and its physical consequence will be examined separately.

As in the previous CTEQ6 analysis, correlated experimental systematic errors are fully incorporated whenever available. The best fit obtained with these inputs shall be called CTEQ6HQ, or C6HQ for short. A broad measure of the quality of this fit is provided by the overall $`\chi ^2`$ of 2033 for a total number of 1950 data points ($`\chi ^2`$/DOF = 1.04). This is to be compared to a $`\chi ^2`$ of 1946 for 1811 points ($`\chi ^2`$/DOF = 1.07) in the case of CTEQ6M (abbreviated to C6M in the following) cteq6 . To gain a better feel of how these fits compare, we show in Table 1 (the 3rd and 4th columns) the overall $`\chi ^2`$ values and, in parentheses, $`\chi ^2`$ per data point, as well as for the applicable individual data sets.

The total number of data points in this head-to-head comparison (not including the charm production data points) is 1925. The new C6HQ fit reduces the overall $`\chi ^2`$ by 29 out of $``$2000 as compared to the C6M fit. The improvement of this generalized $`\overline{MS}`$ result over the zero-mass $`\overline{MS}`$ result is encouraging, since the generalized $`\overline{MS}`$ formalism represents a more accurate formulation of PQCD. However, a difference of $`\chi ^2`$ of 29 is within the current estimated range of uncertainty of PDF analysis. cteq6 ; Hesse ; Lagrange  Therefore, the significance of this difference is arguable. We also note that the improvement in $`\chi ^2`$ is spread over most of the data sets: there is no smoking gun for the overall difference.

The last two columns of Table 1 compare the above results with two possible uses of the PDFs that represent a *mis-use* of PQCD in principle, but occur frequently in the literature in practice, perhaps out of necessity. These involve using PDFs obtained in the general-mass scheme convoluted with hard-scattering cross-sections (Wilson coefficients) defined in the zero-mass scheme, and vice versa.111111For instance, the MRST distributions MRST  are obtained in the General-Mass formalism (using the Robert-Thorne implementation ThorneRob ); they are often used in applications convolving with readily available Zero-Mass hard-scattering cross-sections. For the same data sets, these mis-matched schemes result in a difference of 420$``$490 in the overall $`\chi ^2`$. These are quite large differences relative to the tolerances discussed in Refs. cteq6 ; MRST , and results in clear discrepancies with some of the precision DIS data sets, as can be seen in Table 1. The lesson is clear: *for quantitative applications, it is imperative to maintain consistency between the PDFs and the hard-scattering cross-sections*.

## IV Comparison with data

In Ref. cteq6 , we presented an extensive comparison between the CTEQ6 results and the data sets used in the global analysis, including new ways to explicitly account for the correlated experimental systematic errors. Since the new C6HQ fit is generally similar to the C6M fit, we shall not duplicate the same comparisons for those quantities where the differences are minimal. In subsection IV.1, we compare with the DIS charm production data, which were not used in the previous CTEQ6 analysis. In subsection IV.2, we discuss some implications of the neutrino data.

### IV.1 DIS Charm Production Data from HERA

Because the H1 and ZEUS $`F_2^c`$ structure function are sensitive to the charm and gluon distributions, these will play an important role in extracting the C6HQ PDFs. In Fig. 1, we display comparisons of the charm production data from H1 H1f2c01  and ZEUS ZeusF2c  with the C6HQ theory value. For simplicity, we combine all the data points from each experiment in a single plot, and we scale the plot by the theoretical calculation. On the x-axis, the data points are ordered using $`x`$ and $`Q`$ as the primary and secondary sorting variables, respectively. The error bars represent the statistical and uncorrelated systematic errors added in quadrature. We see that the fits are good with $`\chi ^2`$ per data point of 0.881 (7/8) for H1, and 0.98 (18/18) for ZEUS, respectively.

### IV.2 Neutrino Data in Global Parton Analyses

The NuTeV measurement of the weak mixing angle $`\mathrm{sin}^2\theta _\mathrm{W}`$ (Ref. nuanom ) has recently focused considerable attention on the neutrino induced DIS process. nuint  In this measurement, the weak mixing angle $`\mathrm{sin}^2\theta _\mathrm{W}`$ is extracted from the ratio of neutral current (NC) to charged current (CC) $`\nu `$ and $`\overline{\nu }`$ cross-sections. Even before this $`\mathrm{sin}^2\theta _\mathrm{W}`$ measurement, there have been other long-standing unresolved issues in comparing CC structure functions $`F_2^{CC}`$ (measured in $`\nu `$ and $`\overline{\nu }`$ scattering) with NC structure functions $`F_2^{NC}`$ (measured in $`e^\pm `$ and $`\mu ^\pm `$ scattering) at modestly low-$`x`$ ($`10^2`$). One manifestation of this is that the recently measured structure function $`\mathrm{\Delta }xF_3^{\nu /\overline{\nu }N}`$ is not compatible with the QCD predictions. kosty ; dimuons  As this structure function is particularly sensitive to the heavy quark components, $`\mathrm{\Delta }xF_3^{\nu /\overline{\nu }N}s(x)c(x)`$, it is certainly important to properly treat the heavy quark mass. While the deviations of $`\mathrm{\Delta }xF_3^{\nu /\overline{\nu }N}`$ are at a non-dramatic $`1\sigma `$ level, the pattern has been systematic and persistent.

Therefore, we re-examine the influence of the neutrino structure functions in the current global analysis keeping in mind the following issues. First, the high statistics neutrino structure function data of CCFR used here are from the newer “physics model independent” data analysis pmi , rather than the earlier ones which contained model-dependent corrections. ccfrsf  Secondly, meaningful (i.e., quantitative) comparison of NC and CC structure functions (which represent different hard-scattering processes) can be made only within a common underlying theoretical framework with an accuracy comparable to the experimental uncertainties.121212For example, direct comparisons of the ratio $`F_2^\nu /F_2^\mu `$ to theory (suggested by the “5/18-rule”) are meaningful only at the LO parton model level, which is not appropriate for the precision data from current experiments. The proper treatment of charm mass effects, such as those described in the current study, are an important part of that theoretical formalism, cf. nuc .

First, we ask whether the neutrino structure function, $`F_2^{\nu /\overline{\nu }N}`$, obtained by the newer model-independent analysis is consistent with the muon-induced NC fixed-target structure functions, $`F_2^{\mu ^\pm N}`$, and the HERA collider structure functions, $`F_2^{e^\pm N}`$, in a consistent NLO analysis.131313Ref. pmi  suggests these data are consistent by examining the ratio of the structure functions. However, as explained in the text, it is not appropriate to compare the ratio of two different structure functions with theory in QCD beyond LO; a complete analysis involving a NLO fit is required. Fig. 2 shows the comparison between the data with the CTEQ6HQ fit (using $`m_c=1.3\mathrm{GeV}`$).

We see that, even if this data set is included in the global fit (cf., Table 1), the measured values in the low-$`x`$ bins (left figure) consistently lie above the theory curves. The systematic deviation of data points from the global fit is reflected in the relatively large $`\chi ^2`$’s associated with this data set, as seen in Table 1. The fit is dominated by the much more extensive NC data (from the BCDMS, NMC, H1, and ZEUS experiments), hence one can regard the theory curves as representing the NC data properly “corrected” for NLO QCD effects. The pattern of deviation of the neutrino data from the average NC predictions at low-$`x`$ is similar to earlier comparisons. In other words, the discrepancy between CCFR $`F_2^\nu `$ and the other NC $`F_2^{e,\mu }`$ measurements at small $`x`$ persists in the context of contemporary global QCD analyses, even when charm mass effects are properly treated as in the present analysis. Other possible theoretical sources for this difference (in addition to the NLO contributions included in this comparison) have been discussed in Refs. kosty ; nuint ; nuref , and in the literature quoted therein. On the experimental side, preliminary NuTeV measurements bernstein  seem to yield improved agreements with the NLO QCD predictions at low-$`x`$. It remains to be seen whether the final NuTeV results will resolve this problem.

In our analysis, the strangeness distribution is constrained to be proportional to $`[\overline{u}(x)+\overline{d}(x)]`$; and its normalization at the scale $`Q_0`$ is constrained to be

$$\kappa =\frac{𝑑xx[s(x)+\overline{s}(x)]}{𝑑xx[\overline{u}(x)+\overline{d}(x)]}\frac{1}{2},$$
(8)
as in most current global analysis. In the future, a NLO analysis of the CC charm production data ($`\nu sc\mu \mu ^\pm \mu ^{}X`$) from CCFR and NuTeV charmcc  has the potential to determine $`\{s(x),\overline{s}(x)\}`$ with much more precision. Then, a re-assessment of the apparent discrepancy between the measured $`\mathrm{\Delta }xF_3^{\nu /\overline{\nu }N}`$ and the low-lying QCD predictions kosty  will be warranted. An accurate determination of the strangeness of the proton, as well as the strangeness asymmetry $`[s(x)\overline{s}(x)]`$ will have important implications for these measurements, including the NuTeV anomaly, cf., Refs. nuanom ; nuint .

Clearly, the interpretation of the neutrino DIS data requires a detailed NLO analysis. Global QCD analyses of parton distributions, such as the one performed in this paper, are essential to unraveling the underlying physics. Progress on both the experimental and the theoretical fronts are needed; this advancement is ongoing.

## V Comparison with related PDFS: CTEQ6HQ, CTEQ6M, CTEQ5HQ

The CTEQ6HQ and CTEQ6M fits provide comparable descriptions of the global QCD data in two different schemes. Some of the differences in the PDFs arise purely from the choice of scheme. We are particularly interested in the differences for the charm distribution, and the closely correlated gluon distribution, due to the improved treatment of heavy quark effects in the generalized $`\overline{MS}`$ scheme. It is also interesting to compare the differences between the earlier (previous generation) CTEQ5HQ (C5HQ) distributions with the new CTEQ6HQ distributions; differences between these PDFs are attributable both to new data, and to minor differences in the way the theoretical inputs are implemented. We combine these comparisons in the figures presented below.

### V.1 Light partons at the input scale

Figs. 3a,b show the $`u`$-quark and $`d`$-quark distributions from C5HQ, C6M and C6HQ at the initial $`Q`$-scale of 1.3 GeV (the charm mass). In order to exhibit the behavior of the PDFs at both large and small $`x`$ values, the x-axis is scaled linearly in $`\sqrt{x}`$. To give some additional physical insight to these plots, the y-axis is chosen to be $`x^{1.5}f(x,Q)`$, so that the area under each curve represents the momentum fraction carried by that distribution. This chosen scaling then focuses on the valence peak, and suppresses the small-$`x`$ sea excitations. We see that the valence $`u`$-quark and $`d`$-quark distributions are well constrained by the precision DIS data; the scheme dependence is not pronounced.

Figs. 4a,b show the same PDFs, except the axes are scaled differently in order to highlight the sea quark distributions. The differences between the CTEQ5 and CTEQ6 generations of PDFs are a result of the improved low-$`x`$ data. When we compare the HQ and $`\overline{MS}`$ type PDFs of the CTEQ6 generation, we observe the general pattern that the HQ distributions overshoot the $`\overline{MS}`$ ones at low $`x`$ ($`10^3`$), and that the distributions meet again at ultra-small $`x`$ ($`\stackrel{<}{}\mathrm{\hspace{0.17em}10}^5`$). We can qualitatively understand the larger sea quark distributions for the HQ scheme since these distributions are compensating for the suppression of charm in this scheme; this suppression is absent in the zero-mass $`\overline{MS}`$ scheme. In the region of ultra-small $`x`$, there is no data constraining these distributions; we show this effect graphically in the next Section where we consider the PDF uncertainty band of the the HQ and $`\overline{MS}`$ distributions.

Fig. 5 shows the comparison of the gluon distribution at $`Q_0`$. Here the difference between C5HQ and the CTEQ6 generation of gluon distributions is pronounced. As discussed in Ref. cteq6 , the change in this least-well-determined parton distribution is due to the recent precision DIS data (most influential in the small $`x`$ region) in conjunction with the greatly improved inclusive jet data from the Tevatron (critical for the medium to large $`x`$ regions). The differences between C6HQ and C6M gluons at large $`x`$ are due to a combination of scheme-dependence, and the inherent uncertainty range of the current analysis. cteq6

Fig. 6 shows the comparison of the strange distributions at the same $`Q_0`$. The noticeable difference between the C5HQ curve and the others, in this case, is largely the result of different theoretical inputs: namely, the $`\kappa `$ parameter (cf., Eq. 8). $`\kappa `$ determines the ratio of strange to non-strange sea quarks at the initial scale $`Q_0`$. This $`\kappa `$ factor, known only approximately, was chosen to be $`1/2`$ both in the CTEQ5 and CTEQ6 analyses, but for slightly different values of $`Q_0`$—1.0 GeV for C5HQ and 1.3 GeV for CTEQ6 sets.

### V.2 Charm and gluon distributions at $`Q^2=10`$ GeV2

To compare the differences of the charm distributions, we need to move some distance above the charm mass scale, since all current sets assume a zero charm distribution at $`\mu =m_c`$. Fig. 7a makes this comparison at $`Q^2=10`$ GeV2. We see a substantial difference between the C5HQ and the two CTEQ6 charm distributions. Unlike the case of the strange quark, this difference is physical; it is mainly a reflection of the difference in the gluon distributions (shown already for a lower $`Q`$ value in Fig. 5) because charm is radiatively generated from the gluon. To confirm this hypothesis, Fig. 7b shows the gluon distributions at the same $`Q^2=10`$ GeV2; the similarity is clear.

### V.3 $`Q`$-dependence of heavy quark mass effects

In this subsection, we compare the $`Q`$ dependence of the C6M and C6HQ parton distributions at fixed values of $`x`$ in order to see how the differences between PDFs in the zero-mass and general-mass schemes vary with increasing $`\mu `$-scale. We expect these differences (relative to the general PDF uncertainties) will decrease with increasing $`Q`$. The purpose of this study is two-fold. First, we check the self-consistency of the analysis, which is based on the expectation that power suppressed mass terms vanish asymptotically as $`Q\mathrm{}`$. Secondly, we determine the specific scales where the mass-effects become insignificant, in practice. The latter question is relevant for high-energy applications of PDFs derived from lower-energy data.141414It also serves to validate the self-consistency of analyses where threshold effects have been neglected (such as for collider processes). It can also be seen as a general illustrative example demonstrating how the relative importance of power suppressed terms decrease in the high $`Q^2`$ region where the conventional zero-mass PQCD results becomes dominant.

To examine the effects due to heavy flavor masses in the theoretical formalism, we cannot directly compare physical quantities that are being fit—they are the same (within the errors of the global fit) by construction. On the other hand, parton distribution function of individual flavors are not a good gauge either, because they are not physical—they are scheme dependent. Thus, for this comparison, we compromise and examine the quark singlet combination ($`\mathrm{\Sigma }_q(q+\overline{q})`$) and the gluon distribution.

Figs. 8a,b show, as solid lines, the C6HQ gluon and quark singlet distributions normalized to that of C6M. The shaded bands represent the ranges of uncertainty of the relevant quantities due to experimental sources, as estimated in the CTEQ6 global analysis. cteq6  We observe the following pattern: for most of the parameter space, the deviation of the C6HQ from C6M falls within the estimated uncertainty range. In particular, for the gluon distribution this is true for all values of $`x`$. For the quark distribution, the HQ fit overshoots the zero-mass fit significantly (i.e., beyond the uncertainty band) at moderately low $`x`$ $`[10^2,10^3]`$. At very low $`x`$ $`10^5`$, the differences fall within the uncertainty bands again (primarily due to the increasing uncertainty). Although the HQ quark PDFs at $`x[10^2,10^3]`$ eventually evolve into the error band of the zero-mass fit as $`Q`$ increases, one may wonder why the two fits do not merge faster. The GM-VFNS and ZM-VFNS differ by the manner in which they organize terms of order $`𝒪(m^2/Q^2)`$. Thus, the Wilson coefficients in the two schemes will be comparable when $`(m^2/Q^2)<<1`$. The boundary conditions of the PDFs at $`Q_0`$ will also differ between the two schemes by terms of order $`𝒪(m^2/Q^2)`$; however, since the PDFs evolve logarithmically in the scaling variable, while the power corrections fall off more quickly, this can contribute to scheme differences at scales significantly above the scale $`m`$. This effect is evident in the case of the singlet quark PDF at moderately low $`x`$ values. Thus, the observed behavior is quite natural.

## VI Concluding Remarks

The CTEQ6HQ PDFs presented here complement the previously published CTEQ6 sets by providing distributions which can be used in the generalized MS-bar scheme with non-zero mass partons. This analysis includes the complete set of NLO processes including the real and virtual quark-initiated terms. Additionally, the ACOT($`\chi `$) scheme is used to introduce a generalized scaling variable which provides numerically stable results for the entire energy range—from heavy quark thresholds to the high energy limit.

While the zero-mass parton scheme is sufficient for many purposes, the fully massive scheme can be important when physical quantities are sufficiently sensitive to heavy quark contributions. This is evident when comparing the CTEQ6HQ and CTEQ6M fits to the mis-matched sets (Table 1) where the precise DIS data from HERA highlights the discrepancies.

The CTEQ6HQ fits also provide the basis for a series of further studies involving more quantitative analysis of strange, charm, and bottom quark distributions inside the nucleon. For example, the CTEQ6HQ PDFs are necessary for a consistent analysis of resummed differential distributions for heavy quark production such as in Ref. Nadolsky:2002jr . Using the full range of data from both the charged and neutral current processes, these distributions can reduce the uncertainties in the calculations; hence, they have significant implications for charm and bottom production, and can help resolve questions about intrinsic heavy quark constituents inside the proton, the $`\mathrm{\Delta }xF_3`$ structure function, and the extraction of $`\mathrm{sin}\theta _W`$.

## Acknowledgment

We thank our colleagues J. Huston, P. Nadolsky, J. Pumplin, and D. Stump, for fruitful collaboration on the CTEQ6 project that forms the foundation of this study. F.O. acknowledges the hospitality of MSU and BNL where a portion of this work was performed. S.K. is grateful to RIKEN, Brookhaven National Laboratory and the U.S. Department of Energy (contract No. DE-AC02-98CH10886) for providing the facilities essential for the completion of this work. This research was supported by the National Science Foundation (grant No. 0100677), and by the Lightner-Sams Foundation.


text:
# Amplified imitation in percolation model of stock market

## 1 Introduction

Over the last decade, artificial economic world living in computer made markets of financial assets has become an interest for physicists experienced in complex systems research. See  for review. The reasons for the Cont-Bouchaud proposition is to explain effects of herd behavior which is observable among real traders. Differently from models of markets with agents rational, or heterogeneous, or adaptive, the Cont-Bouchaud model (CB-model, in short) considers noisy traders only, in the sense that traders do not have any rational strategy. Instead the strong communication and information links between traders is proposed. In effect agents form groups- called clusters, within which behavior of agants is common.

Complex system approach provides percolation as a powerful tool in study of connected objects of different type and origin, . The Cont-Bouchaud proposition is to adapt the percolation phenomena to design a stock market organization . The original model uses infinite-range bond percolation. Similar interactions are proposed by considering nearest neighbor and site percolation instead  At a very few assumptions, the model reproduces, so-called stylized, facts known about fluctuations of price of real market assets which are widely accepted, . In particular, the price fluctuations measured by so-called returns, i.e., change of the market price, (i) are not correlated, (ii) average of them is zero and (iii) tails of a distribution of returns decay accordingly to the power-law. These features are achieved only due to exploration of critical properties attributed to the percolation phenomena.

In the CB-model there is no other interaction between investors than instantaneous imitation through a cluster structure (these structure can change in time ). Our proposition is to consider imitation of one cluster by other clusters but taken in the $`next`$ time step. We say the cluster imitates some other cluster with a time step delay. In this sense the herd behavior of the CB-model is amplified by the delay imitation.

The additional interaction introduced seems to be particularly suitable when poorly developed stock markets are considered. In the result of multi century experience the large markets, like London City or New York Stock Exchange, have developed the solid structures and their regulations are firmly executed. One can say that the fundamental assumption of Efficient Market Hypothesis  that the whole information is available for all investors at the same instant, is close to be satisfied. The race and competition between individual investors effects in that the same information is differently understood and in consequence differently transferred into market decisions. There is no place for imitation between investors for at least two reasons. First, the next piece of information arrives at the market and investors are forced to concentrate on this new subject. Second, investors are known to be extremely strong individuals in the sense it is said, they are overconfident about their personal understanding of the market. Such attitude leads to the strong disagreement between investors rather than cooperation. However when a market is far from equilibrium, because, for example, a market is strongly influenced by political decisions or liquidity of a stock is limited then the only strategy to win is to search for ‘well informed’ investors to join their team or if it is impossible to follow their decisions.

In the subsequent sections we provide: setting of the model (Section 2), analysis of results obtained in simulations of the model (Section 3), examples of stocks data from the Warsaw Stock Exchange market for which the hypothesis of amplified imitation as underlying mechanism of driving the price can be proposed (Section 4). In the closing section further possible modifications to the model are proposed.

## 2 Details of the model

In the Stauffer et al. representation  of the CB-model one deals with a number of investors who are distributed on a square lattice of a linear size $`L`$. Each of the investors is assigned to a site at random with the probability $`p`$. A site can be occupied by one investor or left empty. The geometric relations arising from nearest neighboring occupied sites, called clusters, are interpreted as coalitions of investors- groups of traders that act together on a stock market. The coalitions last for a limited period of time after which a new arrangement of stock market participants is introduced. The clusters take decisions randomly and independently of each other to $`buy`$ or $`sell`$, with the probability $`a`$, or just $`wait`$ with a probability $`12a`$. The resulting stock price is proportional to the sum of demand and sell orders from all clusters which have been active during a time step. The parameter $`a`$ is called $`activity`$ and can be seen as a control of the length of a time step . If activity is high, which means that $`a`$ is close to $`0.5`$, then all clusters are active during a session. Such activity characterizes a long time step. The independent and stochastic decisions taken during long sessions lead to the histogram of price changes of a Gaussian type. If the activity is low, $`a<\frac{1}{L^2}`$, what means that a time step is short, then on average one cluster acts during a given market session. In effect the histogram of price changes reflects properties of the cluster size distribution. If $`p`$ is close to the percolation threshold then the power-law decay of tails of the distribution is observed with the exponent $`\tau =2`$ in case of a square lattice .

Let investors occupy randomly a fraction $`p_c=0.592746`$ of the sites of $`L\times L`$ square lattice. This particular value of the parameter $`p_c`$ provides critical properties for a cluster structure on a square lattice: a cluster that spans the lattice occurs with the probability $`\frac{1}{2}`$ . In the following we use the Hoshen-Kopelman algorithm  to identify clusters.

The clusters of size greater than 2 are split into two classes:
— class A contains clusters of well informed agents. These clusters take decisions independently and at random: they buy or sell at the probability $`a`$ and at probability $`12a`$ they do nothing.
— class B contains clusters which imitate other clusters. Each cluster is linked to some other cluster. Imitation means that the next time step decision of a B-class cluster is exactly the same as the present decision of the cluster to which the B-class cluster is linked.
Let us assume that all clusters of size smaller than 3 run the ordinary CB-model dynamics.

The following two types of imitation will be considered: model A \- each cluster from the class B at random chooses a cluster from the class A to imitate, and model AB \- each cluster from the class B is linked to a randomly choosen cluster from any class. So that the pool of clusters is arranged into trees, see Fig. 1. The root of each tree is a leader- an independently acting participant of a market. The rest of clusters grupped in a tree is either of one step depth tree — case of the model A, or many steps depth tree — case of the model AB. Clusters from a tree take the same market decision as their leader but with a delay which equals to the path length to the leader of the tree.

Let us introduce the parameter $`b(0,1]`$ the probability that a cluster of size greater than 2 belongs to the class A. At a given $`b`$ the average number of clusters in one tree is

$$1+\frac{(1b)N}{bN}=\frac{1}{b}$$

where $`N`$ is the average number of clusters of size greater than 2. Hence if $`b=1`$ then both models collapse to the original CB-model.

## 3 Results

Our simulations are performed on a square lattice of the linear size $`L=200`$. On this lattice at $`p_c`$ the average number of clusters of size greater than 2 is $`391\pm 25`$ (total number of all clusters of 1’s is $`1205\pm 50`$). Setting $`b=0.001`$ indicates that at most one tree is present on the lattice while $`b=0.005`$ ($`b=0.01`$) allows two (four, respectively) trees to exist on the lattice. If $`b=0.1`$ then at average about $`40`$ trees take part in the market game. Therefore, very long runs must be considered to obtain satisfactory statistics. Our simulations are performed for (number of lattice rearrangements) $``$ (length of each of arrangement) = $`10^6200`$ time steps and experiments were repeated several times.

Let us start with presenting examples of time series emerging from the two microdynamics studied, Fig. 2. The activity $`a=0.001`$ in all experiments. This value of activity has been found as satisfactory low to be sensitive to observe critical cluster distribution, i.e., the power-law decay of distribution tails. . At average the activity of the whole market is constants, but existing tree structures lead to strongly synchronized actions. In effect we observe huge impacts on a price change after each decision done by a member of the class A if there are only few trees on a lattice, see Fig. 2(a) and Fig. 2(b). In case of the model A isolated high changes of a price are presents. In case of the model AB an avalanche of price changes developes. If the number of leaders grows, see Fig. 2(c), trees of dependencies are smaller and values of returns go down. However we still observe bursts of activity. The volatility, by means the standard deviation, of the recorded data for different values of $`b`$ is as following:

$$\begin{array}{ccc}b\hfill & \hfill modelA& \hfill modelAB\\ 0.001\hfill & \hfill 21.2& \hfill 31.2\\ 0.005\hfill & \hfill 38.5& \hfill 32.6\\ 0.01\hfill & \hfill 39.5& \hfill 30.7\\ 0.1\hfill & \hfill 21.8& \hfill 20.4\\ 0.5\hfill & \hfill 19.5& \hfill 19.4\\ 1.0\hfill & \hfill 19.3& \hfill 19.3\end{array}$$
(1)
So that one can say that for $`b>0.1`$ the volatility in both systems is close to the volatility of the pure CB-model.

It occurs that the two-point autocorrelation function decays sharply in 1 step independently of $`b`$ as far as the model A is considered. However in case of the model AB the presence of delay trees transmits into weak correlation. We observe nonzero correlation but $`<0.2`$ for few time steps if $`b<0.4`$, see Fig. 3.

On Fig. 4 we show distributions of normalized returns. The normalization means dividing each return by the standard deviation. The synchronized decisions of clusters result in huge returns and what follows in fat wings of the return distribution. However the distribution arisen from the model A points at the absence of small returns also. Close to the zero-return peak, symmetrically on its both sides, we observe two minimums. But with $`b0.1`$ both models lead to distributions of returns hardly distinguishable from the CB model distribution. It seems that effect on a price caused by actions of more than 40 independent groups of investors is similar to that ones which are undertaken by the whole 400 independent investors possible to exists on $`200\times 200`$ lattice.

## 4  Polish stock market and imitation

We have examined 73 series of prices of stocks from Polish stock market to find signs of the model A: dips around the zero-return peak, or the model AB: fat tails. We analysed the date which were collected during last six or seven years. The Warsaw Stock Exchange regulations on the minimal price change imply that the size of a histogram bin has to be carefully selected. In the following figures we plot distributions of logarithmic returns, namely, return(t)= log price(t) -log price (t-1), which are normalized by their standard deviation. The minimal price change rule leads to the separation from the return zero less than 0.064. Therefore we decided that the value 0.20 for the size for a histogram bin is proper. Accordingly to the distribution shape the stocks have been divided into two groups: regular - with the distribution similar to the distribution of WIG- the Warsaw Stock Exchange market index, and others, see Fig. 5. Let us point at the minimum on the left side of the zero-return peak. To identify stocks with irregular distribution we searched for the zero-returns properties and fat tails. Series classified as possible examples which exhibit properties of the model A and the model AB are collected in Fig. 6. Together we show time series of these stocks. These stocks represent weak firms: low capital, low liquidity, often at the edge of the bankruptcy. The DFA method  for quantification of long-termed correlation in time-series provides the crossover to the antipersistency behavior for long time distances. 40 stocks of all 73 considered by us are found with the corresponding exponent $`\alpha <0.48`$. Among them 19 shows $`\alpha <0.40`$.

## 5 Closing

Presented results are only preliminary and further investigations should be undertaken. The proposed microdynamics leads to the non-realistic global effect in a sense that each price change is accepted. In real markets there are rules preventing from unlimited large or infinitely small changes. Nevertheless, the observed on the Polish market time series justify the modifications introduced by us to the original CB-model. Specially when emerging markets are considered then the assumption about imitating others in place of undertaking own reasoning seems to be legalized. However one must be careful in a quantitative conclusion suggested when comparing the model A distribution to a real series distribution that a stock is driven by not more than a few trees of investors.

There are some natural possibilities for further model development. The existing extension of the CB model such as adjusting activity with respect to the latest price movements  offers the straightforward way. The next possibility is to inject imitating agents, agents of the class B, into different from the CB model dynamics. For example, it should be interesting to consider the Lux-Marchesi model  with speculative investors divided into imitation trees. Finally, one should consider switching imitating agents into overconfident investors, which are know to be present on well developed markets . Such investors being overconfident about own reasoning, understand better a market than other investors. Therefore each of overconfident investors can predict with certainty what some other investor will do at the next time step. Having such a knowledge he modifies his present decision. Thus, he buys now if the other investor will buy, he sells if the other investor will sell. Simulations with described microdynamics are performed and soon will be presented.

Acknowledgement
We are grateful to Professor Ditrich Stauffer for the source code of Cont-Bouchaud model. This work is supported by Gdańsk University: BW 5400-5-0014-3. The simulations, partially, have been performed in TASK- Academic Computer Center in Gdansk.