text:
# Stochastic Cellular Automata Model for Stock Market Dynamics

## I Introduction

Since the successful application of the Black-Scholes theory for option pricing black73  in 1973 more and more physicists have been attracted by the idea of understanding the behaviour of the market dynamics in terms of complex system theory, where self-organized criticality soc ; bak1997  and stochastic processes mantegna ; paul  play important roles. The aim of the microscopic models proposed so far (for general reviews levy2000 ; feigenbaum2003 ) is to reproduce some stylized facts sornette99  concerning the temporal fluctuations of the price indices, $`P(t)`$. In particular, the logarithmic price returns

$$R(t)=lnP(t+1)lnP(t),$$
(1)
and the volatility, defined in the present work as

$$v(t)=|R(t)|,$$
(2)
have been studied extensively mantegna  from an empirical point of view. The results have shown that while long time correlations are present in the volatility, a phenomenon known as volatility clustering, they cannot be found in the time series of returns. Moreover the latter show an intermittent behaviour that recalls in some aspects hydrodynamic turbulence mantegna97 ; mantegna95 ; frisch , characterized by power law tails in the probability distribution function (pdf). Microsimulations have demonstrated that this kind of behaviour can originate both as a stochastic process with multiplicative noise kaizoji ; krawiecki02 ; takayasu9798  and as a percolation phenomenon stauffer ; stauffer98 ; cont00 .

In order to reproduce these features of real markets we introduce a stochastic cellular automata model, representing an open market. That is, a market where the number of active traders, defined as cells with spin state different from 0, namely $`\sigma _i(t)=\pm 1`$, evolve in time according to a percolation dynamics. The percolation dynamics is chosen in order to simulate the herding behaviour typical of investors sornette03 . According to this, active traders gather in clusters or networks where, following a stochastic exchange of information, they formulate the trading strategy for the next time step. The results obtained by the simulations are then compared with the time series of daily closures of the S&P500 index mantegna ; plerou99 ; cizeau97 ; gopikrishnan99 ; liu99 ; sornette03 ; ausloos2003 ; Tsallis2003 ; canessa2000  over a period of about 50 years.

Moreover, recently, the fractal properties feder  of the price fluctuations have also been investigated for different markets rodrigues01 ; gorski02 ; auloos02 ; dimatteo03 . A common feature found in these studies is the existence of a nonlinear, multifractal spectrum that excludes the possibility of efficient market behaviour bachelier00 . The origin of the multifractality of in the financial time series has also been at the center of discussions lux2001 ; mandelbrot2001 ; lebaron2001 ; stanley2001 . In this paper we consider the multifractal spectrum of the price fluctuations as a stylized fact of the market time series without addressing any question about the underlying process able to generate it. The multifractal spectrum is used as a further test for our model.

A parallel between multifractal and thermodynamical formalism has also been investigated. We found, in agreement with the previous work of Canessa canessa2000 , that the analogue specific heat can provide a good tool to characterize intermittency, that is financial crashes or bubbles, from a thermodynamics-equivalent point of view.

## II The Model

In the present work we simulate the financial market dynamics via a stochastic cellular automata model. The agents of the market are represented by cells on a two dimensional grid, 512x128. The $`i`$th agent at the discrete time step $`t`$ is characterized by three possible states or spin orientations, $`\sigma _i(t)=0,\pm 1`$. The value $`\sigma _i(t)=+1`$ is associated with the purchase of a stock while $`\sigma _i(t)=1`$ with selling. The former states are called active. The cells with spin value $`\sigma _i(t)=0`$ are inactive traders. The active traders herd in networks or clusters via a direct percolation method related to a forest fire model stauffer . The information carried by the active traders, that is their spin state, is shared with the other members of the cluster. The percolation dynamics allow a time dependent herding behaviour and the market can be interpreted as an open system not bounded by conservation laws. The clustering process will be discussed in detail in the next subsection.

The trading dynamics is instead related to the synchronous update of the spins of the active traders, ruled by a stochastic exchange of information between them, similar to a random Ising model kaizoji ; krawiecki02 . A particular feature of the present simulation is that the information is not spread all over the grid, as in other multi-agent scenarios kaizoji ; krawiecki02 , but it is limited by the clusters of interaction previously defined. The mechanism for the spin dynamics is explained in “stochastic trading dynamics” subsection.

### II.1 Percolation Clustering

One of the aims of our cellular automata model is to reproduce the herding behaviour of active traders sornette03 . We refer to herding behaviour as the tendency of people involved in the market to aggregate in networks or clusters of influence. The traders then use the information obtained by their network in order to formulate a market strategy. Even if the topological structure of these networks of information is not important, since several kinds of long range interaction are available nowadays krawiecki02 , the number of connections for each trader must be, in any case, finite and not extended over the whole market. In this framework a direct percolation method is used to simulate herding dynamics between active traders. If we assume that the neighbours of influence are those of von Neumann (up, down , left, right), the percolation is fixed by the following parameters:

$`p_h`$: the probability that an active trader can turn one of his inactive neighbours into an active one at the next time step, $`\sigma _i(t)=0\sigma _i(t+1)=\pm 1`$. This simulates the fact that certain information possessed by a trader may induce a potential trader to join the market dynamics.

$`p_d`$: the probability that an active trader diffuses and so becomes inactive, $`\sigma _i(t)=\pm 1\sigma _i(t+1)=0`$, due to each of his inactive neighbour. This mimics the fact that only traders at the borders of a network, that is the weaker links, can quit the market.

$`p_e`$: the probability that a non trading cell spontaneously decides to enter the market dynamics, $`\sigma _i(t)=0\sigma _i(t+1)=\pm 1`$.

The values of the adimensional parameters, $`p_h`$, $`p_d`$ and $`p_e`$, influence the stability of the system and the percentage of active traders on the grid. In order to test different market activities we fix the values $`p_d=0.05`$ and $`p_e=0.0001`$ while we tune the parameter $`p_h`$. At the beginning of the simulation the grid is loaded randomly with a small percentage of active traders and then the system is permitted to evolve according to the previous rules. If we are in a stable range of the parameter $`p_h`$, after a transient period that depends both on the parameter values and the initial number of active cells, the number of active traders on the grid begins to fluctuate around a certain average, as shown in Fig. 1 (Top). The market can be considered open since the number of agents changes dynamically in time. In this regime, the competition between herding and diffusion produces a power law distribution of the cluster size, as shown in Fig. 1 (Bottom), $`\rho (S)S^\lambda `$, where $`S`$ is the cluster dimension, defined as the number of active cells belonging to the same cluster, and $`\lambda >0`$, creating a hierarchy of networks. This hierarchy is necessary if were are to take into account a real aspect of the market, namely that different traders also have different trading powers. A reasonable assumption is that people having a larger number of sources of information, so belonging to greater clusters, can be associated with professional investors that, most likely, are able to move a greater amount of stocks compared to the occasional investor. Using this assumption we are able to define a proper weight for the trading power of different cells, as we will discuss in the next subsection.

A similar percolation model has also been used to reproduce some statistical and geometrical features of solar activity wentzel ; vlahos .

### II.2 Stochastic Trading Dynamics

The dynamics of the spins of the active traders, $`\sigma _i^k(t)=\pm 1`$ for $`i=1,\mathrm{},N^k(t)`$ (where the superscript $`k`$, from now on, refers to the $`k`$th cluster of the grid configuration at time step $`t`$) follows a stochastic process that mimics the human uncertainty in decision making krawiecki02 . Their values are updated synchronously according to a local probabilistic rule: $`\sigma _i^k(t+1)=+1`$ with probability $`p_i^k`$ and $`\sigma _i^k(t+1)=1`$ with probability $`1p_i^k`$. The probability $`p_i^k`$ is determined, by analogy to heat bath dynamics with formal temperature $`k_bT=1`$, by

$$p_i^k(t)=\frac{1}{1+e^{2I_i^k(t)}},$$
(3)
where the local field, $`I_i^k(t)`$, is

$$I_i^k(t)=\frac{1}{N^k(t)}\underset{j=1}{\overset{N^k(t)}{}}A_{ij}^k\sigma _j^k(t)+h_i^k.$$
(4)
The $`A_{ij}^k(t)`$ are time dependent interaction strengths between agents and $`h_i^k(t)`$ is an external field reflecting the effect of the environment krawiecki02 . The interaction strengths and the external field change randomly in time according to $`A_{ij}^k(t)=A\xi ^k(t)+a\eta _{ij}(t)`$ and $`h_i^k(t)=h\zeta _i^k(t)`$. The variables $`\xi ^k(t),\eta _{ij}(t),\zeta _i^k(t)`$ are random variables uniformly distributed in the interval (-1,1) with no correlation in time or space. The measure of the strengths of the previous terms, $`A`$, $`a`$ and $`h`$, are constant and common for all the grid.

In this contest the dynamics of the price index, $`P(t)`$, can be easily derived if we assume that the index variation is proportional to the difference between demand and supply,

$$\frac{dP}{dt}xP,$$
(5)
and using a weighted average for the orientation of the spins,

$$x(t)=\beta \underset{k=1}{\overset{N_{cl}(t)}{}}\underset{i=1}{\overset{N^k(t)}{}}N^k(t)\sigma _i^k(t),$$
(6)
where $`N_{cl}(t)`$ and $`N^k(t)`$ are, respectively, the number of clusters on the grid and the size of the $`k`$th cluster, while $`\beta `$ is a normalization constant. The relation (6) follows from the assumption of proportionality between the financial power of an active cell and the size of the cluster to which it belongs. A justification of this relation is that professional investors, able to move larger amounts of capital, are more likely to be linked with a large number of other investors than the occasional traders. From (5) we find that the logarithm of the price returns (1), that is the fundamental quantity we aim to model, is proportional to the average orientation defined in (6), $`R(t)x(t)`$.

## III Numerical Results

We compare the results of our model with the Standard & Poor 500 (S&P500) index that is one of the most widely used benchmarks for US equity performance. The database analyzed is composed by of the daily indices from 3/1/1950 to 18/7/2003 for a total of 13468 data. The time series of the index prices, $`P(t)`$, is converted into the logarithmic returns (1) and then normalized over the time interval, $`T`$,

$$r(t)=\frac{R(t)R(t)_T}{\sigma (R(t))}.$$
(7)
where $`\mathrm{}_T`$ denotes the temporal average and $`\sigma `$ is the standard deviation. In this way all the sets of returns will have zero average and $`\sigma =1`$.

Before discussing the results we briefly describe the behaviour of the cellular automata with respect to changes of the parameters. Regarding the percolation, the herding parameter $`p_h`$, as previously seen, determines the concentration of active traders on the grid. For a low concentration the active traders will be distributed in small clusters and the information will be extremely split up on the grid. In this case large coherent events will be more rare. A higher concentration, where clusters of the order of a thousand cells are present, allows large events to occur with higher frequency. In fact, the phase transition of large clusters can easily trigger a crash or a bubble in the market. Time series related to a higher herding parameter are therefore more intermittent. We can then infer that a market made by small groups of traders behaves like a noisy market, while big crashes or bubbles must necessarily be related to the interaction of large networks and so to a kind of crowd behaviour. For most of the simulations we fixed the parameter $`p_h=0.0493`$ in order to have bigger clusters, with $`S^{max}(t)`$ of the order of 2500 cells.

