text:
# Single-electron soliton avalanches in tunnel junction arrays

## Abstract

Numerical modeling of correlated single-electron tunneling in uniform 2D arrays of small conducting islands separated by tunnel junctions shows the possiblility of soliton-antisoliton avalanches. Though the time duration of any avalanche, and the total charge $`\mathrm{\Delta }Q=ne`$ transferred across the array during the avalanche, are always finite, in arrays with length $`N`$ larger than certain critical value $`N_c`$ and large width $`Mf(N)`$, the avalanche magnitude $`n`$ may be exponentially large, resulting in particular in a giant increase of shot noise. Thermal fluctuations and disorder gradually suppress the avalanche effect.

PACS numbers: 72.70+m, 73.23.Hk

During the past few years, there was much interest in correlated single-electron transport  in large arrays of small conducting islands separated by tunnel barriers \- see, e.g., Refs.  and references therein. This attention is due to several interesting properties of the arrays, the most basic of which is the existence of single-electron solitons . Such a soliton consists of an additional single-electron charge placed on an island of the array, surrounded by a group of totally neutral but strongly polarized islands which screen the electron field at large distances. Similarly, removal of an electron from a single island creates an “antisoliton” (i.e., a single-hole soliton).

The concept of single-electron solitons and antisolitons and their interaction allows a natural explanation of all the peculiarities of the array statistics and dynamics, including the Kosteritz-Thouless-like phase transition in the 2D case , single-electron oscillations of frequency $`f=I/e`$ , Coulomb drag , effects of disorder , and shot noise suppression . The goal of this Letter is to report the prediction of a new phenomenon in 2D SET arrays: soliton/antisoliton avalanches which lead to shot noise enhancement rather than suppression.

We have obtained theoretical evidence of this effect during Monte Carlo simulation of 2D array dynamics within the “orthodox” model of single-electron transport . In this model, which is quantitatively valid when the tunnel conductances $`G`$ between the array islands are sufficiently low ($`Ge^2/\mathrm{}`$), single electron tunneling events are treated as incoherent transitions with rates

$$\mathrm{\Gamma }=\frac{G}{e^2}\frac{\mathrm{\Delta }W}{1\mathrm{exp}(\mathrm{\Delta }W/k_BT)},$$
(1)
where $`\mathrm{\Delta }W`$ is the drop of the electrostatic free energy caused by the particular transition. The energy $`W`$ of a charge configuration was calculated within the usual approximation  (strictly correct for an array sufficiently close to a conducting ground plane) in which the capacitance matrix includes only diagonal terms $`C_0`$ and near-diagonal terms $`C`$ which represent island self-capacitances and mutual capacitances between the neighboring islands, respectively. We studied not only the average current $`I`$ through the array, but also the spectral density $`S_I(\omega )`$ of current fluctuations.

The results show that at larger currents the Fano factor $`F=S_I(0)/2eI`$ decreases to 1/$`N`$ ($`N`$ is the longtitudinal number of junctions), just as in 1D arrays , so that the low-frequency shot noise is indeed suppressed in comparison with the Schottky value $`2eI`$. However, rather unexpectedly, we have found that in uniform arrays at low temperature $`T`$ and low currents $`I`$ (near the Coulomb blockade threshold) the Fano factor may be much larger than one - see, e.g., Fig. 1.

This fact invites the suggestion that electrons are transferred through the array in “bundles” with total charge $`\mathrm{\Delta }Qe`$; however, since single-electron solitons repel , the mechanism of the bundle formation was not immediately clear. The solution to this puzzle has turned out to be very simple. It is illustrated by Fig. 2 which shows a few snapshots of a typical Monte-Carlo-simulated process at zero temperature, when the dc voltage $`V`$ applied to the array is just slightly (in this particular case, by $`10^6e/C`$ ) above the Coulomb blockade threshold $`V_t`$. At $`VV_t`$, the soliton entrance into the initially empty array is the main transport sequence bottleneck, and takes place only after a considerable time pause when the array is free of any extra charges. After the entrance from one electrode, the soliton starts to drift, hopping along the electric field applied to the array (Fig. 2a). When the soliton approaches the opposite electrode, its field may induce the entrance of one or more antisolitons, i.e. single-hole solitons (Fig. 2b). After these antisolitons have passed their transport bottleneck at the entrance, they drift mostly along the applied electric field (in the direction opposite to that of the soliton), the attracting field of the initial soliton having little effect on the drift. Because of this, soliton-antisoliton pairs frequently miss their chance to recombine, though lateral tunnel junctions between island rows allow such recombination and we do observe such events in our simulation. In turn, each antisoliton approaching the array electrode may trigger the entrance of one or more solitons (Fig. 2c), $`\mathrm{𝑒𝑡𝑐}.`$ This chain reaction results in a soliton-antisoliton avalanche, very much similar to an electric discharge in a gas due to surface impact ionization . Notice that in 1D arrays the soliton and antisoliton cannot pass each other, so that the recombination always happens; as a result avalanches are absent and $`F1`$ .

Due to the stochastic character of transport within the framework of the orthodox theory, there is always a chance that every soliton and antisoliton leaves the array without triggering the entrance of any solitons of the opposite charge. Because of this, the avalanches always have a finite duration in time and finite “magnitude” $`n`$ (the latter may be defined by the equation $`\mathrm{\Delta }Q=ne`$ for the total charge transferred through the array during one avalanche). Nevertheless, the avalanches may be rather large: in wide arrays we have observed $`n`$’s up to $`4\times 10^6`$, limited only by the computer simulation time.

In the range where avalanches are distinct, computation of $`F`$ may be sped up considerably by using formulas for “charge blocks”, which were derived in Ref.  (in a different context):

$`S_I(0)={\displaystyle \frac{2}{\tau }}\left(e^2n^2+I^2\tau ^22eIn\tau \right),`$ (2)
$`I=en/\tau .`$ (3)
If avalanches do not overlap, $`n`$ is just the magnitude of a single avalanche, while $`\tau `$ is the interval between the beginnings of the adjacent avalanches. In the particular case $`I0`$, Eqs. (2) and (3) yield

$$F=\frac{n^2}{n}.$$
(4)
Figures 1 and 3 show the dependence of $`F`$ on the main parameters of the system: array length $`N`$, width $`M`$, island capacitance ratio $`C_0/C`$, and applied voltage $`V`$ \[the last dependence is presented parametrically via the induced dc current $`I(V)`$\]. Most of these dependences may be readily understood, at least qualitatively, using the avalanche picture discussed above. For example, if the ratio $`C_0/C`$ is increased, the soliton radius decreases , so that the soliton-antisoliton interaction is quenched, and $`F`$ is decreased. Larger $`N`$ gives the soliton more time to induce solitons of the opposite sign, so that $`F`$ grows. Finally, in very narrow arrays (small $`M`$) the soliton-antisoliton recombination suppresses the avalanche magnitude, so that the Fano factor also decreases.

We have found that the avalanche statistics dependence on $`M`$ in the opposite limit of large $`M`$ may be understood as follows. In wide arrays, the probability $`p`$ that a soliton triggers the entrance of an antisoliton, and does not recombine with it, should not depend much on the array width $`M`$. The average number $`f`$ of rows separating the soliton and the neighboring antisoliton, may depend on $`N`$ and $`C_0/C`$, but should be also virtually independent of $`M`$. This is why we may introduce an approximate “macromodel” which is schematically shown in the inset of Fig. 4: an array of width $`M`$ is presented as a parallel connection of $`m=M/f`$ channels of equal width $`f=f(N,C_0/C)`$. We break the avalanche into time stages of equal duration and assume that a passage of a soliton in some channel at stage $`i`$ triggers the antisoliton passage in each of the two neighboring channels at stage $`(i+1)`$ with probability $`p=p(N,C_0/C)`$. The lateral (open) sides of the array are described by setting the corresponding probability to 0. Mathematically, our macromodel is exactly the problem of the directed bond percolation problem on a 2D square lattice  within a stripe of width $`m`$.

Despite the approximate character of the macromodel, it allows the results of the initial problem to be reproduced remarkably well. For example, Fig. 4 shows avalanche amplitude histograms (probability to find an avalanche with particular magnitude $`n`$) for 5 arrays of different widths. The Monte Carlo simulation results virtually coincide with the macromodel results, despite only two parameters ($`p`$ and $`f`$) were available for fitting all 5 curves. (An equally good fit was obtained for all other values of $`N`$ and $`C_0/C`$ we have studied if $`M`$ was large enough, $`m1`$.) At the same time, the macromodel is much faster for simulation, so that with the same computer resources, results may be obtained in a much broader range of array parameters. For example, open points in Fig. 3 show Monte Carlo results for the initial model; for $`N=10`$, we could hardly get acceptable accuracy at $`M>40`$. For the macromodel described above, calculation of $`F`$ in that point (with the same error bars) took 35 times less CPU time, and we could continue calculations all the way up to $`M`$ = 55 (closed triangles).

The macromodel also gives a clear explanation why the growth of the avalanche magnitude with array width $`M`$ saturates for shorter arrays, but is unbound for longer arrays (Fig. 3). The directed bond percolation problem on square lattice has a percolation threshold $`p_c=0.6447`$, beyond which there is a finite probability of having an infinite percolation cluster on an infinite lattice. We have found that for our problem, $`p`$ becomes larger than $`p_c`$ if the array length is above a certain value $`N_c`$. ($`N_c8`$ for $`C_0=0`$, and grows with the increase of the $`C_0/C`$ ratio; for large enough $`C_0/C`$, $`p`$ never reaches $`p_c`$ at all – see inset in Fig. 3.) This means that at $`N>N_c`$ and unlimited $`m`$ there is a finite probability of having an infinite avalanche which, once started, would never end. The large, but finite, array width (and hence a large but finite $`m1`$) stops the avalanche growth in time, and limits its magnitude at a finite (though exponentially large in $`m`$) level.

The avalanche effect is most strongly expressed in uniform arrays at zero temperature, while both thermal fluctuations and array disorder gradually suppress it. Figure 5 shows a typical dependence of the Fano factor on the applied voltage for several values of temperature. For $`T=0`$, transport is possible only above the Coulomb blockade threshold $`V>V_t`$, so that below this point $`F`$ is undetermined. However, even very small temperature fluctuations reveal the second branch of this dependence, since the initial energy barrier for soliton entrance to the array may be overcome by thermal activation, so that dc current becomes finite (though may be very small). Figure 5 shows that for $`V<V_t`$, the Fano factor decreases as applied voltage is decreased, approaching 1 for smaller $`V`$. The reason for this decrease is that as $`V0`$, it is harder and harder for the initial soliton to trigger the antisoliton entrance, so that single-soliton passages become the dominant component of transport. (At very small voltages, $`eV/2<k_BT`$, the Fano factor starts to grow again as $`2k_BT/eV`$ due to quasi-equilibrium thermal fluctuations. However, in our case of low temperatures, $`k_BTeV_t`$, this growth corresponds to exponentially low dc current.) Larger temperature leads to gradual de-correlation of the moments when solitons and antisolitons enter the array, and hence to a gradual suppression of the avalanches which depend on this correlation.

The effect of disorder is qualitatively similar (in qualitative agreement with the results for the directed bond percolation ). For example, in a $`5\times 20`$ junction array with $`C_0=0`$, the introduction of background charges $`Q_0=0.5e`$ to just 6 percent of the islands at random usually suppresses avalanches completely, reducing $`F`$ to 1.

To summarize, we have found strong theoretical evidence for (and a simple explanation of) a new effect of single-electron-soliton and antisoliton avalanches in 2D arrays of small conducting islands. This effect arises due to soliton-assisted entrance of antisolitons into neighboring, but different rows of the array, thus avoiding recombination. The basic physics of the effect is well captured by a simple “macromodel”, equivalent to the directed percolation model on a square-lattice strip of finite width. Thermal fluctuations and array disorder suppress the avalanches. Presently, we are working toward a quantitative characterization of the latter suppression, in particular with the goal to evaluate feasibility of experimental observation of the soliton avalanches using various techniques of array fabrication . Apparently, this may be easier to accomplish in an array with tunnel-uncoupled rows (essentially, a generalization of the 2-row system studied earlier ), where tuning of dc voltages applied to each row may be sufficient to compensate for a moderate disorder.

