text:
# Triangular arbitrage as an interaction among foreign exchange rates

## 1 Introduction

The triangular arbitrage is a financial activity that takes advantage of the three exchange rates among three currencies. Suppose that we exchange one US dollar to some amount of Japanese yen, exchange the amount of Japanese yen to some amount of euro, and finally exchange the amount of euro back to US dollar; then how much US dollar do we have? There are opportunities that we have more than one US dollar. The triangular arbitrage transaction is the trade that takes this type of opportunities. It has been argued that the triangular arbitrage makes the product of the three exchange rates converge to a certain value . In other words, the triangular arbitrage is a form of interaction among currencies.

The purpose of this paper is to show that there is in fact triangular arbitrage opportunities in foreign exchange markets and they generate an interaction among foreign exchange rates. We first analyze real data in Sec. 2, showing that the product of three foreign exchange rates has a narrow distribution with fat tails. In order to explain the behavior, we propose in Sec. 3 a model of the time evolution of exchange rates with an interaction. The simulation explain the real data well.

## 2 Observation of the Triangular Arbitrage Opportunity in High-Frequency Data

In the present paper, we analyze actual data of the yen-dollar rate, the yen-euro rate and the dollar-euro rate, taken from January 25 1999 to March 12 1999 except for weekends. We show in this section that there are actually triangular arbitrage opportunities and that the three exchange rates correlate strongly.

### 2.1 Existence of Triangular Arbitrage Opportunities

In order to quantify the triangular arbitrage opportunities, we define the quantity

$$\mu (t)=\underset{i=1}{\overset{3}{}}r_i(t),$$
(1)
where $`r_i(t)`$ denotes each exchange rate at time $`t`$. We refer to this quantity as the rate product. There is a triangular arbitrage opportunity whenever the rate product is greater than unity.

To be more precise, there are two types of the rate product. One is based on the arbitrage transaction in the direction of dollar to yen to euro to dollar. The other is based on the transaction in the opposite direction of dollar to euro to yen to dollar. Since these two values show similar behavior, we focus on the first type of $`\mu (t)`$ hereafter. Thus, we specifically define each exchange rate as

$`r_1(t)`$ $``$ $`{\displaystyle \frac{1}{\text{yen-dollar ask }(t)}}`$ (2)
$`r_2(t)`$ $``$ $`{\displaystyle \frac{1}{\text{dollar-euro ask }(t)}}`$ (3)
$`r_3(t)`$ $``$ $`\text{yen-euro bid }(t).`$ (4)
(Note the difference between

$$\begin{array}{ccccc}1& [dollar]& & \text{yen-dollar bid}& [yen]\end{array}$$
(5)
and

$$\begin{array}{ccccc}1& [yen]& & \frac{1}{\text{yen-dollar ask}}& [dollar].\text{)}\end{array}$$
(6)
We assume here that an arbitrager can transact instantly at the bid and the ask prices provided by information companies and hence we use the prices at the same time to calculate the rate product.

Figure 1(a)-(c) shows the actual changes of the three rates: the yen-euro ask, the dollar-euro ask and the yen-euro bid. Figure 1(d) shows the behavior of the rate product $`\mu (t)`$. We can see that the rate product $`\mu `$ fluctuates around the average

$$m\mu (t)0.99998.$$
(7)
(The average is less than unity because of the spread; the spread is the difference between the ask and the bid prices and is usually of the order of $`0.05\%`$ of the prices.) The probability density function of the rate product $`\mu `$ (Fig. 2) has a sharp peak and fat tails while those of the three rates (Fig. 3) do not. It means that the fluctuations of the exchange rates have correlation that makes the rate product converge to the average $`m`$.

### 2.2 Feasibility of the Triangular Arbitrage Transaction

We discuss here the feasibility of the triangular arbitrage transaction. We analyze the duration of the triangular arbitrage opportunities and calculate whether an arbitrager can make profit or not.

The shaded area in Fig. 2 represents triangular arbitrage opportunities. We can see that the rate product is grater than unity for about $`6.4\%`$ of the time. It means that triangular arbitrage opportunities exist about ninety minutes a day. The ninety minutes, however, include the cases where the rate product $`\mu `$ is greater than unity very briefly. The triangular arbitrage transaction is not feasible in these cases.

In order to quantify the feasibility, we analyze the duration of the triangular arbitrage opportunities. Figure 4 shows the cumulative distributions of the duration $`\tau _+`$ of the situation $`\mu >1`$ and $`\tau _{}`$ of $`\mu <1`$. It is interesting that the distribution of $`\tau _+`$ shows a power-law behavior while the distribution of $`\tau _{}`$ dose not. This difference may suggest that the triangular arbitrage transaction is carried out indeed.

In order to confirm the feasibility of the triangular arbitrage, we simulate the triangular arbitrage transaction using our time series data. We assume that it takes an arbitrager $`T_{\mathrm{rec}}`$\[sec\] to recognize triangular arbitrage opportunities and $`T_{\mathrm{exe}}`$\[sec\] to execute a triangular arbitrage transaction; see Fig. 5. We also assume that the arbitrager transacts whenever the arbitrager recognizes the opportunities. Figure 6 shows how much profit the arbitrager can make from one US dollar (or Japanese yen or euro) in a day. We can see that the arbitrager can make profit if it takes the arbitrager a few seconds to recognize the triangular arbitrage opportunities and to execute the triangular arbitrage transaction.

## 3 Modeling

We here introduce a new model that takes account of the effect of the triangular arbitrage transaction as an interaction among the three rates. Many models of price change have been introduced so far: for example, the L$`\stackrel{´}{\mathrm{e}}`$vy-stable non-Gaussian model ; the truncated L$`\stackrel{´}{\mathrm{e}}`$vy flight ; the ARCH/GARCH processes . They discuss, however, only the change of one price. They did not consider an interaction among multiple prices. As we discussed in Sec. 2, however, the triangular arbitrage opportunity exists in the market and is presumed to affect price fluctuations in the way the rate product tends to converge to a certain value.

### 3.1 Basic Time Evolution

The basic equation of our model is a time-evolution equation of the logarithm of each rate:

$$\mathrm{ln}r_i(t+\mathrm{\Delta }t)=\mathrm{ln}r_i(t)+f_i(t)+g(\nu (t)),(i=1,2,3)$$
(8)
where $`\nu `$ is the logarithm of the rate product

$$\nu (t)\mathrm{ln}\mu (t)=\underset{i=1}{\overset{3}{}}\mathrm{ln}r_i(t).$$
(9)
Just as $`\mu `$ fluctuates around $`m=\mu 0.99998`$, the logarithm rate product $`\nu `$ fluctuates around

$$ϵ\mathrm{ln}\mu 0.00091$$
(10)
(Fig. 7(a)). In this model, we focus on the logarithm of the rate-change ratio $`\mathrm{ln}(r_i(t+\mathrm{\Delta }t)/r_i(t))`$, because the relative change is presumably more essential than the absolute change. We assumed in Eq. (8) that the change of the logarithm of each rate is given by an independent fluctuation $`f_i(t)`$ and an attractive interaction $`g(\nu )`$. The triangular arbitrage is presumed to make the logarithm rate product $`\nu `$ converge to the average $`ϵ`$; thus, the interaction function $`g(\nu )`$ should be negative for $`\nu `$ greater than $`ϵ`$ and positive for $`\nu `$ less than $`ϵ`$:

$`g(\nu )\{\begin{array}{ccc}<0& ,& \text{for }\nu >ϵ\\ >0& ,& \text{for }\nu <ϵ.\end{array}`$ (13)
As a linear approximation, we define $`g(\nu )`$ as

$$g(\nu )a(\nu ϵ)$$
(14)
where $`a`$ is a positive constant which specifies the interaction strength.

The time-evolution equation of $`\nu `$ is given by summing Eq. (8) over all $`i`$:

$$\nu (t+\mathrm{\Delta }t)ϵ=(13a)(\nu (t)ϵ)+F(t),$$
(15)
where

$$F(t)\underset{i=1}{\overset{3}{}}f_i(t).$$
(16)
This is our basic time-evolution equation of the logarithm rate product.

### 3.2 Estimation of Parameters

The interaction strength $`a`$ is related to the auto-correlation function of $`\nu `$ as follows:

$$13a=c(\mathrm{\Delta }t)\frac{\nu (t+\mathrm{\Delta }t)\nu (t)\nu (t)^2}{\nu ^2(t)\nu (t)^2}.$$
(17)
Using Eq. (17), we can estimate $`a(\mathrm{\Delta }t)`$ from the real data series as a function of the time step $`\mathrm{\Delta }t`$. The auto-correlation function $`c(\mathrm{\Delta }t)`$ is shown in Fig. 8(a). The estimate of $`a(\mathrm{\Delta }t)`$ is shown in Fig. 8(b). Hereafter, we fix the time step at $`\mathrm{\Delta }t=60`$\[sec\] and hence use

$$a(1[\mathrm{min}])=0.17\pm 0.02$$
(18)
for our simulation.

On the other hand, the fluctuation of foreign exchange rates is known to be a fat-tail noise . Here we take $`f_i(t)`$ as the truncated L$`\stackrel{´}{\mathrm{e}}`$vy process :

$`P_{\mathrm{TLF}}(f;\alpha ,\gamma ,l)=qP_\mathrm{L}(f;\alpha ,\gamma )\mathrm{\Theta }(l|f|),`$ (19)
where $`q`$ is the normalization constant, $`\mathrm{\Theta }(x)`$ represents the step function and $`P_\mathrm{L}(x;\alpha ,\gamma )`$ is the symmetric L$`\stackrel{´}{\mathrm{e}}`$vy distribution of index $`\alpha `$ and scale factor $`\gamma `$:

$$P_\mathrm{L}(x;\alpha ,\gamma )=\frac{1}{\pi }_0^{\mathrm{}}e^{\gamma |k|^\alpha }\mathrm{cos}(kx)dk0<\alpha <2.$$
(21)
We determine the parameters $`\alpha `$, $`\gamma `$ and $`l`$ by using the following relations for $`1<\alpha <2`$ :

$`c_2`$ $`=`$ $`{\displaystyle \frac{\alpha (\alpha 1)\gamma }{|\mathrm{cos}(\pi \alpha /2)|}}l^{2\alpha },`$ (22)
$`\kappa `$ $`=`$ $`{\displaystyle \frac{(3\alpha )(2\alpha )|\mathrm{cos}(\pi \alpha /2)|}{\alpha (\alpha 1)\gamma }}l^\alpha ,`$ (23)
where $`c_n`$ denotes the $`n`$th cumulant and $`\kappa `$ is the kurtosis $`\kappa =c_4/c_{2}^{}{}_{}{}^{2}`$. The estimates are shown in Table 1.

The generated noises with the estimated parameters are compared to the actual data in Fig. 9.

We simulated the time evolution (15) with the parameters given in Eqs. (10), (18) and Table 1. The probability density function of the results (Fig. 7(b)) is compared to that of the real data (Fig. 7(a)) with $`\mathrm{\Delta }t=1`$\[min\] in Fig. 10. The fluctuation of the simulation data is consistent with that of the real data. In particular, we see good agreement around $`\nu ϵ`$ as a result of the linear approximation of the interaction function. Figure 7(c) shows $`\nu (t)`$ of the simulation without the interaction, i.e. $`a=0`$. The quantity $`\nu `$ fluctuates freely, which is inconsistent with the real data.

### 3.3 Analytical approach

We can solve the time-evolution equation (15) analytically in some cases. Let us define

$$\omega (t)\nu (t)ϵ$$
(24)
and