The strength $`A`$ also plays an important role in the trading dynamics. With $`p_h`$ fixed, this parameter is related to the intermittency of the system. Both for large values of the activity ($`A>10`$) and for $`A0`$ we observe an approach of the pdf toward a Gaussian-like shape. That is very, large fluctuations become more and more rare, and $`A`$ can be regarded as a temporal scale for the system, similar to the activity parameter in the Cont-Bouchaud model cont00 . In spite of this some large fluctuations can be still identified. This is probably one of the main differences between the Cont-Bouchaud model and the present. In fact Monte Carlo simulations of the former stauffer98  show that an increase of the activity brings a rapid convergence toward a Gaussian distribution because a large number of clusters are trading at the same time following a random procedure of decision making cont00 ; stauffer98 : there are no clusters that can influence the market more than others and so the resulting global interaction is noise-like. In our model fluctuations are always allowed because of the heat bath dynamics. Clusters of active traders can always be subjected to phase transitions, independently of the state of other clusters, creating a displacement between demand and supply.

In order to reproduce the behaviour of financial time series we work with $`1.5<A<2.5`$. In this range there is a nonlinear dependence of the intermittency on $`A`$. The value of the parameter $`a`$ is fixed by the relation $`a=2A`$, already used in the work of Krawiecki et al. krawiecki02 . However, we note that the ratio $`a/A`$ is actually not essential in reproducing the intermittency found in real data. Even the parameter $`h`$ does not have a central role in the simulation, as long as $`h<<1`$. In fact, the percolation dynamics of the active traders already introduces a natural noise. In most of the runs we simply set $`h=0`$.

Now we discuss the results of the cellular automata. In Fig. 2 (Top) and Fig. 2 (Bottom) the normalized logarithmic returns of the S&P500 and of the simulation are shown, respectively. The average number of active traders, in the stable regime, is $`16000`$, as shown in Fig. 1 (Top). The model reproduces the intermittent behaviour of the S&P500 time series, as expressed by the leptokurtic pdf of Fig.3. The tails of the distribution follow a power law decay, reflecting the fact that large coherent events, far from the average, are likely to occur with a frequency higher than expected for a random process (where the shape would be a Gaussian). These large events are related to financial crashes or bubbles of the market and, in our model, to a phase transition in the spin state of large networks of active traders, as we will discuss further on. From a power law fit, $`\rho (r)r^{1\gamma }`$ (for $`|r|>2`$), we find $`\gamma 3`$ for both the S&P500 and the model, confirming the good agreement between the two.

The problem of finding the best distribution describing the price returns is a very important issue from a practical point of view mantegna . The standard Black-Scholes theory for option pricing black73 ; mantegna ; paul  assumes that the returns are normally distributed. This fact has been proven to be empirically false, as shown also in Fig. 3 (see Ref. mantegna  for a general reference). Finding a more appropriate distribution would be an important improvement in this field of research***An intriguing framework of investigation has been provided by the non-extensive statistical mechanics proposed by Tsallis Tsallis1988 ; Tsallis1999 ; tsallis\_market ; ausloos2003 ; Tsallis2003 . A more complete discussion on this important topic is beyond the scope of this paper and will be discussed in future work..

In order to understand the trading dynamics of the automata we also show two snapshots of the grid configuration, the first during a normal session, Fig.4 (Top), and the second during a crash, Fig.4 (Bottom). During the normal session the orientations of the spins are distributed uniformly over the various clusters and there is no sharp difference between demand and supply. The situation is different during a crash. In this case the clusters at the top of the hierarchy, the bigger ones, play a fundamental role. In fact they undergo a phase transition where the greatest part of their spins share the same orientation. The capacity of the clusters to generate a coherent orientation of the spins, and hence of their trading state, can be interpreted in terms of a multiplicative noise process paul ; liu99 , where the collective synchronization arises as a result of the randomly varying interaction strengths between agents. The peculiarity of our model is that crashes or bubbles (sudden price changes) are related not to a phase transition of the whole market krawiecki02 ; kaizoji  but rather to a phase transition in one or more of the larger clusters that have a greater influence on the trading session. This behaviour is probably closer to the real market where the synchronization of trading opinion is more likely to happen between large groups of traders than over the whole market.

The temporal correlations of the logarithmic returns and of the volatility are investigated via the autocorrelation function, defined as

$$c(\tau )=\underset{t=1}{\overset{T\tau }{}}x(t+\tau )x(t),$$
(8)
where T is the length of the time series and $`\tau `$ is a time delay for the normalized variable $`x(t)`$. The results for both the model and the S&P500 are shown in Fig.5 (Top) and Fig.5 (Bottom), respectively. While the temporal correlation for the returns is lost almost immediately, the volatility manifests a slow decay in time, related to the phenomenon of volatility clustering. The previous temporal dependencies have been found in both real data and in the simulation.

## IV Multifractal Analysis

It is also worth pointing out that financial time series present an inherent multifractality feder . In the past few years the work of many authors rodrigues01 ; gorski02 ; auloos02 ; dimatteo03  has been addressed to the characterization of the multifractal properties of financial time series, and nowadays multifractality can be considered as a stylized fact. In order to study the multifractal properties of our model we use the generalized Hurst exponent mandelbrot , $`H(q)`$, derived via the $`q`$order structure function,

$$S_q(\tau )=|x(t+\tau )x(t)|^q_T\tau ^{qH(q)},$$
(9)
where $`x(t)`$ is a stochastic variable over a time interval $`T`$ and $`\tau `$ the time delay. The generalized Hurst exponent, defined in (9), is an extension of the Hurst exponent, $`H`$, introduced in the context of reservoir control on the Nile river dam project, around 1907 feder ; hurst51 . This technique provides a sensitive method for revealing long-term correlations in random processes. If $`H(q)=H`$ for every $`q`$ the process is said to be monofractal and $`H`$ is equivalent to the original definition of the Hurst exponent. This is the case of simple Brownian motion or fractional Brownian motion.

If the spectrum of $`H(q)`$ is not constant with $`q`$ the process is said to be multifractal. From the definition (9) it is easy to see that the function $`H(1)`$ is related to the scaling properties of the volatility. By analogy with the classical Hurst analysis, a phenomenon is said to be persistent if $`H(1)>1/2`$ and antipersistent if $`H(1)<1/2`$. For uncorrelated increments, as in Brownian motion, $`H(1)=1/2`$. In Fig. 6 a comparison is shown between the multifractal spectra of the model and the S&P500 obtained from the prices time series. It is clear that both processes have a multifractal structure and the price fluctuations cannot be associated with a simple random walk as in the classical efficient market hypothesis bachelier00 .

The multifractality of the time series can also be discussed in terms of thermodynamic equivalents, according to multifractal physics lee88 ; canessa1993 ; vollhardt97 ; canessa2000 . In this approach we divide the time series, $`x(t)`$ , for $`t=1,\mathrm{},L`$ into $`𝒩`$ equal sub-intervals. Then we can write the following measure for each of these,

$$\mu _i(\tau )=\frac{|x(t+\tau )x(t)|}{_{t=1}^𝒩|x(t+\tau )x(t)|},$$
(10)
with $`i=1,\mathrm{},𝒩`$ and $`\tau `$ the time delay. The quantity $`\mu _i(\tau )`$ can be viewed as a normalized probability measure. The corresponding generating function is given by

$$Z(q,𝒩)=\underset{i=1}{\overset{𝒩}{}}\mu _i(\tau )^q𝒩^{\chi _q},$$
(11)
that is an analogue of the partition function in thermodynamics. According to Ref. lee88  the scaling exponent $`\chi _q`$ is directly related to the generalized multifractal dimension $`D_q`$ feder ,

$$\chi _q(q1)D_q.$$
(12)
If we consider $`\chi _q`$ as the free energy of our system, Eq.(12) provides a link between the classical thermodynamical formalism and multifractality. Assuming $`q`$ as an equivalent temperature, we can define an analogue specific heat lee88 ; canessa1993 ; canessa2000 ; ivanova02 ,

$$C_q=\frac{^2\chi _q}{q^2}.$$
(13)
The previous equation for $`C_q`$ can also be written in terms of singular measure formalism davis94 . In this framework we define a measure, $`ϵ(1,t)`$, as

$$ϵ(1,t)=\frac{|x(t+1)x(t)|}{\frac{1}{L1}_{t=1}^{L1}|x(t+1)x(t)|}.$$
(14)
We can then generate a series of measures on shorter intervals of length $`\delta `$, $`ϵ(\delta ,l)`$ where $`\delta `$ is an integer power of two and $`l`$ is the index of the first point of the sub segments at that resolution. The average measure in the interval \[$`l,l+\delta `$\] is

$$ϵ(\delta ,l)=\frac{1}{\delta }\underset{l^{}=l}{\overset{l+\delta 1}{}}ϵ(1,l^{}),$$
(15)
for $`l=0,\mathrm{},L\delta `$. In this case we have a scaling property for the ensemble average with respect to the scale $`\delta `$:

$$ϵ(\delta ,l)\delta ^{K_q}.$$
(16)
In a multifractal process the exponent $`K_q`$ is a nonlinear function of $`q`$ related to the intermittency of the time series and to the generalized dimension via hentschel83 ; halsey86

$$(q1)D_q=q1K_q.$$
(17)
From (13) and (17) we have that ivanova02 ,

$$C_q=\frac{^2K_q}{q^2}.$$
(18)
Following Ref. canessa2000  we have found the analogue specific heat for both the S&P500 and for the price time series generated by our model, see Fig.7. For $`\tau =1`$ we observe a double-humped shape for both the model and the empirical data. If we take longer time delays the shoulder on the right hand side disappears, leaving only a sharp peak, similar to a first order phase transition around $`q=1.5`$. The results are in agreement with the analysis of Canessa canessa2000  and recall the Hubbard model for small to intermediate values of the local interaction vollhardt97 . As also suggested in canessa2000 , the second peak is due to the large fluctuations at small scales, that is crashes and bubbles. Increasing the time delay means that the fluctuations tend to be smoothed and the time series of returns approach a noise-like regime. For this reason the analogue specific heat shapes for $`\tau =100`$ are basically indistinguishable. From this argument we can interpret the analogue specific heat, and in particular the second peak for low time delays, as a way to characterize crashes, or in general the degree of intermittency in a time series. The difference in the shapes and heights of the shorter peak for $`\tau =1`$ is due to a slightly different correlation of the fluctuations in the two time series. Moreover, in the model we can link the shorter peak to the physical phase transition in the spin state of a network of traders.

## V Conclusions