The work was supported in part by the Engineering Research Program of the Office of Basic Energy Sciences at the Department of Energy.


format: text

text:
# Simulating Dynamical Features of Escape Panic

One of the most disastrous forms of collective human behaviour is the kind of crowd stampede induced by panic, often leading to fatalities as people are crushed or trampled. Sometimes this behaviour is triggered in life-threatening situations such as fires in crowded buildings; at other times, stampedes can arise from the rush for seats or seemingly without causes. Tragic examples within recent months include the panics in Harare, Zimbabwe, and at the Roskilde rock concert in Denmark. Although engineers are finding ways to alleviate the scale of such disasters, their frequency seems to be increasing with the number and size of mass events. Yet, systematic studies of panic behaviour,- and quantitative theories capable of predicting such crowd dynamics,- are rare. Here we show that simulations based on a model of pedestrian behaviour can provide valuable insights into the mechanisms of and preconditions for panic and jamming by incoordination. Our results suggest practical ways of minimising the harmful consequences of such events and the existence of an optimal escape strategy, corresponding to a suitable mixture of individualistic and collective behaviour.

Up to now, panics as a particular form of collective behaviour occuring in situations of scarce or dwindling resources has been mainly studied from the perspective of social psychology.- Panicking individuals tend to show maladaptive and relentless mass behaviour like jamming and life-threatening overcrowding,- which has often been attributed to social contagion (see Brown for an overview of theories). According to Mintz, the observed jamming is a result of incoordination and depends on the reward structure.

After having carefully studied the related socio-psychological literature,- reports in the media and available video materials (see http://angel.elte. hu/$`\stackrel{~}{}`$panic/), empirical investigations, and engineering handbooks, we can summarise the following characteristic features of escape panics: (i) People move or try to move considerably faster than normal. (ii) Individuals start pushing, and interactions among people become physical in nature. (iii) Moving and, in particular, passing of a bottleneck becomes incoordinated. (iv) At exits, arching and clogging are observed. (v) Jams are building up. (vi) The physical interactions in the jammed crowd add up and cause dangerous pressures up to 4,450 Newtons per meter, which can bend steel barriers or tear down brick walls. (vii) Escape is further slowed down by fallen or injured people turning into “obstacles”. (viii) People show a tendency of mass behaviour, i.e., to do what other people do. (ix) Alternative exits are often overlooked or not efficiently used in escape situations.

These observations have encouraged us to model the collective phenomenon of escape panic in the spirit of self-driven many-particle systems. Our computer simulations of the crowd dynamics of pedestrians are based on a generalised force model, which is particularly suited to describe the fatal build up of pressure observed during panics.- We assume a mixture of socio-psychological and physical forces influencing the behaviour in a crowd: Each of $`N`$ pedestrians $`i`$ of mass $`m_i`$ likes to move with a certain desired speed $`v_i^0`$ into a certain direction $`𝒆_i^0`$, and therefore tends to correspondingly adapt his or her actual velocity $`𝒗_i`$ with a certain characteristic time $`\tau _i`$. Simultaneously, he or she tries to keep a velocity-dependent distance to other pedestrians $`j`$ and walls $`W`$. This can be modelled by “interaction forces” $`𝒇_{ij}`$ and $`𝒇_{iW}`$, respectively. In mathematical terms, the change of velocity in time $`t`$ is then given by the acceleration equation

$$m_i\frac{d𝒗_i}{dt}=m_i\frac{v_i^0(t)𝒆_i^0(t)𝒗_i(t)}{\tau _i}+\underset{j(i)}{}𝒇_{ij}+\underset{W}{}𝒇_{iW},$$
(1)
while the change of position $`𝒓_i(t)`$ is given by the velocity $`𝒗_i(t)=d𝒓_i/dt`$. We describe the psychological tendency of two pedestrians $`i`$ and $`j`$ to stay away from each other by a repulsive interaction force $`A_i\mathrm{exp}[(r_{ij}d_{ij})/B_i]𝒏_{ij}`$, where $`A_i`$ and $`B_i`$ are constants. $`d_{ij}=𝒓_i𝒓_j`$ denotes the distance between the pedestrians’ centers of mass, and $`𝒏_{ij}=(n_{ij}^1,n_{ij}^2)=(𝒓_i𝒓_j)/d_{ij}`$ is the normalised vector pointing from pedestrian $`j`$ to $`i`$. If their distance $`d_{ij}`$ is smaller than the sum $`r_{ij}=(r_i+r_j)`$ of their radii $`r_i`$ and $`r_j`$, the pedestrians touch each other. In this case, we assume two additional forces inspired by granular interactions, which are essential for understanding the particular effects in panicking crowds: a “body force” $`k(r_{ij}d_{ij})𝒏_{ij}`$ counteracting body compression and a “sliding friction force” $`\kappa (r_{ij}d_{ij})\mathrm{\Delta }v_{ji}^t𝒕_{ij}`$ impeding relative tangential motion, if pedestrian $`i`$ comes close to $`j`$. Herein, $`𝒕_{ij}=(n_{ij}^2,n_{ij}^1)`$ means the tangential direction and $`\mathrm{\Delta }v_{ji}^t=(𝒗_j𝒗_i)𝒕_{ij}`$ the tangential velocity difference, while $`k`$ and $`\kappa `$ represent large constants. In summary, we have

$$𝒇_{ij}=\left\{A_i\mathrm{exp}[(r_{ij}d_{ij})/B_i]+kg(r_{ij}d_{ij})\right\}𝒏_{ij}+\kappa g(r_{ij}d_{ij})\mathrm{\Delta }v_{ji}^t𝒕_{ij},$$
(2)
where the function $`g(x)`$ is zero, if the pedestrians do not touch each other ($`d_{ij}>r_{ij}`$), otherwise equal to the argument $`x`$.

The interaction with the walls is treated analogously, i.e., if $`d_{iW}`$ means the distance to wall $`W`$, $`𝒏_{iW}`$ denotes the direction perpendicular to it, and $`𝒕_{iW}`$ the direction tangential to it, the corresponding interaction force with the wall reads

$$𝒇_{iW}=\left\{A_i\mathrm{exp}[(r_id_{iW})/B_i]+kg(r_id_{iW})\right\}𝒏_{iW}\kappa g(r_id_{iW})(𝒗_i𝒕_{iW})𝒕_{iW}.$$
(3)
Probably due to the fact that escape panics are unexpected and dangerous events, which also excludes real-life experiments, we could not find suitable data on escape panics to test our model quantitatively. This scarcity of data calls for reliable models. We have, therefore, specified the parameters as follows: With a mass of $`m_i=80`$ kg, we represent an average soccer fan. The desired velocity $`v_i^0`$ can reach more than 5 m/s (up to 10 m/s), but the observed free velocities for leaving a room correspond to $`v_i^00.6`$ m/s under relaxed, $`v_i^01`$ m/s under normal, and $`v_i^01.5`$ m/s under nervous conditions. A reasonable estimate for the acceleration time is $`\tau _i=0.5`$ s. With $`A_i=210^3`$ N and $`B_i=0.08`$ m one can reflect the distance kept at normal desired velocities and fit the measured flows through bottlenecks, amounting to 0.73 persons per second for an effectively 1 meter wide door under conditions with $`v_i^00.8`$ m/s. The parameters $`k=1.210^5`$ kg s-2 and $`\kappa =2.410^5`$ kg m-1s-1 determine the obstruction effects in cases of physical interactions. Although, in reality, most parameters are varying individually, we chose identical values for all pedestrians to minimise the number of parameters for reasons of calibration and robustness, and to exclude irregular outflows because of parameter variations. However, to avoid model artefacts (gridlocks by exactly balanced forces in symmetrical configurations), a small amount of irregularity of almost arbitrary kind is needed. This irregularity was introduced by uniformly distributed pedestrian diameters $`2r_i`$ in the interval $`[0.5\text{ m},0.7\text{ m}]`$, approximating the distribution of shoulder widths of soccer fans.

Based on the above model assumptions, we will now simulate several important phenomena of escape panic, which are insensitive to reasonable parameter variations, but fortunately become less pronounced for wider exits.

1. Transition to Incoordination due to Clogging. The simulated outflow from a room is well-coordinated and regular, if the desired velocities $`v_i^0=v_0`$ are normal. However, for desired velocities above 1.5 m/s, i.e. for people in a rush, we find an irregular succession of arch-like blockings of the exit and avalanche-like bunches of leaving pedestrians, when the arches break (see Fig. 1a, b). This phenomenon is compatible with the empirical observations mentioned above and comparable to intermittent clogging found in granular flows through funnels or hoppers (although this has been attributed to static friction between particles without remote interactions, and the transition to clogging has been observed for small enough openings rather than for a variation of the driving force).

2. “Faster-Is-Slower Effect” due to Impatience. Since clogging is connected with delays, trying to move faster (i.e., increasing $`v_i^0`$) can cause a smaller average speed of leaving, if the friction parameter $`\kappa `$ is large enough (see Fig. 1c, d). This “faster-is-slower effect” is particularly tragic in the presence of fires, where fleeing people reduce their own chances of survival. The related fatalities can be estimated by the number of pedestrians reached by the fire front (see http://angel.elte.hu/$`\stackrel{~}{}`$panic/).

Since our friction term has, on average, no deceleration effect in the crowd, if the walls are sufficiently remote, the arching underlying the clogging effect requires a combination of several effects: i) slowing down due to a bottleneck such as a door and ii) strong inter-personal friction, which becomes dominant when pedestrians get too close to each other. Consequently, the danger of clogging can be minimised by avoiding bottlenecks in the construction of stadia and public buildings. Notice, however, that jamming can also occur at widenings of escape routes! This surprising result is illustrated in Fig. 2. Improved outflows can be reached by columns placed asymmetrically in front of the exits, which also prevent the build up of fatal pressures (see http://angel.elte.hu/$`\stackrel{~}{}`$panic/).

3. Mass Behaviour. We investigate a situation in which pedestrians are trying to leave a smoky room, but first have to find one of the invisible exits (see Fig. 3a). Each pedestrian $`i`$ may either select an individual direction $`𝒆_i`$ or follow the average direction $`𝒆_j^0(t)_i`$ of his neighbours $`j`$ in a certain radius $`R_i`$, or try a mixture of both. We assume that both options are weighted with some parameter $`p_i`$:

$$𝒆_i^0(t)=𝒩\left[\left(1p_i\right)𝒆_i+p_i𝒆_j^0(t)_i\right],$$
(4)
where $`𝒩(𝒛)=𝒛/𝒛`$ denotes normalisation of a vector $`𝒛`$. As a consequence, we have individualistic behaviour if $`p_i`$ is low, but herding behaviour if $`p_i`$ is high. Therefore, $`p_i`$ reflects the degree of panics of individual $`i`$.

Our model suggests that neither individualistic nor herding behaviour performs well (see Fig. 3b). Pure individualistic behaviour means that each pedestrian finds an exit only accidentally, while pure herding behaviour implies that the complete crowd is eventually moving into the same and probably blocked direction, so that available exits are not efficiently used, in agreement with observations. According to Figs. 3b and c, we expect optimal chances of survival for a certain mixture of individualistic and herding behaviour, where individualism allows some people to detect the exits and herding guarantees that successful solutions are imitated by the others. If pedestrians follow the walls instead of “reflecting” at them, we expect that herd following causes jamming and inefficient use of doors as well (see Fig. 1), while individualists moving in opposite directions obstruct each other.

In summary, we have developed a continuous pedestrian model based on plausible interactions, which is, due to its simplicity, robust with respect to parameter variations and suitable for drawing conclusions about the possible mechanisms beyond escape panic (regarding an increase of the desired velocity, strong friction effects during physical interactions, and herding). After having calibrated the model parameters to available data on pedestrian flows, we managed to reproduce many observed phenomena including i) the build up of pressure, ii) clogging effects at bottlenecks, iii) jamming at widenings, iv) the “faster-is-slower effect”, v) inefficient use of alternative exits (see Fig. 3d), and vi) initiation of panics by counterflows and impatience (i.e., seemingly without any logical reason, see Fig. 1). Moreover, we are able to simulate situations of dwindling resources and estimate the casualties. Therefore, the model can be used to test buildings for their suitability in emergency situations. It accounts for both, the different dynamics in normal and panic situations just by changing a single parameter $`p_i=p`$. In addition, our simulations suggest that the optimal behaviour in escape situations is a suitable mixture of individualistic and herding behaviour.

