text:
# Cohesion induced deepening transition of avalanches

## I Introduction

Granular avalanches have received much attention since sandpile models are used as paradigms of so-called self-organized criticality Bak1987 . However, observations of critical-type distributions of avalanches in real physical systems are still rare, with as a notable exception the recent rice pile experiments by Frette et al. Frette1996 . It was suggested by Christensen et al. Christensen1996  that the anisotropy in the rice grains allows more stable packing configurations in a granular pile, and that this could be responsible for the successful observation of criticality. Some of the recent attention has been drawn to avalanches in cohesive granular materials with the premise that cohesion will also allow the sand more packing configurations and thus increase the likelihood of observing critical scaling behavior. While the goal of finding criticality in cohesive sandpiles remains to be fulfilled even after the experimental work by Quintanilla et al. Quintanilla2001 , the effect of cohesion in granular avalanches represents an interesting direction for a theoretical study.

In this article, we’ll use the discrete-height version of the sandbox (DHSB) model introduced in Ref. Chen2002  for an unloading sandbox (Fig. 1) to understand the effects of cohesion in directed avalanche systems.

In the following section, we’ll discuss how we can model cohesiveness in avalanche systems. In Sec. III, we’ll review the DHSB model and introduce a cohesion parameter. Previous results in Refs. Chen2002 ; Chen2002a  represent a special case of the model where the system is in the deep avalanche phase with the cohesion parameter $`p=1/2`$. In Sec. IV, we describe the step-flow random-deposition (SFRD) interface growth model which underlies the DHSB model and the directed percolation (DP) roughening transition of the SFRD model. In Sec. V, we focus on the two deterministic limits of the model and present the exact solution in one of these limits. In Sec. VI, numerical results for the avalanches in the flat phase of the interface model are presented. In Sec. VII, we investigate the scaling behavior at the transition point where the interface roughness increases logarithmically in time. We show that the avalanche-scarred sand surface, while being rougher than nonscarred ones, retains the same scaling exponent of the roughness in the thermodynamic limit. However, we’ll also show that at the transition point, the violation of mass hyperscaling relation spoils the reduction to two independent exponents established in Ref. Chen2002a . We’ll summarize our results in Sec. VIII.

## II Tunable parameter for cohesion

One interesting characters of cohesion in sand is that it possesses hysteresis behavior. Consider building a sand castle on a beach. It’s common sense that we’ll need to add water to the sand before we can shape it into a standing castle. However, without disturbance, the sand castle can somehow maintain its shape even after it dries out Hornbaker1997 . The moisture in sand increases the cohesion between the sand particles Nase2001  and allows one to manipulate the sand into a stable shape that, while not as attainable, is more or less an equally valid stable shape for dry sand.

In accounting for this standing-sand-castle effect, we’ll use the same stability condition for all cohesiveness of the sandbox. While, in reality, the space of possible stable configurations for wet and dry sand should not be exactly identical, in this article, we shall ignore this distinction to avoid complicating the rules too much.

On the other hand, the way an unstable sand surface topples surely depends on the cohesiveness. In the DHSB model discussed below, there are only two possible final stable states for any toppling site. We’ll call them the minimal stable state and the maximal stable state. These two states are similar to the angle of repose and maximal stable angle in a real sandpile. However, in sandbox model, these states are microscopic while the “angles” of a real sandpile are macroscopic. We’ll use a parameter $`p`$, which is a real number between $`0`$ and $`1`$, to represent the strength of cohesion. In the model, $`p`$ is the probability for a toppling site of the sandpile to settle into the maximal stable state instead of the minimal one. For wet sand, the $`p`$ is large, and for dry sand, the $`p`$ is small.

## III Discrete-height sandbox model

With the discussion of the previous section in mind, let’s review the dynamic rules of the discrete-height sandbox model. The surface of a sandbox (see Fig. 1) is represented by an integer height variable $`h`$ defined on a two-dimensional square lattice which is tilted at $`45^{}`$ with respect to the lowering wall as illustrated in Fig. 2.

This is equivalent to considering only the lattice points whose integer $`x`$ and $`y`$ coordinates satisfy the condition that $`x+y`$ is an even number. The lowering wall that drives the system by creating unstable sites is located at the $`y=0`$ row and the activities in the system propagate only in the positive $`y`$ direction. In our numerical simulations, the system is periodic in the $`x`$ direction, which is parallel to the driving wall. The sizes of the system in the $`x`$ and $`y`$ directions are denoted by the numbers of sites $`L_x`$ in each row and the number of rows $`L_y`$ respectively.

As in most sandpile processes, the dynamics of the sandbox model is defined by a stability condition, a toppling rule, and a driving method. They are as follows. The stability condition of the DHSB is given by

$$h(x,y)\mathrm{min}[h(x1,y1),h(x+1,y1)]+s_c$$
(1)
with $`s_c=1`$, which represents the local maximal stable slope. The unstable sites in the system topple with the rule

$$h(x,y)\mathrm{min}[h(x1,y1),h(x+1,y1)]+\eta ,$$
(2)
where $`\eta =0`$ with probability $`1p`$ and $`\eta =1`$ with probability $`p`$. (In the earlier studies Chen2002 ; Chen2002a , the value of $`p`$ is always $`1/2`$ .) This is the only place in the dynamics of the DHSB that the cohesion parameter $`p`$ comes into play. The lowing wall which drives the system is implemented in the model by randomly picking one of the highest sites $`(x_i,0)`$ on the $`y=0`$ row and by reducing its height by $`1`$:

$$h(x_i,0)h(x_i,0)1,$$
(3)
where $`i`$ is the Monte Carlo time, which also serves as an age index for the avalanches.

A typical configuration of the DHSB before and after an avalanche is shown in Fig. 3. Since the toppling of a site on a given row $`y`$ only affects the stability of the two sites immediately above it at the $`y+1`$ row, we choose to update the system in a row-by-row fashion. For each avalanche, the entire system is stabilized by such a single sweep of topplings from $`y=0`$ to $`y=L_y`$.

## IV Underlying interface dynamics

The underlying interface dynamics of the sandbox models is given by the step-flow random-deposition (SFRD) models with a two-step growth rule Chen2002 ; Chen2002a  as illustrated in Fig. 4.

The mapping between the sandbox system and the interface growth model involves identifying the $`y`$ coordinate of the sandbox model with the time $`t`$ of the interface growth. Each stable sandbox surface thus can be viewed as a space-time world-sheet configuration of the interface growth. Models similar to this generally belong to the Kardar-Parisi-Zhang (KPZ) universality class Kardar1986  with the critical exponents $`\alpha =1/2`$, $`\beta =1/3`$, and $`z=\alpha /\beta =3/2`$ which characterize the scaling of interface roughness

$$W^2\overline{(h\overline{h})^2}.$$
(4)
Starting from a flat interface at $`t=0`$, the interface grows rougher with

$$Wt^\beta .$$
(5)
And, after a characteristic time $`t_cL^z`$, the roughness will saturate with a value

$$WL^\alpha $$
(6)
depending on the system size $`L`$.

From the mapping introduced in Ref. Chen2002 , the avalanche exponents are given by

$$\tau _l=\frac{\sigma 1\alpha }{z}=2,$$
(7)
$$\tau _w=\sigma z\alpha =\frac{5}{2},$$
(8)
and

$$\tau _\delta =\frac{\sigma 1z}{\alpha }=4$$
(9)
for the distribution functions, $`P_l(l)l^{\tau _l}`$, $`P_w(w)w^{\tau _w}`$, and $`P_\delta (\delta )\delta ^{\tau _\delta }`$, of avalanche length $`l`$, width $`w`$, and depth $`\delta `$. As defined in Ref. Chen2002 , the avalanche length $`l`$ (width $`w`$) represents maximum $`y`$ ($`x`$) distance of the toppling sites from the triggering point while the avalanche depth $`\delta `$ is the maximum height change of the toppling sites. The $`\sigma `$ in these expressions was eliminated with the mass hyperscaling relation

$$\sigma =2+z+2\alpha .$$
(10)
obtained from the compactness of the avalanche clusters, i.e., assuming $`mlw\delta `$.

However, the discrete-height version of the SFRD model undergoes a DP roughening transition at $`p=p_c0.294515`$ similar to those studied by Kertész and Wolf Kertesz1989  also Alon et al. Alon1996 . The KPZ scaling behavior only applies when the value of the control parameter $`p`$ is greater than the critical value $`p_c`$. Below this transition point the interface is in a trivial flat state, where, for a stationary interface (interface time $`y\mathrm{}`$), the density of sites at the bottom $`h=h_0`$ layer is finite. The interface is thus *pinned* at this level and its growth rate becomes zero.