In this paper we have introduced a stochastic cellular automata model for the dynamics of the financial markets. The main difference between our model and other stochastic simulations based on spin orientation of agents krawiecki02 ; kaizoji  is the temporal evolution of the networks of interaction and therefore the concept of an “open” market. The active traders follow a direct percolation dynamics in order to aggregate in networks of information. This makes our simulation, even if still a raw approximation, surely closer to the real market, where no conservation rules for the number of agents can be claimed. Crashes and bubbles can be interpreted as a synchronization of the spin orientation of the more influential networks in the market. Moreover, the introduction of a limitation in the number of interacting agents reduces drastically the number of computations on the grid per time step. In a system where all the agents interact with each other this number goes like $`N_a^2`$, being $`N_a`$ the number of active agents, while in our model, considering the distribution of the clusters, it is easy to see that it goes like $`N_a^{2\lambda }`$. The value of $`\lambda `$ found for several herding parameters, $`p_h`$, is $`\lambda 0.6÷1.4`$, so that the computational cost is much lower for the present model. This gives one the possibility to simulate the market using a very large range of agents. The model is able to reproduce most of the stylized aspects of the financial time series, supporting the idea that crashes and bubbles are related to a collective synchronization in the trading behaviour of large networks of traders, where the information is exchanged according stochastic interaction between them.

###### Acknowledgements.

This work was supported by the Australian Research Council.


format: text

text:
# 1 The COMPETE analysis of forward data

## 1 The COMPETE analysis of forward data

Analytic parametrizations of forward $`(t=0)`$ hadron scattering amplitudes is a well-established domain in strong interactions.

However, in the past, the phenomenology of forward scattering had quite a high degree of arbitrariness : i) An excessive focus on $`pp`$ and $`\overline{p}p`$ scattering ; ii) Important physical constraints are often mixed with less general or even ad-hoc properties ; iii) The cut-off in energy, defining the region of applicability of the high-energy models, differs from one author to the other ; iv) The set of data considered by different authors is sometimes not the same ; v) No rigorous connection is made between the number of parameters and the number of data points ; vi) No attention was paid to the necessity of the stability of parameter values ; vii) The experiments were performed in the past in quite a chaotic way : huge gaps are sometimes present between low-energy and high-energy domains or inside the high-energy domain itself.

The COMPETE (COmputerized Models and Parameter Evaluation for Theory and Experiment) project tries to cure as much as possible the above discussed arbitrariness.

The $`\chi ^2/dof`$ criterium is not able, by itself, to cure the above difficulties : new indicators have to be defined. Once these indicators are defined, an appropriate sum of their numerical values is proposed in order to establish the rank of the model under study : the highest the numerical value of the rank the better the model under consideration.

The final aim of the COMPETE project is to provide our community with a periodic cross assessments of data and models via computer-readable files on the Web .

We consider the following exemplar cases of the imaginary part of the scattering amplitudes :

$$ImF^{ab}=s\sigma _{ab}=P_1^{ab}+P_2^{ab}+R_+^{ab}\pm R_{}^{ab}$$
(1)
where :
\- the $`\pm `$ sign in formula (1) corresponds to antiparticle (resp. particle) - particle scattering amplitude.
\- $`R_\pm `$ signify the effective secondary-Reggeon ($`(f,a_2),(\rho ,\omega )`$) contributions to the even (odd)-under-crossing amplitude

$$R_\pm (s)=Y_\pm \left(\frac{s}{s_1}\right)^{\alpha _\pm },$$
(2)
where $`Y`$ is a constant residue, $`\alpha `$ \- the reggeon intercept and $`s_1`$ \- a scale factor fixed at 1 GeV2 ;
\- $`P_1(s)`$ is the contribution of the Pomeron Regge pole

$$P_1^{ab}(s)=C_1^{ab}\left(\frac{s}{s_1}\right)^{\alpha _{P_1}},$$
(3)
$`\alpha _{P_1}`$ is the Pomeron intercept $`\alpha _{P_1}=1,`$ and $`C^{ab}`$ are constant residues.
\- $`P_2^{ab}(s)`$ is the second component of the Pomeron corresponding to three different $`J`$-plane singularities :

* a Regge simple - pole contribution

$$P_2^{ab}(s)=C_2^{ab}\left(\frac{s}{s_1}\right)^{\alpha _{P_2}},$$
(4)
with $`\alpha _{P_2}=1+ϵ,ϵ>0,`$ and $`C_2^{ab}`$ const. ;
* a Regge double-pole contribution

$$P_2^{ab}(s)=s\left[A^{ab}+B^{ab}\mathrm{ln}\left(\frac{s}{s_1}\right)\right],$$
(5)
with $`A^{ab}`$ and $`B^{ab}`$ const. ;
* a Regge triple-pole contribution

$$P_2^{ab}(s)=s\left[A^{ab}+B^{ab}\mathrm{ln}^2\left(\frac{s}{s_0}\right)\right],$$
(6)
where $`A^{ab}`$ and $`B^{ab}`$ are constants and $`s_0`$ is an arbitrary scale factor.

We consider all the existing forward data for $`pp,\overline{p}p,\pi p,Kp,\gamma \gamma `$ and $`\mathrm{\Sigma }p`$ scatterings. The number of data points is : 904, 742, 648, 569, 498, 453, 397, 329 when the cut-off in energy is 3, 4, 5, 6, 7, 8, 9, 10 GeV respectively. A large number of variants were studied and classified. All definitions and numerical details can be found in Ref. 1.

The 2-component Pomeron classes of models are RRPE, RRPL and RRPL2, where by RR we denote the two effective secondary-reggeon contributions, by P - the contribution of the Pomeron Regge-pole located at $`J=1`$, by E - the contribution of the Pomeron Regge-pole located at $`J=1+ϵ`$, by L - the contribution of the component of the Pomeron, located at $`J=1`$ (double pole), and by L2 - the contribution of the component of the Pomeron located at $`J=1`$ (triple pole). We also studied the 1-component Pomeron classes of models RRE, RRL and RRL2.

The highest rank are get by the RRPL2u models (see Table 1 and Figure 1), corresponding to the $`\mathrm{ln}^2s`$ behaviour of total cross sections first proposed by Heisenberg 50 years ago .

The $`u`$ index denotes the universality property : the coupling $`B`$ of the $`\mathrm{ln}^2s`$ term is the same in all hadron-hadron scatterings and $`s_0`$ is the same in all reactions. We note that the coupling $`B`$ is remarkably stable for the different models : 0.3157 mb (RRPL2u(19)), 0.3152 mb (RRPL2u(21)), 0.3117 mb ((RR)cPL2u(15)), etc. This reinforces the validity of the universality property.

We note also that the familiar RRE Donnachie-Landshoff model is rejected at the 98% C.L. when models which achieve a $`\chi ^2/dof`$ less than 1 for $`\sqrt{s}5`$ GeV are considered.

The predictions of the best RRPL2u model, adjusted for $`\sqrt{s}5`$ GeV, are given in Tables 2-4.

The uncertainties on total cross sections, including the systematic errors due to contradictory data points from FNAL (the CDF and E710/E811 experiments, respectively), can reach $`1.9\%`$ at RHIC, $`3.1\%`$ at the Tevatron, and $`4.8\%`$ at the LHC, whereas those on the $`\rho `$ parameter are respectively $`5.4\%`$, $`5.2\%`$, and $`5.4\%`$. The global picture emerging from fits to all data on forward observables supports the CDF data and disfavors the E710/E811 data at $`\sqrt{s}=1.8`$ TeV.

Any significant deviation from the predictions based on model RRPL2u will lead to a re-evaluation of the hierarchy of models and presumably change the preferred parametrisation to another one. A deviation from the “allowed region” would be an indication that strong interactions demand a generalization of the analytic models discussed so far, e.g. by adding Odderon terms, or new Pomeron terms, as suggested by QCD.

## 2 A generalization of the Heisenberg model for the total cross-section

In his remarkable paper of 1952, Heisenberg investigated production of mesons as a problem of shock waves. One of his results was that the total cross section increases like the square of the logarithm of the centre-of-mass energy. It is noteworthy that this result coincides with very recent calculations based on AdS/CFT dual string-gravity theory  or on the Colour Glass Condensate Approach  and, of course, saturates the Froissart-Martin bound . In contradistinction to the latter case however the coefficient of the $`\mathrm{ln}^2`$ term is an estimate at finite energies and not an asymptotic bound as the one obtained by Lukaszuk and Martin .

We have shown that  by modifications of the original model of Heisenberg motivated by the enormous progress of knowledge in the 50 years that passed thence, the model yields some general and even some quantitative results which describe the data very well.

The considerations of Heisenberg concerning the total cross section are essentially geometrical ones, but the crucial ingredient is that the energy density and not the hadronic density is the essential quantity to be taken into account.

Proton-proton collisions are considered in the centre-of-mass system and the energy $`\sqrt{s}`$ is supposed to be high enough that Lorentz contraction allows to view the nucleons as discs.

Interaction takes place only in the overlap region and the crucial assumption is made that a reaction can only occur if the energy density is high enough in order to create at least a meson pair.

The result of Heisenberg is

$$\sigma =\frac{\pi }{m^2}\mathrm{ln}^2\frac{\sqrt{s}}{k_0}.$$
(7)
We see that implicitely the assumption has been made that if a meson production is energetically possible, it will happen (black disk). Of course, Heisenberg was taking the pion mass for the meson mass. For the energy of the produced mesons he deduced, in his dynamical considerations, assuming interactions of maximal strength, that the energy $`k_0`$ (for two produced mesons) increases only slowly with energy, at any rate not by a power of $`s`$. Therefore the asymptotically leading term in the cross section is $`(\pi /4m_\pi ^2)\mathrm{ln}^2s`$, the coefficient $`\pi /4m_\pi ^2`$ being 1/4 of the Lukaszuk-Martin bound. The argument can be extended easily to hadron-hadron scattering in general, and therefore we have the result that the coefficient of the $`\mathrm{ln}^2s`$ term is universal for all hadron reactions.

There are two obvious necessary modifications of the Heisenberg model :

* If we want to apply it to all kind of hadrons, we have to take care of the different hadron sizes, since in the above treatment all sizes are equal to $`1/m`$.
* We have to take into account that direct pion exchange, though being the exchange with the lightest particle, is not relevant at high energies. This is due to the fact that exchanged gluons have spin 1 and pions spin 0. Therefore already in Born approximation gluon exchange dominates at high energies. In Regge theory this is manifested by the fact that intercept of the pion is much lower than that of the Pomeron. For the mass we rather insert a mass $`M`$ in the range of the glueball mass instead of the pion mass $`m`$, since we believe that the high-energy behavior is dominated by gluon exchange.

We then obtain