$$\{\begin{array}{cccccc}\omega \hfill & =\hfill & \omega (t),\hfill & \omega ^{}\hfill & =\hfill & \omega (t+\mathrm{\Delta }t),\hfill \\ F\hfill & =\hfill & F(t),\hfill & F^{}\hfill & =\hfill & F(t+\mathrm{\Delta }t),\hfill \\ c\hfill & =\hfill & c(\mathrm{\Delta }t).\hfill & & & \end{array}$$
(25)
Equation (15) is then reduced to

$$\omega ^{}=c\omega +F$$
(26)
Assume that the probability of $`\omega `$ having a value in $`\omega \omega +\mathrm{d}\omega `$ is $`P_\omega (\omega )`$. The joint probability of $`\omega ^{}`$ having a value in $`\omega ^{}\omega ^{}+\mathrm{d}\omega ^{}`$ and $`F^{}`$ having a value in $`F^{}F^{}+\mathrm{d}F^{}`$ is given by

$`P_{\omega ^{},F^{}}(\omega ^{},F^{})\mathrm{d}\omega ^{}\mathrm{d}F^{}`$ $`=`$ $`P_\omega (\omega )P_F(F^{})\mathrm{d}\omega \mathrm{d}F^{}`$ (27)
$`=`$ $`{\displaystyle \frac{1}{c}}P_\omega ({\displaystyle \frac{\omega ^{}F^{}}{c}})P_F(F^{})\mathrm{d}\omega ^{}\mathrm{d}F^{},`$
where $`P_F(F)`$ is the probability of $`F`$ having a value in $`FF+\mathrm{d}F`$. The probability density function of $`\omega ^{}`$ is thus given by

$$P_\omega ^{}(\omega ^{})=\frac{1}{c}P_\omega (\frac{\omega ^{}F^{}}{c})P_F(F^{})dF^{}.$$
(28)
The characteristic function of $`\omega ^{}`$ is the Fourier transform

$$\stackrel{~}{P}_\omega ^{}(k)=P_\omega ^{}(\omega ^{})e^{i\omega ^{}k}d\omega ^{}=\stackrel{~}{P}_\omega (ck)\stackrel{~}{P}_F(k),$$
(29)
where $`\stackrel{~}{P}_\omega `$ and $`\stackrel{~}{P}_F`$ are the Fourier transforms of $`P_\omega `$ and $`P_F`$, respectively. Then we obtain

$`\stackrel{~}{P}_{\omega (t)}(k)`$ $`=`$ $`\stackrel{~}{P}_{\omega (t\mathrm{\Delta }t)}(ck)\stackrel{~}{P}_{F(t\mathrm{\Delta }t)}(k)`$ (30)
$`=`$ $`\stackrel{~}{P}_{\omega (0)}(c^Nk){\displaystyle \underset{n=0}{\overset{N}{}}}\stackrel{~}{P}_F(c^nk)`$
$`=`$ $`\stackrel{~}{P}_{\omega (0)}(c^Nk){\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle \underset{n=0}{\overset{N}{}}}\stackrel{~}{P}_{f_i}(c^nk)`$
$`=`$ $`{\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle \underset{n=0}{\overset{N}{}}}\stackrel{~}{P}_{f_i}(c^nk),`$
where $`t=N\mathrm{\Delta }t`$. We here assumed $`P_{\omega (0)}(\omega )=\delta (\omega )`$ and hence $`\stackrel{~}{P}_{\omega (0)}(k)=1`$.

The above argument shows the essential reason of the sharp peak and fat tails in Fig. 2. If we had $`c=1`$, or $`a=0`$ (without the interaction), the noise $`F`$ at every time step would accumulate in $`\omega `$ and the probability density function of $`\omega =\nu ϵ`$ would be Gaussian due to the central limit theorem. If we have $`c<1`$, or $`a>0`$ (with the interaction), the noise at the past time steps decay as $`c^n`$. The largest contribution to $`\omega `$ comes from the noise one time step before, which is a fat-tail noise .

As a special case, if the noises $`f_i(t)`$ obey a L$`\stackrel{´}{\mathrm{e}}`$vy distribution of the same index $`\alpha `$ and the same scale factor $`\gamma `$, namely if

$$\stackrel{~}{P}_{f_i}(k)=e^{\gamma |k|^\alpha }\text{for all}i,$$
(31)
the distribution of $`\omega `$ is also a L$`\stackrel{´}{\mathrm{e}}`$vy distribution of the same index $`\alpha `$ and a different scale factor $`\gamma ^{}`$ given by

$$\gamma ^{}=\frac{3}{1\{c(\mathrm{\Delta }t)\}^\alpha }\gamma .$$
(32)
## 4 Summary

We first showed that triangular arbitrage opportunities exist in the foreign exchange market. The rate product $`\mu `$ fluctuates around a value $`m`$. Next, we introduced a model including the interaction caused by the triangular arbitrage transaction. We showed that the interaction is the reason of the sharp peak and the fat-tail property of the distribution of the logarithm rate product $`\nu `$. Finally we showed that our model is solvable analytically in some cases.

## Acknowledgments

One of the authors (Y. A.) is grateful to Mr. Tetsuro Murai and Mr. Takashi Anjo for their interesting comments.


format: text

text:
# Interplay Between Time-Temperature-Transformation and the Liquid-Liquid Phase Transition in Water

## Abstract

We study the TIP5P water model proposed by Mahoney and Jorgensen, which is closer to real water than previously-proposed classical pairwise additive potentials. We simulate the model in a wide range of deeply supercooled states and find (i) the existence of a non-monotonic “nose-shaped” temperature of maximum density line and a non-reentrant spinodal, (ii) the presence of a low temperature phase transition, (iii) the free evolution of bulk water to ice, and (iv) the time-temperature-transformation curves at different densities.

preprint: M. Yamada et al.

Much effort has been invested in exploring the overall phase diagram of water and the connection among its liquid, supercooled and glassy states Debenedetti96 ; Smith99 ; Mishima94 , with particular interest in understanding the origin of the striking anomalies at low temperatures, such as the $`T`$-dependence of the isothermal compressibility $`K_T`$, the constant pressure specific heat $`C_P`$, and the thermal expansivity $`\alpha _P`$.

The “stability limit conjecture” attributes the increase of the response functions upon supercooling to a continuous re-tracing spinodal line bounding the superheated, supercooled and stretched (negative pressure) metastable states Speedy82 . This line at its minimum intersects the temperature of maximum density (TMD) curve tangentially. More recently, a different hypothesis has been developed, for which the spinodal does not re-enter into the positive pressure region, but rather the anomalies are attributed to a critical point below the homogeneous nucleation line Poole92 . The TMD line, which is negatively sloped at positive pressures, becomes positively sloped at sufficiently negative pressures and does not intersect the spinodal. A line of first order phase transitions — interpreted as the liquid state analog of the line separating low and high density amorphous glassy phases Poole92 ; Mishima94  — develops from this critical point.

Simulations of supercooled metastable states are possible because the structural relaxation time at the temperatures of interest is several orders of magnitude shorter than the crystallization time. It is difficult, but not impossible kusalik95 , to observe crystallization in simulations of molecular models Huitema00  because homogeneous nucleation rarely occurs on the time scales reachable by present day computers. Bulk water simulations have been crystallized by applying a homogeneous electric field Svishchev94  or placing liquid water in contact with pre-existing ice Borzsak99 ; Clancy95 , but spontaneous crystallization of deeply supercooled model water has not been observed in simulations.

In contrast, experimental measurements of metastable liquid states are strongly affected by homogeneous nucleation. The nucleation and growth of ice particles from aqueous solution has been extensively studied, and the “nose-shaped” time-temperature-transformation (TTT) curves have been measured Debenedetti96 ; MacFarlane83 ; Kresin91 . The non-monotonic relation between crystallization rate and supercooling depth results from the competition between the thermodynamic driving force for nucleation and the kinetics of growth Debenedetti96 . Crystallization hinders direct experimental investigation of pure metastable liquid water below the homogeneous nucleation line; only indirect measurements can be made by studying the metastable melting lines of ices Mishima94 .

This work attempts to unify the phenomena connected with the existence of a liquid-liquid phase transition and homogeneous nucleation in a single molecular dynamics simulation study. We simulate a system of $`N=343`$ molecules interacting with the TIP5P potential Mahoney00 . TIP5P is a five-site, rigid, non-polarizable water model, not unlike the ST2 model Stillinger74 . The TIP5P potential accurately reproduces the density anomaly at 1 atm and exhibits excellent structural properties when compared with experimental data Mahoney00 ; Head-Gordon . The TMD shows the correct pressure dependence, shifting to lower temperatures as pressure is increased. Under ambient conditions, the diffusion constant is close to the experimental value, with reasonable temperature and pressure dependence away from ambient conditions Mahoney00 .

We perform equilibration runs at constant $`T`$ (Berendsen’s thermostat), while we perform production runs in the microcanonical (NVE) ensemble. After thermalization at $`T=320`$ K we set the thermostat temperature to the temperature of interest. We let the system evolve for a time longer than the structural relaxation time $`\tau _\alpha `$, defined as the time at which $`F_s(Q_0,\tau _\alpha )=1/e`$, where $`F_s(Q_0,t)`$ is the self-intermediate scattering function evaluated at $`Q_0=18`$ nm-1, the location of the first peak of the static structure factor. In the time $`\tau _\alpha `$, each molecule diffuses on average a distance of the order of the nearest neighbor distance. We use the final configuration of the equilibration run to start a production run of length greater than several $`\tau _\alpha `$ and then analyze the calculated trajectory. We check that no drift in any of the studied quantities and no crystallization occurs during the production run.

In Fig. 1 we show results for pressure along isotherms. At lower temperatures an inflection develops, which becomes a “flat” isotherm at the lowest temperature, $`T=215`$ K. The presence of a flat region indicates that a phase separation takes place, and we estimate the critical temperature $`T_C^{}=217\pm 3`$ K, the critical pressure $`P_C^{}=340\pm `$ 20 MPa, and the critical density $`\rho _C^{}=1.13\pm 0.04`$ g/cm3.

In Fig. 2(a) we plot the pressure along isochores. The curves show minima as a function of temperature; the locus of the minima is the TMD line, since $`(P/T)_V=\alpha _P/K_T`$. It can be seen that the pressure exhibits a minimum if the density passes through a maximum ($`\alpha =0`$). It is clear that, as in the case of ST2 water, TIP5P water has a TMD that changes slope from negative to positive as P decreases. Notably, the point of crossover between the two behaviors is located at ambient pressure, $`T4`$ C , and $`\rho 1`$ g/cm3.

We also plot the spinodal line. We calculate the points on the spinodal line fitting the isotherms (for $`T300K`$) of Fig. 1 to the form $`P(T,\rho )=P_s(T)+A\left[\rho \rho _s(T)\right]^2`$, where $`P_s(T)`$ and $`\rho _s(T)`$ denote the pressure and density of the spinodal line. This functional form is the mean field prediction for $`P(\rho )`$ close to a spinodal line. For $`T250K`$, we calculate $`P_s(T)`$ by estimating the location of the minimum of $`P(\rho )`$. The results in Fig.2 show that the liquid spinodal line is not reentrant and does not intersect the TMD line highT .

A supercooled liquid is metastable with respect to the crystal, so it is driven to crystallize Debenedetti96 . However, crystallization of model water has not been found in simulations because the homogeneous nucleation time far exceeds the CPU time. However, for TIP5P, crystallization times lie within a time window accessible to present-day simulations, and we observe crystallization at densities $`\rho =1.15`$ g/cm3, and $`1.20`$ g/cm3 for a wide range of temperatures (Fig. 3).

To quantify the crystallization process, we analyze four independent configurations thermalized at temperature $`T=320`$ K and instantaneously quenched to the temperature of interest. We monitor the potential energy as well as the time evolution of the structure factor $`S_\stackrel{}{Q}(t)\rho _\stackrel{}{Q}(t)\rho _\stackrel{}{Q}^{}(t)/N`$ at all wave vector $`\stackrel{}{Q}`$ values, ranging from the smallest value $`2\pi /L`$ allowed by the side $`L`$ of the simulation box up to $`50`$ nm-1. The oxygen density fluctuation $`\rho _\stackrel{}{Q}(t)_{i=1}^N\mathrm{exp}(i\stackrel{}{Q}\stackrel{}{r}_i)`$, where $`\stackrel{}{r}_i`$ is the oxygen coordinate of molecule $`i`$. The onset of crystallization coincides with the occurrence of (i) a sudden drop in potential energy and (ii) a sharp increase in the density fluctuations at one or more wave vector values. When crystallization occurs, the value of $`S_\stackrel{}{Q}(t)`$ jumps from $`𝒪`$(1) in the liquid to $`𝒪`$(N).

Defining the crystallization time is somewhat arbitrary because of the stochasticity which accompanies the onset of crystallization and the definition of the critical nucleus. We define $`\tau _{cryst}`$ as the time at which any density fluctuation $`S_\stackrel{}{Q}(t)`$ grows above a threshold value $`S^{}=15`$ and remains continuously above the threshold for a time exceeding $`t^{}=40`$ ps. This threshold prevents transient density fluctuations from being attributed to crystallization. We also perform calculations for other definitions of $`S^{}`$ and $`t^{}`$, but the above values are sufficient to unambiguously identify the onset of crystallization without requiring excessive simulation time. Fig. 3 shows the crystallization times $`\tau _{cryst}`$, averaged over the four independent runs, for two different densities and for a broad range of $`T`$. The resulting TTT curve shows a characteristic “nose” shape, arising from the competition between two effects, the thermodynamic driving force for nucleation and the kinetics of growth Debenedetti96 ; kusalik95 . As temperature is lowered, both the thermodynamic driving force and the relaxation time increase, and it becomes more difficult for particles to diffuse to the energetically-preferred crystalline configuration. For both densities, $`\rho =1.15`$, and $`1.20`$ g/cm3, the $`T`$ at which nucleation is fastest is around $`240`$ K. At this $`T`$, the onset of crystallization requires about $`3`$ ns. At the lowest studied $`T`$, the crystallization time has grown to $`30`$ ns.

Fig. 3 also shows the relaxation times $`\tau _\alpha `$. The $`T`$ dependence of $`\tau _\alpha `$ can be described by a power law $`\tau _\alpha (TT_c)^\gamma `$, in agreement with the prediction of mode coupling theory leshouches . Since the relation $`\tau _\alpha \tau _{cryst}`$ holds at each temperature, including in the deeply supercooled region, “equilibrium” studies of metastable water can be achieved before nucleation takes place. The liquid can be connected to the deeply supercooled state via “equilibrium” metastable states if we choose a quench rate larger than the critical cooling rate $`_c(T_mT_n)/\tau _n`$, where $`T_m`$ is the melting $`T`$ and $`T_n`$ and $`\tau _n`$ locate the nose in the TTT curve Debenedetti96 . For TIP5P, $`_c10^{10}`$ K/s at the two studied densities. For $`\rho =1.10`$ g/cm3 and $`T=240`$ K, we observe only one (out of four) crystallization event within a time of 70 ns. For densities smaller than $`\rho =1.10`$ g/cm3, we observe no crystallization events within a time of 60 ns and hence we can only estimate that $`_c`$ is smaller than $`10^9`$ K/s (the experimental value for water at ambient pressure is $`_c10^7`$ K/s uhlmann72 ).

In Fig. 4 we show a typical crystal configuration. The crystal structure, after energy minimization at constant volume, is a proton-ordered structure similar to ice-B, first observed by Baez and ClancyClancy95 . Ice-B is a variant of the ice IX structure, which is the proton-ordered form of ice III. The density of ice IX and ice III is in fact $`1.16`$ g/cm3, close to our value.

We have shown that the liquid-liquid phase separation can be observed in metastable equilibrium (i) if the cooling rate is faster than $`_c`$, and (ii) if the observation time is shorter than the crystallization time at the critical point. While both such conditions can be realized in numerical simulations — as shown here — they cannot be met in experiments. Our simulations also show that a continuity of states between liquid and glassy phases of water exists Hallbrucker89 . Liquid states below the crystallization temperature can be accessed provided the cooling rate exceeds $`_c`$.

We acknowledge useful discussions with D. R. Baker, G. Franzese, W. Kob, E. La Nave, M. Marquez and C. Rebbi. We acknowledge support from the NSF Grant CHE-0096892. M. Y. acknowledges support from NSF Grant GER-9452651 as a Graduate Research Trainee at the Boston University Center for Computational Science. F. S. acknowledges support from MURST COFIN 2000 and INFM Iniziativa Calcolo Parallelo.


text:
# X-ray Synchrotron Spectral Hardenings from Compton and Synchrotron Losses in Extended Chandra Jets

## 1 Introduction

The excellent spatial resolution of the Chandra X-ray Observatory is providing X-ray images of extended radio galaxy jets in the 0.2-8 keV band at resolutions better than $`0.4^{\prime \prime }`$ FWHM, and spectral detail of discrete components at flux levels between $`10^{13}`$ and $`10^{14}`$ ergs cm-2 s-1. Combined with high resolution radio and optical maps, spectral energy distributions (SEDs) of spatially-resolved hot spots and knots in the jets of radio galaxies can be studied from the radio through the X-ray regime.

The broadband SEDs in many of the hot spots and knots show a general behavior whereby the X-ray spectra are harder than the optical spectra, and the X-ray fluxes are above an extrapolation from the optical fluxes. This behavior is observed in the SEDs of at least three of the four regions along the X-ray jet of 3C 273 (Sambruna et al., 2001; Marshall et al., 2001), the western hot spot of Pictor A (Wilson, Young, and Shopbell, 2001), knot WK7.8 of PKS 0637-752 (Schwartz et al., 2000; Chartas et al., 2000), hot spot A of Cygnus A (Wilson, Young, and Shopbell, 2000), and in some of the knots in the jet of M87 (Wilson & Yang, 2001). In other cases, such as the X-ray jet in 3C 66B (Hardcastle, Birkinshaw, & Worrall, 2001), knots A1 and A3 in 3C 273 (Marshall et al. (2001); see, however, Sambruna et al. (2001) for a different spectral analysis), and some other knots in M87’s jet (Wilson & Yang, 2001), the X-ray–optical–radio spectrum is consistent with a single or smoothly broken power-law spectrum, indicating that the broadband emission is entirely due to nonthermal synchrotron radiation.

Three main leptonic processes have been considered to account for the X-ray emission in the extended jets (e.g., Harris and Krawczynski, 2002): nonthermal synchrotron radiation, synchrotron self-Compton (SSC) radiation, and X-ray emission from Compton-scattered external radiation fields. Nonthermal synchrotron emission is favored for smooth spectra and for the western hot spot of Pictor A (Wilson, Young, and Shopbell, 2001), where an SSC origin implies a very weak magnetic field, and a combined SSC/synchrotron origin requires that the two emission components be carefully matched in the X-ray regime. The SSC process can, however, explain the X-ray emission mechanism from the bright radio hot spots of Cygnus A (Wilson, Young, and Shopbell, 2000), the radio hot spots of 3C 295 (Harris et al., 2000), and the eastern hot spot of 3C 123 (Hardcastle, Birkinshaw, & Worrall, 2001a).

Tavecchio et al. (2000) and Celotti, Ghisellini, & Chiaberge (2001) argue that the X-ray emission from knot WK7.8 is due to Compton-scattered cosmic microwave background radiation (CMBR) rather than to SSC emission, because an SSC model requires a system well out of equilibrium and a jet that is significantly debeamed. This model invokes an emitting region in bulk relativistic motion on size scales of several hundred kpc, in accord with observations of one-sidedness in the jets on these size scales. Such a model has been applied to X-ray emitting knots in the jet of 3C 273 by Sambruna et al. (2001). A synchrotron model for the X-ray emission is discounted because it would apparently require two populations of relativistic electrons to explain the radio, optical, and X-ray fluxes. Another possibility is that the X-rays are synchrotron radiation of high-energy protons (Aharonian, 2001).

Here we show that synchrotron emission from a single population of electrons injected with a nonthermal power law spectrum can produce X-ray spectral hardenings. When CMBR cooling in the Thomson regime exceeds synchrotron cooling, a hardening in the electron spectrum is formed at electron energies where KN effects become important. This produces a hardening in the synchrotron spectrum between optical and X-ray frequencies. In the next section, an analytic model for the synchrotron origin of X-ray emission from extended jets is presented, motivated by Chandra observations of 3C 273. Results of a numerical simulation with cascade reprocessing are given in Section 3. Section 4 gives a discussion and summary.

## 2 Analytic Model

In 3C 273, the X-ray energy spectral index in the inner knot region A has $`\alpha _x1.1`$, a value characteristic of a cooling shock-accelerated electron spectrum, and there is an indication of spectral softening along the jet (Sambruna et al., 2001). The peaks of the X-ray emission along the jet of 3C 273 closely track the peaks in the optical waveband, and the overall X-ray emission profile displays a decreasing mean intensity along the jet, whereas the optical intensity is roughly uniform along the jet. The radio profile also displays brightenings in spatial coincidence with the X-ray and optical peaks, but the radio intensity increases dramatically along the jet. The optical emission profiles are narrower than those of the radio emission.

These behaviors suggest that the electron population producing the radio emission accumulates due to weak cooling, whereas the electrons producing the optical and X-ray emission strongly cool. Electrons that produce synchrotron radio emission at observer’s frame frequencies $`\nu `$ have Lorentz factors $`\gamma 5\times 10^4[(\nu /5\mathrm{GHz})(1+z)/\delta B_{\mu \mathrm{G}}]^{1/2}`$, whereas the electrons that Compton scatter the CMBR to X-ray energies $`E_x`$ have $`\gamma 1\times 10^3[(E_x/1\mathrm{keV})/\delta \mathrm{\Gamma }]^{1/2}`$. Consequently the X-ray emitting electrons will cool more slowly than the radio-emitting electrons unless $`B`$ exceeds milli-Gauss intensities. The fast-cooling requirement for the X-ray emitting electrons is readily satisfied assuming a synchrotron origin for the X-rays. In this paper, we neglect complications that involve variations of the magnetic field $`B`$ and the Doppler factor $`\delta =[\mathrm{\Gamma }(1\beta _\mathrm{\Gamma }\mathrm{cos}\theta )]^1`$ along the length of the jet. (Here $`\beta _\mathrm{\Gamma }=(1\mathrm{\Gamma }^2)^{1/2}`$ and $`\theta `$ is the angle between the directions of the jet and the observer.)

The synchrotron energy-loss rate of relativistic nonthermal electrons with a random pitch-angle distribution in a region with comoving mean intensity $`B`$ is $`\dot{\gamma }_{syn}1.3\times 10^{21}B_{\mu \mathrm{G}}^2\gamma ^2`$ s-1 $`k_{syn}\gamma ^2`$ (Blumenthal & Gould, 1970). Under the same conditions, the energy loss rate due to Compton-scattered CMB radiation in the Thomson limit is $`\dot{\gamma }_\mathrm{T}`$ = $`(4c\sigma _\mathrm{T}/3m_ec^2)\gamma ^2(1+z)^4\widehat{u}_{CMB}\mathrm{\Gamma }^2(1+\beta _\mathrm{\Gamma }^2/3)`$ $`1.7\times 10^{20}\gamma ^2(1+z)^4\mathrm{\Gamma }^2\mathrm{s}^1`$ $`k_\mathrm{T}\gamma ^2`$, where $`\widehat{u}_{CMB}4\times 10^{13}`$ ergs cm-3 is the local ($`z=0`$) intensity of the CMBR. This expression holds when $`\gamma \gamma _{\mathrm{KN}}`$, where $`\gamma _{\mathrm{KN}}`$ is the Lorentz factor above which KN effects become important, given through $`4\gamma _{\mathrm{KN}}\widehat{ϵ}^{}=1`$. For the CMB radiation, the dimensionless mean photon energy $`\widehat{ϵ}^{}\widehat{ϵ}_{CMB}\mathrm{\Gamma }(1+z)2.7(k_\mathrm{B}\times 2.72\mathrm{K})\mathrm{\Gamma }(1+z)/m_ec^2`$ and $`\widehat{ϵ}_{CMB}=1.24\times 10^9`$, so that electrons with $`\gamma \gamma _{\mathrm{KN}}=2\times 10^8/\mathrm{\Gamma }(1+z)`$ scatter CMB radiation in the Thomson regime.

In the extreme KN limit, the electron energy loss rate $`\dot{\gamma }_{\mathrm{KN}}4.6\times 10^3(1+z)^2\mathrm{ln}[2.54\gamma _{10}\mathrm{\Gamma }(1+z)]`$ s-1 (Blumenthal & Gould, 1970), where $`\gamma _{10}=\gamma /10^{10}`$, and $`\mathrm{\Theta }=k_\mathrm{B}(2.72)(1+z)/m_ec^2=4.6\times 10^{10}(1+z)`$ is the dimensionless temperature. We approximate the combined energy-loss rate for analytic simplicity by the expression

$$\dot{\gamma }=[k_{syn}+\frac{k_\mathrm{T}}{1+(a\gamma /\gamma _{\mathrm{KN}})^b}]\gamma ^2K_{syn\mathrm{C}}(\gamma )\gamma ^2,$$
(1)
where $`a=0.5`$ and $`b=1.7`$ are chosen to fit roughly the $`\dot{\gamma }\mathrm{ln}\gamma `$ behavior in the extreme KN regime and to provide spectral results within a factor $`2`$-3 of the numerical results. Equation (1) shows that when $`k_\mathrm{T}k_{syn}`$, $`\dot{\gamma }`$ is dominated by Thomson losses at $`\gamma \gamma _{\mathrm{KN}}`$, $`\dot{\gamma }`$ flattens at $`\gamma _{\mathrm{KN}}\gamma \overline{\gamma }=\gamma _{KN}\sqrt{k_\mathrm{T}/k_{syn}}\gamma _{KN}\sqrt{N_{Ts}}`$, and $`\dot{\gamma }`$ is dominated by synchrotron losses at $`\gamma \overline{\gamma }`$. Adiabatic losses are neglected here but can be shown to be negligible for X-ray emitting electrons.

For numerical calculations in the next section, we solve the time-dependent continuity equation for the energy spectrum of electrons $`N_e(\gamma ,t)`$ injected in the source with differential injection rate $`Q(\gamma ,t)`$. For the analytic model treated here, we approximate the injection rate as a stationary function starting from $`t_0=0`$, with a single power-law behavior $`Q(\gamma )=Q_0\gamma ^p\mathrm{exp}(\gamma /\gamma _{max})`$ for $`\gamma \gamma _{min}`$, implying a comoving frame electron power $`L_e=m_ec^2_{\gamma _{min}}^{\gamma _{max}}𝑑\gamma \gamma Q(\gamma )`$. In the stationary injection case with $`p>1`$, the energy distribution of electrons at time $`t`$ can be approximated by $`N_e(\gamma ,t)Q(\gamma )t_{\mathrm{ac}}`$, where the characteristic electron accumulation time $`t_{\mathrm{ac}}=\mathrm{min}(t,t_{cool}=\gamma /\dot{\gamma })`$. In this approximation

$$\gamma ^3N_e(\gamma )Q_0\{\begin{array}{cc}\frac{\gamma ^{3p}}{[K_{syn\mathrm{C}}(\gamma _{cool})\gamma _{cool}]},\hfill & \text{for }\gamma _{min}\gamma \gamma _{cool}\hfill \\ & \\ \frac{\gamma ^{2p}}{K_{syn\mathrm{C}}(\gamma )},\hfill & \text{for }\gamma _{cool}\gamma \gamma _{max}\hfill \end{array}$$
(2)
where $`\gamma _{cool}\gamma _{cool}(t)`$ is found from the equation $`t_{cool}(\gamma _{cool})=t`$. Equation (2) applies when $`\gamma _{min}<\gamma _{cool}`$ and is easily generalized in the opposite case. We consider the case where $`\gamma _{cool}\gamma _{\mathrm{KN}}\gamma _{max}`$, and $`\gamma _{max}=\sqrt{3e_{max}e/\sigma _\mathrm{T}B}4.6\times 10^{10}\sqrt{e_{max}/B_{\mu \mathrm{G}}}`$ (e.g., de Jager et al., 1996), where the parameter $`e_{max}1`$.

In the $`\delta `$-function approximation for the synchrotron and Thomson radiation processes, the $`\nu F_\nu `$ synchrotron radiation spectrum from a uniform blob, assumed spherical in the comoving frame, is

$$f_ϵ^{syn}\delta ^4\left(\frac{c\sigma _\mathrm{T}u_B}{6\pi d_L^2}\right)\gamma _{syn}^3N_e(\gamma _{syn}),$$
(3)
where $`\gamma _{syn}=\sqrt{(1+z)ϵ/(\delta ϵ_B)}`$, $`d_L`$ is the luminosity distance, $`ϵ=h\nu /m_ec^2`$, and $`ϵ_B=B/4.414\times 10^{13}`$G. The $`\nu F_\nu `$ Thomson radiation spectrum for an external isotropic monochromatic radiation field is

$$f_ϵ^\mathrm{T}\delta ^6\left(\frac{c\sigma _\mathrm{T}u_{}}{6\pi d_L^2}\right)\gamma _\mathrm{C}^3N_e(\gamma _\mathrm{C})$$
(4)
(Dermer et al., 1997; Dermer & Schlickeiser, 2001). Here $`\gamma _\mathrm{C}=\delta ^1\sqrt{(1+z)ϵ/(2ϵ_{})}`$, $`ϵ_{}=\widehat{ϵ}_{CMB}(1+z)`$ and $`u_{}=\widehat{u}_{CMB}(1+z)^4`$ are the mean photon energy and radiation energy density, respectively, of the ambient CMBR field in the stationary frame. An accurate representation of the Compton-scattered spectrum in the KN regime is given by Georganopoulos, Kirk, & Mastichiadis (2001).

The KN effects on the synchrotron spectrum are observed at $`\nu _{\mathrm{KN}s1}`$(Hz) $``$ $`2\delta B_{\mu \mathrm{G}}\gamma _{\mathrm{KN}}^2/(1+z)`$ $`8\times 10^{16}\delta B_{\mu \mathrm{G}}/[\mathrm{\Gamma }^2(1+z)^3]`$. The reduction of Compton losses compared to the $`\dot{\gamma }\gamma ^2`$ behavior hardens the electron spectrum until synchrotron losses dominate at $`\gamma >\overline{\gamma }`$, resulting again in a steepening of the synchrotron spectrum in the hard X-ray domain. The KN effects on the Compton-scattered CMBR spectrum are observed at dimensionless photon energy $`ϵ_{\mathrm{KNC}}=h\nu _{\mathrm{KNC}}/m_ec^210^8(\delta /\mathrm{\Gamma })(1+z)^2`$, implying a hardening in the photon spectrum at photon energy $`E>E_{\mathrm{KNC}}50(\delta /\mathrm{\Gamma })(1+z)^2`$ TeV. We set the normalization $`Q_0=(p2)10^{44}L_{44}\mathrm{ergs}\mathrm{s}^1/[(m_ec^2)(\gamma _{min}^{2p}\gamma _{max}^{2p})]`$ and let $`p=2.3`$, corresponding to the likely spectral index of particles accelerated by relativistic shocks (Achterberg et al., 2001).

Figs. 1a and 1b show two suites of models inspired by the Chandra data. In Fig. 1a, we assume that $`\mathrm{\Gamma }=10`$ in a nearby radio galaxy with $`z=0.15`$. We also assume that the jet axis is inclined at an angle $`\theta =4/\mathrm{\Gamma }=23^{}`$. The comoving magnetic field ranges in values from 5 to 50 $`\mu `$G, and effects of different values of $`\gamma _{cool}`$ are illustrated. A field of 44 $`\mu `$G gives equipartition between the magnetic-field and comoving CMBR energy densities. The spectral hardening between the optical and X-ray regimes is apparent. Note that the use of the $`\delta `$-function approximation in the analytic approach enhances spectral features.

The $`\mathrm{\Gamma }=10`$ jet with $`B_{\mu \mathrm{G}}=30`$ and $`\gamma _{cool}=10^5`$ is placed at $`z=1`$ in Fig. 1b, but now with the outer jet oriented at various angles to the observer. The increasingly dominant Compton component at small angles results from the different beaming factors for synchrotron and Compton processes (Dermer, 1995). For sufficiently small values of $`\theta `$ and $`\gamma _{min}`$, Thomson-scattered X-ray CMBR could make a significant or dominant contribution to knot and hot spot emission, as in knot WK7.8 of the $`z=0.651`$ superluminal source PKS 0637-752, with apparent transverse speeds reaching $`18c`$ (Schwartz et al., 2000; Tavecchio et al., 2000). When observing too close within the beaming cone, the direct inner jet radiation may however dominate, which could result in a pattern of $`\gamma `$-ray flares superimposed on a significantly weaker but persistent (on time scales of thousands of years) flux.

## 3 Numerical Model

Numerical calculations are done using the well-known solution for the energy distribution $`N_e(\gamma ,t)`$ of electrons suffering energy losses $`P(\gamma )(\dot{\gamma })`$, which are contributed mostly by Compton and synchrotron losses. Calculations are done in the comoving frame of a source moving relativistically with Lorentz-factor $`\mathrm{\Gamma }`$, and then the emerging radiation spectra are transformed to the observer frame. We neglect possible escape losses of accelerated particles. This is valid for the large characteristic size of the X-ray knots of order $`R_{knot}1`$-10 kpc or so, because for $`\gamma >10^4`$ the electron cooling time $`t_{cool}10^5\mathrm{yr}`$ in the comoving frame, which is smaller than any reasonable escape time from such a large source. Particle escape may be important in much smaller-scale inner jets of blazars (Atoyan & Dermer, 2001, 2002) or jets of microquasars (Atoyan & Aharonian, 1999).

Gamma rays with energies $`10^{15}\mathrm{\Gamma }^1(1+z)^1\mathrm{eV}`$ could be absorbed in $`\gamma \gamma `$ collisions with the CMB photons inside the source. This process provides an injection function $`Q_{1,\gamma \gamma }(\gamma ,t)`$ for the first generation electrons $`(e^+,e^{})`$ of the pair-photon cascade, which we also take into account in our calculations. Note however that the cascade radiation turns out to be insignificant compared with the radiation from the main injection $`Q(\gamma ,t)`$ unless a large fraction of the overall power is injected at electron energies $`100`$ TeV (Atoyan & Dermer, 2002).

We consider an electron injection spectrum

$$Q(\gamma ,t)\gamma ^pe^{\gamma /\gamma _{max}}(1+t/t_{inj})^q$$
(5)
at energies $`\gamma _{min}<\gamma `$ extending with index $`p=2.3`$ to PeV energies, and we set $`\gamma _{max}=2\times 10^9`$. The injection function is approximated by $`Q(\gamma /\gamma _{cool})^2`$ at $`\gamma \gamma _{min}=\gamma _0\mathrm{\Gamma }`$ with $`\gamma _0=300`$. Note that for electron acceleration by relativistic shocks, the values of $`\gamma _0`$ could reach a large fraction of $`(m_p/m_e)\mathrm{\Gamma }2000\mathrm{\Gamma }`$. The time profile of the injection in equation (5) treats both stationary (at $`t0`$) injection if $`t_{inj}t`$, as well as gradually declining injection at $`t_{inj}t`$, when $`q>0`$.

In Fig. 2a we show the synchrotron and Compton fluxes expected from a relativistic knot with $`\mathrm{\Gamma }=10`$ in a distant blazar at $`z=1`$ for three different times in the knot frame: $`t^{}=10^4\mathrm{yr}`$, $`5\times 10^4\mathrm{yr}`$, and $`3\times 10^5`$ yr . If the injection starts relatively close to the core, then these timescales would effectively correspond to knots at distances $`l=ct^{}\mathrm{\Gamma }33\mathrm{kpc}`$, $`166`$ kpc, and $`1`$ Mpc from the core. Injection is assumed to be stationary with $`L_{44}=1`$, and the magnetic field and the jet angle to the observer are $`B=30\mu \mathrm{G}`$ and $`\theta =11.4^{}`$ respectively. Fig. 2b shows the fluxes expected at the same epochs $`t^{}`$ from a closer blazar at $`z=0.15`$, calculated for a smaller magnetic field $`B=6\mu \mathrm{G}`$ and $`\mathrm{\Gamma }=5`$, and with $`\theta =23^{}`$. The electron injection here is assumed to decline with $`q=1/2`$ and with a characteristic timescale $`t_{inj}=5\times 10^3\mathrm{yr}`$, for initial power $`L_{44}=5`$. Fig. 2b demonstrates that the X-ray and optical fluxes may significantly drop while the radio flux continues to rise with increasing distance, as observed in 3C 273 (Sambruna et al., 2001). At the same time, the X-ray fluxes may be systematically harder than the optical spectra, again in qualitative agreement with observations of 3C 273.

## 4 Discussion and Summary

We have shown that the combined effects of Compton and synchrotron losses on a power-law electron injection spectrum can produce a hardening between the optical and X-ray regimes. This model is in accord with observations of X-ray spectral hardenings in the knots of 3C 273 (Sambruna et al., 2001; Marshall et al., 2001) and other sources, suggesting that the X-ray emission is due to a cooling spectrum of electrons accelerated by a strong shock. In a synchrotron model for the X-rays, the X-ray profiles are narrower than radio profiles, because the X-rays are emitted by higher energy electrons and therefore cool more rapidly than the radio-synchrotron emitting electrons.

The Thomson energy-loss rate must exceed the synchrotron energy loss-rate to produce X-ray spectral hardenings, so a large fraction of the energy in nonthermal electrons is radiated as $`\gamma `$-rays with energies $`50(\delta /\mathrm{\Gamma })(1+z)^2`$ TeV. The Compton-scattered CMB radiation from the knots and hot spots make nearly aligned extended jet sources potentially detectable with GLAST and the next generation of ground-based air Cherenkov $`\gamma `$-ray telescopes such as VERITAS, HESS, and MAGIC. The relative intensities of the X-ray and $`\gamma `$-ray fluxes can be used to infer $`\delta `$, but the limited spatial resolution of gamma-ray telescopes will pose difficulties in separating core jet components (excepting those which are highly variable) from steady extended emission components. SIRTF will be important to offer greater spectral detail about the synchrotron component, and to correlate spectral-cooling breaks in the synchrotron component with the break in the Thomson component observed with GLAST.

The CMB radiation scattered in the Klein-Nishina regime will provide a source of ultra-high energy gamma-rays that can pair produce in the diffuse CMB and infrared radiation fields to form pair halos (Aharonian, Coppi, & Völk, 1994), as well as providing a source of synchrotron-emitting electrons. The mean-free-path for $`\gamma \gamma `$ attenuation in the CMB reaches a minimum value of $`8\mathrm{kpc}/(1+z)^3`$ at observed photon energies $`E_\gamma 1`$ PeV, and increases to Mpc scales at lower and higher energies. In the model of Atoyan & Dermer (2001, 2002), the injection of energy into the extended jets of FRII radio galaxies is a consequence of neutral beams composed of ultra-high energy neutrons and $`\gamma `$-rays ($`10^{13.5}`$-$`10^{18}`$eV) formed through photomeson production in the inner jet, explaining the colinearity of the inner and extended jets up to Mpc scales. The highly collimated neutral-beam energy is deposited in the intergalactic medium to drive relativistic outflows and form shocks that accelerate high-energy particles which radiate X-ray emission observed with Chandra.

We thank Hui Li for discussions about the synchrotron radiation from high-energy pair cascades, Andrew Wilson for discussions about the Chandra data, and the referee for a useful report. AA appreciates the hospitality and support of the NRL Gamma and Cosmic Ray (now High Energy Space Environment) Branch during his visit when this work has been done. The work of CD is supported by the Office of Naval Research and NASA grant No. DPR S-13756G.


text:
# Detection of superimposed periodic signals using wavelets

## 1 Introduction

The search for periodic signals is common in many areas of astronomy, and, often, simple Fourier based methods or epoch folding methods are able to detect the periods that are present in the data sets. However, this turns out to be difficult when the periodic signals are not sinusoidal or when the signal-noise ratio is very low.

Specifically, in high energy astrophysics two problems are responsible for reduced signal-noise ratios in the data: the relatively small number of high energy photons when compared to the number of photons available at other wavelengths, and the low efficiency in photon counting of high energy detectors. Hence, the astronomical data supplied by X-ray telescopes, especially flux monitorings, often suffer from poor statistics. This poses many problems when processing and analysing the data in order to determine the flux, detect periodic signals, etc. In addition, some X-ray sources present several periods, which arise from different physical processes such as pulses, eclipses in two-body systems or occultation by a precessing accretion disc. Some examples of this kind of sources are SMC X-1 (Wojdowski et al., 1998), LMC X-4 (Levine et al., 1991), and LS I +65010 (Corbet, Finley & Peele, 1999).

Classical period determination methods are based either on epoch folding techniques or on Fourier decomposition analysis. The first group of methods are based on the analysis of the dispersion of the diferent light curves produced by folding the data over a range of trial periods (see, for example, Lafler & Kinman 1965; Jurkevich 1971 and Stellingwerf 1978). The second type of methods use the Fourier transform in combination with deconvolution techniques to deal with the data sampling function (see, for example, Lomb 1976; Scargle 1982, 1989; Roberts, Lehár & Dreher 1987 and Press & Rybicky 1989).

Epoch folding methods and Fourier based methods work well when applied to data sets that present a unique period. However, when the analysed data set contains several periodic signals, the behavior of these algorithms highly depends on intrinsic signal characteristics. Every individual periodic signal may present different spectral behavior (amplitude and profile) and depending on them, the algorithm will detect some kinds of signal better than others. Fourier based techniques will, in general, successfully detect two combined sinusoidal signals, even if one of them has a low signal-noise ratio. However, they have difficulties in detecting the non-sinusoidal signals that might be present in a data set (as in the astronomical case of pulsed emission superimposed on an orbital period). On the contrary, epoch folding methods are well suited for detecting non-sinusoidal signals but they tend to fail when two or more periods are present in the data set, especially when the signal-noise ratio is low.

In general, we can say that the greater the difference between the spectral characteristics of each signal, the more difficult it is to detect each period with classical methods, and the less statistically significant are the results given by them. Therefore, a new algorithm becomes immediately necessary for the detection of superimposed periodic signals.

In this paper, we show how the methodology of wavelet theory is very well suited to this problem, since it is completely oriented to decompose functions into their several spectral characteristics. This allows us to isolate every signal present in our data and to analyse them separately, avoiding their mutual influences. In Section 2 we outline some concepts in wavelet theory relevant to the stated problem. In Section 3 we propose an algorithm to detect each of the periodic signals present in a data set by combining wavelet decomposition with classical period determination methods. In Sections 4 and 5 we present some examples of synthetic data we used to test the algorithm and the results we obtained. We summarise our conclusions in Section 6.

## 2 Outline of the wavelet transform

Multiresolution analysis based on the wavelet theory introduces the concept of details between successive levels of scale or resolutions (Chui 1992; Daubechies 1992; Meyer 1993; Young 1993; Kaiser 1994; Vetterli & Kovacevic 1995).

Wavelet decomposition is increasingly being used for astronomical signal processing (Starck & Murtagh 1994; Starck, Murtagh & Bijaoui 1998) and remotely sensed data (Yocky 1995; Datcu, Luca & Seidel 1996). The method is based on the decomposition of the signal into multiple channels based on their local frequency content. The wavelet transform is an intermediate representation between the Fourier and the temporal one. In the Fourier transform we know the global frequency content of our signal, but we have no information about the temporal localization of these frequencies. However, the wavelet transform gives us an idea of both the local frequency content and the temporal distribution of these frequencies. Since in Fourier space the basis functions are sinusoidal, they extend along all space and do not have temporal concentration. However, wavelets are concentrated around a central point, thus they have a high degree of temporal localization.

### 2.1 Some theory

Next, we will briefly present an outline of the wavelet transform, relying on the multiresolution signal representation concept.

Given a signal $`f(t)`$, we construct a sequence $`F_m(f(t))`$ of approximations of $`f(t)`$. Each $`F_m(f(t))`$ is specific for the representation of the signal at a given scale (resolution). $`F_m(f(t))`$ represents the projection of the signal $`f(t)`$ from the signal space $`S`$ onto subspace $`S_m`$. In this representation, $`F_m(f(t))`$ is the ‘closest’ approximation of $`f(t)`$ with resolution $`2^m`$.

The differences between two consecutive scales $`m`$ and $`m+1`$ are the multiresolution wavelet planes or ‘detail’ signal at resolution $`2^m`$:

$$w_m(f(t))=F_m(f(t))F_{m+1}(f(t)).$$
(1)
This detail signal can be expressed as:

$$w_m(f(t))=\underset{l}{}W_{m,l}(f)\psi _{m,l}(t),$$
(2)
where coefficients $`W_{m,l}(f)`$ are given by the direct wavelet transform of the signal f(t):

$$W_{m,l}(f)=_{\mathrm{}}^+\mathrm{}f(t)\psi _{m,l}(t)𝑑t.$$
(3)
The coefficients $`W_{m,l}(f)`$ are called wavelet coefficients of $`f(t)`$. Such coefficients correspond to the fluctuations of the signal $`f(t)`$ near the point $`l`$ at resolution level $`m`$. Thus, the wavelet transform (3) represents the expansion of signal $`f(t)`$ in the set of basis functions $`\psi _{m,l}(t)`$. These basis functions are scaled and translated versions of a general function $`\psi (t)`$ called Mother Wavelet. To construct the basis functions $`\psi _{m,l}(t)`$, the Mother Wavelet is dilated and translated according to parameters $`m`$ and $`l`$ as follows:

$$\psi _{m,l}(t)=2^{m/2}\psi (2^mtl).$$
(4)
Therefore, all the basis functions $`\psi _{m,l}(t)`$ have the same profile, that is, the Mother Wavelet profile. Using (4) we obtain an orthonormal wavelet basis (Daubechies, 1992).

The inverse discrete wavelet transform is given by the reconstruction formula:

$$f(t)=\underset{m}{}\underset{l}{}W_{m,l}(f)\psi _{m,l}(t).$$
(5)
In summary, the wavelet transform describes at each resolution step the difference between the previous and the current resolution representation. By iterating the process from the highest to the lowest available resolution level we obtain a pyramidal representation of the signal. This usually includes a decimation process, i.e., in each iteration 1 out of 2 points is discarded, which implies that the number of data points at lower frequencies is highly reduced.

In the Fourier transform, the noise contribution is spread along all frequencies. In contrast, one of the interesting properties of the wavelet transform (which is also frequency based) is that noise contribution is more important on the higher frequency wavelet planes.

### 2.2 The ‘à trous’ algorithm

In order to obtain a discrete wavelet decomposition for signals, we follow Starck & Murtagh (1994) and we use the discrete wavelet transform algorithm known as ‘à trous’ (‘with holes’) to decompose the signal into wavelet planes. Given a signal $`p`$ we construct the sequence of approximations:

$$F_1(p)=p_1,F_2(p_1)=p_2,F_3(p_2)=p_3,\mathrm{}.$$
(6)
To construct the sequence, this algorithm performs successive convolutions with a filter obtained from an auxiliary function named scaling function. We use a scaling function which has a $`B_3`$ cubic spline profile. The use of a $`B_3`$ cubic spline leads to a convolution with a mask of 5 elements, all elements being scaled up by 16: (1,4,6,4,1).

As stated above, the wavelet planes are computed as the differences between two consecutive approximations $`p_{i1}`$ and $`p_i`$. Letting $`w_i=p_{i1}p_i(i=1,\mathrm{},n)`$, in which $`p_0=p`$, we can write the reconstruction formula:

$$p=\underset{i=1}{\overset{n}{}}w_i+p_r.$$
(7)
In this representation, the signals $`p_i(i=0,\mathrm{},n)`$ are versions of the original signal $`p`$ at increasing scales (decreasing resolution levels), $`w_i(i=1,\mathrm{},n)`$ are the multiresolution wavelet planes and $`p_r`$ is a residual signal (in fact $`n=r`$, but we explicitly substitute $`n`$ by $`r`$ to clearly express the concept of $`residual`$). In our case, we are using a dyadic decomposition scheme. Thus, the original signal $`p_0`$ has double resolution than $`p_1`$, the signal $`p_1`$ double resolution compared to $`p_2`$, and so on. If the resolution of signal $`p_0`$ is, for example, $`10\mathrm{\Delta }t`$ (being $`\mathrm{\Delta }t`$ the sampling of the data), the resolution of $`p_1`$ would be $`20\mathrm{\Delta }t`$, the resolution of $`p_2`$ would be $`40\mathrm{\Delta }t`$, etc. All these $`p_i(i=0,\mathrm{},n)`$ signals have the same number of data points, in contrast to some more usual wavelet decomposition algorithms which reduce the number of points by a factor 2 when going through increasing scales (process known as decimation step). Later we will see why we are interested in maintaining the number of points.

## 3 Proposed algorithm

We propose to apply the wavelet decomposition using the ‘à trous’ algorithm to solve our initial stated problem: to isolate each of the periodic signals contained in a set of data and study them separately.

Hence, the algorithm we propose is as follows:

1. Choose a value for $`n`$ and decompose the original signal $`p`$ into its wavelet planes $`w_i(i=0,\mathrm{},n)`$.
2. Detect periods in each of the $`n`$ obtained wavelet plane $`w_i`$.

Any method can be used to detect the periods present in each $`w_i`$. In our tests, we used Phase Dispersion Minimization (PDM) (Stellingwerf, 1978) and CLEAN (Roberts et al., 1987) methods. PDM belongs to the group of epoch folding methods mentioned in the introduction. Therefore it folds the data points over a set of trial periods and analyses the dispersion of every light curve in the phase space. To do that, the data points are grouped in phase bins and the most probable period is then found by minimizing the sum of the bin variances. On the other hand, CLEAN belongs to the group of Fourier transform decomposition methods. It works in the frequency domain and it is based on the fact that the spectrum obtained from the Fourier transformation of the data set contains artefacts caused by the imcompleteness of the sampling (time gaps) and the finite time span of it. CLEAN uses a non-linear deconvolution algorithm to clean (hence its name) these artefacts from the original (dirty) spectrum. This is done by interatively substracting from the dirty spectrum the expected response of a signal composed by a unique harmonical function (clean component). From the set of clean components, a clean spectrum is obtained which is at the end corrected to preserve the noise level and the frequency resolution.

PDM is well suited to identify signals with any profile, not just a sinusoidal one. For instance, if we have a burst-like periodic signal, it is more successfully detected by the PDM algorithm than by CLEAN. On the other hand, if we knew that we are trying to estimate the period of a sinusoidal signal, it would be better to use CLEAN.

The algorithm we present here could be seen as an improvement over the usual period detection algorithms, because prior to using them we decompose our signal into its $`w_i`$ wavelet planes, and afterwards we apply these algorithms to each $`w_i`$.

In PDM and CLEAN methods there is only one data set to study and to estimate its period. But in the wavelet based method there are several data sets to study (several wavelet planes), and this improves the detection probabilities. We can search for a certain period in several wavelet planes, which helps to discard some of the usual marginal artefact detections present in PDM and CLEAN and to detect the true periods present in our original data.

We can use other wavelet based algorithms (Szatmáry, Vinkó & Gál, 1994), which are based on approximations of a continuous wavelet transform and the study of wavelet space coefficients. However, these algorithms present a non direct inverse wavelet transform: the search for periodicities is based on the fit between the wavelet base function profile and the signal one, and therefore on the values of the wavelet transform coefficients, which highly depend on the wavelet base used. For example, if we took a signal constructed by a Gaussian with 1 day duration which repeats every 30 days, it would be very difficult to find a wavelet base that fits well to this sparse signal. On the contrary, if we take a suitable wavelet base and a decomposition scheme which allows us to perform easily its inverse transform, this base will fit the profile of every Gaussian. Hence, performing the inverse transform of every wavelet scale, we can work in the temporal representation to find directly the periodicities. This is the case of the wavelet decomposition algorithm we use in the method we present here.

Another reason to use the ‘à trous’ algorithm is that the sampling of the wavelet planes is the same as the original data. This allows us to work directly in the temporal space with the frequency content of the corresponding wavelet plane. In other wavelet decomposition algorithms, we are usually forced to work on the decimated wavelet space.

Hereafter, and for notational convenience, the wavelet based PDM and CLEAN methods will be called WPDM and WCLEAN respectively.

## 4 Simulated data

In order to check the benefit of applying WPDM versus PDM, or WCLEAN versus CLEAN, we generated several sets of simulated data containing two superimposed periodic signals. Each data set is composed of a high amplitude primary sinusoidal function and a secondary low amplitude one. For the latter we used two kinds of functions: sinusoidal and Gaussian. The first are intended to simulate variable sources with pure sinusoidal intensity profile (like precession of accretion discs), and the second burst-like events (like pulses or eclipses). Finally, we added a white-Gaussian noise to this combination of signals. We increased the value of the noise standard deviation, $`\sigma `$, up to the value where detection of periodic signals became statistically insignificant in both the classical and the wavelet based methods.

Here we list the characteristics of the signals:

1. Each signal is generated as an evenly spaced data set composed of 1000 points, separated by a unit between them.
2. The high amplitude sinusoidal function has an amplitude equal to 1 and a period of 108.5 units. The period selected should satisfy two conditions: not be an exact divisor of the number of points in the data set, and allow the presence of more than 5 complete periods in the simulated signal to guarantee the possibility of detection.
3. The amplitudes of the low amplitude periodic functions (sinusoidal or Gaussian) are 0.1 and 0.5.
4. The periods used for the secondary functions are 13.13 and 23.11 units. Others could be used, satisfying the condition that they are not harmonics of the high amplitude sinusoidal function period.
5. In the case of the Gaussian signal, we have used two values for the Full Width Half Maximum (FWHM), corresponding to 2 and 6 units, respectively.

In the first three columns of Table 1 we present the parameters used to generate each simulated data set in the case of sine+sine periodic signals. In Table 2 the first four columns show the parameters used to generate the sine+Gaussian signals.

## 5 Results

The simultaneous use of two independent methods, such as CLEAN and PDM, is usually applied to discriminate false period detections from the true ones. A similar procedure can be used with each one of the wavelet planes in the WPDM and WCLEAN methods. Therefore, when comparing the behavior of these different methods, we have to compare the usual PDM-CLEAN method combination for period estimation prior to the new WPDM-WCLEAN combination.

The primary period (108.5 units) is always detected by all methods, and does not appear in Table 1 and Table 2. In these tables, the four last columns show the detected low amplitude periods for each data set using PDM, CLEAN, WPDM and WCLEAN methods, respectively. A dash is shown when a period is not detected, and a question mark when the detection is difficult or doubtful. When a period is indicated in the wavelet based methods, we also show in parentheses the wavelet planes where it is detected.

### 5.1 Sine+sine case

To illustrate the performances of WPDM and WCLEAN methods, we show in Fig. 1 a simulated data set which contains two pure sinusoidal functions without noise, and both their PDM and CLEAN outputs (in PDM we look for minima with subharmonics, while in CLEAN we look for maxima). The 108.5 units period is detected by both methods. In the case of PDM we have zoomed into the range of short trial periods and dispersion to try to detect the 23.11 units minimum. It is clear that even in the case of signals without noise, PDM method is unable to find the secondary signal (period=23.11 units). In contrast, CLEAN is able to find it, although with a low normalized power (less than 0.1, which is the ratio between the amplitudes of both periodic signals). The output of WPDM is shown in Fig. 2 and the one of WCLEAN in Fig. 3. As one would expect, the secondary period is only seen in the lower wavelet planes ($`w_1`$ to $`w_4`$), while the primary long-term period of 108.5 units is better detected when working with higher wavelet planes ($`w_4`$ to $`w_6`$). We have worked until $`w_6`$, because in this plane the 108.5 units period is clearly detected, and there is no reason to continue with higher values of $`n`$. This will be the case in all the analysed data we present here.

Some remarks can be made after inspection of Table 1:

1. There is a better performance of CLEAN relative to PDM, because the latter has problems when dealing with superimposed low amplitude periodic signals. Only when the amplitude of the secondary signal is close enough to the primary one (amplitude=0.5) is PDM able to detect it.
2. In all methods, with high noise-signal ratios the detected periods are slightly different from the simulated ones.
3. WPDM and WCLEAN perform always better than PDM and CLEAN methods. When CLEAN marginally detects the secondary period, WPDM and WCLEAN have no problems to detect it, and they work properly even with higher noise. In the WPDM case, the results are always much better than with PDM.
4. As the noise increases, the detection starts to fail in the lower wavelet planes (higher frequencies), and only the higher ones (lower frequencies) are noise-free enough to allow period detection.
5. The maximum ratio between the noise and the signal amplitude to detect a period with WPDM and WCLEAN, is nearly constant for a given period, being around 12 for the 13.13 units period and around 8 for the 23.11 units period. The reason for this difference is that there are 76 complete periods in the case of 13.13 units period, and only 43 complete periods in the case of 23.11 units period.

### 5.2 Sine+Gaussian case

For illustrative purposes, we show in Fig. 4 the wavelet decomposition of a sine+Gaussian signal. In the lower wavelet planes (planes $`w_1`$ and $`w_2`$) we have isolated the Gaussian contribution. In the $`w_3`$ we have a mixed contribution from both periods, but in the higher planes we only have the longest period contribution.

We must note that the use of two different FWHM for the Gaussian, combined with two different periods (13.13 and 23.11 units), gives 4 different profiles. Hence, the phase duration of the burst-like event ranges from very low to relatively high values in the following order: FWHM=2 and period=23.11, FWHM=2 and period=13.13, FWHM=6 and period=23.11, and finally FWHM=6 and period=13.13.

In view of the results displayed in Table 2 we can make the following comments:

1. Again, there is a better performance of CLEAN over PDM, for the reason explained above. However, we must note that when the FWHM is only 2 units, PDM never detects the secondary period. Only with FWHM=6 units and an amplitude of 0.5, can PDM detect the low amplitude periodic signals.
2. In all methods, with high noise-signal ratios the detected periods are slightly different from the simulated ones.
3. WPDM and WCLEAN perform always better than PDM and CLEAN methods. When CLEAN marginally detects the secondary period, WPDM and WCLEAN have no problems to detect it, and they work properly even with higher noise. In the WPDM case, the results are always much better than with PDM.
4. As the noise increases, the detection starts to fail in the lower wavelet planes (higher frequencies), and only the higher ones (lower frequencies) are noise-free enough to allow period detection.
5. In all cases with amplitude=0.1, the WPDM performance is very similar to WCLEAN. In the amplitude=0.5 cases, WPDM is always better than WCLEAN because the signal is not sinusoidal and the amplitude is high enough to allow a good detection.
6. For a given amplitude of the Gaussian signal, the maximum noise-signal ratio achieved with WPDM and WCLEAN increases with the phase duration of the FWHM.

Finally, and for illustrative purposes, we show in Fig. 5 and Fig. 8 two simulated data sets generated with the following common parameters: 13.13 units period, FWHM=6.0 units and amplitude=0.5. The only difference is that in Fig. 5, $`\sigma _{\mathrm{n}oise}=0.5`$, while in Fig. 8, $`\sigma _{\mathrm{n}oise}=1.5`$. The outputs of PDM, CLEAN, WPDM and WCLEAN are shown from Fig. 5 to Fig. 10.

As can be seen in the figures, in the case $`\sigma _{\mathrm{n}oise}=0.5`$, the 13.13 units period is clearly detected by all methods. In WPDM and WCLEAN the period is detected in the wavelet planes $`w_2`$ and $`w_3`$. These detections are much better than those of PDM and CLEAN.

On the other hand, in the $`\sigma _{\mathrm{n}oise}=1.5`$ case, the period is not detected by PDM nor by CLEAN. WCLEAN also fails to detect it: the maximum in the wavelet plane $`w_3`$ is not indicative of a period detection and, moreover, its position is in the 15.5 units trial period. The only clear detection with a subharmonic is in $`w_3`$ of WPDM.

## 6 Conclusions

When dealing with two superimposed periodic signals, the performance of classical methods, like PDM and CLEAN, depends on the characteristics of these signals. These methods often fail to detect simultaneously both periods, even with high or moderate signal-noise ratios.

With the wavelet based methods we present here, WPDM and WCLEAN, it is possible to reach lower values of signal-noise ratios and still detect both periodic signals. This means that the wavelet based methods are less affected by noise than PDM and CLEAN, and they allow us to detect periodic signals with noisier data.

Another advantage of the proposed methods is the simultaneous detection of periods in several wavelet planes. If a period is marginally detected in more than one wavelet plane, we are probably detecting a true periodic signal. This allows us to improve the detection confidence.

These facts prove the major ability of WPDM and WCLEAN over PDM and CLEAN to deal with noisy signals and to detect superimposed periods even with low signal-noise ratios.

## Acknowledgments

We acknowledge partial support from DGI of the Ministerio de Ciencia y Tecnología (Spain) under grant AYA2001-3092. This research has made use of facilities of CESCA and CEPBA, coordinated by C4 (Centre de Computació i Comunicacions de Catalunya). During this work, M. Ribó has been supported by a fellowship from CIRIT (Generalitat de Catalunya, ref. 1999 FI 00199).


text:
# In-plane Tunneling Spectrum into a [110]-Oriented High-𝑇_𝑐 Superconductor in the Pseudogap Regime

## Abstract

Both the differential tunneling conductance and the surface local density of states (LDOS) of a -oriented high-temperature superconductor in the pseudogap (PG) regime are studied theoretically. As a competing candidate for the mechanism of PG state, the charge-density wave (CDW), spin-density wave (SDW), $`d`$-density wave (DDW), and $`d`$-wave superconducting (DSC) orderings show distinct features in the tunneling conductance. For the CDW, SDW, and DSC orderings, the tunneling conductance approaches the surface LDOS as the barrier potential is increased. For the DDW ordering, we show for the first time that there exist midgap states at the  surface, manifesting themselves as a sharp zero-energy peak in the LDOS, as in the case of DSC ordering. However, due to the particle-hole pair nature of the DDW state, these states do not carry current, and consequently the one-to-one correspondence between the tunneling conductance and the surface LDOS is absent.

One of the most intriguing properties of high temperature superconductors is that, in the underdoped regime, the density of low-lying excitations is suppressed below a characteristic temperature $`T^{}`$, which could be well above the superconducting transition temperature $`T_c`$ Timusk99 . The origin for this so called pseudogap (PG) phenomenon remains unsettled. The existing scenarios for the mechanism of PG state include: the preformed pair model, where the Cooper pairs are formed above $`T_c`$ but either the phase associated with them fluctuates so strongly Emery95  or the coherence length (essentially the size of a Cooper pair) is so small Chen98  that the superconducting transition temperature is reduced from $`T^{}`$ to $`T_c`$; stripes Zaanen89 ; Poilblanc89 ; Machida89 ; Emery90 ; Martin00 ; Castellani00  or antiferromagnetic fluctuations Pines97 , relying on the spin and/or charge fluctuations; spin-charge separation Lee99 , where the PG comes from the spinon sector; quantum criticality Sachdev99 , of which the $`d`$-density wave (DDW) ordering breaking both time-reversal and translational symmetry Chakravarty01  and the circulating currents breaking time-reversal but preserving translational symmetry Varma99 , are the most representative models. Recently, several tunneling experiments have been proposed for both the preformed pair Choi00 ; Sheehy00  and antiferromagnetic correlation scenarios Bang00 . In this paper, based on a lattice model, we calculate the tunneling conductance, within the framework of a scattering approach, of a -oriented normal-metal/underdoped high-$`T_c`$ superconductor junction. We typically consider the possible orderings: (i) the charge-density wave (CDW) to model the charge stripes; (ii) the spin-density wave (SDW) to model the spin stripes or antiferromagnetic correlations; (iii) the DDW state; (iv) the $`d`$-wave superconducting (DSC) state. The exposure of distinct features among these orderings may help us pin down the mechanism for the PG. On the other hand, it has been predicted Hu94  that a sizable areal density of midgap states (i.e., states with essentially zero energy relative to the Fermi surface) exists on the -oriented surface of a $`d_{x^2y^2}`$-wave superconductor. These states gives rise to a narrow surface local density of states (LDOS) peak at the Fermi energy, where the bulk density of states dips to zero. One of the consequences of these midgap states is a sharp zero bias conductance peak in single particle tunneling Tanaka95 ; Xu96 , which has been confirmed by several carefully controlled experiments Covington97 ; Alff97 ; Sinha98 ; Wei98 . A question arises naturally: Whether such kind of midgap states also exist in those normal-state orderings, and if so, how to distinguish them from those in the DSC ordering. This question constitutes another motivation for the present work. In light of this, we also calculate the LDOS near the -oriented surface of the superconductor with all orderings listed above. We choose the  orientation, where the crystalline $`a`$ axis of the superconductor has a $`\frac{\pi }{4}`$ angle with respect to the direction normal to the surface, because the surface sensitivity of the electronic structure of some (e.g., the DDW and DSC) orderings can be manifested most pronouncedly in this setup.

The equation of motion for quasiparticles in a superconductor with a variety of orderings can be described in a unified manner by the generalized Bogoliubov-de Gennes equations defined in a two-dimensional tight-binding model:

$$\underset{j}{}\left(\begin{array}{cc}_{ij,\sigma }& \mathrm{\Delta }_{ij}\\ \mathrm{\Delta }_{ij}^{}& _{ij,\sigma }^{}\end{array}\right)\left(\begin{array}{c}u_{j,\sigma }\\ v_{j,\sigma }\end{array}\right)=E\left(\begin{array}{c}u_{i,\sigma }\\ v_{i,\sigma }\end{array}\right).$$
(1)
Here the quasiparticle wavefunction $`\left(\begin{array}{c}u_{i,\sigma }\\ v_{i,\sigma }\end{array}\right)`$ comprises the electrons at energy $`E`$ with spin $`\sigma `$ coupled with the holes of the same energy but opposite spin $`\sigma `$; $`i=(i_a,i_b)`$ are the lattice indices in the coordinate system defined by the crystalline $`a`$ and $`b`$ axes (Hereafter the lattice constant is chosen to be unity); the single particle Hamiltonian $`_{ij}=t\delta _{i+\tau ,j}+W_{ij}+(U_i+S_{i,\sigma })\delta _{ij}`$, where $`t`$ is the nearest-neighbor hopping integral and $`\tau =(\tau _a,\tau _b)=(\pm 1,0)`$ or $`(0,\pm 1)`$ represent the relative position of those sites nearest neighboring to the $`i`$th site; $`W_{ij}=(1)^{i_a+i_b}W_0\delta _{ij}`$ for the conventional CDW ordering while $`W_{ij}=(1)^{i_a+i_b}iW_d\text{sgn}(\tau _a^2\tau _b^2)`$ for the DDW ordering, $`S_{i,\sigma }=\sigma (1)^{i_a+i_b}S_0`$ represents the SDW order parameter; $`U_i`$ is the scattering potential to model the impurities, surface or interface. The quantity $`\mathrm{\Delta }_{ij}=\mathrm{\Delta }_s\delta _{ij}`$ for the $`s`$-wave superconducting order parameter while $`\mathrm{\Delta }_{ij}=\mathrm{\Delta }_d\text{sgn}(x_{ij}^2y_{ij}^2)`$ for the DSC order parameter. We note that the CDW, SDW, and DDW fields participate in the diagonal elements of the $`2\times 2`$ BdG matrix equations in spin space while the superconducting order parameter appears as off-diagonal elements. Of most interest to us in the present work, as explained above, we consider the tunnel junction comprising a normal metal on the left hand side and a superconductor on right hand side with a  oriented surface. We choose that the interface runs along the $`y`$ direction while the normal to the interface is parallel to the $`x`$ direction, along which the transport takes place. The insulating barrier at the interface is modeled by a scattering potential, $`U_i=U_0\delta _{i_x,0}`$. In addition, we assume that the order parameter is identically zero in the normal metal and constant in the superconductor. In general, due to the proximity effect, the order parameter is degraded when the interface is approached from the interface within a characteristic length. However, the discontinuity of the pair field remains at the interface and the physics discussed here is not changed qualitatively. By using the translation invariance of the system along the specular interface with spacing $`\sqrt{2}`$, the differential conductance (in units of the conductance quantum $`e^2/h`$) is calculated within the framework of the Blonder-Tinkham-Klapwijk (BTK) BTK82  theory:

$$G(E)=\frac{1}{N_y}\underset{k_y,\sigma }{}[1|r_{n,\sigma }(E,k_y)|^2+|r_{a,\sigma }(E,k_y)|^2]$$
(2)
where $`Ny\sqrt{2}`$ is the linear dimension of the junction along the $`y`$ direction, $`r_n`$ and $`r_a`$ are respectively the normal and Andreev reflection amplitude for a beam of electrons with spin $`\sigma `$ incident from the normal metal at a fixed energy $`E`$ and the transverse wavevector $`k_y=\sqrt{2}\pi l_y/N_y`$ with $`l_y=0,\mathrm{},N_y1`$. We determine the values of $`r_n`$ and $`r_a`$ by matching the wavefunction. Away from the interface, the solution to the BdG equations consists of the incoming electron, reflected and transmitted electron and hole waves. The normalized zero-temperature tunneling conductance is defined as: $`g(E)=G(E)/G_{nn}(E)`$, where $`G_{nn}(E)=\frac{1}{N_y}_{k_y,\sigma }\frac{4}{4+U_0^2/t_{}^2}`$ with $`t_{}=2t\mathrm{cos}(k_y/\sqrt{2})`$ is the tunneling conductance through a normal junction of the identical interface condition. We remark that the BTK formula can be considered as an extended version of the Landauer-Büttiker Landauer70  formula for normal metal/superconductor junctions. The latter has found spectacular success in the description of quantum coherent transport in mesoscopic normal conductors. At the same time, we exactly diagonalize Eq. (1) to calculate the LDOS for a -oriented surface of the superconductor, which is defined as:

$$\rho _{i_x}(E)=\frac{1}{N_y}\underset{n,k_y}{}[|u_{i_x,k_y}^n|^2\delta (EE_{k_y}^n)+|v_{i_x,k_y}^n|^2\delta (E+E_{k_y}^n)],$$
(3)
where for a fixed value of $`k_y`$, the summation is over the complete set of eigenstates.

Tunneling into CDW. In the case of normal-metal/CDW junction, the reflection amplitudes are obtained as:

$`r_{n,\sigma }`$ $`=`$ $`{\displaystyle \frac{\left(\frac{U_0}{t_{}}+\frac{i(EW_0)}{\sqrt{E^2W_0^2}}\right)+i\text{sgn}(E)}{\left(\frac{U_0}{t_{}}+\frac{i(EW_0)}{\sqrt{E^2W_0^2}}\right)+i\text{sgn}(E)}},`$ (4a)
$`r_{a,\sigma }`$ $`=`$ $`0.`$ (4b)
Notice that for the CDW, SDW, and DDW orderings (as considered below), the electron correlation is through the particle-hole channel so that the charge is conserved and the Andreev reflection process is absent. In Fig. 1(a), we show the tunneling conductance $`g(E)`$ as a function of energy for various strength of the tunneling barrier. The surface LDOS $`\rho (E)`$ is displayed in Fig. 1(b). As shown in Fig. 1, both the tunnel conductance and the LDOS vanish when the energy is below the energy gap $`W_0`$. As the tunnel barrier becomes more opaque, the profile of the energy dependence of the tunneling conductance resembles that of the surface LDOS. Both $`g(E)`$ and $`\rho (E)`$ differ dramatically from those of a normal-metal/$`s`$-wave superconductor junction BTK82 , by exhibiting no coherent peak at one of the gap edges. They are also asymmetric with respect to the Fermi energy $`E=0`$. We also find the property that $`g(E,W_0)=g(E,W_0)`$ and $`\rho (E,W_0)=\rho (E,W_0)`$, which show the strong dependence of both the tunneling conductance and LDOS on the macroscopic phase of the CDW.

Tunneling into SDW. In the case of normal-metal/SDW junction, we arrive at the reflection amplitudes as:

$`r_{n,\sigma }`$ $`=`$ $`{\displaystyle \frac{\left(\frac{U_0}{t_{}}+\frac{i(E\sigma S_0)}{\sqrt{E^2S_0^2}}\right)+i\text{sgn}(E)}{\left(\frac{U_0}{t_{}}+\frac{i(E\sigma S_0)}{\sqrt{E^2S_0^2}}\right)+i\text{sgn}(E)}},`$ (5a)
$`r_{a,\sigma }`$ $`=`$ $`0.`$ (5b)
In contrary to the CDW case, where the spin is degenerate for the tunnel conductance, the reflection amplitude for the SDW is spin dependent. In particular, there holds the following symmetry property of the tunneling conductance between different spin bands: $`g_\sigma (E,S_0)=g_\sigma (E,S_0)`$. It follows immediately that the total tunneling conductance is symmetric with respect to the Fermi energy, i.e., $`g(E,S_0)=g(E,S_0)`$. As shown in Fig. 2, both the tunneling conductance and the LDOS vanishes when the energy below the energy gap $`S_0`$. The profile of $`g(E)`$ also resembles to that of $`\rho (E)`$ as the strength of the tunnel barrier is increased. In addition, the coherent peak shows up at the gap edge $`E=\pm S_0`$ in both $`g(E)`$ and $`\rho (E)`$. Therefore, the overall behavior of the tunneling conductance of the N/SDW junction is qualitatively similar to the normal-metal/$`s`$-wave superconductor junction.

Tunneling into DDW. In the case of normal-metal/DDW junction, the reflection amplitudes are obtained as:

$`r_{n,\sigma }`$ $`=`$ $`{\displaystyle \frac{\left(\frac{U_0}{t_{}}+\frac{iE}{\sqrt{E^2W_𝐤^2}+W_𝐤}\right)+i\text{sgn}(E)}{\left(\frac{U_0}{t_{}}+\frac{iE}{\sqrt{E^2W_𝐤^2}+W_𝐤}\right)+i\text{sgn}(E)}},`$ (6a)
$`r_{a,\sigma }`$ $`=`$ $`0,`$ (6b)
where $`W_𝐤=4W_d\mathrm{sin}(k_y/\sqrt{2})`$. The energy gap for both the CDW and SDW ordering, as discussed above, is constant on the Fermi surface in the momentum space. The bulk energy gap of the DDW ordering is momentum-dependent, and of the form $`\mathrm{cos}k_x\mathrm{cos}k_y`$, which shows clearly that the DDW energy is closed along the nodal directions on the Fermi surface. In Fig. 3, we show the tunneling conductance of the N/DDW junction and the surface LDOS of the superconductor at the DDW state. Noticeably, a sharp peak appears at $`E=0`$ in the surface LDOS, which indicates the existence of the midgap states at the -oriented surface of the DDW ordering. However, the tunnel conductance exhibits a V-shaped feature, which persists even when the tunnel barrier is opaque. The V-shape originates from the momentum dependence of the energy gap of the DDW ordering. In particular, the tunneling conductance does not show a zero-bias anomaly, it instead dips to zero. It can also be seen clearly that when $`E=0`$, $`|r_n|^2=1`$, which gives $`g(E)=0`$, regardless of the strength of the tunnel barrier. Therefore, in the tunneling limit, the profile of the conductance is strikingly different from the surface LDOS. Physically, the energy gap of the DDW ordering comes from the particle-hole pairs. In contrast to the superconducting state, there exists no Andreev reflection process to convert, at the interface, the single particle current into the supercurrent through the condensate, so that all quasiparticles with energy below the momentum-dependent DDW gap are completely reflected. Consequently, although there exist midgap states at the -oriented surface of the superconductor in the DDW ordering, these states do not carry charge. However, as in the DSC case Zhu00 , the midgap states in the DDW state can also lead to a giant magnetic moment when a strong in-plane magnetic field is applied, the study of which is beyond the scope of current work. At this stage, we can conclude that no midgap states exist at the surface of the superconductor with either CDW or SDW ordering, but they really show up in the DDW ordering.

Tunneling into DSC. The tunneling conductance of normal-metal/DSC junctions, has been studied very intensively in the continuum model. Considering that the Andreev reflection process should take place even above $`T_c`$ in the preformed pair scenario, or in the spin-charge separation picture (where the spinons form $`d`$-wave resonant-valence-bond state), we believe it is still helpful to present, within the lattice model, the results for the DSC case. We arrive at the reflection amplitudes:

$`r_{n,\sigma }`$ $`=`$ $`{\displaystyle \frac{\left(2i\frac{U_0}{t_{}}+\frac{U_0^2}{t_{}^2}\right)|E|}{2(|E|+\sqrt{E^2|\mathrm{\Delta }_𝐤|^2})+\frac{U_0^2}{t_{}^2}|E|}},`$ (7a)
$`r_{a,\sigma }`$ $`=`$ $`{\displaystyle \frac{2\text{sgn}(E)|\mathrm{\Delta }_𝐤|}{2(|E|+\sqrt{E^2|\mathrm{\Delta }_𝐤|^2})+\frac{U_0^2}{t_{}^2}|E|}},`$ (7b)
where $`\mathrm{\Delta }_𝐤=4\mathrm{\Delta }_d\mathrm{sin}(k_y/\sqrt{2})`$. As shown in Fig. 4, the conductance at zero energy does not vanish. In particular, when the strength of the tunneling barrier is increased, a sharp peak occurs at zero energy in the conductance. In the opaque limit of the tunneling barrier, the profile of the conductance resembles that of the surface LDOS, which also exhibits a zero-energy peak, indicating the existence of midgap states. As we have emphasized above, the superconducting energy gap originates from the formation of Cooper pairs due to the electronic correlation in the particle-particle channel. Therefore, the charge is not conserved in a superconductor, and the Andreev reflection can happen in the interface between the normal metal and superconductor. It can be found from Eq. (7) that, at $`E=0`$, $`r_{n,\sigma }=0`$ and $`r_{a,\sigma }=1`$, regardless of the strength of the tunneling barrier. Therefore, the existence of the midgap states in the superconducting state corresponds to the perfect Andreev reflection process. In the phase fluctuation scenario, the long-range phase coherence is absent but the short-range correlation is still retained, the localized nature of the midgap states should be robust against the phase fluctuations. Therefore, if the phase fluctuation is the mechanism for the PG, we expect that the zero-bias conductance anomaly is still observable.

In summary, we have studied the tunneling conductance and the surface LDOS of a -oriented high-$`T_c`$ superconductor in the underdoped regime. The distinct features among the different kinds of ordering—CDW, SDW, DDW, and DSC, are investigated. We have shown, in the opaque limit of the tunneling barrier, that: (i) for the CDW ordering, asymmetric tunneling conductance and surface LDOS with vanishing intensity within the gap; (ii) for the SDW ordering, symmetric tunneling conductance and surface LDOS with vanishing intensity within the gap; (iii) for the DDW ordering, V-shaped tunneling conductance but zero-energy-peaked surface LDOS; (iv) for the DSC ordering, zero-energy peaked tunneling conductance and surface LDOS. These different features can serve to distinguish the competing scenarios for the mechanism of the PG. We have also found the existence of surface midgap states in the DDW ordering, which is, to best of our knowledge, the first example occurring in a gapped normal system. Also in this case, the correspondence between the tunneling conductance and the surface LDOS is broken down.

Acknowledgments: The author thanks S. A. Trugman for helpful discussions. This work was supported by the Department of Energy through the Los Alamos National Laboratory.