text:
# Relaxation and Diffusion in a Globally Coupled Hamiltonian System

## I Introduction

Relaxation to thermal equilibrium has been studied in Hamiltonian systems with long-range interactions ruffo-94 ; antoni-95 ; tsuchiya-96 ; yamaguchi-96 ; torcini-99 ; latora-02 . One of the characteristic phenomena in the relaxation process is anomalous diffusion, since normal diffusion is expected at equilibrium. Anomalous diffusion was firstly investigated in a one-dimensional chaotic map to describe enhanced diffusion in Josephson junctions geisel-85 , and is observed in many systems both numerically kaneko-89 ; klafter-94 ; latora-98 ; latora-99b ; torcini-99  and experimentally solomon-93 .

Anomalous diffusion is also observed in Hamiltonian dynamical systems. It is explained as due to power-type distribution functions klafter-94 ; meiss-85 ; zaslavsky-93  of trapping and untrapping times of the orbit in the self-similar hierarchy of cylindrical cantori geisel-95 . Self-similarity is expected to be one of the important concepts to understand statistics and motion in Hamiltonian systems, but cannot be the main feature in systems with many degrees of freedom. Then, as the first step of approaching the study of self-similarity, we have to clarify when anomalous diffusion appears, and what is the origin of the anomaly.

Latora et al. latora-99b  discussed the relation between the process of relaxation to equilibrium and anomalous diffusion in a globally coupled rotator system, by comparing the time series of the temperature and of the mean squared displacement of the phases of the rotators. They showed that anomalous diffusion changes to a normal diffusion after a crossover time, and that the crossover time coincides with the time when the canonical temperature is reached. They also claim that anomalous diffusion occurs in the quasi-stationary states, which appear before the system goes towards equilibrium.

The crossover from anomalous to normal diffusion determines the time when the anomalous diffusion finishes. However, it is not clearly pointed out when the anomalous diffusion starts, and hence the study of the relation between the relaxation process and anomalous diffusion is still not complete. Moreover, in Ref.latora-99b , the numerical calculations were performed by using only one type of initial condition, but different types of initial condition may change the conclusion koyama-01 .

In this article, we study the globally coupled rotator system considered in Ref.latora-99b , and we exhibit the relation between relaxation to equilibrium and anomalous diffusion with a different type of initial condition from one used in Ref.latora-99b . Then we show that diffusion becomes anomalous if and only if the state is neither stationary nor quasi-stationary. In other words, diffusion is shown to be normal in quasi-stationary states, although a stretched exponential correlation function is present, contrary to previous claims that report power law type function montemurro-02 . Simple scaling laws of the correlation function imply that the result holds irrespective of degrees of freedom.

This paper is organized as follows. The model, initial condition and observed quantities are described in Sec.II. In Sec.III, we study relaxation process, which we divide into three stages, quasi-stationary, relaxational and equilibrium stages. Diffusion process in each stage is investigated in Sec.IV by using stretched exponential correlation functions of momenta. Dependence on degrees of freedom is also reported both in Secs.III and IV. The last section V is devoted to summary.

## II Model, Initial Condition and Observed Quantities

The model considered in this paper has $`N`$ classical and identical rotators confined to move on the unit circle, and the Hamiltonian is composed of a kinetic and a potential part ruffo-94 ; antoni-95 ; yamaguchi-96 ; latora-99b ; montemurro-02  ,

$$H=K+V=\underset{j=1}{\overset{N}{}}\frac{p_j^2}{2}+\frac{1}{2N}\underset{i,j=1}{\overset{N}{}}[1\mathrm{cos}(\theta _i\theta _j)].$$
(1)
The $`N`$ particles are globally coupled through the mean field defined as

$$𝑴=\frac{1}{N}\underset{j=1}{\overset{N}{}}(\mathrm{cos}\theta _j,\mathrm{sin}\theta _j)=M(\mathrm{cos}\varphi ,\mathrm{sin}\varphi ),$$
(2)
where the modulus $`M(0M1)`$ represents the magnetization of this system. We remark that the potential $`V`$ and the kinetic energy $`K`$ are related to the magnetization $`M`$ as follows:

$$2V/N=1M^2,2K/N=2U1+M^2,$$
(3)
where $`U`$ is the energy per particle, i.e. $`U=E/N`$, and $`E`$ is the total energy. The free energy of this system has been obtained in the canonical ensemble ruffo-94 ; antoni-95 ; latora-99a , and it has been shown that system (1) has a second-order phase transition at the critical energy $`U_c=0.75`$. If the energy $`U`$ is greater than the critical energy, the largest Lyapunov exponent goes to zero in the thermodynamic limit ($`N\mathrm{}`$) firpo-98 . Then, all rotators freely rotate, and diffusion becomes ballistic. On the contrary, if $`U`$ is small compared to $`U_c`$, all rotators are trapped in the potential well and no diffusion occurs. We are therefore interested in a value of the energy which is near but less than the critical energy in order to allow some particle diffusion. Hereafter, we set $`U=0.69`$ (a value studied also in Refs.latora-99a ; latora-99b ).

The canonical equations of motion for system (1) can be cast in a form that uses the mean field (2) as follows:

$$\frac{\text{d}\theta _j}{\text{d}t}=p_j,\frac{\text{d}p_j}{\text{d}t}=M(t)\mathrm{sin}(\theta _j\varphi (t)),(j=1,\mathrm{},N).$$
(4)
We numerically integrate Eq.(4) by using 4th order symplectic integrators yoshida-93 ; mclachlan-92 . The time slice of the integrator is set at $`\mathrm{\Delta }t=0.2`$ or $`0.4`$, and it suppresses the relative energy error down to $`|\mathrm{\Delta }E/E|<5\times 10^7`$.

We have performed the integrations starting from $`M(0)=0`$. To prepare these initial conditions numerically, we set $`q_j(0)=2\pi j/N`$, and $`p_j(0)`$ is taken from a uniformly random distribution whose support is $`[\overline{p},\overline{p}]`$, where the value $`\overline{p}`$ is chosen to get the energy density $`U`$. The total momentum $`_{j=1}^Np_j`$ is an integral of the motion and we initially set it to zero. This initial state corresponds to a local entropy minimum antoni-02 , and to a stationary stable solution to the Vlasov-Poisson equation antoni-95 , although the system goes towards Gibbs equilibrium due to finite size effects latora-02 . With respect to the $`M(0)=1`$ initial condition chosen in Ref.latora-99b , the one we choose has the advantage of being a quasi-stationary state from the start.

We numerically observe the time series of two quantities. One is for the relaxation process, and the other is for the diffusion process.

To observe the relaxation process, we use the magnetization $`M(t)`$. Note that observing $`M(t)`$ corresponds to observing $`2K(t)/N`$ by using Eq.(3), and $`2K(t)/N`$ is the time series of the temperature, since the canonical average of $`2K/N`$ coincides with the canonical temperature.

To observe the diffusion process, we introduce the mean square displacement of phases $`\sigma _\theta ^2(t)`$ defined as

$$\sigma _\theta ^2(t)=\frac{1}{N}\underset{j=1}{\overset{N}{}}[\theta _j(t)\theta _j(0)]^2=[\theta _j(t)\theta _j(0)]^2_N.$$
(5)
The symbol $`_N`$ represents the average over all the $`N`$ rotators. The quantity $`\sigma _\theta ^2(t)`$ typically scales as $`\sigma _\theta ^2(t)t^\alpha `$, and the diffusion is anomalous when $`\alpha 1,2`$, while it is normal when $`\alpha =1`$ and ballistic for $`\alpha =2`$. The quantity $`\sigma _\theta ^2(t)`$ can be rewritten by using the correlation function of momenta $`C_p(t;\tau )`$ as

$$\begin{array}{cc}\hfill \sigma _\theta ^2(t)& =_0^t\text{d}t_1_0^t\text{d}t_2p_j(t_1)p_j(t_2)_N\hfill \\ & =2_0^t\text{d}s_0^{ts}\text{d}\tau C_p(s;\tau ),\hfill \end{array}$$
(6)
where $`C_p(t;\tau )`$ is defined as

$$C_p(t;\tau )=p_j(t+\tau )p_j(\tau )_N.$$
(7)
Moreover, if the system is stationary and $`C_p(t;\tau )`$ does not depend on $`\tau `$ accordingly,

$$C_p(t;\tau )=C_p(t;0),(\tau >0)$$
(8)
then Eq.(6) is simplified as

$$\sigma _\theta ^2(t)=2_0^t(ts)C_p(s;0)\text{d}s.$$
(9)
## III Relaxation Process

Temporal evolutions of $`M(t)`$ are shown in Fig.1. In order to suppress fluctuations, we have calculated averages over realizations. Throughout this article, unless no comments appear, the numbers of realizations are $`n=1000,100`$ and $`8`$ for $`N=100,1000`$ and $`10000`$, respectively. We divide the temporal evolutions into three stages, Stage-I,II and III. In Stage-I, the value of magnetization is almost constant but smaller than the canonical value. After Stage-I, magnetization rapidly increases towards its equilibrium value $`M_{\text{eq}}`$, and we call this time interval Stage-II. Finally the system reaches equilibrium, during Stage-III.

Let us define boundary times between Stage-I and II, $`t_{\text{I/II}}`$, and between Stage-II and III, $`t_{\text{II/III}}`$, as follows. The magnetization takes the local minimum at $`t_{\text{min}}`$, and we adopt $`t_{\text{I/II}}=t_{\text{min}}`$. We define the other boundary time $`t_{\text{II/III}}`$ as the first passage time which satisfies $`M(t)=0.99M_{\text{eq}}`$. Values of the two boundary times are reported in Fig.2 as functions of degrees of freedom. The local minimum time is proportional to $`N^{1.7}`$ for $`N100`$ with our initial condition $`M(0)=0`$, as with another initial condition $`M(0)=1`$ zanette-02 . For small $`N`$, we cannot neglect the initial time region $`t<6`$ in which the level of $`M(t)`$ goes to $`O(1/\sqrt{N})`$ coming from the law of large numbers (see Fig.3(d)), and hence the power law breaks. The power law recovers by subtracting the initial increasing time $`6`$ from $`t_{\text{I/II}}`$ as shown in Fig.2, i.e. $`t_{\text{I/II}}6N^{1.7}(N10)`$.

A theoretical prediction of $`t_{\text{II/III}}`$, the upper curve in Fig.2, is obtained by fitting the magnetization $`M(t)`$ as hyperbolic tangent function,

$$M(t)=[1+\mathrm{tanh}(a(\mathrm{log}_{10}tb))]c+d.$$
(10)
The parameter $`d`$ represents the initial level of $`M(t)`$, and $`c`$ the half width between initial and equilibrium levels of $`M(t)`$. The product $`ac`$ is the slope at $`\mathrm{log}_{10}t=b`$, i.e. $`ac=\text{d}M/\text{d}(\mathrm{log}_{10}t)|_{\mathrm{log}_{10}t=b}`$, and $`10^b`$ is the time scale. As shown in Fig.3, these four parameters are fitted as

$$\begin{array}{cc}& a(N)=\frac{\sqrt{N}}{100c(N)},10^{b(N)}=\frac{1}{9}N^{1.7},\hfill \\ & c(N)=\frac{(M_{\text{eq}}d(N))}{2},d(N)=\frac{1.7}{\sqrt{N}}.\hfill \end{array}$$
(11)
By using the scaling law (11), we can predict when $`M(t)`$ reaches a given threshold level, $`M_{\text{th}}`$, as a function of $`N`$. Let $`t_{\text{th}}`$ be the threshold time, which satisfies $`M(t_{\text{th}})=M_{\text{th}}`$, then $`t_{\text{th}}`$ is expressed as

