text:
# Ordering phase transition in the one-dimensional Axelrod model

## I Introduction

In recent years cellular automata and lattice models have been introduced in many fields in order to investigate the emergence of nontrivial collective phenomena out of simple microscopic interactions between a large number of elementary objects Marrobook ; Jensenbook ; Axelrodbook . In the context of social phenomena, numerical simulations of simple automata have been used to study the formation and spreading of cultural domains Axelrod97 ; Axtell96 ; Castellano00 ; Klemm02 . Culture is intended here as the set of behaviors, beliefs, technical standards or values that individuals may mutually exchange. The cultural features are modeled as a set of $`F`$ variables, which can assume integer values between 1 and $`Q`$, specifying the various traits that each feature can assume Axelrod97 . In general, the environment where exchanges take place is a two-dimensional square grid, where each site represents an individual or a group of individuals, whose culture is defined by the specific values of the $`F`$ variables. Starting from an initial random configuration, interactions take place between nearest neighbors, leading to the formation and evolution of regions of homogeneous culture. It is obvious that the assessment of how realistic the model is or to which extent its results can be useful in the interpretation of empirical data in the social science context is a task left to sociologists. From the point of view of statistical physics this model (Axelrod model) is appealing, because it has a non-trivial out of equilibrium dynamics, similar but not equivalent to other well studied models Liggett85 ; Dornic01 .

In previous works Axelrod97 ; Axtell96 ; Castellano00 ; Klemm02 , it has been investigated how the model behavior changes with respect to the values of the parameters $`F`$ and $`Q`$. In particular it has been shown by large-scale simulations that, on a two-dimensional lattice, the model exhibits a non-equilibrium transition between a phase with all the system sharing a common culture and a disordered phase of coexisting regions with different cultures. This is a very interesting example of ordering phase transition in the context of far from equilibrium statistical physics. In this paper we perform a thorough analysis of the model defined in Ref. Axelrod97  on a one-dimensional lattice. This geometry acquires a special role in the description of social systems, since it can be used as starting point for complex topologies such as small-world networks Wattsbook ; Strogatz01 . Noticeably, in dimension $`d=1`$ it is possible to provide an analytical treatment of the model by mapping the evolution dynamics into a multi-species reaction-diffusion system, that we study within a mean-field approximation. We recover also in one dimension the existence of an ordering phase transition, occurring at a critical value of the parameter $`Q`$. When the initial value of the cultural variables on each site is chosen randomly with uniform distribution between 1 and $`Q`$, we find that the system undergoes a transition between a fully ordered state for $`Q<Q_c`$ and a disordered one for $`Q>Q_c`$. When the distribution of initial variables is instead Poissonian of parameter $`q`$, the system exhibits a phase transition between order and disorder, but only for $`F>2`$. For $`F=2`$ any finite value of $`q`$ leads, in the thermodynamic limit, to an asymptotic disordered state. In order to test the analytical predictions we perform numerical simulations of the original dynamics, finding a good agreement with the theoretical results.

The paper is organized as follows. In Section II the model is introduced, the dynamics qualitatively described and the results of the two-dimensional case briefly summarized. In Sec. III we present the mapping to the reaction-diffusion system and the mean-field approach to it. In Sec. IV we report the results of numerical simulations of the model in one dimension, both for the case of uniform and Poisson initial distributions, and compare them with the analytical predictions. Finally in Sec. V we draw some conclusions.

## II The model

We consider a one-dimensional lattice of linear size $`L`$. On each site $`i`$ there are $`F`$ integer variables $`\sigma _{i,f}`$ which define the cultural features of the individuals living on that site. In the original model Axelrod97 , each feature $`f=1,\mathrm{},F`$ on each site $`i`$ is initially drawn randomly from a uniform distribution on the integers between $`1`$ and $`Q`$. The parameter $`Q`$ is a measure of the initial cultural variability (i.e. disorder) in the system. In Ref. Castellano00  instead, the initial values were chosen according to a discrete Poisson distribution of parameter $`q`$ $`\mathrm{Prob}(\sigma _{i,f}=k)=q^ke^q/k!`$, so that the positive real $`q`$ was the average of the values extracted. Results were found to be qualitatively similar independently from the initial distribution, and only the Poissonian case was studied in detail. Here we will consider also the case with uniform distribution and we will discuss analogies and difference between the two cases.

The dynamics is defined as follows. At each time step, a pair of nearest neighbor sites $`i`$ and $`j`$ is randomly chosen. A feature $`f`$ is selected and if $`\sigma _{i,f}\sigma _{j,f}`$ nothing happens. If instead $`\sigma _{i,f}=\sigma _{j,f}`$ then an additional feature $`f^{}`$ is randomly chosen among those taking different values across the bond, $`\sigma _{i,f^{}}\sigma _{j,f^{}}`$. Such a feature is then set equal in the two sites: $`\sigma _{i,f^{}}\sigma _{i,f^{}}^{}=\sigma _{j,f^{}}`$. Time is measured as the total number of attempts divided by the number of sites $`L`$. Axelrod model can be seen as $`F`$ coupled voter models Liggett85 . The evolution of the system is basically a coarsening dynamics of ordered regions. The initial state is totally disordered since the choice of variables in each site is uncorrelated to any other one. As the dynamics starts, the interaction between sites tends to make neighbors similar, introducing correlations in the system. Ordered regions start to form with all sites sharing exactly the same value of all variables. Clearly there are many different possible types of ordered regions. As time proceeds boundaries between regions move, and the average domain size increases. On a finite lattice the ordering process continues until the system reaches an absorbing state such that no further dynamics is possible. Absorbing states can be of two types. A perfectly ordered absorbing state is made by a single region covering the whole system: all sites have the same set of variables $`\sigma _{i,f}`$. Otherwise an absorbing state can also be made by different ordered regions such that adjacent regions have all variables different. In that case all sites are either totally equal or totally different from their neighbors and no interaction may take place.

The nature of the final absorbing state depends on the full evolution of the system, i. e. it is the result of the competition between the disorder of the initial configuration and the drive toward order due to local interactions. Intuitively one expects that for small initial disorder (small $`Q`$ or $`q`$) local interactions will prevail and the final state will be perfectly uniform, while in the opposite limit the ordering process will not last long enough to fully overcome the disorder of the initial condition. This picture turns out to be correct in two dimensions Castellano00 . In that case extensive numerical simulations have shown that for small values of $`q`$, i. e. small initial disorder, the model converges toward an ordered state where a single culture occupies a macroscopically large fraction of the system. On the contrary for large values of $`q`$ many different cultures coexist in the final absorbing state. In the thermodynamic limit of an infinitely large system one sees a well defined non-equilibrium phase transition for a precise value $`q_c`$, separating an ordered phase (for $`q<q_c`$) from a disordered one (for $`q>q_c`$). Interestingly, the nature of the transition is different depending on the number $`F`$ of variables used to define the culture of a site. If $`F=2`$ the transition is continuous, with the order parameter (defined below) vanishing for $`qq_c`$. For $`F>2`$ instead the order parameter exhibits a jump at the transition. The different behavior is reflected also in the distribution of region sizes at the transition. For $`F=2`$ it is a power law with an exponent $`\tau <2`$, while for $`F>2`$ the exponent is universal and larger than 2, so that the transition does not involve the divergence of a correlation length Castellano00 . It is important to stress that the kind of nonequilibrium phase transition discussed here is not an absorbing phase transition Marrobook . Here the system always reaches a frozen absorbing state. The transition occurs in the properties of such absorbing state and hence is static in nature, although it is clearly determined by the dynamical evolution of the problem.

In this paper we shall investigate whether the same kind of transition occurs also in a one-dimensional geometry. To monitor the evolution of the model we mainly focus on the density (i. e. the number divided by the system size $`L`$) of active bonds $`n_A(t)`$ and of fixed bonds $`n_F(t)`$. We define as active a bond between two adjacent sites such that they have more than 0 and less than $`F`$ equal features. We define as fixed a bond between two adjacent sites such that they have no equal variables. Clearly, active bonds are where the dynamics take place and form the boundary between regions that may interact, leading to further ordering. Boundaries formed by fixed bonds are instead frozen. All bonds connecting sites belonging to the same region are neither active nor fixed. We denote their density as $`n_0`$. Clearly, for all times $`n_0+n_A+n_F=1`$. Both the density $`n_A`$ of active bonds and the density $`n_F`$ of fixed bonds measure correlations in the system: $`1/n_F`$ and $`1/n_A`$ are correlation lengths growing in time. Since for any finite lattice an absorbing state is always reached by the dynamics, $`n_A(t)`$ goes asymptotically to zero in all cases. The value of $`n_F(t=\mathrm{})`$ measures the order in the system. If $`n_F(\mathrm{})`$ keeps finite as $`L`$ diverges, the system remains disordered in the thermodynamic limit, while $`n_F(\mathrm{})0`$ indicates full order.

Another useful quantity that will be considered below, is $`S_{max}/L`$, where $`S_{max}`$ is the average over different realizations of the size of the largest region existing in the final frozen state Castellano00 . $`S_{max}/L`$ is the fraction of the whole lattice occupied by the largest region and is an alternative order parameter, being zero for a disordered state and finite when a macroscopically large region exists in the system.

