text:
# Fearless versus Fearful Speculative Financial Bubbles

## 1 Introduction on the detectability of bubbles

Looking back at the financial bubble of the late 1990’s, questions are now increasingly raised as to why people did not realise in time the existence of a bubble and why nobody acted early enough to prevent the maniacal increase in stock prices and the investment frenzy in technological businesses that happened during the bubble. Contradictory statements abound. Greenspan’s 1996 catchphrase, “irrational exhuberance,” interpreted to refer to an overvalued U.S. stock market is often remembered as a premonitory insight . In a private meeting, Greenspan told Lawrence Lindsey (former Fed Governor) that he saw “all sorts of parallels to the late 1920s” (quote from Investor’s Business Daily, 8/24/1999). Similarly, “The Greenspan model unveiled by the Federal Reserve chairman in 1997, calculates that the S&P is a whopping 67% above its appropriate level.” (quote of Andrew Bary, Barron’s, 1/24/2000). Robert Rubin, the chairman of the executive committee of Citigroup Inc., the biggest US financial services company, said “the world faces serious and continuing danger of another economic crisis like Asia’s in 1997 and 1998. Threats to wider world trade, the growing gap between rich and poor nations and unwise investments based on the idea that new technology will ensure prosperity all could set off such a crisis… The idea that new technologies can erase business cycles is a myth,” Rubin said. “New technologies are of profound importance, but they are not the first new technologies of significance. Autos, electricity, railroads and medicine led earlier productivity booms, and none of them, separately or together, produced one-way prosperity.” (Bloomberg News, quoted by William Fleckenstein, The Contrarian, 2/7/2000). Michel Camdessus, outgoing head of the International Monetary Fund, said “Overvalued US share prices and low US savings rates are two risks to an otherwise rosy world economic picture… you have this extreme risk of dramatically insufficient domestic savings translating into a current account deficit of a very great size… And there is indeed a risk of an overvalued financial market – by any conventional standards – which create a situation of vulnerability of which we must be mindful.” (Michel Camdessus to Janet Gittsman, Reuters, 2/8/2000). All these statements were issued before the crash on the Nasdaq in March-April 2000 and before the start of the descent of all major US markets (as well as most western markets).

In contrast, consider the following statements by A. Greenspan, issued after the burst of the last two notworthy bubbles. “We did not foresee such a breakdown in Asia. I suspect that the very nature of the process may make it virtually impossible to anticipate. It is like water pressing against a dam. Everything appears normal until a crack brings a deluge” (Fall 1998) . “The struggle to understand developments in the economy and financial markets since the mid-1990s has been particularly challenging for monetary policymakers. We were confronted with forces that none of us had personally experienced. Aside from the then recent experience of Japan, only history books and musty archives gave us clues to the appropriate stance for policy. We at the Federal Reserve considered a number of issues related to asset bubbles–that is, surges in prices of assets to unsustainable levels. As events evolved, we recognized that, despite our suspicions, it was very difficult to definitively identify a bubble until after the fact–that is, when its bursting confirmed its existence. Moreover, it was far from obvious that bubbles, even if identified early, could be preempted short of the central bank inducing a substantial contraction in economic activity–the very outcome we would be seeking to avoid” . In the same conference in August 2002, A. Greenspan added “The notion that a well-timed incremental timing could have been calibrated to prevent the late 1990’s bubble is surely an illusion.” However if money supply growth had been kept in check, interest rates would have risen, and people would have been more careful with their investment choices, which would have had a dampening effect on the frantic rising stock markets. But, before one can use money supply as a financial instrument, it requires a knowledge that one is actually in the middle of a bubble phase. In hindsight, this could sound like an easy task, but the above experiences suggest that even the most skilled professionals have a very hard time.

The situation on the detection of bubbles is perhaps even murkier in the academic literature. Empirical research has largely concentrated on testing for “explosive” exponential trends in the time series of asset prices and foreign exchange rates , with however limited success. The first reason lies in the absence of a general definition, as bubbles are model specific and generally defined from a rather restrictive framework. The concept of a fundamental price reference does not necessarily exist, nor is it necessarily unique. A major problem is that apparent evidence for bubbles can be reinterpreted in terms of market fundamentals that are unobserved by the researcher. Another suggestion is that, if stock prices are not more explosive than dividends, then it can be concluded that bubbles are not present, since bubbles are by definition the explosive component to stock prices not explained by dividends. However, periodically collapsing bubbles are not detectable by using standard tests to determine whether stock prices are more explosive or less stationary than dividends . Some recent attempts use for instance Markov-switching models to distinguish a high-return stable state from a low-return volatile state  or variance bounds tests, bubble specification tests, and cointegration tests on equity prices and dividends , or measures of deviations of future from expected returns .

Here, we present a different approach based on the idea that a bubble is more than just an exponential growth, it is a super-exponential growth. We thus use a recently introduced rational expectation model of bubbles, capturing the interplay between stochasticity and nonlinear positive feedback of prices onto future returns and volatility. We use this model to develop a methodology for the detection of bubbles. The next section recalls the formulation of the model and formulates the hypothesis and the associated null hypothesis. Section 3 describes the methodology of the test and its implementation. Section 4 describes the results obtained on 9 time series that have been previously argued to be speculative bubbles ending in crashes. Section summarizes our results and concludes.

## 2 Formulation of the nonlinear rational expectation bubble model

One can summarize the formulations of previous models of bubbles $`B(t)`$ as

$$\frac{dB(t)}{B(t)}=\mu dt+\sigma dW_t\kappa dj,$$
(1)
where $`B(t)`$ is the price above the fundamental price, $`\mu `$ is the abnormal return rate above the fundamental return, $`\sigma `$ is the volatility of the bubble and the jump term $`dj`$ describes a correction or a crash that may occur with amplitude $`\kappa `$. As we already mentioned, the problem is that (1) is of the same form as the standard geometrical Brownian motion often used to model the fundamental price. The qualification of a bubble in this context boils down to measuring an abnormal exponential rate $`\mu >0`$ of growth, a rather uncontrolled procedure since the fundamental growth rate acting as a reference is itself not known.

In Ref. , we have proposed to generalize (1) by allowing $`\mu (B(t))`$ and $`\sigma (B(t))`$ to become nonlinearly dependent on the instantaneous realization of the price. Specifically,

$`\mu (B)B`$ $`=`$ $`{\displaystyle \frac{m}{2B}}[B\sigma (B)]^2+\mu _0[B(t)/B_0]^m,`$ (2)
$`\sigma (B)B`$ $`=`$ $`\sigma _0[B(t)/B_0]^m,`$ (3)
where $`B_0`$, $`\mu _0`$, $`m>0`$ and $`\sigma _0`$ are four parameters of the model, setting respectively a reference scale, an effective drift, the strength of the nonlinearity and the magnitude of stochastic component which sets the scale of the volatility. The first term in the r.h.s. (2) is added as a convenient device to simplify the Ito calculation of these stochastic differential equations. In Ref. , we discussed several mechanisms that justify the positive feedback of $`\mu (B(t))`$ and $`\sigma (B(t))`$ on stock prices. Herding is perhaps the most obvious mechanism that leads to a positive nonlinear feedback, since the mere fact that the price of a stock rises/falls attracts more buyers/sellers respectively. Reputation herding, information cascades, the wealth effect and hedging strategies are other sources of nonlinear positive feedbacks of the type described by (2,3) (see  and references therein).

The model (1) with (2,3) is a rational expectation model of bubbles . The feature that distinguishes it from ordinary bubble models is that the solutions may be at some times growing “super-exponentially,” i.e., with a growth rate growing itself with time. This is fundamentally different from previous bubble models based on exponential growth (with constant average growth rate). This suggests to qualify bubbles by testing for the existence of nonlinear super-exponential growth. Similar ideas have been previously developed by combining the rational expectation theory of bubbles together with models of collective agent behavior  (see also  and references therein).

The model (1) with (2,3) can be reformulated in the Stratonovich interpretation given by the expression:

$$\frac{dB}{dt}=(a\mu _0+b\eta )B^m\kappa Bdj,$$
(4)
where $`a`$ and $`b`$ are two positive constants and $`\eta `$ is a delta-correlated Gaussian white noise. The form (4) examplifies the fundamental ingredient of our theory based on the interplay between nonlinearity and multiplicative noise.

The solution of (1) with (2) and (3) is derived in  and reads

$$B(t)=\alpha ^\alpha \frac{1}{\left(\mu _0[t_ct]\frac{\sigma _0}{B_0^m}W(t)\right)^\alpha },\mathrm{where}\alpha \frac{1}{m1}$$
(5)
with $`t_c=y_0/(m1)\mu _0`$ is a constant determined by the initial condition with $`y_0=1/[B(t=0)]^{m1}`$. Expression (5) is correct as long as a crash $`dj=1`$ has not occurred, which may happen at any instant according to the crash hazard rate

$$h(t)=\frac{\mu (B(t))}{\kappa },$$
(6)
determined from the no-arbitrage condition. Here, $`\kappa `$ is the average amplitude calculated over some pre-determined distribution of $`\kappa `$. When a crash occurs, $`B(t)`$ is transformed into $`B(t)(1\kappa )`$ which provides a novel initial condition which is plugged in (1) with (2) and (3).

In the deterministic case $`\sigma _0=0`$, (5) reduces to $`B(t)1/[t_ct]^{\frac{1}{m1}}`$, i.e., the bubble follows a hyperbolic growth path which would diverge in finite time, if not checked by crashes according to (6). Note that this hyperbolic growth is the signature of the positive feedback characterized by $`m>1`$ of the price $`B(t)`$ on the return rate $`\mu `$. Reintroducing the stochastic component $`\sigma _00`$, we see from (5) that a ghost finite-time singularity still exists but its visit is controlled by the first passage of a biased random walk at the position $`\mu _0t_c`$ such that the denominator $`\mu _0[t_ct]\frac{\sigma _0}{B_0^m}W(t)`$ vanishes. Due to the crash hazard rate which grows even faster than does the bubble price, this ghost-like singularity is never attained. For extensions on the concept of noisy finite-time singularities similar to (5), see Ref. .

Taking the limit $`1/\alpha 0`$ ($`m1`$) in (5) recovers the standard Black-Scholes solution

$$B_{\mathrm{BS}}(t)=\mathrm{exp}\left(\mu _0t+\sigma _0W(t)\right).$$
(7)
This is expected from (1) with (2,3) for $`m=1`$ which has the standard differential form of a geometric random walk.