$$t_{\text{th}}=10^b\left(\frac{M_{\text{th}}d}{M_{\text{eq}}M_{\text{th}}}\right)^{{\displaystyle \frac{\mathrm{ln}10}{2a}}}.$$
(12)
In Fig.2, $`t_{\text{th}}`$ is reported for $`M_{\text{th}}=0.99M_{\text{eq}}`$, and the prediction is in good agreement with numerical results. We remark that, roughly speaking, $`t_{\text{I/II}}`$ is asymptotically proportional to $`N^{1.7}`$.

The system seems quasi-stationary in Stage-I. The existence of quasi-stationary states for sufficiently long time has been questioned in Ref.zanette-02 . We will answer to the question by observing dependence on $`\tau `$ of the correlation function $`C_p(t;\tau )`$ in Sec.IV.2.

## IV Diffusion Process

As described in Eq.(6), the mean square displacement $`\sigma _\theta ^2(t)`$ is obtained from correlation function of momenta $`C_p(t;\tau )`$, and hence we study diffusion process by observing the correlation function. We start from the simplest stage, Stage-III, because we may use the simple expression (9). Next, we progress to Stage-I, where we expect that the system is quasi-stationary, and that we may use Eq.(9) again. In non-stationary stage, Stage-II, we check whether diffusion is of a power type. Finally we investigate dependence on degrees of freedom for some important parameters.

### IV.1 Diffusion at Equilibrium

Assuming the system has reached equilibrium at $`t=t_{\text{eq}}`$, we observe $`C_p(t;t_{\text{eq}})`$ and $`\sigma _\theta ^2(t;t_{\text{eq}})`$, where

$$\sigma _\theta ^2(t;t_{\text{eq}})=[\theta _j(t+t_{\text{eq}})\theta _j(t_{\text{eq}})]^2_N.$$
(13)
At equilibrium we may assume that the system is stationary,

$$C_p(t;t_{\text{eq}}+\tau )=C_p(t;t_{\text{eq}}),(\tau >0)$$
(14)
and hence

$$\sigma _\theta ^2(t;t_{\text{eq}})=2_0^t(ts)C_p(s;t_{\text{eq}})\text{d}s.$$
(15)
Now let us consider the correlation function for $`N=1000`$. We adopt $`t_{\text{eq}}=2^{20}10^6`$ which is long enough to reach equilibrium imaging from Fig.1. The correlation function $`C_p(t;t_{\text{eq}})`$ is reported in Fig.4(a), and is well approximated by the stretched exponential function phillips-96

$$C_p(t;t_{\text{eq}})=0.47\mathrm{exp}[(t/410)^{0.32}],$$
(16)
rather than by a pure exponential (see the inset of Fig.4(a) which is a log-linear plot of $`C_p(t;t_{\text{eq}})`$).

We remark that a stretched exponential function $`\mathrm{exp}[x^\beta ]`$ with a small exponent $`|\beta |<<1`$ is indistinguishable from a power-type function in the region $`|\beta \mathrm{ln}x|<<1`$:

$$\mathrm{exp}[x^\beta ]=\mathrm{exp}[\mathrm{exp}(\beta \mathrm{ln}x)]\mathrm{exp}[1\beta \mathrm{ln}x]=x^\beta /\text{e}.$$

However the fitting function (16) well agree with numerical result even around $`|0.32\mathrm{ln}(t/410)|=1`$, whose two solutions are $`t18,9330`$. We therefore adopt a stretched exponential function as an approximation of $`C_p(t;t_{\text{eq}})`$.

By using the fitting function (16) and Eq.(15), we numerically reproduce $`\sigma _\theta ^2(t;t_{\text{eq}})`$, and the reproduced curve well approximates the numerical result as shown in Fig.4(b). Note that $`\sigma _\theta ^2(t;t_{\text{eq}})`$ is proportional to $`t^2`$ in the limit of $`t0`$, since $`C_p(s;t_{\text{eq}})`$ in Eq.(15) goes to the constant $`C_p(0;t_{\text{eq}})`$. On the other hand, in the limit of $`t\mathrm{}`$, $`\sigma _\theta ^2(t;t_{\text{eq}})`$ is proportional to $`t`$, because both $`C_p(s;t_{\text{eq}})`$ and $`sC_p(s;t_{\text{eq}})`$ are almost zeros in long time region, and hence their integrals become constants. The crossover from $`t^2`$ to $`t`$ is also observed if we assume an exponential correlation function, and hence we conclude that diffusion at equilibrium is normal as expected although a stretched exponential is present.

### IV.2 Diffusion in Quasi-Stationary State

Except for Stage-III, we cannot expect stationarity to hold (8) any more. However, from the temporal evolutions of $`M(t)`$, Fig.1, we may expect quasi-stationarity in Stage-I,

$$C_p(t;\tau )=C_p(t;0)+ϵ(t;\tau ),$$
(17)
where $`\tau `$ belongs to Stage-I and $`ϵ(t)`$ is suitably small.

The correlation function $`C_p(t;\tau )`$ for various values of $`\tau `$ is reported in Figs.5(a) and (b) for $`N=1000`$, and the relative error of correlation function, defined as

$$R(\tau )=\underset{t}{\mathrm{max}}\frac{|ϵ(t;\tau )|}{C_p(0;0)},$$
(18)
is also reported in Fig.5(c) as a function of $`\tau `$. The error $`R(\tau )`$ stays small up to the end of Stage-I, and hence we conclude that the system is quasi-stationary in Stage-I. We believe that the quasi-stationary states correspond to stationary stable states of the Vlasov equation yamaguchi-03 . We remark that $`R(\tau )`$ is constant in Stage-III again due to stationarity at equilibrium.

It seems natural that we regard $`C_p(t;\tau )`$ as a series of stretched exponential functions of $`t`$ rather than power-type functions, since this function fits $`C_p(t;\tau )`$ in more than two decades of time (power law fits of the correlation functions hold in one decade). Moreover, at equilibrium, $`C_p(t;t_{\text{eq}})`$ is also a stretched exponential rather than a pure exponential, as shown in Fig.4.

In the quasi-stationary region, Stage-I, the mean square displacement $`\sigma _\theta ^2(t)`$ can be derived by the correlation function $`C_p(t;0)`$, which is reported in Fig.6 for $`N=100,1000`$ and $`10000`$. We approximate $`C_p(t;0)`$ by a stretched exponential function as

$$\begin{array}{cc}\hfill N=100:& C_p(t;0)=0.38\mathrm{exp}[(t/20)^{0.68}],\hfill \\ \hfill N=1000:& C_p(t;0)=0.38\mathrm{exp}[(t/180)^{0.91}],\hfill \\ \hfill N=10000:& C_p(t;0)=0.38\mathrm{exp}[(t/2200)^{0.90}].\hfill \end{array}$$
(19)
The prefactor $`0.38`$ comes from $`C_p(0;0)=2K(0)/N`$.

Using the approximate functions (19) and Eq.(9), we are able to reproduce $`\sigma _\theta ^2(t)`$, as shown in Fig.7. The approximation is good in Stage-I, i.e. in the quasi-stationary time region, irrespective of the value of $`N`$. Consequently, there is no anomaly in diffusion in Stage-I, since the diffusion is explained by stretched exponential correlation function.

### IV.3 Diffusion in Non-Stationary State

After the quasi-stationary region, diffusion becomes anomalous, which is faster than normal diffusion, in Stage-II. If we fit $`\sigma _\theta ^2(t)`$ by a power-type function $`t^\alpha `$ in Stage-II, the exponent $`\alpha `$ is estimated as $`1.54,1.59`$ and $`1.74`$ for $`N=100,1000`$ and $`10000`$ respectively. The values of exponent tend to increase as $`N`$ increases as reported for the system having the so-called $`2`$-dimensional egg-crate potential torcini-99 . On the other hand, the duration in which diffusion is anomalous becomes shorter and shorter in logarithmic time scale as $`N`$ increases , in accordance with the sharper change of $`M(t)`$. Moreover, $`\sigma _\theta ^2/t^\alpha `$ is not constant, but has a wave in Stage-II (see Fig.8). Hence we guess that the anomaly in diffusion is not anomalous diffusion taking a power-type function, but a transient anomaly due to non-stationarity of Stage-II.

Let us proceed to investigate the origin of anomaly in diffusion. We focus on the behavior for $`N=1000`$. The mean square displacement $`\sigma _\theta ^2(t)`$ is perfectly determined by the correlation function $`C_p(t;\tau )`$ using Eq.(6), once we assume that $`C_p(t;\tau )`$ is a series of stretched exponential functions. We introduce three parameters, $`C_p(0;\tau )`$, $`t_{\text{corr}}(\tau )`$ and $`\beta (\tau )`$, to describe the stretched exponential function as

$$C_p(t;\tau )=C_p(0;\tau )\mathrm{exp}[\{t/t_{\text{corr}}(\tau )\}^{\beta (\tau )}].$$
(20)
We investigate which of the three parameters is the most important to yield anomaly in diffusion.

The strategy is as follows. We reproduce $`\text{d}\sigma _\theta ^2(t)/\text{d}t`$ by using the three parameters and the formula

$$\frac{\text{d}\sigma _\theta ^2}{\text{d}t}(t)=2_0^t\text{d}\tau C_p(0;\tau )\mathrm{exp}[\{(t\tau )/t_{\text{corr}}(\tau )\}^{\beta (\tau )}].$$
(21)
We consider the first derivative of $`\sigma _\theta ^2`$ instead of $`\sigma _\theta ^2`$ itself, because the former requires only single integration while the latter requires double integrations (6). We first omit the dependence on $`\tau `$ of the parameter $`C_p(0;\tau )`$ and fix it to a constant value to observe how it affects the anomaly in diffusion. We then fix the two other parameters $`t_{\text{corr}}(\tau )`$ and $`\beta (\tau )`$ to determine their effect on the mean square displacement.

From the numerical results of $`C_p(t;\tau )`$, Fig.5(b), we determine the values of three parameters $`C_p(0;\tau ),t_{\text{corr}}(\tau )`$ and $`\beta (\tau )`$ at some value of $`\tau `$ by using the least square method. The discrete values of the parameters are not enough to reproduce $`d\sigma _\theta ^2(t)/dt`$ accurately, and then we approximate the parameters by hyperbolic tangent functions as follows:

$$\begin{array}{cc}\hfill C_p(0;\tau )& =0.046[1+\mathrm{tanh}(2.5(\mathrm{log}_{10}\tau 4.35))]+0.385,\hfill \\ \hfill t_{\text{corr}}(\tau )& =80[1+\mathrm{tanh}(1.5(\mathrm{log}_{10}\tau 3.4))]+170,\hfill \\ \hfill \beta (\tau )& =0.31[1+\mathrm{tanh}(1.5(\mathrm{log}_{10}\tau 3.8))]+0.29.\hfill \end{array}$$
(22)
The hyperbolic tangent functions are in good agreement with numerical results, as shown in Fig.9. To confirm the validity of the approximation, we reproduced $`\text{d}\sigma _\theta ^2/\text{d}t`$ using Eqs.(21) and (22), and the reproduced one is in good agreement with numerical results, as shown in Fig.10(a).

