text:
# Nonlinear magneto-optical rotation with frequency-modulated light

## Abstract

A magnetometric technique is demonstrated that may be suitable for precision measurements of fields ranging from the sub-microgauss level to above the Earth field. It is based on resonant nonlinear magneto-optical rotation caused by atoms contained in a vapor cell with anti-relaxation wall coating. Linearly polarized, frequency-modulated laser light is used for optical pumping and probing. If the time-dependent optical rotation is measured at the first harmonic of the modulation frequency, ultra-narrow ($``$ a few Hz) resonances are observed at near-zero magnetic fields, and at fields where the Larmor frequency coincides with half the light modulation frequency. Upon optimization, the sensitivity of the technique is expected to exceed $`10^{11}G/\sqrt{Hz}`$.

When light near-resonant to an atomic transition propagates through an atomic medium immersed in a magnetic field, the light polarization can be affected. For example, when a magnetic field is applied along the direction of light propagation, there is light-power-dependent rotation of the polarization plane known as nonlinear magneto-optical (Faraday) rotation (NMOR). Recently Budker et al. (1998), we observed ultra-narrow ($`1\mathrm{Hz}`$) zero-field resonances in NMOR with rubidium atoms contained in vapor cells with high-quality anti-relaxation coating Alexandrov et al. (1996). These resonances arise due to preservation of atomic polarization over thousands of collisions with the walls of the cell. The sensitivity of an NMOR-based magnetometer for sub-microgauss fields could, in principle, exceed $`10^{11}G/\sqrt{Hz}`$ Budker et al. (2000), approaching the fundamental shot-noise limit (given the number of atoms in the cell and the polarization relaxation rate). In the present work, we describe a new technique that involves both inducing and measuring NMOR with a single frequency-modulated light beam. Additional ultra-narrow resonances are observed when the light modulation frequency coincides with twice the Larmor frequency, 2$`\mathrm{\Omega }_L`$ (and also with $`\mathrm{\Omega }_L`$, etc.), the use of which can extend the dynamic range of an NMOR-based magnetometer to above the Earth field range ($`1`$G) of interest in many applications.

The high sensitivity of NMOR magnetometry is based in part on the ability of the polarimeter to measure small optical rotation angles. The sensitivity of polarimeters without modulation is limited by low-frequency noise. Therefore, the most successful polarimetric techniques always involve some type of fast modulation. In our previous work Budker et al. (1998, 2000), we used polarization-modulation polarimetry, where in addition to the vapor cell under investigation, a Faraday rotator is inserted between the crossed polarizer and analyzer. The Faraday rotator modulates the polarization direction at a frequency of $`1`$kHz, and the rotation angle is determined by measuring the first harmonic of the signal from a photodetector in the dark output of the analyzer. Such a method allows detection of signals at high frequencies, eliminating sensitivity to excess low-frequency noise. However, the polarization modulation technique has a number of shortcomings. Importantly, it is not immune to the problem of “drifting zero” – any relative rotation of the polarizer and analyzer, or a change in the birefringence of optical elements, etc. is detected as a rotation signal indistinguishable from rotation produced by the atoms. In order to avoid such problems, a light frequency modulation technique was introduced in Refs. Barkov and Zolotorev (1978); Barkov et al. (1988) and was applied to measurements of *linear* magneto-optical rotation and parity-violating optical rotation. This technique allowed detection of optical rotation at the light modulation frequency without introducing additional optical elements (such as a Faraday modulator) between the polarizer and analyzer, considerably reducing spurious rotations. Sensitivity to most spurious rotations (including the “drift of zero”) is further reduced because such effects generally do not lead to spectral features as sharp as an atomic resonance. Optical pumping with frequency-modulated light was studied and applied to $`{}_{}{}^{4}\mathrm{He}`$ resonance magnetometers Cheron et al. (1995, 1996); Gilles et al. (2001), but in that work transmission was monitored, in contrast to the measurements of optical rotation performed in the present work. The performance of the $`{}_{}{}^{4}\mathrm{He}`$ magnetometer is limited by laser intensity noise Gilles et al. (2001), which is not expected to limit the performance of NMOR-based magnetometers.

A simplified diagram of the experimental apparatus is shown in Fig. 1. A paraffin-coated vapor cell (diameter 10 cm) containing an isotopically enriched sample of Rb atoms ($`94\%`$ of $`{}_{}{}^{87}\mathrm{Rb}`$; atomic density $`7\times 10^9\mathrm{cm}^3`$ at 20C) is placed in a multi-layer magnetic shield equipped with a system of magnetic coils within the inner shield that are used to compensate residual magnetic fields (to a level of $`0.1\mu `$G), and to apply well-controlled, arbitrarily-directed fields to the cell Budker et al. (1998, 2000). The cell is placed between a polarizer and an analyzer oriented at $`45^{}`$ with respect to each other (a balanced polarimeter). The frequency of the laser (tuned near the D1-line) is modulated with an amplitude typically of a few hundred MHz. For the measurements presented here, we have chosen modulation frequencies up to $`\mathrm{\Omega }_m=2\pi \times 3\mathrm{kHz}`$ achieved by modulating a piezo actuator EN (1) in the diode laser (New Focus, Vortex 6000). The amplitudes at the first and second harmonics of $`\mathrm{\Omega }_m`$ of the difference signal between the two photodiodes at the outputs of the analyzer are detected with a lock-in amplifier (Stanford Research Systems, SR810DSP). Upon normalization by the time-averaged sum of the photodiode outputs, this represents the signal measured in this experiment (i.e. the effective amplitude of the first or second harmonic of the time-dependent polarization rotation angle).

The signals at the first and second harmonics of $`\mathrm{\Omega }_m`$ as a function of longitudinal (along the light propagation direction) magnetic field are shown in Fig. 2. The laser frequency is modulated at $`\mathrm{\Omega }_m=2\pi \times 1\mathrm{kHz}`$ with modulation amplitude $`\mathrm{\Delta }\omega =2\pi \times 220\mathrm{MHz}`$. For the first and the second harmonic signals, the laser is tuned to where the respective signals are of maximal size: the low-frequency slope of the $`F=2F^{}=1`$ resonance ($`F,F^{}`$ are the total angular momenta of the lower and upper state) for the first harmonic (Fig. 4(a)), and to the center of the $`F=2F^{}=1`$ resonance for the second harmonic.

For both harmonics, a narrow, dispersively-shaped resonance is seen in Fig. 2 near $`B=0`$. Its origin is similar to that of the ultra-narrow resonances observed in the previous work Budker et al. (1998): the atoms are pumped by linearly-polarized light into a state aligned EN (1) along the light polarization; optically-pumped atomic alignment evolves due to the applied magnetic field; the resulting atomic polarization EN (2) is “probed” by the light, leading to a modified light polarization at the output of the cell. For the zero-field resonances in Fig. 2, the modulation frequency is much faster than the Larmor frequency and the optical pumping rate for the cell, so that frequency modulation is equivalent to spectral broadening of the pumping light, which does not significantly affect the pumping process. On the other hand, modulation is essential at the probing stage. As the laser frequency is scanned through the atomic resonance, there is a time-dependent optical rotation, so the signal contains various harmonics of $`\mathrm{\Omega }_m`$. The peaks in the magnetic field dependence of the signal occur at $`\mathrm{\Omega }_Lg_F\mu B=\pm \gamma _{\mathrm{rel}}/2`$, where $`g_F`$ is the Land$`\stackrel{´}{\mathrm{e}}`$ factor ($`g_F=1/2`$ for the $`F=2`$ state of 87Rb), $`\mu 2\pi 1.40\mathrm{MHz}/\mathrm{G}`$ is the Bohr magneton, and $`\gamma _{\mathrm{rel}}`$ ($`2\pi 1\mathrm{Hz}`$ for the data in Fig. 2) is the relaxation rate of atomic polarization. The zero-field resonances can be applied to low-field NMOR-magnetometry Budker et al. (2000), with all of the advantages of the frequency-modulation technique as in the case of linear optical rotation Barkov and Zolotorev (1978); Barkov et al. (1988).