In the next Section we describe an analytical treatment of the problem based on the mapping of the dynamics into a multi-species reaction-diffusion system, which is then analyzed in terms of mean-field equations. This approach turns out to capture remarkably well the qualitative behavior of the system. In some cases it gives predictions that are quantitatively accurate.

## III Theoretical approach

### III.1 Mapping to a reaction-diffusion system

The Axelrod model considers the sites as elementary objects, which interact according to the rules described above. This dynamics can also be thought in terms of bonds between sites whose dynamics is induced by the updating of sites. More precisely, let us define as bond of type $`k`$ a bond connecting two sites with $`k`$ different features. Clearly $`k`$ can be any number between $`0`$ and $`F`$. Bonds of type $`0`$ connect sites with all variables equal, i. e. belonging to the same region. Bonds of type $`F`$ instead connect sites completely different and are therefore the fixed bonds introduced above. Bonds of type $`1,\mathrm{},F1`$ are what we called active bonds in the previous section.

To understand how the update of sites induces a dynamics for bonds, let us consider three neighbor sites (A), (B), (C).

Site (B) is updated depending on the configuration of its neighbor (A). If the bond between sites (A) and (B) is of type $`0`$ or $`F`$ nothing occurs. Otherwise, if the bond is of type $`n`$ one variable on site (B) is set equal to the corresponding variable on site (A). Hence the bond (A)-(B) becomes of type $`n1`$. The update of site (B) modifies in general also the state of bond (B)-(C), which initially was of type $`m`$. Notice that for this bond the number of features which are different may increase, decrease or even remain unchanged, so that after the update the bond may be of type $`m1`$, $`m`$ or $`m+1`$. In summary, if we denote a bond of type $`k`$ as $`A_k`$, the update of a site is then equivalent to a two-particle reaction

$$A_n+A_mA_{n1}+A_{m+j}(j=0,\pm 1).$$
(1)
A specific example is presented in Fig. 1. Bonds of type $`0`$ can be considered as a sort of vacuum state where other types of bonds move. Reactions of the type $`A_1+A_0A_0+A_1`$ are then diffusion processes for bonds of type $`1`$. These are the only diffusion processes allowed.

In order to fully specify the dynamics we must provide the transition rates for the reactions (1). Let us consider the three sites of Fig. 1. One of the features of site (B) is selected and compared with the corresponding one of site (A). In order for a reaction to occur the two variables must be the same. This occurs with a probability $`R_n^{(1)}=(1n/F)(1\delta _{n,0})`$ where we have taken into account that nothing happens if the bond is of type $`0`$.

Hence with probability $`R_n^{(1)}`$, another variable $`\sigma _{B,f^{}}`$ is set equal to $`\sigma _{A,f^{}}`$. What is the effect of this change on the bond (B)-(C)? If $`\sigma _{C,f^{}}`$ was equal to $`\sigma _{B,f^{}}`$, then the update increases by one the number of features different in the bond (B)-(C). The probability that $`\sigma _{C,f^{}}=\sigma _{B,f^{}}`$ is the probability that we randomly select one of the $`Fm`$ equal variables out of the $`F`$ total variables

$$R_{m,+1}^{(2)}=\frac{Fm}{F}=1\frac{m}{F}.$$
(2)
With probability $`1R_{m,1}^{(2)}`$ instead, the variable selected will be different in site (B) and (C), $`\sigma _{C,f^{}}\sigma _{B,f^{}}`$. In this case, since $`\sigma _{B,f^{}}`$ is set equal to $`\sigma _{A,f^{}}`$ the final state of bond (B)-(C) depends on whether $`\sigma _{A,f^{}}`$ and $`\sigma _{C,f^{}}`$ are equal or not. Let us denote as $`\lambda `$ the probability that $`\sigma _{A,f^{}}=\sigma _{C,f^{}}`$. Then the probability that bond (B)-(C) goes from type $`m`$ to type $`m1`$ is

$$R_{m,1}^{(2)}=\frac{m\lambda }{F}.$$
(3)
The probability $`R_{m,0}^{(2)}`$ is obtained via the normalization condition

$$R_{m,0}^{(2)}=\frac{m(1\lambda )}{F}.$$
(4)
In summary the reaction $`A_n+A_mA_{n1}+A_{m+j}`$ occurs with rate

$$R_n^{(1)}R_{m,j}^{(2)}.$$
(5)
The reactions (1) and the rates (5) fully specify the dynamics of the model. Provided one is able to write down the appropriate form of $`\lambda `$ the mapping from the original Axelrod model is exact and Eqs. (1) and (5) are an alternative formulation of the same dynamics.

### III.2 Mean-field approximation

As discussed above, $`\lambda `$ is the probability that two independent random variables (A and C), different from a third one (B), are equal. Here we will determine it by introducing a mean-field approximation, i. e. assuming spatial uniformity.

Let us call $`k`$ the value of variable (B) and $`p_B(k)`$ its distribution. The $`\lambda `$ can be written as $`_kp_B(k)\lambda _k`$ where $`\lambda _k`$ is the probability that A and C assume the same value, different from $`k`$. Clearly

$$\lambda _k=\underset{ik,jk}{}\stackrel{~}{p}_A(i)\stackrel{~}{p}_C(j)\delta _{i,j}$$
(6)
where $`\stackrel{~}{p}_A(i)=p_A(i)/_{ik}p(i)`$ is the probability distribution of the variable $`A`$ constrained to be different from $`k`$. We assume that the distributions of all variables A, B and C remain equal to the initial ones at all times. For uniform initial distribution $`p(i)=1/Q`$, we recover the expected value

$$\lambda =\frac{1}{Q1}.$$
(7)
For variables distributed according to a Poissonian, the plot of $`\lambda `$ as a function of $`q`$, is presented in Fig. 2.

The two behaviors are qualitatively similar. However, one has to remark that in the Poissonian case $`\lambda =1`$ for $`q=0`$, in the trivial situation when all variables are the same. For uniform initial condition instead $`\lambda =1`$ for $`Q=2`$, i. e. when the variables may assume two values.

From the knowledge of the possible reactions and their associated transition rates it is straightforward to write down the mean-field equations for the densities $`n_k`$ of bonds of type $`k`$.

We obtain the set of equations ($`k=0,\mathrm{},F`$)

$`\dot{n}_k`$ $`=`$ $`R_k^{(1)}n_k+R_{k+1}^{(1)}n_{k+1}`$
$`+`$ $`{\displaystyle \underset{N=1}{\overset{F1}{}}}R_N^{(1)}n_N\left[n_k\left(1R_{k,0}^{(2)}\right)+n_{k1}R_{k1,1}^{(2)}+n_{k+1}R_{k+1,1}^{(2)}\right]`$
where clearly $`n_10n_{F+1}`$. More explicitly

$`\dot{n}_k`$ $`=`$ $`\left(1{\displaystyle \frac{k}{F}}\right)(1\delta _{k,0})n_k+\left(1{\displaystyle \frac{k+1}{F}}\right)n_{k+1}`$
$`+`$ $`{\displaystyle \underset{N=1}{\overset{F1}{}}}\left(1{\displaystyle \frac{N}{F}}\right)n_N\left[n_k\left(1{\displaystyle \frac{k}{F}}(1\lambda )\right)+n_{k1}\left(1{\displaystyle \frac{k1}{F}}\right)+n_{k+1}{\displaystyle \frac{k+1}{F}}\lambda \right].`$
It is easily verified that $`_{k=0}^F\dot{n}_k=0`$. Eqs. (III.2) are mean-field equations, since they neglect spatial fluctuations and noise in the system. They are complemented by the initial conditions, which are

$$n_k(t=0)=P_{eq}^{Fk}\left(1P_{eq}\right)^k\left(\begin{array}{c}F\\ k\end{array}\right)$$
(10)
where $`P_{eq}`$ is the probability that a feature takes the same value in two adjacent sites. Clearly $`P_{eq}=1/Q`$ in the uniform case, while

$$P_{eq}=\underset{h=0}{\overset{\mathrm{}}{}}\frac{q^{2h}e^{2q}}{(h!)^2}$$
(11)
for a Poisson initial distribution.

The mean-field system of equations of motion (III.2) consists of $`F`$ coupled nonlinear equations. A detailed study of the temporal evolution is necessarily performed by integrating them numerically. However some important conclusions on the asymptotic state of the system can be drawn analytically by considering the fixed points of Eqs. (III.2). A first crucial observation is that, for any $`F`$, a state with all $`n_k=0`$ ($`k=1,\mathrm{},F1`$) is a fixed point of the equations. Therefore the system always admits a family of frozen asymptotic states, parametrized by the density of fixed bonds $`n_F`$ ($`n_0=1n_F`$ because of normalization). To understand which of these fixed points is selected by the dynamics we perform a stability analysis, by linearizing the dynamical system (III.2) around the generic fixed point $`n_k=0`$.