At the transition point $`p=p_c`$, we find the roughness of the SFRD interface diverges only logarithmically in time

$$W^2(\mathrm{ln}t)^\gamma ,$$
(11)
with the exponent $`\gamma 1`$ similar to that of the Kertész and Wolf’s model as well as the restricted version of the models by Alon et al..

## V Deterministic limits

In the two limits, $`p=1`$ and $`p=0`$, the toppling process of the avalanches becomes deterministic and the sand surface topples down layer by layer. The only randomness in the process comes from the driving method (3), i.e., that we randomly lower one of the highest sites at the $`y=0`$ row to trigger an avalanche. The typical avalanche scar configurations at these two limits are shown in Fig. 5. These are the edges of avalanche clusters left on the surface, some of which are partially erased by newer avalanches.

### V.1 Domain walls at $`p=0`$

The $`p=0`$ limit runs into the complication that in the bulk of the system ($`y>0`$) the sand surface goes down by $`2`$ units at a time. Since $`\mathrm{\Delta }hh(x,y)\mathrm{min}[h(x1,y1),h(x+1,y1)]=1`$ is stable according to the stability condition (1), and the sites on the $`y=0`$ row always goes down by $`1`$ unit each time according to the driving method (3), the sites on the $`y=1`$ row will only topple when their heights are $`2`$ units higher than the triggering sites and they alway go down by $`2`$ units to the same height of the triggering site according to the toppling rule (2). All the sites at higher rows will be locked into the same even-oddness as the sites triggering their toppling. Therefore, after all sites have participated in at least one avalanche, their even-oddness will be fixed for all subsequent topplings. This means the even-oddness of a site is preserved by the toppling process, and that the lines separating the even and odd sites thus form impenetrable domain walls for the avalanches (see Fig. 6). This hinders the applicability of the same type of analysis as presented below for the $`p=1`$ limit.

However, the numerical results in Sec. VI will show that the same scaling exponents as those of $`p=1`$ case control this limit, too.

### V.2 Exact solution at $`p=1`$

The $`p=1`$ limit has a nice solution. Since the sites in the bulk topple from $`\mathrm{\Delta }h=2`$ to $`\mathrm{\Delta }h=1`$, the sand surface indeed goes down only one layer at a time without the complications as the $`p=0`$ case. An exact solution can be obtained by considering the avalanches taking place in such one single layer. For a brand-new layer, the two boundaries of the first avalanche open up linearly until the avalanche spans the system in the $`x`$ direction and leaves two scar lines on the surface. The two boundaries of the second avalanche expand until they meet the scar lines created by the first avalanche. Then, they turn and follow those scar lines until they meet with each other and terminate the avalanche. Subsequent avalanches follow the same scenario. The maximum distance an avalanche cluster can expand from its triggering point to each side in the $`x`$ direction is exactly half the distance from the nearest triggering point of the previous avalanches in the same layer on that side. As the triggering points are chosen in an uncorrelated manner, the maximum width $`w`$ of an avalanche should follow the Poisson distribution

$$P_w(w)=\frac{\mu ^we^\mu }{w!}$$
(12)
if $`\mu `$ is the average distance between the triggering points of the previous avalanches in the same layer in the stationary state. The avalanche under consideration could be any one of the avalanches happening in the same layer. Thus, we need to average over the number of avalanches $`n`$ taking place before this one in the same layer. For a system of transverse size $`L_x`$, $`n=L_x/\mu `$, the integral can be carried out explicitly and gives

$$_0^{\mathrm{}}\frac{\mu ^we^\mu }{w!}𝑑\frac{1}{\mu }=\frac{(w2)!}{w!}w^2,$$
(13)
which results in

$$\tau _l=\tau _w=2.$$
(14)
The same results can also be derived from Eq. (7) and (8) by assuming $`z=1`$ and $`\alpha =0`$. Since the avalanches are compact, the hyperscaling relation (10) and other exponent relations (7)–(9) from Ref. Chen2002  hold.

## VI Shallow-avalanche phase

Below the transition point, the underlying interface model is in a flat phase where the bottom layer percolates with finite density. All the information of the initial configuration of the interface (the $`y=0`$ row next to the wall) is wiped out at a time scale proportional to the sizes of the islands higher than the bottom layer in the initial state. (Without deposition, the sizes of these islands decrease linearly in time.) While the underlying interface model is in a trivial phase, much like the uncorrelated stationary state in Dhar and Ramaswamy’s directed sandpile model Dhar1989 , the avalanche distributions of the system may still exhibit power-law scaling. The numerical values of the scaling exponents shown in Fig. 7 confirm the power-law scaling of the distributions and they are similar to those values found at the $`p=1`$ fixed point following $`z=1`$ and $`\alpha =0`$.

While an exact solution is not available in this phase, we can understand the scaling exponent $`z=1`$ from the perspective that the DP clusters triggered from single seeds in the percolating phase open up linearly $`lw`$; and also that roughness exponent $`\alpha =0`$ comes from that the interface is flat. However, a difference is that while $`p<p_c`$ represents an entire phase of shallow avalanches which should be controlled by an attractive fixed point, the $`p=1`$ fixed point is unstable in the sense that the scaling behavior falls back to the KPZ universality class for any small deficiency in the cohesiveness $`p`$ from the value $`1`$.

## VII DP roughening transition

At the transition point $`p=p_c`$, the interface roughness diverges logarithmically thus the $`\beta `$ and $`\alpha `$ exponents, defined by Eqs. (5) and (6), are both zero. Nonetheless, the dynamic exponent $`z`$ has a nontrivial value $`z_{\mathrm{DP}}1.582`$ originating from the DP nature of the bottom-layer dynamics. Moreover, at the transition point, the avalanche clusters lose their compact shapes (see Fig. 8) and we should not expect the exponent relations (7)–(10), nor the calculation in Ref. Chen2002a  for the corrections to scaling to remain valid. In this section we will demonstrate the break down of mass hyperscaling relation (10) and how the avalanches affect the roughness of the sand surface.

### VII.1 Breakdown of mass hyperscaling

At the transition point, the bottom layer of an avalanche cluster follows the critical DP dynamics. Therefore, we should expect from the fractal DP cluster shape that the density of sites at the lowest $`h=h_0`$ level goes to zero in the thermodynamic limit for large avalanches. However, the overall shape of an avalanche consists, in addition, of sites at $`h_0+1,h_0+2,`$ …levels. The higher-level sites that participate in the avalanche fill into the holes and voids next to the bottom layer cluster and more or less bring the avalanche cluster back to a compact shape. We can verify this compactness of the avalanche cluster by a direct measurement of the ratio $`a/(lw)`$, with $`a`$ being the area of (or, the number of sites participating in) an avalanche. The result is shown as the solid line in Fig. 9.

The approach to a finite value on the vertical axis demonstrates the compactness of the avalanche clusters by the existence of a finite area density $`0.2`$ in the thermodynamic limit. The FSS estimates are plotted against $`1/\mathrm{ln}y`$ instead of $`1/y`$ since the roughness of the surface diverges only logarithmically in $`y`$, which will be elaborated later.

Contrary to a finite area density, as also shown in Fig. 9, the mass density $`m/(lw\delta )`$ (the dashed line) goes to zero in the thermodynamic limit. The absence of a finite mass density breaks the scaling

$$mlwd,$$
(15)
which lead, in Ref. Chen2002 , to the mass hyperscaling relation (10). The plot of the combined exponent $`\sigma z2\alpha `$ in Fig. 10 shows the violation of Eq. (10) as the FSS estimates approach $`1.72`$ which is much lower than the expected value $`2`$ for compact avalanches obeying Eq. (15).

Also shown in Fig. 10 are the plots for the $`\alpha `$, $`z`$, and $`\sigma `$ exponents. They are consistent with $`z=z_{\mathrm{DP}}`$ and more or less with $`\alpha =0`$. This confirms that the scaling behavior of the avalanches follows those of the SFRD interface. The slow convergence of $`\alpha `$ is to be expected from the logarithmic divergence of the interface roughness.

### VII.2 Interface roughness