Additional resonances are seen at magnetic field magnitudes of $`714.4\mu `$G ($`\mathrm{\Omega }_m=2\mathrm{\Omega }_L`$) and $`1428.9\mu `$G ($`\mathrm{\Omega }_m=\mathrm{\Omega }_L`$) for the first and second harmonics, respectively (Fig. 2). For these resonances, there are both dispersively-shaped in-phase signals (Fig. 2(a,c)) and $`\pi /2`$ out of phase (quadrature) components peaked at the centers of these resonances (Fig. 2(b,d)). The origin of the resonances can be understood as follows (Fig. 3). As the laser frequency is modulated, the optical pumping rate changes depending on the instantaneous detuning of the laser from the atomic transition. This occurs with periodicity dictated by the modulation frequency $`\mathrm{\Omega }_m`$. On the other hand, the state of atomic polarization also changes periodically due to the presence of the magnetic field. In the case of alignment transverse to the magnetic field, this change occurs at frequency $`2\mathrm{\Omega }_L`$ (alignment precesses at $`\mathrm{\Omega }_L`$, but it returns to the same state after a half of the corresponding period). If the optical pumping rate is synchronized with Larmor precession, a resonance occurs, and the atomic medium is optically pumped into an aligned state whose axis rotates at $`\mathrm{\Omega }_L`$. The optical properties of the medium are thus being modulated at $`2\mathrm{\Omega }_L`$, causing rotation of the light polarization at this frequency. This leads to dispersively-shaped resonances in optical rotation that are detected as the signal at the appropriate harmonic of $`\mathrm{\Omega }_m`$. The quadrature components (Fig. 2(b,d)) arise because exactly on resonance, the aligned atoms produce maximum optical rotation when the alignment axis is at an angle of $`\pi /4`$ to the direction of the light polarization. (There is no quadrature component at zero magnetic field because atomic alignment does not rotate at all in this case.) Since the widths of the additional resonances at $`\mathrm{\Omega }_m=2\mathrm{\Omega }_L`$ and $`\mathrm{\Omega }_m=\mathrm{\Omega }_L`$ are determined by $`\gamma _{\mathrm{rel}}`$ as for the zero-field resonance, they are similarly narrow. In addition to these resonances, there are also smaller resonances in the first- and second-harmonic signals, generally occurring for commensurate $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_L`$. There are several mechanisms that lead to these resonances, which will be described in detail elsewhere.

The overall slope of the curves in Fig. 2 is due to the so-called transit effect (see e.g. Budker et al. (1998) and references therein), for which the peaks occur at $`B=\pm \gamma _{tr}/2g_F\mu 50`$mG, where $`\gamma _{tr}`$ is the rate of atoms’ transit through the laser beam.

The spectral dependences of various signals are shown in Fig. 4. Fig. 4(a) depicts the spectrum of the first harmonic signal obtained when the magnetic field was fixed at $`2.2\mu `$G (corresponding to maximum signal), and the central frequency of the laser was scanned across the hyperfine structure of the Rb D1-line (the transmission spectrum is shown for comparison in Fig. 4(d)). Comparing this spectrum with the NMOR spectrum obtained without frequency modulation (as in earlier work Budker et al. (1998, 2000)) shown in Fig. 4(c), it is seen that the spectral profile of the signal with frequency modulation is similar to $`d\phi (\nu )/d\nu `$, where $`\phi (\nu )`$ is the NMOR angle and $`\nu `$ is the laser frequency. Fig. 4(b) shows the spectrum obtained with the magnetic field of $`B=(714.4+2.2)\mu `$G, corresponding to an $`\mathrm{\Omega }_m=2\mathrm{\Omega }_L`$ resonance (Fig. 2). The spectral dependence of the quadrature component of the first harmonic for $`B=714.4\mu `$G (not shown) is similar to that of the in-phase component. The signals acquired at the second harmonic of $`\mathrm{\Omega }_m`$ (not shown) have spectra resembling $`d^2\phi (\nu )/d\nu ^2`$.

We have investigated the dependence of the signals on the modulation amplitude $`\mathrm{\Delta }\omega `$ and on the light power, and found similar behavior for the zero-field and the $`\mathrm{\Omega }_m=2\mathrm{\Omega }_L`$ resonances. The maximum first harmonic amplitude is achieved for modulation amplitude $`\mathrm{\Delta }\omega `$ comparable to the Doppler width of the transition $`300`$MHz. The signals are proportional to light power at low powers, peak at around $`30\mu `$W, and decrease with further power increase due to light broadening. Subsequent work will include systematic optimization of all parameters for highest magnetometric sensitivity, as it was done earlier Budker et al. (2000) for low-field NMOR with conventional polarization-modulation polarimetry. Additional improvement may be expected if the pump and probe laser beams are separated, and the modulation method is optimized (e.g., a combination of frequency and amplitude modulation). With the present unoptimized parameters and set-up, we estimate that the shot-noise-limited sensitivity of the device is $`10^{11}G/\sqrt{Hz}`$ (comparable or superior, e.g., to the most sensitive superconducting quantum interference (SQUID) sensors). Also note that shot-noise-limited operation in the Earth-field range will be facilitated by the reduced laser noise at high frequencies, and by the fact that this method uses resonant optical rotation which allows suppression of the common mode noise of the two photodetectors and spurious rotations. Due to the absence of inertia in magnetic spin-precession, a magnetometer based on this technique (which employs transverse polarization) could, in principle, also have a high band width ultimately limited by electronics.

In conclusion, we have demonstrated a novel magnetometric technique based on nonlinear optical rotation with frequency-modulated light. This technique is useful for near zero-field magnetometry and it also opens a possibility for ultra-sensitive magnetometers with broad dynamic range that includes the Earth field range, of importance to many applications, such as geophysics, magnetic prospecting and navigation. Possible new applications also include nuclear magnetic resonance spectroscopy and imaging using weakly-polarized samples that would not require any additional magnet, except the Earth itself Trabesinger et al. (2001). Narrow magneto-optical resonances may also find applications in atomic tests of fundamental symmetries Yashchuk et al. (1999). For example, it may be advantageous to carry out a search for a parity- and time-reversal invariance violating electric dipole moment of Cs in the presence of a Rb co-magnetometer in a bias magnetic field. Due to the difference between the Land$`\stackrel{´}{\mathrm{e}}`$ factors for Rb and Cs, this avoids “locking” of their polarizations to each other. In addition, presence of a bias field eliminates the effects of stray transverse fields.

The authors are grateful to E. B. Alexandrov for stimulating discussions. This work is supported by the Office of Naval Research (grant N00014-97-1-0214). D. B. also wishes to acknowledge support by an NSF CAREER grant.


format: text

text:
# ADP-01-55/T487 JLAB-THY-01-40 Baryon resonances from a novel fat-link fermion action

## 1 INTRODUCTION

Understanding the dynamics responsible for baryon excitations provides valuable insight into the forces which confine quarks inside baryons and into the nature of QCD in the nonperturbative regime. One of the long-standing puzzles in spectroscopy has been the low mass of the first positive parity excitation of the nucleon (the $`J^P=\frac{1}{2}^+`$ $`N^{}(1440)`$ Roper resonance) compared with the lowest lying odd parity excitation. Another challenge for spectroscopy is presented by the $`\mathrm{\Lambda }^{1/2}(1405)`$, whose anomalously small mass has been interpreted as indicating strong coupled channel effects involving $`\mathrm{\Sigma }\pi `$, $`\overline{K}N`$, $`\mathrm{}`$ , and a weak overlap with a three valence constituent quark state.

In this paper we present the first results of excited octet baryon mass simulations using an $`𝒪(a^2)`$ improved gluon action and an improved Fat Link Irrelevant Clover (FLIC)  quark action in which only the irrelevant operators are constructed using fat links. Configurations are generated on the new Orion computer cluster dedicated to lattice gauge theories at the CSSM at Adelaide University. After reviewing in Section 2 the main elements of lattice calculations of excited hadron masses, we describe in Section 3 various features of interpolating fields used in this analysis. In Section 4 we present results for $`J^P=\frac{1}{2}^\pm `$ nucleons and hyperons. Finally, in Section 5 we make concluding remarks and discuss some future extensions of this work.

## 2 BARYONS ON THE LATTICE

The history of excited baryons on the lattice is quite brief, although recently there has been growing interest in finding new techniques to isolate excited baryons, motivated partly by the experimental $`N^{}`$ program at Jefferson Lab. Previous work on excited baryons on the lattice can be found in Refs. .

Following standard notation, we define a two-point correlation function for a baryon $`B`$ as:

$$G_B(t,\stackrel{}{p})\underset{\stackrel{}{x}}{}e^{i\stackrel{}{p}\stackrel{}{x}}0\left|T\chi _B(x)\overline{\chi }_B(0)\right|0,$$
(1)
where $`\chi _B`$ is a (positive parity) baryon interpolating field, and we have suppressed Dirac indices. The choice of interpolating field $`\chi _B`$ is discussed in Section 3 below. For large Euclidean time, the correlation function can be written as a sum of the lowest energy positive and negative parity contributions:

$`G_B(t,\stackrel{}{p})`$ $``$ $`\lambda _{B^+}^2{\displaystyle \frac{\left(\gamma p+M_{B^+}\right)}{2E_{B^+}}}e^{E_{B^+}t}`$ (2)
$`+`$ $`\lambda _B^{}^2{\displaystyle \frac{\left(\gamma pM_B^{}\right)}{2E_B^{}}}e^{E_B^{}t},`$
where a fixed boundary condition in the time direction is used to remove backward propagating states, and where the overlap of the field $`\chi _B`$ with positive or negative parity states $`|B^\pm `$ is parameterized by a coupling strength $`\lambda _{B^\pm }`$, with $`E_{B^\pm }=\sqrt{M_{B^\pm }^2+\stackrel{}{p}^2}`$ the energy. The energies of the positive and negative parity states are obtained by taking the trace of $`G_B`$ with the operator $`\mathrm{\Gamma }_\pm `$, where

$`\mathrm{\Gamma }_\pm `$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1\pm {\displaystyle \frac{M_{B^\pm }}{E_{B^\pm }}}\gamma _4\right).`$ (3)
For $`\stackrel{}{p}=0`$, $`E_{B^\pm }=M_{B^\pm }`$ and the operator $`\mathrm{\Gamma }_\pm `$ projects out the mass, $`M_{B^\pm }`$, of the baryon $`B^\pm `$. In this case, positive parity states propagate in the 1, 1 and 2, 2 elements of the Dirac matrix of Eq. (2), while negative parity states propagate in the 3, 3 and 4, 4 elements.

## 3 INTERPOLATING FIELDS

In this analysis we consider two types of interpolating fields which have been used in the literature. The notation adopted follows that of Leinweber et al. . For the positive parity proton we use as interpolating fields:

$`\chi _1^{p+}(x)`$ $`=`$ $`ϵ_{abc}\left(u_a^T(x)C\gamma _5d_b(x)\right)u_c(x),`$ (4)
and

$`\chi _2^{p+}(x)`$ $`=`$ $`ϵ_{abc}\left(u_a^T(x)Cd_b(x)\right)\gamma _5u_c(x),`$ (5)
where the fields $`u`$, $`d`$ are evaluated at Euclidean space-time point $`x`$, $`C`$ is the charge conjugation matrix, $`a,b`$ and $`c`$ are color labels, and the superscript $`T`$ denotes the transpose. As pointed out by Leinweber , because of the Dirac structure of the “diquark” in the parentheses in Eq. (4), the field $`\chi _1^{p+}`$ involves both products of upper $`\times `$ upper $`\times `$ upper and lower $`\times `$ lower $`\times `$ upper components of spinors for positive parity baryons, so that in the nonrelativistic limit $`\chi _1^{p+}=𝒪(1)`$. Furthermore, since the “diquark” couples to a total spin 0, one expects an attractive force between the two quarks, and hence a lower energy state than for a state in which two quarks do not couple to spin 0.

The $`\chi _2^{p+}`$ interpolating field, on the other hand, is known to have little overlap with the ground state . Inspection of the structure of the Dirac matrices in Eq. (5) reveals that it involves products of upper $`\times `$ lower $`\times `$ lower components only for positive parity baryons, so that $`\chi _2^{p+}=𝒪(p^2/E^2)`$ vanishes in the nonrelativistic limit. As a result of the mixing, the “diquark” term contains a factor $`\stackrel{}{\sigma }\stackrel{}{p}`$, meaning that the quarks no longer couple to spin 0, but are in a relative $`L=1`$ state. One expects therefore that two-point correlation functions constructed from the interpolating field $`\chi _2^{p+}`$ are dominated by larger mass states than those arising from $`\chi _1^{p+}`$.