$`\sigma `$ $`=`$ $`{\displaystyle \frac{\pi }{4M^2}}\mathrm{ln}^2s`$ (8)
$`+{\displaystyle \frac{\pi }{M}}\mathrm{ln}s\left\{(R_1+R_2)+{\displaystyle \frac{1}{M}}\mathrm{ln}{\displaystyle \frac{\alpha }{k_0}}\right\}`$
$`+\pi (R_1+R_2)^2+{\displaystyle \frac{\pi }{M^2}}\mathrm{ln}^2{\displaystyle \frac{\alpha }{k_0}}`$
$`+{\displaystyle \frac{2\pi }{M}}(R_1+R_2)\mathrm{ln}{\displaystyle \frac{\alpha }{k_0}},`$
where $`0\alpha 1`$.

We see that the leading $`\mathrm{ln}^2s`$ term is still universal, but now dominated rather by a glueball than by the pion mass. Since $`R_1`$ and $`R_2`$ are supposed to be of the size of the electromagnetic radii, the second term in Eq. (8) will dominate over the $`{\displaystyle \frac{\pi }{M^2}}`$ term except at high energies, $`sk_0^2/\alpha ^2`$.

In order to perform a rough numerical estimate, we may insert for the glueball mass a value between 1.4 and 1.7 GeV, yielding

$$\frac{\pi }{4M^2}=0.110.16\text{mb}.$$
(9)
For $`R_1`$ and $`R_2`$ we may insert the electromagnetic radii. In contrast to Heisenberg, we insert for $`k_0`$ the minimal energy of two produced particles. Since production seems to occur in clusters with mass around 1.3 GeV , we can put $`k_0=2.6`$ GeV. The value of $`\alpha (0\alpha 1)`$ might be process dependent. For very small objects (”onia”) it might be very small.

In the past, application of the Heisenberg model to the global analyses of the forward hadronic data were performed in , but the universality of the leading term was not discussed there. This universality was treated by Gershtein and Logunov , who made the assumption, as in the present paper, that the growth of the hadron-hadron total cross-sections is related to resonance production of glueballs.

The COMPETE value of $`B`$ (see the previous section) corresponds to a mass $`M`$ of 1 GeV, a bit small for a glueball, but not unreasonable given the crude approximations. $`Z^{HH}`$ (where $`Z^{HH}=C_1^{HH}+A^{HH}`$) are in the right order of magnitude of $`R^2`$.

A consequence of the universal $`\mathrm{ln}^2s`$ term is that at asymptotic energies all hadron cross sections become equal. At finite but high energies the pion and kaon proton cross sections are therefore expected to rise somewhat faster than the nucleon-nucleon cross sections. This seems indeed to be indicated by the data.


text:
# Phantom Cosmologies

## I  INTRODUCTION

The combined analysis of SNe Ia observations 1 ; 2 , galaxy cluster measurements 3  and the latest CMB data 4  provides compelling evidence for the existence of dark energy which dominates the present day universe and accelerates the cosmic expansion. The recent detection of Integrated Sachs-Wolfe effect 5  also gives a strong and independent support to dark energy. In principle, any physical field with positive energy density and negative pressure, which violates the strong energy condition, may cause the dark energy effect of repulsive gravitation. Of late, phantom fields 6  have emerged as potential candidates for dark energy. Scalar fields with super-negative equation of state $`(p=w\rho ,w<1)`$ are called phantom fields as their energy density increases with the expansion of the universe in contrast to quintessence energy density $`(w>1)`$ which scales down with the cosmic expansion. The phantom models violate the dominant energy condition $`(p+\rho )<0`$ as such they may not be physically stable models of dark energy; but, strangely enough, phantom energy is found to be compatible with most of the classical tests of cosmology 6  based on current data from SNe Ia observations, CMB anisotropy and mass power spectrum.

The peculiar nature of phantom energy, violation of dominant energy condition and its strange consequences, possible rip-off of the large and small scale structures of matter, occurrence of future singularity and probable decay of phantom e nergy have attracted many cosmologists 7 -24 52 ; 53 ; 54  and made ‘phantom cosmologies’ a hot topic of research.

In section 2, we have discussed the dynamics of minimally coupled scalar fields with special reference to phantom fields. There is extensive literature in cosmology on rolling scalar fields 22 ; 23 , quintessence fields and tracker fields 24 -40 50 ; 51 , the cosmological constant $`\mathrm{\Lambda }`$ 25 ; 41 ; 42 ; 49  and other forms of dark energy. The major problem in cosmology is to identify the form of dark energy that dominates the universe today whether it is phantom energy, quintessence, simply $`\mathrm{\Lambda }`$ or something else. Maor, Brustein and Steinhardt 43  have discussed the degeneracy in the measurement of the dark energy parameter $`w`$ from SNe Ia data, its time variation and pitfalls in taking $`w`$ to be a constant.

In section 3, we have computed the present age $`t_0`$ of the universe in two steps. First we calculate the expansion age up to the end of matter-dominated era, denoted by $`t_m`$. In order to supplement it with the expansion age during dark energy dominated era, we express $`t_0`$ in terms of $`t_m`$ and thereby we calculate $`t_0`$. The advantage of this method is that the expansion age in the two segments can be expressed separately in terms of the redshift $`z_m`$ at the end of matter-dominated era which is again a function of the parameter $`w`$.

In section 4, assuming $`w`$ to be constant, we have discussed a kind of degeneracy in the value of $`w(z_m)`$ which leads to duality in the behavior of phantom and quintessence models with respect to transition redshift from deceleration to accelerating phase of expansion. In fact two distinct values of parameter $`w`$, usually one lying in the range of quintessence field and another in the range of phantom field, lead to the same transition redshift $`z_m`$.

In section 5, we have tried to constrain the range of the dark energy parameter $`w`$ on the basis of data analysis of Kiselev 44 , Freedman and Turner 45  Schubnell 46  and the precise observational data from WMAP 47  in combination with SDSS 48 . In section 6, we conclude with some remarks on phantom energy.

## II DYNAMICS OF PHANTOM COSMOLOGY

Consider a 2-component cosmic fluid in a Friedmann universe comprising (i) pressure-free matter of energy density $`\rho _m`$ and (ii) a minimally coupled scalar field of energy density $`\rho _x`$ and equation of state $`p=w\rho `$ which contributes to dark energy in the universe.The energy densities $`\rho _ma^3`$ and $`\rho _xa^{3(1+w)}`$ evolve independently in the expanding universe.

The dark energy might be due to

(a) quintessence field if $`1<w<\frac{1}{3}`$

(b) cosmological constant $`\mathrm{\Lambda }`$ if $`w=1`$

or (c) phantom field if $`w<1`$.

In the above classification, the equation of state parameter $`w`$ plays the role of dark energy parameter.

The Friedmann equations are

$$\frac{\dot{a}^2}{a^2}=\frac{8\pi G}{3}\left[\rho _m+\rho _x\right]=H_0^2\left[\mathrm{\Omega }_m^0(a_0/a)^3+\mathrm{\Omega }_x^0(a_0/a)^{3(1+w)}\right]$$
(1)
and

$`{\displaystyle \frac{\ddot{a}}{a}}`$ $`=`$ $`{\displaystyle \frac{4\pi G}{3}}[\rho _m+\rho _x(1+3w)]={\displaystyle \frac{4\pi G}{3}}\rho _x[\mathrm{\Omega }_x^1+3w]`$ (2)
$`=`$ $`{\displaystyle \frac{4\pi G}{3}}\rho _x\left[{\displaystyle \frac{\mathrm{\Omega }_m^0}{\mathrm{\Omega }_x^0}}\left({\displaystyle \frac{a_0}{a}}\right)^{3w}+1+3w\right]`$
The cosmic expansion decelerates as long as

$$\mathrm{\Omega }_x^1+3w>0$$
(3)
With the growth of phantom energy density parameter, $`\mathrm{\Omega }_x^1`$ goes on decreasing with time until the transition to accelerating phase takes place at cosmic time $`t=t_m`$. The transition epoch $`t_m`$ corresponds to the red-shift $`z_m`$ given by

$$1+z_m=\left[\frac{(3w+1)\mathrm{\Omega }_x^0}{\mathrm{\Omega }_m^0}\right]^{\frac{1}{3w}}$$
(4)
It may be emphasized that the transition epoch $`t_m`$ marks the end of the ’effective matter dominated’ era or the beginning of the accelerating phase in the cosmic expansion after which the large scale structure formation in the universe must cease although $`\mathrm{\Omega }_m`$ may still be greater than $`\mathrm{\Omega }_x`$. It is evident from Eq.(4) that the transition epoch $`t_m`$ depends upon the choice of $`w`$. FIG 3 shows the variation of $`z_m`$ with $`w`$.

As $`\mathrm{\Omega }_x^1`$ decreases further after the transition epoch

$$\mathrm{\Omega }_x^1+3w<\mathrm{\hspace{0.17em}0}$$
(5)
and the Hubble expansion in the accelerating phase of the universe is given by Eq.(1)

$$\frac{\dot{a}}{a}=H_0\sqrt{\mathrm{\Omega }_x^0}\left(\frac{a_0}{a}\right)^{3(1+w)/2}\left[1+\frac{\mathrm{\Omega }_m^0}{\mathrm{\Omega }_x^0}\left(\frac{a}{a_0}\right)^{3w}\right]^{\frac{1}{2}}$$
(6)
Expanding binomially and integrating, we get

$$\chi t+c=\left(\frac{a}{a_0}\right)^{3(1+w)/2}\left[\frac{2}{3(1+w)}\frac{\mathrm{\Omega }_m^0(a)^{3w}}{\mathrm{\Omega }_x^0(a_0)^{3w}(9w+3)}+\mathrm{}\right]$$
(7)
where $`\chi =H_0\sqrt{\mathrm{\Omega }_x^0}`$. The binomial expansion of the right hand side of Eq.(6) is valid under the condition $`\frac{\mathrm{\Omega }_m^0}{\mathrm{\Omega }_x^0}(\frac{a}{a_0})^{3w}<1<(3w+1)`$ which holds during the accelerating phase over the range (i)$`a<a_0`$, $`w<0`$ and $`\mathrm{\Omega }_m^0<\mathrm{\Omega }_x^0`$ and (ii) $`a<a_0`$, $`w<0`$ and $`(1+z)^{3w}<\frac{(3w+1)\mathrm{\Omega }_x^0}{\mathrm{\Omega }_m^0}`$. This ensures that the accelerated expansion during the regime $`\mathrm{\Omega }_m>\mathrm{\Omega }_x`$ continues as long as $`(3w+1)\mathrm{\Omega }_x^0>\mathrm{\Omega }_m^0`$.

Since $`w<0`$, the successive terms on the right hand side of Eq.(7) decrease by O(3w) of magnitude and the scale factor is effectively given by

$$a^{3(1+w)/2}(t)=\frac{3(1+w)}{2}\chi t+c$$
(8)
It shows that the contribution of the matter density is almost negligible during phantom dominated universe since the contribution of $`\mathrm{\Omega }_m^0`$ falls down steeply by 3 orders of magnitude or more with each successive term in Eq.(7). According to Eq.(1)also, the cosmic expansion in the accelerating phase is essentially driven by the phantom field since its energy density scales up as $`\rho _xa^{3(1+w)}`$ whereas $`\rho _m`$ scales down as $`a^3`$.