The remaining question is how the scaling behavior of the roughness is changed by the iterated avalanche process. We approach this by looking at the change of the global surface roughness itself and by comparing the scaling of this change to the scaling of the original interface roughness. The same analysis was performed in Ref. Chen2002a  which concerns only the $`p=1/2`$ case of the DHSB, and it was found that the change in the global roughness by the avalanche process only represents large corrections to the KPZ scaling behavior of the surface. However, at the DP transition point, the interface roughness diverges only logarithmically. This makes the scaling of interface roughness more likely to be overwhelmed by the change in the roughness due to the avalanche process, and we generally would not expect the values of the scaling exponents to remain the same. In the following, we’ll show the scaling of the interface does follow the same logarithmic divergence.

We perform a direct measurement of the global interface roughness at the transition point. The results are shown in Fig. 11(a).

As in the $`p=1/2`$ case, the surface is made rougher by the iterated avalanches. The increase in the roughness $`\mathrm{\Delta }W^2`$ scale as $`(\mathrm{ln}t)^{\gamma _\mathrm{\Delta }}`$ with the exponent $`\gamma _\mathrm{\Delta }0.4`$ which is shown as the dashed line in Fig. 11(b). Since the interface roughness itself scales as $`W^2(\mathrm{ln}t)^\gamma `$ with $`\gamma 1`$ which is shown as the solid line in Fig. 11(b), the change in the roughness is irrelevant comparing to the interface scaling. We can thus conclude that in the thermodynamics limit, the stationary surfaces of DHSB have the same $`\gamma `$ exponent as the SFRD interfaces. Just as the in the deep phase ($`p=1/2`$) of the the avalanche, the iterated avalanche process only gives rise to sizable corrections to the interface scaling behavior.

## VIII Summary

In this article, we introduced the DHSB as a model for avalanches in granular materials with variable cohesiveness. This model exhibits a deepening transition from a shallow-avalanche phase where avalanches only involve a couple of surface layers of the granular material, into a deep-avalanche phase where the depths of avalanches increase as power laws in their lengths or widths. In the deep-avalanche phase, the scaling behavior of the avalanches belongs to the KPZ universality class: The avalanche clusters scale anisotropically with $`lw^{3/2}`$ and depth increase as $`\delta w^{1/2}`$. In the flat phase, the avalanche clusters scale isotropically $`lw`$ with finite depths.

In both phases, the mass hyperscaling relation (10) based on compactness (15) of the avalanches holds. On the other hand, at the transition point, the hierarchical DP structure, pointed out by Täuber et al. Tauber1998 , for each height level breaks this scaling in a subtle way. While the mass density $`m/(lw\delta )`$ of the avalanche clusters goes to zero in the thermodynamic limit, the area density $`a/(lw)`$ remains finite. However, the exact scaling behavior of the systems at this DP roughening transition point remains unclear even without the iterated avalanche in the DHSB model Alon1998 ; Lopez1998 ; Goldschmidt1999 .

While we are not aware of any experimental study on how the avalanche behavior of a system will vary with a gradual change in the cohesiveness of the grains, the cohesiveness in granular system is known to vary with moisture Nase2001  and grain sizes Valverde2000 ; Quintanilla2001 . We thus expect experimental studies in this direction to be feasible. The DHSB model represents a system with a layered structure where the heights are discrete, and the DP nature of the deepening transition relies heavily on a well-defined bottom layer or minimal stable configuration of the system. It thus wouldn’t be a surprise if exact DP scaling were not to be observed in the avalanches of most experimental sandpiles. Nonetheless, the breakdown of the mass hyperscaling relation (10) comes from the fractal aspect of the hierarchical DP clusters and is a more fundamental property. It would serve as a hallmark of such a transition if it’s to be observed experimentally.

The author thanks Marcel den Nijs for many helpful discussions and critical reading of the manuscript. This research is supported by the National Science Foundation under Grant No. DMR-9985806.


format: text

text:
# Rabi oscillations, Ramsey fringes and spin echoes in an electrical circuit

## Abstract

We present a superconducting tunnel junction circuit which behaves as a controllable atom, and whose ground and first excited state form an effective spin 1/2. By applying microwave pulses, we have performed on this circuit experiments demonstrating the controlled manipulation of the spin : Rabi precession, Ramsey interferences, and spin echoes.

The state variables of an electrical circuit, like voltages and currents can be made to behave quantum mechanically by minimizing the coupling to external degrees of freedom through a proper design. Circuits based on superconducting tunnel junctions have displayed signatures of macroscopic quantum behavior , but the level of coherence of the quantum states remained until now much smaller than for isolated atoms or ions. In the “quantronium” circuit presented here, a coherence quality factor of more than $`10^4`$ has been obtained, thus allowing the coherent manipulation of the state of the system like in atomic physics and NMR experiments .

The quantronium consists of a superconducting loop interrupted by two adjacent small Josephson tunnel junctions with capacitance $`C_j/2`$ and Josephson energy $`E_J/2`$ each, which define a low capacitance superconducting electrode called the “island”, and by a large Josephson junction with large Josephson energy $`(E_{J0}20E_J)`$ (see Fig. 1). The island is biased by a voltage source $`U`$ through a gate capacitance $`C_g`$. In addition to $`E_J`$, the quantronium has a second energy scale which is the Cooper pair Coulomb energy $`E_{CP}=(2e)^2/2\left(C_g+C_j\right)`$. The temperature $`T`$ and the superconducting gap $`\mathrm{\Delta }`$ satisfy $`k_BT\mathrm{\Delta }/\mathrm{ln}𝒩`$ and $`E_{CP}<\mathrm{\Delta }k_BT\mathrm{ln}𝒩`$, where $`𝒩`$ is the total number of paired electrons in the island. The number of excess electrons is then even , and the system has discrete quantum states which are in general quantum superpositions of several charge states with different number $`\widehat{N}`$ of excess Cooper pairs in the island. Neglecting the loop inductance and the charging energy of the large junction, the Hamiltonian of the circuit is

$$\widehat{H}=E_{CP}\left(\widehat{N}N_g\right)^2E_J\mathrm{cos}(\frac{\widehat{\gamma }+\varphi }{2})\mathrm{cos}\widehat{\theta }E_{J0}\mathrm{cos}\widehat{\gamma }I_b\phi _0\widehat{\gamma }\text{,}$$
(1)
where $`N_g=C_gU/2e`$ is the dimensionless gate charge, $`\varphi =\mathrm{\Phi }/\phi _0`$ is a phase bias, with $`\mathrm{\Phi }`$ the external flux imposed through the loop and $`\phi _0=\mathrm{}/2e`$, $`\widehat{\gamma }`$ is the phase across the large junction, and $`\widehat{\theta }`$ is the phase operator conjugate to the Cooper pair number $`\widehat{N}`$. The bias current $`I_b`$ is zero except during readout of the state .

In our experiment $`E_JE_{CP}`$ and neither $`\widehat{N}`$ nor $`\widehat{\theta }`$ is a good quantum number. In contrast, the large junction is shunted by a large capacitor $`C`$ so that $`\widehat{\gamma }`$ is almost a classical variable. In this regime, the energy spectrum is sufficiently anharmonic for the two lowest energy states $`|0`$ and $`|1`$ to form a two-level system. This system corresponds to an effective spin one-half $`\stackrel{}{s}`$ with eigenstates $`|0|s_z=+1/2`$ and $`|1|s_z=1/2.`$ At $`N_g=1/2,`$ $`I_b=0`$ and $`\varphi =0,`$ its “Zeeman energy” $`h\nu _{01}`$, of the order of $`E_J,`$ is stationary with respect to $`N_g`$, $`I_b`$ and $`\varphi `$ , making the system immune to first order fluctuations of the control parameters. Manipulation of the quantum state is thus performed at this optimal point by applying microwave pulses $`u(t)`$ with frequency $`\nu \nu _{01}`$ to the gate, and any superposition $`|\mathrm{\Psi }=\alpha |0+\beta |1`$ can be prepared, starting from $`|0`$. In a frame rotating around the quantization axis $`z`$ at frequency $`\nu ,`$ the microwave voltage acts on $`\stackrel{}{s}`$ as an effective dc magnetic field in the $`xz`$ plane, with $`x`$ and $`z`$ components proportional to the microwave amplitude and to the detuning $`\mathrm{\Delta }\nu =\nu \nu _{01},`$ respectively.