## 3 Testing for nonlinear bubbles

### 3.1 Formulation of the bubble detection test

Using the two expressions (5) and (7), we propose the following testing procedure in order to qualify a given time series as a bubble. We want a test that distinguishes between a standard biased geometric Brownian motion (7), called from now on the Black-Scholes (BS) model taken as the benchmark for “normal” times, and the stochastic hyperbolic or “singular bubble” (SB) regime (5). Since the BS solution (7) is embedded in the more general SB expression (5) as the special case $`1/\alpha 0`$, testing for the presence of a bubble amounts to test whether $`1/\alpha `$ is different from $`0`$ with a sufficient high confidence. A first step is to resort to the general method of testing for embedded hypotheses and use the Likelihood Ratio test (Wilks’ chi-square statistics) .

In our empirical tests, we shall focus on time series which have been argued in the past to be bubbles because prices did appreciate by large margins before crashing precipitously over a short time interval. They are shown in figures 1-9. The table below each figure gives the studied time period, the SB and BS acceptance percentages defined below, and the values of the fit parameters for the plots shown in the figures. The crash happens right after the final date for each time series. By focusing on the part of the price trajectory before the crash, we do not need to worry about the crash jump process and do not need to address the impact of the crash hazard rate.

We assume that the observed price $`P(t)`$ evolves according to

$$\mathrm{ln}P(t)=rt+\mathrm{ln}[F]+\mathrm{ln}[B(t)].$$
(8)
Expression (8) accounts for the presence of a fundamental component $`F`$ and for the existence of an exogenous growth rate $`r`$ of the economy. This formulation (8) is different from our previous model in which $`F`$ was taken additive to $`B`$ . We propose that (8) is more appropriate since $`\mathrm{ln}[B(t)]`$ follows a geometric Brownian motion according to the null BS hypothesis and $`\mathrm{ln}(F)`$ thereby gives the constant level of the random walk proportional to $`W(t)`$.

The SB hypothesis is defined by the model (8) with (5) which has six parameters $`(\alpha ,F,r,t_c,\mu _0,\sigma _0)`$. The BS null hypothesis is specified by the model (8) with (7) and has three parameters $`(F,r+\mu _0,\sigma _0)`$. Note that the two parameters $`r`$ and $`\mu _0`$ which play distinct roles in the SB model can be fused into a single parameter $`r+\mu _0`$ for the BS model as seen from equations (7,8), making it effectively a three-parameter model.

What are we really testing when comparing the SB with the BS model? As can be seen from expression (5) and from (2), the first characteristics of a SB is its potential for faster-than-exponential accelerated growth, as we already stressed. There is another important undissociated characteristics, as shown by expression (3), that the volatility of a SB also tends to increase and to accelerate as the bubble increases. Thus, our tests for the detection of SB are probing the existence of simultaneous stochastic accelerations of the price and of the volatility. Our model has thus a behavior superficially in contradiction with the so-called leverage effect, in which volatility seems to rise when a stock price drops, and fall when the stock goes up (see the early reference  and  and references therein). The leverage effect is usually identified by a careful average over many different time series, usually over many decades of data and is thus believed to characterize the stationary behavior of prices. The leverage effect does not exclude the possible existence of isolated periods of positive correlation between past returns and future volatility. Indeed, while the leverage effect describes a short-lived correlation between past returns and future volatility, our emphasis is the cumulative trends and accelerations both in prices and volatility occurring on the longer time scales. Our SB model thus emphasizes the possibility for two anomalous behaviors occurring simultaneously: (i) super-exponential price growth and (ii) accompanying volatility growth. We shall refer to this behavior as the “fearful SB,” which is characterized by an increase of risk aversion as the bubble develops. It can thus be expected that some market phases will exhibit only one of the two, for instance only an acceleration of price without increasing volatility or even with decreasing volatility. In this case, our model may reject the SB hypothesis. The coexistence of the two characteristics (i) and (ii) can be summarized pictorially by the sentence: “Bull markets climb a wall of worry, while Bear markets decline on a slope of hope” (citation attributed to Robert Pretcher). We are going to show that this is not always the case.

### 3.2 Implementation

We describe the procedure for the SB case (the BS case follows a similar procedure). We first use a Maximum Likelihood procedure to estimate the parameters $`(\mu _0,\sigma _0)`$. From the initial value of the financial data, $`t_c`$ can in turn be determined. The remaining three parameters $`(\alpha ,F,r)`$ are then determined by a goodness-of-fit to the financial data. This improves significantly on the inversion method used in our previous work , based a Kolmogorov-Smirnov test of the Gaussianity of the increments $`dW`$ obtained by inverting the analog of (8) with an additive fundamental contribution $`F`$ to the price $`P`$. Our previous calibration of the model used only the information on the one-point statistics (distribution) of increments. We have also tried to extend this methodology by developing tests probing the residual correlation structure of the inverted increments $`dW`$ in the spirit of  but found this approach worse that our first method based on ensuring the Gaussianity of the distribution of the increments $`dW`$. Our present procedure performs even better than our first one by taking into account simultaneously the information in all multi-point statistics of the increments.

Define the variable $`X(t)`$ as:

$$X(t)=\left[P(t)e^{rt}/F\right]^{\frac{1}{\alpha }}.$$
(9)
If the observable price $`P(t)`$ was completely specified by our positive feedback model (5), the increments $`dX(t)`$ would be distributed according to a Gaussian distribution with variance $`\frac{\sigma _0}{\alpha }dt`$ and mean $`\frac{\mu _0}{\alpha }dt`$. The likelihood of $`N`$ increments is known explicitely from the model

$$L=\underset{i=1}{\overset{N}{}}p_L(dX_i),$$
(10)
where

$$p_L(dX)=\frac{1}{\sqrt{2\pi (\sigma /\alpha )^2}}\mathrm{exp}\left(\frac{\left(dX(\mu _0/\alpha )dt\right)^2}{2(\sigma _0/\alpha )^2dt}\right).$$
(11)
The maximum likelihood value of $`\mu _0`$ is thus

$$\mu _0/\alpha =\frac{1}{Ndt}\underset{1}{\overset{N}{}}dX(i),$$
(12)
where $`N`$ is the length of the time series of increments and $`T=Ndt`$ its total duration. Given fixed values of the parameters $`(\alpha ,F,r)`$, $`\mu _0`$ can then be determined by (9) and (12).

Similarly, let us define a variable $`Z`$ as

$$Z(t)=X(t)\mu _0\left(t_ct\right).$$
(13)
The probability density distribution of the increment $`dZ`$ is

$$p_L(dZ)=\frac{1}{\sqrt{2\pi (\sigma _0/\alpha )^2}}\mathrm{exp}\left(\frac{\left(dZ\right)^2}{2\left(\sigma _0/\alpha \right)^2dt}\right).$$
(14)
Constructing the likelihood as $`_{i=1}^Np_L(dZ_i)`$, its maximization with respect to $`\sigma `$ gives

$$(\sigma _0/\alpha )^2=\frac{1}{Ndt}\underset{1}{\overset{N}{}}[dZ(i)]^2.$$
(15)
Given a fixed parameter set $`(\alpha ,F,r)`$, one first determines $`\mu `$ from (12). $`\sigma `$ can then be determined by (13) and (15).

Using the initial value of the observable price $`P(t=0)`$, $`t_c`$ is then determined by:

$$t_c=\frac{\alpha }{\mu _0}[F/P(t=0)]^{1/\alpha }$$
(16)
The maximum likelihood method applied to the three other parameters $`(\alpha ,F,r)`$ leads to a set of nonlinear equations that have no analytic solutions and are not conveniently dealt with numerically. We have resorted to a non-standard method which consists literally to perform an OLS (ordinary least-square) fit of the price time series by the model (8) with (5). Specifically, we have to find the minimum in the three dimensional parameter space of $`(F,r,\alpha )`$ of the expression:

$$Q_{\mathrm{sb}}\underset{i=1}{\overset{N}{}}\left[P_{\mathrm{data}}(i)P(i)\right]^2/N.$$
(17)
In practice, this search is performed by scanning of the parameters $`(\alpha ,r)`$ and then searching for the minimum of $`Q_{\mathrm{sb}}`$ in (17) as a function of $`F`$. For each fix value of $`(F,r,\alpha )`$, $`\mu _0`$ is first calculated from (12) after which $`\sigma _0`$ can be determined from (15) and $`t_c`$ from (16). Subsequently, $`Q_{\mathrm{sb}}`$ can then be calculated from (17).

The non-standard nature of this OLS fit stems from the fact that the model $`P(i)`$ depends on the realization of the stochastic component $`W(t)`$. For each fixed values of the parameters $`(\mu _0,\sigma _0,t_c,F,r)`$, we generated $`n`$ different configurations of $`W`$. $`Q_{\mathrm{sb}}`$ was then determined as the minimum over the $`n`$ realizations. The choice $`n=5`$ is a compromise between a sufficient exploration of the space of realizations and numerical feasibility. Increasing the number $`n`$ of stochastic realizations too much is hindered by the slow convergence in the three-dimensional search space $`(\alpha ,F,r)`$ due precisely to the existence of the stochastic component of the model.

By applying the same procedure to the BS model (8) with (7), we obtain $`Q_{\mathrm{bs}}`$, defined as the minimum of the expression equivalent to (17) for the BS model over the $`n`$ stochastic realizations. In principle, we can then use Wilks’ theorem  on nested hypotheses and calculate the log-likelihood-ratio

$$T2N(\mathrm{log}\sqrt{Q_{\mathrm{sb}}}\mathrm{log}\sqrt{Q_{\mathrm{bs}}}).$$
(18)
The log-likelihood-ratio allows us to test the SB model against the BS model. A large value of $`T`$ means that the SB model gives a strong improvement over the BS model in fitting the data. To assess the statistical significance of a given value $`T`$, we can use the standard result that $`T`$ is asymptotically a $`\chi ^2`$ variable with three degrees of freedom, because the BS model has three less parameters than the SB model. The Wilks stastitics then gives the probability $`p(T)`$ that the observed value of $`T`$ can be surpassed by chance alone, given the confidence level $`1p(T)`$ at which the SB model is prefered over the null BS hypothesis.

Actually, we find that $`T`$ can be negative, that is, the BS model may give a better fit that the SB model. This is shown in the tables after each figure, where we give the $`T`$-value for the example shown in the figure and the corresponding probability for acceptance. This is a priori surprising since, by construction, the BS model is embedded in the SB model. The later cannot in principle be worse that one of its special case. In all cases, we find that the absolute value of $`T`$ is very large so as to give an acceptance of SB/BS of unity. The origins of these results are the following.