Interpolating fields for a negative parity proton can be constructed by multiplying the positive parity fields by $`\gamma _5`$, $`\chi ^B\gamma _5\chi ^{B+}`$, which reverses the role of the terms in Eq. (2). While the masses of negative parity baryons are obtained directly from the (positive parity) interpolating fields in (4) and (5) by using the parity projectors $`\mathrm{\Gamma }_\pm `$, it is instructive nevertheless to examine the general properties of the negative parity interpolating fields. In contrast to the positive parity case, both the interpolating fields $`\chi _1^p`$ and $`\chi _2^p`$ mix upper and lower components, and consequently both $`\chi _1^p`$ and $`\chi _2^p`$ are $`𝒪(p/E)`$. Physically, two nearby $`J^P=\frac{1}{2}^{}`$ states are observed in the nucleon spectrum. In simple quark models, the splitting of these two orthogonal states is largely attributed to the extent to which scalar diquark configurations compose the wave function . It is reasonable to expect $`\chi _1^p`$ to have better overlap with scalar diquark dominated states, and thus provide a lower effective mass in the large Euclidean time regime explored in lattice simulations. If the effective mass associated with the $`\chi _2^p`$ correlator is larger, then this would be evidence of significant overlap of $`\chi _2^p`$ with the higher lying $`N^{\frac{1}{2}}`$ states. In this event, further analysis directed at resolving these two states is warranted.

## 4 RESULTS

In this paper we report results of calculations of octet excited baryon masses performed on a $`16^3\times 32`$ lattice at $`\beta =4.60`$ with a lattice spacing of $`a=0.125(2)`$ fm. The analysis is based on a preliminary sample of 50 configurations generated on the new Orion computer cluster at the CSSM, Adelaide. For the gauge fields, a mean-field improved plaquette plus rectangle action is used, while for the quark fields, the FLIC  action is implemented. In the present analysis we have imposed fixed boundary conditions in the time direction ($`U_t(\stackrel{}{x},nt)=0\stackrel{}{x}`$), and periodic boundary conditions in spatial directions. Although the simulations were performed with both $`n=4`$ and 12 fattening sweeps, the improved gauge fields were found to be smooth after only 4 sweeps. Since the results with $`n=4`$ sweeps exhibit slightly better scaling than those with $`n=12`$ , we shall focus on the results with 4 smearing sweeps. The 12 sweep results lead to the same conclusions as presented in the following. Further details of the simulations are given in Ref. .

In Fig. 1 we show the $`N`$ and $`N^{}(\frac{1}{2}^{})`$ masses as a function of the squared pseudoscalar meson mass, $`m_\pi ^2`$. The results of the new simulations are indicated by the filled symbols (filled circles are FLIC; filled diamonds are Wilson). For comparison, we also show results from earlier simulations with domain wall fermions (DWF)  (open triangles), and a nonperturbatively (NP) improved clover action at $`\beta =6.2`$ . The scatter of the different NP improved results is due to different source smearing and volume effects: the open squares are obtained by using fuzzed sources and local sinks, the open circles use Jacobi smearing at both the source and sink, while the open diamonds, which extend to smaller quark masses, are obtained from a larger lattice ($`32^3\times 64`$) using Jacobi smearing. The empirical masses of the nucleon and the lowest $`\frac{1}{2}^{}`$ excitation are indicated by the asterisks along the ordinate. There is excellent agreement between the different improved actions for the nucleon mass, in particular between the FLIC, NP improved clover  and DWF  results. On the other hand, the Wilson results lie systematically low in comparison to these due to large $`𝒪(a)`$ errors in this action . A similar pattern is repeated for the $`N^{}(\frac{1}{2}^{})`$ masses. Namely, the FLIC, NP improved clover and DWF masses are in agreement with each other, while the Wilson results again lie systematically low. A mass splitting of around 400 MeV is clearly visible between the $`N`$ and $`N^{}`$ for all actions, including the Wilson action, contrary to previous claims .

Figure 2 shows the ratio of the masses of the $`N^{}(\frac{1}{2}^{})`$ and the nucleon. Once again, there is good agreement between the FLIC and DWF actions. However, the results for the Wilson action lie above the others, as do those for the anisotropic D234 action . The D234 action has been mean-field improved, and uses an anisotropic lattice which is relatively coarse in the spatial direction ($`a0.24`$ fm). This is an indication of the need for nonperturbative or fat link improvement.

The mass splitting between the two lightest nearby $`N^{}(\frac{1}{2}^{})`$ states ($`N^{}(1535)`$ and $`N^{}(1650)`$) can be studied by considering the $`\chi _1`$ and $`\chi _2`$ interpolating fields in Eqs.(4) and (5). Recall that the “diquarks” in $`\chi _1`$ and $`\chi _2`$ couple differently to spin, so that even though the correlation functions built up from the $`\chi _1`$ and $`\chi _2`$ fields will be made up of a mixture of many excited states, they will have dominant overlap with different states yielding different masses . The results, shown in Fig. 3 for the FLIC action, indicate that indeed the $`N^{}(\frac{1}{2}^{})`$ corresponding primarily to the $`\chi _2`$ field (labeled “$`N_2^{}`$”) lies systematically above the $`N^{}(\frac{1}{2}^{})`$ associated primarily with the $`\chi _1`$ field (“$`N_1^{}`$”).

As has long been known, the positive parity $`\chi _2`$ interpolating field (“$`N_2`$”, which is also sometimes denoted by “$`N^{}(\frac{1}{2}^+)`$”) does not have good overlap with the nucleon ground state , and corresponds to a state which lies around 400 MeV above the negative parity excitation of $`\chi _1`$. There is little evidence that this state is the $`N^{}(1440)`$ Roper resonance (first $`\frac{1}{2}^+`$ excitation of the nucleon). While it is possible that the Roper resonance may have a strong nonlinear dependence on the quark mass at $`m_\pi ^2<0.2`$ GeV2, arising from pion loop corrections, it is unlikely that this behavior would be so dramatically different from that of the $`N^{}(1535)`$ so as to reverse the level ordering obtained from the lattice. A more likely explanation is that the $`\chi _2`$ interpolating field does not have good overlap with either the nucleon or the $`N^{}(1440)`$, but rather a (combination of) excited $`\frac{1}{2}^+`$ state(s).

Recall that in a constituent quark model in a harmonic oscillator basis the mass of the Roper is higher than the mass of the lowest $`P`$-wave excitation. The lattice data thus appear to be consistent with the naive quark model expectation at large values of $`m_q`$. Better overlap with the Roper resonance is likely to require more exotic interpolating fields.

## 5 CONCLUSION

We have presented the first results for the excited baryon spectrum from lattice QCD using an $`𝒪(a^2)`$ improved gauge action and an improved Fat Link Irrelevant Clover (FLIC) quark action in which only the links of the irrelevant dimension five operators are smeared. The simulations have been performed on a $`16^3\times 32`$ lattice at $`\beta =4.60`$, providing a lattice spacing of $`a=0.125(2)`$ fm. The analysis is based on a set of 50 configurations generated on the new Orion computer cluster at the CSSM, Adelaide.

Good agreement is obtained between the FLIC and other improved actions, such as the nonperturbatively improved clover  and domain wall fermion  actions, for the nucleon and its chiral partner, with a mass splitting of $`400`$ MeV. Our results for the $`N^{}(\frac{1}{2}^{})`$ improve on those using the D234  and Wilson actions. Despite strong chiral symmetry breaking, the results with the Wilson action are still quite reasonable rendering earlier conjectures invalid. Using the two standard nucleon interpolating fields, we also confirm earlier observations  of a mass splitting between the two nearby $`\frac{1}{2}^{}`$ states. We find no evidence of overlap with the $`\frac{1}{2}^+`$ Roper resonance.

We have not attempted to extrapolate the lattice results to the physical region of light quarks, since the nonanalytic behavior of $`N^{}`$’s near the chiral limit is not as well understood yet as that of the nucleon . It is vital that future lattice $`N^{}`$ simulations push closer towards the chiral limit. On a promising note, our simulations with the 4 sweep FLIC action are able to reach relatively low quark masses ($`m_q60`$–70 MeV) already. We have also not addressed the question of to what extent quenching may affect any of our results. We naturally expect the effects of quark loops to be relatively unimportant at the currently large quark masses, although quenching may well produce some artifacts as one nears the chiral limit.

For future work, we intend to use variational techniques to better resolve individual excited states, for instance, those corresponding to the $`N_1^{}`$ and $`N_2^{}`$ fields (using a $`2\times 2`$ correlator matrix). In order to further explore the origin of the Roper resonances, more exotic interpolating fields involving higher Fock states, or other nonlocal operators should be investigated. Finally, the present $`N^{}`$ mass analysis will be extended in future to include $`NN^{}`$ transition form factors through the calculation of three-point correlation functions.

We thank D.G. Richards for providing the data points from Ref. . This work was supported by the Australian Research Council, and the U.S. Department of Energy contract DE-AC05-84ER40150, under which the Southeastern Universities Research Association (SURA) operates the Thomas Jefferson National Accelerator Facility (Jefferson Lab).


text:
# Two-dimensional gauge theories of the symmetric group 𝑆_𝑛 and branched n-coverings of Riemann surfaces in the large-n limit.

## 1 Introduction

The problem of counting in how many ways a Riemann surface can be covered $`n`$-times without allowing folds but allowing branch points is equivalent to computing the partition function of a lattice gauge theory defined on a discretization of the surface, with the symmetric group $`S_n`$ as the gauge group. The choice of the action for the gauge theory determines the branch-point structure of the coverings one wants to count.

Therefore 2d $`S_n`$ gauge theory is of interest not only beacuse it is a simple but non-trivial non-abelian gauge theory, but also because it has a natural interpretation in terms of a 2d string theory, that is a theory of maps from a 2d world sheet to a 2d target space. Moreover the strings described by $`S_n`$ gauge theory, that can have branch points but not folds, are precisely the ones that enter the string-theoretic formulation of $`SU(N)`$ 2d gauge theory, discovered in . Finally, branched coverings of Riemann surfaces appear naturally when one quantizes 2d Yang Mills theory in the unitary gauge. This gauge choice was shown in Refs. to actually define a different theory, with a natural intepretation in terms of 2d strings.

In this work we will consider 2d $`S_n`$ gauge theory on a genus-zero Riemann surface, and in particular the theory that corresponds to coverings with quadratic branch points only, in the large-$`n`$ limit, and study its phase diagram. It turns out that one can obtain a non-trivial phase diagrams by letting the number of plaquettes $`p`$ of the discretized target Riemann surface scale with the order $`n`$ of the coverings as