We are now calling for complementary data and additional video material on escape panics to test our model quantitatively and compare it with alternative ones which, for example, include direction- and velocity-dependent interpersonal interactions, specify the individual variation of parameters, study the effect of fluctuations, consider falling people, integrate acoustic information exchange, implement more complex strategies and interactions (also three-dimensional ones), or allow for switching of strategies. A superior theory would have to reproduce the empirical findings equally well with less parameters, reach a better quantitative agreement with data with the same number, or reflect additional observations.

REFERENCES

Acknowledgments: D.H. thanks the German Research Foundation (DFG) for financial support by a Heisenberg scholarship. T.V. and I.F. are grateful for partial support by OTKA and FKFP.


text:
# Where the Slopes Are

## 1. Setting up the Problem

Let $`p`$ be a prime number, $`k2`$ an even integer, and $`N`$ a positive integer not divisible by $`p`$. Let $`\mathrm{ord}_p`$ be the $`p`$-adic valuation mapping, normalized by $`\mathrm{ord}_p(p)=1`$. For any field $`F`$ of characteristic zero, we write $`S_k(N,F)`$ to denote the $`F`$-vector space of cuspidal modular forms of weight $`k`$ for $`\mathrm{\Gamma }_0(N)`$ (with trivial character) whose Fourier coefficients all belong to $`F`$. We will essentially be concerned only with $`F=_p`$, since the Newton polygon (and therefore the slopes) can be computed already in this context, though the eigenforms themselves may only be defined over some extension of $`_p`$. Our computations will be restricted to the case $`N=1`$ (and hence $`k12`$), but it seems reasonable to set up the problem for general level.

There are two natural inclusions of $`S_k(N,F)`$ into $`S_k(Np,F)`$; on $`q`$-expansions the first is the identity mapping and the second is the Atkin-Lehner $`V`$ operator, which sends $`q`$ to $`q^p`$. The subspace spanned by the images of both maps is called the space of *oldforms* in $`S_k(Np,F)`$; it has a natural complement called the space of *newforms*.

The Atkin-Lehner $`U`$ operator maps $`S_k(Np,F)`$ to itself, acting on $`q`$-expansions by

$$U(a_nq^n)=a_{np}q^n.$$

It follows from the Atkin-Lehner theory of change of level (see \[AL70\]) that the action of $`U`$ on newforms can be diagonalized (possibly after extending the base field), and that all the eigenvalues are equal to $`\pm p^{(k2)/2}`$, and hence have slope equal to $`(k2)/2`$. Thus, as far as the slopes are concerned, the interesting questions have to do with the action of $`U`$ on the oldforms. This is best understood by relating it to the action of the Hecke operator $`\mathrm{T}_p`$ on forms of level $`N`$; this yields the theory of “twin eigenforms” discussed in \[GM92\].

The Hecke operator $`\mathrm{T}_p`$ can be diagonalized on $`S_k(N,_p)`$. Let $`fS_k(N,_p)`$ be a normalized cuspidal eigenform, and let $`a_p`$ be the eigenvalue of $`\mathrm{T}_p`$ acting on $`f`$. Finally, let $`f_1,f_2S_k(Np,_p)`$ be the two images of $`f`$ under the maps described above. The $`U`$ operator stabilizes the two-dimensional space generated by $`f_1`$ and $`f_2`$, and its characteristic polynomial is $`x^2a_px+p^{k1}`$. If this polynomial has two distinct roots, the action of $`U`$ on this two-dimensional subspace can be diagonalized, and the slopes of the two resulting eigenforms can be easily determined:

* If $`\mathrm{ord}_p(a_p)<(k1)/2`$, the two eigenvalues have $`p`$-adic valuation equal to $`\mathrm{ord}_p(a_p)`$ and $`k1\mathrm{ord}_p(a_p)`$.
* If $`\mathrm{ord}_p(a_p)(k1)/2`$, then both eigenvalues have $`p`$-adic valuation $`(k1)/2`$.

It has been conjectured by Ulmer that the polynomial $`x^2a_px+p^{k1}`$ always has two distinct roots. Specifically:

###### Conjecture  (Ulmer).

The action of $`U_p`$ on $`S_k(\mathrm{\Gamma }_0(Np),_p)`$ is semisimple. In particular, the polynomial $`x^2a_px+p^{k1}`$ always has distinct roots.

Coleman and Edixhoven have shown that this is true for $`k=2`$ and that for general $`k`$ it follows from the Tate Conjecture (see \[CE98\]).

It is easy to see that if the polynomial has a double root then we must have $`\mathrm{ord}(a_p)=(k1)/2`$ (unless $`p=2`$, in which case we would have $`\mathrm{ord}(a_p)=(k+1)/2`$). In our computations (which were restricted to the case $`N=1`$), *we never found a form for which $`\mathrm{ord}(a_p)(k1)/2`$*, so that Ulmer’s conjecture is true in every case we considered. This shows that the computations below are correct independent of Ulmer’s conjecture.

For what follows, *we assume Ulmer’s conjecture holds*. In that case, we can determine the slopes of $`U`$ acting on the oldforms in $`S_k(Np,_p)`$ by determining the slopes of $`\mathrm{T}_p`$ acting on $`S_k(N,_p)`$. For each slope obtained in level $`N`$ one obtains a pair of slopes $`\alpha ^{}`$ and $`\alpha ^{\prime \prime }`$ in level $`Np`$, satisfying

* $`0\alpha ^{}\alpha ^{\prime \prime }k1`$
* $`\alpha ^{}+\alpha ^{\prime \prime }=k1`$

with $`\alpha ^{}<\alpha ^{\prime \prime }`$ unless they are both equal to $`(k1)/2`$. We define the *slope sequence for level $`N`$, weight $`k`$, and prime $`p`$* to be the ordered list of slopes

$$(\alpha _1,\alpha _2,\mathrm{},(k1)\alpha _2,(k1)\alpha _1)$$

for $`U`$ acting on the oldforms in $`S_k(Np,_p)`$, where we repeat slopes that occur with multiplicity. The number of elements in this sequence is equal to twice the dimension of $`S_k(N,_p)`$. Since the slope sequence is symmetric under $`\alpha (k1)\alpha `$, we will usually specify it by giving only the first half of the slope sequence. The discussion above shows that this first half is the same as the slope sequence for $`\mathrm{T}_p`$ acting on $`S_k(N,_p)`$, except that all slopes larger than $`(k1)/2`$ are replaced by $`(k1)/2`$. (As mentioned above, we found no example in which the latter case occurs.)

Since we know that the slopes are in the interval $`[0,k1]`$ (and we want to vary $`k`$), it makes sense to normalize the slopes by dividing them by $`k1`$.

###### Definition .

Suppose $`f`$ is either a $`\mathrm{T}_p`$-eigenform of level $`N`$ or a $`U`$-eigenform $`f`$ of level $`Np`$. Let $`k`$ be the weight of $`f`$ and let $`a_p(f)`$ be the eigenvalue (of $`\mathrm{T}_p`$ or of $`U`$). We define the *supersingularity* of $`f`$ by

$$\mathrm{ss}(f)=\frac{\mathrm{ord}_p(a_p(f))}{k1}.$$

Let $`fS_k(N,_p)`$ is an eigenform for $`\mathrm{T}_p`$ and (still assuming Ulmer’s Conjecture is true) $`f^{},f^{\prime \prime }`$ are the two old $`u`$-eigenforms corresponding to it as above. Then, provided that $`\mathrm{ss}(f)1/2`$, we have

$$\mathrm{ss}(f^{})=\mathrm{ss}(f)$$

and

$$\mathrm{ss}(f^{\prime \prime })=1\mathrm{ss}(f),$$

and both numbers are in the interval $`[0,1]`$. Thus, the sequence of supersingularities corresponding to old eigenforms of weight $`k`$ and level $`Np`$ is a normalized version of the slope sequence, and can be computed via the supersingularities of forms of level $`N`$, provided these are small enough.

(One can think of $`\mathrm{ss}`$ as a function on the *eigencurve* studied by Coleman and Mazur in \[CM98\]. It will be an continuous function on the eigencurve, except along the $`k=1`$ locus. Notice, however, that classical eigenforms of weight $`1`$ will always have slope zero; defining $`\mathrm{ss}(f)=0`$ for such forms gives a continuous extension of $`\mathrm{ss}`$ to classical forms of weight $`1`$. No such continuous extension is possible at points corresponding to non-ordinary forms of weight $`1`$.)

We define the *supersingularity sequence in weight $`k`$*

$$(\eta _1,\eta _2,\mathrm{},1\eta _2,1\eta _1)$$

by

$$\eta _i=\mathrm{ss}(f_i)=\frac{\mathrm{slope}(f_i)}{k1}$$

as $`f_i`$ runs through the old eigenforms of weight $`k`$ on $`\mathrm{\Gamma }_0(Np)`$. The supersingularity sequence is contained in the interval $`[0,1]`$ and is symmetric under $`\eta 1\eta `$.

The main problem we consider is to understand the distribution of the supersingularities in the interval $`[0,1]`$ when we fix the level $`N`$ and let $`k\mathrm{}`$. This problem can be expressed in measure-theoretic terms, as in \[Ser97\]: considering $`N`$ and $`p`$ as fixed, for each $`k`$ we define a probability measure $`\mu _k`$ on the interval $`[0,1]`$ by putting a point mass at each supersingularity $`\eta _i`$: let $`d_k=dimS_k(N,_p)`$, and set

$$\mu _k=\frac{1}{2d_k}\underset{i=1}{\overset{d_k}{}}\left(\delta _{\eta _i}+\delta _{1\eta _i}\right),$$

where $`\delta _x`$ is the Dirac measure at $`x`$. The question then is whether the measures $`\mu _k`$ tend to a limit as $`k\mathrm{}`$, and if so to determine that limit measure.

## 2. Computations

For our computations, we restricted to the case $`N=1`$, which then means that one only gets non-trivial results for even weights $`k12`$. For each prime number $`p100`$, we computed the Newton polygon of $`\mathrm{T}_p`$ acting on forms of weight $`k`$ and level $`1`$ for weights $`k500`$. Since in every case the slopes were less than $`(k1)/2`$, the slopes we obtained are exactly the first half of the slope sequence for the $`U`$ operator acting on oldforms of level $`p`$, as described above.

The method used for computation was straightforward: the space of modular forms of weight $`k`$ and level $`1`$ has a basis consisting of forms $`E_4^aE_6^b\mathrm{\Delta }`$, where $`E_4`$ and $`E_6`$ are the Eisenstein series of weight $`4`$ and $`6`$ respectively, $`\mathrm{\Delta }`$ is the unique cuspform of weight $`12`$, and $`4a+6b+12=k`$. Using this explicit basis we determined the characteristic polynomial of $`\mathrm{T}_p`$ and computed its Newton slopes, then produced supersingularities by dividing by $`k1`$. The computation was done with the GP calculator \[BBCO\]; the basic GP functions we needed were based on a script originally written by Robert Coleman. The main constraint on the computation was the memory required for computing the characteristic polynomial: larger $`k`$ meant working with a larger basis, and larger $`p`$ meant that we needed to use more terms from the $`q`$-expansion of the modular forms. The full output of the computations can be found on the web at http://www.colby.edu/personal/fqgouvea/slopes/.

As already mentioned above, in every case we found that every slope in the Newton polygon of $`\mathrm{T}_p`$ acting on forms of level $`N`$ was smaller than $`(k1)/2`$, from which it follows that $`U`$ acts semisimply on the space of oldforms and that the slope sequence for $`\mathrm{T}_p`$ is indeed the same as the first half of the slope sequence of $`U`$ acting on oldforms of level $`Np`$.

###### Question .

Fix a prime number $`p`$ and a level $`N`$. Let $`fS_k(N,_p)`$ be an eigenform for $`\mathrm{T}_p`$. Is it true that

$$\mathrm{slope}(f_i)<\frac{k1}{2}$$

always?