During the deceleration phase, the Hubble expansion is given by

$$\frac{\dot{a}}{a}=H_0\sqrt{\mathrm{\Omega }_m^0}\left(\frac{a_0}{a}\right)^{3/2}\left[1+\frac{\mathrm{\Omega }_x^0}{\mathrm{\Omega }_m^0}\left(\frac{a_0}{a}\right)^{3w}\right]^{1/2}$$
(9)
The deceleration condition (3) implies that $`\frac{\mathrm{\Omega }_x^0}{\mathrm{\Omega }_m^0}(\frac{a_0}{a})^{3w}<\frac{1}{3w+1}<1`$ since $`w<1`$.

Therefore expanding Eq.(9) binomially and integrating we get

$$\xi t=a^{3/2}\left[\frac{2}{3}\frac{(1+z)^{3w}\mathrm{\Omega }_x^0}{3(12w)\mathrm{\Omega }_m^0}+\frac{(1+z)^{6w}}{4(14w)}\left(\frac{\mathrm{\Omega }_x^0}{\mathrm{\Omega }_m^0}\right)^2+\mathrm{}\right]$$
(10)
Since $`w<0`$, the first term on the right hand side of Eq.(10) dominates while the remaining terms decrease for high redshifts. Therefore during the matter dominated era, the scale factor is given by

$$a^{3/2}(t)=\frac{3}{2}\xi t$$
(11)
Eq.(11) holds at the epoch $`t=t_m`$, as such

$$a^{3/2}(t_m)=\frac{3}{2}\xi t_m$$
(12)
At the beginning of the dark energy dominated phase, Eq.(8) gives

$$a^{3(1+w)/2}(t_m)=\frac{3}{2}\chi t_m+c$$
(13)
Matching the junction conditions at $`t=t_m`$, Eqs,(8),(12) and (13) yield

$$\left[\frac{a(t)}{a(t_m)}\right]^{3(1+w)/2}=\mathrm{\hspace{0.17em}1}+\frac{3/2(w+1)\chi (tt_m)}{a^{3(1+w)/2}(t_m)}=1+(w+1)\frac{tt_m}{t_m}$$
(14)
Therefore the scale factor in the phantom (dark energy) dominated universe is given by

$$a(t)=\frac{a(t_m)}{[w+(1+w)t/t_m]^{\frac{2}{3(1+w)}}}fort>t_m$$
(15)
Since $`1+w<0`$ for phantom fields, $`a(t)`$ diverges to infinity at $`t^{}=\frac{w}{1+w}t_m`$. Prior to blow over time $`t^{}`$, $`H>0`$ and the deceleration parameter $`q=1+3(1+w)/2`$ in contrast to $`q_0=1/2+3/2w\mathrm{\Omega }_x`$ at the present epoch. On the contrary in the quintessence dominated universe$`(1+w>0)`$, the cosmic expansion is singularity-free with the scale factor $`a(t)=a(t_m)\left[1+(1+w)(tt_m)/t_m\right]^{\frac{2}{3(1+w)}}`$.

Using Eq.(15), the energy density of the phantom universe $`(t>t_m)`$ is given by

$$\rho _x(t)=\frac{\rho (t_m)}{[w+(1+w)t/t_m]^2}$$
(16)
Accordingly $`\rho _x`$ goes on increasing gradually, followed by steep rise to infinite value at $`t^{}=\frac{w}{1+w}t_m`$. Therefore the phantom models have a finite life time ending in a singularity. On the other hand, $`\rho _x`$ scales down with time in the quintessence universe whereas $`\rho _\mathrm{\Lambda }`$ remains stationary. For $`t<t_m`$, the Hubble expansion is dominated by matter density; accordingly $`at^{2/3}`$, , $`\rho _mt^2`$ but $`\rho _x`$ varies independently as $`t^{2(1+w)}`$ as shown in FIG. 2.

The expansion and density singularity in the phantom universe corresponds to the curvature singularity as the Ricci scalar $`R_i^i`$ also tends to infinity at this epoch. The total age $`t^{}`$ of the phantom universe depends upon the choice of $`w`$ and is larger for values of $`w`$ closer to $`1`$. For example, $`t^{}=21t_m`$ for $`w=1.05`$ whereas $`t^{}=6t_m`$ for $`w=1.2`$

## III PRESENT AGE OF THE UNIVERSE

The present age $`t_0`$ of the universe depends on the cosmological density parameters and the equation of state parameter $`w`$ which determines the redshift $`z_m`$ at the transition epoch $`t_m`$ (at the end of matter dominated phase). We calculate $`t_0`$ in two steps. During the matter-dominated phase, the Hubble expansion is given by Eq.(9). Integrating the first term in the binomial expansion of Eq.(9) over the redshift range from infinity to $`z_m`$, we calculate the age of the universe from the beginning to the end of matter-dominated era as given by the expression

$$t_m=(H_0\sqrt{\mathrm{\Omega }_m^0})^1\left[\frac{2}{3}(1+z_m)^{3/2}\right]$$
(17)
where $`\mathrm{\Omega }_m^0=0.27,\mathrm{\Omega }_x^0=0.73`$.

Using Eq.(15), the present age $`t_o`$ is given by the equation

$$t_0=\left[1+\frac{(1+z_m)^{3(1+w)/2}1}{1+w}\right]\times t_m$$
(18)
Combining Eq.(17) and Eq.(18), the present age of the universe is calculated for a wide range of $`w`$ as shown in the Table 1.

## IV DUALITY IN QUINTESSENCE AND PHANTOM BEHAVIOR

We have investigated the correlation between the transition redshift $`z_m`$ (corresponding to the end of matter-dominated era) and the dark energy parameter $`w`$ and found a sort of duality in the behavior of quintessence fields (Q) and phantom fields (P) \[see FIG. 3\]. For every value of $`z_m`$ in FIG. 3, the corresponding parameter $`w`$ has,in general, two values, one lying in the range of Q fields and the other in the range of P fields, both leading to cosmological parameters (like the present age $`t_0`$ and the deceleration parameter $`q_0`$) which seem to be compatible with the observational data. For example, if the transition from the decelerating phase to accelerating phase occurs at $`z_m=0.739`$, the corresponding dark energy parameter may be either $`w=0.8`$ (Quintessence field) or $`w=1.1`$ (phantom field) which yield $`14.37`$ Gyr and $`13.58`$ Gyr respectively for the present age of the universe. Hence both the Q model and the P model seem to be compatible with the observational value 45  of $`t_0`$. This duality poses a question whether the phantom fields really exist or they are merely ghost fields which replicate the quintessence-like behavior for super-negative equation of state, violating the dominant energy condition. It might explain the concordance of SNIa and galaxy cluster abundance observations in the extended $`w\mathrm{\Omega }_m`$ parameter space for $`w<1`$ (phantom models) as shown by Caldwell et al 7 .

The above-mentioned duality is essentially different from the form-invariance symmetry between standard cosmology and phantom cosmology pointed out independently by Dabrowski et al 24  and Chimento et al 55 . They have shown that in the case of a single component cosmological model, the scale factors of the standard and phantom models for a given energy density have reciprocal relationship with the equation of state parameter $`w_{ph}=w2,(w>1)`$. This formalism seems inadequate to describe the evolution of the scale factor in dark energy models in conjunction with pressure-free matter.

## V CONSTRAINTS OVER PHANTOM COSMOLOGIES

One of the greatest challenges in cosmology is to understand the nature of the dark energy. Dark energy models are characterized by two parameters $`\mathrm{\Omega }_x`$ and $`w`$. From the analysis of WMAP data 47 , $`\mathrm{\Omega }_x=0.73\pm 0.04`$ is known up to high precision but $`w<0.8`$ (95 % cl) leaves the field open to speculation whether the dark energy is phantom energy or quintessence energy. According to Melchiorri’s analysis 13 , $`1.38<w<0.82`$. We have examined the possibility of constraining the range of $`w`$ by comparison of the theoretically calculated age of the universe in Table 1 with the age derived from the observational data from WMAP, SDSS and data analysis of Tegmark et al 48 . By assuming $`H_0=\mathrm{\hspace{0.17em}71}_{3.0}^{+4.0}`$ and $`t_0=\mathrm{\hspace{0.17em}13.7}\pm 0.2`$ Gyr from WMAP data, we can find a narrow range $`1.18<w<0.93`$ for dark energy parameter and the corresponding range $`0.8<q_0<0.52`$ for the deceleration parameter (consistent with Kisilev’s analysis 44 ) but it would be more realistic to allow for errors in the measurement of Hubble constant and take $`H_0=\mathrm{\hspace{0.17em}72}\pm 7.0andt_0=\mathrm{\hspace{0.17em}13.0}\pm 1.5Gyr`$ 45  but this yields a wide range $`1.5<w<0.75`$ for variation of dark energy parameter and a very loose constraint on $`w`$ and $`q_0`$.

However the degeneracy in the determination of of $`w(z_m)`$ for a chosen transition redshift may be used to constrain $`w`$ and discriminate between the quintessence and phantom dark energy. Since the formation of the galactic clusters ceases with the end of matter dominated era at redshift $`z_m`$, the lowest redshift observations of the galactic clusters can indicate the maximum value of the redshift at which large scale structure formation in the universe would stop. In general there might be two values of $`w`$ corresponding to this particular value of $`z_m`$ out of which the one satisfying the WMAP range for the present age of the universe $`(t_0=13.7\pm 0.2)`$, may be taken as the correct value of $`w`$ to determine the genuine candidate of dark energy filling the universe.

## VI CONCLUDING REMARKS

The combined analysis of the latest cosmological observations provides a definite clue of the existence of dark energy in the universe but it is difficult to distinguish between the various forms of dark energy at present. So far as phantom energy is concerned, it is found to be compatible with SNe Ia observations and CMB anisotropy measurements but the violation of the ‘dominant Energy condition’ makes phantom models physically unstable; however, phantom models may be considered to be phenomenologically viable provided their age happens to be less than the time scale of the singularity. In case the instability occurs earlier and the dark energy decays into gravitons, as discussed by Carroll et al 14 , the universe might escape the ordeal of ‘Rip-Off’ 7  and phantom singularity. With the large number of SNe Ia observations expected from SNAP, LOSS and other surveys in the coming years, a clear picture of the dark energy profile is likely to emerge which would reveal whether $`\mathrm{\Lambda }`$CDM, QCDM or PCDM is the correct cosmology of the universe.

###### Acknowledgements.

This work was done under the SERC Project, submitted to DST, India. The author gratefully acknowledges the hospitality extended by W. I. Fine Theoretical Physics Insitute, University of Minnesota, useful discussions with Keith Olive and valuable help of Vassilis Spanos and Roman Zwicky in the preparation of this manuscript.


text:
# Light transmission assisted by Brewster-Zennek modes in chromium films carrying a subwavelength hole array

## I Introduction