$$p=An\mathrm{log}n$$
(1)
For example if one considers the theory in which exactly one quadratic branch point is placed on each plaquette, one finds that at $`A=\frac{1}{2}`$ there is a phase transition separating a phase in which the covering has more than one connected component from one in which the covering surface is connected. By allowing each plaquette to have either no branch point or one quadratic branch point with a given probability, the phase diagram becomes more complex, with lines of first and second order phase transitions.

These results can be obtained by writing the partition function of the theory, as it is customary in 2d lattice gauge theories, as a sum over irreducible representations of $`S_n`$, and then looking for the representations that give the highest contribution in the large-$`n`$ limit. This procedure is outlined in Sec. 2.

Some insight into the nature of the various phase transitions can be gained from an exact correspondence, that we describe in Sec. 3, between 2d gauge theories on genus-zero Riemann surfaces and random walks on the gauge group. The correspondence is completely general with respect to the choice of the gauge group and the action. If one applies it to our gauge theory, one can map it into a particular random walk on $`S_n`$ that has been studied in the mathematical literature as an example of cutoff phenomenon. In Sec. 4 we will show that one of the phase transition lines we find corresponds to the cutoff phenomenon in this random walk.

Finally, another precise connection can be made with the theory of random graphs, and precisely with what are called phase transitions in random graphs. In Sec. 5 we will describe this connection, and use it to argue that the phases of our model are distinguished by the connectivity properties of the world-sheet.

Sec. 6 will be devoted to some final remarks. In this contribution we just present our results: for more details and derivations, we refer the reader to Ref..

## 2 The model and the phase diagram on the sphere

The equivalence between 2d $`S_n`$ gauge theories on Riemann surfaces and the problem of counting branched coverings of such surfaces is described in detail in Ref.. Such equivalence can be established for all possible branch point structure of the coverings. Here we will consider quadratic branch points only: On each plaquette of the target Riemann surface there is, with probability $`x`$, a quadratic branch point, or, with probability $`1x`$, no branch point.

The partition function on a disk depends on the branch point type on the boundary, that is, in the language of the gauge theory, on the holonomy $`QS_n`$ at the boundary of the disk, and can be written as

$`Z_{n,\mathrm{disk},p}(Q)`$
$`={\displaystyle \frac{1}{n!}}{\displaystyle \underset{r}{}}d_r\mathrm{ch}_r(Q)\left[(1x)+x{\displaystyle \frac{\mathrm{ch}_r(\mathrm{𝟐})}{d_r}}\right]^p`$ (2)
where $`p`$ is the number of plaquettes, the sum runs over all the irreducible representations $`r`$ of $`S_n`$, $`d_r`$ is the dimension of the representation $`r`$ and $`\mathrm{ch}_r(\mathrm{𝟐})`$ is the character of a transposition in $`r`$. The problem of computing the partition function Eq.(2) in the large-$`n`$ limit consists in finding which representations give the largest contribution to the sum when $`n\mathrm{}`$.

The irreducible representations of $`S_n`$ are in one-to-one correspondence with the Young tableaux with $`n`$ boxes. The results of the analysis can be expressed in terms of the shape of the Young tableaux that dominate the sum in the large $`n`$-limit. We will restrict ourselves to the case of the sphere, where $`Q`$ is the identical permutation.

If one lets $`p`$ scale with $`n`$ as in Eq.(1) one finds three phases in the $`(x,A)`$ plane, shown in Fig. 1:

* In phase I, the Young tableaux that dominate the sum have both rows and columns scaling as $`n^{1/2}`$.
* In phase II the sum is dominated by Young tableaux made of a single row of length $`2Ax`$ attached to a part of area $`(12Ax)n`$ whose rows and columns scale like $`n^{1/2}`$.
* In phase III the sum is dominated by a Young tableau made of a single row (or column) of length $`n`$.

The free energy in the large-$`n`$ limit can be written exactly, so that one can determine the order of the various transitions. It turns out that the (I,II) and (I,III) transitions are first order, while (II,III) is second order, characterized by a discontinuity of the second derivative of the free energy with respect to $`A`$.

What this type of analysis does not easily give is a characterization of the three phases in terms of some order parameter of clear physical meaning. This problem will be solved in part by mapping the theory first into a random walk on $`S_n`$, then into a classical problem of random graph theory. This will be done in the next two sections.

## 3 2d gauge theories are random walks on the gauge group

A two-dimensional gauge theory on a disk is equivalent to a random walk on the gauge group manifold, the area of the disk being identified with the number of steps and the gauge theory action with the transition probability at each step.

This result is completely general with respect to the choice of the gauge group and the action, as is proven in Ref.. Here we show how this result emerges in the $`S_n`$ gauge theory defined by Eq.(2) for $`x=1`$ (the extension to general $`x`$ is trivial), that is where all plaquette variables are forced to be equal to a transposition.

As shown by Eq.(2) the partition function depends only on the holonomy $`Q`$ and the area of the disk. It follows that we can freely choose any cell decomposition of the disk made of $`p`$ plaquettes, e.g. the one shown in Fig. 2.

To compute the partition function is to count the ways in which we can place permutations $`P`$ on all the links in such a way that

* the ordered product of the links around each plaquette is a transposition
* the ordered product of the links around the boundary of the disk is a permutation in the same conjugacy class as $`Q`$

Now we can use the gauge invariance of the theory to fix all the radial links to contain the identical permutation. At this point the links on the boundary are forced to contain transpositions: therefore the partition function with holonomy $`Q`$ is the number of ways in which one can write the permutation $`Q`$ as an ordered product of $`p`$ transpositions. This in turn can be seen as a random walk on $`S_n`$, in which, at each step, the permutation is multiplied by a transposition chosen at random: the gauge theory partition function for area $`p`$ and holonomy $`Q`$ is the probability that after $`p`$ steps the walker is in $`Q`$.

The cutoff phenomenon in random walks was discovered in Ref. by studying this random walk on $`S_n`$ (more precisely the one corresponding to our model with $`x=1/n`$). It was shown in Ref. is that if the number of steps scales as $`An\mathrm{log}n`$, then in the large-$`n`$ limit for $`A>1/2`$ the probability of finding the walker in any given element $`QS_n`$ is just $`1/n!`$ for all $`Q`$: complete randomization has been achieved and all memory of the initial position of the walker has been erased.

In terms of the corresponding gauge theory, this statement translates into the following: for $`A>1/2`$, the partition function does not depend on the holonomy $`Q`$, and at fixed holonomy $`Q`$, does not depend on $`A`$. This is true in particular for $`Q=1`$, corresponding to the partition function on the sphere. Therefore the transition between phases II and III, found in Sec. 2, has a natural interpretation as a cutoff phenomenon in the corresponding random walk (see Ref. for a discussion about generalizing the result of Ref. to arbitrary values of $`x`$).

## 4 A connection with the theory of random graphs

The connection with the theory of random graphs is a consequence of the possibility, established in the previous section, of mapping one of the phase transitions of our theory of coverings into a cutoff phenomenon of a random walk. In fact it was shown in Ref.  that some phase transitions in random graphs can be interpreted as cutoff phenomena in certain random walks that naturally generate ensembles of random graphs.

Consider the random walk corresponding to our model, taken now with free boundary conditions, and take for simplicity $`x=1`$: at each step the permutation is multiplied by a random transposition. From the point of view of coverings, a step in which the transposition $`(ij)`$ is used corresponds to adding a simple branch point that connects the two sheets $`i`$ and $`j`$ of the covering surface. One can think of the process as the construction of a graph of $`n`$ sites, in which at each step one link, chosen at random, is added. After $`p=An\mathrm{log}n`$ steps the expected number of links is equal to the number of steps (since the number of available links is $`O(n^2)`$ the fact that the same link can be added more than once can be neglected in the large $`n`$ limit; see Ref. ).

It is a classic result in the theory of random graphs  that if the number of links $`p`$ is smaller than $`1/2n\mathrm{log}n`$ then the graph is almost certainly disconnected while for $`p>1/2n\mathrm{log}n`$ the graph is almost certainly connected, where “almost certainly” means that the probability is one in the limit $`n\mathrm{}`$. Therefore we conclude that the transition at $`A=1/2`$, for the model with free boundary conditions, is the point where the covering surface becomes connected.

Notice that the free boundary conditions are crucial for this argument to work: in the case, say, of the sphere, the corresponding random walk is forced to go back to the initial position in $`p`$ steps, so that links in the graph are not added independently and the result of Ref. do not apply. However there are strong arguments (see Ref.) suggesting that connectedness of the covering surface can be used to distinguish the phases in the case of the sphere as well: phase III should correspond to a connected covering surface while phases I and II should not.

## 5 Conclusions

We have shown that the two dimensional gauge theory of the symmetric group, that describes the statistics of branched coverings of a Riemann surface, has a remarkably rich phase structure in the large-$`n`$ limit.

The theory on the sphere can be studied by a variational approach that identifies which representations give the largest contribution in the large-$`n`$ limit. If the number of branch points is taken to scale as $`n\mathrm{log}n`$, three phases can be found, corresponding to different shapes of the dominant Young tableaux. The transitions between these lines can be located analytically and their order determined.

All two-dimensional gauge theories on a genus-0 surface can be mapped into random walks in the corresponding group manifold. In our case, this allows us to interpret one of the transition lines as a cutoff phenomenon in the corresponding random walk.

The theory on a disk, with free boundary conditions, can be studied with methods of the theory of random graphs: this allows one to show that there is a phase transition on a disk from a disconnected to a connected covering surface. From this one can argue, and with some limitations prove, that the connectness of the covering is what characterizes the different phases also on the sphere.

This last result is particularly interesting since it provides a description of the various phases in terms of geometric properties of the covering surfaces, that is in terms of quantities that are natural when one considers the model as a two-dimensional string theory rather than as a Yang-Mills theory. Along this line, it was recently proposed in  that the phase transition we find is a close analogue of the Hagedorn phase transition in Matrix String Theory described in Refs. .


text:
# On the scent of the knee — air shower measurements with KASCADE

## 1 INTRODUCTION