In fact, one sees much more. Even a cursory observation of the tables suggests that the slopes are much smaller than one might expect. In fact, we found that *in almost every case* the supersingularities for weight $`k`$ and prime $`p`$ are smaller than $`1/(p+1)`$. In other words, the inequality

$$\mathrm{ss}(f)\frac{1}{p+1}$$

holds almost always for forms of level $`N`$. It follows that the sequence of supersingularities for weight $`k`$ is almost always contained in $`[0,\frac{1}{p+1}][\frac{p}{p+1},1]`$.

###### Question .

Fix a prime number $`p`$ and a level $`N`$. Let $`fS_k(N,_p)`$ be an eigenform for $`\mathrm{T}_p`$. Is it true that

$$\mathrm{ss}(f)\frac{1}{p+1}$$

almost always as $`k\mathrm{}`$?

To be more explicit about “almost always,” in our computations exceptions to this inequality occurred only for $`p=59`$ and $`p=79`$; for each of these primes, the inequality fails to hold for the highest-slope form in certain weights. See Table 1 for the list of weights at which exceptional slopes appear; we discuss this list of weights further below. Other exceptions to the inequality, outside the range of this computation, can be read off from the results in \[Gou97\]; they correspond to forms of weights $`k=12,16,20`$ that are non-ordinary with respect to large primes. The final entry in the table comes from a computation by Atkin. The full list of primes and weights for which we know of a slope that does not satisfy the inequality is given in Table 1. For each $`(p,k)`$ pair, we found that exactly one slope in the first half of the slope sequence violates the inequality.

If we focus on the exceptional slopes and compute the corresponding supersingularities, we see that in our examples $`\mathrm{ss}(f)`$ seems to get closer to $`1/(p+1)`$ as $`k`$ grows. Let $`p=59`$, for example; the sequence of supersingularities corresponding to the exceptional slopes in Table 1 is

$`0.066,0.022,0.026,0.019,0.022,0.018,0.020,0.017,0.019,`$
$`0.017,0.019,0.017,0.018,0.017,0.018,0.017,0.018`$
Here, of course, $`1/(p+1)=1/60=0.01666\mathrm{}`$, and the exceptional values of the supersingularity seem to be approaching this value as the weight grows.

Encouraged by this, we can try to make “almost always” precise by using the measure-theoretic formulation:

###### Question .

Is it true that the sequence of measures $`\{\mu _k\}`$ converges, as $`k\mathrm{}`$, to a measure supported on the set $`[0,\frac{1}{p+1}][\frac{p}{p+1},1]`$?

One way to think about these results is by analogy with the results in \[Gou97\]. There, we fixed a modular form and computed its slopes with respect to varying primes $`p`$, and found that the slope was almost always zero. In this case, we fix a prime and consider the full slope sequences for varying weights $`k`$, and we find that the slopes are almost always bounded by $`(k1)/(p+1)`$. The two exceptional sets are connected, of course. For example, a form of weight $`k`$ which has slope equal to $`1`$ at a prime $`p>k+2`$ will be exceptional from both points of view.

We know of no general results which suggest that “exceptional” slopes are rare. For specific primes, there are some hints. For $`p=3`$ and $`k=23^a`$, Lawren Smithline has shown in \[Smi00\] that for every form of level $`1`$ we have

$$\mathrm{slope}(f)<\frac{k}{4}.$$

For integral slopes (as we point out below, the slopes seem to be almost always integral), this is in fact equivalent to our inequality

$$\mathrm{slope}(f)\frac{k1}{4}.$$

For $`p=2`$, Kevin Buzzard has formulated a conjectural description of *all* the slopes that implies, in particular, that the inequality

$$\mathrm{slope}(f)\frac{k1}{3}$$

always holds for forms of level $`1`$, i.e., that there are no exceptional forms for $`p=2`$. (Buzzard’s conjectural description also implies that all the $`2`$-adic slopes for level $`1`$ are integral.)

These hints suggest that something is going on. It seems to us that the fact that the inequality is so often true demands some explanation. In particular, one would like to know whether there is something special about the cases where it fails.

###### Question .

Can one identify a specific property of the modular forms or the Galois representations corresponding to exceptional pairs $`(p,k)`$?

The list of exceptional cases is itself suggestive. Consider the case of $`p=59`$. For $`k=16`$, one finds a form of slope $`1`$, which is therefore an exception to the inequality. There is no clear explanation for this first exceptional slope, but the fact that the other counterexamples occur at weights $`k=46,76,106,136,\mathrm{}`$ suggests a systematic pattern. A possible interpretation of this pattern via the theory of $`\mathrm{\Theta }`$-cycles developed in \[Joc82\] was suggested by Kevin Buzzard and will be discussed in the next section.

Finally, we note that *almost all the slopes we obtained are integers*. Of all the observations we make, this is the one that is most likely to be merely an effect of the fact that we work only with small primes. The location of the exceptions, however, suggests that something else may be going on. Specifically, non-integral slopes occur in our computations only for $`p=59`$ and $`p=79`$, the same primes for which exceptional slopes occur. Furthermore, the fractional slopes we observe are connected to the exceptional slopes, in the following remarkable way:

1. A weight $`k`$ for which there exists exceptional form with slope equal to $`2`$ is preceded by a weight $`k2`$ for which the slope sequence contains two slopes equal to $`1/2`$. For $`p=59`$, this happens for the pairs of weights $`(74,76)`$ and $`(104,106)`$; for $`p=79`$, weights $`(116,118)`$, $`(122,124)`$.
2. A weight $`k`$ for which there exists exceptional form with slope equal to $`3`$ is preceded by a weight $`k2`$ for which the slope sequence contains two slopes equal to $`3/2`$, and that weight is preceded by a weight $`k4`$ whose slope sequence contains two slopes equal to $`1/2`$. This happens for $`p=59`$ and weights $`(132,134,136)`$ and $`(162,164,166)`$; for $`p=79`$, weights $`(194,196,198)`$ and $`(200,202,204)`$.
3. A weight $`k`$ for which there exists exceptional form with slope equal to $`4`$ is preceded by a weight $`k2`$ for which the slope sequence contains two slopes equal to $`5/2`$, that weight is preceded by a weight $`k4`$ whose slope sequence contains two slopes equal to $`3/2`$, and that weight is preceded by a weight $`k6`$ whose slope sequence contains two slopes equal to $`1/2`$. This happens for $`p=59`$ and weights $`(190,192,194,196)`$ and $`(220,222,224,226)`$; for $`p=79`$ and weights $`(272,274,276,278)`$ and $`(278,280,282,284)`$.
4. And so on. An exceptional form of slope $`n`$ is and weight $`k`$ is accompanied by a “train” of pairs of forms of slope

$$\frac{2n3}{2},\frac{2n5}{2},\frac{2n7}{2},\mathrm{},\frac{1}{2}$$

and weight

$$k2,k4,k6,\mathrm{},k2(n1)$$

(one pair for each weight). See the tables in section 4 for the actual slope sequences.

These patterns can overlap without interfering. For example, there are exceptional forms of slope $`6`$ for $`p=79`$ and weights $`438`$ and $`444`$. Each has its trail of weights for which fractional slopes exist. Because of the exceptional form at weight $`438`$, there are forms of slope $`1/2`$ at weight $`428`$ and at each subsequent weight, up to forms of slope $`9/2`$ at weight $`436`$. Because of the exceptional form at weight $`444`$, there are forms of slope $`1/2`$ at weight $`434`$ and at each subsequent weight, up to forms of slope $`9/2`$ at weight $`442`$. Hence, for example, at weight $`436`$ we have *both* a pair of forms of slope $`3/2`$ and a pair of forms of slope $`9/2`$.

This suggests a final question:

###### Question .

How often are the slopes integral? What is the connection between non-integral slopes and exceptional slopes?

## 3. $`\mathrm{\Theta }`$-cycles?

In this section we look a little more closely at the exceptional cases that occur for $`p=59`$ and $`p=79`$. The relevant slope sequences for $`p=59`$ and for $`p=79`$ are listed in the tables in section 4 The exceptional and the fractional slopes are printed in boldface.

The list of weights at which exceptional slopes occur strongly suggest a connection to the theory of $`\mathrm{\Theta }`$-cycles, as in \[Joc82\]. We recall the basic ideas. As above, we restrict to the case of level $`N=1`$. The $`\mathrm{\Theta }`$ operator is the operator that acts on $`q`$-expansions as $`q\frac{d}{dq}`$, so that

$$\mathrm{\Theta }\left(a_nq^n\right)=na_nq^n.$$

(If we think of modular forms as functions on the complex upper half-plane, then $`q=e^{2\pi iz}`$ and $`\mathrm{\Theta }`$ is just $`\frac{d}{dz}`$.) As is well known (see \[Gou97\]), if $`f`$ is a modular form then $`\mathrm{\Theta }f`$ is *not* a modular form, though it is “almost” modular in some sense (the $`p`$-adic story is a little different: see \[Kat73, Gou88, CGJ95\]). On the other hand, $`\mathrm{\Theta }`$ does define an operator on modular forms modulo $`p`$, in which case it maps forms of weight $`k`$ to forms of weight $`k+p+1`$ (see \[Joc82, Kat77\]).

When one considers modular forms modulo $`p`$, it is possible for forms whose weight differs by a multiple of $`p1`$ to have identical $`q`$-expansions. (Basically, this is because the $`q`$-expansion of the Eisenstein series $`E_{p1}`$ is congruent to $`1`$ modulo $`p`$.) Thus, if $`f`$ is a modular form, we define its *filtration* $`w(f)`$ to be the minimal weight for which there is a modular form whose $`q`$-expansion is congruent modulo $`p`$ to the $`q`$-expansion of $`f`$. The theory of $`\mathrm{\Theta }`$-cycles describes what happens to the filtration under the $`\mathrm{\Theta }`$ operator.

The basic facts are the following. Let $`f`$ be a modular form of weight $`k`$ on $`\mathrm{\Gamma }_0(N)`$. If $`f`$ is not ordinary (i.e., $`U(f)0(modp)`$, which is the case in all the examples we will consider), then $`\mathrm{\Theta }^{p1}ff(modp)`$, and hence $`w(\mathrm{\Theta }^{p1}f)=w(f)`$. This is why one speaks of a “$`\mathrm{\Theta }`$-*cycle*.” The filtration normally increases by $`p+1`$ each time we compute $`\mathrm{\Theta }`$; specifically, if $`w(f)`$ is not divisible by $`p`$, then $`w(\mathrm{\Theta }f)=w(f)+p+1`$. If, on the other hand, $`w(f)`$ is divisible by $`p`$, the filtration goes down, so that $`w(\mathrm{\Theta }f)=w(f)+p+1n(p1)`$, where $`n`$ is some integer. What values $`n`$ can assume is completely determined in \[Joc82\]. The case that is relevant here is the one where $`w(f)=k`$, $`4kp1`$, and $`f`$ is not ordinary. In this case, the cycle looks as follows. First,

$$w(f)=k,w(\mathrm{\Theta }f)=k+p+1,\mathrm{},w(\mathrm{\Theta }^{pk}f)=p^2(k1)p.$$

Then, since this is divisible by $`p`$, the filtration falls, with $`n=p+2k`$, so

$$w(\mathrm{\Theta }^{pk+1}f)=p+3k,\mathrm{},w(\mathrm{\Theta }^{p2}f)=(k2)p,$$

and finally $`w(\mathrm{\Theta }^{p1}f)=k`$, closing the cycle.

Notice that this whole theory refers *only* to modular forms modulo $`p`$. If $`f`$ is an eigenform, then all of the $`\mathrm{\Theta }^if`$ are eigenforms modulo $`p`$; by the Deligne-Serre Lemma (see \[DS74, Lemma 6.11\] or \[AS86, Prop. 1.2.2\]), they lift to eigenforms in characteristic zero (but recall that the $`\mathrm{\Theta }^if`$ are not themselves modular forms, so the lifts will only be congruent to them). We know these eigenforms will have positive slope (a property which is “visible” modulo $`p`$), but there seems to no reason to predict anything further about their slope.