If we fix $`C_p(0;\tau )`$ at its middle value, $`0.431`$, we find that the dependence on $`\tau `$ of $`C_p(0;\tau )`$ does not affect significantly $`\text{d}\sigma _\theta ^2/\text{d}t`$, as shown in Fig.10(b). By fixing $`t_{\text{corr}}(\tau )`$ at its middle value, $`250`$, we obtain the same conclusion for $`t_{\text{corr}}(\tau )`$ as for $`C_p(0;\tau )`$, particularly in Stage-II (see Fig.10(c)). On the contrary, if we fix $`\beta (\tau )`$ at $`0.6`$ or $`0.9`$, we observe no anomaly in diffusion as shown in Fig.10(d), because $`\text{d}\sigma _\theta ^2(t)/\text{d}t`$ is proportional to $`t`$ and is constant in short and long time regions respectively, and the same behavior is obtained at equilibrium (see Fig.4(b)). Consequently, among the three parameters, $`\beta (\tau )`$ plays a crucial role to produce anomaly in diffusion.

### IV.4 Dependence on Degrees of Freedom

In Stage-I (and III), we fit $`C_p(t;0)`$ (resp. $`C_p(t;t_{\text{eq}})`$) by a stretched exponential function, which has three parameters: $`C_p(0;0)`$, $`t_{\text{corr}}(0)`$ and $`\beta (0)`$ (resp. $`C_p(0;t_{\text{eq}})`$, $`t_{\text{corr}}(t_{\text{eq}})`$ and $`\beta (t_{\text{eq}})`$). In order to obtain scaling laws for the parameters, we show them as functions of degrees of freedom $`N`$. The parameters $`C_p(0;0)`$ and $`C_p(0;t_{\text{eq}})`$ represent temperature at $`t=0`$ and at equilibrium respectively, and hence they do not depend on $`N`$. We therefore focus on the other four parameters, $`t_{\text{corr}}(0)`$, $`\beta (0)`$, $`t_{\text{corr}}(t_{\text{eq}})`$ and $`\beta (t_{\text{eq}})`$. The correlation functions, $`C_p(t;0)`$ and $`C_p(t;t_{\text{eq}})`$, are shown in Fig.11, and values of the four parameters are reported as functions of $`N`$ in Fig.12.

For large $`N`$, $`N200`$, the correlation times are proportional to $`N`$, that is $`t_{\text{corr}}(0)=N/5`$ and $`t_{\text{corr}}(t_{\text{eq}})=N/2`$, and the stretching exponents $`\beta (0)`$ and $`\beta (t_{\text{eq}})`$ are almost constants. We expect that these scaling lows for the four quantities are kept even in the thermodynamic limit, although they break for small $`N`$ where $`t_{\text{corr}}(0)`$ is larger than $`N/5`$ and $`\beta (0)`$ is smaller than the constant. The duration of Stage-I, $`t_{\text{I/II}}`$, is around $`23`$ for $`N=100`$, and hence $`t_{\text{corr}}(0)`$ and $`\beta (0)`$ are estimated mainly not in Stage-I, but in Stage-II from Fig.11(a). In Stage-II, $`t_{\text{corr}}(\tau )`$ and $`\beta (\tau )`$ are increasing and decreasing functions of $`\tau `$ respectively (see Fig.9), and hence $`t_{\text{corr}}(0)`$ and $`\beta (0)`$ are larger and smaller than expected values, respectively.

## V Summary

As a summary, we have investigated the relation between relaxation and diffusion in a Hamiltonian system with long-range interactions. The relaxation process is divided into three stages: quasi-stationary, relaxational and equilibrium. We showed that diffusion becomes anomalous only in the second non-stationary stage, where magnetization is increasing and goes towards to the canonical value. The result mentioned above does not depend on the number of degrees of freedom, at least from $`N=100`$ to $`10000`$.

The interval where the anomaly in diffusion appears becomes shorter and shorter in logarithmic time scale as $`N`$ increase corresponding to a sharper change of magnetization. Moreover, a detailed investigation exhibits the absence of power-type diffusion even in the non-stationary stage. We guess that anomaly in diffusion is a transient anomaly due to non-stationarity.

Diffusion is obtained by integrating the correlation function of momenta $`C_p(t;\tau )`$ and the correlation function is approximated by a series of stretched exponential functions $`C_p(t;\tau )=C_p(0;\tau )\mathrm{exp}[(t/t_{\text{corr}}(\tau ))^{\beta (\tau )}]`$. Among the three parameters, $`C_p(0;\tau )`$, $`t_{\text{corr}}(\tau )`$ and $`\beta (\tau )`$, the stretching exponent $`\beta (\tau )`$ plays a crucial role to yield anomaly in diffusion. If we assume that $`\beta (\tau )`$ is a constant, we never observe anomaly in diffusion. This result is consistent with the fact that anomaly in diffusion does not appear in (quasi-)stationary state, because correlation function $`C_p(t;\tau )`$, and $`\beta (\tau )`$ accordingly, are almost invariant with respect to $`\tau `$.

We also investigated scaling laws concerning degrees of freedom $`N`$. The duration of quasi-stationary stage is proportional to $`N^{1.7}`$, and relaxation time, at which the system reaches at equilibrium, is also proportional to $`N^{1.7}`$ asymptotically although some corrections must be added. In both quasi-stationary and equilibrium stages, $`t_{\text{corr}}`$ is proportional to $`N`$, and $`\beta `$ is almost constant. These simple scaling laws imply that fitting by stretched exponential functions is valid irrespective of degrees of freedom.

We have not understood the theoretical reason of the appearance of a stretched exponential function. If we assume that several time scales with exponential correlation function, $`\mathrm{exp}(t/t_{\text{corr}})`$, are present, and we assume probability distribution function of $`t_{\text{corr}}`$, $`P(t_{\text{corr}})`$, then we obtain a stretched exponential function $`P(t_{\text{corr}})\mathrm{exp}(t/t_{\text{corr}})\text{d}t_{\text{corr}}`$ by choosing suitable forms for $`P(t_{\text{corr}})`$ palmer-84 ; pelcovits-83 . In our model, $`P(t_{\text{corr}})`$ corresponds to the distribution of time scales of individual rotators. The investigation of the macrovariable $`C(t;\tau )`$ in relation with the microvariables of the individual particle correlation functions will be a subject of future work.

###### Acknowledgements.

I thank Stefano Ruffo for a careful reading of the manuscript and useful comments. I acknowledge valuable discussions with Alessandro Torcini, Freddy Bouchet and Julien Barré.


format: text

text:
# Orbital periods of twenty two subdwarf B stars

## 1 Introduction

Subdwarf B (sdB) stars can be identified with models for Extreme Horizontal Branch (EHB) stars. The surface gravities and temperatures of sdB stars suggest that they have helium cores of mass $`0.5\text{M}_{}`$ and thin hydrogen envelopes of mass $`0.02\text{M}_{}`$ \[Heber et al. 1984, Saffer et al. 1994\]. A recent asteroseismological study of an sdB star results in a value for its mass of 0.49$`\pm `$0.02 $`\text{M}_{}`$ \[Brassard et al. 2001\]. Several evolutionary scenarios have been proposed to explain the formation of sdB stars, in particular the loss of the hydrogen envelope. Evolution within a binary star is an effective method for envelope removal, and yet it is hard to see why this should have happened to a horizontal branch star since it would have been much larger during its preceding red giant stage. A solution to this problem was presented by D’Cruz et al. \[D’Cruz et al. 1996\] who found that if a red giant star with a degenerate helium core loses its hydrogen envelope when it is within $``$0.4 magnitudes of the tip of the red giant branch, the core can go on to ignite helium, despite the dramatic mass loss, and may then appear as an sdB star. The advantage of this model is that it very nicely explains the masses of sdB stars as a consequence of the core mass at the helium flash. D’Cruz et al. \[D’Cruz et al. 1996\] supposed that mass loss occurred because of an enhancement of the stellar wind, but it could as well have been driven by binary interaction.

If sdB stars do form within binary systems and if they still have their companions, then the companions must be low-mass main-sequence stars or compact stellar remnants to avoid outshining the sdB star. If so, it is probable that in many cases the companions were unable to cope with the mass transferred from the sdB progenitor and a single “common” envelope formed around the two stars. Driving off such envelopes drains energy and angular momentum from the binary orbit, which as a result becomes much smaller than it was at the start of mass transfer \[Webbink 1984\]. It is therefore possible that many sdB stars are now members of close binary systems. Maxted et al. \[Maxted et al. 2001\] found exactly this, discovering 21 binary sdB stars in a sample of 36, suggesting, after allowance for detection efficiency, that some two-thirds of all sdB stars are in short period binary systems ($`P10\mathrm{d}`$). The other third seems to be made up a combination of long period binary stars that avoided a common envelope phase \[Green, Liebert & Saffer 2000\] and apparently single sdB stars.

If D’Cruz et al.’s \[D’Cruz et al. 1996\] model of the formation of sdB stars is correct, then the stage immediately prior to mass transfer in binary sdBs is well defined. This, together with the fact that the binary does not have enough time to change its orbital period significantly following its emergence from the common envelope, makes the sdB stars a superb population for testing models of the common envelope phase. Moreover, the detection of sdB binary stars is not compromised by the strong and poorly understood selection effects that plague other populations of close binary stars, such as the cataclysmic variable stars. The properties of sdB binaries (e.g. their orbital period distribution) can be compared fairly directly with the results of binary population synthesis codes and are therefore a strong test of population synthesis models for binary stars.

Following on from the detection of many binary stars by Maxted et al. \[Maxted et al. 2001\], we started a project to measure their orbits. The orbit of one of the new binary stars has been presented in Maxted et al. \[Maxted et al. 2002a\]. In this paper we present the orbits of a further 22 systems. We then consider the implications of the known sdB binary stars for their evolution. It should be noted that our sample is biased against sdB stars with G/K-type companions as the majority of our stars were selected from the PG survey which excludes most stars that show a Caii H-line.

## 2 Observations and Reduction

The data used in this study were taken with the Intermediate Dispersion Spectrograph (IDS) at the 2.5m Isaac Newton Telescope (INT) on the island of La Palma. Two different configurations of the IDS were used for the observations. The first setup consisted of the 500 mm camera with the R1200R grating centred in $`\text{H}\alpha `$ and the TEK (1kx1k) Charge couple device (CCD) giving a dispersion of 0.37 Å/pix and a resolution of 0.9 Å. The second setup used the 235 mm camera with the R1200B grating centred at 4350Å and the thinned EEV10 (2kx4k) CCD covering the Bamer lines from $`\text{H}\beta `$ to H$`ϵ`$, giving a dispersion of 0.48 Å/pix and a resolution of 1.4 Å. We carried out the observations during six different runs. The dates when the observations were taken and the setup used during each campaign are given in Table 1. We took two consecutive observations of each object and bracketed them with CuAr plus CuNe frames to calibrate the spectra in wavelength. We subtracted from each image a constant bias level determined from the mean value in its over-scan region. Tungsten flatfield frames were obtained each night to correct for the pixel to pixel response variations of the chip. Sky flatfields were also obtained to correct for the pixel to pixel variations of the chip along the slit. After debiasing and flatfielding the frames, spectral extraction proceeded according to the optimal algorithm of Marsh \[Marsh 1989\]. The arcs were extracted using the profile associated with their corresponding target to avoid systematic errors caused by the spectra being tilted. Uncertainties on every point were propagated through every stage of the data reduction.

## 3 Results