The earth’s atmosphere is continuously bombarded by highly relativistic ionized particles, first discovered and named ”cosmic rays” by V. Hess in 1912. Present–day experiments show the cosmic–ray energy spectrum extending up to more than $`10^{20}`$ eV. The flux spectrum follows a power law $`dN/dEE^\gamma `$ over many decades in energy. The most prominent feature is the knee in the spectrum around 3 PeV where the spectrum steepens from $`\gamma 2.7`$ to $`\gamma 3.1`$. The origin of cosmic rays is still under debate. Strong, relativistic shock fronts expanding from supernova explosions are favoured by popular models for the acceleration of ionized particles. Such models explain the particle acceleration up to energies of about $`Z10^{15}`$ eV, with the nuclear charge $`Z`$ of the particle. This upper limit coincides for primary protons with the mentioned steepening of the spectrum, and the origin of the knee is related to the upper limit of acceleration in several models.

Since the charged particles are deflected in the interstellar magnetic fields, the only hint of their sources are their energy spectrum and the mass composition, or more preferable, the energy spectra of individual elements. Cosmic rays at energies below 1 PeV have been directly observed by balloon–borne instruments at the top of the atmosphere or in outer space. At higher energies, the steep falling flux spectrum requires large detection areas or long observation periods, presently only possible in ground–based installations. These detector systems measure the secondary particles produced by cosmic rays in the atmosphere, the extensive air showers (EAS).

To investigate the cosmic rays from several $`10^{13}`$ eV up to $`10^{17}`$ eV the air shower experiment KASCADE (”Karlsruhe Shower Core and Array DEtector”)  has been built on–site at the Forschungszentrum Karlsruhe in Germany. The experiment detects the three main components of EAS simultaneously. A $`200\times 200`$ m2 scintillator array measures the electromagnetic and muonic components. The 320 m2 central detector system combines a large hadron calorimeter  with several muon detection systems . In addition, high energetic muons are measured by an underground muon tracking detector .

## 2 INVESTIGATION OF HADRONIC INTERACTIONS

In order to obtain the mass and energy of the shower–inducing particle, measured EAS observables are compared with predictions of simulations. These simulations describe the development of an EAS in the atmosphere and, thereafter, the signal response for all individual particles hitting the detectors. A critical part in the simulation chain is the model used to describe the high–energy hadronic interactions, since the model has to extrapolate into kinematical and energy regions not covered by present–day collider experiments. The Karlsruhe EAS simulation program CORSIKA  provides several high–energy hadronic interaction models — HDPM, DPMJET, NEXUS, QGSJET, SIBYLL, VENUS — based on different phenomenological descriptions. Below 80 GeV the codes GHEISHA and URQMD are available. One objective of the KASCADE experiment is to evaluate these models and to provide criteria for their improvement.

The hadron calorimeter is a valuable detector for testing the interaction models. The structure of the hadronic component is examined in energy and spatial coordinates. Observables used include the number of hadrons as well as their energy sum, their lateral distribution and their energy spectrum, the energy of each individual hadron relative to the most energetic hadron in an EAS, the maximum hadron energy, and the spatial distribution of the hadrons. All observables are investigated as functions of the number of electrons and muons as well as of the hadronic energy sum.

An example of such investigations is given in figure 1. Shown are frequency distributions for the number of reconstructed hadrons above 50 GeV and the hadronic energy sum for two shower size intervals for EAS with shower core inside the calorimeter. KASCADE findings are compared with predictions of the models QGSJET, DPMJET, and NEXUS for two extreme mass scenarios, the primaries being only protons or iron nuclei . Combining all observables, it turned out that the model QGSJET delivers the most reliable description . This model is used for the analyses described in the following sections.

## 3 ENERGY SPECTRA AND MASS COMPOSITION

Detailed investigations have been performed with the KASCADE data to study the influence of different observables and hadronic interaction models on the inferred cosmic–ray mass composition.

An example is shown in figure 2. In a Bayesian analysis , the mean logarithmic mass $`\mathrm{ln}A=_ir_iA_i`$, where $`r_i`$ is the relative abundance of elements with mass $`A_i`$, is calculated for different combinations of observables. Several electromagnetic, muonic, and hadronic observables have been used, including number of electrons $`N_e`$ and muons $`N_\mu ^{tr}`$, number of hadrons above 100 GeV $`N_h^{100}`$, their energy sum $`E_h`$, the maximum energy of a hadron per shower $`E_h^{max}`$, the number of high energy muons ($`E_\mu >2`$ GeV) $`N_\mu ^{}`$, and their spatial distribution $`D_6`$. The left graph summarizes results obtained including the electromagnetic shower size $`N_e`$, whereas the results in the right panel are obtained without $`N_e`$. The combinations including $`N_e`$ lead to smaller $`\mathrm{ln}A`$ values as compared to those without. The total spread between the results is $`\mathrm{\Delta }\mathrm{ln}A0.8`$. Similar systematic effects have been observed earlier for different hadronic observables . These discrepancies are most probably caused by inconsistencies in the hadronic interaction models used to interpret the EAS data.

Several analysis methods have been carried out to obtain the cosmic–ray energy spectrum. The KASCADE data show a knee in the electromagnetic, muonic, and hadronic shower size spectra . From the size spectra, the primary energy spectrum has been derived (see e.g. ), including the absolute flux spectrum obtained for the first time using hadronic observables .

A recent analysis  uses electromagnetic and muonic shower size spectra in three different zenith angle bins. With a four component assumption for the mass composition of primary cosmic rays (protons, helium, CNO-group, and iron group), an unfolding algorithm is applied, taking into account shower fluctuations and experimental effects. The hadronic interactions have been simulated using the model QGSJET above 80 GeV and the GHEISHA code below. Individual energy spectra for the four mass groups are obtained as shown in figure 3. Each spectrum shows a knee–like structure. The energy dependence of these cut–offs suggests a rigidity–dependent behaviour.

The mean logarithmic mass calculated from the individual energy spectra is shown in figure 4, indicating an increase of the average mass above the knee.

## 4 CONCLUSION

The EAS observables obtained by KASCADE are sensitive to hadronic interaction models used in simulations to interpret EAS data. At present, the combination CORSIKA/QGSJET best describes the measurements. The systematic dependence of the mean logarithmic mass on different observables and models has been investigated. The systematic error for the model QGSJET is in the order of $`\mathrm{\Delta }\mathrm{ln}A0.8`$. The reconstruction of individual energy spectra for groups of elements indicate a rigidity dependent cut–off, which explains the knee in the all particle cosmic–ray flux spectrum.

The KASCADE experiment is supported by the Ministry for Research of the German government and embedded in collaborative WTZ projects between Germany and Romania (RUM 97/014), Poland (POL 99/005), and Armenia (ARM 98/002). The Polish group acknowledges the support by KBN grant no. 5 P03B 133 20.


text:
# Galactic Halo Cusp versus Core: Tidal Effects in Mergers

## 1. Introduction

The ‘standard’ model of cosmology, CDM, which assumes hierarchical buildup of structure, is facing difficulties in explaining observed properties of galaxies, such as the number density of dwarfs (e.g., Klypin et al. 1999b), the angular-momentum crisis (e.g., Navarro & Steinmetz 2000), and the cusp/core problem. Our approach in addressing these problems within CDM is to isolate and model in simple physical terms the key relevant processes, as a guide for possible solutions. We first model the buildup of dark-matter (DM) halos in N-body simulations based on tidal effects, and then incorporate the inevitable baryonic feedback processes in an attempt to explain the apparent discrepancies. We address the angular momentum problem in Maller & Dekel (2002) and Dekel & Maller (2002), and summarize here our progress in the cusp/core problem (Dekel & Devor 2002; Dekel et al. 2002).

Cosmological N-body simulations have revealed that the density profiles of DM halos scatter about a universal shape, $`\rho (r)=\rho _\mathrm{c}(r/r_\mathrm{c})^\alpha (1+r/r_\mathrm{c})^{\alpha 3}`$, with an inner cusp of slope $`\alpha `$. Navarro, Frenk & White (1995, NFW) found this function, with $`\alpha 1`$, to be a good fit in the range $`(0.011)R_{\mathrm{vir}}`$ for different hierarchical cosmological scenarios. High-resolution simulations of a few individual halos (Moore et al. 2001; Klypin et al. 2001) found that the cusp could be as steep as $`\alpha 1.5`$, though it may flatten towards $`\alpha 1`$ at $`r<0.01R_{\mathrm{vir}}`$ (private comm. with Navarro, Frenk, Springel & White). While the formation of a cusp with $`1\alpha 1.5`$ has been established in the simulations, a basic understanding of its origin is still lacking. An even more intriguing puzzle is introduced by observations of low surface-brightness galaxies, whose centers are dominated by their DM halos, which indicate flatter inner cores with $`\alpha 0`$ (de Block et al. 2001). This seems to introduce a challenge to the CDM paradigm.

In the first two sections we develop a toy model for tidal stripping and test it against an N-body simulation. In §2 we describe the compression and rapid steepening in a core. In §3 we address stripping where $`\alpha >1`$. In §4 we analyze the convergence to an asymptotic profile. In §5 we discuss our results.

## 2. Core into Cusp due to Tidal Compression

We consider a fixed spherical halo of mean density profile $`\overline{\rho }(r)M(r)/r^3`$. Denote $`\alpha (r)d\mathrm{ln}\overline{\rho }/d\mathrm{ln}r`$ such that locally $`\overline{\rho }r^\alpha `$, with $`\alpha (r)`$ constant or increasing in the range $`0\alpha 3`$. Consider a satellite at $`r`$, spiraling into the halo under gravity and dynamical friction. The maximum tidal force by the halo on a unit mass at satellite radius $`\mathrm{}`$$`\mathrm{}`$$`\mathrm{}`$ ($`\mathrm{}r`$) is

$$𝐅_{\mathrm{tide}}=\frac{\mu (r)\mathrm{}}{r^3}\widehat{\mathrm{}\mathrm{}\mathrm{}},\mu (r)=2M(r)r\frac{dM}{dr}=(\alpha 1)M(r).$$
(1)
The maximum is obtained along the line connecting the halo centers, but while the radial component is smaller in all other directions, it is of the same sign. The familiar case is of a point mass, $`\alpha =3`$, for which the tidal force is pulling outwards from the satellite. For flatter halo slopes it becomes weaker $`(\alpha 1)`$. An important feature that is often overlooked is that when the halo density profile is flat enough, $`\alpha <1`$, the tidal force reverses direction into compression, resulting in accretion rather than stripping. Note that the critical slope of unity coincides with the cusp of the NFW profile, and that in a core, $`\alpha 0`$, the tides induce strong compression.