For readout, we have implemented a strategy reminiscent of the Stern and Gerlach experiment , in which the information about the spin of silver atoms is transferred onto their transverse position. In our experiment, the information on $`\stackrel{}{s}`$ is transferred onto the phase $`\widehat{\gamma },`$ and the two states are discriminated through the supercurrent in the loop $`\widehat{I}=\widehat{H}/\widehat{\delta }/\phi _0,`$ where $`\widehat{\delta }=\widehat{\gamma }+\varphi `$ is the phase difference across the series combination of the two small junctions . For this purpose, a trapezoidal readout pulse $`I_b(t)`$ with a peak value slightly below the critical current $`I_0=E_{J0}/\phi _0`$ is applied to the circuit (see Fig. 2).

When starting from $`\widehat{\delta }0`$, the phases $`\widehat{\gamma }`$ and $`\widehat{\delta }`$ grow during the current pulse, and consequently an $`\stackrel{}{s}`$-dependent supercurrent develops in the loop. The loop current is the analog of the transverse acceleration experienced by a silver atom in the magnetic field gradient. This current adds to the bias-current in the large junction, and by precisely adjusting the amplitude and duration of the $`I_b(t)`$ pulse, the large junction switches during the pulse to a finite voltage state with a large probability $`p_1`$ for state $`|1`$ and with a small probability $`p_0`$ for state $`|0`$. A switching event corresponds to the impact of a silver atom on the top spot of the screen. The absence of switching corresponds to the impact of a silver atom on the bottom spot of the screen. For the parameters of the experiment, the efficiency of this projective measurement should be $`\eta =p_1p_0=0.95`$ for optimum readout conditions. Large ratios $`E_{J0}/E_J`$ and $`C/C_j`$ provide further protection from the environment.

An actual “quantronium” sample is shown on the right side of Fig. 1. It was fabricated with standard technique of aluminum evaporation through a shadow-mask obtained by e-beam lithography. With an external microwave capacitor $`C=1\mathrm{p}F`$, the plasma frequency of the large junction with $`I_0=0.77\mu A`$ is $`\omega _p/2\pi 8\mathrm{G}Hz`$. The sample and last filtering stage were anchored to the mixing chamber of a dilution refrigerator with 15 mK base temperature. The switching of the large junction to the voltage state is detected by measuring the voltage across it with a room temperature preamplifier followed by a discriminator with a threshold voltage $`V_{th}`$ well above the noise level (Fig. 2). By repeating the experiment (typically a few 104 times), we can determine the switching probability, and hence, the occupation probabilities $`\left|\alpha \right|^2`$ and $`\left|\beta \right|^2`$.

The readout part of the circuit was tested by measuring the switching probability $`p`$ as a function of the pulse height $`I_p`$, for a current pulse duration of $`\tau _r=100\mathrm{n}s`$, at thermal equilibrium . The discrimination between the currents corresponding to the $`|0`$ and $`|1`$ states was found to have an efficiency of $`\eta =0.6`$, which is lower than the expected $`\eta =0.95`$. Measurements of the switching probability as a function of temperature and repetition rate indicate that the discrepancy between the theoretical and experimental readout efficiency could be due to an incomplete thermalization of our last filtering stage in the bias current line.

Spectroscopic measurements of $`\nu _{01}`$ were performed by applying to the gate a weak continuous microwave irradiation suppressed just before the readout current pulse. The variations of the switching probability as a function of the irradiation frequency display a resonance whose center frequency evolves as a function of the dc gate voltage and flux as the Hamiltonian (1) predicts, reaching $`\nu _{01}16.5\mathrm{G}Hz`$ at the optimal working point (see Fig. 3). The small discrepancy between theoretical and experimental values of the transition frequency at nonzero magnetic flux is attributed to flux penetration in the small junctions not taken into account in the model. We have used these spectroscopic data to precisely determine the relevant circuit parameters and found $`i_0=18.1\mathrm{n}A`$ and $`E_J/E_{CP}=1.27`$. The linewidth $`\mathrm{\Delta }\nu _{01}`$ is given in the bottom panels in Fig. 3 as a function of $`\varphi `$ and $`N_g`$. At the optimal working point, the linewidth was found to be minimal with a $`0.8`$MHz full width at half-maximum, corresponding to a quality factor $`Q=2\times 10^4`$. The lineshape was found to be irreproducible, probably because of slight shifts of the resonance frequency during the measurement, related to low frequency charge or noise on the phase $`\delta `$. If one considers the narrowest lines recorded, the linewidth varies linearly when departing from the optimal point $`(N_g=1/2,\varphi =0,I_b=0)`$, the proportionality coefficients being $`\mathrm{\Delta }\nu _{01}/N_g250\mathrm{M}Hz`$ and $`\mathrm{\Delta }\nu _{01}/(\varphi /2\pi )430\mathrm{M}Hz`$. These values can be translated into RMS deviations $`\mathrm{\Delta }N_g=0.004`$ and $`\mathrm{\Delta }(\delta /2\pi )=0.002`$ during the time needed to record the resonance. The residual linewidth at the optimal working point can be explained by the second order contribution of these noises. The amplitude of the charge noise is in agreement with measurements of $`1/f`$ charge noise , and its effect could be minimized by increasing the $`E_J/E_C`$ ratio. By contrast, the amplitude of the phase noise corresponds to a large flux noise , but it could be also attributed to bias current noise.

When varying the delay between the end of a resonant irradiation and the measurement pulse at the optimal working point, the switching probability decays with a time constant $`T_1=1.8\mu s`$ (see Fig. 4). Supposing that the energy relaxation of the system is only due to the bias circuitry, a calculation along the lines of Ref.  predicts that $`T_110\mu s`$ for a crude discrete element model. This result shows that no detrimental sources of dissipation have been seriously overlooked in our circuit design.

We have then performed controlled rotations of $`\stackrel{}{s}`$ with large amplitude microwave pulses. Prior to readout, a single pulse at the transition frequency with variable amplitude $`U_{\mu w}`$ and duration $`\tau `$ was applied. The resulting change in switching probability is an oscillatory function of the product $`U_{\mu w}\tau `$ (see Fig. 5), in agreement with the theory of Rabi oscillations . It provides direct evidence that the resonance indeed corresponds to an effective spin rather than to a spurious harmonic oscillator resonance in the circuit. The proportionality ratio between the Rabi period and $`U_{\mu w}\tau `$ was used to calibrate microwave pulses for the application of controlled rotations of $`\stackrel{}{s}`$.

The measurement of the coherence time of $`\stackrel{}{s}`$ during free evolution was obtained by performing a Ramsey-fringes-like experiment . One applies on the gate two phase coherent microwave pulses corresponding each to a $`\pi /2`$ rotation around $`x`$  and separated by a delay $`\mathrm{\Delta }t`$ during which the spin precesses freely around $`z`$. For a given detuning $`\mathrm{\Delta }\nu `$ of the microwave frequency, the switching probability displays decaying oscillations of frequency $`\mathrm{\Delta }\nu `$ (see Fig. 6), which correspond to the “beating” of the spin precession with the external microwave field. The envelope of the oscillations yields the coherence time $`T_\phi 0.5\mu s`$. Given the transition period $`1/\nu _{01}60\mathrm{p}s`$, this means that $`\stackrel{}{s}`$ can perform on average 8000 coherent free precession turns.

When the circuit is biased away from the optimal point, the coherence time $`T_\phi `$ of the oscillation is strongly reduced, as shown in the top panel of Fig. 7 for $`N_g=1/2+0.02`$, $`\varphi =0`$. In order to determine the contribution to dephasing of low frequency charge noise, we have performed spin echo experiments: an intermediate $`\pi `$ pulse is inserted between the two $`\pi /2`$ pulses of the Ramsey sequence (see Fig. 7, right side). The effect of the $`\pi `$ pulse is to make the phases accumulated during the two free evolution time intervals $`t_1`$ and $`t_2`$ to subtract one from the other (see Fig. 7). By symmetry, when $`t_1=t_2,`$ the total phase accumulated is independent of $`\mathrm{\Delta }\nu `$ if $`\mathrm{\Delta }\nu `$ is constant over the complete sequence. Compared to the Ramsey fringes experiment, where coherence during $`\mathrm{\Delta }t`$ is revealed by the periodic evolution of $`\mathrm{cos}\left[2\pi \mathrm{\Delta }\nu \mathrm{\Delta }t\right],`$ the echo signal varies as $`(1\mathrm{cos}\left[2\pi \mathrm{\Delta }\nu (t_2t_1)\right])/2`$ and is therefore less sensitive to fluctuations of $`\mathrm{\Delta }\nu `$ from sequence to sequence when $`t_1t_2`$. In the experiment, we have recorded the switching probability at fixed values of $`\mathrm{\Delta }t,`$ as a function of the delay $`t_1`$ (left panels of Fig. 7). Up to $`\mathrm{\Delta }t1`$ $`\mu `$s, fringes emerge around $`t_1=t_2=(\mathrm{\Delta }t\tau )/2`$ (here, $`\tau 15\mathrm{n}s`$), indicating that during pulse sequences of this duration, coherence was at least partly conserved. As expected, the period of the oscillations is twice as short in the echo experiment than in the Ramsey experiment. The observation of spin echoes at time scales much larger than the decay time of the Ramsey fringes indicates that in this situation decoherence was essentially due to charge fluctuations at frequencies lower than $`1/\mathrm{\Delta }t1\mathrm{M}Hz`$. No echo was seen in experiments performed at $`\varphi 0`$, suggesting that the relevant phase noise was at higher frequencies.