The optical properties of thin metallic films containing periodic arrays of subwavelength holes have been actively studied in relation with the possibility of transmission engineering. Originally, the remarkable optical transmission found in these structures was pointed out by T.W. Ebbesen et al back in 1998T.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff (1998). After these early experimental observations, the role of the thin metallic film surface plasmons (SPs) was put forward in order to explain the peculiar wavelength dependence of the transmissionT.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff (1998); U. Schröter, D. Heitmann (1998); T. Thio, H.F. Ghaemi, H.J. Lezec, P.A. Wolff, T.W. Ebbesen (1999); H.F. Ghaemi, T. Thio, D.E. Grupp, T.W. Ebbesen, H.J. Lezec (1998); J.A. Porto, F.J. Garcia-Vidal, J.B. Pendry (1999); D.E. Grupp, H.J. Lezec, T.W. Ebbesen, K.M. Pellerin, T. Thio (2000); T. Thio, H.J. Lezec, T.W. Ebbesen (2000); Krishnan et al. (2001); Q. Cao, Ph. Lalanne (2002). As a consequence of this suggestion, it was generally believed that the typical transmission features could only be obtained with metallic filmsT.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff (1998); D.E. Grupp, H.J. Lezec, T.W. Ebbesen, K.M. Pellerin, T. Thio (2000). However it is important to note that other explanations had been developed which could compete with the SPs modelJ.-M. Vigoureux (2001); M.M.J. Treacy (2002); Q. Cao, Ph. Lalanne (2002); H.J. Lezec, T. Thio (2004); M. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003); C. Genet, M. P. van Exter and J. P. Woerdman (2003). For instance, it was suggested that these phenomena could also be described in terms of the short range diffraction of evanescent wavesJ.-M. Vigoureux (2001), or in terms of dynamical diffraction effectsM.M.J. Treacy (2002). Another explanation suggests to emphasize the role of cavity resonances taking place in the holes to explain the transmission enhancementQ. Cao, Ph. Lalanne (2002).

Many questions remain to completely clarify the scattering processes involved in these experiments. Nevertheless, even if the exact role of SPs is not clearly assessedT.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff (1998); U. Schröter, D. Heitmann (1998); T. Thio, H.F. Ghaemi, H.J. Lezec, P.A. Wolff, T.W. Ebbesen (1999); H.F. Ghaemi, T. Thio, D.E. Grupp, T.W. Ebbesen, H.J. Lezec (1998); J.A. Porto, F.J. Garcia-Vidal, J.B. Pendry (1999); D.E. Grupp, H.J. Lezec, T.W. Ebbesen, K.M. Pellerin, T. Thio (2000); T. Thio, H.J. Lezec, T.W. Ebbesen (2000); Krishnan et al. (2001); Q. Cao, Ph. Lalanne (2002), most interpretations tend to admit that SPs should play a key role in the observed transmission spectrum features. The observed transmission exhibits a set of peaks and dips. Many authors suggested that transmission peaks are created by SPs resonancesT.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff (1998); U. Schröter, D. Heitmann (1998); T. Thio, H.F. Ghaemi, H.J. Lezec, P.A. Wolff, T.W. Ebbesen (1999); H.F. Ghaemi, T. Thio, D.E. Grupp, T.W. Ebbesen, H.J. Lezec (1998); J.A. Porto, F.J. Garcia-Vidal, J.B. Pendry (1999); D.E. Grupp, H.J. Lezec, T.W. Ebbesen, K.M. Pellerin, T. Thio (2000); T. Thio, H.J. Lezec, T.W. Ebbesen (2000); Krishnan et al. (2001) but others rather saw the SPs resonances at the spectral dipsQ. Cao, Ph. Lalanne (2002). In a recent workM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003) we suggested to redefine the role of SPs in Ebbesen experiments in the required context of resonant Wood anomaliesWood (1935). In doing so, we have shown that the transmission spectrum could be better depicted as a series of Fano profilesFano (1938). These recognizable lineshapes result from the interference of non resonant transfers with resonant transfers which involve the film SPs and evanescent diffraction orders. We could then point out that each transmission peak-dip pair is nothing else than a Fano profileFano (1938). The appearance of an assymetric Fano lineshape does not necessarily locates the peak or the dip at the SP resonanceM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003). However, we have shown that the existence of the SPs is a condition for the presence of the Fano lineshapesM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003). It can be noted that recent results by C. Genet et al confirm this descriptionC. Genet, M. P. van Exter and J. P. Woerdman (2003).

As a further outcome of this workM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003) it became clear this kind of transmission spectrum, with its Fano lineshapes, could also be obtained in a more general context and that the generic concept of eigenmodes could substitute that of SPs. These results implied two important new ideas. First, according to the contributions of resonant and non-resonant processes in the Fano profile, eigenmodes can be associated with wavelengths closer to the peaks or to the dips. Second, it is possible to obtain transmission curves similar to those of metal films, by substituting SPs with guided modes or other types of polaritons. Many examples can be found, including highly refractive materials defining guided modes or ionic crystals in the restrahlen band defining phonon-polaritons.

In a previous paperM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003), it was made clear that a strong hypothesis on the origin of the eigenmodes was not required to account for the general features of the membrane transmission. As metal films were considered there, it seemed natural to involve SPs in the transmission mechanisms. By contrast, in another recent paperM. Sarrazin, J.-P. Vigneron (2003), simulations of light transmission through an array of subwavelength cylindrical holes in a tungsten layer deposited on glass substrate was given an interpretation based on dielectric guided modes. Indeed, in the wavelength domain under examination, tungsten exhibits a positive permittivity real partD.W. Lynch, W.R. Hunter (1991). Though SPs cannot exist, the transmission pattern is found very similar to that obtained with a metallic film. These theoretical results have recently been confirmed experimentallyH.J. Lezec, T. Thio (2004). Surprisingly, the very fact that the typical transmission pattern can be observed even in non-metallic systems convinced some authorsH.J. Lezec, T. Thio (2004) to fully reject the SPs hypothesis and rather consider models involving non-resonant evanescent waves diffraction. To our knowledge, this point of view is not supported by recent resultsM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003); C. Genet, M. P. van Exter and J. P. Woerdman (2003); M. Sarrazin, J.-P. Vigneron (2003) and the SPs must be replaced by other kinds of eigenmodesM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003); M. Sarrazin, J.-P. Vigneron (2003) in the non-metallic cases.

The purpose of the present paper is to demonstrate that the transmission profile of a chromium film, in the restricted wavelength domain ($`11121292`$ $`nm`$) where the dielectric constant is positiveD.W. Lynch, W.R. Hunter (1991), should involve eigenmodes which are not SPs nor guided modes. In this case, SPs can be substituted by Brewster-Zennek modes (BZ modes). This confirms that the observation of the effects described by T.W. Ebbesen et al in systems where SPs do not exist do not preclude the existence of a mechanism involving the excitation of other eigenmodes.

## II Brewster-Zennek modes

Let us recall some properties of surface modes, including the concept of BZ modesG. S. Agarwal (1973); F. Yang, J.R. Sambles, G.W. Bradberry (1991). We consider an interface between two media, $`1`$ and $`2`$, described by their respective permittivities $`\epsilon _1`$ and $`\epsilon _2`$. It is well know that for an isotropic material, only $`p`$ polarized modes can occur as a surface modes. We define $`k`$ as the wave vector component parallel to the interface and we set $`k_0=\omega /c`$, where $`\omega `$ is the angular frequency of the mode. Then, for the surface mode, the normal component of the wave vector is given in each medium by

$$\beta _i=\sqrt{k^2k_0^2\epsilon _i}.$$
(1)
The subscript $`i`$ denotes one of the media $`1`$ or $`2`$.

In the following we enforce a real positive value of $`\epsilon _1`$, which means that medium $`1`$ is a lossless dielectric material, while we admit the form $`\epsilon _2=\epsilon _2^{}+i\epsilon _2^{\prime \prime }`$ ($`\epsilon _2^{\prime \prime }0`$) for medium $`2`$, which means that this medium is allowed to generate losses. Applying electromagnetic boundary conditions pertinent to complex-response discontinuities, we are lead to the following well-known expression which controls the appearance of surface modesRaether (1988)

$$\epsilon _1\beta _2+\epsilon _2\beta _1=0$$
(2)
From this, and eq. (1), we can develop the following expression, out of which a dispersion relation can be drawn

$$k^2=k_0^2\frac{\epsilon _1\epsilon _2}{\epsilon _1+\epsilon _2}$$
(3)
However, it is important to realize that in this expression, the factors $`\beta _i`$ have been squared, so that eq. (3) and (2) are not equivalent. The solutions of eq. (2) can be found among those of eq. (3), but some of the solutions of eq. (3) are spurious and should be discarded. One way to do this is to take the possible values of $`k^2`$ from eq. (3), inject them in eq. (2), and check for the correct matching. This actually generates a second condition, which can be used together with eq. (3), to make precise the acceptable eigenmodes branches. This condition turns out to be

$$\sqrt{\frac{z^2}{1+z}}+z\sqrt{\frac{1}{1+z}}=0$$
(4)
where, somewhat surprisingly, the only variable that needs to be examined is the ratio $`z=\epsilon _2/\epsilon _1`$. The function of $`z`$, at the left-hand side of this equation, when carefully calculated, defines a border line which splits the complex plane into two extended connected regions. On one side of this border, the left-hand side of eq. (4) cannot vanish, while on the other side, it will vanish everywhere. The surface states can only show up when the ratio $`z`$ lies in this latter region.

The border separating the two regions can be obtained as follows. Both terms added in eq. 4 are complex numbers with the same modulus. In one of the regions mentioned above, they are opposite complex numbers and they cancel each other; in the rest of the complex plane, they turn out to be identical and they cannot add to zero. The borderline will then be found at the complex points $`z`$ where the sign determination of one of the square roots in eq. 4 changes, and this requires to look for the cut of the complex square-root function $`\sqrt{Z}`$. It is convenient to take this cut as the straight segment $`Z=i\gamma `$ ($`\gamma `$ real and positive) and this choice lead us to identify two distinct implicit paths,

$$\frac{1}{1+z}=i\gamma $$
(5)
and

$$\frac{z^2}{1+z}=i\gamma .$$
(6)
which form the border. This choice of the cut is dictated in part by the choice of a time-dependent oscillation written as $`\mathrm{exp}(i\omega t)`$ and the requirement that $`\sqrt{1}=1`$. Once the borderline is known, a quick check at the case $`z=1`$ shows that no cancellation can occur there, so that the active region where eigenmodes appear is, unambiguously, that region which does not contain $`z=1`$.

On the basis of this analysis, it is easily shown that, whatever $`\epsilon _2=\epsilon _2^{}+i\epsilon _2^{\prime \prime }`$ such that $`\epsilon _2^{}<0`$ and $`\epsilon _2^{\prime \prime }0`$ surface eigenmodes exist. These modes are SPs and they will appear at metal/dielectric interfaces regardless of the metal component of the heterojunction.