This effect can be demonstrated using an N-body simulation of a merger (following Mihos & Hernquist 1996). We use a total of $`10^5`$ DM particles, the mass ratio is 1:10, the satellite spirals in on a quasi-circular orbit, and the initial profile is a truncated isothermal sphere with a core, where $`\alpha `$ ranges in practice from 0.6 to 2.8. In Fig. 1 (left), the stripping of each satellite shell, marked by the onset of a steep rise, can be identified with a halo radius $`r`$, with $`\alpha (r)`$ decreasing from 3 downwards as one moves from outer shells inwards. However, as the stripping point approaches $`\alpha 1`$, inwards to the 30% satellite mass shell, the stripping stops and some shrinking can be seen instead.

The compression at $`\alpha <1`$ implies that any part of the satellite which makes it intact into the halo core would sink towards the center without further stripping. This should cause a rapid steepening of the profile to $`\alpha >1`$, as seen in Fig. 1 (right). The NFW inner slope of $`\alpha =1`$ is thus a robust lower bound, as we know from cosmological simulations; a flatter density core cannot survive as long as satellites deposit enough mass in the inner halo.

## 3. Tidal Puffing-Up and Stripping at $`\alpha >1`$

At $`\alpha >1`$, the effects are more subtle. Let $`\mathrm{}`$ mark shell radii within the unperturbed satellite, whose mean-density profile is $`\overline{\sigma }(\mathrm{})m(\mathrm{})/\mathrm{}^3`$. Assume that when it is at halo radius $`r`$, mass is lost beyond a momentary stripping radius $`\mathrm{}(r)`$ and is added to the halo at $`r`$ (on average). We wish to determine the correspondence between $`\mathrm{}`$ and $`r`$ at stripping. Define $`\psi (r,\mathrm{})\overline{\rho }(r)/\overline{\sigma }(\mathrm{})`$. Traditionally, the stripping radius is assumed to be determined by the resonance condition $`\psi (r,\mathrm{})=1`$, but this ignores the earlier effects of tides on the satellite structure. A key new feature in our analysis is that, as $`r`$ decreases, the tides stretch the satellite orbits, which can be modeled as an effective puffing up of the relevant shells before they are being torn away. We define for every shell $`\mathrm{}`$ when the satellite is at $`r`$ a momentary puffing factor by $`p(r,\mathrm{})\mathrm{}_\mathrm{p}/\mathrm{}`$, where $`\mathrm{}_\mathrm{p}(r)`$ is the momentary shell radius. One can then show that the resonance condition, for $`r`$ and $`\mathrm{}`$ at stripping, becomes

$$\psi =\alpha ^1p^3.$$
(2)
In the regime where $`\alpha >1`$, we expect puffing, $`p>1`$, so the corrected resonance condition implies $`\psi <1`$, and for large $`\alpha `$ even $`\psi 1`$. This means more efficient stripping compared to the old condition ignoring puffing. To obtain an explicit stripping condition we wish to express $`\psi `$ as a function of $`\alpha `$, so we need to estimate how $`p(\mathrm{})`$ evolves as the satellite falls into smaller $`r`$ positions. By applying an adiabatic invariant, we obtain a puffing equation for any shell $`\mathrm{}`$ when it is at $`r`$:

$$p(\alpha 1)\psi p^4=1.$$
(3)
In the two equations above we have omitted geometrical factors of order unity, to be calibrated later using simulations (see Dekel & Devor 2002). When we combine the above equations we obtain a new stripping condition:

$$\frac{\overline{\rho }(r)}{\overline{\sigma }(\mathrm{})}=\psi [\alpha (r)]=\{\begin{array}{cc}\alpha ^4\hfill & 1<\alpha <\alpha _\mathrm{c}1.4\hfill \\ 0.1/(\alpha 1)\hfill & \alpha >\alpha _\mathrm{c}\hfill \end{array}.$$
(4)
Given the slope profile $`\alpha (r)`$, it relates every satellite shell $`\mathrm{}`$ to the position $`r`$ where it should be stripped. The function $`\psi (\alpha )`$ is thus predicted to decrease monotonically as a function of $`\alpha `$, towards values of order 0.1-0.2 at $`\alpha 2`$ and below 0.1 at $`\alpha =3`$, Fig. 2 (left). This implies that the stripping process is more efficient for steeper halo profiles. We also implement in the deposit prescription additional effects near and below $`\alpha =1`$, due to the finite size of the satellite and the tidal compression in the inner halo. These effects may cause deposit of satellite material before the stripping condition is fulfilled.

The puffing before stripping can be seen in the merger simulation, e.g., for the outer and intermediate shells in Fig. 1 (left). The magnitude of the puffing is 30-50%, as expected, corresponding to a factor of 2-3 in density. A similar effect has been qualitatively noticed in simulations before (e.g. Klypin et al. 1999a, Fig. 6). In the merger simulation, we measure the deposit radius $`r(\mathrm{})`$ from the final distribution of stripped satellite mass about the halo center \[can be read from Fig. 1 Fig. 1 (left)\]. The corresponding values of $`\psi (r,\mathrm{})`$ are plotted against $`\alpha (r)`$ in Fig. 2 (left). The qualitative agreement between the simulation result and the model predictions indicates that despite the crude approximations made, our very simplified model mimics the main features of the tidal stripping and deposit process. When puffing is ignored, $`\psi =1`$, the model clearly fails.

## 4. Halo Asymptotic Profile

If the stripping is described by a condition similar to eq. 4, with $`\psi (\alpha )`$ a decreasing function, the profile evolves slowly towards an asymptotic stable power law with $`\alpha _\mathrm{a}\stackrel{>}{}1`$. We assume that the halo and satellite are drawn from a cosmological distribution; they are homologous, with their characteristic radii and densities scaling like $`\mathrm{}_\mathrm{c}/r_\mathrm{c}m^{(1+\nu )/3}`$ and $`\sigma _\mathrm{c}/\rho _\mathrm{c}m^\nu `$, where $`\nu 0.33`$ for $`\mathrm{\Lambda }`$CDM. Fig. 2 (right) helps understanding the origin of an asymptotic slope. We write $`\overline{\rho }_{\mathrm{final}}(r)=\overline{\rho }(r)+\overline{\sigma }(\mathrm{})\mathrm{}^3/r^3`$, and obtain for the change of $`\alpha `$ in a merger

$$\mathrm{\Delta }\alpha (r)\frac{d}{dr}\left(\frac{\overline{\sigma }(\mathrm{})}{\overline{\rho }(r)}\frac{\mathrm{}^3}{r^3}\right).$$
(5)
One can see that every power law is a self-similar solution, $`\mathrm{\Delta }\alpha (r)=0`$, but not necessarily a stable one. When $`\alpha `$ is increasing with $`r`$, $`\mathrm{}/r`$ is decreasing with $`r`$. Thus, when puffing is ignored, $`\psi =`$const., one has continuous steepening, $`\mathrm{\Delta }\alpha (r)>0`$. With realistic puffing, $`1/\psi `$ is increasing with $`r`$, which can produce a stable fixed point at a certain asymptotic value $`\alpha _\mathrm{a}`$, where $`\mathrm{\Delta }\alpha =0`$ and the second derivative is negative. A rigorous linear perturbation analysis determines the rate of convergence to $`\alpha _\mathrm{a}`$ and yields an equation for its value for a sequence of mergers with the same mass ratio:

$$\mathrm{\Delta }\alpha \alpha (\alpha 3)\psi ^{}(\alpha )/\psi (\alpha )+3\mathrm{ln}[(m/M)^\nu \psi (\alpha )]=0.$$

The solutions are typically in the range $`1<\alpha _\mathrm{a}1.5`$.

In order to test the linear analysis, we perform toy simulations of the profile buildup by mergers, where we implement the stripping recipe of §3 (replacing the crude stripping recipe used in earlier work, e.g., Syre & White 1998). For given halo and satellite profiles, we solve the stripping equation for $`\mathrm{}(r)`$ and add the stripped satellite mass to the halo accordingly. Among other numerical complications, we implement a smoothing scheme to ensure that $`\alpha (r)`$ remains monotonic. We then follow a sequence of cosmological mergers and study the evolution towards an asymptotic slope. Fig. 3 (left) shows the convergence of $`\alpha `$ at a fixed $`r`$ to the asymptotic value. Fig. 3 (right) shows how the profile evolves through momentary profiles which are probably more relevant for comparison with real halos at different times during their buildup process. These figures are for certain given values of the geometrical factors, the mass ratio and $`\nu `$. In Dekel et al. (2002) we address a sequence of mergers with a cosmological distribution of mass ratios, and the robustness to the cosmological model.

## 5. Discussion

Our analysis demonstrates that the way to maintain a flat core is by disrupting satellites outside the core. This may be achieved if the cores of satellite halos are puffed up due to gas processes. As an example, the following speculative scenario is based on enhanced feedback due to tidal compression. Consider a satellite made of DM and $`10\%`$ baryons passing through the halo core towards a turn-around on the other side. Assume that the baryons have already cooled and contracted into the satellite center. The satellite loses its outer DM layers in the outer halo such that when it enters the halo core it is baryon-rich. In the core, the tides compress the satellite, creating shocks and an efficient burst of star formation. Before turn-around on the other side, there is time for the resulting supernovae to blow out the satellite gas. If the remaining satellite loses half its mass in this blow-out, its density drops by a factor between 8 to infinity, depending on whether the gas expulsion is adiabatic or impulsive. Thus, the remaining satellite becomes much more susceptible to tidal stripping, which disrupts it completely before it re-enters the halo core.

In our recent work, we address the different problems within the successful cosmological framework of CDM by appealing to inevitable feedback effects. In Maller & Dekel (2002) we address the angular-momentum catastrophe, where simulations including gas produce disks significantly smaller than the galactic disks observed, and with a different internal distribution of angular momentum. We first construct a toy model for the angular-momentum buildup by mergers based on tidal stripping and dynamical friction, which helps us understand the origin of the spin problem as a result of over-cooling in satellites. We then incorporate a simple model of feedback, motivated by Dekel & Silk (1986), and find that it can remedy the discrepancies, and in particular explain the low baryon fraction and angular-momentum profiles in dwarf disk galaxies. Feedback effects may also provide the cure to the missing dwarf problem, where the predicted number of dwarf halos in CDM is much larger than the observed number of dwarf galaxies (Bullock, Kravtsov & Weinberg 2000).

Another approach (e.g., Hogan & Dalcanton 2000) is to appeal to a Warm Dark Matter (WDM) scenario, despite the fact that it requires fine-tuning of the particle mass to $`1keV`$. The main feature of WDM is the partial suppression of small halos, which should help maintaining a core, but simulations of halos in WDM seem to still show inner cusps (Bullock, Kravtsov & Colin 2001). While the explicit merger picture modeled above may be invalid in this case, the gravitational processes involved in the halo buildup still mimic a similar behavior. We note that the tidal compression in the core may amplify density perturbations and make them behave like merging satellites.

The suppression of small halos in WDM may remedy the spin catastrophe by avoiding over-cooling (Sommer-Larsen & Dolgov 2001), but it is harder to see how it would explain the angular-momentum profile in galaxies. Furthermore, while the number of dwarfs is already suppressed in WDM, the addition of minimum feedback effects is likely to cause an overkill, where the number of dwarfs is predicted to be much lower than observed (J. Bullock, private comm.).

The success of our toy model in matching several independent observations indicates that it indeed captures the relevant elements of the complex processes involved, and in particular that feedback effects may indeed provide the cure to all three problems of galaxy formation in CDM.

## References

Bullock J.S., Kravtsov A.V. & Colin P. 2001, ApJL, in press

Bullock J.S., Kravtsov A.V. & Weinberg D.H. 2000, ApJ, 539, 517

deBlok, W.J.G., McGaugh, S.S., Bosma, A. & Rubin, V.C. 2001, ApJL, 552, L23

Dekel, A. & Devor, J. 2002, in preparation

Dekel, A., Arad, I., et al. 2002, in preparation

Dekel, A. & Maller, A.H. 2002, in The Mass of Galaxies at Low and High Redshift, eds. R. Bender & A. Renzini (Springer-Verlag, ESO Astrophysics Symposia)

Dekel A. & Silk J. 1986, ApJ, 303, 39

Hogan C.J. & Dalcanton J.J. 2000, PRD, 62

Klypin, A., Gottl ber, S., Kravtsov, A.V. & Khokhlov, A.M. 1999a, ApJ, 516, 530

Klypin, A.A., Kravtsov, A.V., Bullock, J.S. & Primack, J.R. 2001, ApJ, 554, 903

Klypin, A., Kravtsov, A.V., Valenzuela, O. & Prada, F. 1999b, ApJ, 522, 82

Maller A.H. & Dekel A. 2002, MNRAS, submitted

Moore, B., Quinn, T., Governato, F., Stadel, J. & Lake, G. 1999, MNRAS, 310, 1147

Navarro J.F., Frenk C.S. & White S.D.M. 1995, MNRAS, 275, 56

Navarro J.F. & Steinmetz M. 2000, ApJ, 538, 477

Sommer-Larsen J. & Dolgov A. 2001, ApJ, 551, 608

Syre, D. & White, S.D.M. 1998, MNRAS, 293, 337


text:
# FINDING REGULATORY SITES FROM STATISTICAL ANALYSIS OF NUCLEOTIDE FREQUENCIES IN THE UPSTREAM REGION OF EUKARYOTIC GENES.

## 1 Introduction

As more and more complete genomic sequences are being decoded it is becoming of crucial importance to understand how the gene expression is regulated. A central role in our present understanding of gene expression is played by the notion of “regulatory network”. It is by now clear that a particular expression pattern in the cell is the result of an intricate network of interactions among genes and proteins which cooperate to enhance (or depress) the expression rate of the various genes. It is thus important to address the problem of gene expression at the level of the whole regulatory network and not at the level of the single gene.

In particular, most of the available information about such interactions concerns the transcriptional regulation of protein coding genes. Even if this is not the only regulatory mechanism of gene expression in eukaryotes it is certainly the most widespread one.

In these last years, thanks to the impressive progress in the DNA array technology several results on these regulatory networks have been obtained. Various transcription factors (TF’s in the following) have been identified and their binding motifs in the DNA chain (see below for a discussion) have been characterized. However it is clear that we are only at the very beginning of such a program and that much more work still has to be done in order to reach a satisfactory understanding of the regulatory network in eukaryotes (the situation is somehow better for the prokaryotes whose regulatory network is much simpler).

In this contribution we want to discuss a new method which allows to reconstruct these interactions by comparing existing biological information with the statistical properties of the sequence data. This is a line of research which has been pursued in the last few years, with remarkable results, by several groups in the world. For a (unfortunately largely incomplete) list of references see . In particular the biological input that we shall use is the fact that some genes, being involved in the same biological process, are likely to be “coregulated” i.e. they should show the same expression pattern. The simplest way for this to happen is that they are all regulated by the same set of TF’s. If this is the case we should find in the upstream111With this term we denote the portion of the DNA chain which is immediately before the starting point of the open reading frame (ORF). We shall characterize this region more precisely in sect.3 below. region of these genes the same TF binding sequences. This is a highly non trivial occurrence from a statistical point of view and could in principle be recognized by simple statistical analysis.

As a matter of fact the situation is much more complex than what appears from this idealized picture. TF’s not necessarily bind only to the upstream region. They often recognize more than one sequence (even if there is usually a “core” sequence which is highly conserved). Coregulation could be achieved by a complex interaction of several TF’s instead than following the simple pattern suggested above. Notwithstanding this, we think that it is worthwhile to explore this simplified picture of coregulation, for at least three reasons.

* Even if in this way we only find a subset of the TF’s involved in the coregulation, this would be all the same an important piece of information: It would add a new link in the regulatory network that we are studying.
* Analyses based on this picture, being very simple, can be easily performed on any gene set, from the few genes involved in the Glycolysis (the first example that we shall discuss below) up to the whole genome (this will be the case of the second example that we shall discuss). This feature is going to be more and more important as more and more DNA array experiment appear in the literature. As the quantity of available data increases, so does the need of analytical tools to analyze it.
* Such analyses could be easily improved to include some of the features outlined above, keeping into account, say, the sequence variability or the synergic interaction of different TF’s.

To this end we have developed two different (and complementary) approaches. The first one (which we shall discuss in detail in sect.3 below) follows a more traditional line of reasoning: we start from a set of genes which are known to be coregulated (this is our “biological input”) and then try to recognize the possible binding sites for the TF’s. We call this approach the “direct search” for coregulating TF’s.

The second approach (which we shall briefly sketch in sect.4 below and is discussed in full detail in ) is completely different and is particularly suitable for the study of genome-wide DNA array experiments. In this case the biological input is taken into account only at the end of the analysis. We start by organizing all the genes in sets on the basis of the overrepresented common sequences and then filter them with expression patterns of some DNA array experiment. We call this second approach the “inverse search” for coregulating TF’s.

It is clear that all the candidate gene interactions which we identify with our two methods have to be tested experimentally. However our results may help selecting among the huge number of possible candidates and could be used as a preliminary test to guide the experiments.

This contribution is organized as follows. In sect.2 we shall briefly introduce the reader to the main features of the regulatory network (this introduction will necessarily be very short, the interested reader can find a thorough discussion for instance in ). We shall then devote sect.3 and 4 to explain our “direct” and “inverse” search methods respectively. Then we shall discuss two instances which well exemplify the two strategies. First in sect.5 we shall study the coregulation of glycolysis in Drosophila melanogaster. Second, in sect.6 we shall discuss the diauxic shift in Saccharomyces cerevisiae. The last section will be devoted to some concluding remarks.

## 2 Transcription factors.

As mentioned in the introduction, a major role in the regulatory network is played by the Transcription Factors, which may have in general a twofold action on gene transcription. They can activate it by recruiting the transcription machinery to the transcription starting site by binding enhancer sequences in the upstream noncoding region, or by modifying chromatine structure, but they can also repress it by negatively interfering with the transcriptional control mechanisms.

The main point is that in both cases TFs act by binding to specific, often short DNA sequences in the upstream noncoding region. It is exactly this feature which allows TF’s to perform a specific regulatory functions. These binding sequences can be considered somehow as the fingerprints of the various TF’s. The main goal of our statistical analysis will be the identification and characterization of such binding sites.

### 2.1 Classification

Even if TF’s show a wide variability it is possible to try a (very rough) classification. Let us see it in some more detail, since it will help understanding the examples which we shall discuss in the following sections. There are four main classes of binding sites in eukaryotes.

* Promoters

These are localized in the region immediately upstream of the coding region (often within 200 bp from the transcription starting point). They can be of two types:

+ short sequences like the well known CCAAT-box, TATA-box, GC-box which are not tissue specific and are recognized by ubiquitous TFs
    + tissue specific sequences which are only recognized by tissue specific TFs
* Response Elements

These appear only in those genes whose expression is controlled by an external factor (like hormones or growing factors). These are usually within 1kb from the transcription starting point. Binding of a response element with the appropriate factor may induce a relevant enhancement in the expression of the corresponding gene
* Enhancers

these are regulatory elements which, differently from the promoters, can act in both orientations and (to a large extent) at any distance from the transcription starting point (there are examples of enhancers located even 50-60kb upstream). They enhance the expression of the corresponding gene.
* Silencers

Same as the enhancers, but their effect is to repress the expression of the gene.

### 2.2 Combinatorial regulation.

The main feature of TF’s activity is its “combinatorial” nature. This means that:

* a single gene is usually regulated by many independent TF’s which bind to sites which may be very far from each other in the upstream region.
* it often happens that several TF’s must be simultaneously present in order to perform their regulatory function. This phenomenon is usually referred to as the “Recruitment model for gene activation” (for a review see ) and represents the common pattern of action of the TF’s. It is so important that it has been recently adopted as guiding principle for various computer based approaches to detect regulatory sites (see for instance ).
* the regulatory activity of a particular TF is enhanced if it can bind to several (instead of only one) binding sites in the upstream region. This “overrepresentation” of a given binding sequence is also used in some algorithms which aim to identify TF’s. It will also play a major role in our approach.

## 3 The “direct” search method.

In this case the starting point is the selection of a set of genes which are known to be involved in the same biological process (see example of sect. 5).

Let us start by fixing a few notations:

* Let us denote with $`M`$ the number of genes in the coregulated set and with $`g_i`$, $`i=1\mathrm{}M`$ the genes belonging to the set
* Let us denote with $`L`$ the number of base pairs (bp) of the upstream non coding region on which we shall perform our analysis. It is important to define precisely what we mean by “upstream region” With this term we denote the non coding portion of the DNA chain which is immediately before the transcription start site. This means that we do not consider as part of this region the UTR5 part of the ORF of the gene in which we are interested. If we choose $`L`$ large enough it may happen that other ORFs are present in the upstream region. In this case we consider as upstream region only the non coding part of the DNA chain up to the nearest ORF (even if it appears in the opposite strand). Thus $`L`$ should be thought of as an upper cutoff. In most cases the length of the upstream region is much smaller and is gene dependent. We shall denote it in the following as $`L(g)`$.
* In this upstream region we shall be interested in studying short sequences of nucleotides which we shall call words. Let $`n`$ be the length of such a word. For each value of $`n`$ we have $`N4^n`$ possible words $`w_i`$, $`i=1\mathrm{}N`$. The optimal choice of $`n`$ (i.e. the one which optimize the statistical significance of our analysis) is a function of $`L`$ and $`M`$. We shall see some typical values in the example of sect.5 In the following we shall have to deal with words of varying size. When needed, in order to avoid confusion, we shall call $`k`$-word a word made of $`k`$ nucleotides.

Let us call $`U`$ the collection of upstream regions of the $`M`$ genes $`g_1,\mathrm{}g_M`$. Our goal is to see if the number of occurrences of a given word $`w_i`$ in each of the upstream regions belonging to $`U`$ shows a “statistically significant” deviation (to be better defined below) from what expected on the basis of pure chance. To this end we perform two types of analyses.

First level of analysis.

This first type of analysis is organized in three steps

* Construction of the “Reference samples”. The first step is the construction of a set of $`p`$ “reference samples” which we call $`R_i,i=1,\mathrm{}p`$. The $`R_i`$ are nonoverlapping sequences of $`L_R`$ nucleotides each, extracted from a noncoding portion of the DNA sequence in the same region of the genome to which the genes that we study belong but “far” from any ORF. From these reference samples we then extract for each word the “background occurrence probability” that we shall then use as input of the second step of our analysis. The rationale behind this approach is the idea that the coding and regulating parts of the genome are immersed in a large background sea of “silent” DNA and that we may recognize that a portion of DNA has a biological function by looking at statistical deviations in the word occurrences with respect to the background. However it is clear that this is a rather crude description of the genome, in particular there are some obvious objections to this approach:

+ There is no clear notion of what “far” means. As we mentioned in the introduction one can sometimes find TF’s which keep their regulatory function even if they bind to sites which are as far as $`50kb`$ from the ORF
    + It is possible that in the reference samples the nucleotide frequencies reflect some unknown biological function thus inducing a bias in the results
    + It is not clear how should one deal with the long repeated sequences which very often appear in the genome of eukaryotes

We shall discuss below how to overcome these objections.
* Background probabilities. For each word $`w`$ we study the number of occurrences $`n(w,i)`$ in the $`i^{th}`$ sample. They will follow a Poisson distribution from which we extract the background occurrence probability of the word. This method works only if $`p`$ and $`L_R`$ are large enough with respect to the number of possible words $`N`$ (we shall see in the example below some typical values for $`p`$ and $`L_R`$). However we have checked that our results are robust with respect to different choices of these background probabilities.
* Significant words. From these probabilities we can immediately construct for each $`n`$-word the expected number of occurrences in each of the upstream sequences of $`U`$ and from them the probabilities $`p(n,s)`$ of finding at least one $`n`$-word simultaneously present in the upstream region of $`s`$ (out of the $`M`$) genes. By suitably tuning $`L`$, $`s`$ and $`n`$ we may reach very low probabilities. If notwithstanding such a low probability we indeed find a $`n`$-word which appears in the upstream region of $`s`$ genes then we consider this fact as a strong indication of its role as binding sequence for a TF’s. We may use the probability $`p(n,s)`$ as an estimate of the significance of such a candidate binding sequence.

As we have seen the critical point of this analysis is in the choice of the reference sample. We try to avoid the bias induced by this choice by crossing the above procedure with a second level of analysis

Second level of analysis.

The main change with respect to the previous analysis is that in this case we extract the reference probabilities for the $`n`$-words from an artificial reference sample constructed with a Markov chain algorithm based on the frequencies of $`k`$-words with $`k<<n`$ (usually $`k=1,2`$ or 3) extracted from the upstream regions themselves. Then the second and third step of the previous analysis follow unchanged. The rationale behind this second approach is that we want to see if in the upstream region there are some $`n`$-words (with $`n=7`$ or $`8`$, say) that occur much more often than what one would expect based on the frequency of the $`k`$-words in the same region.

These two levels of analysis are both likely to give results that are biased according to the different choices of reference probabilities that define them. However, since these biases are likely to be very different from each other, it is reasonable to expect that by comparing the results of the two methods one can minimize the number of false positives found.

## 4 The “inverse” search method.

A major drawback of the analysis discussed in the previous section is that it requires a precise knowledge of the function of the genes examined. As a matter of fact a large part of the genes of eukaryotes have no precisely know biological function and could not be studied with our direct method. Moreover in these last years the richest source of biological information on gene expression comes form microarray experiments, thus it would be important to have a tool to study gene coregulation starting from the output of such experiments. These two observations suggested us the inverse search method that we shall briefly discuss in this section. We shall outline here only the main ideas of the method, a detailed account can be found in .

The method we propose has two main steps: first the ORFs of an eukaryote genome are grouped in (overlapping) sets based on words that are overrepresented in their upstream region, with respect to their frequencies in a reference sample which is made of all the upstream regions of the whole genome. Each set is labelled by a word. Then for each of these sets the average expression in one ore more microarray experiments are compared to the genome-wide average: if a statistically significant difference is found, the word that labels the set is a candidate regulatory site for the genes in the set, either enhancing or inhibiting their expression.

An important feature is that the grouping of the genes into sets depends only on the upstream sequences and not on the microarray experiment considered: It needs to be done only once for each organism, and can then be used to analyse an arbitrary number of microarray experiments.

We refer to  for a detailed description of how the sets are constructed, we only stress here that this construction only requires three external parameters which must be fixed by the user: the length $`L`$ of the upstream region (see sect.3 for a discussion of this parameter), the length $`n`$ of the words that we use to group the sets and a cutoff probability $`P`$ which quantifies the notion of “overrepresentation” mentioned above.

## 5 Example: glycolysis in Drosophila melanogaster.

As an example of the analysis of sect.3, we studied the 7 genes of Drosophila melanogaster involved in glycolysis. These genes are listed in Tab.1.

We performed our analysis with two choices of the parameters:

1\] Promoter region. In this first test we decided to concentrate in the promoter region. Thus we chose $`L100`$. With this choice, and since $`M=7`$, we are bound to study $`n`$-words with $`n=3,4,5`$ in order to have a reasonable statistical significance. In particular we concentrate on $`n=3`$ In the first level of analysis we chose $`L_R=100`$ and $`p=1000`$ ($`p`$ is the number of reference samples). In the second level of analysis we chose $`k=1`$ ($`k`$ being the number of nucleotides of the $`k`$-words used to construct the Markov chain). We found (among few other motifs which we do not discuss here for brevity) that a statistically relevant signal is reached by the sequence GAG. This result has a clear biological interpretation since it is the binding site of an ubiquitous TF known as GAGA factor which belongs to the class of the so called “zinc finger” TF’s222The commonly assumed binding site for the GAGA factor is the sequence GAGAG, however it has been recently realized that the minimal binding sequence is actually the 3-word GAG .. We consider this finding as a good validation test of the whole procedure.