In all our time domain experiments, the oscillation period of the switching probability closely agrees with theory, meaning a precise control of the preparation of $`\stackrel{}{s}`$ and of its evolution. However, the amplitude of the oscillations is smaller than expected by a factor of three to four. This loss of contrast is likely to be due to a relaxation of the level population during the measurement itself. In principle the current pulse, whose rise time is 50 ns, is sufficiently adiabatic not to induce transitions directly between the two levels. Nevertheless, it is possible that the readout or even the preparation pulses excite resonances in the bias circuitry which in turn could induce transitions in our two-level manifold. Experiments using better shaped readout pulses and a bias circuitry with better controlled high-frequency impedance are needed to clarify this point.

In conclusion we have designed and operated a superconducting tunnel junction circuit which behaves as a tunable two-level atom that can be decoupled from its environment. When the readout is off, the coherence of this “quantronium” atom is of sufficient quality ($`Q_\phi =2.5\times 10^4`$) that an arbitrary quantum evolution can be programmed with a series of microwaves pulses. Coupling several of these circuits can be achieved using on-chip capacitors. The ability to tune and address them individually would allow to produce entangled states and probe their quantum correlations. These fundamental physics experiments could lead to the realization of quantum logic gates, an important step towards the practical implementation of solid-state quantum processors .

Acknowledgements: The indispensable technical work of Pief Orfila is gratefully acknowledged. This work has greatly benefited from direct inputs from J. M. Martinis and Y. Nakamura. The authors acknowledge discussions with P. Delsing, G. Falci, D. Haviland, H. Mooij, R. Schoelkopf, G. Schön and G. Wendin. This work is partly supported by the European Union through contract IST-10673 SQUBIT and by the Conseil Général de l’Essonne (EQUM project).


text:
# SMC X-1 As An Intermediate-Stage Flaring X-ray Pulsar

## 1 Introduction

Neutron star X-ray binaries are generally categorized into two groups: low mass X-ray binaries (LMXBs) and X-ray pulsars. The surface magnetic field of the central neutron star in an LMXB is thought to be $``$108 G. The mass accretion with this magnetic field is most likely spherical, so that no significant inhomogeneity in the X-ray emission over the neutron star surface is expected – i.e., no persistent coherent pulsations are observed. The strong magnetic field, $``$1012 G, of an X-ray pulsar, on the other hand, can funnel the accretion matter onto the magnetic pole, which makes the central neutron star appear as a pulsar.

Of particular interest are the so-called “intermediate-stage sources” speculated to lie between the LMXBs and X-ray pulsars, including “the Rapid Burster” (MXB 1730–355), GRO J1744–28 (“the Bursting Pulsar”), and SAX J1808.4–3658 (“the accreting millisecond pulsar”) (Lewin et al. 1976; Fishman et al. 1995; in ’t Zand et al. 1998). Both the Rapid Burster and GRO J1744–28 exhibit Type II X-ray bursts; however, only the former shows Type I X-ray bursts while only the latter has apparent coherent pulsations (Lewin et al. 1996). SAX J1808.4–3658, on the other hand, shows both Type I bursts and coherent pulsations (Wijnands & van der Klis 1998; Chakarbarty & Morgan 1998), but not Type II bursts. The magnetic field strengths of these sources have been inferred to be $``$108-11 G, between those of LMXBs and X-ray pulsars.

Another possible intermediate-stage source is the X-ray pulsar SMC X-1, which has similar properties to GRO J1744–28, including its fast spin period ($``$0.72 s for SMC X-1; $``$0.47 s for GRO J1744–28), steady spin-up, and inferred magnetic field ($``$1011 G ) (Bildsten & Brown 1997; Li & van den Heuvel 1997). In addition, once SMC X-1 was observed with an X-ray burst that resembles Type II bursts (Angelini, Stella, & White 1991). It may be possible, therefore, that SMC X-1 and GRO J1744–28 belong to a distinctive group of X-ray binaries, “bursting pulsars”, which show both coherent pulsations and Type II X-ray bursts (Li & van den Heuvel 1997). To investigate this important possibility, we analyze all publicly available RXTE data for SMC X-1, searching for phenomena that may be related to X-ray bursts. We report that SMC X-1 in fact exhibits active flares resembling Type II bursts from GRO J1744–28.

## 2 Data Analysis and a Flare Search

We analyzed all publicly available RXTE Proportional Counter Array observations toward SMC X-1. Photon arrival times from the Good Xenon mode were transformed to the solar system barycenter using the JPL DE400 ephemeris. The Very Large Event Models were used to subtract backgrounds, and only the Standard 2 data obtained from the top xenon layers of Proportional Counter Units 0, 1, and 4 were considered for spectral analysis. Both the data within 30 minutes after passages through South Atlantic Anomaly and/or with a high ($`>`$ 0.1) electron ratio were ignored.

The data toward SMC X-1 are occasionally contaminated by outbursts of the nearby ($``$27 away) transient source XTE J0111.2–7317 (Chakrabarty et al. 1998). The contamination was easily identified by the source’s $``$31-s pulsations, with powers reach up to $`>`$ 5000 in Leahy-normalized PDSs. We searched for weak contamination in PDSs as narrow statistically significant ($`>`$ 10 in Leahy-nomalized PDSs) peaks at $``$0.032 Hz (and/or at its harmonic frequencies). We excluded any data contaminated by XTE J0111.2–7317 from our analyses, as well as data containing dip-like features (e.g., Her X-1; Moon & Eikenberry 2001) or close to the eclipse (i.e., binary orbital phase between 0.9 and 0.1) of SMC X-1. We obtained a total of 150 data segments with an average length of $``$1360 s, giving a net observational duration of $``$204 ks. For spectral analysis, we considered only the spectrum between 2.5–25 keV range due to a PCA responsivity problem (R. Remillard 2001, private communication), and assumed a systematic uncertainty of 1 %.

In our analysis, we define “flares” to be the part of a light curve that has three or more consecutive data bins with photon counts larger than a threshold value of 3$`\sigma `$ Poisson noise above the mean photon count in a given light curve. We analyzed all the 150 light curves as follows searching for the flares. First, we binned each light curve to 4-s time resolution and made a flare list. We extended the search with 2- and 8-s time resolutions, excluding the flares already found with different time resolutions. We found a total of 323 flares, and fit them a Gaussian function to estimate width and peak intensity. We obtained $``$18 s and $``$15 s for the mean and rms standard deviation of FWHM; the reduced $`\chi ^2`$ (= $`\chi _\nu ^2`$) of the fits range from 0.5 to 1.9. The total duration above half-maximum is $``$5.8 ks, indicating that SMC X-1 spends $``$3 % of its time on flares.

## 3 Flare Examples and Correlations with Other Parameters

The rms variability of the 150 light curves ranges between 6 and 16 %, with a mean of $``$11 %. Figure 1 presents three light curves (and their PDSs) with very low (7 %), average (11 %), and very high (16 %) rms variability as examples representing three different variability levels. The mean photon count rates are 120.5 $`\pm `$ 8.5, 116.1 $`\pm `$ 12.9, and 146.4 $`\pm `$ 23.3, for Figure 1a, 1b, and 1c, respectively. No flare is found in Figure 1a, which has a low rms variability; three flares are found Figure 1b, which has an average rms variability. Strong flaring activity is well illustrated in Figure 1c, which has a high rms variability, with several apparent flares lasting for a few tens of seconds. At the flare peak, the photon count rate rises up to $``$2.5 times of that outside the flares. The inset shows a Gaussian fit to the intense flare at $`t`$ $``$ 900 s, with an FWHM estimated to be $``$23 s. Leahy-normalized PDSs of the three light curves (Figure 1d, 1e, and 1f) show strong QPO-like peaks around 10 mHz, as well as the peaks caused by the source’s coherent pulsations at $`\nu `$ $``$ 1.4 Hz (and their harmonics), indicating the existence of $``$100-s aperiodic variability independent of the flaring activity. The fractional rms (FRMS) amplitudes in the 2–50 mHz range are 7 $`\pm `$ 1, 9 $`\pm `$ 1, and 14 $`\pm `$ 1 % for 1d, 1e, and 1f, respectively.