Consider now the exceptional slopes for $`p=59`$. For the case of weight $`16`$, the occurrence of a large slope seems to be “accidental,” but its “propagation” to higher weights seems to be linked a $`\mathrm{\Theta }`$-cycle. Let $`f_0`$ be the (unique) form of weight $`16`$ and level $`1`$; its $`59`$-adic slope is $`1`$, because its $`59`$-th Fourier coefficient is divisible (once) by $`59`$. For each $`i`$, let $`f_i`$ denote the minimal weight modular form which is congruent modulo $`p`$ to $`\mathrm{\Theta }^if`$. Then we have:

* $`f_1`$ is of weight $`76`$,
* $`f_2`$ is of weight $`136`$,
* $`f_3`$ is of weight $`196`$,
*
* $`f_{43}`$ is of weight $`2596=22\times 59`$,
* $`f_{44}`$ is of weight $`46`$,
* $`f_{45}`$ is of weight $`106`$,
* $`f_{46}`$ is of weight $`166`$,
*
* $`f_{57}`$ is of weight $`826=14\times 59`$,
* finally, $`f_{58}=f_0`$.

Note that the smaller weights (up to $`496`$) on this list are precisely the weights for which we found exceptional slopes.

Because the reductions modulo $`p`$ of these forms are in the image of $`\mathrm{\Theta }`$, all of the liftings of forms in the $`\mathrm{\Theta }`$-cycle must have positive slope, but the lifting theory suggests no reason to expect the slope to increase along the cycle. Nor, in fact, is there any reason to expect the slopes to be integers. What we observe, however, is the following:

* The exceptional form of weight $`76`$ is congruent modulo $`59`$ (but not modulo $`59^2`$) to $`\mathrm{\Theta }f_0`$ *and has slope $`2`$*.
* The exceptional form of weight $`136`$ is congruent modulo $`59`$ (but not modulo $`59^2`$) to $`\mathrm{\Theta }^2f_0`$ *and has slope $`3`$*.
*
* The exceptional form of weight $`46`$ is congruent modulo $`59`$ (but not modulo $`59^2`$) to $`\mathrm{\Theta }^{44}f_0`$ *and has slope $`1`$*.
* The exceptional form of weight $`106`$ is congruent modulo $`59`$ (but not modulo $`59^2`$) to $`\mathrm{\Theta }^{45}f_0`$ *and has slope $`2`$*.
*

In other words, the forms we get are only connected via $`\mathrm{\Theta }`$ modulo $`59`$, but their slopes increase by $`1`$ as one might expect under $`\mathrm{\Theta }`$, except that at the point where the filtration falls back down to $`46`$ the slope falls back down to $`1`$. If the pattern continues to hold for the full cycle, every one of the resulting forms will be exceptional (i.e., their supersingularities will be bigger than $`1/60`$), though their supersingularities get closer and closer to $`1/60`$ as the weight gets larger.

Also interesting is to ask whether one still gets exceptional slopes beyond the $`\mathrm{\Theta }`$-cycle; for example, is there an exceptional form for $`p=59`$ of weight $`886`$? \[Computation still in progress!\]

Similarly, we can examine the exceptional slopes for $`p=79`$. Again, the exceptional slopes seem to fit into a $`\mathrm{\Theta }`$-cycle starting from a form of weight $`48`$ (the weights in the cycle would be $`48`$, $`128`$, $`208`$, …, $`2528`$, $`34`$, $`114`$, $`194`$, …, $`3634`$, $`48`$). Again, all we see is the lower end of the cycle, and again, if the pattern persists throughout the cycle every one of these forms would be exceptional.

If the existence of such “$`\mathrm{\Theta }`$-cycles” of slopes reflects a general phenomenon, then this would suggest that the set of possible slopes for a given $`p`$ and varying $`k`$ has far more structure than predicted, say, in \[GM92\]. This agrees with Kevin Buzzard’s conjectures regarding the case $`p=2`$.

## 4. Tables of Exceptional and Fractional Slopes

The tables that follow give the lower halves of the slope sequences for $`p=59`$ and $`p=79`$ and all the weights for which there is either a fractional or an exceptional slope. To make it easier for the reader, we have printed the fractional and exceptional slopes in bold. The full tables, containing all the slopes sequences for $`p100`$ and $`k500`$, can be found on the web at http://www.colby.edu/personal/fqgouvea/slopes/.

First the table for $`p=59`$:

Now the table for $`p=79`$:


text:
# SrCu2(BO3)2 – a Two Dimensional Spin Liquid

## Abstract

We study an extended Shastry-Sutherland model for SrCu2(BO3)2 and analyze the low lying parts of the energy spectrum by means of a perturbative unitary transformation based on flow equations. The derivation of the 1-magnon dispersion (elementary triplets) is discussed. Additionally, we give a quantitative description (symmetries and energies) of bound states made from two elementary triplets. Our high order results allow to fix the model parameters for SrCu2(BO3)2 precisely: $`J_1=6.16(10)`$meV, $`x:=J_2/J_1=0.603(3)`$, $`J_{}=1.3(2)`$meV. To our knowledge this is the first quantitative treatment of bound states in a true 2d model.

Presently low dimensional quantum antiferromagnets are investigated intensively both in experiment and theory. Systems that do not show a long ranged ordered ground state, so called spin liquids, are particularly interesting. Besides low spin and low coordination number spin liquids are more likely to form in strongly frustrated geometries. The recently synthesized antiferromagnetic material SrCu2(BO3)2  is a nice two dimensional example since it is a realization of the Shastry-Sutherland model  (cf. Fig. 1).

In this model frustration plays an essential role. Each spin is coupled to pairs of spins on dimers. If these pairs form singlets the couplings between dimers is without effect and the singlet-on-dimers product state is always an eigen state and for certain parameters the ground state . In this dimer phase the system is gaped. The Hamiltonian reads

$$H=J_1\underset{H_0}{\underset{}{\underset{}{}\stackrel{}{S}_i\stackrel{}{S}_j}}+J_2\underset{H_1}{\underset{}{\underset{}{}\stackrel{}{S}_i\stackrel{}{S}_k}}+J_3\underset{H_2}{\underset{}{\underset{}{}\stackrel{}{S}_i\stackrel{}{S}_l}},$$
(1)
where the bonds corresponding to interactions $`J_1`$, $`J_2`$ and $`J_3`$ are shown in Fig. 1.
For larger ratios $`x=J_2/J_1`$ and $`y=J_3/J_1`$ the system transits into a Néel phase, competing with a helical phase for certain values of $`y`$. Finally there exists a ferromagnetic phase in the region of larger negative $`x`$ and $`y`$. A profound analysis of these scenarios can be found in Ref. . Along the $`x`$-axis Koga and Kawakami  found an additional plaquette-singlet phase, which we see as beeing related to the helical phase.
In this article we start from the dimer phase where an elementary excitation is given by breaking up $`one`$ singlet and substituting a triplet instead, which acquires a dispersion by hopping from dimer to dimer. For the next higher excitations we focus on bound states formed from pairs of the elementary triplets. We will perturbatively calculate the low lying excitations of the model about the limit of strong dimerization. The problem can be stated as

$$\frac{H}{J_1}=H_0+xH_\mathrm{S},\text{with}H_\mathrm{S}=H_1+\frac{y}{x}H_2.$$
(2)
In the limit of isolated dimers ($`x=0`$, with $`y/x`$ finite) $`H`$ is bounded from below and has an equidistant energy spectrum, since $`H_0`$ simply counts the number of excited dimers (up to a trivial constant). Furthermore $`H_\mathrm{S}`$ can be decomposed

$$H_\mathrm{S}=T_1+T_0+T_1,$$
(3)
where $`T_i`$ creates $`i`$ (destroys for $`i<0`$) elementary triplets (energy quanta). These properties of $`H_0`$ and $`H_\mathrm{S}`$ allow us to use the perturbative unitary transformation  based on flow equations . This technique enables us to link smoothly and uniquely $`H`$ at $`x0`$ to an effective $`H_{\mathrm{eff}}`$ conserving the number of triplets on dimers $`[H_{\mathrm{eff}},H_0]=0`$. This permits a clear distinction between the ground state (without triplets), the 1-triplet sector, the 2-triplet sector etc.. Thus the effective Hamiltonian is block diagonal. The explicit form is given by

$$_{\mathrm{eff}}=_0+\underset{k=1}{\overset{\mathrm{}}{}}x^k\underset{|\underset{¯}{m}|=k,M(\underset{¯}{m})=0}{}C(\underset{¯}{m})T(\underset{¯}{m}),$$
(4)
where $`\underset{¯}{m}`$ is a vector of dimension $`k`$ of which the components are in $`\{\pm 1,0\}`$; $`M(\underset{¯}{m})=0`$ signifies that the sum of the components vanishes which reflects the conservation of the number of energy quanta. The operator product $`T(\underset{¯}{m})`$ is defined by $`T_{\underset{¯}{m}}=T_{m_1}T_{m_2}\mathrm{}T_{m_k}`$. The coefficients $`C(\underset{¯}{m})`$ are generally valid fractions, which we computed up to order $`k=15`$.
In terms of $`H_{\mathrm{eff}}`$ it is easy to show that a hopping, displacing an elementary triplet by a finite distance, starts in sixth order in x, while in the 2-triplet sector $`correlated`$ hopping starts in second order (see  and Refs. therein), explaining the rather flat 1-triplet dispersion in contrast to the much stronger pronounced 2-triplet dispersion as found by Kageyama et al. by INS measurements . By means of degenerate perturbation theory Momoi and Totsuka  also derive an effective Hamiltonian for the Shastry-Sutherland model. Their third order result also shows the significance of correlated hopping in the multi-triplet dynamics of the model.
Let $`|𝐫=|r_1,r_2`$ denote the state of the system with one triplet on the dimer at $`𝐫`$ and singlets on all other sites. The amplitude $`t_{𝐫^{}𝐫}^{o(𝐫)}`$ for a triplet-hopping from site $`𝐫`$ to site $`𝐫^{}`$ is then given by

$$t_{𝐫^{}𝐫}^{o(𝐫)}=𝐫^{}|_{\mathrm{eff}}|𝐫,$$
(5)
where the upper index $`o(𝐫)\{v,h\}`$ allows to distinguish whether the hopping started on a vertically oriented ($`v`$) or a horizontally oriented dimer ($`h`$). Further we choose to split the hopping amplitudes into a net part $`\overline{t}_𝐬`$ and a deviation part $`dt_𝐬`$ ($`𝐬=𝐫^{}𝐫`$)

$$t_𝐬^{o(𝐫)}=\overline{t}_𝐬+e^{iQ𝐫}dt_𝐬,$$
(6)
with $`Q=(\pi ,\pi )`$.
Since $`H_{\mathrm{eff}}`$ conserves the number of triplets one has

$$H_{\mathrm{eff}}|𝐫=\underset{𝐫^{}}{}t_𝐫^{}^{o(𝐫)}|𝐫+𝐫^{}.$$
(7)
We introduce the Fourier transformed states

$$|\sigma ,𝐤=\frac{1}{\sqrt{N}}\underset{𝐫}{}|𝐫e^{i(𝐤+\sigma Q)𝐫}$$
(8)
with the total number of dimers $`N`$, the new quantum number $`\sigma \{0,1\}`$ reflecting the sub lattice structure and $`𝐤`$ a vector in the magnetic Brillouin zone (MBZ). $`H_{\mathrm{eff}}`$ acts as a $`2\times 2`$ matrix on the states $`|\sigma ,𝐤`$ and $`|1\sigma ,𝐤`$. Its diagonalization yields

$$\omega _{1/2}(𝐤)=\underset{\omega _0(𝐤)}{\underset{}{\frac{a_0+a_1}{2}}}\pm \frac{1}{2}\sqrt{(a_0a_1)^2+4b^2}.$$
(9)
Here we have defined