The stability matrix has $`F`$ eigenvalues. One eigenvalue is always zero, corresponding to the direction parametrized by $`n_F`$ and all made of fixed points. $`F2`$ other eigenvalues are real negative numbers or complex ones with negative real part, corresponding to stable directions. The last one, that we denote as $`\mu `$, is always real and it crucially depends on the value of $`\lambda `$ note . Its value as a function of $`\lambda `$, computed analytically for $`F=2`$ and $`F=3`$, is plotted in Fig. 3 for two different values of $`n_F`$.

Clearly $`\mu `$ depends from the value of $`n_F`$, but for any $`n_F>0`$ it is negative for $`\lambda <1/(F1)`$ and positive for $`\lambda >1/(F1)`$. For $`n_F=0`$ (homogeneous absorbing state) $`\mu =0`$ for any value of $`\lambda `$, so that the stability of this fixed point cannot be assessed directly within a linear analysis. However its stability can be inferred indirectly from the stability of other fixed points: if all other fixed points are unstable the fixed point with $`n_F=0`$ must be stable; if all others are stable, the dynamical system will not go to $`n_F=0`$. This allows us to conclude that for values of the parameter $`Q`$ or $`q`$ of the model such that $`\lambda <1/(F1)`$ the equations will converge exponentially fast to a finite value $`n_F>0`$ implying that the system remains disordered. If instead $`\lambda >1/(F1)`$ the only possible attractive fixed point is $`n_F=0`$, which will attract the solution more slowly than exponentially. In this case the system exhibits a power-law convergence to an ordered state. Whether this transition actually takes place and for which value of the parameter of the initial distribution, it depends on the precise value of $`\lambda `$, that is plotted in Fig. 2.

For uniform initial distribution, since $`\lambda =1/(Q1)`$ it is immediate to see that a phase transition occurs for

$$Q_c=F$$
(12)
so that for $`Q<Q_c`$ the system is ordered, while the system remains asymptotically disordered for $`Q>Q_c`$.

For Poisson initial distribution and $`F>2`$, the scenario is the same, with an ordered phase for small $`q`$ separated by a sharp transition from an ordered phase for large $`q`$. For $`F=2`$ all finite values of $`q`$ imply a positive eigenvalue $`\mu `$ and hence a disordered state. Only for $`q=0`$ order is recovered, but of a very trivial type, since $`q=0`$ means that all variables are the same already in the initial state. Notice that, since $`\lambda `$ decays as $`q^{1/2}`$ for large $`q`$, the transition point diverges for large $`F`$ as $`q_cF^2`$, differently from the uniform case.

For easier insight into the behavior of the model we plot in Fig. 4 the mean-field phase diagram.

### III.3 Uniform initial distribution with $`F=2`$

Let us now give a closer look at the case of uniform initial distribution with $`F=2`$. In this case the set of mean-field equations (III.2) is reduced to two coupled equations

$`\dot{n}_1`$ $`=`$ $`\left[{\displaystyle \frac{3}{4}}+{\displaystyle \frac{1}{4(Q1)}}\right]n_1^2+{\displaystyle \frac{2Q}{2(Q1)}}n_1n_F`$ (13)
$`\dot{n}_F`$ $`=`$ $`{\displaystyle \frac{1}{4}}n_1^2{\displaystyle \frac{1}{2(Q1)}}n_1n_F`$ (14)
where $`n_0`$ has been eliminated using normalization. In the case $`Q=2`$ it is easy to solve (14) explicitly, since the equation for $`n_1`$ does not depend on $`n_F`$. One obtains

$`n_1(t)`$ $`=`$ $`{\displaystyle \frac{1}{2+t}}`$ (15)
$`n_F(t)`$ $`=`$ $`{\displaystyle \frac{1}{(2+t)^{1/2}}}{\displaystyle \frac{1}{2(2+t)}}.`$ (16)
Hence we find, as expected, that both the density of active bonds and of fixed bonds vanish asymptotically as power-laws.

For $`Q>2`$ the two equations (14) become fully coupled. The asymptotic value of $`n_F`$ is found to be

$$n_F(t=\mathrm{})=\frac{Q2}{Q}$$
(17)
for all values of $`Q`$, which will be shown below to agree very well with numerical results.

## IV Numerical results

In order to validate the results of the theoretical approach we compare them with extensive numerical simulations of the model. For each value of the parameters, we average over 10000 or more realizations. For both types of initial distributions, all values of $`F`$ and of $`q`$ or $`Q`$, we find that the system reaches for long times an absorbing state with no active bonds ($`n_A=0`$). Extrapolation to the limit $`L\mathrm{}`$ indicates that also in the thermodynamic limit the system freezes asymptotically. This is the first prediction of the mean-field approach that is verified numerically. We now turn to a more detailed analysis of the asymptotic state as parameters are changed.

### IV.1 Uniform initial distribution

Let us consider first the case with $`F=2`$ variables on each site, which in two dimensions gives a continuous transition. The mean-field approach predicts the existence of a transition between an ordered and a disordered phase for $`Q_c=2`$. Such a prediction is confirmed by the behavior of the order parameter $`S_{max}/L`$. As Fig. 5 illustrates, $`S_{max}/L`$ keeps for $`Q=2`$ a finite value close to 0.71 in the thermodynamic limit, while it goes to zero for $`Q>2`$.

Hence one macroscopic region occupies a finite fraction of the system for $`Q=2`$, while for larger values of $`Q`$, no macroscopic domain is formed.

We remark that in the ordered case the asymptotic value of $`S_{max}/L`$ is not 1. This should not be taken as evidence that the system reaches a partially ordered state with the largest region spanning only a finite part of the whole lattice. In fact the situation is more complex. $`S_{max}`$ is the average over a large number of realizations. In many of them the system becomes fully ordered and $`S_{max}/L=1`$. However in a finite fraction of the realizations the system remains partially disordered so that the order parameter has a value smaller than 1. In practice the distribution of the values of $`S_{max}/L`$ for $`L\mathrm{}`$ is a $`\delta `$-function for $`S_{max}/L=1`$ plus a finite nontrivial contribution for $`S_{max}/L<1`$ This leads to $`S_{max}/L0.71`$.

For $`Q>2`$ the system remains disordered and a finite fraction $`n_F`$ of the bonds is fixed. In the previous Section, Eq. (17) is the mean-field result for $`n_F`$ in the asymptotic state as a function of $`Q`$. Such a formula is compared with numerical results in Fig. 6.

The agreement is excellent, indicating that the mean-field approach is quantitatively correct in this case.

Let us now analyze the whole dynamical evolution. In the main part of Fig. 7 we plot the temporal behavior of the density of active bonds $`n_A`$ in the case $`Q=2`$, for different system sizes.

All curves collapse on the same one when time is rescaled by $`L^2`$, i. e. $`n_A`$ obeys the scaling form

$$n_A(t,L)=t^{1/2}g_A(t/L^2).$$
(18)
The scaling function $`g_A(x)`$ goes to a constant as $`x1`$ and goes to zero exponentially for $`x1`$. In the inset of the same figure it is shown that also the density of fixed bonds obeys a scaling form

$$n_F(t,L)=t^{1/3}g_F(t/L^2)$$
(19)
where $`g_F(x)=\text{const}`$ for $`x1`$ and $`g_F(x)x^{1/3}`$ for $`x1`$.

The physics behind Fig. 7 is clear. In the initial state a finite number of active bonds is present. Notice that for $`F=2`$ and $`Q=2`$ a single type of active bond exists, connecting sites with one equal and one different feature. Interaction between sites implies diffusion of active bonds and annihilation whenever two active bonds meet. Fixed bonds are affected by active bonds but the converse is not true: As far as active bonds are concerned, the system is equivalent to the reaction-diffusion system $`A+A0`$. The decay as $`t^{1/2}`$ of the density of active bonds is easily related to the combination of diffusive wandering and mutual annihilation Redner . The inverse of $`n_A`$ is a correlation length that grows in time. In an infinite system such length grows indefinitely. In a finite lattice, growth stops when $`1/n_A`$ reaches the system size $`L`$, and this yields a cutoff time scaling as $`L^2`$. The evolution of fixed bonds is enslaved by the dynamics active bonds but it is slower. $`1/n_Ft^{1/3}`$ is a correlation length non equivalent to $`1/n_At^{1/2}`$. It grows with time but remains finite when the cutoff time is reached, so that the in the frozen state $`n_F`$ reaches a finite value decreasing as $`L^{2/3}`$, which vanishes in the thermodynamic limit.

A comparison of the exponents found numerically with the mean-field values \[Eq. (16)\] indicates that the mean-field approximation correctly reproduces that $`n_F`$ decays more slowly than $`n_A`$, but it does not predict the correct exponents. This is not a surprising result: we have computed the mean-field rates assuming the system to be homogeneous, while instead ordered regions grow indefinitely. The problem is here the same of the mean-field treatment for the reaction-diffusion system $`A+A0`$ Redner . For higher values of $`Q`$ the disagreement between the mean-field and numerical results is similar. It is somewhat striking that despite the dynamics is not captured quantitatively by the mean-field approximation, the asymptotic value for $`n_F`$ is remarkably accurate (Fig. 6).