### 3.1 Radial velocity measurements

To measure the radial velocities we used least squares fitting of a model line profile. The model line profile is the summation of three Gaussian profiles with different widths and depths. For any given star, the widths and depths of the Gaussians are optimised and then held fixed while their velocity offsets from the rest wavelengths of the lines in question are fitted separately for each spectrum; see Maxted, Marsh & Moran \[Maxted et al. 2000c\] for further details of this procedure. For the data taken on the April 2000, the February, March, and May 2001 observing runs, the fitting was performed to the $`\text{H}\alpha `$ line whereas for the August and September 2001 observing campaigns, the fitting was carried out simultaneously to all the Balmer lines present in the spectra. Table 6 gives a list of the radial velocity measurements for the 22 objects presented in this paper.

Using the measured radial velocities of the lines, we determined the orbital periods of our targets. The description of the orbital period determination is given in Section 3.2. The results of folding the radial velocities of each object on its orbital period are plotted in Fig. 1. The error bars on the radial velocity points are, in most cases, smaller than the size of the symbol used to display them. For that reason we also display the residuals of the fit on a scale 10 times larger than the scale of the radial velocity curves. We find that there is no sign of ellipticity in any of the radial velocity curves, not even at long periods where departures from circular orbits might be expected. This is good evidence for the action of the common envelope. The values of the orbital periods, systemic velocities and radial velocity semi-amplitudes for each system are given in Table 2.

### 3.2 Determination of orbital periods

We use a “floating mean” periodogram to determine the periods of our targets (e.g. Cumming, Marcy & Butler 1999). The method consists in fitting the data with a model composed of a sinusoid plus a constant of the form:

$$\gamma +\mathrm{K}\mathrm{sin}(2\pi f(tt_0)),$$

where $`f`$ is the frequency and $`t`$ is the observation time. The key point is that the systemic velocity is fitted at the same time as $`K`$ and $`t_0`$. This corrects a failing of the well-known Lomb-Scargle \[Lomb 1976, Scargle 1982\] periodogram which starts by subtracting the mean of the data and then fits a plain sinusoid; this is incorrect for small numbers of points.

We obtain the $`\chi ^2`$ of the fit as a function of $`f`$ and then identify minima in this function. Section 3.3 gives a detailed explanation on the probability of the periods obtained being incorrect.

Table 2 gives a list of the orbital parameters derived for each sdB binary star. The orbital period of the second best period is also given, along with the difference in $`\chi ^2`$ between the two best periods found. The resulting periodograms ($`\chi ^2`$ versus orbital frequency) are given in Figs. 2 and 3. Each panel includes an blow up of the region in frequency where the minimum $`\chi ^2`$ is found. It is clear from the figures that in all cases, apart from PG1043+760, the difference in $`\chi ^2`$ between the best and the second alias is at least 10. We have made an exception for PG1043+760, because in this case the competing aliases are so close (owing to 1 cycle/year aliasing) that it makes more sense to consider the nearest competing group of aliases.

### 3.3 The probability that our periods are incorrect

A sometimes frustrating characteristic of radial velocity work is that while one can very soon know for sure that a star is binary – perhaps after just two measurements – it can take much longer to pin down the orbital period. One somehow has to know when the orbital period is “correct”. To compound this problem, another feature of radial velocity orbits is that if one picks an incorrect alias, the period can be *completely* wrong, even when the quoted uncertainty on the best-fit period is tiny. This is simply because the statistics are not Gaussian so that an error of 100 or even 1000 times quoted uncertainty can happen. Perhaps the most common way around this issue is the “method of overkill” where one takes so much data as to put the issue beyond any doubt, but this is necessarily inefficient. This is a particular problem with the sdB stars where there is no shortage of potential targets, but always a shortage of telescope time in which to observe them. Our approach whilst observing was to use the rule-of-thumb that an orbit was determined once the best-fit orbit improved upon the next-best by at least 10 in $`\chi ^2`$. The basis for this is that the probability of a period in the Bayesian sense is dominated by the term $`\mathrm{exp}\chi ^2/2`$ (see the appendix), and so a difference of more than 10 shows that the second-best alias is at least $`\mathrm{exp}5150`$ times less probable than the best. This rule-of-thumb was propagated into the submitted version of the paper, but the referee made two cogent criticisms of this approach. First, while the peak of the second alias may be $`>150`$ times less probable than the peak of the best alias, there is no guarantee that the total probability of *any* other period was as low. Second, some of our $`\chi ^2`$ values were larger than expected given the number of data points, suggesting some extra source of uncertainty that if included would reduce one’s faith in the best alias.

Both of these criticisms were justified, and the first in particular prompted us to develop a more rigorous approach which we describe below. To address the second point we have computed the level of systematic uncertainty that when added in quadrature to our raw error estimates gives a reduced $`\chi ^2=1`$. The reasoning behind this is that there may be an un-accounted source of error such as true variability of the star or slit-filling errors causing the poor fits of a few stars. It seems unlikely that such errors will be either correlated with the orbit or with the statistical errors we estimate, and therefore we add a fixed quantity in quadrature with our statistical errors as opposed to applying a simple multiplicative scaling to them. In all cases we use a minimum value of $`2\mathrm{km}\mathrm{s}^1`$ corresponding to 1/10th of a pixel which we believe to be a fair estimate of the true limits of our data.

The end result is that in two cases there is indeed a relatively large probability that our periods are in error (KPD0025+5402; PG1032+406), nevertheless we continue to list them in this paper on the bases that (a) the best periods remain strongly favoured in all cases and therefore are probably correct and (b) the IDS spectrograph is being withdrawn from service and it is not clear when we will get the opportunity to gather more data. These systems (KPD0025+5402, PG1032+406) therefore come with a health warning: there is a larger than desirable probability that we may not have picked the best alias for them.

In most cases the probability of the orbital period being further than 1 and 10 per cent from our favoured value is the same. This is because all the significant probability lies within a very small range around the best period, with all the significant competition (i.e. next best alias) placed outside the 10 per cent region around the best alias. In the case of HD171858, the probability of the true orbital period being further than 1 per cent from the one given in this paper is very high. This is due to the short baseline of our observations (1 day) which means that the period is not determined anything like as well as the others and in particular is not confined to within 1 per cent of the best fit. For the purpose of comparing observations to evolutionary models there is no need for better than 10 per cent measurements and we are quite certain that the period is correct within those limits. The probabilities for all systems are listed in Table 3 where we give the logarithm (base 10) of the chance that the true period lies further than 1 and 10 per cent of our favoured value.

To calculate these probabilities, we integrated equation A4 from the appendix of Marsh, Dhillon & Duck \[Marsh, Dhillon & Duck 1995\]. Full details of our method are given in the appendix to this paper. The calculation is Bayesian, and therefore involves prior probabilities over all parameters, i.e. the systemic velocity, semi-amplitude, phase and period. These “priors” are uncertain and therefore the final probabilities listed in Table 3 are similarly uncertain. This we believe is the main reason why one would like the probability of being incorrect to be *very* small, so as to overcome any plausible uncertainty in the prior probabilities; this is one justification for the method of overkill. For most of our sample, this is in fact the case, but we feel that the uncertainties are large enough to leave a grain of doubt when the probability listed in Table 3 rises above $`0.1`$%, or $`3`$ in the log, as it does for the two stars discussed.

## 4 Discussion

### 4.1 Effective temperature, surface gravity and helium abundance

For those sdBs that we observed in the blue, we measured the effective temperature, Teff, the surface gravity, $`\mathrm{log}g`$, and the helium abundance, $`\mathrm{log}(\mathrm{He}/\mathrm{H})`$. We used Saffer et al.’s \[Saffer et al. 1994\] procedure to fit the profiles of the Balmer, the He i and the He ii lines present in the spectra by a grid of synthetic spectra. The synthetic spectra obtained from hydrogen and helium line blanketed NLTE atmospheres \[Napiwoztki 1997\] were matched to the data simultaneously. For stars cooler than 27 000 K we used the metal line-blanketed LTE model atmospheres of Heber, Reid, & Werner \[Heber, Reid & Werner 2000\]. Before the fitting was carried out, we convolved the synthetic spectra with a Gaussian function to account for the instrumental profile. See Heber et al. \[Heber, Reid & Werner 2000\] for details of the models. The values of Teff, $`\mathrm{log}g`$, and $`\mathrm{log}(\mathrm{He}/\mathrm{H})`$ for some of the objects are the same as those presented by Maxted et al. \[Maxted et al. 2001\]. Those that are different from the ones obtained by Maxted et al. \[Maxted et al. 2001\] and those for which Teff and $`\mathrm{log}g`$ have not been calculated previously (KPD0025+5402 and KPD1946+4340) are given in Table 4.

All the objects from Table 4 apart from two lie in or near the band defined by the zero-age extreme horizontal branch (ZAEHB), the terminal-age extreme horizontal branch (TAEHB) and the He main sequence (HeMS) and are therefore extreme horizontal branch stars (EHB). (See Fig. 2 of Maxted et al. \[Maxted et al. 2001\] for a $`\mathrm{log}g`$ versus Teff plot for the objects not included in Table 4.) The two objects that lie outside the EHB band are post-EHB stars. PG0839+399 appears in Saffer et al. \[Saffer et al. 1994\] as an EHB star, and KPD1946+4340 is a new post-EHB star.

### 4.2 Orbital parameters known up to now

Table 5 gives a list of all of the orbital parameters of sdB binary stars known to date. In each case we require that both the orbital period and the radial velocity semi-amplitude are measured. We combine the orbital periods and the radial velocity semi-amplitudes to calculate the mass function, $`f_m`$, of the system according to the well-known relation:

$$f_m=\frac{M_2^3\mathrm{sin}^3i}{(M_1+M_2)^2}=\frac{PK_1^3}{2\pi G},$$

where the subscript “$`1`$” refers to the sdB star and “$`2`$” to its companion.

If we take a canonical mass of $`0.5\text{M}_{}`$ for the sdB star, we can also calculate the minimum mass of its companion, $`M_{2\mathrm{m}\mathrm{i}\mathrm{n}}`$. The values for $`f_m`$ and $`M_{2\mathrm{m}\mathrm{i}\mathrm{n}}`$ obtained in each case are given in Table 5. We have also added in the table a column indicating the nature of the companion, where known, i.e. whether it is a compact object (most likely a white dwarf, indicated by “WD”) or a non-degenerate object (a main sequence star or a brown dwarf, indicated by “MS”). The evolution of each type is fundamentally different since in the first case the system must go through at least two mass transfer episodes whereas in the second case the system suffers only one mass transfer episode. It is notable that all the sdB stars with “MS” companions have very short orbital periods.

In Fig. 4 we present a histogram of all the known orbital periods of sdBs. The dashed line shows the systems with previously published periods and the solid line shows the combination of the previously published periods and the ones measured in this work. There is no sign of any fine structure such as the well-known “period gap” of the cataclysmic variable stars. The main points to take from this plot are the large dynamic range and that our survey has substantially increased the numbers of long period systems ($`P>`$ a few days). This is probably a consequence of the large time base of our data.