If $`\epsilon _2^{}0`$, surface eigenmodes exist only when $`0\epsilon _2^{}<\epsilon _1`$ if the condition

$$\epsilon _2^{\prime \prime }>\epsilon _2^{}\sqrt{\frac{\epsilon _1+\epsilon _2^{}}{\epsilon _1\epsilon _2^{}}}$$
(7)
is fulfilled. Obviously such modes are not SPs and it is well known that a dielectric/dielectric interface, without loss, cannot support a surface mode if $`\beta _i`$ is real. Nevertheless, surface modes can exist for a dielectric/dielectric interface if one of the media is a lossy material such that the imaginary part of its permittivity verify the above conditions. Such modes are known as Brewster-Zenneck modes. They verify $`Re\left\{k\right\}<\frac{\omega }{c}\sqrt{\epsilon _1}`$ and are necessarily related, by definition, to radiative surface modes (contrasting SPs in metals such that iron, gold or copper). Fundamentally it is well known that SPs consist in a collective motion of electronsRaether (1988). In the context of an electromagnetic model, the electronic character of SPs is then contained in the permittivity $`\epsilon _2`$. Indeed, the permittivity properties take their origin in the electronic characteristics of the mediumRaether (1988). This is true whatever the value of $`\epsilon _2^{}`$, positive or negative. So, despite the above-mentioned differences between SPs and BZ modes, both are related intrinsically to a collective electronic phenomenon. In this way, the transition from the SPs towards the BZ modes and vice versa occurs without discontinuities. Note that the present paper deals with metals only. In the case of ionic crystals in the far infrared, it would be natural to expect for the BZ modes counterpart associated with phonon-polaritons instead of SPs.

## III Brewster-Zennek-based transmission

Now we study the properties of the surface modes in the case of a chromium grating similar to those submitted to experiments by Thio et alT. Thio, H.F. Ghaemi, H.J. Lezec, P.A. Wolff, T.W. Ebbesen (1999) (see fig. 1). In a large part of the spectrum of interest here, the real part of the chromium permittivity is negative (fig. 2), as for simple metals below the plasma frequency. However, for wavelengths in the interval from $`1112`$ to $`1292`$ nm, it becomes positive, as, specifically, we observe $`0<\epsilon _2^{}<0,4`$ and $`25,54<\epsilon _2^{\prime \prime }<27,3`$. With vacuum and substrate (glass) permittivities equal to $`1`$ and $`2,24`$ respectively, the conditions for the formation of BZ modes are verified at both interfaces of the chromium film. This is a peculiar situation : with tungsten films in the dielectric domainD.W. Lynch, W.R. Hunter (1991); M. Sarrazin, J.-P. Vigneron (2003), for instance, the same conditions are not fulfilled and the BZ modes do not appear.

The transmission spectrum can easily be understood if we first localize the BZ modes in the spectrum of a flat, homogeneous film. For normal incident light falling on a such a homogeneous film considered periodic on a square lattice of parameter $`a`$ (the “empty” lattice case), surface modes resonance will appear at the wavelength $`\lambda _{i,j}`$ provided by the following equationRaether (1988)

$$\lambda _{i,j}=\frac{a}{(i^2+j^2)^{1/2}}Re\sqrt{\frac{\epsilon _1\epsilon _2}{\epsilon _1+\epsilon _2}}$$
(8)
where $`(i,j)`$ denotes the related vector $`𝐠`$ of the reciprocal lattice, such that $`𝐠=\left(2\pi \mathrm{/}a\right)(i,j)`$. These wavelengths are indicative only, as they are based on a model where the cylindrical holes have not been accounted for.

More detailed simulations which do not contain this simplification have been carried out. The calculations which will now be described are based on a coupled-modes method which combines a scattering matrix formalism with a plane wave representation of the fields. This technique provides a computation scheme for the amplitude and polarization ($`s`$ or $`p`$) of reflected and transmitted fields in any diffracted order. A brief account of the method has been presented elsewhereM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003); M. Sarrazin, J.-P. Vigneron (2004); J.P. Vigneron, F. Forati, D. André, A. Castiaux, I. Derycke, A. Dereux (1995) so that we will not recall any technicality about these computations.

Following the experimentsT. Thio, H.F. Ghaemi, H.J. Lezec, P.A. Wolff, T.W. Ebbesen (1999), fig. 3a shows the calculated zeroth order transmission of a grating prepared by perforating a $`100`$ nm thick chromium film deposited on glass with circular holes of $`250`$ nm radius (fig. 1). The two-dimensional square lattice parameter is $`a=1000`$ nm (solid line) or $`a=1200`$ nm (dashed line). The incident beam is normal to the grating surface and polarized so that the electromagnetic wave electric field is oriented along one of the holes nearest-neighbor directions ($`Oy`$ axis). A and B on Fig. 3a indicate the main transmission maxima. Fig. 3b represents the amplitude of the calculated resonant $`p`$-polarized diffraction orders $`(0,\pm 1)`$ for the substrate/metal interface (peak 1) and for the vacuum/metal interface (peak 2). Again, in Fig. 3b, the dashed line refers to a lattice parameter of $`a=1200`$ nm. It can be noted that, as justified by the above discussion, these peaks do not correspond to the transmission maxima. We have shown that the spectral lineshape of the resonance can be interpreted as a Fano’s profileM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003); C. Genet, M. P. van Exter and J. P. Woerdman (2003). What is remarkable here is that, for a lattice parameter of $`a=1200`$ nm, the resonance has moved to a spectral region ($`1200`$ nm) where the real part of the chromium permittivity is positive, and that clearly SPs do not exist there. In this situation the resonance is related to the appearance of a BZ mode. Fig 4 shows the location of both peaks (1) and (2) as a function of the grating parameter $`a`$ (solid lines (1) and (2), respectively). We could compare these values with the wavelengths obtained from eq. 8 (dashed lines (1) and (2) respectively) for the appropriate diffraction orders. A good agreement is found between these calculated quantities. We show that as the lattice parameter increases, the resonance wavelength increases linearly as predicted by eq. 8. The hatched domain in fig. 4 is characterized by the presence of BZ modes. The SP domain lies everywhere else. One interesting observation is that both kinds of surface modes (BZ modes and SPs modes) substitute each others in a continuous way in accordance with our previous explanation in section II. The location of transmission maxima A (circle dots) and B (square dots) are also shown on fig. 4. These, as already underlined, do not indicate directly the eigenmodes resonance wavelengths though these transmission maxima are known to be related to the surface modes via a Fano profile.

Note that, in a previous workH.F. Ghaemi, T. Thio, D.E. Grupp, T.W. Ebbesen, H.J. Lezec (1998), it was inferred that the SPs wavelength were red-shifted, compared to those calculated for an empty lattice. It was assumed that such a shift was needed to make the transmission maxima coincide with SPs wavelengths and the justification suggested for this shift was based on the presence of holes. This is not observed in actual detailed calculations. The present results show that the exact resonant wavelengths of the surface modes, as computed numerically, are not shifted enough to justify their location on transmission maxima. Nevertheless, as mentioned before, an eigenmode resonance does not necessarily coincide with a transmission maximum, but a maximum-minimum pair indicates the presence of a nearby eigenmode resonanceM. Sarrazin, J.-P. Vigneron, J.-M. Vigoureux (2003). In fact, in this case, the surface modes resonances are incidentally located closer to the minima of the transmission. It is interesting to note that Thio et al have compared experimentally the transmission from gratings similar to ours, with lattice parameters $`a=1000`$ nm and $`a=1200`$ nmT. Thio, H.F. Ghaemi, H.J. Lezec, P.A. Wolff, T.W. Ebbesen (1999). The present computations of the transmission matched very well with these dataT. Thio, H.F. Ghaemi, H.J. Lezec, P.A. Wolff, T.W. Ebbesen (1999). We take this as a strong evidence of the relevance of the Brewster-Zennek modes in the description of these experiments.

## IV Conclusion

In the context of studies of the optical transmission properties of subwavelength hole arrays, we have shown that Brewster-Zennek modes can substitude surface plasmons when considering a layer made from a lossy dielectric medium, for specific lattice parameters. When modelling a grating etched in a chromium film we can use the concept of Brewster-Zennek modes to explain the spectra found in previous experiments. This study recalls the existence requirement of eigenmodes for observing the typical transmission lineshapes and we have shown that, contrasting a frequent believe, surface modes wavelengths based on an empty lattice model are not significatively red-shifted to justify their location on transmission maxima.

###### Acknowledgements.

We acknowledge the use of the Namur Interuniversity Scientific Computing Facility (Namur-ISCF), a joint project between the Belgian National Fund for Scientific Research (FNRS), and the Facultés Universitaires Notre-Dame de la Paix (FUNDP). This work was partially supported by the EU Belgian-French INTERREG III project “PREMIO”, the EU5 Centre of Excellence ICAI-CT-2000-70029 and the Inter-University Attraction Pole (IUAP P5/1) on “Quantum-size effects in nanostructured materials” of the Belgian Office for Scientific, Technical, and Cultural Affairs.


text:
# STM observation of the quantum interference effect in finite-sized graphite

## 1 Introduction

A superperiodic pattern on metal surfaces has been observed on Cu(111) and Ag(111) surfaces by scanning tunneling microscopy (STM).ref1 ; ref2  This is interpreted as an interference wave pattern of free electrons that move on the Fermi surface originating from the surface and the bulk states, and are scattered by adatoms, step edges, and defects. Recently, a superperiodic pattern was also reported on semiconductor surfaces such as InAs/GaAs(111)A.ref3  This pattern is also explained by interference in the two-dimensional (2D) electron gas, which are electrons generated due to the band bending caused by a surface reconstruction and confined in a very thin layer near the surface. As for 2D structures, layered materials such as graphite are other candidates in generating the electronic wave interference effect due to a large 2D anisotropy. As an indicative fact of a large anisotropy of graphite, the foliation is often generated by heat treatment at a high temperature or by the direct collision of a STM tip against a graphite surface.ref4  The foliated nanographene sheets are occasionally folded back after tearing a graphite surface, and in one case they are adhered to the graphite surface with its edge structure rolled like the carbon nanotube. In this case, the nanographene sheets except for a nanotube-like part interact with the graphite substrate and sometimes a superperiodic pattern appears as a consequence of a defaulted stacking. Actually, if the foliated nanographene sheets are not folded completely and not adhered tightly to the graphite substrate surface, the interaction between the nanographene sheets and the substrate becomes very week, in the nanographene, where the standing wave is expected to generate due to the electron confinement in the 2D structured nanographene. From its anisotropy, other types of superperiodic patterns are also generated by stacking faults. One is a moire pattern related to a relative rotation between adjacent layers and another is a pattern originating from the dislocation-network structures due to abrupt transitions of stacking with the lattice distortion in plane.ref5 ; ref6 ; ref7  Though the periodicity of a moire pattern and a pattern originating from the dislocation-network structures can be determined by the angle of the relative rotation and the periodicity of glides, respectively, the corrugation amplitudes observed by STM cannot be explained just by the calculated local density of states (LDOS) of different stacked graphite.ref5 ; ref8  Instead, the scattering potential at the interface between the overlayer (a faulted part) and the substrate should give rise to the interference of the perturbed and the unperturbed waves in the overlayer, which affects the corrugation amplitudes of patterns dependent on the bias voltage of STM and the overlayer height from the substrate.ref9
In this paper, we report observed superperiodic patterns, whose periodicities vary gradually in plane, on nanographene sheets inclined with respect to a HOPG substrate, and superperiodic patterns originating from the dislocation-network structures, where the corrugation amplitudes depend on the bias voltage and the overlayer height from the substrate.