Let us now discuss the case $`F>2`$ corresponding to a discontinuous transition in $`d=2`$. For illustrating the behavior common to all values $`F>2`$ let us focus on the case $`F=10`$. As in the two-dimensional case Castellano00 , the system exhibits a discontinuous transition between a fully ordered phase for small $`Q`$ and a disordered phase for $`Q>Q_c`$. This can be inferred from Fig. 8, where $`1S_{max}/L`$ is plotted versus $`Q`$.

By observing that the convergence of $`S_{max}/L`$ to 1 is exponential for $`Q6`$, while curves are bent upward for $`Q>7`$ one is lead to conclude that the transition occurs for $`Q_c7`$.

### IV.2 Poissonian initial distribution

Let us consider now a system such that the initial values of the variables are chosen according to a discrete Poisson distribution of parameter $`q`$. $`q`$ plays here a role similar to the parameter $`Q`$ for the uniform initial distribution, since it fixes the average value of the variables in the initial state. However, no matter how small is $`q`$, arbitrarily large values can be extracted with a Poisson distribution.

Let us consider first the case $`F=2`$ (Fig. 9).

The behavior differs qualitatively from the uniform case. For any value $`q>0`$ it is clear that the order parameter goes asymptotically to zero, even if for small $`q`$, $`S_{max}/L`$ remains close to 1 for relatively large $`L`$. We conclude that there is strong evidence that for $`F=2`$ an initial distribution of Poisson type destroys long range order for any value of $`q>0`$, as predicted by the mean-field approach.

We finally turn to the case $`F=10`$, as an example of what happens in general for $`F>2`$. The value of $`1S_{max}/L`$ is plotted in Fig. 10.

For small $`q`$ it is clear that the system becomes ordered as $`L`$ grows while for large $`q`$, the curves tend to go toward 1, indicating the existence of a discontinuous phase transition for intermediate $`q`$. A rough estimate of the transition point is between $`q=4`$ and $`q=5`$.

These numerical data corroborate the picture emerging from the mean-field approach. For $`F=2`$ we actually find no evidence of long range order in the system, no matter how small $`q`$ is taken. For $`F=10`$ we observe instead the pattern already found with uniform initial distribution: a transition between a disordered phase for large $`q`$ and an ordered one for smaller values of the parameter $`q`$. Also the location of the transition seems to be captured with good accuracy.

## V Conclusions

In this paper we have analyzed in detail the behavior in one dimension of the Axelrod model, a cellular automaton introduced in Ref. Axelrod97  for the description of social influence. Interestingly, the one-dimensional version of the model displays a phenomenology as rich as the two-dimensional case. More precisely, we find by numerical simulations a non-equilibrium phase transition between a disordered and an ordered phase. This phenomenology is recovered accurately by a simple mean-field treatment of the problem, which exploits a mapping onto a reaction-diffusion problem involving $`F`$ different species. The mean-field approach is successful in several aspects. First, it correctly predicts that the system always reaches an asymptotically absorbing state. Secondly, it captures the existence of a transition in the final absorbing state or the absence of it (for Poissonian initial conditions and $`F=2`$). Finally, it provides reasonably accurate estimates of the critical value of the control parameter $`Q`$ or $`q`$. Mean-field fails instead in the prediction of the exponents of the temporal behaviors close to the transition. This is a typical scenario for mean-field treatments: transition points can be computed with good accuracy while exponents are not captured correctly. The theoretical approach presented in this paper is the simplest suitable for this kind of problem. For more precise values of the exponents one should go beyond mean-field, by including diffusion terms and noise in Eq. (III.2) and applying dynamic renormalization group techniques Peliti85 ; Lee95 .


format: text

text:
# Lattice gas models of coherent strained epitaxy

## I Introduction

The phenomena of self-assembly and self-organization of coherent (i. e., dislocation-free) size calibrated nano- and atomic-scale structures observed during the heteroepitaxial growth in some systems Eaglesham and Cerullo (1990); Notzel et al. (1994) are considered to be promising tools for fabrication of microelectronic devices Orlov et al. (1996).

A major factor influencing the phenomenon of self-assembly is the lattice size misfit (LSM) between the substrate and the growing overlayer commonly encountered in heteroepitaxial systems van der Merwe et al. (1994). The LSM-induced strain is believed to be the driving force behind the above phenomena van der Merwe et al. (1994); Priester and Lannoo (1995). So an adequate account of the strain effects should lie at the basis of any theory of strained epitaxy. Because strained systems exhibit complicated kinetics and morphologies, analytic approach is difficult, so a major technique in theoretical studies of strained epitaxy is the kinetic Monte Carlo simulation. The application of this technique, however, is severely hampered by the necessity to simulate the continuous atomic displacements. Because of this, atomistic models in such simulations are currently restricted to rather small systems consisting of only a few thousand atoms Khor and Sarma (2000) while experimentally observed 3D quantum dots sometimes consist of several tens of thousand atoms each Notzel et al. (1994).

Our research is based on the observation that as long as we are interested only in the coherent structures, there is a possibility to map the system onto a purely lattice model because in the absence of dislocations the relaxed atom can only either deviate from its symmetric position inside the same cell or to be displaced to another cell but there always exists a lattice site of a regular lattice to which this atom can be ascribed. So our first goal is to develop a formalism which would allow to map a coherent heteroepitaxial system with the continuous variables onto a lattice gas model with only discrete variables.

One of important goals of the heteroepitaxial studies is the development of techniques of growth of 1D quantum wires (QW) which can be used, e. g., for experimental investigation of the Luttinger model of interacting 1D electrons Segovia et al. (1999). Furthermore, the QW may find application in microelectronics circuitry Orlov et al. (1996), as well as in magnetic memory devices Dorantes-Dávila and Pastor (1998). The latter applications would require the 1D structures of finite length. In the case of the memory devices it is also desirable that these atomic clusters (or chains) were of similar length and that they were arranged into periodic arrays in order to simplify the memory access. All these requirements can be satisfied by self-organized size calibrated structures similar to quantum dots of the 2D epitaxy Notzel et al. (1994). A phenomenological theory explaining the mechanism of formation of quantum dots was proposed in Ref. Priester and Lannoo (1995). The theory is quite general and can be applied to objects in any number of dimensions. So to illustrate the techniques developed in the present paper we will study the conditions of formation of the 1D size calibrated monatomic chains.

## II The model

A simple approach to theoretical description of the misfit is provided by the 1D Frenkel-Kontorova model (FKM) which is frequently been used in qualitative Snyman and van der Merwe (1974); Korutcheva et al. (2000) and semi-quantitative studies Erwin et al. (1999) of strained epitaxy. To illustrate the generality of our approach we first present the formalism appropriate to the 2D epitaxy but in concrete examples we will restrict ourselves to the 1D case.

We consider an ensemble of a fixed number ($`N`$) of atoms coherently deposited on a surface with a square lattice of deposition sites. The rectangular lattice geometry was chosen because it allows for the separation of $`x`$ and $`y`$ variables (see Ref. Snyman and van der Merwe (1974) and below). Let us first consider more general model with the energy of the epilayer of the form

$$U=\underset{𝐢}{}n_𝐢V_s(𝐑_𝐢+𝐮_𝐢)+\frac{1}{2}\underset{\mathrm{𝐢𝐣}}{}n_𝐢n_𝐣V(𝐮_𝐢+𝐑_𝐢𝐮_𝐣𝐑_𝐣),$$

where $`𝐑_𝐢=(a^xi_x,a^yi_y)`$ ($`a^\gamma `$ the lattice constants in the two directions, $`i_x`$ and $`i_y`$ integers ), $`n_𝐢=0,1`$ is the occupation number of site $`𝐢`$, $`𝐮_𝐢`$ the atomic displacement, $`V_s`$ the potential of interaction with the substrate, and $`V`$ the interatomic potential.

The analysis of the system considerably simplifies at low temperatures where it can be approximately reduced to a lattice gas model. This is achieved by exploiting the fact that the residence time of atoms at the lattice sites can be arbitrarily large due to the Arrhenius law obeyed by the probability of activated hopping over the potential barriers separating neighboring sites Uebing and Homer (1991). The dynamics of the variables $`𝐮_𝐢`$, on the other hand, do not have any energy barriers. So at sufficiently low temperature these variables are capable of reaching their thermal equilibrium distribution during the time intervals between the atomic hops, i. e., with the atomic configuration remaining unchanged. Averaging over $`𝐮_𝐢`$ will leave us with an effective non-equilibrium free energy function $`F_{\text{eff}}`$ of variables $`n_𝐢`$ only:

$$\mathrm{exp}(F_{\text{eff}}/k_BT)=\underset{\{n_𝐢=1\}}{}d𝐮_𝐢\mathrm{exp}(U/k_BT).$$
(1)
This purely lattice model can be further used in both equilibrium and kinetic studies.

For the purposes of qualitative analysis it will suffice to average out the displacement variables in Eq. (1) in the harmonic approximation, i. e., by expanding $`V_s`$ and $`V_p`$ in the above equation up to the second order in the displacement variables $`𝐮_𝐢`$ Erwin et al. (1999); Orr et al. (1992) which is valid for small $`|u_i^\gamma |/a^\gamma `$. This allows to perform the calculation exactly (see below). However, the integration can be extended to approximately account also for anharmonic terms. This extension, besides making the approach more accurate, can also account for some qualitative phenomena Korutcheva et al. (2000).