Fig. 5 shows $`M_{2\mathrm{m}\mathrm{i}\mathrm{n}}`$ versus orbital period, assuming $`M_1=0.5\text{M}_{}`$. All the systems given in Table 5 are included in the plot. The systems whose periods are measured in this paper are plotted as asterisks whereas sdB binaries with previously published orbital periods are plotted as plus signs. In addition, to test the idea of mass loss followed by a common envelope shortly before the helium flash, we present limits based upon a common envelope phase initiated when the progenitor of the sdB star was at the tip of the red giant branch. We calculate the effect of the common envelope in the standard fashion \[Webbink 1984\] with a fraction $`\alpha _{\mathrm{CE}}`$ of the loss of orbital energy set equal to the binding energy of the envelope. We apply the standard mass transfer equations instead of those calculated by Nelemans et al. \[Nelemans et al. 2001\] because in most cases we are treating the second phase of mass transfer instead of the first phase. The change in orbital energy is

$$\mathrm{\Delta }E_{\mathrm{orb}}=\frac{GM_1M_2}{2a_i}\frac{GM_{\mathrm{sdB}}M_2}{2a_f},$$

where $`M_1`$ is the mass of the sdB star’s progenitor and $`a_i`$ and $`a_f`$ are the initial and final orbital separations. Parameterising the envelope binding energy as

$$E_{\mathrm{env}}=\frac{GM_1(M_1M_{\mathrm{sdB}})}{\lambda R_1},$$

where $`\lambda `$ depends upon the structure of the envelope and the contribution of internal as well as gravitational energy \[Dewi & Tauris 2000\]. We then have

$$\alpha _{\mathrm{CE}}\lambda \left(\frac{GM_{\mathrm{sdB}}M_2}{2a_f}\frac{GM_1M_2}{2a_i}\right)=\frac{GM_1(M_1M_{\mathrm{sdB}})}{R_1}.$$

Thus given $`M_1`$, $`M_2`$, $`M_{\mathrm{sdB}}`$, $`R_1`$ and the combination $`\alpha _{\mathrm{CE}}\lambda `$, the final separation can be computed for any given initial separation.

Following this formalism, the different lines in the plots represent the following constraints:

1. the leftmost curve represents the limit imposed by the lack of mass transfer in these binaries. The limit has been calculated assuming that the companion is a main sequence star and therefore does not apply for white dwarf companions. We have used Eggleton’s \[Eggleton 1983\] formula for the Roche lobe radius to obtain this limit.
2. the curve marked with $`M_1>1.9\text{M}_{}`$ marks the limit beyond which there would be no helium flash. We have employed Hurley et al.’s \[Hurley et al. 2000\] analytic formulae to calculate how large the progenitor of the sdB was, assuming that mass transfer occurred near the tip of the RGB. For a given initial mass, metallicity and common envelope efficiency one can then predict the current period as a function of $`M_2`$, assuming that the sdB star has a mass of $`0.5\text{M}_{}`$ at the end. Massive companions get rid of the envelope easily and therefore end up at long periods. To get short periods and high companion masses one needs some combination of inefficient envelope ejection, a small red giant and a massive envelope. The latter two conditions occur for large mass progenitors since the radius at the tip of the giant branch decreases with mass for low mass stars. This sets a lower limit on the progenitor mass; $`M_1=1.9\text{M}_{}`$ is taken as the maximum since above this the core is not degenerate and there is then no natural explanation for the $`\mathrm{log}g`$/$`T_{\mathrm{eff}}`$ distribution of sdB stars.
3. the vertical line on the right hand side is the upper limit on the final period given that the sdB progenitor mass had to be $`1\text{M}_{}`$ for it to have evolved to this stage within the lifetime of the Galaxy.
4. the maximum value for $`M_2`$ has been chosen to be $`1.4\text{M}_{}`$, corresponding to the Chandrasekhar stellar limit. Any main sequence star this massive would be easily visible.
5. the dashed line that appears at the bottom of the figures indicates the minimum $`M_2`$ required to get the smallest amplitude that we can detect ($`10\mathrm{km}\mathrm{s}^1`$). This crudely shows that selection effects are likely to be weak in our data.

In addition to show the sensitivity to variations of these limits, we show the case for $`M_1<1.75\text{M}_{}`$ and $`M_2<0.65\text{M}_{}`$. Most of the data are contained within these more restrictive limits. It should be remembered when looking at these plots that the companion masses are lower limits owing to the unknown orbital inclinations of most targets.

Three panels are shown to show the effect of varying either the metallicity \[Fe/H\] (which affects the radius of the star) or the common envelope efficiency/envelope structure parameter, $`\alpha _{\mathrm{CE}}\lambda `$. Increasing the metallicity increases the radius at mass transfer and results in longer period systems being formed. This is also the result if the efficiency with which orbital energy is used to expel the common envelope increases. The largest value of $`\alpha _{\mathrm{CE}}\lambda `$ used in our calculations is 1.1. If the value of $`\alpha _{\mathrm{CE}}\lambda =`$ 2, favoured by Nelemans et al. \[Nelemans et al. 2001\], is considered, this results in even longer period systems being formed.

The conclusion to be drawn from all this is that although a common envelope phase as a result of mass loss immediately prior to the helium flash covers the range of parameters of most systems, it does not explain simultaneously systems of long orbital period and systems of short orbital period with massive companions. The two systems KPD1930+2752 ($`P=0.0951\mathrm{d}`$, $`M_2>0.97\text{M}_{}`$) and KPD0422+5421 ($`P=0.0902\mathrm{d}`$, $`M_2>0.50\text{M}_{}`$) look particularly discrepant and may require an alternative formation path, such as a descent from progenitors with $`M_1>2\text{M}_{}`$.

### 4.3 sdB stars with main-sequence companions

Five of the short period sdB binary stars are known to have main-sequence or brown dwarf companions. All of them have very short periods indeed, with PG1329+159 having the longest period at $`4.0\mathrm{hr}`$. This fits with the low masses of their companions, $`0.1\text{M}_{}`$ compared to $`0.5\text{M}_{}`$ for a typical white dwarf. It also suggests an interesting possibility for the $`10`$$`20`$ percent of sdB stars that are apparently single (Maxted et al. 2001; Green et al. 2000): perhaps these stars could be the result of merging of even lower mass companions, unable to survive the common envelope phase. This would at least avoid the need to have a single star route for sdB formation to supplement the binary star route which already seem to account for 80 to 90 percent of these stars. The long period sdB/main-sequence binaries identified by Green et al. \[Green, Liebert & Saffer 2000\] must have avoided a common envelope. Perhaps they could do so because their relatively high mass main-sequence components were able to accrete at a high enough rate to avoid the common envelope. It would be of great interest to measure element abundances in the main-sequence components of these stars, e.g. abundance ratios of 12C/13C \[Sarna et al. 1995\].

## 5 Conclusions

We have confirmed the binary nature of 22 subdwarf B stars and have measured their orbital parameters. This work increases the sample of sdB binaries with known orbital parameters to 38. The observations extend over several months allowing us to detect orbital periods of the order of tens of days, longer than any previously measured.

We have measured Teff, $`\mathrm{log}g`$, $`\mathrm{log}(\mathrm{He}/\mathrm{H})`$ for the sdBs where previous measurements of these quantities did not exist. When we place the results in the Teff-$`\mathrm{log}g`$ plot we find that two out of the 22 sdBs are post-EHB stars.

The large range of orbital periods and companion masses is a challenge to simple theories for the formation of sdB stars. Although binary-induced mass-loss at the tip of the red giant branch is able to explain most systems, it appears unlikely to be the only formation route. Full population synthesis will be needed to establish the viability of alternative paths.

Amongst sdB stars with known orbits, those with low-mass main-sequence or brown dwarf companions have particularly short periods. It seems likely that a fraction of such systems, particularly those of very low companion mass, may not have survived the common envelope phase. We suggest that these could now be the single sdB stars.

## Acknowledgements

LMR was supported by a PPARC post-doctoral grant. The reduction and analysis of the data were carried out on the Southampton node of the STARLINK network. The Isaac Newton Telescope is operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. We thank PATT for their support of this program.

## Appendix A The probability of being wrong

Marsh, Dhillon & Duck \[Marsh, Dhillon & Duck 1995\] derived the following equation (their A4) for the probability of a binary star of orbital frequency $`f`$ versus a single star given a set of radial velocity data, $`D`$

$`{\displaystyle \frac{P(B,f|D)}{P(S|D)}}`$ $`=`$ $`{\displaystyle \frac{P(B,f)}{P(S)}}{\displaystyle \frac{4\left(2\pi \right)^{3/2}}{R_K^2R_\gamma \left(det𝐀\right)^{1/2}}}{\displaystyle \frac{R_\gamma \left(_iw_i\right)^{1/2}}{\left(2\pi \right)^{1/2}}}`$
$`\times \mathrm{exp}{\displaystyle \frac{1}{2}}\left(𝐛^t𝐀^1𝐛{\displaystyle \frac{\left(_iw_iV_i\right)^2}{_iw_i}}\right),`$
where the vector $`𝐛`$ is given by

$$𝐛=\left(\begin{array}{c}_iw_i\\ _iw_ic_i\\ _iw_is_i\end{array}\right),$$
(1)
and the matrix $`𝐀`$ is

$$𝐀=\left(\begin{array}{ccc}_iw_i& _iw_ic_i& _iw_is_i\\ _iw_ic_i& _iw_ic_i^2& _iw_ic_is_i\\ _iw_is_i& _iw_is_ic_i& _iw_is_i^2\end{array}\right),$$
(2)
with $`w_i=1/\sigma _i^2`$, where $`\sigma _i`$ is the uncertainty on the $`i`$-th point and $`c_i=\mathrm{cos}(2\pi ft_i)`$, $`s_i=\mathrm{sin}(2\pi ft_i)`$ where $`t_i`$ is the time of the $`i`$-th point; $`V_i`$ is the measured radial velocity for the $`i`$-th point. The term inside the brackets of the exponential is equal to the difference in $`\chi ^2`$ between adopting a constant velocity model versus a model of a constant plus a sinusoid, $`\chi _c^2\chi _s^2`$. That is, it is a measure of the gain one makes by adding a sinusoid to a constant model. Since the constant model $`\chi _c^2`$ is independent of frequency, this is the justification for our statement in the main text that the probability of a certain period is dominated by the term $`\mathrm{exp}\chi ^2/2`$.

Already some assumptions have gone into this result which is the end product of integrating over the systemic velocity and the sine and cosine amplitudes for a given frequency. In particular we took the prior probabilities over these parameters to be uniform in a box extending from $`R_\gamma `$ to $`+R_\gamma `$ in systemic velocity, and $`R_K`$ to $`+R_K`$ in each of the semi-amplitudes (a factor 4 has entered into Eq. A compared to the original version from Marsh, Dhillon & Duck \[Marsh, Dhillon & Duck 1995\] to account for a factor 2 change in our definition of these ranges). We continue to adopt this prior, largely because it then leads to a tractable integral, but also because it is hard to justify any particular choice of prior, and so one might as well adopt the simplest (although see below for more on this).

Eq. A is correct, given our assumptions, as long as the volume of integration wholly encloses the integrand in the systemic velocity/semi-amplitude space. The integrand is a 3D Gaussian centred upon the best-fit values. For good enough data, this will occupy a small region, comparable in dimension to the uncertainties quoted on these parameters. At least this is the case for the best-fit frequency, but it is less clear that it will always be the case. To account for the case where a significant fraction of the integrand lies outside the integration region, we wrote a program that calculated the value of Eq. A and, in addition, a correction to it from a Monte Carlo integration. Essentially this computes the fraction of the integrand lying outside the volume of integration. A detail was that we changed the volume of integration to a cylinder with axis along the systemic velocity axis and radius $`R_K`$, which gives a uniform prior in orbital phase as one expects. The Monte Carlo work was carried out by computing the eigenvalues and eigenvectors of the matrix $`𝐀`$ and then generating points scattered around the maximum using Gaussian random number generators to define the multiple of each eigenvector to use. The correction factor was then the ratio of the number of points that fell within the volume of integration to the total number of points generated. In practice our data was good enough that none of this made much difference to the results.