Figure 2 presents the pulse profiles of the three light curves, all showing the double-peaked, smooth profile typical of SMC X-1 (e.g., Levine et al. 1993). They have the same pulsed fractions of $``$39 %, with their second peaks at phase 0.56 with respect to the first ones at phase 0. The ratio of the second peak to the first peak is $``$0.78 for Figure 2a, while it is $``$0.96 for 2b and 2c. No significant variation in the pulse peak ratio has been found within a given light curve, regardless whether it is obtained inside or outside flares.

Figure 3 compares the phase-averaged softness ratio, defined to be the ratio of the soft X-ray (2–5 keV) photon count rates to those of the hard X-ray (5–13 keV), of the three light curves in Figure 1 with the total photon count rates at 2–25 keV with 4-s resolution. The average softness ratio of the three light curves is invariant: 0.47 $`\pm `$ 0.03, 0.47 $`\pm `$ 0.02 , and 0.46 $`\pm `$ 0.04 for Figure 1a, 1b, and 1c, respectively. This contrasts with the softness ratio of the LMC X-4 flares (Moon, Eikenberry, & Wasserman 2002), which shows a strong linear correlation with the total photon count rate (inset in Figure 3).

We examined the spectral invariance of SMC X-1 over the flaring activity indicated by the constant softness ratio distribution more thoroughly via fitting the 32-s spectra obtained from the largest peak in each of the three light curves (i.e., the peak at $`t`$ $``$ 900 s in Figure 1a, at $`t`$ $``$ 1050 s in 1b, and at $`t`$ $``$ 1110 s in 1c) to a model spectrum. The model spectrum consists of a power-law component with a high-energy cutoff (for non-thermal magnetospheric emission) and a Gaussian component (for iron line emission), together with a component for photoelectric absorption by intervening interstellar matter. We fixed the central energy of the Gaussian component to be 6.7 keV based on previous results Angelini et al. (1991). Figure 4 compares the observed spectra with the best-fit model spectra, and Table 1 lists the best-fit parameters and $`\chi _\nu ^2`$ of the fits. Although all parameters are poorly constrained, the power-law index does not change significantly over the flaring activity – consistent with the softness ratio distribution (Figure 3). The difference in $`N_\mathrm{H}`$ may be caused by the motion of the precessing, tilted accretion disk suggested to be responsible for the super-orbital motion of SMC X-1, but the low energy limit of the spectral fits (i.e., 2.5 keV) makes it difficult to constrain $`N_\mathrm{H}`$ propertly, because the photoelectric absorption by intervening interstellar matter is expected to be most significant in the soft energy band.

In order to perform statistical analyses, we calculated the flare fraction, which we define to be the ratio of the integrated time that SMC X-1 is flaring (i.e., within FWHMs of the Gaussian fits) to the total observing time of a given parameter, and investigate its correlation with the parameter. Figure 5 shows the distribution of the flare fraction as functions of the binary orbital phase (5a), rms variability of the light curve (5b), and the pulse peak ratio (5c). Some important results are worth noticing: the flare fraction (1) is larger than 1.4 % over all the orbital phases (with its minimum at the orbital phase of $``$0.35), (2) shows a significant variation from phase to phase, and (3) increases with the rms variability of the light curve, as well as with the pulse peak ratio. The correlation between flare fraction and rms variability identified in Figure 5b remains very similar when we use the flare-subtracted rms variability.

## 4 Discussion

SMC X-1 shows active flares that occupy $``$3 % of the total observing time. The flares seem to occur over all binary orbital phases, and the flaring activity is proportional to the rms variability of the light curve (Figure 5). Except for the small change in the pulse peak ratio (Figure 2 and 5), no significant change is found along with the flaring activity. All these suggest that the SMC X-1 flares may be simple extensions of the persistent emission of a normal state with increased accretion rates, but without significant changes in the geometry of the accretion flows and the magnetosphere. Under the hypothesis that the double peaks in the pulse profile (Figure 2) are due to the two magnetic poles of SMC X-1, one simple explanation for the change in the pulse peak ratio may be that the increase in the accretion rate onto the fainter pole is higher than that onto the brighter pole during the flares. The invariant pulse peak ratio in a given light curve (regardless of flares), together with the correlation between the flare-subtracted rms variability and the flare fraction, on the other hand, may indicate the existence of variability whose time scale is longer than the flare recurrence time, $``$100 s.

The SMC X-1 flares differ from Type I X-ray bursts from LMXBs for various reasons, including the shape of profiles and spectral properties. While the profile of Type I bursts shows an abrupt increase with an exponential decay in most cases, the SMC X-1 flares have symmetric Gaussian shapes. The X-ray spectrum of the SMC X-1 flare is far from thermal spectrum of Type I bursts (although the thermal spectrum of an X-ray pulsar is not well constrained, so it is not completely excluded that the SMC X-1 flare spectrum is thermal). In addition, the SMC X-1 spectrum does not show any apparent variation within a flare, while Type I burst shows spectral cooling as the burst continues. The changes in the softness ratio and pulse profile of the SMC X-1 flares differ from those of the LMC X-4 flares (Figure 3; Moon et al. 2002), which is one of the most well known and regular flaring X-ray binaries.

The SMC X-1 flares recall the Type II X-ray bursts found in GRO J1744–28 mainly due to the spectral invarince over the flaring activity. In fact, the burst spectrum of GRO J1744–28 was found to be very similar to that of SMC X-1, with a similar photon index of $``$1.2 and e-folding energy of $``$14 keV (Lewin et al. 1996). One difference is that the SMC X-1 flares lack the post-flare dip which often follows the Type II bursts from GRO J1744–28. The rather gradual rise of the SMC X-1 flares, however, may be responsible for it via offering sufficient time to replenish the material in the accretion disk. This is consistent with the interpretation that the SMC X-1 flares are simple extensions of a normal state with increased accretion rates, which is probably the most compelling argument for attributing the origin of the SMC X-1 flares to an accretion disk instability. Given the difference between GRO J1744–28 and SMC X-1 (i.e., a transient low-mass system versus a persistent high-mass one), we consider that the magnetic field strength, suggested to be comparable for the two sources, is critical to understanding the Type II bursts from GRO J1744–28 and the X-ray flares from SMC X-1.

We note that SMC X-1 may be capable of experiencing a viscous instability, namely the Lightman-Eardley instability (Lightman & Eardley 1974), due to its $``$1029 G cm3 magnetic moment Li & van den Heuvel (1997). The radiation pressure in this case is comparable to the gas pressure around the inner disk radius (i.e., the magnetospheric radius), resulting in a viscous instability with slightly increased mass accretion (e.g., Cannizzo 1996, 1997). Because the instability develops near the inner edge of the accretion disk, it has an advantage in explaining bursts/flares with a short recurrence time. The viscosity parameter ($`\alpha `$) of the classical $`\alpha `$-disk (Shakura & Sunyaev 1973) at the transition radius between the “inner region” and the “middle region” around an 1.4 $`M_{}`$ neutron star is $`\alpha 216t_{\mathrm{vis}}^{3/2}\dot{M}_{17}`$, where $`t_{\mathrm{vis}}`$ is the viscous time scale in seconds and $`\dot{M}_{17}`$ is the mass accretion rate in units of 1017 g s-1. For a typical value $`\dot{M}_{17}`$ $``$ 20 for SMC X-1 (e.g., Wojdowski et al. 2000), the viscosity parameter is $`\alpha `$ $``$ 0.14 and 4.4 when $`t_{\mathrm{vis}}`$ is 1000 and 100 s, respectively. The value 4.4 for $`t_{\mathrm{vis}}`$ = 100 s is somewhat larger than generally expected, $`\alpha `$ $``$ 1. However, even larger values ($`\alpha `$ $``$ 10–100) were obtained for the dwarf nova HT Cas, possibly due to the patch nature of the accretion disk (Vrielmann, Hessman, & Horne 2002). Alternately, the strong magnetic field of SMC X-1 may help increase the viscosity parameter (R. Lovelace 2002, private communication). We need detailed numerical studies to investigate this scenario more thoroughly.