To facilitate comparison with other studies based on the Frenkel-Kontorova model Snyman and van der Merwe (1974); Ratsch and Zangwill (1993); Korutcheva et al. (2000); Erwin et al. (1999) we write the harmonic approximation as the second order power series expansion in $`𝐮_𝐢`$ for $`V_s`$ and in $`𝐮_𝐢𝐮_{𝐢\gamma }`$ for the pair potential $`V_p`$:

$`U`$ $``$ $`V_sN+{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐢\gamma }{}}n_𝐢k_s^\gamma (u_𝐢^\gamma )^2+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mathrm{𝐢𝐣}}{}}n_𝐢n_𝐣V_{\mathrm{𝐢𝐣}}`$ (2)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐢\gamma }{}}k_p^\gamma [(u_𝐢^\gamma u_{𝐢+\gamma }^\gamma f_\gamma )^2f_\gamma ^2]n_𝐢n_{𝐢+\gamma }`$
$`=`$ $`V_sN+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mathrm{𝐢𝐣}}{}}V_{\mathrm{𝐢𝐣}}n_𝐢n_𝐣`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mathrm{𝐢𝐣}\gamma }{}}k_p^\gamma D_{\mathrm{𝐢𝐣}}^\gamma u_𝐢^\gamma u_𝐣^\gamma +{\displaystyle \underset{𝐢\gamma }{}}k_p^\gamma f_\gamma u_𝐢^\gamma _𝐢^\gamma n_𝐢n_{𝐢+\gamma }`$
where $`V_s`$ and $`V_{\mathrm{𝐢𝐣}}`$ are the values of the substrate and the pair potentials at the lattice sites for the atom in the symmetric position (the zero order approximation), $`\gamma `$ denotes $`x`$ or $`y`$ component; $`f_\gamma =V_\gamma ^{}/V_{\gamma \gamma }^{\prime \prime }`$ (where $`V`$ is the pair potential) is interpreted as the misfit parameter in the direction $`\gamma `$, $`k_s^\gamma `$ and $`k_p^\gamma `$ are the second derivatives of the corresponding potentials, $`D_{\mathrm{𝐢𝐣}}^\gamma `$ is the dimensionless dynamical matrix defined by this equality, and $`_𝐢^\gamma F_𝐢F_𝐢F_{𝐢\gamma }`$. A major simplification is achieved under the nearest neighbor (NN) approximation for the relaxation because, as is seen from the last line of Eq. (2), in this case the variables $`u_𝐢^x`$ and $`u_𝐢^y`$ separate, and the relaxations along the two directions are independent.

With approximation (2) the statistical average amounts to the Gaussian integration to give

$`F_{\text{eff}}=`$ $`\overline{V}_sN+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mathrm{𝐢𝐣}}{}}n_𝐢n_𝐣V_{\mathrm{𝐢𝐣}}{\displaystyle \frac{k_BT}{2}}{\displaystyle \underset{\gamma }{}}\mathrm{ln}detG^\gamma `$ (3)
$`{\displaystyle \frac{k_pf^2}{2}}{\displaystyle \underset{\mathrm{𝐢𝐣}\gamma }{}}G_{\mathrm{𝐢𝐣}}^\gamma _𝐢^\gamma n_𝐢n_{𝐢+\gamma }_𝐣^\gamma n_𝐣n_{𝐣+\gamma },`$
where $`G^\gamma =1/D^\gamma `$ and in $`\overline{V}_s`$ we gathered all terms which are proportional to the total particle number, such as the normalization of the determinant coming from the in-plane Gaussian integration and a similar term from the integration along $`z`$-direction. We do not discuss these contributions here because in the present study we restrict ourselves to systems with a fixed number of deposited atoms $`N`$. In case of necessity these terms can be easily recovered. We also note the entropic contribution (the second term in the first line) which naturally appears in our formalism and which was shown to be crucial for proper description of the processes of deposition Venables (1987) as well as for the correct prediction of the shape of the atomic clusters Doye and Calvo (2000).

Because of the gradient factors, only the ends of contiguous chains of atoms contribute into the last term of Eq. (3) leaving two matrix elements of $`G_{\mathrm{𝐢𝐣}}^\gamma `$ in the sum over $`\mathrm{𝐢𝐣}`$ inside every chain: the diagonal one which we denote as $`G_l^{(0)}`$ ($`l`$ the length of the chain; we omit the superscript $`\gamma `$ to simplify notation) and the matrix element $`G_l^{(l1)}`$ connecting the two ends of the chain. Furthermore, in the NN approximation the matrix $`D`$ is block-diagonal because the atoms belonging to different chains do not couple. Therefore, the determinant factorizes and the relaxation part of the free energy which consists of the terms in Eq. (3) containing $`G`$ takes the form

$`F_{\text{relax}}`$ $`=`$ $`{\displaystyle \underset{\text{chains}}{}}\{{\displaystyle \frac{k_BT}{2}}\mathrm{ln}det{\displaystyle \frac{G_l}{G_1}}+k_pf^2[G_l^{(0)}G_l^{(l1)}]\}`$ (4)
$``$ $`{\displaystyle \underset{\text{chains}}{}}(TS_l+W_l){\displaystyle \underset{\text{chains}}{}}E_l,`$
where the summation is over the chains consisting of two or more atoms. Besides, in the first term on the right hand side (r. h. s.) we subtracted the single atom entropy term by assuming that it is accounted for in $`\overline{V}_s`$ which we discarded to simplify notation. It is easy to show that $`F_{\text{relax}}`$ can be expanded into an infinite sum of multiatom interactions of the form $`V_ln_𝐢n_{𝐢+\widehat{\gamma }}\mathrm{}n_{𝐢+(l1)\widehat{\gamma }}`$, where $`\widehat{\gamma }`$ is the unit vector in the direction $`\widehat{\gamma }`$. The expansion coefficients are given by the expression

$$V_l=E_l2E_{l1}+E_{l2}$$
(5)
valid for all $`l2`$ with $`E_0=E_1=0`$.

## III The relaxation free energy

Omitting the superscript $`\gamma `$ to simplify notation, the matrix $`D_l`$ for chain of length $`l`$ is defined according to Eq. (2) as

$$k_p\underset{ij}{}(D_l)_{ij}u_iu_j=k_su_1^2+k_p(u_1u_2)^2+k_su_2^2+\mathrm{}+k_p(u_{l1}u_l)^2+k_su_l^2.$$
(6)
Hence, the square matrix $`k_pD_l`$ have the following structure:

$$k_pD_l=\left(\underset{l}{\underset{}{\begin{array}{ccccccc}k_p+k_s& k_p& 0& \mathrm{}& 0& 0& 0\\ k_p& 2k_p+k_s& k_p& 0& \mathrm{}& 0& 0\\ 0& k_p& 2k_p+k_s& k_p& 0& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& \mathrm{}& 0& k_p& 2k_p+k_s& k_p& 0\\ 0& 0& \mathrm{}& 0& k_p& 2k_p+k_s& k_p\\ 0& 0& 0& \mathrm{}& 0& k_p& k_p+k_s\end{array}}}\right)$$
(7)
In Appendix A we have shown that the matrix elements of $`G_l=1/D_l`$ entering the expression for the relaxation free energy (4) can be calculated with the use of recursion formulas for the tridiagonal matrices as \[see Eqs. (A) and (20)\]