We adopted a value of $`R_\gamma =500\mathrm{km}\mathrm{s}^1`$ and set $`R_K`$ by taking the maximum companion mass to be $`1.5\mathrm{M}_{}`$. There is certainly a case for having a prior in the systemic velocity that falls off more gradually from a peak close to zero. A Gaussian for example could still be accommodated analytically by a small modification of the matrix $`𝐀`$, however we feel that the greater ease of understanding of a uniform prior is an advantage that should not be wasted. What this means in practice is that we will over- or under-estimate the probabilities listed in Table 3 according to whether the best-fit systemic velocity is far from or close to zero velocity. This is a quantitative manifestation of the lack of confidence one would have in an orbit with a systemic velocity of, say, $`2000\mathrm{km}\mathrm{s}^1`$ compared to a more moderate one. Similar remarks apply to the prior for $`K`$. The uncertainty over the priors is why the final probabilities are themselves subject to uncertainty.

The final step was to integrate the probability as a function of frequency over ranges of frequency. The only technical difficulty here is that the integrand over frequency is extremely peaky with large ranges of very low probability pierced by narrow spikes of high probability. There is little one can do about this apart from computing the integrand over a very finely spaced series of frequencies. We applied straightforward trapezoidal integration, halving the spacing of the frequency grid until the integral converged. Typically this required of the order of one million points, but did not take overly long on any one object. The same points about the systemic velocity and semi-amplitude priors applies to the prior over frequency. We took a prior uniform in frequency, but this has no particular justification. The final probability is uncertain to the extent that this prior deviates from reality. A factor of 10 would not especially surprise us here.


text:
# Untitled Document

On the Minimum Energy Configuration of a Rotating Barotropic Fluid:

A Response to Narayan & Pringle astro-ph 0208161

Sébastien Fromang (IAP) & Steven A. Balbus (UVa)

In a recent posting, Narayan & Pringle (2002, astro-ph/0208161) criticized a paper written by the current authors (Fromang & Balbus 2002, astro-ph/0207561 hereafter FB02). FB02 presented a second order variational calculation and concluded that supersonically rotating barotropic fluids had lower energy states neighboring that of equilibrium. (Note that this is not the same as claiming instability, a point noted explicitly in FB02.) Narayan & Pringle challenged this, presented a specific class of equilibrium solution, and claimed to prove that it represented a global energy minimum. They concluded that FB02 was incorrect.

The calculation presented in FB02 is indeed incorrect. The difficulty is that the second order changes in the energy cannot be calculated directly from the energy itself, even when the variations are restricted to satisfy the constraints of mass and angular momentum conservation. Rather, the constraints themselves change the effective functional form of the energy, by forcing additional second order correlations that would otherwise be absent. Narayan & Pringle (2002) is therefore correct in stating that our second order variational calculation is in error.

Nevertheless, we are compelled to respond to the critique of Narayan & Pringle (2002). Our reasons are twofold: (1) the critique unfairly misrepresents statements in FB02; and (2) the main discussion presented in the critique appears itself to have flaws. In particular, the conclusion of Narayan & Pringle that their chosen example is a global minimum with respect to arbitrary variations from equilibrium does not follow from the argument presented in their paper. We very briefly discuss these two points in turn.

1. Lower energy states and instability. The existence of lower energy states neighboring the equilibrium solution is not a sufficient criterion for instability. This point is stated explicitly in FB02:

“Note that our procedure does not ensure instability in the latter case \[supersonic rotation\], it simply shows that neighboring states of lower energy exist.”

Unfortunately, Narayan & Pringle misrepresent FB02, stating, for example:

“Since FB02 claim that the only stable barotropic configurations are those that rotate uniformly and are subsonic at all radii, $`\mathrm{}`$

Throughout the critique, the existence of lower energy states is referred to as a “stability criterion.” No such claim was ever made in FB02. A stability analysis of a uniformly rotating polytrope in an external potential need not show any unstable modes, any more than a stability analysis of, say, a Keplerian disk is bound to show linear instability, even though such a profile is not an energy minimum. When a weak magnetic field is added to a Keplerian disk, no matter how small the field energy, pathways are opened to states of lower disk energy and the MRI ensues, even for purely axisymmetric disturbances. In the case of a rotating polytrope, the question is not whether the system itself is directly unstable per se, it is whether subsequent evolution drives the initial state further away from its original equilibrium when one makes energetically negligible modifications — say a weak magnetic field and a tiny angular velocity gradient.

2. Global Energy Minima. Narayan & Pringle analyze rotating polytropes in the fixed external potential

$$\mathrm{\Phi }=\frac{1}{2}\mathrm{\Omega }_K^2(R^2+z^2)$$

where $`\mathrm{\Omega }_K`$ is a constant and $`R`$ and $`z`$ are cylindrical radius and vertical coordinate respectively. For a given mass $`M`$ and angular momentum $`J`$, it is shown that there is a unique equilibrium solution. Combining this result with the positive-definite character of the energy $`E`$, Narayan & Pringle conclude that the energy extremum represented by their solution must be a global minimum.

Clearly something more is required to establish a global minimum. This line of reasoning is insufficient to establish the existence of a local minimum, let alone a global one. (A gaussian function is everywhere positive-definite and possesses a unique extremum, which is, of course, a maximum.) We do not wish to belabor this issue or other technical difficulties with the presentation, since the point of discrediting the variational calculation is moot. The point of physics, however, is not without interest.

Consider the gravitational potential arising from a constant density sphere of large but finite radius. Within the sphere, the potential is exactly of the above form (perhaps with an inconsequential additive constant), while outside the sphere the potential is Keplerian. There will be values of $`M`$ and $`J`$ for which the unique equilibrium solution lies entirely in the sphere, and is of the precise form considered by Narayan and Pringle. The total energy is bounded from below (and can thus be made positive by adding a trivial constant to the potential). But for any equilibrium solution lying entirely within the sphere, there exists another state with the same $`M`$ and $`J`$ that has a lower energy: one in which the mass distribution is essentially identical, and the angular momentum is removed to a ring of arbitrarily small mass and arbitrarily small total energy at arbitrarily large distances. (Supersonic disks do in fact avail themselves of these lower energy solutions.) Note that if the potential rises sufficiently steeply at large distances, the argument need not hold. Thus, the nature of the potential at large $`r`$ plays a role in establishing the existence of global minima. If the solutions of Narayan & Pringle (2002) are in fact true global minima, this is dependent upon the rapid growth of the harmonic potential at large distances, and is not a simple consequence of the existence of a unique extremum and positivity of the energy.


text:
# Atom Lithography using MRI-type feature placement

## Abstract

We demonstrate the use of frequency-encoded light masks in neutral atom lithography. We demonstrate that multiple features can be patterned across a monotonic potential gradient. Features as narrow as 0.9 $`\mu `$m are fabricated on silicon substrates with a metastable argon beam. Internal state manipulation with such a mask enables continuously adjustable feature positions and feature densities not limited by the optical wavelength, unlike previous light masks.

Atom lithography has been used to make lines as small as 13 nm BellSmall ; CrSmall , to fabricate features with aspect ratios greater than 2:1 rhese ; RIE , and to pattern areas as large as 38 cm2 kjapl . Atoms patterned with light masks are uniquely well-suited for structured doping doping  and for fabrication demanding long-range coherence nanoruler . However, feature spacings in light masks have been limited by the length scale of the optical wavelength used timp-prl ; jabez-science ; siuau-al ; drodofsky ; lison ; SWQ ; 2D ; engels ; mutzel ; brezger ; jabezover8 . In this Letter, we demonstrate a technique that extends light masks beyond the feature density of optical lithography and to more complex patterns than standing wave interferences.

Research in neutral atom lithography has been motivated by several considerations. Neutral atomic beams with thermal kinetic energies have short ($`<.1`$ nm) de Broglie wavelengths, minimizing diffractive resolution limitations. Neutral atoms are insensitive to stray electric and magnetic fields, and long-range inter-particle interactions are weak. A wealth of atom optics (mirrors, guides, etc.) have been developed, including light masks, which, unlike physical masks, do not sag and cannot be clogged or damaged. Laser-accessible internal structure of atoms allows efficient cooling, the storage of internal energy, and complex manipulations. State-sensitive atom lithography with metastable atoms has been demonstrated with both positive- and negative-tone resists capable of sub-10-nm resolution karl-sam ; kjapl ; passH .

In this work we encode the desired position distribution using the frequency, instead of the intensity, of the light. In the presence of an appropriate potential gradient, the optical resonant frequency of the atom will be position-dependent. Thus each spectral component of a laser beam interacts can interact with the atoms at a specific but separate location. The frequency encoding of spatial information is commonly known from magnetic resonance imaging (MRI), where protons (instead of atoms) are imaged (instead of patterned) in a potential gradient by radio- (instead of optical-) frequency radiation. We will refer to atom lithography with frequency-encoded light masks as Atomic Resonance Lithography (ARL). Thomas and coworkers have developed neutral atom measurement techniques based on a more direct analogue to MRI thomas\_th ; thomas\_exp .

Frequency encoding has several promising features. Just as sub-100-$`\mu `$m features can be resolved with 100-meter wavelength radio waves in MRI, far sub-optical (ie, sub-micron) features and spacings can be created with optical resonances in atoms thomas\_th ; thomas\_exp ; olshanii . As shown in Ref. olshanii , arbitrary patterns can be frequency-encoded and generated. Furthermore, ARL can create two-dimensional patterns, just as MRI can image two- lauterbur  or three-dimensional distributions. Finally, frequency-encoded masks do not demand high transverse coherence of the incident atoms. By contrast, although intricate holographic patterns fujita ; morinaga ; shimizu  have been detected, the required coherent sources lack the flux to pattern a surface: experiments have used integrated surface densities on the order of $`10^9`$ below a typical lithographic dose.

Figure 1 describes a realization of ARL based on quenching directly within the potential gradient. First, metastable argon atoms are optically pumped to the $`|J=2,m_J=2`$ sublevel of the $`4s[3/2]_2^o`$ or $`1s_5`$ state. Not shown is earlier optical quenching of the atoms in the $`4s^{}[1/2]_0^o`$ or 1s3 state. If unquenched, these “J=0” atoms would contribute a background rate of resist formation. Next, the atoms pass through a zone of two overlapping beams. A resonant $`\sigma ^+`$ “shift” beam (S) on the 812 nm $`|J=2,m_J=2`$$`|J^{}=3,m_J^{}=3`$ cycling transition overlaps with a $`\sigma ^{}`$ “quench” beam (Q) containing four frequency components, $`\pm \delta _1`$ and $`\pm \delta _2`$. As shown in Fig. 2a, these frequencies are resonant at two stark shifts (i.e., four positions) in the resonant shift beam. At these four positions, atoms are optically pumped by Q to the true atomic ground state $`|g`$. Note that two frequencies (e.g., $`\pm \delta _1`$) were required for complete optical pumping at each resonant position because atomic population existed at both the strong- and weak-field-seeking dressed states. This internal state manipulation creates a pattern on a surface by exploiting the internal energy difference between the the ground state and the $`4s[3/2]_2^o`$ metastable state of argon: atoms in $`|m`$ release 11.5 eV of internal energy to activate the formation of a resist on the surface (as described below), while atoms in $`|g`$ do not affect the surface. Thus, the detunings $`\delta _1`$ and $`\delta _2`$ control where features – here, stripes of unprotected surface – will appear.