1. For numerical reasons, we cannot really take the limit $`\alpha \mathrm{}`$ but only scan up to $`\alpha _{max}=5`$. Taking larger $`\alpha `$ values more often run into overflow problems (when the dominator in the SB fit is close to $`0`$). The SB-fit can therefore never truely fit a BS model even in the deterministic case ($`W=0`$). Hence, the practical implementation of the SB model does not fully encompass the BS model. As a consequence, it is natural that the BS fit may be superior to the SB fit.
2. Due to the fluctuations in different realizations of $`W(t)`$, even exploring the limit $`\alpha \mathrm{}`$ in the SB fits, randomness of $`W`$ still makes possible that $`Q_{SB}>Q_{BS}`$.

As a consequence, a true SB model can be “mistaken” for a BS model and vice versa for a true BS model.

Consequently, we can not use the Likelihood-Ratio test to conclude reliably on what is the correct model for a given time series. The above procedure only tells us which one of the two models (BS versus SB) fits better the data, given the implementation constraints. To obtain a better sense of the statistical confidence with which one model is better than the other one for a given time series, we need an additional test that we know describe.

### 3.3 Synthetic tests and bootstrapping

In order to understand the meaning of the results presented below, it is useful to test what our procedure described in the previous section 3.2 gives for controlled cases.

As a starting point, let us take some typical parameters values $`(\alpha =2.0,F=170,r=1.0\times 10^4,t_c=5000,\mu _0=5\times 10^3,\sigma _0=1.0\times 10^4)`$ for the Dow Jones Industrial Average (DJIA<sub>1929</sub>) from 1/6 1927 to 30/9 1929, obtained by the procedure of section 3.2. Let us generate a synthetic SB time series from equations (8,5) using one specific realisation of the random process $`W(t)`$. We call this synthetic generated SB bubble for $`P_{\mathrm{SBsyn}}^1`$. We then consider $`P_{\mathrm{SBsyn}}^1`$ as an input data time series (called $`P_{\mathrm{data}}`$ in (17)) and we want to see if our detection procedure can recognize this data series as a SB bubble. We therefore launch a search for a time series $`P(t)`$ describing respectively a SB bubble (8,5) and a BS bubble (8,7). Each search for a $`P`$ described in section 3.2 depends on the initial configuration of $`W`$, and we therefore launch 10 searches with different initial configurations of $`W`$. For the case SB as well as for BS we calculate $`P_{\mathrm{SB}}^1`$, $`P_{\mathrm{BS}}^1`$ from the 10 searches with different inital configuration. Another realisation $`P_{\mathrm{SBsyn}}^2`$ of the Dow bubble is generated and again we launch 10 searches so as to calculate $`P_{\mathrm{SB}}^2`$ and $`P_{\mathrm{BS}}^2`$. Making 10 realisations of the Dow bubble, and defining $`P_{\mathrm{SB}/\mathrm{BS}}=(1/10)_{i=1}^{10}P_{\mathrm{SB}/\mathrm{BS}}^i`$ we finally get an estimate of $`P_{\mathrm{SB}}=0.71`$ and $`P_{\mathrm{BS}}=0.28`$. As explained above, we find always that one of the two models works overwhelmingly better, corresponding to a large positive or negative $`T`$-value.

This means that the SB hypothesis is not necessary to describe 28% of the realizations generated in the SB regime and the BS description appears sufficient for them. In other words, a market phase which is in a genuine SB regime by construction, characterized by strong nonlinear positive feedbacks according to the values of its parameters (especially $`1/\alpha `$), can be mistaken 28% of the time as a “normal” BS regime! Symmetrically we apply the same test in the BS regime with parameters $`(F=900,r=0.0002,\mu _0=0.001,\sigma _0=0.0001983)`$ found for the Hang Seng index from 1/6 1984 to 19/10 1987 before the crash in October 1987, we find that 37% of the BS replicas are detected as being in the BS regime and 62% are qualified as being in the SB regime. In contrast, comparing the SB and BS regimes without stochasticity ($`\sigma _0=0`$), we find that a 100% success rate in qualifying a SB as a SB and a BS as a BS.

Thus, it is clear that the failures to detect come from the stochastic character of the realizations. And this should not be too surprising: indeed, a pure BS random walk can give large spontaneous deviations that look like nonlinear stochastic accelerations. Reciprocally, a SB realization may have its $`W(t)`$ wandering such that the denominator in the right-hand-side of (5) remains far from $`0`$ with $`W(t)`$ being only a relatively minor correction to the other terms such that $`B(t)`$ has only a weak nonlinear correction compared with the standard random walk. This is related to the general problem that there is no certainty in assessing the random character of a finite string of numbers : for instance, the string head-head-head-head has exactly the same probability as head-tail-tail-head or any other for an unbiased coin toss. As the theory of algorithm complexity has taught us, saying that the later is more random that the former relies on the possibility to describe it with a shorter description or algorithm. Similarly, a SB is detected when it has a sufficiently strong linear structure that makes it unprobable to result from a pure random realization. The unfortunate but unescapable consequence is that the present stochastic rational expectation bubble formulation makes the detection of a SB not certain, with only the possibility of a probabilistic statement.

Another natural test starts with the generation of many synthetic BS models with fixed parameters obtained from the fit of the BS model to the data. Then, fitting each of these synthetic BS replicas with the SB model, we obtain a distribution of $`1/\alpha `$, which should be centered around $`0`$ and whose width should give us information on the uncertainty in recognizing a BS bubble using a SB fit. We have found no significant difference with this procedure, with a significant probability of taking a BS time series for a SB model, as found with the other approach described above.

It is also interesting to note that, whenever we get a predominant BS fit ($`P_{BS}>P_{SB}`$), the range of $`\alpha `$ for SB is restricted to relatively small values. In contrast, for predominant SB fits, we get a broad range of $`\alpha `$. Specifically, for bubbles clearly identified as BS, like e.g. HSI 1987, the range of $`\alpha `$ for the SB fits is “small” ($`0.5<\alpha <2.5`$ in that case), whereas for bubbles identified as SB like e.g. DOW 1929 the range of $`\alpha `$ is “large” ($`0.5<\alpha <5`$ in that case). This again reflects the large fluctuations inherent in the stochastic bubble model. The fact that $`\alpha `$ is smaller and clusters to smaller values for predominant BS time series should give additional credit to the belief that those markets are indeed BS-like since the well determined SB fits are inferior. In contrast, the larger values of $`\alpha `$ for the predominant SB fits give credit to the belief that those markets are more SB-like. The associated larger range of values $`\alpha `$ is a nuisance coming again from the strongly stochastic character of the bubbles and the existence of large deviations in the realizations of $`W(t)`$.

## 4 Results of empirical tests on putative bubbles

Figure 1 shows the Dow Jones Industrial Average (DJIA) Index from 1/6 1927 to 1929 (open circles), a period preceding the famous crash of October 1929, together with the best fits with the SB (solid line) and BS (dashed line) models. The parameters of the fits are given in the table. From the value of the $`T`$-statistics (18), we obtained a confidence level for the SB hypothesis of $`p(T)=1.0`$. Following the procedure described in section 3.3, consisting in making 100 searches with the DJIA as the data time series we find that 66% of the searches are qualified as SB while 33% are qualified as BS.

These two statistics confirm that the DJIA index from 1/6 1927 to 30/9 1929 exhibited a pattern consistent with a “fearful SB,” to use the terminology defined in section 3.1. The lower panel showing the volatility of the DJIA index with time confirms that the increase of the price was in fact accompanied by an overall increase in volatility.

Figure 2 shows the Nasdaq Future 100 Index from 18/6 1999 to 27/3 2000 (open circles), a period often refered to as the new economy information technology bubble, which preceded the crash in March-April 2000. The best fits of one search with the SB (solid line) and BS (dashed line) models are also shown. The parameters of the fits are given in the table. From the value of the $`T`$-statistics (18), we obtained a confidence level for the SB hypothesis of $`p(T)=0`$. As seen from the figure, the BS fit is clearly superior to the SB fit, especially at late times of the bubble, where the increasing volatility of the SB model does not match the approximately constant volatility observed for the index. Furthermore, using the procedure described in section 3.3, we find that 3% of the searches are qualified as SB while a whooping 97% are qualified as BS. The new economy information technology bubble cannot be classified as a “fearful SB.”

Figure 3 shows the S&P500 Index from 1/7 1985 to 31/8 1987 (open circles) before the crash on Black Monday, October 19, 1987, together with the best fits with the SB (solid line) and BS (dashed line) models. The parameters of the fits are given in the table. From the value of the $`T`$-statistics (18), we obtained a confidence level for the SB hypothesis of $`p(T)=1`$. Using the procedure described in section 3.3, we find that 70% of the synthetic bubbles are qualified as SB while 28% are qualified as BS. The Black monday crash can thus be classified as the end of a “fearful SB.”

Figure 4 shows the Hong Kong Hang Seng Index from 1/6 1984 to 19/10 1987 (open circles) before its crash in October, 1987, together with the best fits with the SB (solid line) and BS (dashed line) models. The parameters of the fits are given in the table. From the value of the $`T`$-statistics (18), we obtained a confidence level for the SB hypothesis of $`p(T)=0`$. Using the procedure described in section 3.3, we find that 0% of the synthetic bubbles are qualified as SB while 100% are qualified as BS. The accelerating price before this crash cannot be classified as a “fearful SB.”

Figure 5 and figure 6 give the same information as for the previous figures for the period from 3/1 1995 to 3/10 1997 preceding the crash of the Hang Seng index in October 1997 and for the period from 2/1 1991 to 4/9 1997 preceding the approximately $`7\%`$ one-day drop of the S&P500 index on October 27, 1997. From the value of the $`T`$-statistics (18), we obtain a confidence level for the SB hypothesis of $`p(T)=0`$ for the Hang Seng index and of $`p(T)=1`$ for S&P500 index. Using the procedure described in section 3.3, we find that 27% of the searches are qualified as SB while 69% are qualified as BS for the Hang Seng index, rejecting it as a “fearful SB.” For the S&P500 index, we find that 75% of the synthetic bubbles are qualified as SB while 25% are qualified as BS, describing it as a “fearful SB.”