## 5 Summary and Conclusions

Through the analysis of the all publicly available RXTE data toward the X-ray pulsar SMC X-1, we find that the source is highly variable on short time scales ($`<`$ 1 h), and that Gaussian flares occur over all orbital phases with a $``$100 s recurrence time scale. The flares occupy $``$3 % of the total observing time, and the flaring activity is proportional to the overall variability of the source. While the pulse peak ratio shows a small systematic change along with the flaring activity, the PDSs, pulse profiles, softness ratios, and X-ray spectra during the flares are very similar to those outside the flares, indicating that the flares are probably extensions of a normal state just with increased accretion rates. This supports the interpretation that the SMC X-1 flares have their origin in an accretion disk instability and the suggestion that it may belong to a distinctive group of “bursting pulsars” with the “bursting pulsar” GRO J1744–28, owing to its $``$1011 G surface magnetic field. A viscous instability near the inner-edge of the accretion disk might be responsible for the SMC X-1 flares, although detailed studies on this scenario need to be done in the future.

We would like to thank the anonymous referee for the useful comments and suggestions. D.-S. M. acknowledges Akiko Shirakawa, Dong Lai, and Richard Lovelace for their comments. This research has made use of data obtained from the High Energy Astrophysics Science Archive Research Center provided by NASA’s Goddard Space Flight Center. D.-S. M. is supported by NSF grant AST-9986898. S. S. E. is supported in part by an NSF Faculty Early Careeer Development award (NSF-9983830). I. M. W. is supported by the NASA grant NAG-5-8356.


text:
# Fluctuations of self-flattening surfaces

## Abstract

We study the scaling properties of self-flattening surfaces under global suppression on surface fluctuations. Evolution of self-flattening surfaces is described by restricted solid-on-solid type monomer deposition-evaporation model with reduced deposition (evaporation) at the globally highest (lowest) site. We find numerically that equilibrium surface fluctuations are anomalous with roughness exponent $`\alpha 1/3`$ and dynamic exponent $`z_W3/2`$ in one dimension (1D) and $`\alpha =0(\mathrm{log})`$ and $`z_W5/2`$ in 2D. Stationary roughness can be understood analytically by relating our model to the static self-attracting random walk model and the dissociative dimer type deposition-evaporation model. In case of nonequilibrium growing/eroding surfaces, self-flattening dynamics turns out to be irrelevant and the normal Kardar-Parisi-Zhang universality is recovered in all dimensions.

Structural properties for fluctuating surfaces under thermal noise have been studied extensively. Equilibrium surfaces with proper surface tension are always rough in one dimension (1D) and display a roughening transition in two dimensions (2D). Higher dimensional surfaces are always smooth. Surface roughness is well documented and classified as the Edwards-Wilkinson (EW) universality class. The EW class is generic and robust for equilibrium surfaces with local surface tension. Only specific nonlinear contributions in nonequilibrium growth processes may become relevant and drive the system into other universality classes, e.g., the Kardar-Parisi-Zhang (KPZ) universality class .

In this paper, we introduce a new global mechanism to suppress surface fluctuations, besides ordinary local surface tension. We call it self-flattening mechanism to reduce growing (eroding) probability at the globally highest (lowest) point on the surface. This global type suppression makes the surface less rough, which may bring forth new universality classes for equilibrium and nonequilibirum surfaces. Inclusion of suppression at all local extremal points leads to less interesting layer-by-layer growth processes and the steady state surfaces are always smooth with finite fluctuation width.

We describe surface configurations in terms of integer height variables $`\{h(\stackrel{}{r})\}`$ at site $`\stackrel{}{r}`$ on a $`D`$-dimensional hypercubic lattice. They are subject to the restricted solid-on-solid (RSOS) constraint, $`h(\stackrel{}{r}+\widehat{e_i})h(\stackrel{}{r})=0,\pm 1`$ with $`\widehat{e_i}`$ a primitive lattice vector in the $`i`$-th direction ($`i=1,\mathrm{},D`$). The RSOS constraint effectively generates local surface tension which prevents indefinite growth of surface fluctuations for finite systems.

Evolution rule for the ordinary RSOS type monomer deposition-evaporation model is given as follows. First, select a site $`\stackrel{}{r}`$ randomly. Next, deposit a particle, $`h(\stackrel{}{r})h(\stackrel{}{r})+1`$, with probability $`p`$ or evaporate a particle, $`h(\stackrel{}{r})h(\stackrel{}{r})1`$, with probability $`q=1p`$. Any deposition/evaporation attempt is rejected if it would result in violating the RSOS constraint. Equilibrium surfaces at $`p=q`$ belongs to the EW class, while nonequilibrium growing/eroding surfaces at $`pq`$ the KPZ class.

For self-flattening surfaces, we need a slight variation of the evolution rule to incorporate the global suppression : only when deposition (evaporation) is attempted at the globally highest (lowest) site, the attempt is accepted with probability $`u`$ and rejected with probability $`1u`$. At $`u=1`$, the ordinary RSOS model is recovered. The $`u=0`$ case is trivial, because the surface is confined within initial surface width.

We perform numerical simulations, starting from a flat surface of linear size $`L`$ with periodic boundary conditions. We measure the surface fluctuation width $`W`$ as

$$W^2(L,t)=\frac{1}{L^D}\underset{\stackrel{}{r}}{}\left[h(\stackrel{}{r},t)\frac{1}{L^D}\underset{\stackrel{}{r}}{}h(\stackrel{}{r},t)\right]^2,$$
(1)
where $`\mathrm{}`$ represents the ensemble average with equal weights. Therefore, our simulations at $`p=q`$ correspond to the infinite temperature limit of equilibrium RSOS surfaces. The surface width satisfies the dynamic scaling relation

$$W(L,t)=L^\alpha f\left(t/L^{z_W}\right),$$
(2)
where the scaling function $`f(x)\mathrm{const}.`$ for $`x1`$ and $`f(x)x^\beta `$ $`(\beta =\alpha /z_W)`$ for $`x1`$.

First, we report the numerical results for equilibrium surfaces ($`p=q`$). For 1D, we run simulations for $`L=2^5,\mathrm{},2^{11}`$ at $`u=0.1`$, 0.3, 0.6, and 0.8, and average over at least 300 independent samples. In early time regime ($`tL^{z_W}`$), the surface width grows with time, $`Wt^\beta `$, and saturates to a finite value which increases with size, $`W_sL^\alpha `$.

In Fig. 1, we show the plot of $`\mathrm{ln}W`$ versus $`\mathrm{ln}t`$ at $`u=0.6`$ for various system sizes. The growth exponent $`\beta `$ is estimated by a simple straight line fitting of early time data for the largest system size $`L=2^{11}`$. Our estimate is $`\beta =0.22(1)2/9`$. We also analyze the data at other values of $`u`$ and find that $`\beta `$ does not vary with $`u`$.

In order to extract the stationary property, we average over data in the saturated regime ($`tL^{z_W}`$) for given $`L`$ to estimate $`W_s(L)`$. For efficient estimation of $`\alpha `$, we introduce effective exponents

$$\alpha _{eff}(L)=\mathrm{ln}\left[W_s(mL)/W_s(L)\right]/\mathrm{ln}m,$$
(3)
where $`m`$ is an arbitrary constant (here, we set $`m=2`$).

Effective exponents at various values of $`u`$ are plotted in Fig. 2. Close to $`u=1`$, our data show large corrections to scaling as expected, due to the presence of the EW fixed point ($`\alpha =1/2,\beta =1/4`$) at $`u=1`$. However, the asymptotic estimates seem to be independent of $`u`$. We estimate that $`\alpha =0.33(1)1/3`$ for all $`u`$. We check the dynamic scaling relation directly by plotting $`W/L^\alpha `$ versus $`t/L^{z_W}`$ in the inset of Fig. 1. Our data collapse very well with $`\alpha =1/3`$ and $`z_W=3/2`$ for all $`u`$, which are consistent with the above results.

This set of scaling exponents form a new universality class, distinct from the EW and any previously known growth-type universality class. It implies that the self-flattening dynamics is a relevant perturbation to the EW fixed point in 1D. Therefore, the continuum equation to describe self-flattening surfaces must contain a global-type nonlinear term. Further study in this direction is left for future research.

In case of 2D EW surfaces, it is well known that the surface width grows logarithmically with time and its saturated value also increases logarithmically with size. Especially, the saturated width $`W_s`$ scales for large $`L`$ as