$$\{\begin{array}{ccc}G_{l+1}^{(0)}\hfill & =& 1/(1+\alpha d_l)\hfill \\ G_{l+1}^{(l)}\hfill & =& detG_{l+1}=G_{l+1}^{(0)}b_l,\hfill \end{array}$$

where $`\alpha =k_s/k_p`$; $`d_l`$ and $`b_l`$ are generated by the recursion relations

$$d_{l+1}=1/(2+\alpha d_l)\text{and}b_{l+1}=d_{l+1}b_l.$$
(8)
The recursion relations (8) at large $`l`$ drive $`d_l`$ and $`b_l`$ to the fixed point $`d_{},b_{}`$:

$`d_{}`$ $`=`$ $`1+{\displaystyle \frac{\alpha }{2}}\sqrt{\alpha +{\displaystyle \frac{\alpha ^2}{4}}}<1`$ (9)
$`b_{}`$ $`=`$ $`0.`$ (10)
This means that at large $`l`$ $`G_l^{(0)}`$ saturates to a constant value while $`G_l^{(l)}=detG_ld_{}^l`$ \[see Eqs. (III)\]. This in particular means that the entropic contribution into $`F_{relax}`$ which is proportional to $`\mathrm{ln}(detG_l)`$ is asymptotically linear in $`l`$ (see Fig. 1).

Two contributions into $`F_{\text{relax}}`$ calculated with these formulas are shown in Fig. 1. The entropic contribution is negative because in our canonical ensemble formalism we discarded the entropy corresponding to $`N`$ isolated atoms which is reflected in the denominator $`G_1`$ in the expression (4) for the free energy. This means, in particular, that when $`k_p0`$ the atomic displacements within the cell became mutually independent (the on-site pair interactions do not depend on the atomic positions inside the cell), so $`G_lG_1\widehat{I}`$, where $`\widehat{I}`$ is the unit matrix, and the entropy term tends to zero \[see Eq. (4) and Fig. 1\]. When $`k_p`$ is not equal to zero, the atomic displacements became restricted by its neighbors, so their entropy diminishes. For large chains the interior atoms are practically in translationally invariant environment, so for long chains the entropy loss is proportional to $`l`$. The above restrictions on the atomic displacements are the more stringent the larger $`k_p`$. This is reflected in the slope of $`S_l`$ which is larger for smaller $`\alpha =k_s/k_p`$ (see Fig. 1).

Because the entropic contribution in Eq. (4) is multiplied by $`T`$, at low temperature it can be neglected. The remaining term $`W_l`$ is negative and saturates to a finite value at large $`l`$ (see Fig. 1b). From Eq. (2) it follows that the NN interatomic interaction $`V_{\text{NN}}=V_{\text{NN}}^{\text{ch}}+k_pf^2`$, where the first term ($`V_{\text{NN}}^{\text{ch}}`$) is some “chemical” interaction and the second term is the unrelaxed repulsion due to LSM. Thus, the relaxation energy in Eq. (4) is the difference between the relaxed elastic energy of the chain $`E_{rel}^{el}(l)`$ and the unrelaxed one $`E_0^{el}(l)=(l1)k_pf^2`$ with both energies being positive, as it should be for elastic energies. Thus,

$$W_l=E_{rel}^{el}(l)E_0^{el}(l)$$

and is negative because relaxation lowers the energy. Now, if $`V_{\text{NN}}^{\text{ch}}`$ is negative but such that $`V_{\text{NN}}`$ is not too large, then the reduced chain energy $`[W_l+V_{\text{NN}}(l1)]/l`$ will have a minimum corresponding to an equilibrium chain length. In two dimensions these chains will unify in square platelets. Thus, our model contains the mechanism for the size calibration of atomic clusters which was phenomenologically described in Ref. Priester and Lannoo (1995). In 2D this conclusion as well as the formula of Ref. Priester and Lannoo (1995) were verified with the use of the Monte Carlo simulation Tokar and Dreyssé (2002).

Because the entropic contribution $`S_l`$ (see Fig. 1a) is practically linear in $`l`$ according to Eq. (5) which has the form of the discrete second derivative it essentially contributes only into the pair interaction $`V_2`$ (the multiatom contributions were found to be $`10\%`$). Thus, the relaxation entropy formally amounts to effective NN interatomic repulsion which grows linearly with temperature. The entropic forces of this kind were earlier discovered in alloys (see Appendix E of Ref. van de Walle and Ceder (2002) and references to earlier literature therein).

At low temperature the main contributions into $`V_l`$ for $`l3`$ come from $`W_l`$ (Fig.1b). Because the r. h. s. of Eq. (5) has the form of the discrete second derivative and the curve $`W_l(l)`$ is concave, all multiatom interactions are repulsive. On the other hand, in the case under consideration when atoms are supposed to assemble into clusters, there should exist attractive interactions which in our model are necessarily the pair ones \[see Eq. (2)\].

## IV Examples from metallic heteroepitaxial systems

To illustrate the above formalism with realistic examples of strained epitaxy we consider two metallic heteroepitaxial systems—Ag/Pt and Pt/Co—which currently are being actively studied both experimentally Becker et al. (1993); Gambardella et al. (2000a, b) and theoretically Goyhenex and Tréglia (2000); Goyhenex et al. (1999). For simplicity we considered 1D case and used the geometry of Ref. Picaud et al. (2000) where the growth on the steps of the closed packed (111) vicinal surface was studied. The position of a deposited atom was relaxed to its equilibrium value in order to find the value $`k_s`$ as the second derivative of the potential near equilibrium. The many-body “potentials” and corresponding parameters were taken from Ref. Picaud et al. (2000) for the Ag/Pt system and from Ref. Goyhenex et al. (1999) for the Pt/Co system. These potentials were devised specifically to application in the heteroepitaxy and for the (111) surfaces so we expect our results below are quite reliable. Other parameters listed in Table 1 were calculated for the atomic pairs relaxed only in vertical position because in our approach we need the first and second derivatives \[see discussion after Eq. (2)\] calculated in the center of of symmetry of the 1D cell. In both calculations only the NN interactions were taken into account because the NNN ones were found to be negligeable.

According to our calculations the Ag/Pt system has the following parameters: (the energy unit is eV, the length unit Å): $`V_{\text{NN}}^p0.57`$, $`V_{\text{NNN}}^p8.810^3`$, $`k_pf^20.72`$, and $`\alpha 3.7`$. The large value of $`\alpha `$ means that the relaxation of the strain is very weak (see Fig. 1) so there is no size calibration with the above parameters. The crucial parameter $`\alpha `$, however, can strongly vary in different systems. For example, according to our estimates based on the potentials of Ref. Goyhenex et al. (1999), in Pt/Co $`\alpha 0.32`$ is an order of magnitude smaller. Because of this the system is quite close to the self-assembly but the misfit strain is still too small. To enhance the misfit, in our model calculations we assumed that the Co underlayer is further compressed (e. g., by means of deposition on an appropriate substrate) by the factor $`1ϵ`$. (see Table 1 and Fig. 2). But for $`ϵ`$ below the critical value $`ϵ_c2\%`$ there is no self-assembly at $`T=0`$ because there is no minimum in $`E_l(T=0)/l`$ (see Ref. Priester and Lannoo (1995)). However, because of the $`T`$-dependent entropic contribution in the effective energy (3), slightly below $`ϵ_c`$ the system does exhibit the self-assembly at intermediate temperatures which, however, disappears as $`T0`$ (see Fig. 2). Thus, our model predicts a new phenomenon which may be called the entropy driven transient self-assembly. For higher values of strain the $`E_l(T=0)/l`$ curve does have a minimum (see Fig. 2) in which case the qualitative analysis of Ref. Priester and Lannoo (1995) fully applies.

## V conclusion

In this paper we considered a simple model of strained 1D epitaxy and have shown that similarly to the 2D case the mechanism of self-assembly and size calibration proposed in Ref. Priester and Lannoo (1995) is operative also in this case. We considered only two explicit examples for which there exist in literature reliable interatomic potentials and in one system already found the size calibration for stressed substrate. In our opinion, this shows that the above phenomena should be as common in 1D as they are in 2D heteroepitaxy. Further argument in favor of this conclusion is that for the size calibration to be plausible the crucial parameter $`\alpha =k_s/k_p`$ should be as small as possible. This favors small values of $`k_s`$. However, the geometry considered in the present paper is not quite favorable because of the high coordination of the atoms deposited at the steps with 5 NN atoms of the substrate (see Ref. Picaud et al. (2000)) which enhances $`k_s`$. It may be hoped that with lower coordination of the deposited atoms (as in the case of 1D structures of Ref. Segovia et al. (1999)) the conditions for the size calibration will be more favorable.

###### Acknowledgements.

One of us (VT) expresses his gratitude to University Louis Pasteur de Strasbourg and IPCMS for their hospitality.

## Appendix A

According to Eq. (7)

$$D_l=\left(\underset{l}{\underset{}{\begin{array}{ccccccc}1+\alpha & 1& 0& \mathrm{}& 0& 0& 0\\ 1& 2+\alpha & 1& 0& \mathrm{}& 0& 0\\ 0& 1& 2+\alpha & 1& 0& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& \mathrm{}& 0& 1& 2+\alpha & 1& 0\\ 0& 0& \mathrm{}& 0& 1& 2+\alpha & 1\\ 0& 0& 0& \mathrm{}& 0& 1& 1+\alpha \end{array}}}\right),$$
(11)
where $`\alpha =k_s/k_p`$. Because matrices $`D_l`$ are tri-diagonal, their determinants satisfy recurrence relations which can be used to calculate all quantities entering Eqs. (4) and (5). The diagonal element of the matrix $`G_l=1/D_l`$ is

$$G_l^{(0)}=r_{l1}/detD_l,$$
(12)
where $`r_{l1}`$ is the determinant of the matrix obtained from $`D_l`$ by crossing out its first row and the first column:

$$r_{l1}=\left|\underset{l1}{\underset{}{\begin{array}{cccccc}2+\alpha & 1& 0& \mathrm{}& 0& 0\\ 1& 2+\alpha & 1& 0& \mathrm{}& 0\\ 0& 1& 2+\alpha & 1& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& \mathrm{}& 0& 1& 2+\alpha & 1\\ 0& 0& \mathrm{}& 0& 1& 1+\alpha \end{array}}}\right|.$$
(13)
Expanding $`detD_l`$ with respect to the first row we get

$$detD_l=(1+\alpha )r_{l1}r_{l2}.$$
(14)
Comparing this with Eq. (12) we get

$$G_l^{(0)}=1/(1+\alpha d_{l1}),$$
(15)
where $`d_{l1}=r_{l2}/r_{l1}`$. Now, expanding determinant Eq. (13) with respect to the elements of the first row we get the following three term recurrence relation

$$r_{l1}=(2+\alpha )r_{l2}r_{l3}.$$
(16)
But because Eq. (12) includes only the ratio $`r_{l2}/r_{l1}`$ by dividing Eq. (16) by $`r_{l2}`$ we can transform it to a simpler form

$$\frac{1}{d_{l1}}=2+\alpha d_{l2}.$$
(17)
The off-diagonal matrix element $`G_l^{(l1)}`$ of the inverse matrix $`G_l`$ is equal to the ratio of the determinant of the matrix obtained from $`D_l`$ by crossing out its first row and the last column multiplied by $`(1)^{l+1}`$ and divided by $`detD_l`$. As is easy to see from Eq. (11) the matrix thus obtained is a triangular matrix with the diagonal elements being all equal to -1. Hence its determinant is equal to $`(1)^{l1}`$ and

$$G_l^{(l1)}=1/detD_l=detG_l.$$
(18)
From Eq. (14) we get

$$G_l^{(l1)}=\frac{1}{detD_l}=\frac{1}{r_{l1}(1+\alpha d_{l1})}=\frac{b_{l1}}{1+\alpha d_{l1}}=b_{l1}G_l^{(0)},$$
(19)
where $`b_{l1}=1/r_{l1}=d_{l1}b_{l2}`$ (see the definition of $`d_{l1}`$ above). Finally, making the shift $`l1l+1`$ the above formulas can be summarized as

$$\{\begin{array}{ccc}G_{l+1}^{(0)}\hfill & =& 1/(1+\alpha d_l)\hfill \\ G_{l+1}^{(l)}\hfill & =& detG_{l+1}=G_{l+1}^{(0)}b_l,\hfill \end{array}$$

where $`d_l`$ and $`b_l`$ are generated by the recursion relations

$$d_{l+1}=1/(2+\alpha d_l)\text{and}b_{l+1}=d_{l+1}b_l$$
(20)
initialized by $`d_0=b_0=1`$.


text:
# Ellipsoidal variability of symbiotic giants

## 1. Introduction

Amongst the evidence for the predominance of wind-accretion is the fact that ellipsoidal variations, characteristic of tidally distorted stars are rarely observed for symbiotic binaries. However, the general absence of the tidally distorted giants in the symbiotic binaries can be in fact due to the lack of systematic searches for the ellipsoidal variations in the red and near-IR range where the cool giant dominate the continuum light. This can be best illustrated by the eclipsing symbiotic systems BF Cyg (Fig. 1) and CI Cyg, (Fig. 2; Mikołajewska 2001). In both systems, during the hot component outburst and its decline, eclipses in the $`UBV`$ bands are narrow with well-defined eclipse contacts whereas at quiescence very broad minima and continuous nearly sinusoidal variation are observed. In addition, the quiescent $`VRI`$ and near-IR light curves show a modulation with half-orbital period as expected for ellipsoidal variability of the red giant.

## 2. Analysis of the quiescent light curves of YY Her

Similar behaviour has been recently found in YY Her. In particular, the photometric and spectroscopic study of YY Her during the major outburst in 1993 and its decline (Tatarnikova et al. 2001) has shown that the minima are due to obscuration of the hot component and the ionized nebula. The return to quiescence in 1997/98 was accompanied by significant changes in the shape of light curves, and in particular by the appearance of a secondary minimum in the $`VRI`$ light curves (Fig. 2). Tatarnikova et al. (2001) interpreted the secondary minimum in terms of ellipsoidal changes of the red giant, whereas Hric et al. (2001, hereafter H01) argued that the secondary minimum is caused by an eclipse of the red giant by a circumstellar envelope around the hot component.

The salient characteristic of the quiescent light curves, namely: (i) nearly sinusoidal shape of the $`U`$ light curves; (ii) a flat toped $`B`$ light curve; and (iii) the two minima of almost comparable depth in $`I^{}`$ light, combined with the fact that the quiescent, ultraviolet + optical, spectra of YY Her can be satisfactorily matched with a three-component model composed of a cool giant dominating the red spectral range, a gaseous nebula dominating the near UV range, and the hot component predominant in the short UV flux, led us to a simple phenomenological model in which the light changes are described by a combination of variable nebular emission and ellipsoidal changes of the red giant (Mikołajewska et al. 2002a, hereafter M02).

Each light curve has been converted to fluxes, $`F_\lambda =10^{0.4m_\lambda }`$, and then independently fitted by the Fourier cosine series

$$F_\lambda =F_{\lambda ,0}A_{\mathrm{n},\lambda }\mathrm{cos}\varphi A_{\mathrm{g},\lambda }\mathrm{cos}2\varphi .$$
(1)
In the case of $`U`$ data, the best fit is obtained with the $`\mathrm{cos}\varphi `$ term only. This indicates that the contribution of the giant to the $`U`$ light is negligible, and that all changes can be attributed to the nebular emission. Assuming that the amplitude of the nebular emission changes is the same at all wavelengths, $`(F_{\mathrm{n},\lambda }^0+A_{\mathrm{n},\lambda })/(F_{\mathrm{n},\lambda }^0A_{\mathrm{n},\lambda })(F_{\mathrm{U},0}+A_\mathrm{U})/(F_{\mathrm{U},0}A_\mathrm{U})=C4.5`$, we can separate the contributions of the nebular, $`F_{\mathrm{n},\lambda }^0`$, and the giant, $`F_{\mathrm{g},\lambda }^0`$, components to the constant $`F_{\lambda ,0}`$. Our best fit models are plotted in Fig. 2.

Using the maximum, $`F_{\lambda ,\mathrm{max}}(\varphi =0.5)=F_{\mathrm{n},\lambda }^0+A_{\mathrm{n},\lambda }`$, fluxes of the nebular component, corrected for the contribution from emission lines and the interstellar reddening the maximum emission measure of the nebula, $`EM`$, has been derived for each band (see M02 for details). These $`EM`$ values, agree with each other to within a factor of 2, and they cluster around $`EM2\times 10^{14}`$ and $`3\times 10^{14}\mathrm{cm}^5`$ for $`T_\mathrm{e}=10^4`$ and $`2\times 10^4\mathrm{K}`$, respectively. The nebular emission measure derived from our light curve analysis requires a Lyman continuum luminosity of $`6`$$`5\times 10^{45}(d/1\mathrm{kpc})^2\mathrm{phot}\mathrm{s}^1`$ for $`T_\mathrm{e}=10^4`$$`2\times 10^4\mathrm{K}`$, which agrees with independent estimates for the quiescent $`L_\mathrm{h}`$ to within a factor $`<2`$ (M02).

The magnitudes, and colours of the cool component derived from our light curve analysis, $`BV=1.73`$–1.77, $`VR=1.73`$–1.92 and $`VI=3.33`$–3.61, are in excellent agreement with those expected for a moderately reddened M3–M4 giant. Moreover, they indicate somewhat higher temperature (earlier spectral type) at maxima ($`\varphi =0.25`$ and 0.75) than that at minima ($`\varphi =0`$ and 0.5) as expected for a tidally distorted giant.

## 3. Concluding remarks

One of fundamental questions in relation of symbiotic binaries is the mechanism that powers the Z And-type multiple outburst activity observed in many classical symbiotic systems. The quiescent spectra of these systems can mostly be fitted by a hot stellar source (most likely a white dwarf powered by thermonuclear burning of the material accreted from its companion) and its ionizing effect on a nebula. Their outburst activity, however, with time scales of a few/several years cannot be simply accounted by the thermonuclear models. A promising interpretation of this activity involves changes in mass transfer and/or accretion disc instabilities. Detection of an ellipsoidal hot continuum source during outbursts of CI Cyg, AX Per, YY Her, AS 338 and other Z And-type systems suggests the presence of an optically thick accretion disc which strongly supports this interpretation (Mikołajewska et al. 2002b, and references therein) Related to this problem is the process of mass transfer – Roche lobe overflow or stellar wind – and possibility of an accretion disc formation. The presence of secondary minima in the light curves of YY Her, apparently due to ellipsoidal changes of the red giant, as well as in a few other similar systems provides here important clues. Although it is very premature to claim that all symbiotic systems with the Z And-type activity do have tidally distorted giants, and – at least during active phase – accretion discs, whereas the non-eruptive systems (such as RW Hya) do not, the former is certainly the case for YY Her, CI Cyg, BF Cyg and a few other active systems.

Systematic searches for the ellipsoidal variations in both active and non-eruptive symbiotic stars are necessary to address the problem, and confirm or exclude tidally distorted giants. The observations must be, however, carried in the red and near-IR range where the cool giant dominated the continuum light. We also note that any optical/red light curve analysis must be supported by spectroscopic information about contamination of the broad-band photometry by the nebular component.

### Acknowledgments.

AAT and SYuS would like to thank LOC and SOC for hospitality and financial support. This research was partly founded by KBN Research Grant No. 5 P03D 019 20, and RFBR grants Nos. 00-15-96553, 02-02-16462, 02-02-26635 and 02-02-16235.

## References

Hric, L., Petrik, K., Niarchos, P. et al. 2001, IBVS 5046 (H01)

Mikołajewska, J. 2001, in ASP Conf. Ser., 246, Small-Telescope Astronomy on Global Scales, ed. B. Paczyński, W.P. Chen, & C. Lemme, 167

Mikołajewska, J., Kolotilov, E.A., Shugarov, S.Yu., Yudin, B.F., 2002a, A&A, in press (M02)

Mikołajewska, J., Kolotilov, E.A., Shenavrin, V.I., Yudin, B.F. 2002b, in ASP Conf. Ser., 261, The Physics of Cataclysmic Variables and Related Objects, ed. B.T. Gansicke, K. Bauermann & K. Reinsch, 645

Tatarnikova, A.A., Esipov, V.F., Kholotilov, J.A. et al. 2001, Astron. Letters, 27, 703


text:
# VLBI Techniques in Pulsar Astronomy

## 1. Introduction

The case for measuring pulsar parallaxes is strong – many models based on other measurements of pulsars require estimates of the distances to the pulsars in the sample. Applications of pulsar distances are not within the scope of this paper as they are discussed in detail by Chatterjee in the following paper. This paper concentrates on making useful measurements of pulsars with VLBI and focuses mainly on high-precision astrometry.

## 2. VLBI

Interferometry is the radio astronomers answer to the optical astronomer’s resolution. Building large single centimeter-wave telescopes with resolutions that come close to those of large optical telescopes is impossible. At 1.5 GHz, the dish would need to be 20 km in diameter – quite a mechanical feat. Even if such a marvel could be constructed, the ionosphere and troposphere could not be expected to be coherent over the entirety of the aperture. Fortunately phase-preserving amplifiers and stable electronics are available in the gigaHertz frequencies which allow the electrical combination of signals from many small antennas. Also important is that the radio photons are in the classical regime with occupation numbers much greater than unity, which allows amplification and duplication of the signals without a large sensitivity penalty. Once one can combine the signals from many small antennas, one can use aperture synthesis to probe features smaller than a milliarcsecond.

VLBI is probably best distinguished from connected-element interferometry by its scale – thousands of kilometers, typically. The distance and geographic barriers require the sampled electric field measurements to be recorded onto (usually magnetic) media to be shipped from the antennas to a central correlator and the use different clocks at each station which must be stable and synchronized to better than microsecond precision. It also requires use of very detailed correlator models, especially for astrometric projects. These models include effects such as bending of light around the sun and Jupiter, tidal deformation of the earth, excursions from Earth’s average spin, and plate tectonic motion. This model is used to calculate, to 10 ps or better accuracy, the expected geometric delay between the arrivals of the radiation from an astronomical source at the antennas. The input parameters to the model include station positions, the time of the observation, and coordinates of the source being observed (the phase center). The data products of the correlator are the cross-correlations of the signals from each pair of stations after taking into account the modeled delay. These are called visibilities.

## 3. Scattering

This paper will not focus on scattering, however for completeness I include a brief discussion. VLBI with its high angular resolution can measure the scattering disc size, constraining the distribution of scattering material along the line of sight to the source. Scattering increases strongly (as $`\lambda ^2`$) with wavelength.

While an image that exhibits the scattering disc can be reconstructed with the set of visibilities measured for a source, a more accurate measurement can be made by modeling the visibilities directly. The amplitude of a visibility is the flux density of the source being observed multiplied by the correlation coefficient, which for a point source would be unity. Earth rotation causes the projected antenna separations (baseline lengths) to change, probing a wide range of baselines. A detection of a scattering disc manifests itself as a drop in the correlation coefficient as a function of baseline length. The scattering-based decorrelation phenomenon can be useful for probing the ISM, but is generally a nuisance for VLBI at low frequencies as it reduces the effective resolution.

## 4. VLBI astrometry

Astrometry is the science of making precision measurements of positions of astronomical sources. Multi-epoch astrometry leads to measurements of proper motions and perhaps parallaxes. The phases of the complex visibilities contain astrometric information. In the absence of correlator model errors or other deviations from perfect calibration, the visibility phase for a point source is given by

$$\varphi (\nu )=\frac{\nu }{c}\left(ul+vm\right),$$
(1)
where $`c`$ is the speed of light, $`\nu `$ is the observing frequency, $`(u,v)`$ is the projected baseline vector, and $`(l,m)`$ is the displacement vector of the source relative to the correlator model phase center. For a source that is exactly at the correlator phase center, this phase is zero. With multiple measurements of $`\varphi `$ made with different projected baselines, $`(u,v)`$, the source position, $`(l,m)`$, can be obtained.

Pulsar VLBI measurements are made almost exclusively in the 20 cm observing band. Pulsars have steep spectral indexes, making observing at higher frequencies more difficult. The ionosphere, scattering, interference, and lower angular resolution all contrive to make observing at frequencies below 1 GHz difficult. Pulsar gating is often employed to improve the signal-to-noise of the pulsar measurements. Gating disables accumulation in the correlator during a pulsar’s off-pulse phase.

### 4.1. Phase-referencing

The ionosphere, just like the ISM, contains a non-homogeneous distribution of electrons which impart a dispersive delay to the incoming radiation. The ionosphere above each station can change substantially on timescales less than 1 minute, causing the relative phases to change, corrupting the astrometry. The bulk of this can be removed by a technique known as phase-referencing (Shapiro et al. 1979). In this process, a nearby point-like calibrator source with a well known position is observed alternately with the target with a cycle time of 2 to 4 minutes. Since the calibrator’s position is well known, its visibility phases can be assumed to be purely due to calibration problems. These phases are then interpolated across the pulsar’s observation and subtracted from the pulsar’s corresponding visibility phases. This standard observing technique calibrates away most of the ionosphere’s effect, but the uncorrected ionosphere gradients often leave significant and systematic phase errors.

### 4.2. Ionospheric calibration

After phase-referencing at frequencies lower than about 5 GHz, the remaining phase errors are dominated by the effects of the ionosphere. The dispersive nature of the ionosphere can be exploited in its removal. The unwanted phase as a function of frequency for a given visibility can be expressed as

$$\varphi ^{\mathrm{iono}}(\nu )=\frac{B}{\nu },$$
(2)
where B is the strength of the uncalibrated ionosphere. The phase-referenced visibility phase then has the form

$$\varphi ^{\mathrm{vis}}(\nu )=\frac{\nu }{c}\left(ul+vm\right)+\frac{B}{\nu }.$$
(3)
The difference in frequency dependence of the first, wanted term and the second, unwanted term make it possible to measure $`B`$ and remove its effect from the visibility phase, allowing precise astrometry. Additional details of this technique and a discussion of its applicability can be found in Brisken et al. (2000 & 2002). A project to measure the parallaxes of 10 pulsars was started in 1999 after success of this technique was demonstrated on pulsar B0950+08 (Brisken et al. 2000). In this experiment, the NRAO111The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Very Long Baseline Array (VLBA) was used in the 1400–1730 MHz band. Eight frequency sub-bands, each 8 MHz wide, were placed across 330 MHz. Five epochs of observation were planned for each pulsar over the course of one year. One pulsar, J2145-0750, was not detected in either of its first two epochs and was subsequently abandoned. Successful parallax measurements were made for the nine remaining pulsars with a typical precision of 100 microarcseconds (Brisken et al. 2002). See Table 1 for the results.

### 4.3. In-beam calibration

A second method to improve phase-referencing is to use a phase-reference calibrator source that is very close to the target. In the case that a compact source of sufficient flux density (about 5 mJy for the VLBA) is within the same primary beam as the target (about 25 arcminutes for the VLBA at $`\lambda 20`$ cm), both can be observed simultaneously. At correlation time, each source is correlated separately with its own phase center. The simultaneous observation of the calibrator and target mean that no temporal interpolation is needed and the spatial proximity implies that ionospheric gradients contribute little to the phase-referenced visibility phase. Finding a suitable in-beam calibrator is often difficult or impossible. Once a good calibrator is found, a successful parallax measurement is almost guaranteed.

### 4.4. Disc-based recording

VLBI has traditionally relied on large (in physical dimension and data capacity) magnetic tapes for storage and transportation of data from the antennas to the correlator. Market pressure over the last decade has forced the price of commodity IDE discs to fall below that of magnetic tape. Discs have additional benefits over tapes for reasons including: random access to data is possible, much less maintenance is required, recording rates can be greater and capacity scales with the consumer market. The pulsar community may greatly benefit from the deployment of disc-based VLBI recorders at most large radio telescopes in the near future.

Haystack Observatory is leading the Mark 5 disc-based recorder project. To date, roughly 20 prototype units are deployed. The VLBA achieved first fringes with its first Mark 5 prototype on 13 Oct. 2002. Pulsar B0950+08 and its phase reference calibrator were observed with magnetic tape at the St. Croix and Hancock, NH VLBA stations and onto a Mark 5 recorder at the Pie Town, NM VLBA station. This data point augments 8 previous epochs and is shown in Fig. 1. The updated parallax for B0950+08 is shown in Table 1.