2\] large scale analysis

In this second test we chose $`L=5000`$. This allowed us to address $`n`$-words with $`n=6,7,8`$. For the reference samples we used $`L_R=5000`$ and $`p=21`$ As a result of our analysis we obtained the probabilities $`p(n,s)`$ of finding at least one $`n`$-word in the upstream region of $`s`$ out of the 7 genes that we are studying. As an example we list in Tab.2 the values of $`p(n,s)`$ for $`s=7`$ and $`n=6,7,8`$. For the Markov chain analysis we used $`k=1,2`$.

In this case we found a 7-word which appeared in the upstream region of all the seven genes: A fact that, looking at the probabilities listed in tab.3 certainly deserves more attention. The word is TTTAAAT. A survey in the literature shows that this is indeed one of the binding sequences of a TF known as “even-skipped” which is known to regulate segmentation (and also the development of certain neurons) in Drosophila. This TF has been widely studied due to its crucial role in the early stages of embryo development, but it was not directly related up to now to the regulation of glycolysis.

## 6 Example: diauxic shift in S. cerevisiae

As an example of the analysis of sect.4, we studied the so called diauxic shift, (i.e. the metabolic shift from fermentation to respiration), in S. cerevisiae the pattern of gene expression during the shift was measured with DNA microarrays techniques in Ref. . In the experiment gene expression levels were measured for virtually all the genes of at seven time-points while glucose in the medium was progressively depleted. As a result of our analysis we found 29 significant words, that can be grouped into 6 motifs (i.e. groups of similar words). Five of them correspond to known regulatory motifs (for a database of known and putative TF’s binding sites in S. cerevisiae see ref. ). In particular three of them: STRE, MIG1 and UME6 (for the meaning of these abbreviations see again ) were previously known to be involved in glucose-induced regulation process, while for the two other known motifs: PAC and RRPE this was a new result. We consider the fact of having found known regulatory motifs a strong validation of our method.

Finally we also found a new binding sequence: ATAAGGG, which we could not associate to any known regulatory motif.

## 7 Conclusions

We have proposed two new methods to extract biological information on the Transcription Factors (and more generally on the mutual interactions among genes) from the statistical distribution of oligonucleotides in the upstream region of the genes. Both are based on the notion of a “regulatory network” responsible for the various expression patterns of the genes, and aim to find common binding sites for TFs in families of coregulated genes.

* The methods can be applied both to selected sets of genes of known biological functions (direct search method) or to the genome wide microarrays experiments (inverse search method).
* They require a complete knowledge of the upstream oligonucleotide sequences and thus they can be applied for the moment only to those organisms for which the complete genoma has been sequenced.
* In the direct method, once the set of coregulated genes has been chosen, no further external input is needed. The significance criterion of our candidates binding sites only depends on the statistical distribution of oligonucleotides in the upstream region (or in nearby regions used as test samples)
* Both can be easily implemented and could be used as standard preliminary tests, to guide a more refined analysis.

Even if they already give interesting results, both our methods are far from being optimized. In particular there are three natural directions of improvement.

a\] Taking into account the variability of the binding sequences,

b\] Recognizing dyad like binding sequences (see for instance ) which are rather common in eukaryotes,

c\] Recognizing synergic interactions between TF’s.

Work is in progress along these lines. Needless to say the candidate binding sequences that we find with our method will have to be tested experimentally. However our method could help to greatly reduce the number of possible candidates and could be used as a guiding line for experiments.

## References