Figure 7 is the same for the Hang Seng index from 2/1 1992 to 6/1 1994, preceding the crash on February 1994. From the value of the $`T`$-statistics (18), we obtained a confidence level for the SB hypothesis of $`p(T)=0`$. Using the procedure described in section 3.3, we find that 27% of the synthetic bubbles are qualified as SB while 73% are qualified as BS for the Hang Seng index, rejecting it as a “fearful SB.”

Figures 8 and 9 are the same for two exchange rate time series preceding two crashes, namely the US dollar expressed in Swiss Franc and in German Mark, respectively, for the period from March 1983 to March 1985. From the value of the $`T`$-statistics (18), we obtained a confidence level for the SB hypothesis of $`p(T)=1`$ for US$/CHF exchange rate and of $`p(T)=1`$ for US$/DEM echange rate. Using the procedure described in section 3.3, we find that 75% of the searches are qualified as SB while 24% are qualified as BS for the S$/CHF exchange rate, describing it as a “fearful SB.” For the US$/DEM echange rate, we find that 74% of the synthetic bubbles are qualified as SB while 26% are qualified as BS, describing it also as a “fearful SB.”

## 5 Synthesis and concluding remarks

We have studied 9 time series which all ended in a crash or severe correction. Many analysists have classified these time series as bubbles, that is, exagerated price increase over fundamental value (see for instance  and references therein). We have examined these 9 events through the lenses of our rational expectation bubble model. The two fundamental ingredients of nonlinear positive feedback and of stochasticity lead to the existence of two anomalous behaviors occurring simultaneously: (i) super-exponential price growth and (ii) volatility growth. We refer to this behavior as the “fearful singular bubble” regime. Out of the 9 time series, we find that 5 pass our tests and can be characterized as ‘fearful singular bubbles.” The 4 other cases which belong to a different class are the information technology Nasdaq bubble and the three bubbles of the Hang Seng index. Our model is not able to conclude whether or not these 4 cases are bubbles. Many previous analyses have classified these 4 cases as bubbles, but here we identify a distinguishing feature. According to our model, these four bubbles (if we accept this classification) have developed with essentially no significant increase of their volatility. They have been really impressive price appreciation without detectable “fear” or “worry” as expressed by the volatility. In contrast, the other 5 bubbles exhibited characteristics which make them recognizable in a framework emphasizing the effect of positive feedback of prices acting both on returns and volatility. This means that bubble markets do not always “climb a wall of worry,” as long as the “worry” is quantified by price volatilities. The information technology Nasdaq bubble and the three bubbles of the Hang Seng index have been characterized by a remarkable ingenuous behavior on the part of investors, apparently lacking any concern for possible associated risks or any sense of fear. In other words, the large price appreciations during the bubbles preceding the four crashes were not coming with an increasing perception of risk on the part of investors. One could tentatively conclude that the market participants were not mindful of the possibility of a crash risk looming. In contrast, the five other bubbles are more in line with a rational expectation bubble model, in which investors are at least partly aware of an increasing looming risk.

Acknowledgements: We thank Sam Salama for pointing out the references of several of the quotes in the introduction.

| Dow Jones Industrial Average from 1/6 1927 to 30/9 1929 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
| Percentage of 100 searches $`P_{SB}=0.66`$, $`P_{BS}=0.33`$ | | | | | | |
| Parameters of the curves in the figure | $`F`$ | $`r`$ | $`\sigma `$ | $`\mu `$ | $`\alpha `$ | $`T_c`$ |
| SB | 167.93 | 4.0 $`\times 10^5`$ | 6.49 $`\times 10^5`$ | 7.65 $`\times 10^4`$ | 1.5 | 1929 |
| BS | 167.47 | 3.2 $`\times 10^4`$ | 1.09 $`\times 10^4`$ | 6.82 $`\times 10^3`$ | - | - |
| T=119.42 Prob.(SB) $`>`$ 0.9999 | | | | | | |

| Nasdaq Future 100 Index from 18/6 1999 to 27/3 2000 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
| Percentage of 100 searches $`P_{SB}=0.03`$, $`P_{BS}=0.97`$ | | | | | | |
| Parameters of the curves in the figure | $`F`$ | $`r`$ | $`\sigma `$ | $`\mu `$ | $`\alpha `$ | $`T_c`$ |
| SB | 1262.08 | 2.0 $`\times 10^5`$ | 8.89 $`\times 10^5`$ | 1.55 $`\times 10^3`$ | 1.0 | 368.34 |
| BS | 1845.20 | 1.60 $`\times 10^4`$ | 5.10 $`\times 10^4`$ | 3.73 $`\times 10^3`$ | - | - |
| T=-45.30 Prob.(BS) $`>`$ 0.9999 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |

| S&P 500 Index from 1/7 1985 to 31/8 1987 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
| Percentage of 100 searches $`P_{SB}=0.70`$, $`P_{BS}=0.28`$ | | | | | | |
| Parameters of the curves in the figure | $`F`$ | $`r`$ | $`\sigma `$ | $`\mu `$ | $`\alpha `$ | $`T_c`$ |
| SB | 110.64 | 2.0 $`\times 10^5`$ | 4.59 $`\times 10^5`$ | 7.36 $`\times 10^4`$ | 3.0 | 3390.54 |
| BS | 185.06 | 2.0 $`\times 10^4`$ | 7.81 $`\times 10^5`$ | 9.65 $`\times 10^4`$ | - | - |
| T=183.71 Prob.(SB) $`>`$ 0.9999 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |

| Hang Seng Index from 1/6 1984 to 19/10 1987 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
| Percentage of 100 searches $`P_{SB}=0.01`$, $`P_{BS}=0.99`$ | | | | | | |
| Parameters of the curves in the figure | $`F`$ | $`r`$ | $`\sigma `$ | $`\mu `$ | $`\alpha `$ | $`T_c`$ |
| SB | 448.53 | 3.2 $`\times 10^4`$ | 8.17 $`\times 10^5`$ | 7.59 $`\times 10^4`$ | 2.5 | 2473.03 |
| BS | 732.50 | 6.4 $`\times 10^4`$ | 1.98 $`\times 10^4`$ | 9.14 $`\times 10^4`$ | - | - |
| T=-535.38 Prob.(BS) $`>`$ 0.9999 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |

| Hang Seng Index from 3/1 1995 to 3/10 1997 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
| Percentage of 100 searches $`P_{SB}=0.27`$, $`P_{BS}=0.69`$ | | | | | | |
| Parameters of the curves in the figure | $`F`$ | $`r`$ | $`\sigma `$ | $`\mu `$ | $`\alpha `$ | $`T_c`$ |
| SB | 7075.35 | 8.0 $`\times 10^5`$ | 9.50 $`\times 10^5`$ | 6.81 $`\times 10^4`$ | 1.5 | 2055.85 |
| BS | 8637.04 | 1.6 $`\times 10^4`$ | 1.60 $`\times 10^4`$ | 8.06 $`\times 10^3`$ | - | - |
| T=-151.21 Prob.(BS) $`>`$ 0.9999 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |

| S&P 500 Index from 2/1 1991 to 4/9 1997 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
| Percentage of 100 searches $`P_{SB}=0.75`$, $`P_{BS}=0.25`$ | | | | | | |
| Parameters of the curves in the figure | $`F`$ | $`r`$ | $`\sigma `$ | $`\mu `$ | $`\alpha `$ | $`T_c`$ |
| SB | 436.18 | 3.2 $`\times 10^4`$ | 5.18 $`\times 10^5`$ | 3.03 $`\times 10^4`$ | 5 | 17481.7 |
| BS | 329.97 | 3.2 $`\times 10^4`$ | 4.98 $`\times 10^4`$ | 3.01 $`\times 10^4`$ | - | - |
| T=696.78 Prob.(SB) $`>`$ 0.9999 | | | | | | |

| Hang Seng Index from 2/1 1992 to 6/1 1994 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
| Percentage of 100 searches $`P_{SB}=0.26`$, $`P_{BS}=0.74`$ | | | | | | |
| Parameters of the curves in the figure | $`F`$ | $`r`$ | $`\sigma `$ | $`\mu `$ | $`\alpha `$ | $`T_c`$ |
| SB | 3454.43 | 3.2 $`\times 10^4`$ | 7.30 $`\times 10^5`$ | 8.90 $`\times 10^4`$ | 1.0 | 902.14 |
| BS | 4507.98 | 3.2 $`\times 10^4`$ | 2.12 $`\times 10^4`$ | 1.62 $`\times 10^3`$ | - | - |
| T=-283.28 Prob.(BS) $`>`$ 0.9999 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |

| Swiss Francs (CHF) from 3/1 1983 to 19/3 1985 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
| Percentage of 100 searches $`P_{SB}=0.75`$, $`P_{BS}=0.24`$ | | | | | | |
| Parameters of the curves in the figure | $`F`$ | $`r`$ | $`\sigma `$ | $`\mu `$ | $`\alpha `$ | $`T_c`$ |
| SB | 2.63 | 4.0 $`\times 10^5`$ | 4.37 $`\times 10^5`$ | 5.84 $`\times 10^4`$ | 3.5 | 6486.15 |
| BS | 1.89 | 3.1.6 $`\times 10^3`$ | 4.02 $`\times 10^5`$ | 4.44 $`\times 10^4`$ | - | - |
| T=43.02 Prob.(SB) $`>`$ 0.9999 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |

| Deutch Marks (DEM) from 3/1 1983 to 8/3 1985 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
| Percentage of 100 searches $`P_{SB}=0.74`$, $`P_{BS}=0.26`$ | | | | | | |
| Parameters of the curves in the figure | $`F`$ | $`r`$ | $`\sigma `$ | $`\mu `$ | $`\alpha `$ | $`T_c`$ |
| SB | 2.17 | 3.2 $`\times 10^4`$ | 3.37 $`\times 10^5`$ | 3.16 $`\times 10^4`$ | 2.0 | 6036.89 |
| BS | 2.29 | 1.6 $`\times 10^4`$ | 3.95 $`\times 10^5`$ | 5.07 $`\times 10^3`$ | - | - |
| T=256.898987 Prob.(SB) $`>`$ 0.9999 | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |


format: text

text:
# Two-Dimensional Quantum XY Model with Ring Exchange and External Field

## Abstract

We present the zero-temperature phase diagram of a square lattice quantum spin 1/2 XY model with four-site ring exchange in a uniform external magnetic field. Using quantum Monte Carlo techniques, we identify various quantum phase transitions between the XY-order, striped or valence bond solid, staggered Néel antiferromagnet and fully polarized ground states of the model. We find no evidence for a quantum spin liquid phase.