$`a_\sigma `$ $`=`$ $`\left[\overline{t}_\mathrm{𝟎}+2{\displaystyle \underset{𝐫>0}{}}\overline{t}_𝐫\mathrm{cos}((𝐤+\sigma Q)𝐫)\right]\text{ and}`$ (10)
$`b`$ $`=`$ $`2{\displaystyle \underset{\genfrac{}{}{0pt}{}{𝐫>0}{r_1+r_2\text{ even}}}{}}dt_𝐫\mathrm{cos}(\mathrm{𝐤𝐫}),`$ (11)
with $`𝐫>0`$ if $`r_1>0,`$ or $`r_1=0`$ but $`r_2>0`$. Thus the 1-triplet dispersion splits into two branches. We want to point out, however, that at $`𝐤=\mathrm{𝟎}`$ and on the borders of the MBZ the two branches fall onto each other leading to a 2-fold degenerate dispersion at these points. This can be derived by showing that the square root in Eq. (9) vanishes at these points . An analogous degeneracy is noticed in the 2-triplet sector. In Ref.  we give a detailed analysis of the symmetries of the model (2d space group p4mm with underlying point group 4mm) and show that the degeneracies are due to glide line symmetries.
We calculated the amplitudes $`\overline{t}_𝐫`$ and $`dt_𝐫`$ (and therefore the dispersion) as exact polynomials in $`J_1`$, $`x`$ and $`y`$ up to and including 15th order.
Expanding the square root in Eq. (9) about the limit of vanishing $`x`$ and $`y`$ produces terms $`x^\alpha y^\beta `$ with $`\alpha +\beta 10`$. Hence the energy splitting starts in 10th order only. It is negligible for all reasonable values of $`x`$ and $`y`$.
By substituting $`y=0`$ in $`\omega _0(𝐤)`$ we verify the decimal numbers Weihong et al.  obtained previously.
For three different sets of $`J_1`$, $`x`$ and $`y`$ the 1-triplet dispersion is plotted in the inset of Fig. 2. ESR , FIR  and INS  data suggest a value of $`\omega (\mathrm{𝟎})=2.98`$meV. At finite $`𝐤`$ we have to rely on the INS measurement, which have rather large errors. We get a very good agreement. But it is not possible to fix the model parameters unambiguously from information in the 1-triplet sector alone. Thus we continue our analysis in the 2-triplet sector.
The dynamics of two triplets at large distances is governed by 1-triplet hopping. At smaller distances a 2-particle interaction occurs additionally given by $`W_{h;𝐝;𝐫,𝐝^{}}`$ ($`W_{v;𝐝;𝐫,𝐝^{}}`$) starting with one triplet on a horizontal (vertical) dimer and another at distance $`𝐝`$. The action of $`H_{\mathrm{eff}}`$ is to shift the triplets to $`𝐫`$ and to $`𝐫+𝐝^{}`$. Nothing else is possible due to triplet number conservation. Since the total spin is conserved ($`S\{0,1,2\}`$) the distances are restricted to $`𝐝,𝐝^{}>\mathrm{𝟎}`$, because the exchange parity is fixed.
The coefficients $`W`$ for $`S\{0,1,2\}`$ and the 1-triplet hopping yield the complete 2-particle dynamics. We compute $`W`$ up to $`x^{12}`$, the coefficients for the lowest-lying states even up to $`x^{14}`$.
Analogously to the 1-triplet problem we use the following basis for the 2-triplet states

$$|𝐤,𝐝,\sigma :=N^{1/2}\underset{𝐫}{}e^{\left(i(𝐤+\sigma 𝐐)(𝐫+𝐝/2)\right)}|𝐫,𝐫+𝐝,$$
(12)
where $`𝐤`$ is the conserved total momentum in the MBZ applying due to the two sub lattices; $`|𝐫,𝐫+𝐝`$ denotes the state with triplets at $`𝐫`$ and at $`𝐫+𝐝`$ ($`𝐝>0`$). Since 1-triplet hopping is of higher order than the interaction an analytic expansion for the energies of the bound states is possible. At finite order in $`x`$ only configurations contribute where the two triplets are not too far away from each other. Of course, higher orders imply larger, but still finite distances. In particular, the energies of the four states which evolve from neighboring triplets can be computed very well since their leading interaction is linear. Investigating the matrix elements shows that it is sufficient to study the distances $`𝐝\{(0,1),(1,0),(1,\pm 1)\}`$ for order 5. To $`x^{14}`$ only $`𝐝\{(1,\pm 2),(2,\pm 1),(0,2),(2,0),(2,\pm 2)\}`$ must be added. So, for given total momentum only a finite $`8\times 8`$ or $`24\times 24`$ matrix has to be analyzed. Furthermore, the elements connecting shorter distances to longer distances and the elements among longer distances do not need to be known to very high order.
We have analyzed the dispersions in $`x^5`$ of the four states bound linearly in $`x`$ in the MBZ. Fukumoto’s results are mostly confirmed . The dispersion of bound states starts only in $`x^3`$ (contrary to $`x^4`$ claimed in Ref. ). At particular points of high point group symmetry ($`(0,0)`$,$`(0,\pi )`$) the Hamiltonian splits into six blocks corresponding to different representations of the square point group 4mm. At these points the analysis up to $`x^{14}`$ is carried out. The symmetries are classified according to the irreducible representations (four 1D, one 2D) of the point group 4mm $`\mathrm{\Gamma }_1(1),\mathrm{\Gamma }_2(x^2y^2),\mathrm{\Gamma }_3(xy),\mathrm{\Gamma }_4(xy(x^2y^2)),\mathrm{\Gamma }_5(x,y)`$ where simple polynomials are given in brackets to show the transformation behavior.

The extrapolated energies are depicted in Figs. 2 ($`S=0`$) and 3 ($`S=1`$) as functions of $`x`$. The double degeneracy for $`𝐤=(0,\pi )`$ results from the same symmetry reasons as described for the 1-magnon dispersion. The dashed-dotted curve at $`(0,\pi )`$ has to be compared to the solid and the long-dashed curve to assess the dispersion of these two modes from $`\mathrm{𝟎}`$ to $`(0,\pi )`$. While for $`S=0`$ this dispersion always has the expected behavior with $`\omega (\mathrm{𝟎})<\omega ((0,\pi ))`$ the energies for $`S=1`$ are reversed for small values of $`x`$ (cf. ). Only above $`x0.55`$ the relation $`\omega (\mathrm{𝟎})<\omega ((0,\pi ))`$ holds for $`S=1`$.

For $`S=0`$, the lowest mode vanishes at the same $`x`$ as does the elementary triplet gap $`\mathrm{\Delta }`$. So no additional instability occurs for $`S=0`$. This provides evidence against the competing singlet phase as presumed in Ref. . There is, however, a salient instability for $`S=1`$ (Fig. 3) at $`x=0.63`$. This comes as a surprise since one expects in antiferromagnets binding effects to be largest for $`S=0`$ (discussed in more detail in Ref. ).
Assuming $`\mathrm{\Delta }=2.98`$meV as above and $`\omega |_{S=1}=4.7`$meV from ESR , FIR  and INS  we can determine the model parameters for SrCu2(BO3)2 precisely ($`y=0`$): $`x=0.603(3)`$ and $`J_1=6.16(10)`$meV.
The inset in Fig. 2 shows that the 1-magnon dispersion agrees very nicely with experiment for these values. Further, the energy of the $`\mathrm{\Gamma }_3`$ singlet matches the 30cm-1 peak in Raman scattering  perfectly. The $`\mathrm{\Gamma }_4`$ singlet at 25cm-1 is forbidden by symmetry, since the Raman operator is $`\mathrm{\Gamma }_3`$ at T=0 . Calculating the next $`\mathrm{\Gamma }_3,S=0`$ bound state (not shown) yields 45cm-1 in good agreement with the experimental 46cm-1 line, too.
We conclude that the 2D model (Fig. 1) explains the low-lying excitations of SrCu2(BO3)2 perfectly. Thermodynamic quantities like the susceptibility $`\chi (T)`$ require the inclusion of the interplane coupling $`J_{}`$ which is fully frustrated not changing the dimer spins . The 3D $`\chi _{3\mathrm{D}}`$ is computed from $`\chi _{2\mathrm{D}}`$ on the mean-field level as $`\chi _{3\mathrm{D}}^1=\chi _{2\mathrm{D}}^1+4J_{}`$ (inset in Fig. 3). Fits to the experiment  give $`J_{}=1.3(2)`$meV leading to a curve agreeing without flaw above 40K. Our value for $`J_{}`$ is significantly higher than the one in Ref.  due to different values of $`x`$ and $`J_1`$.
Summarizing, we presented the first quantitative description of 2-particle bound states in 2D.
An unexpected instability for the $`S=1`$ 2-triplet bound state is found at $`x0.63`$ indicating a transition to a triplet condensate probably related to the helical phase found previously . We conjecture that this transition is first order occuring at lower $`x`$ than assumed so far.
The symmetries of experimentally relevant states were determined. The reliability of the high order results allows to fix the experimental coupling constants very precisely ($`J_1=6.16(10)`$mev, $`J_2/J_1=0.603(3)`$, $`J_{}=1.3(2)`$meV). Thereby, various experiments (ESR, FIR, INS, Raman, $`\chi (T)`$) are explained consistently.
The work is supported by the DFG in SFB 341 and in SP 1073.


text:
# A Submillimeter Selected Quasar in the Field of Abell 478

## 1 Introduction

We are performing a survey of gravitationally lensed submm sources — a project aimed at studying the cosmic star formation history. Our survey covers 11 clusters of galaxies and has so far resulted in the detection of at least 25 submm sources. Here we present results on the brightest submm source detected.

## 2 Observations and Results

In the field of Abell478 we have detected with SCUBA (JCMT, Hawaii) a very bright submm source, $`F_{850\mu \mathrm{m}}=25\pm 3\mathrm{mJy}`$ and $`F_{450\mu \mathrm{m}}=63\pm 20\mathrm{mJy}`$, and possibly a fainter source of $`F_{850\mu \mathrm{m}}=9\pm 2\mathrm{mJy}`$. The errors contain the detection uncertainty and the calibration error. To identify the sources we have obtained a deep $`I`$–band image with FORS1 (VLT, Paranal). Our optical identification of the bright submm source, which we henceforth refer to as SMMJ04135+1027, is an $`I=20.5\pm 0.1\mathrm{mag}`$ point source. There are no obvious optical counterparts for the faint submm source. In Fig. 1 the $`I`$–band image is shown overlayed by the SCUBA contours.

The Galactic extinction of A478 is $`E(BV)=0.52\mathrm{mag}`$ (Schlegel et al 1998). Hence, the corrected $`I`$ magnitude is $`19.5\pm 0.1\mathrm{mag}`$ assuming a Milky Way type extinction law with $`R_V=3.1`$. All the optical observations presented in the rest of this paper have been corrected accordingly.

To determine the redshift and the nature of SMMJ04135+1027 we have obtained an optical spectrum also with FORS1. The spectrum was observed in MOS mode with resolution $`R=150`$. The spectrum exhibits the characteristica of high redshift quasars, such as broad emission lines (e.g. Ly$`\alpha `$, CIV, CIII and SiIV), Lyman Forest absorption bluewards of the Ly$`\alpha `$ emission line, and a power law continuum. The redshift has been measured using the CIII line to be $`2.83\pm 0.01`$, which is consistent with the other emission lines. The spectrum is shown in Fig. 2. If omitting the correction for Galactic extinction, the quasar appears very red. When corrected, as mentioned above, the quasar has a colour similar to optically selected quasars (Francis et al. 1992).

## 3 Discussion

### 3.1 Optical spectrum

The spectrum shown in Fig. 2 shows a number of unusual features when compared to the QSO sample of Francis et al. (1992). In the first place, Ly$`\alpha `$ is strongly suppressed. The origin of the suppression may be continuum absorption by dust, or Ly$`\alpha `$ absorption associated with the QSO. A higher resolution spectrum is needed to distinguish between these possibilities. Secondly, the CIV emission line is remarkably strong with an equivalent width of $`110\mathrm{\AA }`$.

### 3.2 The Far–IR SED

The fluxes quoted here have not been corrected for the weak magnification arising from the gravitational lensing caused by the intervening galaxy cluster. In this case, no arcs are detected in the vicinity of the quasar, hence the gravitational magnification is expected to be small and achromatic.