Metastable argon (Ar\*) atoms exposed the substrate through the light mask with a flux density of $`3\times 10^{12}`$ cm-2 s-1. Because the system is evacuated with diffusion pumps, surfaces within the system are covered with several monolayers of siloxane hydrocarbons. When the metastable atoms hit the surface, they transfer their internal energy to the physisorbed hydrocarbons and induce a chemical change that results in the formation of a durable material. The resist material remains on the surface even after exposure to air and/or solution, and can thus be used as a mask for etching. The beam apparatus and lithographic process used in this work have been described in more detail elsewhere kjapl .

A gold substrate was exposed for five hours through the optical mask of Fig. 1. The light mask parameters were as follows. The S beam had 5.9 mW cylindrically focused to a profile of 46 $`\mu `$m by 925 $`\mu `$m, creating a maximum stark shift that was spectroscopically measured to be 150 MHz. Cylindrically symmetric 110-$`\mu `$m-wide quench beams Q1 and Q2 each had 70 $`\mu `$W of power, split between the pair of detunings, $`\pm \delta _1=\pm 2\pi \times 40`$ MHz and $`\pm \delta _2=\pm 2\pi \times 110`$ MHz. After exposure, the substrate was removed from the vacuum system and etched for 7 minutes in a ferricyanide solution.

Figure 2b shows an optical micrograph of the etched substrate. At positions where either Q1 or Q2 was resonant, metastable atoms were quenched, leaving the underlying substrate unprotected against the subsequent gold etch. At all other positions, a resist material was formed as described above, and gold remained on the silicon substrate after etching. A pattern was formed across 2 cm of the substrate. Figure 2b was taken in the center of that pattern, where atoms entered perpendicular to the light mask. The comparison (indicated by dotted arrows) to Fig. 2a shows that, as desired, there are four features: two on each side of the Gaussian intensity profile of the S beam. We fit the sum of four Lorentzians to the integrated reflectivity profile of the substrate to determine more precise feature locations and positions. The lower (upper) two lines have a half-width of 10$`\pm `$1 $`\mu `$m (14$`\pm `$1$`\mu `$m) and a center-to-center separation of 26$`\pm `$2$`\mu `$m (34$`\pm `$2 $`\mu `$m). The asymmetry between the upper and lower set of lines was not expected, but could be explained by aberrations in the cylindrical optics that formed the shift beam. The lower set of lines are within experimental error of the expected linewidth of atomic density, 9 $`\mu `$m (as discussed below), and the expected feature separation, 27 $`\mu `$m. Note that the transfer function from atomic flux density to reflectivity of the surface after etching is nonlinear, such that purely atomic calculations may vary systematically from the observed features.

The significance of the result shown in Fig. 2 is that two lithographic features were created across each monotonic part of the potential. Unlike previous atom lithographic results with light masks, the wavelength of light does not determine the minimum feature separation. The ultimate resolution of ARL can be limited by a variety of effects, discussed Refs. thomas\_th ; olshanii  and summarized in the following paragraphs.

The most fundamental limitation is the spectroscopic precision with which each resonant transfer position is defined, $`\delta x_{\mathrm{sp}}=\mathrm{}\mathrm{\Gamma }^{}/F`$, where $`\mathrm{\Gamma }^{}`$ is the spectroscopic width of the transition, and $`F`$ is the potential gradient. In these experiments we used laser light to create the potential gradient, but gradients could also be created with magnetic thomas\_exp  or electric fields. For the realization of ARL discussed above, $`\mathrm{\Gamma }^{}`$ is the power-broadened and time-broadened quenching transition, with a natural linewidth $`\mathrm{\Gamma }=2\pi \times 5`$ MHz. Experimentally, $`\mathrm{\Gamma }^{}=2\pi \times 20`$ MHz. For the gradient $`F/\mathrm{}=2\pi \times 2.3`$ MHz/$`\mu `$m at the feature locations in Fig. 2, we find $`\delta x_{\mathrm{sp}}=9`$$`\mu `$m. Note that for a Raman transition, as will be considered below, $`\mathrm{\Gamma }^{}=v_L/w_0=1/\tau `$, where $`v_L=850`$ m/s is the longitudinal velocity of the Ar\* atoms, $`w_0`$ is the 1/e2 intensity waist of the Raman beams, and $`\tau `$ is the atom-light interaction time.

Atomic motion during the transfer can also limit the resolution. While transferring an atom of mass $`M`$ to a state with a gradient $`F`$, acceleration and diffraction limit the resolution to the order of $`F\tau ^2/2M`$, where $`\tau `$ is again the transfer time thomas\_th . However, the numerical factors depend on which states are in a gradient and for how long. In the case where $`\mathrm{\Gamma }^{}`$ is determined by the interaction time of a Raman transfer, and the pattern is formed by the atoms transferred to the gradient, one can choose an optimal $`\tau `$ to give a resolution of $`\delta x_{\mathrm{}}=2\left(\mathrm{}^2/2MF\right)^{1/3}`$. For 40Ar and $`F/\mathrm{}=2\pi \times 1`$ MHz/$`\mu `$m, the optimal size is $`\delta x_{\mathrm{}}=100`$ nm. This limit is not, however, fundamental: Olshanii et al. show that an appropriate choice of probe frequency variation during atomic transfer (a “magic phase”) can correct for later evolution, such that the spectroscopic resolution $`\delta x_{\mathrm{sp}}`$ is recovered olshanii .

Free flight between the pattern formation and the substrate will allow the atoms to diffract further. Given a free-flight time $`T_{FF}`$, the diffraction limit is $`\delta x_{\mathrm{dif}}=\sqrt{\mathrm{}T_{FF}/2M}`$. For instance, for $`T_{FF}=30\mathrm{\Gamma }^1`$, this limit is $`\delta x_{\mathrm{dif}}=27`$ nm, smaller than the other limits given here. Note that this is a quantitative restatement of one of the desirable qualities of atoms as a patterning constituent: their short de Broglie wavelength $`\lambda _{dB}`$ minimizes the diffraction, $`\sqrt{\lambda _{dB}D}`$, where $`D`$ is the mask-substrate separation.

The collimation of the atomic beam can limit the resolution of the pattern. Given an rms velocity spread $`v_{rms}`$, $`\delta x_\mathrm{v}v_{rms}(\tau +T_{FF})`$. For $`\tau +T_{FF}=30\mathrm{\Gamma }^1`$ and $`v_{\mathrm{rms}}=0.5`$m/s, $`\delta x_\mathrm{v}=0.4\mu `$m. Thus with a small (20 nm) transverse coherence length, sub-micron lithography can be performed.

In order to improve our resolution, we performed a second experiment with several experimental modifications. First, we improved the quality of the optics forming the S beam to reduce the waist by a factor of 3. Next, as shown in Fig. 3, we used a coherent transfer scheme by detuning the S beam by 2 GHz and introducing Raman beams. At resonant positions in the light mask, the Raman beams transfer atoms from $`|m1`$ to $`|m2`$. In the third interaction zone, only the atoms transferred to $`|m2`$ are quenched statefn . With such a scheme, we realized a factor of 10 improvement in resolution, as described below. Furthermore, such a coherent ARL scheme avoids momentum diffusion by spontaneous emission, and can implement further resolution improvements such as the coherent manipulations described in Ref. olshanii .

A silicon $`100`$ substrate with a 15 nm oxide layer was exposed for six hours through the ARL mask depicted in Fig. 3. The R beam of 25 $`\mu `$m by 103 $`\mu `$m had $`\sigma ^{}`$ and $`\sigma ^+`$ components with powers of 65 $`\mu `$W and 60 $`\mu `$W and detunings of $`+61.1`$ MHz and $`+104.2`$ MHz, respectively. A 1.3 G quantization field reduced the differential detuning by 5.4 MHz, such that the net Raman detuning was $`\mathrm{\Delta }/2\pi =37.7`$ MHz. The shift beam S was detuned by $`\delta _S/2\pi =2.00`$ GHz above resonance, had $`6.9\pm 0.3`$ mW of power, and a profile of 15$`\pm 1`$$`\mu `$m by 250 $`\mu `$m. The $`\sigma ^+`$ polarized S beam interacted with both $`|m1`$ and $`|m2`$, but gave a differential gradient $`F/\mathrm{}=2\pi \times 6.4`$ MHz/$`\mu `$m due to the difference in Clebsh-Gordan coefficients: 1 and $`\sqrt{2/5}`$, respectively (note that Fig. 3 shows only the transition $`|m1|e`$, for simplicity). The third interaction zone had a 350 $`\mu `$W quench beam (Q) with a waist of 210 $`\mu `$m, such that about 12 photons are scattered from atoms in $`|m2`$. Because $`|m1`$ was in the $`|2,+2`$, it was dark with respect to the $`\sigma ^+`$ polarized quench. Untransferred atoms arrived at the substrate in the metastable state $`|m1`$, and could therefore activate resist formation on the substrate.

After exposure, the pattern was transferred to the substrate using both dry and wet etches. First, the sample was etched for 35 s in a CHF3 reactive ion etch with parameters 42 mTorr chamber pressure, 100 W forward power, 10 W reflected power, -842 V self-bias, and 25.0 sccm flow rate. Second, the sample was etched for 50 s in 2% HF to prepare an unoxidized silicon face. Finally, the sample was etched for 14 min in 20% KOH at 22 C RIE .

Figure 4 shows a scanning electron micrograph of the etched sample. The features created by ARL are indicated by arrows. Two features are created with a single Raman detuning because the Stark shift is equal to $`\mathrm{\Delta }`$ at a location to either side of its maximum. The dark lines against a bright background indicate places at which silicon was preferentially etched. The sample was patterned across several millimeters of length. At center, for the run parameters above (gradient $`F/\mathrm{}`$=$`2\pi \times 6.4\pm 0.5`$ MHz$`/\mu `$m and $`\tau =29\pm 2`$ ns), the expected feature separation is $`19.0\pm 1.5`$$`\mu `$m, and the spectroscopic limit is $`\delta x_{\mathrm{sp}}=0.86\pm 0.07`$$`\mu `$m. This agrees well with our measurements: a minimum spacing of $`(20.0\pm 0.4)\mu `$m and a rms width of $`(0.9\pm 0.1)\mu `$m.

These fabricated features are thirty times smaller than those that have been detected in holographic masking experiments morinaga , but still two orders of magnitude larger than the smallest features fabricated with atom lithography BellSmall ; CrSmall ; rhese ; RIE . Were the light gradient to be increased, feature sizes in the second ARL scheme demonstrated would be limited by the velocity spread of the atomic beam, at $`\delta x_\mathrm{v}=0.4`$$`\mu `$m. Therefore, the spectroscopically limited feature size observed in Fig. 4 is a factor of two above the minimum observable linewidth for the collimation of our atomic beam. By comparison, the resolution limit of an optimized light mask with a gradient of $`F/\mathrm{}=2\pi \times 6.4`$ MHz/$`\mu `$m is $`\delta x_{\mathrm{}}=54`$ nm for Ar\*.