## 2 Electronic wave interference patterns
in nanographene sheets

A superperiodic pattern, whose periodicity changed gradually in plane, was observed by STM as shown in fig.1(a)-(b). The sample preparation is given elsewhere.ref5 ; ref11  The nanographene was fabricated by turning back after the heat-treatment-induced foliation at 1600 ˚C, which is ascertained from the fact that the shape of a nanographene edge resembles that of a step edge of the substrate as indicated by arrows in fig.1(a). The apparent height of the nanographene sheets from the HOPG substrate was about 0.8 nm from a cross-sectional analysis of the upper square region in fig.1(b). The nanographene probably consists of a stacking of two graphene layers, which interact weakly with the HOPG substrate, by comparison with the interlayer distance of bulk graphite (0.335 nm). The height of the lower square region in fig.1(b) was about 0.9 nm, therefore, the nanographene inclined from the substrate with an average slope of $`\mathrm{\Delta }z/\mathrm{\Delta }y`$$``$$`2\times 10^4`$. Note that the periodicity of the observed superperiodic pattern decreases in increasing the $`y`$ value. The corrugation amplitude of pattern also decreases in increasing the $`y`$ value from the cross-sectional analyses of the STM image in fig.1(b). Interestingly, the pattern in fig.1(b) changed into that in fig.1(c) by attachment of a nanographite flake on the nanographene sheets. By comparing two patterns in figs.1(b) and (c), it is evident that the periodicity changed remarkably at the upper part, though there was only a subtle change in the periodicity in the lower part. Such an effect on the patterns cannot be explained by a simple structural modulation. Instead, we can take the long-distance slope of nanographene sheets as the origin of the pattern and the interference effects of an electronic wave function due to a 2D structure of the nanographene interacting weakly with the substrate. In this case, the spatially varying periodicity can be explained by the electronic confinement effect in 2D-structured nanographnene with a potential gradient generated from the slope of the nanographene. To reproduce the observed patterns, we will characterize the patterns theoretically on the basis of the interference effect in a free electron model at first.ref12  By assuming a linear potential $`Fy`$ generated from the slope of the nanographene along the $`y`$ axis and a confinement effect in the $`x`$ axis due to the square-well potential within $`d/2`$$`<`$$`x`$$`<`$$`d/2`$, where $`d`$ is the width of the nanographene sheets along the $`x`$ axis, the Schrödinger equation is obtained:

$$[\frac{\mathrm{}^2}{2m}(\frac{^2}{x^2}+\frac{^2}{y^2})+V_xFy]\psi (x,y)=E\psi (x,y),$$
(1)
where $`V_x`$ is a square-well potential with infinite depth. The solution in the well is a standing wave in the $`x`$ axis and an oscillatory wave that is given by the Airy function in the $`y`$ axis. The potential gradient constant $`F`$, which was derived from fitting the LDOS to the observed pattern shown in fig.1(b) in the $`y`$ axis, was estimated to be about $`5\times 10^6`$ eV/nm, using the free electron mass. The potential variation for 200 nm along the $`y`$ axis was calculated to be about $`1\times 10^3`$ eV, which is unrealistically small because it is one order of magnitude smaller than thermal energy $`k_\text{B}T`$ at room temperature, where $`k_\text{B}T`$ is the Boltzmann constant. This inconsistency comes from the difference of electronic structures between the free electron model and actual 2D graphite. Then we use the $`kp`$ perturbation model, taking the band structure of 2D graphite with the linear dispersion near the Fermi energy into consideration for making a more realistic model. The Hamiltonian around the $`K`$-point with the linear potential is given as,

$$H=\left(\begin{array}{cc}Fy& i\gamma \frac{}{x}\gamma \frac{}{y}\\ i\gamma \frac{}{x}+\gamma \frac{}{y}& Fy\end{array}\right),$$
(2)
where $`\gamma `$ is $`(\sqrt{3}/2)a\gamma _0`$, $`a`$ and $`\gamma _0`$(=2.7 eV) are the nearest-neighbor distance and the resonance integral between nearest-neighbor carbon atoms, respectively. From the Schrödinger equation, $`H\psi (x,y)`$=$`E\psi (x,y)`$, the solution is a standing wave in the $`x`$ axis and sinusoidal wave of squared $`y`$ in the $`y`$ axis. The potential gradient constant $`F`$ is obtained as $`5\times 10^3`$ eV/nm from the LDOS of eigen wave functions with the Hamiltonian around the $`K`$\- and $`K^{}`$\- points, where the variation of the potential over 200 nm along the $`y`$ axis is about 1 eV. This model satisfies conditions of the experimental results.

## 3 Electronic wave interference effec
in ultra thin graphene layers

Superperiodic patterns that extended over several $`\mu `$m2 were observed by STM on ultra thin graphite layers. Figure 2(a) shows a part of those patterns in an 500$`\times `$500 nm2 image near a step edge observed at a sample bias voltage of 0.3 V. In this image, the patterns are sharp-shaped and the periodicity is very large ($``$100 nm). We consider that the patterns come from the dislocation-network structures due to the similarity of the patterns to those reported by STM and transmission electron microscopy observations.ref6 ; ref7  Cross-sectional analyses of STM images indicate that apparent heights of the lower and upper terraces in fig.2 correspond to two and three graphene layers from the substrate, respectively. Figure 2(b) shows an image that is shifted to bottom direction of fig.2(a) by 200 nm close to the Fermi energy ($`V_s`$=0.02 V). In this image, we can classify patterns into three shapes. Triangular-, rhombic- and net-shaped patterns were observed in domains A, B and C, respectively. The geometric patterns in domains B and C appears to have an inversed contrast though those patterns come from an array of dislocations as is evident from the fact that the patterns are continuous at the step edge. The corrugation amplitude of the patterns in domain A and B in fig.2(b) is larger than that in fig.2(a). In contrast, the corrugation amplitude of the patterns in domain C in fig.2(b) is smaller than that in fig.2(a). These phenomena cannot be explained just by the LDOS of stacking faulted parts and the tunneling probability dependent on a gap between a tip and a sample. By increasing the bias voltage, the patterns on the upper terrace in fig.2(b) changed into the net-shaped pattern although the pattern shape on the lower terrace does not change, as shown in fig.2(c). This phenomenon cannot be explained by a change of the dislocation-network structures driven by the applied bias voltages because the periodicity did not change. We can consistently explain those phenomena as the interference effect originating from the perturbed and unperturbed waves in the overlayer, by assuming the scattering potential at the interface to generate the patterns.
Here, we theoretically treat the interference effect for explaining the bias dependence of superperiodic patterns by a free electron model assuming an effective mass for electrons.ref9 ; ref11  An in-plane square-patterned potential with a periodicity of 2$`L`$ at the interface is employed for a calculation of the LDOS of the wave function confined in the plane to generate an abrupt potential change associated with the dislocation-network structures. Although the square-patterned potential is not consistent with the experimentally observed triangular-shaped pattern, the model makes the analysis considerably simple without any loss of validity. We place a square potential with $`L`$/3 in width and 4$`v_0\delta (z)`$ in height at the lines dividing the pattern into individual geometric unit as a simple model for reproducing the patterns in domain B and C. By the Fourier analysis, this potential is expressed as,

$`V(x,y,z)=(\mathrm{}^2/m_{})v_0{\displaystyle \underset{n}{}}a_n\delta (z)`$
$`\times (e^{iq_{xn}x}+e^{iq_{xn}x}+e^{iq_{yn}y}+e^{iq_{yn}y}),`$ (3)
where $`\mathrm{}`$ is the Planck constant over 2$`\pi `$, $`m_{}`$ is the effective mass along the $`z`$ axis, $`a_n`$ is the $`n`$th component which equals to $`\{2/(n\pi )\}\{\mathrm{sin}(n\pi )\mathrm{sin}(5n\pi /6)\}`$ for the square potential, and $`q_{xn}`$ and $`q_{yn}`$, which take discrete values $`(n\pi /L)`$($`n`$=1,2,…), are the $`n`$th wave vectors along the $`x`$ and $`y`$ axes, respectively. We can find a wave function in the form of

$$\mathrm{\Psi }(x,y,z)=\underset{q_x,q_y}{}A_{q_x,q_y}(z)e^{i(q_xx+q_yy)},$$
(4)
where $`A_q(z)`$ is a wave function for the $`z`$ component. The perturbed wave is obtained as a solution of the Schrödinger equation with the discontinuity of $`d^2/dz^2`$ at the interface. Provided a perturbed wave number and a unperturbed wave number are $`k^{}`$ and $`k`$, respectively, the LDOS of a superperiodic pattern at the surface is roughly propotional to the sum of $`[a_nk\mathrm{sin}(kl)\mathrm{cos}(k^{}l)\{\mathrm{cos}(q_{xn}x)+\mathrm{cos}(q_{yn}y)\}]`$ with index $`n`$, where $`l`$ is the overlayer height from the interface. The effective mass of electrons was estimated from the band dispersion of graphite. The calculated LDOS suggested that the corrugation amplitude of the superperiodic pattern could increase, decrease or invert, dependent on the overlayer thickness and the bias voltage, and consistently reproduced the experimental results.

## 4 Summary

A superperiodic pattern, whose periodicity changes in plane, was observed by STM on nanographene sheets that inclined from the HOPG substrate with weak interaction from the substrate. The periodicity changed remarkably by attachment of a nanographite flake on the nanographene sheets. The pattern can be explained as the interference effect of an electronic wave in the nanographene sheets on the basis of the $`kp`$ perturbation model that assumes a linear potential in the direction of slope and the electron confinement effect.
Superperiodic patterns, which come from the dislocation-network structure, were also observed by STM. Individual shapes of the observed patterns were different at terraces with different number of layers from the substrate and their contrasts changed depending on the bias voltage. The patterns can be explained in terms of the interference effect by perturbed and unperturbed waves in the ultra thin graphite overlayer.

The authors are grateful to Prof. Katsuyoshi Kobayashi (Ochanomizu University), for fruitful discussion. They also thank to Dr. A. Moore for his generous supply of HOPG sample. The present work was supported partly by the Grant-in-Aid for ‘Research for the future’ Program, Nano-Carbon, and 15105005 from Japan Society for the promotion of Science. One of the authors (K.H.) acknowledges the financial support from NEDO via Synthetic Nano Functional Materials Project, AIST, Japan.