Studies of two-dimensional spin-1/2 quantum magnet and boson models have provided insight into novel quantum phases and quantum critical points Gen1 . Recently, interest has focused on models which have multi-site ring exchange Ring1 ; Arun ; LMGH . The ring exchange interaction, either alone or in competition with the usual spin or boson near-neighbor exchange, has been shown to promote a variety of exotic quantum ground states Arun , including in some cases a spin-liquid state LMGH . Of particular importance is the class of two-dimensional model Hamiltonians that contain quantum spin-1/2 or boson operators interacting with ring exchange that can be simulated using quantum Monte Carlo (QMC) techniques without a negative sign problem. With modern algorithms, such models can be studied numerically on large lattices without approximation, providing a laboratory for surveying the critical behavior that separates various quantum phases.

One important model in this respect is the easy-plane $`J`$-$`K`$ model AD1  that has quantum $`S=1/2`$ spins on a square lattice with a near-neighbor exchange $`J`$ and a four-site ring exchange $`K`$. This Hamiltonian is partially motivated by the undoped cuprate materials CupRing , where ring-exchange processes are believed to contribute to experimental signatures beyond those explained by the near-neighbor Heisenberg model. The two-parameter $`J`$-$`K`$ model, despite its simplicity, displays a surprisingly rich and complex phase diagram AD1 ; AD2 , with three distinct zero-temperature phases. These are an XY-ordered or superfluid phase for large $`J`$, a staggered Néel or boson charge density wave (CDW) phase for large $`K`$, and a striped or valence bond solid (VBS) phase for intermediate $`K/J`$. The zero-temperature phase transition between the VBS and Néel phases is first order, however previous numerical results AD1 ; AD2  indicate the existence of a continuous quantum critical point (QCP) at the zero-temperature superfluid-VBS boundary.

The question naturally arises as to the behavior of the easy-plane $`J`$-$`K`$ model under the influence of a magnetic field. This is interesting both as a study of the evolution of the QCP, as well as the behavior of the ground state phases away from half-filling. Using stochastic series expansion (SSE) QMC, we present here the basic features of the zero-temperature phase diagram of the easy-plane model described by the Hamiltonian

$$H=J\underset{ij}{}B_{ij}K\underset{ijkl}{}P_{ijkl}h\underset{i}{}S_i^z,$$
(1)
where $`S_i^z`$ is the $`z`$ component of a quantum spin-1/2, $`B_{ij}=S_i^+S_j^{}+S_i^{}S_j^+`$ is a near-neighbor exchange, and $`P_{ijkl}=S_i^+S_j^{}S_k^+S_l^{}+S_i^{}S_j^+S_k^{}S_l^+`$ generates a four-site ring exchange. Here, $`ij`$ denotes a pair of nearest-neighbor sites and $`ijkl`$ are sites on the corners of a square plaquette on the $`L\times L`$ lattice. For $`K=0`$, this is the standard XY-model in a uniform magnetic field, or alternatively hard-core bosons with a chemical potential. For $`h=0`$, this model is in a XY-ordered or superfluid phase for temperatures less than the Kosterlitz-Thouless transition temperature of $`T_{KT}/J0.68`$ HaradaMelko . With the application of a uniform magnetic field the average magnetization $`m=S^z`$ of the XY-superfluid increases from zero ($`m=0`$) until it saturates into a fully spin polarized ($`m=1/2`$) state at $`h/J=4`$ Bern1 . For $`J=0`$ and $`h=0`$, the ground state of the system is Néel antiferromagnetic order Arun ; AD1 . For $`h=0`$, it was found AD1 ; AD2  that an intermediate VBS phase exists for $`7.9K/J14.5`$, in which the expectation value $`P_{ijkl}`$ alternates in strength with a period of two lattice spacings in one of the lattice directions, suggesting the term “striped” order.

To study the effect of the uniform magnetic field $`h`$ on the ground state properties of the easy-plane $`J`$-$`K`$ model, we use the SSE quantum Monte Carlo simulation method SSE1  that was previously applied to the $`h=0`$, $`J`$-$`K`$ model AD1 ; AD2 . In order to implement the SSE method, the operators in the Hamiltonian (1) are represented as four-spin plaquette operators. Diagonal operators involving $`h`$ terms are added to or removed from the SSE basis state expansion using a simple Metropolis probability algorithm. Off-diagonal ($`J`$ or $`K`$ term) operators are sampled using the directed-loop algorithm SSE1 ; AD2 , which becomes increasingly important for simulation efficiency with increasing magnetic field strength. The directed-loop equations SSE1  for the $`J`$-$`K`$-$`h`$ model are only slightly more complicated than for the pure $`J`$-$`K`$ model AD2 , and are presented elsewhere RAD1 . The QMC algorithms were tested on $`L=4`$ lattice sizes against exact diagonalization results and previous QMC simulations on the pure XY and $`J`$-$`K`$ models. In this paper, simulations were carried out on square lattices of linear dimension $`L`$ (number of spins $`N=L^2`$), at temperatures $`T=1/\beta `$ low enough to ensure convergence into the ground state.

A variety of physical observables of direct relevance to the ground states of the model are accessible through the SSE method. It is straightforward to calculate the internal energy SSE1  since its statistical estimator is just the number $`n`$ of plaquette operators in the SSE basis-expansion operator sequence multiplied by $`T`$: $`E=n/\beta `$. The spin stiffness (or superfluid density in the boson representation) is defined in terms of the energy response to a twist $`\varphi `$ in the periodic boundary of the lattice by

$$\rho _s=\frac{^2E}{\varphi ^2},$$
(2)
and is directly estimated using the winding number fluctuations in the SSE simulation CeperlyAWS . In addition we calculate the plaquette structure factor

$$S_p(q_x,q_y)=\frac{1}{L^2}\underset{a,b}{}\mathrm{e}^{i(𝐫_a𝐫_b)𝐪}P_{a_1a_2a_3a_4}P_{b_1b_2b_3b_4}.$$
(3)
Here, $`a_1,\mathrm{},a_4`$ are the sites belonging to plaquette $`a`$, located at $`𝐫_a`$. In the VBS phase, the square of the magnitude of the order parameter per site is $`M_P^2=[S_p(\pi ,0)+S_p(0,\pi )]/2L^2`$. Similarly, the square of the order parameter $`M_S`$ of the Néel ordered phase is obtained from the $`S^z`$ structure factor

$$S_s(q_x,q_y)=\frac{1}{L^2}\underset{j,k}{}\mathrm{e}^{i(𝐫_j𝐫_k)𝐪}S_j^zS_k^z,$$
(4)
with $`M_S^2=S_s(\pi ,\pi )/L^2`$. Here, $`j`$ and $`k`$ are lattice sites located at lattice coordinate $`𝐫_j`$. The quantities $`M_P^2`$ and $`M_S^2`$ are expected to decrease as $`1/L^2`$ (signifying short-range correlations) in phases without the respective order, but tend to a finite value for large $`L`$ in phases where long-range order occurs.

By directly observing the behavior of the spin stiffness (superfluid density) and the VBS and Néel order parameters, we are able to map out the phase boundaries of the $`J`$-$`K`$-$`h`$ model as illustrated in Fig. 1. In general, we find no persistent regions of quantum disorder (i.e. a spin liquid state) in the vicinity of the $`h=0`$ quantum critical point. Rather, the QCP appears to evolve smoothly into a quantum phase transition between the superfluid and VBS regions for $`0h6`$. The $`J`$-$`K`$-$`h`$ model also exhibits a direct superfluid to Néel order transition for $`6h11`$, a feature not contained in the $`h=0`$ phase diagram. Finally, for large $`h`$, the model finds a fully polarized spin state with $`m=1/2`$. This latter phase transition is strongly first order for $`K/J5`$, displaying pronounced metastability and hysteresis effects in the simulation (see Fig. 2). Renormalization Group treatments of two-dimensional bosons FishSub , as well as spin-wave corrected mean-field theory and simulations of a hard-core boson Hamiltonian Bern1  indicate that at $`K=0`$, the pure XY model exhibits a continuous transition to the fully polarized state at $`h=4J`$. This suggests that a tricritical point (TCP in Fig. 1) exists on the phase boundary somewhere between $`0<K/J5`$, above which the transition to the fully polarized state becomes first order.

The energy crossover and density hysteresis of Fig. 2 provide one indicator of a first-order transition. Alternatively, one may look for an abrupt discontinuity in the order parameter (for large system sizes) or for double-peaked probability histograms for data in the transition region. To illustrate this we turn now to a detailed set of simulation results for the superfluid-Néel phase boundary along cut (B) in Fig. 1. As illustrated in Fig. 3, the boson and superfluid densities develop significant discontinuities for larger systems as the phase boundary is traversed. This abrupt discontinuity does not appear for $`L<20`$, illustrating that the transition is caused by an avoiding level crossing and that large lattices sizes are necessary to quantify the behavior of this model. The first-order nature is apparent in double-peaked magnetization histograms, that were observed for data in the “discontinuity” regions for $`L=16`$, indicating a phase coexistence. For $`K/J16.3`$, the spin spin structure factor (Eq. (4)) develops Bragg peaks at $`(\pi ,\pi )`$ (not illustrated), indicating Néel order. It is interesting to note that a similar phase transition between a superfluid and a $`(\pi ,\pi )`$ staggered solid is found in hardcore boson Hubbard models with nearest and next-nearest neighbor repulsion Schmid1 ; Hebert .

Finally we examine the XY superfluid-VBS transition along cut (C). As illustrated in Fig. 4, simulation data for system size $`L=24`$ does not display an obvious sharp discontinuity as in Fig. 3. However, the presence of a small discontinuity in $`m`$ and $`\rho _s`$ for $`L=32`$ to 48 is suggested by the data. The inset of Fig. 4 displays a double-peaked magnetization probability histogram in the transition region, which indicates the presence of a first-order phase coexistence. This clearly precludes the existence of a continuous quantum phase transition, at least for the field value $`h/J=4`$ that was studied in cut (C) (see Fig. 1). The most immediate conclusion to draw is that the superfluid-VBS phase transition is weakly first order, either along its entirety (excluding the $`h=0`$ QCP), or up to a tricritical point at a field $`0<h<4`$. In this case, the difficulty in seeing a large discontinuity in the superfluid density or plaquette structure factor is due to the small $`h`$ value and the closeness of the magnetization to zero. The persistence of a small region of superfluid density in apparent coexistence with a finite VBS order parameter (for example, the two data points for $`L=48`$, $`K/J=11.60`$ and 11.65 in Fig. 4) is due to the first-order metastability between the superfluid phase and the VBS phase that is obscured by statistical averaging. As a check, we observed the Monte Carlo time correlation between $`S_p(\pi ,0)`$, $`S_p(0,\pi )`$, and the superfluid density $`\rho _s`$ in the $`x`$ and $`y`$-directions at these points. In fact, we find that $`\rho _s`$ show no preference for the $`x`$ or $`y`$ directions when stripe order is present. Rather, both $`\rho _s^y`$ and $`\rho _s^x`$ show very strong anti-correlations whenever $`M_P^2`$ develops Bragg peaks.