We have plotted the two SCUBA data points together with four known far–IR SEDs (Klaas et al 1999, Haas et al 1998, Sodroski et al 1997). The four SEDs have been redshifted to the redshift of the quasar and scaled to the quasar flux at $`\lambda _{obs}=850\mu \mathrm{m}`$. We see that the data points correspond well to the SEDs of the two starburst galaxies, Arp220 and the Antennae, and to the low–z QSO PG0050+124. The cool Milky Way type of dust can easily be ruled out. Shorter wavelength data would be needed to see if the hot dust component that is found in low–z QSOs such as PG0050+124 is also present in SMMJ04135+1027. For an Arp220 SED, the total luminosity of SMMJ04135+1027 is $`3\times 10^{13}\mathrm{L}_{}`$.

### 3.3 Contribution of AGNs to the submm samples

A number of well-studied submm sources exhibit AGN features in their optical spectra (e.g. Ivison et al. 1998). These are also typically the brightest submm sources in the samples. However, the fraction of AGNs in submm samples is still a debated issue. Optical spectroscopy of the fainter submm sources has turned out to be very challenging, therefore it will be difficult to fully characterize the submm population with optical spectroscopy alone. As shown by Fig. 3, shorter wavelength IR data may be used to find QSO–like dust emission. The most promising method for finding AGNs in submm samples is by hard X–ray observations. A recent study with Chandra and SCUBA of the relation between the resolved X-ray background and the resolved FIR background finds almost no overlap between the two populations, which suggest that SCUBA sources primarily are powered by starburst (Fabian et al. 2000). This result then suggests that the fainter submm sources contain relatively few AGNs, and that the AGNs, that are present in the submm samples, are to be found at the bright submm flux levels.


text:
# Phase-separation of binary fluids in shear flow: a numerical study

## I Introduction

Binary fluids quenched below the demixing temperature exhibit an interesting off-equilibrium phenomenology that is nowadays reasonably well understood . Tipically, after an early stage when domains of the two phases are formed, the kinetics proceeds by coarsening of these structures maintaining their morphology invariant. For viscous fluids, such as some polymer blends, hydrodynamic effects can be neglected in a preasymptotic time domain and a description of the dynamics in terms of diffusive processes alone is appropriate. This is the situation we focus on in this paper. In this case, the typical domains size $`R(t)`$ grows according to the Lifshitz-Slyozov law $`R(t)t^\alpha t^{1/3}`$ . This whole phenomenology is known to be correctly described in terms of the time dependent Ginzburg-Landau model . This continuum approach, beside being often more efficient for numerical investigations, is better suited for analytical purposes. In fact, despite an exact solution is presently lacking, successful approximate theories have been developed. Among these a prominent role is played by the large $`N`$ expansion. Indeed, when the model is generalized to a vectorial order parameter with an arbitrary number $`N`$ of components, the limit $`N\mathrm{}`$ turns out to be analytically tractable and perturbation theories around it can be applied. This special limit has proven to provide qualitatively correct informations about the kinetic evolution contributing, in some cases, to clarify the nature of the process .

When segregating systems are driven mechanically, by applying an external field or making the liquid flow, the comprehension we have of the kinetics is much poorer. This is at odd with the wide technological interest of these systems in many application areas . In the case of an applied shear flow considered here, it is known that the evolution is profoundly changed with respect to that of a static liquid . The most noticeable effect is the alignment of the domains of the two species along the flow direction, chosen in the following as the $`x`$ axis. This produces an anisotropy in the system and, in principle, one has to consider the typical domains sizes $`R_x,R_y`$ and $`R_z`$ along the coordinate axis separately. Modified growth exponents for $`R_x,R_y`$ and $`R_z`$ with respect to the case of an immobile fluid are expected and a difference $`\alpha _x\alpha _y=1`$ between the power-law exponents of $`R_x`$ and $`R_y`$ has been reported . In experiments, depending on the system considered, it has also been observed that, after an initial growth, the shear can stabilize the system into a stationary state characterized by a very large $`R_x`$ with $`R_y`$ finite, namely the phase separation process is interrupted .

Stretching of domains requires work against surface tension. Macroscopically the viscoelastic response of the fluid shows up through a non-vanishing stress tensor . On microscopic scales it has been recognized  that the strain exerted by the fluid can produce different effects and in particular the breakup and the recombination of domains. The fragmentation process is accompanied by the formation of small bubbles and more isotropic clusters. Then these structures, which grow by diffusion and join each other favoured by the flow, are successively stretched again and broken in a cyclical way. In an extended system, with many such domains, if the break-up events occur incoherently along the fluid we do not expect any appreciable effect, apart possibly from a global slowing of the growth process. Differently, if the break-up of the domains with the successive more isotropic growth is a coherent process in the network occurring when the stress exceeds a sort of threshold reached syncronously in all the system, one could expect a periodicity superimposed on the net behavior of the typical domains size and of other observables.

In a previous paper  we have shown that a pattern similar to the latter case outlined above is indeed exhibited by the large-$`N`$ model: Physical observables such as the domains radii or the rheological indicators behave for long times as power-laws decorated by log-time periodic oscillations. In that framework it was possible to identify this oscillatory pattern with a corresponding behavior of the four peaks observed in the structure factor. Moreover in the large-$`N`$ limit it was possible to show the existence of a dynamical scaling regime even in the presence of shear with the difference between the power law growth exponents in the flow and in the transverse directions equal to 1. This last observation was confirmed by a successive exact asymptotic solution of the large-$`N`$ model . In the $`t\mathrm{}`$ limit, the only accessible in this solution, however, the periodic oscillations do not survive. Then it must believed that these are long lasting preasymptotic effects. In experiments, on the other hand, a four peaked structure factor has been occasionally reported  while the so called double overshoot of the excess viscosity observed sometimes  can be interpreted as being due to the oscillations of this quantity.

Given this scenario the question of the accuracy of the large-$`N`$ limit for the description of real systems, where $`N=1`$, must be addressed. In fact, the absence of topological defects for $`N>d`$ makes concepts like domains or interfaces inappropriate and the reliability of the approximation, even qualitatively, is not for granted. Concerning this issue important topics to be addressed are, in our opinion, i) the existence of the oscillatory pattern together with its counterpart of a four-fold peaked structure factor and ii) the existence of a scaling regime with specific growth exponents. This we do in this paper by a complete numerical investigation of the time dependent Ginzburg-Landau model for sheared binary fluids. A preliminary account of our results at $`T=0`$ in $`d=2`$ was already published in a previous Letter . There the existence of structure factors with four peaks was confirmed and interpreted in terms of domains distributed in the system with two typical thicknesses. Here we overview and detail further these results. Moreover we study the effects of thermal fluctuations and the role of spatial dimensionality by comparing results for two and three dimensional systems.

Regarding issue i) our results are confirmative of the large-$`N`$ prevision: four peaks of the structure factor are resolved and an initial oscillation of the typical domains radii and other quantities is observed. Domains with different thickness and structure factors with four peaks are observed also in quenching at finite temperatures and in three-dimensional systems. With respect to point ii) despite $`R_x`$ keeps growing, CPU limitations do not allow a determination of the the behavior of $`R_y`$ and $`R_z`$ on sufficiently long timescales to detect a possible power-law behavior and the actual value of the exponents. Our data only provide a rough semi-quantitative agreement with the expected scaling laws which, therefore, cannot be inferred presently from numerical simulations alone.

This paper is organized as follows. In Sec. (II) we introduce the time dependent Ginzburg-Landau model used in the simulations. In Sec. (III) and (IV) we present the results of the numerical simulations in $`d=2`$ and $`d=3`$ respectively. Finally we summarize, discuss some open problem and draw our conclusions.

## II The model

We consider the Ginzburg-Landau free-energy

$$\{\phi \}=d^dx\{\frac{a}{2}\phi ^2+\frac{b}{4}\phi ^4+\frac{\kappa }{2}\phi ^2\}$$
(1)
where $`\phi `$ is the order parameter representing the concentration difference between the two components. The polynomial terms in the free-energy density have a single-well structure when $`a>0`$, $`b>0`$ and describe the disordered state of the mixture with $`\phi =0`$. In the ordered state, for $`a<0`$, $`b>0`$, two symmetric minima are located at $`\phi =\sqrt{a/b}`$. These are the equilibrium values of the order parameter in the small temperature limit. The gradient-squared term in (1) with $`\kappa >0`$ takes into account the energy cost of interfaces between domains of different composition. The kinetics is described by the convection-diffusion equation

$$\frac{\phi }{t}+\stackrel{}{}(\phi \stackrel{}{v})=\mathrm{\Gamma }^2\frac{\delta }{\delta \phi }+\eta $$
(2)
where the order parameter is coupled to an external velocity field $`\stackrel{}{v}`$. Here $`\mathrm{\Gamma }`$ is a mobility coefficient and $`\stackrel{}{v}`$ will be assumed to be a plane shear flow profile, namely

$$\stackrel{}{v}=\gamma y\stackrel{}{e}_x$$
(3)
$`\gamma `$ and $`\stackrel{}{e}_x`$ being a spatially uniform shear rate , and a unit vector in the $`x`$ direction, respectively. The thermal noise $`\eta `$, which describes thermal fluctuations , has zero mean and satisfies the fluctuation-dissipation relationship

$$\eta (\stackrel{}{r},t)\eta (\stackrel{}{r}^{},t^{})=2T\mathrm{\Gamma }^2\delta (\stackrel{}{r}\stackrel{}{r}^{})\delta (tt^{})$$
(4)
where $`T`$ is the temperature of the fluid and $`\mathrm{}`$ denotes the ensemble average. The validity of the present approach is restricted to systems where hydrodynamic effects can be neglected. For weakly sheared polymer blends with large polymerization index and similar mechanical properties of the two species, however, the present model is expected to be satisfactory in a preasymptotic time domain when velocity fluctuations are small .

Equation (2) can be cast in a dimensionless form after a redefinition of time, space and field scales . Then we have chosen $`\mathrm{\Gamma }=|a|=b=\kappa =1`$.

The main observable for the study of the growth kinetics is the structure factor $`C(\stackrel{}{k},t)=\phi (\stackrel{}{k},t)\phi (\stackrel{}{k},t)`$, namely the Fourier transform of the real-space equal time correlation function. From the knowledge of the structure factor one computes the average size of domains in the different directions. These quantities can be defined as

$$R_x(t)=\pi \frac{𝑑\stackrel{}{k}C(\stackrel{}{k},t)}{𝑑\stackrel{}{k}|k_x|C(\stackrel{}{k},t)}$$
(5)
and analogously for the other directions. Of experimental interest are also the rheological indicators, among which the excess viscosity defined as

$$\mathrm{\Delta }\eta =\frac{1}{\gamma }\frac{d\stackrel{}{k}}{(2\pi )^d}k_xk_yC(\stackrel{}{k},t).$$
(6)
We have simulated equation (2) in $`d=2`$ and $`d=3`$ by a first-order Euler discretization scheme. Periodic boundary conditions have been implemented in the $`x`$ and $`z`$ directions (for the $`d=3`$ case); Lees-Edwards boundary conditions  were used in the $`y`$ direction. These boundary conditions, originally developed for molecular dynamic simulation of fluids in shear, require the identification of a point at $`(x,0,z)`$ with the one located at $`(x+\gamma L\mathrm{\Delta }t,L,z)`$, where $`L`$ is the size of the lattice (the same in all the directions) and $`\mathrm{\Delta }t`$ is the time discretization interval. The system was initialized in a high temperature disordered state and the evolution was studied with $`a<0`$. Simulations were run using lattices of size $`L=1024`$ in $`d=2`$ and $`L=256`$ in $`d=3`$ with $`\mathrm{\Delta }x=0.5,1`$. The width of interfaces is given by $`\sqrt{2\kappa }`$ . We do not observe significative differences between these two choices of $`\mathrm{\Delta }x`$. The results here shown were obtained with $`\gamma =0.0488`$, $`\mathrm{\Delta }x=1`$, $`\mathrm{\Delta }t=0.01`$, $`\phi =0`$. The two dimensional system has been studied considering temperatures in the range $`0T5`$, while only the $`T=0`$ case has been studied in $`d=3`$.

## III Numerical simulations: d=2

### A T=0