$$W_s^2(L)\frac{1}{2\pi K_G}\mathrm{ln}L,$$
(4)
where $`K_G`$ is the effective coupling constant of the Gaussian model where equilibrium surface models flow into by renormalization group transformations. The ordinary RSOS model at the infinite temperature (our model at $`u=0`$) is known to take $`K_G=K_G^00.916`$.

Assume the dynamic scaling relation similar to Eq.(2) as

$$W^2(L,t)=\frac{1}{2\pi K_G}\mathrm{ln}\left[Lg\left(t/L^{z_W}\right)\right],$$
(5)
where the scaling function $`g(x)\mathrm{const}.`$ for $`x1`$ and $`g(x)x^{1/z_W}`$ for $`x1`$. Then, in early time regime ($`tL^{z_W}`$), the surface width grows as

$$W^2(t)\frac{1}{2\pi K_Gz_W}\mathrm{ln}t.$$
(6)
The amplitude ratio in Eqs.(4) and (6) yields a value of the dynamic exponent $`z_W`$. The EW surfaces take $`z_W=2`$ in all dimensions.

We run simulations on $`L\times L`$ lattices with $`L=2^3,\mathrm{},2^7`$ at $`u=0.1`$ and 0.5 and average over at least 300 independent samples. In Fig. 3(a), we plot $`W^2`$ against $`\mathrm{ln}t`$ at $`u=0.5`$. It shows a nice linear behavior in the early time regime. In Fig. 3(b), we plot $`W_s^2`$ against $`\mathrm{ln}L`$, which also shows a very nice linear behavior. We measure its slope and find that $`K_G0.92(1)`$ for all $`u`$, which is very close to $`K_G^0`$. In contrast to the 1D surfaces, the global suppression does not seem to change the asymptotic behavior of the stationary surface roughness. As can be seen in Fig. 3, it seems to shift $`W_s`$ only by a constant.

We measure the amplitude ratio by comparing two slopes in Figs. 3(a) and (b). We estimate $`z_W=2.5(1)5/2`$ for all $`u`$, which is clearly distinct from the EW value of 2. We also check the dynamic scaling relation of Eq.(5) by plotting $`W^2W_s^2`$ versus $`t/L^{z_W}`$ in Fig. 3(a). Our data collapse reasonably well with $`z_W=5/2`$ for all $`u`$. Together with our 1D results, we conclude that the self-flattening surfaces display a new type of scaling behavior and form a novel universality class.

The partition function for equilibrium self-flattening surfaces can be written as

$$Z=\underset{\mathrm{RSOS}\mathrm{conf}.}{}e^{\beta \left(h_{\mathrm{max}}h_{\mathrm{min}}\right)},$$
(7)
where the summation is over all height configurations satisfying the RSOS condition, $`\beta `$ a temperaturelike parameter, and $`h_{\mathrm{max}}`$ ($`h_{\mathrm{min}}`$) the globally maximum (minimum) height for a given configuration.

Global suppression for self-flattening dynamics is simply Metropolis type evolution algorithm with this partition function to reach the equilibrium. Deposition (erosion) at the globally highest (lowest) site increases the energylike term $`h_{\mathrm{max}}h_{\mathrm{min}}`$ by one unit and these attempts are accepted with Boltzmann type probability $`e^\beta `$. Any other deposition (erosion) attempts are always accepted, because they do not increase the energylike term. Of course, all attempts resulting in violation of the RSOS constraint are rejected. By identifying $`u=e^\beta `$, our model for self-flattening surfaces is exactly the same as Metropolis evolution with the above partition function.

Stationary property of this system can be understood analytically. In 1D, this system is equivalent to the so-called static self-attracting (timid) random walks . The surface can be mapped to the time trajectory of a random walker by identifying the height $`h(x)`$ at site $`x`$ with the walker position after $`x`$ steps. The system size $`L`$ becomes the total number of steps and the RSOS constraint limits one-step hopping distance to 0 or $`\pm 1`$.

In 1D, the energylike term is simply the number of distinct sites visited by the random walker up to $`L`$ steps. Random walk configurations with less visited sites are preferred. Such a random walker tends to visit previously visited sites, so the walk is self-attractive. Its typical displacements are known rigorously to scale as $`L^{1/(D+2)}`$  under the assumption that the visited sites form a compact cluster. In 1D, the cluster is obviously compact, so the roughness exponent in our model should be $`\alpha =1/3`$ in 1D.

In 2D, the self-flattening surfaces are completely different from the self-attracting walks. The former deals with membrane fluctuations, while the latter polymer fluctuations. In order to understand the scaling behavior of the self-flattening surfaces, we investigate the intricate relation between our model and the dissociative dimer deposition-evaporation model in equilibrium .

In the dimer model, we deposit or evaporate particles only in a dimer form aligned along the surface. There is a global evenness conservation law that the number of particles at each height level must be conserved modulo 2. This leads to a Boltzmann type factor in the partition function as

$$Z=\underset{\mathrm{RSOS}\mathrm{conf}.}{}\underset{h}{}\frac{1}{2}\left(1+z^{v_h}\right),$$
(8)
where the product is over all possible height levels and $`v_h`$ the number of particles at height level $`h`$. The dimer model corresponds to the $`z=1`$ case where only configurations obeying the evenness conservation law (all $`v_h`$ are even) survive in the partition function. At $`z=1`$, the model reduces to the ordinary monomer model.

The self-flattening surfaces correspond to the $`z=0`$ limit. Each term inside the product picks up a factor of $`\frac{1}{2}`$ if $`v_h0`$, otherwise a factor of unity. The number of height levels with nonzero $`v_h`$ (at least one particle) is $`h_{\mathrm{max}}h_{\mathrm{min}}`$. Therefore, the $`z=0`$ case is equivalent to the self-flattening surfaces at $`\beta =\mathrm{ln}2`$. In fact, the $`Q`$-mer generalization corresponds to the $`\beta =\mathrm{ln}Q`$ case .

¿From the Gaussian model type renormalization group argument, one can show that the 2D surface roughness is always logarithmic in the dimer model for $`1z<1`$ (see Eq.(4)) and its amplitude remains unchanged . Our numerical results for all $`u`$ are consistent with this. The dimer characteristics show up only in a form of corrections to scaling. Recently, it is suggested that the corrections to scaling should scale as $`\mathrm{ln}(\mathrm{ln}L)`$, which is confirmed for the dimer model at $`z=1`$. We find no evidence of this type of corrections to scaling in our model ($`z=0`$) and the leading corrections are constants. The origin of this discrepancy between the $`z=0`$ and $`z=1`$ case is not fully understood as yet.

Next, we consider nonequilibrium growing/eroding surfaces $`(pq)`$. We run simulations for $`L=2^5,\mathrm{},2^{11}`$ for 1D and $`L=2^3,\mathrm{},2^7`$ for 2D at $`p=1`$ with $`u=0.5`$ and $`u=1`$ (ordinary RSOS). In Fig.4, we plot $`\mathrm{ln}W_s`$ against $`\mathrm{ln}L`$ and, in the inset, $`\mathrm{ln}W`$ against $`\mathrm{ln}t`$ for the largest system size in 1D and 2D, respectively. We do not find any noticeable change of $`W`$ ascribed to the global suppression. We estimate that $`\alpha 0.50(1)`$ and $`\beta 0.32(1)`$ for 1D and $`\alpha 0.40(1)`$ and $`\beta 0.24(1)`$ for 2D, which are consistent with the results for the ordinary RSOS model . We conclude that the global suppression is irrelevant to nonequilibrium growing/eroding surfaces.

In summary, we studied the scaling properties of the self-flattening surfaces in 1D and 2D. Equilibrium surfaces display dynamic scaling behavior distinct from the EW class and form a new universality class. We show that stationary roughness can be understood through mapping our model to self-attracting random walks in 1D and dissociative dimer type deposition-evaporation model in 2D. In higher dimensions, the surfaces are always smooth. In contrast, nonequilibrium self-flattening surfaces belong to the ordinary KPZ universality class. This implies that the self-flattening dynamics is strong enough to dominate over the EW type local surface tension term, but weaker than the KPZ type nonlinear term. It would be very interesting to find a continuum-type equation to govern the self-flattening dynamics.

We thank Marcel den Nijs, Deok-Sun Lee, and Lucian Anton for useful discussions. This research is supported in part by by grant No. R01-2001-00025 (YK) and by grant No. 2000-2-11200-002-3 \[HP\] from the Basic Research Program of KOSEF.