In summary, using SSE QMC techniques, we have determined the ground state phase diagram (Fig. 1) of the easy-plane $`J`$-$`K`$-$`h`$ model. In addition to the XY superfluid, VBS, and Néel ordered phases observed for $`h=0`$ AD1 , we observe a large region of fully polarized order, which dominates the phase diagram for large $`h`$. The phase transition to the polarized state is continuous at small $`K/J`$ and strongly first order for large $`K/J`$, suggesting the existence of a tricritical point somewhere on the phase boundary for intermediate $`K/J`$. Two other phase transitions were studied in detail, the superfluid-Néel and superfluid-VBS transitions. Both were first order for the parameter values investigated in detail here.

As indicated by our data, the $`J`$-$`K`$-$`h`$ model does not appear to support a region of superfluid-VBS coexistence (i.e. a supersolid), which is observed near a similar transition between a superfluid and $`(\pi ,0)`$ striped solid phase in a hardcore boson Hubbard model Hebert ; HubBosRG . No additional ordered phases were observed in this model, in particular, incommensurate VBS stripes (or striped order away from half filling) which would have been indicated by Bragg peaks in the q-dependent structure factor $`S_p(q_x,q_y)`$ away from $`(\pi ,0)`$.

In the context of the $`h=0`$ superfluid-VBS transition at $`T=0`$ AD1 ; AD2 , the existence of a continuous QCP does not require a continuous phase transition to develop smoothly as $`h`$ is increased from zero. Conversely, the existence of a true continuous phase transition would provide additional supporting evidence for the existence of the $`h=0`$ QCP RogNote , as well as a further region in which to explore the nature of the critical behavior associated with the transition from XY superfluid to VBS order. Ultimately, one would like to determine whether this QCP is an example of the “deconfined” quantum criticality recently discussed by Senthil et al. LeonMatt .

Finally, the inability of any significant region of a spin-liquid phase to develop in the vicinity of the QCP motivates further searches on related models. Of particular interest is the square lattice $`J`$-$`K`$ ring model in a staggered magnetic field, which could conceivable destabilize the superfluid or VBS order near the QCP and promote the development of an extended region of disorder. Work on this model is in progress RAD1 .

The authors would like to thank L. Balents and M. P. A. Fisher for insightful discussions and a critical reading of the manuscript. This work was supported by the Department of Energy, Grant No. DE-FG02-03ER46048 (DJS), and by The Academy of Finland, project No. 26175 (AWS). Supercomputer time was provided by NCSA under grant number DMR020029N. AWS would like to thank the Department of Physics at UCSB for hospitality and support during a visit.


text:
# 1 Conductivity vs. temperature for 𝑛_𝑠=2.46 (a) and 2.23×10¹¹ cm-2 (b) along with the slopes 𝑑𝜎/𝑑𝑇 recalculated from 𝑑𝜌/𝑑𝑇 of Ref. (). Also shown are the fits to the data using the traditional screening theory () with the effective mass values taken from Ref. ().

Comment on “Interaction Effects in Conductivity of Si Inversion Layers at Intermediate Temperatures”

In a recent Letter pudalov03 , Pudalov et al. have claimed that they have found an excellent agreement between their experimental data and theory zala01 . According to them, the anomalous (by an order of magnitude kravchenko03 ) increase of resistance with temperature can be quantitatively described by the theory of small corrections arising from quantum coherent interaction effects. Notably, this claim negates previous publication brunthaler01  entitled “Exclusion of Quantum Coherence as the Origin of the 2D Metallic State in High-Mobility Silicon Inversion Layers”, in which the same samples were studied in the same temperature and density ranges. The purpose of this Comment is to show that when analyzed correctly, the data of Pudalov et al. are not appropriate for accurate comparison with theory and do not allow one to distinguish between the interaction-based theory zala01  and the traditional screening theory gold86 .

Despite theory zala01  yields a linear-in-$`T`$ correction to conductivity ($`\sigma `$), Pudalov et al. choose to analyze the temperature-dependent resistivity ($`\rho `$). This serves to mask the fact that there is essentially no linear-in-$`T`$ interval on their $`\sigma (T)`$ dependences temp . To demonstrate this, in Fig. 1 we plot $`\sigma (T)`$ recalculated from the published $`\rho (T)`$ data for two lowest electron densities lowest  shown in Fig. 1 (b) of Ref. pudalov03 . The solid lines correspond to the slopes $`d\sigma /dT`$, which the authors claim to be in “excellent agreement” with their experimental data (circles). The obvious reason for the dramatic discrepancy between presentations of the data in the conductivity and resistivity forms are large changes of $`\rho `$ with temperature, reaching a factor of two.

Another point of importance is that the data of Ref. pudalov03  are not related to the “anomalous increase of $`\rho `$ with temperature” as they are taken at (relatively) high electron densities $`n_s>2\times 10^{11}`$ cm-2, where interaction effects are weak. Metallic corrections to conductivity in this regime have been observed dolgopolov88  and described by the traditional screening theory gold86  long time ago. The same theory provides reasonably good fits to the current data from Ref. pudalov03  (see, e.g., the dashed lines in Fig. 1). Note that at electron densities $`5\times 10^{11}`$ cm-2, theory zala01  can hardly be applied at all to Si inversion layers because surface roughness scattering becomes dominant.

Summarizing, the comparison performed by Pudalov et alpudalov03  between experimental data (including those for the magnetoresistance) and theory zala01  is not valid. Moreover, their data do not support the interaction-related origin of the metallic $`\rho (T)`$ at the relatively high electron densities used once both theories zala01 ; gold86  are treated on equal footing. As follows from accurate measurements in best samples shashkin ; kravchenko03 , the interactions become important at lower densities close to the metal-insulator transition, where they lead to a strong renormalization of the effective mass which corresponds to strongly temperature-dependent metallic resistivity.

AAS and VTD gratefully acknowledge financial support by the RFBR and the programs “Nanostructures” and “Mesoscopics” from the Russian Ministry of Sciences. SVK is supported by NSF grant DMR-0129652 and the Sloan Foundation.
A. A. Shashkin and V. T. Dolgopolov
Institute of Solid State Physics, Chernogolovka, Moscow District 142432, Russia
S. V. Kravchenko
Physics Department, Northeastern University, Boston, Massachusetts 02115, U.S.A.


text:
# Effects of Fluid Instabilities on Accretion Disk Spectra

## 1. Introduction

The standard assumption central to most accretion disk models is that the turbulence responsible for angular momentum transport and gravitational energy release can be parameterized by an “$`\alpha `$-viscosity” prescription (Shakura & Sunyaev 1973). The hypothesis that turbulence drives accretion has been borne out by the identification of the magnetorotational instability (MRI, Balbus & Hawley 1991) as a generic source for magnetohydrodynamic turbulence in differentially rotating flows. However, it is not yet known how accurately $`\alpha `$-disk models describe real accretion disks. It is possible that both the structure and emitted spectrum of a disk in which the MRI is active may deviate significantly from the predictions of ad hoc $`\alpha `$-disk models.

Numerical simulations of the radiation magnetohydrodynamic (RMHD) equations (Turner et al. 2003) show that MRI turbulence can produce a factor of 20 variation in density. In addition to the MRI, local linear analyses of the RMHD equations demonstrate that a radiatively driven photon bubble instability exists in radiation dominated accretion disks (Gammie 1998). This instability is present as long as the magnetic field is finite, but the growth rate for the instability is greatest when the magnetic pressure exceeds the gas pressure (Blaes & Socrates 2003). Numerical RMHD calculations of radiation pressure supported disks which neglect shear (Turner et al., in preparation) show that the non-linear development of the photon bubble instability does indeed produce large density variations. In section 2 we discuss the effects of these density inhomogeneities on the emergent spectrum.

Fluid motions associated with the MRI may also influence the emission. MRI turbulence is fundamentally magnetic in nature with a characteristic velocity $`v_w`$ given by the Alfvén speed $`v_A`$ which is less than the local sound speed $`c_s`$. Close to the hole, where the radiation pressure is dominant, $`c_s`$ is given by the radiation sound speed which is larger than the sound speed of the gas. Thus, it is possible for the MRI turbulence to be highly supersonic relative to the gas and Compton scattering off of the turbulent eddies may produce a larger spectral modification than thermal Comptonization (Socrates, Davis, & Blaes 2003, hereafter SDB03). In section 3, we summarize these results.

## 2. Monte Carlo Results

We calculate a full-disk-thickness numerical solution of the RMHD equations in a stratified, radiation supported disk using the ZEUS-2D code with flux limited diffusion. Rapid photon diffusion creates large density inhomogeneities with relatively small temperature gradients. We then perform 2-D MC and 1-D Feautrier radiative transfer calculations through the uppermost quarter of the simulation domain. We horizontally average the domain to create a 1-D density and temperature input profile for the Feautrier calculation. Both methods model the effects of electron scattering and free-free emission and absorption. The optical depth at the base of the MC calculation was chosen to be several times the unit effective optical depth.

In figure 1, we compare the emergent spectrum from the MC calculation with the Feautrier spectrum and blackbody emission at a representative domain surface temperature. Though both of the calculated spectra resemble modified blackbody curves, the 2-D calculation’s spectrum is more thermal (closer to blackbody) than the 1-D result. The stronger density dependence of absorption and emission relative to scattering is responsible for this difference. In the inhomogeneous 2-D MC domain, the surface flux is largest in the low density regions but is dominated by photons emitted in the densest regions where absorption is increasingly probable relative to scattering. Thus, on average, the photons escaping through the surface in the 2-D model will undergo fewer scatterings than those photons reaching the surface in the 1-D model resulting in a “less-modified” blackbody spectrum. This enhanced flux exceeds the local Eddington limit by a factor of order ten, consistent with predictions for highly inhomogeneous atmospheres (Begelman 2001).