In this section we discuss the results of the numerical simulations of Eq. (2) without the thermal noise $`\eta `$. A sequence of configurations of the order parameter are shown in Fig. (1) at different values of the strain $`\gamma t`$. After the usual early stage, when domains are forming from the mixed initial state, a bicontinuous structure is observed. The distortions produced by the flow appear evident from $`\gamma t1`$ onwards. At $`\gamma t=11`$ striped domains aligned with a tilt angle $`\theta `$ with respect to the flow are observed. As time elapses $`\theta `$ decreases and the structures align with the flow. In the meanwhile nonuniformities are formed in the system: Regions with domains of different thickness can be clearly observed. Small bubbles are also present, originated from the fragmentation of strained domains.

A more accurate analysis of the spatial properties can be achieved by means of the structure factor, shown in Fig. (2). Initially $`C(\stackrel{}{k},t)`$ exhibits an almost circular shape, corresponding to the early stage. Then the presence of shear deforms $`C(\stackrel{}{k},t)`$ into an ellipse and modifies the profile of the edge until, for $`\gamma t1`$, four peaks can be clearly observed. The presence of a peak in the structure factor is generally interpreted as the signature of a characteristic length in the fluid. Here, due to the anisotropy, for each peak one associates one length for each spatial dimension. Since the peaks are related by the symmetry $`\stackrel{}{k}\stackrel{}{k}`$ one concludes that there are two distinct characteristic lengths for each space direction. This corresponds to the observation of thinner and thicker domains in Fig. (1). The relative height of the peaks in one of the two specular foils in which $`C(\stackrel{}{k},t)`$ results to be separated starting from $`\gamma t4`$, is better seen in Fig. (3). For a better view of $`C(\stackrel{}{k},t)`$ we have enlarged differently the scales on the $`k_x`$ and $`k_y`$ axes. Here one observes that the two peaks dominate alternatively at the times $`\gamma t=11`$ and $`\gamma t=20`$. In a situation like that at $`\gamma t=11`$, the peak with the larger $`k_y`$ dominates, indicating that the stretched thin domains are more abundant. When the strain is further increased their rupture makes the contribution of the thick regions more important; this causes the other peak of $`C(\stackrel{}{k},t)`$ to dominate. This mechanism is very reminiscent of what happens in the large-$`N`$ model, where the recurrent prevalence of the peaks is shown to reproduce periodically in time up to much longer times then those reachable in the simulations here presented.

The mechanism of stretching and break-up which characterizes the domains evolution produces an oscillatory pattern in the typical size of ordered regions $`R_x`$ and $`R_y`$. This is shown in the first picture of Fig. (4) . $`R_y`$ reaches a local maximum when $`C(\stackrel{}{k},t)`$ is of the form of Fig. (3) at $`\gamma t=20`$. As already observed thick domains are more abundant at this time, as one can see in Fig. (1) at $`\gamma t=20`$. In the large-$`N`$ model this oscillatory behavior is periodic in the logarithm of time and is superimposed on an average power-law increase of the length $`R_x,R_y`$. Here CPU limitations prevent a clear check of both these predictions. This behavior reflects itself on most physical observables, among which the excess viscosity. This is shown in the first picture of the panel of Fig. (5) where $`\mathrm{\Delta }\eta `$ is plotted against the shear strain $`\gamma t`$. Starting from zero in correspondence of the isotropic initial condition, $`\mathrm{\Delta }\eta `$ shows a net increase up to a global maximum at $`\gamma t8`$. This is because the stretching of the domains caused by the shear require work against surface tension. This increase, however, is not monotonous and $`\mathrm{\Delta }\eta `$ shows an oscillating pattern with a first pronounced local maximum at $`\gamma t1.5`$. For larger times $`\mathrm{\Delta }\eta `$ decreases due to the ordering process in the fluid. Fig. (5) suggests that oscillations decorates the behavior of the excess viscosity even in this regime. Another oscillation is completed at $`\gamma t40`$. This is further supported by the analogous feature in the large-$`N`$ model. We mention that a double overshoot of $`\mathrm{\Delta }\eta `$ is observed in experiments too .

### B Role of temperature fluctuations

We have performed numerical simulations of Eq. (2) for several temperatures in the range $`05`$. Generally speaking we observe a similar qualitative behavior as in the zero temperature case. In particular from the analysis of the field configurations presented in Fig. (6) for $`T=5`$ one again observes the anisotropic deformation of the bicontinuous pattern emerging after the linear regime and the inhomogeneities in the domains thickness. Thermal fluctuations are responsible for the increased roughness of the interfaces with respect to $`T=0`$ and for the thermal excitations inside the ordered domains, similarly to the case without shear.

The structure factor, shown in Fig. (7), starting from $`\gamma t=4`$ exhibits four peaks whose relative heights prevail recurrently. However, the difference of the heights of the peaks on the same foil is less than in the case $`T=0`$. At $`\gamma t=1`$ instead, one has only two broad maxima along the diagonal $`k_xk_y`$ that is splitted into four peaks only at later times. The same behavior is also observed at $`T=0`$ but for earlier times. This suggests that thermal fluctuations slow the evolution of $`C(\stackrel{}{k},t)`$.

Quantitative differences with the case $`T=0`$, possibly related to the temperature-induced slower evolution, are exhibited by the characteristic sizes of the domains in the $`x`$ and $`y`$ directions. These quantities are compared for several temperatures in Fig. (4). On the basis of renormalization group  or scaling arguments  $`R_y`$ and $`R_x`$ are expected to obey an asymptotic power law growth $`t^\alpha `$ with $`\alpha _y=1/3`$ and $`\alpha _x=\alpha _y+1`$. The exponent $`\alpha _y`$ is the same as without shear and is therefore presumably caused by the evaporation-condensation mechanism responsible for growth in the absence of flow. On the other hand the growth law in the flow direction is induced by the shear. An accurate determination of these exponents is not possible on the basis of our numerical simulations. Finite-size effects become quite soon relevant in the phase separation process. Moreover the possible power-law behaviors are decorated by oscillations and it would be necessary to access much longer times in order to average over this periodic contribution. Due to finite size effects this, in turn, would require even larger lattices, which is too numerically demanding for standard workstations. Despite these limitations, by looking at Fig. (4) one can still make some qualitative observation. At $`T=0`$ one observes for $`\gamma t>30`$ a regime roughly consistent with $`\alpha _x\alpha _y=1`$ even if $`\alpha _x`$ is smaller then the expected value $`4/3`$ and $`\alpha _y`$ is not even well defined because of large oscillations. By raising the temperature one sees that this regime is gradually changed and the effective exponent $`\alpha _x`$ is lowered until, at $`T=5`$, it becomes comparable to the value $`1/3`$ observed without shear. If one believes that the asymptotic exponents are not changed by the strength of thermal fluctuations, as suggested by the renormalization group argument  above mentioned and by the analogy with the case $`\gamma =0`$, then the role of the temperature is that of producing a preasymptotic behavior characterized by a slower growth in the direction of the flow.

In binary mixtures without shear the temperature is known to help the segregation process increasing the value of the amplitude in the expression of $`R(t)At^{1/3}`$ . Here, as an effect of the increased temperature, the contribution to the phase separation of the usual evaporation-condensation mechanism, associated to the power exponent $`1/3`$, becomes more relevant than the effects of the shear induced mechanism. This results in the observed preasymptotic behaviour $`R_x=At^{1/3}`$. At later times, even if not observed in our simulations, we expect a growth with a larger exponent close to 4/3.

The above trend is reflected in the behavior of the excess viscosity shown in Fig.(5). If scaling is obeyed asymptotically, so that $`C(\stackrel{}{k},t)R_xR_yf(k_xR_x,k_yR_y)`$, then, from Eq. (6) one would expect $`\mathrm{\Delta }\eta `$ to scale as the inverse domain volume, namely $`\mathrm{\Delta }\eta R_x^1R_y^1`$. Hence, given the above discussed behavior of $`R_x,R_y`$ as the temperature is changed, one expects $`\mathrm{\Delta }\eta `$ to gradually cross-over from a power law $`t^{5/3}`$ at low temperatures to $`t^{2/3}`$ at higher temperatures, which is rougly consistent with our data.

Finally we make some observations about the influence of the temperature on the oscillating behavior. The role of the oscillations is progressively suppressed increasing $`T`$. This is observed in Fig. (5). The most apparent effect is the disapparence for $`T>1`$ of the initial double overshoot of $`\mathrm{\Delta }\eta `$. A gradual suppression of oscillations due to thermal fluctuations is also observed in the plot of $`R_x`$ and $`R_y`$, Fig. (4). This is consistent with the above discussion about the effects of raising the temperature.

## IV Numerical simulations: d=3

In this section we present the results of numerical simulations of the zero temperature quench process in a three dimensional system. The typical structures produced in the coarsening process are shown in Fig. (8). Domains appear with an interconnected cylindrical shape aligned with a tilt angle $`\theta `$ that shrinks in time with the flow direction. Contour plots of the structure factor along the coordinate planes are presented in Fig. (9). On the plane $`k_x=0`$ the shear term drops out from Eq. (2) as it can be easily checked by transforming Eq. (2) into reciprocal space. Then shear is ineffective along this direction and the structure factor displays the usual circular form as for a static fluid. In the plane $`k_y=0`$ $`C(\stackrel{}{k},t)`$ has an elliptical form with axes along $`k_x`$ and $`k_y`$. The white spots reveal that four peaks with the same height are present. In the plane $`k_z`$ a behaviour reminiscent of the two-dimensional case is observed. These results can be compared and better describe the experimental measures of the structure factor of . As in $`d=2`$, the periodic modulation of the relative heights of the different peaks produces an oscillatory behavior of most observables, as shown in Fig. (10) for the characteristic sizes $`R_x,R_y,R_z`$ and in Fig. (11) for the excess viscosity $`\mathrm{\Delta }\eta `$. While the oscillations in the $`x`$ and $`y`$ directions are out of phase, due to the recurrent prevalence of thin and thick domains described in the two-dimensional system, it is interesting to note that $`R_z`$ has a different phase with respect to $`R_y`$, indicating an even reacher behavior. These results provide at least an indication of a power-law increase of $`R_x`$ and $`R_z`$, apart from the oscillations, with an exponent larger than 1.2 for $`R_x`$. In the $`y`$ direction, however, there is not even any evidence of growth.

## V Conclusions

In this article we have considered the phase separation process of a fluid binary mixture in the presence of an applied plane shear flow. The binary mixture is described by a continuum free-energy functional and we have studied numerically the corresponding time-dependent Ginzburg-Landau equation with a convective term introduced by the flow. Results have been presented for two-dimensional systems quenched to different final temperatures and for three-dimensional fluids quenched to $`T=0`$. Simulations were carried out on lattices larger than those usually considered in phase separation studies in order to get the highest possible resolution for the structure factor.

In real binary fluids the existence of a four-peaked structure factor has been reported  together with a double overshoot of the excess viscosity. To our knowledge, however, a coherent picture indicating a possible connection between these apparently different phenomena with a clear physical interpretation of their origin was not available. The results of this paper strongly suggest that the interplay between the four peaks of the structure factor, with their recurrent prevalence, gives rise to an oscillatory phenomenon which reflects itself on the main observables. On these bases we conjecture that the double overshoot observed in experiment can be interpreted as the first part of an oscillatory pattern superimposed on the global trend of the excess viscosity and that other local maxima of this quantity could be observed in an experiment on longer time scales.

Our simulations suggest that the physical origin of the oscillating behavior is related to the presence of two kind of domains, characterized by different sizes. The kinetics proceeds by stretching the thicker and breaking the thinner ones , so that a cyclical prevalence of one of the two types occurs. This is reflected on the properties of the structure factor and other quantities. These results agree at a semi-quantitative level with the pattern observed in the large-$`N`$ approximation , indicating that the global picture put forward by the $`N=\mathrm{}`$ theory is reliable, despite the different physical character of the large-$`N`$ system without stable topological defects. This predictiveness can be used to infer that the observed oscillations are periodic in the logarithm of time, and that the phase separation proceeds via dynamical scaling, so that the characteristic sizes of the domains grow according to power-laws whose expected exponents have been computed in . However these conjectures could not be clearly addressed by our simulations with the available CPU resources, because much larger system should be considered on much longer times.

###### Acknowledgements.

F.C. is grateful to M.Cirillo, R.Del Sole and M.Palummo for hospitality in the University of Rome. F.C. and G.G. acknowledge support by the TMR network contract ERBFMRXCT980183 and by PRA-HOP 1999 INFM.

FIGURES