In conclusion, we have demonstrated the use of frequency encoding in optical masks for atom lithography. Using a two-photon internal state transfer, we create features in silicon with an rms width of $`0.9\pm 0.1`$$`\mu `$m. With a mask based on optical pumping, we show that feature position can be chosen continuously throughout monotonic potential gradients. Such a method could be used to form intricate, sub-micron, two-dimensional patterns with currently available sources. Lithography using internal manipulation schemes could also be extended to alkali atoms by state-selective ionization and deflection. Unionized alkali atoms, like unquenched metastable atoms, would be left to pattern an atom resist (e.g., CsSAMs ).

###### Acknowledgements.

The authors thank N. Dekker and K. S. Johnson for contributions to the early stages of this project. We also thank F. Altarelli, G. Zabow, and S. Coutreau for assistance. This work was supported in part by the NSF Grant No. PHY-9876929, MRSEC Grant No. DMR-9809363, and the Hertz Foundation.


text:
# New Results in the Analysis of the 16O+28Si Elastic Scattering by Modifying the Optical Potential

## I Introduction

The elastic and inelastic scattering between heavy ions have been amongst the main sources of information about complex nuclei over the last 40 years.. There have been numerous experimental investigations for the systems with a combined mass number of $`A_P+A_T60`$, and these investigations have displayed a common unexpected feature near $`\theta _{CM}`$=180 for the elastic and inelastic scattering cross-sections.

The physical origin of the observed structure is not yet fully understood  and presents a challenge to different approaches that have been proposed to explain it. These approaches range from the occurrence of possibly overlapping shape resonances  and the scattering from surface-transparent Optical potentials  to more exotic effects like explicit parity-dependence of the ion-ion potential . The first approach that has been more popular with researchers so far attempts to describe the data by invoking the properties of the average Optical potential . On the other hand, in the second approach, the structure in the excitation function is associated with isolated or nearly isolated partial-wave resonances superimposed on the scattering properties of a standard Optical potential . At present, none of these approaches provides a consistent explanation for all the existing data for this system.

Consequently, the following problems continue to exist for this reaction : $`(1)`$ The explanation of anomalous large angle scattering data; $`(2)`$ the reproduction of the oscillatory structure near the Coulomb barrier; $`(3)`$ the out-of-phase problem between theoretical predictions and experimental data; $`(4)`$ the consistent description of angular distributions together with the excitation functions data.

Within the framework of the coupled-channels method, we have analyzed such light heavy-ion reactions that pose the above-mentioned problems . The new approach that we proposed within the coupled-channels method has successfully explained the experimental data over wide energy ranges for different systems such as 12C+12C, 16O+28Si, 16O+24Mg and 12C+24Mg. One feature we observed in these analyzes was that these reactions were extremely sensitive to the shape of the nuclear potential in the surface region. In this paper, by taking this feature into account, we consider an extensive investigation of the elastic scattering of this system at numerous energies by using a modified potential. Similar to the previously conducted coupled-channels analyzes for such systems , this new technique modifies the shape of the potential at the surface. Thus, we aim to address the above-mentioned problems within the framework of the Optical model and to obtain results that are as good as the results of the coupled-channels method. Accordingly, we analyze the experimental data from 29.0 to 142.5 MeV in the laboratory system over the whole angular range up to 180. The 180 elastic scattering excitation function has also been studied over this energy range. Extensively modified version of the code CHUCK  has been used for the all calculations.

In the next section, we introduce our Optical model and potential parameters to explain the observed experimental data. Then, we show the results of these analyzes in section III from ELab=29.0 MeV to 142.5 MeV. Our conclusion is given in Section IV.

## II The Model

The standard Optical model with folding model potentials or with similar phenomenological potentials such as the square of the Woods-Saxon has failed to describe certain aspects of the experimental data. Therefore, similar to Mackintosh, Kobos and Satchler’s work , in the present calculations, our total real potential consists of the nuclear potential, $`V_{Nuclear}`$, with two small additional potentials, $`\left[U(r)=U_1(r)+U_2(r)\right]`$:

$$V_{total}(r)=\underset{Realpotential}{\underset{}{V_{Nuclear}(r)+U(r)}}+V_{Coulomb}(r)+V_{Centrifugal}(r)$$
(1)
The nuclear potential is assumed to have the square of a Woods-Saxon shape and the parameters are fixed to reproduce the folding model potential of Mackintosh, Kobos and Satchler :

$$V_{Nuclear}(r)=\frac{V_0}{(1+exp(rR)/a)^2}.$$
(2)
where $`V_0`$=761.5 MeV and R=$`r_0`$($`A_a`$1/3+$`A_A`$1/3) with $`r_0`$=0.75 fm and a=1.425 fm. The parameters of the nuclear potential are fixed as a function of energy and kept constant in the present calculations although small changes were observed to improve the quality of the fits.

The Coulomb potential  due to a charge $`Z_ae`$ interacting with a charge $`Z_Ae`$, distributed uniformly over a sphere of radius $`R_c`$, is also added.

$`V_{Coulomb}(r)`$ $`=`$ $`{\displaystyle \frac{1}{4\pi ϵ_{}}}{\displaystyle \frac{Z_aZ_Ae^2}{r}},rR_c`$ (3)
$`=`$ $`{\displaystyle \frac{1}{4\pi ϵ_{}}}{\displaystyle \frac{Z_aZ_Ae^2}{2R_c}}(3{\displaystyle \frac{r^2}{R_c^2}}),r<R_c`$ (4)
where $`R_c`$=7.78 fm is the Coulomb radius, and $`Z_a`$ and $`Z_A`$ denote the charges of the projectile $`a`$ and the target nuclei $`A`$ respectively.

The sum of the nuclear, Coulomb and the centrifugal potentials is shown in figure 1 for various values of the orbital angular momentum. The superposition of the attractive and repulsive potentials results in the formation of a potential pocket. The width and depth of the pocket depend on the orbital angular momentum quantum number for a given nuclear potential. This pocket is very important for the interference of the barrier and internal waves, which creates the oscillatory structure observed in the cross-section.

The significance of the two small additional potentials should be emphasized here, since without them, it would be impossible to fit the experimental data in the elastic scattering calculations. These two small additional potentials are the derivatives of the Woods-Saxon shape and the parameters are shown in table I where

$$U(r)=4U_1a_1df(r,R_1,a_1)/dr+4U_2a_2df(r,R_2,a_2)/dr$$
(5)
$$f(r,R_i,a_i)=\frac{1}{(1+exp((rR_i)/a_i))}$$
(6)
As shown in the insert of figure 2, they create two minima in the nuclear potential between $``$5.5 fm and $``$9.0 fm. The effect of these potentials can be understood in terms of the interference between the internal and barrier waves that correspond to a decomposition  of the scattering amplitude into two components, the inner and external waves. The inner wave comes from the reflection at the inner face of the total real potential pocket and the external wave comes from the reflection at the outer barrier (see figure 1 for the pocket in the total real potential). The presence of the two small potentials affects the phases and magnitudes of these internal and external components. We observe from the parameters in table I that the two small additional potentials are not strong enough to produce pockets in the total real potential although they have a very significant effect on the scattering. This is demonstrated in figure 4 at 41.17 MeV. In that figure, (d) shows the cross-sections obtained when both potentials are omitted. As a result, the calculated oscillations are out-of-phase with the measured ones even at the intermediate angles. Without these potentials, we were unable to refit the data by merely varying the parameters of the real and imaginary potentials. Thus, the oscillatory structure could not be reproduced correctly.

The imaginary part of the potential was taken as the sum of a Woods-Saxon volume and the surface potential :

$$W(r)=W_Vf(r,R_V,a_V)+4W_Sa_Sdf(r,R_S,a_S)/dr$$
(7)
$$f(r,R,a)=\frac{1}{(1+exp((rR)/a))}$$
(8)
with $`W_V`$=59.9 MeV, $`a_V`$=0.127 fm and $`W_S`$=50.0 MeV, $`a_S`$=0.250 fm. These parameters were also fixed in the calculations and only their radii were increased linearly with increasing energy according to the following formulae.

$$R_V=0.06084E_{CM}0.544$$
(9)
$$R_S=0.2406E_{CM}2.191$$
(10)
The imaginary potentials are shown in figure 3 for $`E_{Lab}`$=41.17 MeV.

Moreover, the relative significance of the volume and surface components of the imaginary potential was evaluated for all the energies. For higher energies, omitting the volume term predominantly affected the amplitude of the cross-section at large angles. However, this effect was small and negligible at lower energies. Omitting the surface term increased the cross-sections at large angles which were as much as two orders of magnitude. It was observed that this term had a significant effect at all the considered energies.

## III The Results

Using the above-described Optical model with fixed real and linearly increasing imaginary potential parameters, we have analyzed the experimental data from 29.0 MeV to 142.5 MeV in the laboratory system. The comparisons between experimental data and the Optical model fits are shown in figures 5, 6 and 7. The $`\chi ^2`$ values for some of the energies studied are also shown in table III. As it can be seen from these figures and the $`\chi ^2`$ table, we have obtained excellent agreement with the experimental data over the whole energy range considered.

We have also analyzed the averaged value of the excitation function over the angular range 180$`\pm `$5 with these potentials. The radii of the imaginary potential increased linearly with energy up to 54.0 MeV according to equations (9) and (10). Beyond this energy, we do not have the available experimental data for the large angles. Therefore, as shown in table II, the radii were kept constant at that energy and were not changed for the higher energies. The result is compared with the measured values in figure 8. The overall agreement is quite good up to about 50 MeV where experimental data are available at large angles. In particular, the maxima and minima in the excitation function are obtained correctly. However, it should be emphasized that the calculated excitation function is much more sensitive to the details of the Optical potential than are the fits to the angular distributions. For example, omitting two small potentials completely destroys the agreement with the excitation function.

We were able to obtain almost the same agreement with the experimental data above 50 MeV regardless of these two small potentials. This shows that for higher energies, the barrier/internal wave decomposition is not very important and thus we may remove these two small potentials. By using the same nuclear potential and the same type of imaginary potentials, but by readjusting the depth and diffuseness of its surface component, we were able to fit the experimental data from $`E_{Lab}`$=50 MeV to 142.5 MeV. These results are shown in figure 7 and the relevant parameters are shown in table II.

## IV Summary

We have shown a consistent description of the elastic scattering of the 16O+28Si system from 29.0 MeV to 142.5 MeV in the laboratory system by using the Optical model calculations. In the introduction, we presented the problems that this reaction manifests. We attempted to find a consistent solution to these problems. However, within the standard Optical model calculations, we failed, as others did, to describe certain aspects of the data.

Our previous elastic and inelastic analyzes with light-heavy ion reactions using the coupled-channels formalism had already shown that these reactions are extremely sensitive to the shape of the potential in the surface region. By taking this feature into account, we have used two small additional potentials that modify the shape of the potential in the surface region. As a result, we have obtained excellent agreement with the experimental data over a wide energy range similar to that of the coupled-channels method. The comparison of the results indicates that a global solution to the problems relating to the scattering observables of this reaction over a wide energy range has been provided by this method. Finally, it should be noted that although these additional two small potentials we used are very small and do not create a pocket in the total nuclear potential, they are very effective for the interference of the barrier and internal waves, which creates the oscillatory structure observed in the cross-section. Further work in order to derive this term from a microscopic viewpoint is still under-progress.

## Acknowledgments

Author wishes to thank W.D.M. Rae and N. Ayşe Boztosun for useful comments. He is also grateful to the members of the Nuclear Physics Laboratory in Oxford University.