## 3. Turbulent Comptonization

In order for turbulent Comptonization to be important in accretion disks, the bulk velocity of the turbulent eddies must exceed the thermal velocity of the electrons. To compare the turbulent and thermal velocities, it is convenient to define a turbulent wave temperature $`T_w`$ in terms of the turbulent velocity on the outer scale. The wave temperature can then be related to thin disk parameters (SDB03)

$$T_w\frac{m_e}{3K_B}v_w^29.0\times 10^9\frac{\alpha }{ϵ^2}\left(\frac{L}{L_{edd}}\right)^2r^3\mathrm{K}$$
(1)
Here, L is bolometric luminosity, $`\mathrm{L}_{\mathrm{edd}}`$ is the Eddington luminosity, $`ϵ`$ is the radiative efficiency, $`r`$ is the cylindrical radius in units of gravitational radii $`r_g=GM/c^2`$, $`\alpha `$ is the standard viscosity parameter, and the effects of general relativity and no-torque inner boundary conditions have been neglected. It can be seen from eq. (1) that the wave temperature dominates the thermal temperature near the inner radius of the disk for sufficiently high accretion rates. The photons can only sample the velocity field on scales of order the photon mean free path such that the relevant wave temperature is generally less than that given by eq. (1). Our lack of understanding of turbulent dissipation in accretion disks hinders our ability to model their vertical structure. This uncertainty is further compounded by the possibility that turbulent Comptonization may itself be a dominant means of energy release. Therefore, we adopt a simple homogeneous model to estimate the reduction in wave temperature appropriate for the reduced velocities.

1-D MC calculations show that for sufficiently high turbulent stresses, turbulent Comptonization may account for a significant amount of the far-UV and X-ray power of accreting sources. For example, the spectrum shown in figure 2 has significant power up to $``$1 keV and resembles soft X-ray excesses observed in many Seyfert 1 galaxies and QSO’s. If soft excesses in these sources are due to turbulent Comptonization, they provide a direct means for probing the turbulence.

## 4. Conclusions

Large density fluctuations due to the photon bubbles instability may alter the spectrum emitted in radiation pressure dominated accretion disks. These results are based on calculations which neglect shear. If shear is not neglected the MRI may affect the photon bubble instability, possibly altering the disk structure. However, for sufficiently large $`v_A`$, the MRI may be compressible in radiation dominated regions and may produce significant inhomogeneities on its own (Turner et al. 2003).

We also analyze the interactions of photons with turbulent eddies and show that for sufficiently large turbulent stresses, a prominent spectral component in the far UV and X-ray bands emerges. Observations of accreting systems may provide a means to directly probe the mechanism responsible for angular momentum transport in these sources.

## References

Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214

Begelman, M. C. 2001, ApJ, 551, 897

Blaes, O.,& Socrates, A. 2003, ApJ, in press (astro-ph/0304348)

Gammie, C. F. 1998, MNRAS, 297, 929

Shakura, N. I., & Sunyaev, R. A., 1973, A&A, 24, 337

Socrates, A., Davis S.W. & Blaes, O., 2003, ApJ, in press (astro-ph/0307158)

Turner, N.J., Stone,J.M., Krolik, J.H., & Sano, T. 2003, ApJ, 593, 922

Turner, N.J., Socrates, A., Blaes, O., & Begelman M.C., in preparation


text:
# Anisotropy in the Distribution of Satellite Galaxy Orbits

## 1. Introduction

Observations of the distribution of bright elliptical galaxies within the Virgo galaxy cluster show a remarkable collinear arrangement (West & Blakeslee 2000). Moreover, this axis also appears to be part of a filamentary bridge connecting Virgo and the rich cluster Abell 1367. This phenomenon has already been recognized (Arp 1968; Bingegeli, Tammann & Sandage 1987) but only with accurate measurements of distance did it become possible to discriminate between this being a genuine three-dimensional structure or merely a chance alignment of galaxies. West & Blakeslee (2000) based their distances upon the surface brightness fluctuations method (cf. Blakeslee, Ajhar & Tonry 1999 for a recent review), concluding that not only were the brightest cluster members distributed anisotropically over the cluster, but the distribution of dwarf ellipticals was elongated in the same direction (Bingelli 1999). Plionis et al. (2003) investigated this substructure-cluster correlation statistically using 303 Abell clusters. They affirm that there is indeed such a signal and also that this signal is related to the dynamical state (and the environment) of the cluster.

Moreover, even on smaller scales – in galactic systems – there are clear observational indications that the distribution of the orbits of satellite galaxies is biased towards the galactic pole (Holmberg 1969) giving an anisotropic distribution. Zaritsky et al. (1997) found that satellites of (isolated) disk galaxies are scattered asymmetrically about the parent galaxy and aligned preferentially with the disk minor axis. At first this signal was only prominent for distances out to $``$50 kpc from the host (Holmberg 1969), but an analysis based on a much larger sample (cf. Zaritsky et al. 1993) of satellites extends it to 200 kpc (Zaritsky et al. 1997). This result is confirmed by a study of the satellites orbiting M31 (Hartwick 2000; Grebel, Kolatt & Brandner 1999). In addition, for the one galaxy where individual satellite orbits are known, the Milky Way, there is also evidence that that the orbits are preferentially polar (Zaritsky & Gonzalez 1999). This is derived from information based upon the alignment of satellites on the sky (Kunkel & Demers 1976; Lynden-Bell 1982), the orientation of the Magellanic Stream (Mathewson, Clearly & Murray 1974), the three-dimensional distribution of satellites (Majewski 1994; Hartwick 1996), and their actual velocities (Scholz & Irwin 1994).

All these observations clearly indicate that on scales spanning from galaxy clusters down to galaxies there is a signal indicating a correlation between the alignment of substructure and the shape of the gravitational potential this substructure moves in. However, the source of this alignment of satellites has been debated for a number of years, with two potential solutions suggested to explain this puzzling arrangement. If the distribution of infalling satellites is initially spherical, then dynamical selection may preferentially destroy or suppress the star-formation in those on more equatorial orbits111With equatorial we mean ”perpendicular to the major axis”. For instance, Peñarubbia, Kroupa & Boily (2002) have argued that in flattened dark matter halos, the timescale for dynamical friction is substantially longer for satellites that are on polar orbits. Alternatively, the present non-isotropic distribution of satellite systems may reflect the fact that it has always been non-isotropic, with the satellites being accreted along preferred directions. These hypotheses can be tested in numerical simulations of structure formation. However, in previous attempts to decipher this signal from cosmological $`N`$-body simulations there still remains a certain amount of uncertainty (Zaritsky et al. 1997). The first fully self-consistent simulations targeting the subject were performed by Tormen (1997). They though lacked resolution in time, space and mass to perform a detailed analysis of the satellite population for differing environments, and hence more refined simulations are required. Tormen (1997) could not follow the satellite distribution within the host’s virial radius but rather tracked all progenitors prior to accretion. This allowed him to analyze the infall pattern rather than the orbital evolution of the satellites. The aim of this study is to investigate this subject with the latest state-of-the-art high-resolution $`N`$-body-simulations. We focus upon a detailed analysis of the temporal and spatial properties of satellite galaxies residing within host galaxy clusters that formed fully self-consistently within a cosmological framework.

## 2. The Simulations

There is mounting evidence that, despite its problems, the Cold Dark Matter structure formation scenario provides the most accurate description of our Universe. Observations point towards the standard $`\mathrm{\Lambda }`$CDM according to which the universe is comprised of about 28% dark matter, 68% dark energy, and luminous baryonic matter (i.e. galaxies, stars, gas, and dust) at a mere 4% level (cf. Spergel et al. 2003). This so-called ”concordance model” induces hierarchical structure formation in which small objects form first and subsequently merge to form larger objects.

The analysis presented in this study is based upon a suite of eight high-resolution $`N`$-body simulations within the $`\mathrm{\Lambda }`$CDM concordance cosmology. They were carried out using the publicly available adaptive mesh refinement code MLAPM (Knebe, Green & Binney 2001). MLAPM reaches high force resolution by refining all high-density regions with an automated refinement algorithm. The refinements are recursive: the refined regions can also be refined, each subsequent refinement having cells that are half the size of the cells in the previous level. This creates an hierarchy of refinement meshes of different resolutions covering regions of interest. The refinement is done cell-by-cell (individual cells can be refined or de-refined) and meshes are not constrained to have a rectangular (or any other) shape. The criterion for (de-)refining a cell is simply the number of particles within that cell and a detailed study of the appropriate choice for this number as well as more details about the particulars of the code can be found elsewhere (Knebe, Green & Binney, 2001).

Each run focuses on the formation and evolution of one particular dark matter halo containing of order greater one million particles. The force resolution reached by MLAPM in the high-density regions is 2$`h^1\mathrm{kpc}`$ corresponding to roughly 0.25% of the host’s virial radii $`R_{\mathrm{vir}}`$. Following Lacey & Cole (1994) we define the formation redshift $`z_{\mathrm{form}}`$ as the time when the halo contains half of its mass today. The host halos were carefully chosen to sample a variety of different triaxialities and merger histories as summarized in Table 1. The virial radii given in that Table correspond to the point where the density of the host (measured in terms of the cosmological background density $`\rho _b`$) drops below $`\mathrm{\Delta }_{\mathrm{vir}}=340`$ with $`M_{\mathrm{vir}}`$ being the mass enclosed by that sphere. All satellite galaxies orbiting in and around the host halo are identified using the MLAPM-based halo finder described in a companion paper (Gill, Knebe & Gibson 2004). This new halo finder is based upon the MLAPM grid structure and works at the same resolution level as the $`N`$-body code itself thus assuring not to miss or include any objects not resolved in the actual simulation. The satellites are then individually traced from $`z_{\mathrm{form}}`$ onwards until redshift $`z=0`$. The mass spectrum of those satellites accounted for in the following analysis can be described by a declining power-law $`dn/dMM^\alpha `$ with $`\alpha 1.71.9`$ in the range from $`2\times 10^{10}`$$`h^1\mathrm{M}_{}`$ (applied mass-cut corresponding to 100 particles and explaining the rather ’low’ number for $`N_{\mathrm{sat}}(<r_{\mathrm{vir}})`$ in Table 1) up to $`10^{13}`$$`h^1\mathrm{M}_{}`$. The total fraction of mass locked up in satellites never exceeds $``$15% though.

The high temporal sampling of the outputs (i.e. $`\mathrm{\Delta }t0.2`$ Gyrs) allows us to accurately measure the orbital parameters of the satellite galaxies and a detailed analysis of their disruption and survival history will be presented elsewhere (Gill, Knebe & Gibson 2004). Gill et al. show that the distribution of orbital eccentricities for the satellite population peaks about the value $`ϵ0.6`$ where $`ϵ=1p/a`$ is defined using the last pericenter $`p`$ and apocenter $`a`$ of its orbit. The pericenter distribution itself peaks about 35% of $`R_{\mathrm{vir}}`$ for all host halos and most of the satellites (i.e. more than 75%) had had one full orbit with the maximum number of orbits found to be 3–5 depending on the host halo. The angular momentum vector of the satellite $`\stackrel{}{L}_{\mathrm{sat}}^{\mathrm{apo}}`$=$`\stackrel{}{R}_{\mathrm{sat}}^{\mathrm{apo}}`$$`\times \stackrel{}{V}_{\mathrm{sat}}^{\mathrm{apo}}`$ was calculated at its last apocenter passage $`\stackrel{}{R}_{\mathrm{sat}}^{\mathrm{apo}}`$.

In order to probe the alignment with the shape and orientation of the host halo we calculated the eigenvectors $`\stackrel{}{E}_{1,2,3}`$ (with $`\stackrel{}{E}_1`$ being the major axis) of its inertia tensor using only the ”core” region as defined by the $`6^{\mathrm{th}}`$ refinement level in MLAPM, i.e. the boundary of this refinement level is an isodensity contour. According to the refinement criterion adopted in the simulations the $`6^{\mathrm{th}}`$ level surrounds material about 3000 times denser than $`\rho _b`$ or in other words 9 times denser than the material at the virial radius. The host halos are well described by the density profiles advocated by Navarro, Frenk & White (1997) with concentration in the range $`c=57`$. Therefore, a density of roughly $`9\times \rho (r_{\mathrm{vir}})`$ corresponds to about the half-mass radius of the host. The eigenvectors now define the orientation of the host halo and the coordinate system used to measure the orbits of the satellites, respectively. Moreover, its eigenvalues $`a>b>c`$ can be used to construct the triaxiality parameter $`T=(a^2b^2)/(a^2c^2)`$ (Franx, Illingworth & Zeeuw 1991) presented in Table 1, too.

The left panel of Fig. 1 presents the (normalized) distribution of angles between $`\stackrel{}{R}_{\mathrm{sat}}^{\mathrm{apo}}`$ and $`\stackrel{}{E}_{1,\mathrm{host}}`$. This graph shows that there is a clear trend in at least six of the eight halos for the two vectors to be aligned (the distribution peaks at $`0^{}`$ and $`180^{}`$, respectively) meaning that the orbits of the satellites are preferentially found along the major axis of the host. The right panel of Fig. 1 shows the (normalized) distribution of angles between the angular momentum vector $`\stackrel{}{L}_{\mathrm{sat}}^{\mathrm{apo}}`$ of the satellite at its last apocenter and $`\stackrel{}{E}_{1,\mathrm{host}}`$, clearly revealing evidence for these two vectors to be parallel (distribution peaks at $`90^{}`$). We like to stress that only satellite galaxies that at least had one or more complete orbits were taken into account in Fig. 1; the figure is not based upon the infall pattern of satellites as investigated by, for instance, Tormen (1997). We also need to stress that the distributions presented in that Figure are normalized, i.e. they are corrected for the bias introduced by plotting them as a function of the angle $`\theta `$ rather than $`\mathrm{cos}\theta `$. The data is binned equally spaced in angle ranging from 0 to 180, and therefore the area probed on the sphere varies with $`\theta `$. Therefore the distribution needed to be normalized by the the respective area $`A`$ specified by the actual range of angles \[$`\theta \mathrm{\Delta }\theta /2,\theta +\mathrm{\Delta }\theta /2`$\]. This area is proportional to $`A\mathrm{cos}(\theta \mathrm{\Delta }\theta /2)\mathrm{cos}(\theta +\mathrm{\Delta }\theta /2)`$.

## 3. Conclusions and Discussion

There are two possible scenarios that explain the substructure-cluster alignment found in observations of galaxy clusters (Plionis et al. 2003; West & Blakeslee 2000) and the analysis of $`N`$-body-simulations presented in this study. Either some dynamical process can be held responsible, meaning that initially the orbits were randomly distributed and only some satellites survive to give rise to the observed correlation in Fig. 1, or alternatively, their orbital parameters are imprinted upon them at the time they enter the host halo as already pointed out by Tormen (1997).

“Dynamical destruction”, whereby satellites on more equatorial orbits suffer more tidal disruption, can be rejected because a comparable (although not as pronounced) alignment signal is found when restricting the analysis to those satellites that are in fact disrupted: We placed mass-less tracer particles at the last centre of the satellites before they were classified ”disrupted”. These “disrupted systems” are presented as the thin lines in Fig. 1.

We therefore conclude that it must be the initial distribution of the satellite systems that is responsible for the present day alignment. For example, galaxies are ”funneling along the filaments” (Kitzbichler & Saurer 2003) which has also been confirmed by X-ray observations of the spatial distribution of substructure in galaxy clusters (West, Jones & Forman 1995) and biases the types of possible orbits. It would now be reassuring to confirm a link between the host halo shape and the surrounding environment, i.e. are the filaments that feed the halo with material preferentially angled with respect to the orientation of the host halo? When projecting the positions of the satellites onto a sphere at the time they enter the virial radius of the host it is clear that they are not randomly distributed; they cluster in directions linked to the filamentary structure surrounding the host halo, as found before observationally (e.g. West & Blakeslee 2003; Plionis & Basilikos 2002, and references therein) and in $`N`$-body-simulations (e.g. Faltenbacher 2002; Hatton & Ninin 2001; Onuora & Thomas 2000; Colberg et al. 1999; Splinter et al. 1997; Tormen 1997). Tormen (1997) already pointed out that there is a strong alignment between the distribution of infalling satellite galaxies and the shape of the dark matter host. However, Tormen’s simulations did not have the spatial and mass resolution to investigate the orbits of satellites within the host halo. We are, for the first time, tracing extremely well-resolved dynamics of the satellites within the dark matter hosts confirming that the alignment for both, disrupted and surviving, satellites is maintained for several orbital periods. Fig. 2 now presents a visual impression of the trajectories of all surviving satellites identified at formation time of host halo #3 until redshift $`z=0`$, qualitatively supporting the scenario that these orbits are linked to the filamentary structure which will leave its imprint in observed anisotropy.

Plionis et al. (2003) also established a link between the alignment and the dynamical age of the clusters. Our simulations though show that this finding can not be generalized: both our oldest halos, in which the satellites had as many as 4–5 orbits, do show the correlation signal. We do not observe a trend for a ”randomization” of the orbits in older halos. The satellites actually preserve the alignment with the host they had when they first fell into the cluster.

In Section 1 we showed that the (observational) substructure-alignment signal spans from galaxy clusters down to galactic scales. Our results are mostly applicable to the former, and any extrapolation to smaller scales has to be handled with care. In the case for galaxies, for instance, the interaction between the galactic disk and the incoming satellite could enhance destruction in the galaxy plane, i.e. satellite on prograde orbits decay faster than the ones on retrograde (or polar) orbits due to orbital resonances between the disk and the satellites (Peñarubbia, Kroupa & Boily 2002). Based on theoretical and numerical studies (cf. Lacey & Cole 1993 and Moore et al. 1999) it should be possible to re-scale our simulations to a Milky Way sized object by requiring that the maximum of the circular velocity curves of our halos equals 220 $`\mathrm{km}\mathrm{s}^1`$. The scaling factor lies in the range of about 4–5 (cf. $`v_{\mathrm{circ}}^{\mathrm{max}}`$ values in Table 1) and hence our “rescaled dark matter halos” would correspond to the ones in the observational data with virial radii in the range 260–295$`h^1\mathrm{kpc}`$. This scaling factor in length entails a scaling in mass of 64–125 (simply the length scale to the power of three). This brings the mass of our galaxy clusters down to $`10^{12}`$$`h^1\mathrm{M}_{}`$ which agrees with the dark matter mass inferred for our Milky Way (Freeman 1996). More problematic, however, are the ages of our systems: our halos are only $``$8.3 Gyrs old opposed to 12 Gyrs for the Milky Way. Satellites in the Milky Way had the chance to complete nearly twice as many orbits leaving more space for an explanation based upon the dynamical destruction scenario.

When re-scaling our data, putting aside the age issue, and trying to explain the Holmberg effect, another uncertainty comes into play: the orientation of the stellar disk. In a triaxial potential there are in general stable closed orbits about both the major and the minor axis (e.g., Binney & Tremaine 1987). So, in principal the disk plane could either be perpendicular to the major or perpendicular to the minor axis. Under the assumption it lies perpendicular to the major axis of the dark matter halo the satellites in our simulations will be on polar orbits. There are indications that this configuration results in the most stable disk configuration within a triaxial halo (Hayashi et al. 2003). Moreover, even though there are clear indications that the angular momentum of the dark matter is well aligned with the minor axis of the halo (e.g. Warren et al. 1992), van den Bosch et al. (2002) showed that the angular momentum of the baryonic component (i.e. gas) not necessarily follows that of the dark matter distribution. They found an average misalignment between $`\stackrel{}{L}_{\mathrm{gas}}`$ and $`\stackrel{}{L}_{\mathrm{DM}}`$ of the order of 40 in their numerical simulations. Hence the orientation of the galactic disk with respects to the dark matter halo is not well determined. Turning to observations does not resolve this question, with studies of polar rings indicating strongly oblate halos (i.e. Iodice et al. 2003), while others suggest dark matter halos are more spherical, or even oblate (Olling & Merrifield 2000; Ibata et al. 2001). The situation is not at all clear, but the assumption that disks are perpendicular to the major axis (which is in agreement with results presented by Hayashi et al. 2003) would provide an explanation for the Holmberg effect, even though this is a very speculative interpretation.

The simulations presented in this paper were carried out on the Beowulf cluster at the Centre for Astrophysics & Supercomputing, Swinburne University. The support of the Australian Research Council and the Swinburne Research Development Grants Scheme is gratefully acknowledged. GFL thanks Suede for “Introducing the band”.