text:
# Interfacial phase transitions in a far-from-equilibrium magnetic growth model

## 1 Introduction

The interaction of a saturated gas in contact with a wall or a substrate may result in the occurrence of very interesting wetting phenomena, where a macroscopically thick liquid layer condenses at the wall, while the bulk fluid may remain in the gaseous phase . The wetting of solid surfaces by a fluid is a phenomenon of primary importance in many fields of practical technological applications (lubrication, efficiency of detergents, oil recovery in porous material, stability of paint coatings, interaction of macromolecules with interfaces, etc. ). Surface enrichment or wetting layers have been observed experimentally in a great variety of systems, such as e.g. polymer mixtures , adsorption of simple gases on alkali metal surfaces , with the recent addition of Hg to the adsorption species exhibiting this kind of transition phenomena , also hydrocarbons on mica , etc.

The study of wetting transitions at interfaces has also attracted considerable theoretical interest , involving, among others, different approaches such as the mean field Ginzburg-Landau method , transfer matrix and Pfaffian techniques , density matrix renormalization group methods , solving the Cahn-Hilliard equation , using Molecular Dynamic simulations , solving self-consistent field equations , and by means of extensive Monte Carlo simulations .

So far, the considerable progress due to all these studies of wetting transitions have been achieved for systems under equilibrium conditions. In contrast, the study of wetting phenomena under nonequilibrium conditions has received much less attention. For instance, Hinrichsen et al. have recently introduced a nonequilibrium growth model of a one-dimensional interface interacting with a substrate. The interface evolves via adsorption-desorption processes, which depart from detailed balance. Then, changing the relative rates of these processes a transition from a binding to a nonbinding phase has been reported .

Within this context, the aim of this work is to perform an extensive numerical study of the irreversible growth of a magnetic material confined between parallel walls where competing surface magnetic fields act. For this purpose, a variant of the irreversible Eden growth model , in which particles are replaced by spins that can adopt two different orientations, is investigated. It is shown that the interplay between confinement and growth mode leads to a physically rich phase diagram. It should also be remarked that, although the discussion is presented here in terms of a magnetic language, the relevant physical concepts can rather straightforwardly be extended to other systems such as fluids, polymers, and binary mixtures. Apart from the fundamental interest to understand this complex physical situation, it may well play a key role in the development of technologies such as micromagnetic materials, micro-fluidics, self-assembly of three-dimensional structures, adhesives, lubricants, and coatings, among others. Indeed, wetting phenomena under far-from-equilibrium conditions are expected to be of wide application to describe a great variety of processes actually encountered in practice.

Furthermore, the proposed study establishes a link with recent investigations of irreversible growth processes. In fact, the study of growth systems under far-from-equilibrium conditions is a subject that has attracted great attention during the last decades. Nowadays, this interdisciplinary field has experienced a rapid progress due to both, its interest in many subfields of physics, chemistry and biology, as well as by its relevance in numerous technological applications. Recent progress in our understanding of growth phenomena, with special emphasis on the properties of rough interfaces, has extensively been reviewed .

Also, the study of wetting phenomena in far-from-equilibrium systems under confinement has an extra ingredient of theoretical interest due to the delicate interplay between surface and bulk properties. Indeed, from the experience gained studying equilibrium systems, it is well known that, using confined geometries with restricted dimensionality, the effects of statistical fluctuations are more pronounced , leading to a new and rich physical behavior, which eventually may be the precursor of the actual critical behavior only observed in the thermodynamic limit. Within this context, in the present work it is shown that in far-from-equilibrium systems, the subtle interplay between finite-size effects, wetting, and interface growth mechanisms leads to more rich and complex physical features than in the equilibrium counterpart. In fact, a complex phase diagram that exhibits a localization-delocalization transition in the interface that runs along the walls and a change of the curvature of the growing interface running perpendicularly to the walls, is evaluated and discussed, firstly for finite-size systems, and subsequently for the extrapolated infinite system.

This manuscript is organized as follows: in Section 2 we give details on the simulation method, Section 3 is devoted to the presentation and discussion of the results, while the conclusions are finally stated in Section 4.

## 2 The model and the simulation method

In the classical Eden model  on the square lattice, the growth process starts by adding particles to the immediate neighborhood (the perimeter) of a seed particle. Subsequently, particles are sticked at random to perimeter sites. This growth process leads to the formation of compact clusters with a self-affine interface . The magnetic Eden model (MEM)  considers an additional degree of freedom due to the spin of the growing particles. Early studies of the MEM have been performed using a single seed placed at the center of the sample , but some subsequent investigations  have adopted instead $`(d+1)`$dimensional rectangular geometries. Following the latter approach, in the present work the MEM in $`(2+1)`$dimensions is studied using a rectangular geometry $`L\times L\times M`$ (with $`ML`$). Figure 1 illustrates the general setup assumed. The location of each site on the lattice is specified through its rectangular coordinates $`(i,j,k)`$, ($`1i,jL`$, $`1kM`$). The starting seed for the growing cluster is a plane of $`L\times L`$ parallel-oriented spins placed at $`k=1`$ and cluster growth takes place along the positive longitudinal direction (i.e., $`k2`$). Periodic boundary conditions are chosen along one of the transverse directions (say in the $`i`$direction), while open boundary conditions are adopted along the remaining transverse direction. Competing surface magnetic fields $`H>0`$ ($`H^{}=H`$) acting on the sites placed at $`j=1`$ ($`j=L`$) are considered. Then, assuming that each spin $`S_{ijk}`$ may adopt two possible orientations, namely up and down (i.e. $`S_{ijk}=\pm 1`$), clusters are grown by selectively adding spins to perimeter sites, which are defined as the nearest-neighbor (NN) empty sites of the already occupied ones. Considering a ferromagnetic interaction of strength $`J>0`$ between NN spins, the energy $`E`$ of a given configuration of spins is given by

$$E=\frac{J}{2}\left(\underset{ijk,i^{^{}}j^{^{}}k^{^{}}}{}S_{ijk}S_{i^{^{}}j^{^{}}k^{^{}}}\right)H\left(\underset{ik,\mathrm{\Sigma }_1}{}S_{i1k}\underset{ik,\mathrm{\Sigma }_L}{}S_{iLk}\right),$$
(1)
where $`ijk,i^{^{}}j^{^{}}k^{^{}}`$ means that the summation in the first term is taken over all occupied NN sites, while $`ik,\mathrm{\Sigma }_1`$, $`ik,\mathrm{\Sigma }_L`$ denote summations carried over occupied sites on the surfaces $`\mathrm{\Sigma }_1`$, $`\mathrm{\Sigma }_L`$ (defined as the $`j=1`$ and $`j=L`$ planes, respectively). Thus, setting the Boltzmann constant equal to unity ($`k_B1`$) and measuring absolute temperature, energy, and magnetic fields in units of $`J`$, the change of energy $`\mathrm{\Delta }E`$ involved in the addition of a spin $`S_{ijk}`$ to the system is given by

$$\mathrm{\Delta }E=S_{ijk}\left(\underset{ijk,i^{^{}}j^{^{}}k^{^{}}}{}S_{i^{^{}}j^{^{}}k^{^{}}}+H\left(\delta _{j1}\delta _{jL}\right)\right),$$
(2)
where the summation $`ijk,i^{^{}}j^{^{}}k^{^{}}`$ is taken over occupied NN sites keeping $`i,j,k`$ fixed, and $`\delta _{j1}`$, $`\delta _{jL}`$ are standard Kronecker delta symbols. Therefore, the probability for a perimeter site to be occupied by a spin $`S_{ijk}`$ is proportional to the Boltzmann factor $`\mathrm{exp}(\frac{\mathrm{\Delta }E}{T})`$, where $`\mathrm{\Delta }E`$ is given by Eq.(2). At each step, the probabilities of adding up and down spins to a given site have to be evaluated for all perimeter sites.

After proper normalization of the probabilities, the growing site and the orientation of the spin are determined through standard Monte Carlo techniques. Although both the interaction energy and the Boltzmann probability distribution considered for the MEM are similar to those used for the Ising model with surface magnetic fields , it must be stressed that these two models operate under extremely different conditions, namely the MEM describes the irreversible growth of a magnetic material and the Ising model is suitable for the study of a magnetic system under equilibrium conditions. In the MEM, the position and orientation of all deposited spins remain fixed. The nonequilibrium nature of the MEM is clear from the fact that the extensive thermodynamic variables (such as energy, entropy, and volume, for instance) grow monotonically with time and tend to diverge. Furthermore, during the growth process, the system develops a rough growth interface and evolves mainly along the longitudinal direction $`k`$ (see figure 1). However, some lattice sites could remain empty even well within the system’s bulk. Since at each growth step all perimeter sites are potential candidates for being occupied by the next spin to be added, these holes become gradually filled. It may appear that the $`k`$-coordinate is something like the “time” in a kinetic Ising model. However, this is not strictly true because already deposited spins at position $`k`$ effectively affect spin growth at position $`k^{}<k`$, and this would mean causality violation.

Far behind the active growth interface, the system is compact and frozen. When the growing cluster interface is close to reach the limit of the sample ($`k=M`$) one can compute the relevant properties of the irreversibly frozen cluster’s bulk (in the region where the growing process has definitively stopped), thereafter erase the useless frozen bulk, and finally shift the growing interface towards the lowest possible coordinate $`k`$. Hence, by repeatedly applying this procedure the growth process is not limited by the lattice length $`M`$.

It should be noticed that this paper involves a large computational effort. On the one hand, as will be seen below, the observables of interest (e.g. the susceptibility) are averaged over many transverse planes of size $`L\times L`$. In order to obtain acceptably small statistical errors, averages over $`10^{56}`$ planes in the stationary regime are typically required. So, in the present work clusters having up to $`10^9`$ spins have been grown. On the other hand, the update algorithm is quite slow as compared with standard Ising simulations, since the growing probability has to be computed after each deposition event.

## 3 Results and discussion

Recent investigations  have shown that the magnetic Eden growth process in a stripped $`(d+1)`$dimensional geometry (with $`d=1,2`$) is characterized by an initial transient of average length $`l_{Tr}`$, followed by a nonequilibrium stationary state that is independent of the starting seed. It has also been shown that the MEM in $`d=1`$ is noncritical (i.e., it only exhibits an ordered phase at $`T=0`$), while the MEM in $`d=2`$ undergoes an order-disorder thermal transition of second order at the temperature $`T_c=0.69\pm 0.01`$. Moreover, the critical exponents associated with the continuous phase transition have been found to agree, within error bars, with those of the Ising model in two dimensions. Hence, the reported findings have led to the conjecture that the $`(d+1)`$dimensional MEM and the $`d`$dimensional Ising model behave identically (unless finite-size differences that vanish in the thermodynamic limit) at criticality for all $`d`$ .

At this stage it is appropriate to briefly recall that a confined Ising film with competing surface fields undergoes an equilibrium wetting transition. Indeed, when an Ising film is confined between two competing walls a distance $`L`$ apart from each other, so that the surface magnetic fields ($`H`$) are of the same magnitude but opposite direction, it is found that the competing fields cause the emergence of an interface that undergoes a localization-delocalization transition. This transition shows up at an $`L`$dependent temperature $`T_w(L,H)`$ that is the precursor of the true wetting transition temperature $`T_w(H)`$ of the infinite system .

In view of the nontrivial correspondence established between the $`(d+1)`$dimensional MEM and the $`d`$dimensional Ising model, it should also be expected an Ising-like wetting transition for the MEM. In fact, by applying surface magnetic fields of opposite signs to the MEM, it should be expected to obtain a well defined phase transition curve between wet and nonwet states on the $`HT`$ plane. In order to deal with a phase transition that will remain in the thermodynamic limit ($`L\mathrm{}`$), one should devote attention to the $`(d+1)`$dimensional MEM with $`d2`$ since, as already pointed out, the MEM is noncritical for $`d=1`$. For this purpose, we have studied the $`(2+1)`$dimensional MEM with magnetic fields $`H`$ and $`H^{}=H`$ applied to the surfaces $`\mathrm{\Sigma }_1`$ and $`\mathrm{\Sigma }_L`$, respectively (recall Section II). As in previous investigations , the mean transverse magnetization is defined as

$$m(k,L,T,H)=\frac{1}{L^2}\underset{i,j=1}{\overset{L}{}}S_{ijk}$$
(3)
for $`k>l_{Tr}`$, in order to exclude the initial transient. Furthermore, it is assumed that the finite-size ($`L`$-dependent) susceptibility can be defined in terms of order parameter fluctuations in the same manner as for equilibrium systems, namely

$$\chi =\frac{L^2}{T}\left(m^2|m|^2\right),$$
(4)
where $`\mathrm{}`$ means the average taken over a sufficiently large number of transverse planes in the stationary regime. Then, using a standard procedure , the localization-delocalization transition curve (on the $`HT`$ plane) corresponding to the up-down interface running along the walls can be obtained considering that a point with coordinates $`(H_w,T_w)`$ on this curve maximizes $`\chi (H,T)`$. Figure 2 shows plots of $`\chi `$ vs $`T`$ for several values of $`H`$ and the fixed lattice size $`L=12`$, illustrating the method used to trace the size-dependent localization-delocalization transition curve, which is shown in figure 3 (open squares). As in the case of the Ising model, this quasi-wetting transition refers to a transition between a nonwet state that corresponds to a localized interface bound to one of the confinement walls, and a wet state associated to a delocalized domain interface centered between roughly equal domains of up and down spins. The localization-delocalization transition in a confined system is indeed the precursor of the true wetting transition that occurs in the thermodynamic limit . In fact, it is observed a finite jump in the wetting layer thickness that takes place as a result of the finite size of the system. As the lattice size is increased, the magnitude of the jump grows and diverges in the $`L\mathrm{}`$ limit, as expected for a continuous wetting transition.

Let us now discuss the critical temperature associated to the bulk order-disorder phase transition. As well known from finite-size scaling theory, there is some degree of arbitrariness in locating the $`L`$dependent critical temperature $`T_c(L)`$ of a finite system. However, the critical point $`T_c`$ of the infinite system, obtained by extrapolating $`T_c(L)`$ to the $`L\mathrm{}`$ limit, is unique and independent of any particular choice for the finite-size critical point.

Let us first consider the case with $`H=0`$, defining the $`L`$dependent critical temperature as given by the peak of the susceptibility at zero surface field, to assure consistence with the evaluation of the localization-delocalization quasi-wetting transition. Indeed, under this assumption, the quasi-wetting curve $`T_w(L,H)`$ and the critical point $`T_c(L)`$ coincide by definition at $`H=0`$. For $`L=12`$, the critical point so defined is $`T_c(L=12)=0.84`$, and is shown in figure 2 by a vertical straight line. Using larger and larger lattices, correspondingly smaller finite-size critical points are found, which tend to the actual critical point of the $`(2+1)`$-dimensional MEM, namely $`T_c=0.69\pm 0.01`$ . Before exploring the case $`H>0`$, it is appropriate to recall a long-standing discussion on this topic generated in the field of equilibrium critical phenomena. Parry and Evans  claim that the critical temperature for a finite system depends on the surface magnetic fields, and only differs from the wetting temperature of the infinite system ($`T_w(H)`$) by a term that vanishes in the thermodynamic limit. Indeed, they suggest a scaling Ansatz such that $`T_w(H)T_c(L,H)`$ is of order $`L^{1/\beta _s}`$, where $`\beta _s`$ is the exponent that describes the growth of the wetting layer. However, Swift et al.  and Indekeu et al.  propose that $`T_c(L,H)`$ is actually a shifted wetting transition (hence called quasi-wetting transition $`T_w(L,H)`$ ), which is different from the bulk critical point $`T_c(L)`$, and such that $`T_w(L,H)`$ tends for $`L\mathrm{}`$ to the actual wetting temperature $`T_w(H)`$. In the case of the MEM, general considerations, supported by our numerical results, appear to favor the latter proposal, so that the bulk critical point is independent of $`H`$ and hence clearly different from the quasi-wetting temperature. However, it should be remarked that the controversy has been established for systems under equilibrium and the present study of the MEM corresponds to far from equilibrium conditions.

Let us first consider an increase in the surface fields from $`H=0`$ for a system within the ordered phase (i.e., for $`T<T_c(L)`$). Since fluctuations in the bulk are governed only by the temperature, it turns clearly out that the bulk will remain in its ordered state irrespective of the applied field. Indeed, as will be discussed below in detail, an increase in the fields favors the formation of a stable longitudinal interface between domains of up and down spins. However, within each domain, it is clear that the state of order will depend only in fluctuations driven by the temperature. Hence, surface fields applied on an ordered system below $`T_c(L)`$ will eventually favor the coexistence of oppositely oriented ordered domains, but are not capable of generating disorder within each domain. These arguments are strongly supported by our simulations. For instance, figures 4(a)-(c) show typical snapshot configurations that correspond to nearly the same temperature (below $`T_c(L=32)=0.76`$) and several different magnetic fields. As expected from our considerations, the fields appear to support the formation of the longitudinal interface between opposite spin domains, but do not affect the bulk ordered state within each domain. In particular, it should be noticed that figure 4(c) corresponds to a field well above the corresponding one on the quasi-wetting curve (note that the snapshots correspond to $`L=32`$, and the associated transition curves are shifted to the left with respect to those for $`L=12`$, shown in the phase diagram of figure 3). If $`T_w(L,H)`$ would be the system’s critical point, it should be expected a system beyond the quasi-wetting curve to be disordered, in remarkable contrast with the ordered configuration shown by figure 4(c). Moreover, as a check of consistency, one can compare the configurations shown in figures 4(c) and 4(e), that correspond to nearly the same fields and differ in temperature. It turns thus evident that the bulk’s order-disorder phase transition occurs at a temperature far away from the quasi-wetting transition curve, and consistent with $`T_c(L=32)=0.76`$.

Further insights on the role of $`H`$ acting within the bulk ordered phase $`T<T_c(L)`$ can be gained by means of the following procedure. Let us point our attention to the stationary regime considering all completely filled columns directed along the $`j`$direction, which are formed by $`L`$ spins and are identified through the values of the remaining coordinates $`i`$ and $`k`$.

For any given column $`(i,k)`$, a bond to each pair of nearest-neighbor sites occupied by oppositely oriented spins is assigned. Summing over the whole column, $`n_b(i,k)`$ is defined as the total number of bonds for that column, so that $`n_b=0`$ for parallel-oriented spins and $`n_b=L1`$ for alternating up-down nearest-neighbor spins. Since all columns are statistically independent, the system is allow to grow for a sufficiently long time and averages are taken over all filled columns. In this way, the normalized bond probability distribution $`P(n_b)`$ can be computed as a function of temperature, surface magnetic fields, and lattice size. For the purposes of the present discussion, it suffices to fix $`L=12`$ and consider the effects of increasing the fields for a given value of temperature below $`T_c(L=12)=0.84`$. Figure 5 shows the bond probability distribution $`P(n_b)`$ vs $`n_b`$ for $`T=0.6`$ and several values of $`H`$. It is observed that $`P(n_b)0`$ for $`n_b2`$, irrespective of the field. Hence, it is concluded that the system remains in its ordered state independently of $`H`$, and that the role of the magnetic field is that of driving the system from a state constituted by a single domain ($`P(n_b=0)1`$) to a state formed by two oppositely oriented ordered domains ($`P(n_b=1)1`$). Figure 6 shows the mean number of bonds per column $`n_b`$ as a function of $`H`$, for three different temperatures. As expected, in all cases the field drives the crossover from a single ordered domain to two opposite ordered domains.

So, our discussion concerning the location of the critical temperature associated to the bulk order-disorder phase transition can be summarized by stating that it is found a compelling evidence for interpreting $`T_w(L,H)`$ as a quasi-wetting transition, clearly different from the finite-size critical temperature $`T_c(L)`$. Furthermore, our results are consistent with the assumption of a field-independent critical point $`T_c(L)`$, since the magnetic field appears to play no role in the state of order within each magnetic domain.

Since the MEM is a nonequilibrium kinetic growth model, it also allows the identification of another kind of phase transition, namely a morphological transition associated with the curvature of the growing interface of the system . To avoid confusion, it shall be remarked that the term growing interface is used here for the transverse interface between occupied and empty lattice sites, while it was used above for the longitudinal interface between up and down spin domains.

Firstly, let us consider a longitudinal slice with a fixed value of $`i`$ in the range $`1iL`$. In order to define the location of the growing interface at time $`t`$, it is assumed that each row contributes to the growing interface with the outermost perimeter site (i.e., the site with the largest value of the longitudinal coordinate $`k`$, for a given row number $`j`$) and that number is called $`I_j(t)`$. Then, the growing interface center of mass, that is taken as the location of the growing interface at time $`t`$, $`I(t)`$, is given by

$$I(t)=\frac{1}{L}\underset{j=1}{\overset{L}{}}I_j(t).$$
(5)
Subsequently, one can evaluate the coordinates of the growing interface relative to its center of mass location at time $`t`$, namely $`IR_j(t)I_j(t)I(t)`$, for $`j=1,2,\mathrm{},L`$. In this way, it is possible to describe the growing interface at any time $`t`$ during the growing process just by evaluating the set $`\{IR_j(t)\}`$. However, one should be cautious at this point. In fact, since the applied surface fields are of equal intensity but have opposite direction, it turns out that the probability of occurrence of a given growing interface $`\{IR_j\}`$ must equal the one corresponding to $`\{IR_j^{^{}}\}`$, where $`j^{^{}}=L+1j`$. But then, unless the growth profile happens to be symmetric (i.e., invariant under $`jL+1j`$, for all $`j`$), the time average of equally probable growing interfaces $`\{IR_j\}`$, $`\{IR_j^{^{}}\}`$ will lead to an unphysical symmetrized profile that is not representative of the actual shape of the growing interface. To avoid this problem, the following procedure is used. First, the largest value of the longitudinal coordinate $`k`$ that corresponds to a completely filled column is located. Then, by means of the sign of the total magnetization of that column, i.e. $`Ssign\left(_jS_{ijk}\right)`$, the orientation of the dominant spin domain in the active growing interface is identified. Supposing that, following the recipe given above, a given profile $`\{IR_j(t)\}`$ is obtained, then the growing interface location is redefined as $`\{IR_j^{}(t)\}`$, where $`j^{}j`$ if $`S=+1`$ and $`j^{}j^{^{}}=L+1j`$ if $`S=1`$, for all $`j`$. Notice that $`j^{}=1`$ ($`j^{}=L`$) corresponds to the side of dominant (non-dominant) spin domain, while $`j=1`$ ($`j=L`$) is the side of positive (negative) magnetic field.

Then, it is possible to evaluate the average relative growing interface $`IR_j^{}`$ by taking into account interface coordinates measured at different times between $`t_i`$ and $`t_f`$, and also by averaging all longitudinal $`i`$fixed slices, i.e.

$$IR_j^{}=\frac{1}{L}\frac{1}{(t_ft_i+1)}\underset{i=1}{\overset{L}{}}\underset{t=t_i}{\overset{t_f}{}}IR_{i,j^{}}(t).$$
(6)
Figure 7 shows $`IR_j^{}`$ versus $`j^{}`$ for different values of the surface magnetic field $`H`$, for a fixed temperature $`T=0.6`$ and a fixed lattice size $`L=32`$. From the figure it follows that three qualitatively distinct growth regimes can clearly be distinguished. Indeed, it is observed that, while for small fields the system grows with convex curvature, increasing the fields the growth process enters into a regime of non-defined curvature, since the dominant spin domain partially wets the confinement wall, while the non-dominant domain does not. But then, further increasing the fields, a point is reached where the non-dominant spin domain also (partially) wets the wall and the growing interface turns concave. This qualitative behavior has been observed for all temperatures and lattice sizes within the range of interest of this work.

To explore this phenomenon quantitatively, the behavior of the contact angles between the growing interface and the confinement walls (as functions of temperature and magnetic field) has to be studied thoroughly. Clearly, two different contact angles must be defined, namely $`\theta _D`$ for the angle corresponding to the dominant spin cluster, and $`\theta _{ND}`$ for the one that corresponds to the non-dominant spin cluster. Figure 8 shows plots of $`\mathrm{cot}(\theta )`$ vs $`H`$ for $`T=0.9`$ and $`L=12`$. The vertical dashed lines indicate the fields that separate a given growth regime from another one. One observes that, increasing the surface fields, the growth regime changes from convex to non-defined to concave, in agreement with the interface profiles plotted in figure 7. Analogously, figures 9(a)-(d) show plots of $`\mathrm{cot}(\theta )`$ vs $`T`$ for $`L=12`$ and several different values for the magnetic field $`H`$. Again, vertical dashed lines correspond to transition temperatures between different growth regimes. Figure 9(a) corresponds to $`H=0.2`$ and displays the characteristic behavior for very small magnetic fields, that is, a convex growing interface irrespective of temperature. For $`H=0.4`$ one observes a single transition from the growth regime of non-defined curvature to the convex growth regime, which shows up by increasing the temperature, as shown in figure 9(b). It should be noticed that the concave growth regime is prevented, since for small enough magnetic fields $`\mathrm{cot}(\theta _{ND})<0`$ for all $`T`$. As the fields are increased, $`\mathrm{cot}(\theta _{ND})`$ moves upwards and crosses $`\mathrm{cot}(\theta _{ND})=0`$, as expected from the plot of figure 8. For instance, the plots of $`\mathrm{cot}(\theta )`$ vs $`T`$ for $`H=0.6`$, shown in figure 9(c), exhibit this behavior. Hence, here one has to deal with three transition temperatures. Finally, by further increasing the fields, the whole low-temperature region is dominated by the concave growth regime and two transition temperatures remain, as shown in figure 9(d) for $`H=1.5`$. All these features are compactly shown on the $`HT`$ phase diagram of figure 3, where open (filled) circles refer to the transition between non-defined and concave (convex) growth regimes.

As shown in figure 3, the phase diagram of the MEM in a confined geometry with competing surface fields is very rich and exhibits eight regions. In order to gain some insight into the physics involved in this phase diagram, some typical snapshot configurations characteristic of the various different growth regimes observed are shown in figure 4, as obtained using lattices of size $`L=32`$.

To begin with, let us analyze Region $`I`$ (see figure 3), that corresponds to the Ising-like nonwet state and the convex growth regime. In this region, temperature is low and the system grows in an ordered state, i.e. the dominant spin domain prevails and the deposited particles tend to have their spins all pointing in the same direction. Small clusters with the opposite orientation may appear preferably on the surface where the non-dominant orientation field is applied. These “drops” might grow and drive a magnetization reversal, thus changing the sign of the dominant domain. In fact, the formation of sequences of well-ordered domains are characteristic of the ordered phase of confined (finite-size) spin systems. For instance, this phenomenon has already been observed in finite Ising strips  and magnetic Eden thin films . Due to the fact that open boundary conditions are imposed at $`j=1`$ and $`j=L`$, perimeter sites at the confinement walls experience a missing neighbor effect, that is, the number of NN sites is lower than for the case of perimeter sites on the bulk. Since the surface magnetic fields in this region are too weak to compensate this effect, the system grows preferentially along the center of the sample as compared to the walls, and the resulting growth interface exhibits a convex shape. A typical snapshot configuration characteristic of Region $`I`$ is shown in figure 4(a).

Let us now consider an increase in the fields, such as the system may enter into Region $`II`$ (see figure 3). Since the temperature is kept low, the system is still in its ordered phase and neighboring spins grow preferably parallel-oriented. The surface fields in this region are stronger and thus capable of compensating the missing NN sites on the surfaces. But, since the fields on both surfaces have opposite signs, it is found that, on the one hand, the field that has the same orientation as the dominant spin cluster favors the growth of surface spins, while on the other hand, the sites on the surface with opposite field have a lower probability to be chosen during the Monte Carlo growth process. Hence, the contact angle corresponding to the dominant spin cluster is then $`\theta _D<\frac{\pi }{2}`$, while the non-dominant is $`\theta _{ND}>\frac{\pi }{2}`$. Thus, on the disfavored side the growing interface becomes pinned and the curvature of the growing interface is not defined. Figure 4(b) shows a typical snapshot corresponding to Region $`II`$.

Keeping $`H`$ fixed within Region $`II`$ but increasing the temperature, thermal noise will enable the formation of drops on the disfavored side that eventually may nucleate into larger clusters as the temperature is increased even further. This process may lead to the emergence of an up-down interface, separating oppositely oriented domains, running in the longitudinal direction (i.e. parallel to the walls). Since sites along the up-down interface are surrounded by oppositely oriented NN spins, they have a low growing probability. So, in this case the system grows preferably along the confinement walls and the growing interface is concave (figure 4(c)). Then, as the temperature is increased, the system crosses to Region $`A`$ (see figure 3) and the onset of two competitive growth regimes is observed, namely: (i) one exhibiting a non-defined growing curvature that appears when a dominant spin orientation is present, as in the case shown in figure 4(b); (ii) another that appears when an up-down interface is established and the system has a concave growth interface, as is shown in figure 4(c). Further increasing the temperature and for large enough fields, the formation of a stable longitudinal up-down interface that pushes back the growing interface is observed. So, the system adopts the concave growth regime (see figure 4(c) corresponding to Region $`IV`$ in figure 3). Increasing the temperature beyond $`T_c(L)`$, a transition from a low-temperature ordered state (Region $`IV`$) to a high-temperature disordered state (Region $`VI`$, see figure 4(e)), both within the concave growth regime, is observed. Analogously, for small enough fields, a temperature increase drives the system from the ordered convex growth regime (Region $`I`$) to the disordered convex growth regime (Region $`V`$, see figure 4(d)). As shown in figure 3, there is also an intermediate fluctuating state (Region $`B`$) between Regions $`V`$ and $`VI`$, characterized by the competition between the disordered convex growth regime and the disordered concave one.

Finally, a quite unstable and small region (Region $`III`$) that exhibits the interplay among the growth regimes of the contiguous regions, can also be identified. Since the width of Region $`III`$ is of the order of the rounding observed in $`T_c(L)`$, large fluctuations between ordered and disordered states are observed, as well as from growth regimes of non-defined curvature to convex ones. However, figure 4(f) shows a snapshot configuration that is the fingerprint of Region $`III`$, that may prevail in the thermodynamic limit, namely a well defined spin up-down interface with an almost flat growing interface.

Let us now extrapolate our results to show that the rich variety of phenomena found in a confined geometry is still present in the thermodynamic limit ($`L\mathrm{}`$), leading to the phase diagram shown in the inset of figure 3. As clearly seen by comparison with the finite-size results, the crossover Regions $`A`$ and $`B`$ collapse in this limit, so that only the six regions that correspond to well identified growth regimes (as illustrated by the snapshot configurations of figure 4) appear to remain.

In order to illustrate the extrapolation procedure, the following seven representative points of the finite-size phase diagram are discussed in detail: (i) the points labeled $`P_1`$, $`P_1^{}`$, $`P_2`$, and $`P_2^{}`$, that correspond to the intersections of the $`H=0.6`$ line with the various transition curves shown in figure 3, and (ii) the points labeled $`P_3`$, $`P_3^{}`$, and $`P_4`$, that refer to the intersection point between Regions $`I,II,III`$, and $`A`$, the minimum of the limiting curve between Regions $`IV`$-$`VI`$ and $`A`$-$`B`$, and the zero-field transition point, respectively.

Figure 10 shows plots of $`T`$ versus $`L^1`$ for $`12L48`$ corresponding to the points $`P_1,P_1^{},P_2,`$ and $`P_2^{}`$. Also shown in the figure are the fits to the data extrapolated to $`L^1=0`$. The results from the extrapolations are: $`T_1=0.67\pm 0.01`$, $`T_1^{}=0.66\pm 0.01,`$ and $`T_2=1.30\pm 0.02`$, $`T_2^{}=1.29\pm 0.01,`$ pointing out that, within error bars, $`P_iP_i^{}`$ $`(i=1,2)`$ in the $`L\mathrm{}`$ limit. Using the same procedure, the extrapolations of $`P_3`$ and $`P_3^{}`$ (not shown here) give: $`H_3=0.30\pm 0.01`$, $`H_3^{}=0.31\pm 0.02,`$ and $`T_3=0.69\pm 0.01`$, $`T_3^{}=0.71\pm 0.03`$. So, one has $`P_3P_3^{}`$ for $`L\mathrm{}`$ within error bars. Finally, the extrapolation of $`P_4`$ is $`T_4=T_c=0.69\pm 0.01`$.

Using the above-mentioned extrapolation procedure, the phase diagram in the thermodynamic limit can be drawn, as shown in the inset of figure 3. By comparison with the finite-size phase diagram of figure 3, one can notice that, as anticipated, the crossover Regions $`A`$ and $`B`$ appear in the phase diagram just as a consequence of the finite-size nature of confined geometries, since they collapse in the $`L\mathrm{}`$ limit. Moreover, we conjecture that Region $`III`$ may remain in the thermodynamic limit. Although this (very tiny!) region corresponds to a physically well characterized growth regime, since one expects that the system in this region may grow in an ordered phase with a delocalized up-down domain interface and a convex growing interface, statistical errors due to large fluctuations close to criticality hinder a more accurate location of this region. The unambiguous clarification of our conjecture remains as an open question that will require a huge computational effort.

Besides an Ising-like continuous wetting transition, coupled morphological transitions in the growing interface, which arise from the MEM’s kinetic growth process, have also been identified. Comparing the equilibrium wetting phase diagram of the Ising model  and that of the MEM, it follows that the nonequilibrium nature of the latter introduces new and rich physical features of interest: the nonwet (wet) Ising phase splits out into Regions $`I`$ and $`II`$ (Regions $`III`$ and $`IV`$), both within the ordered regime ($`T<T_c`$) but showing an additional transition in the interface growth mode. Also, the disordered state of the Ising system ($`T>T_c`$) splits out into Regions $`V`$ and $`VI`$ exhibiting a transition in the interface growth mode.

It should be noticed that we have restricted ourselves to temperatures above $`T=0.5`$ throughout, since the lower the temperature in the ordered phase, the greater the computational effort needed to reliably sample the whole configuration space (indeed, ergodicity is broken in the $`T0`$ limit). Right at $`T=0.5`$ the wetting curve of the phase diagram (inset of figure 3) intercepts the $`H`$-axis close to $`H=0.9`$. On physical grounds no particular features of interest are expected to arise in the $`T0`$ limit, and the critical field $`H=1`$ for $`T=0`$ can be inferred by energetic considerations, as e.g. in the case of the Ising model.

## 4 Conclusions

The growth of magnetic Eden clusters with ferromagnetic interactions between nearest-neighbor spins has been studied in a $`(2+1)`$dimensional geometry with competing surface magnetic fields. Extensive Monte Carlo simulations allow us to locate, on the one hand, an Ising-like localization-delocalization wetting transition, and, on the other hand, a morphological transition associated with the curvature of the growing interface. In this way, eight different regions on the $`HT`$ phase diagram for a finite-size lattice are identified. Moreover, the characteristic behavior of typical growth processes within each region are discussed, and qualitative explanations that account for the observed features are provided. Finally, extrapolating the results to the $`L\mathrm{}`$ limit the phase diagram is obtained. It is composed of six different regions, since two crossover regions identified in the finite-size phase diagram appear to collapse in the thermodynamic limit. The obtained phase diagram shows new and rich physical features of interest, which arise as a consequence of the nonequilibrium nature of the investigated model.

We hope that the presented results will, on the one hand, contribute to the understanding of the rich and complex physical phenomena exhibited by the irreversible growth of binary mixtures in confined geometries, and on the other hand, stimulate further experimental and theoretical work.

ACKNOWLEDGMENTS. This work is supported financially by CONICET, UNLP, and ANPCyT (Argentina). The authors thank the referees for fruitful comments and suggestions.


format: text

text:
# A Current Induced Transition in atomic-sized contacts of metallic alloys

## I Introduction

The scale of electronic devices will soon be reduced to a level where material properties will be very different from the bulk, necessitating research into very small scale devices. A metallic point contact is perhaps the most simple example of an atomic-sized device. It consists of a connection of a few atoms, typically between 1 and 1000, between two macroscopic electrodes. The size of the contact is therefore of the same order of magnitude as the Fermi wavelength of the electrons. In this limit, the conductance is related to the transmission of an electron wave through the system. For a free electron gas, the transmission only depends on the size of the constriction. As a result, the conductance in a 2-dimensional electron gas has been seen to increase by steps as a function of the width of the constriction. The situation is more complicated in a metallic contact, where the transmission depends on the geometry of the contact as well as on the atomic structure of the atoms forming it. In the simplest case the contact consists of a single atom, where the conductance has been shown to be determined by the valence orbitals of the atom forming the contact.

Several ways exist to create a point contact. The method used in this research is the so-called mechanically controllable break junction (MCBJ) technique. Its principle is very simple; it consists of breaking a metallic wire thereby creating a clean fracture surface. After rupture, an atomic-sized contact can be made by indenting the broken ends of the wire into each other. This contact is subjected to a repeated cycle of breaking and indenting, during which the conductance of the contact is measured. As the contact is elongated, the diameter of the constriction reduces and consequently the conductance decreases.

For most metals at small contacts, the conductance decreases in steps of the order of $`G_0=2e^2/h`$, the quantum of conductance. It has been shown that these steps are due to rearrangements of atoms in the contact. The last plateau is generally assumed to correspond to a single atom contact.

The electrical and mechanical properties of atomic point contacts made from pure metals have been extensively studied. For noble metals (Cu, Ag, Au) the last plateau of a conductance trace, corresponding to a contact of one atom, has a conductance value around 1 $`G_0`$. On the other hand, transition metals (of which we will consider Ni, Pd, Pt) have a smallest conductance value around 1.5 $`G_0`$, or higher. These values are explained by the fact that noble metal atoms have a single $`s`$ valence orbital and transition metal atoms have five $`d`$ valence orbitals in addition. Very few studies have been made of the properties of alloys at the atomic scale.

Recently, point contact studies were made of random alloys of a transition metal and a noble metal, namely gold and paladium and copper and nickel, for different concentration ratios. For AuPd, the addition of Pd only leads to a decrease in height of the Au conductance peak, and no shift is observed. We propose below a different explanation for the persistence of the Au peak up to a high concentration of Pt or Pd. The CuNi study was aimed at an investigation of the influence of impurity scatterers on the point contact conductance. For low Ni concentrations the nickel atoms can be considered as impurity scatterers for the electrons, that lead to a decrease of the conductance resulting in a shift towards lower values of the peaks in the conductance histogram. This shift grows linearly with the Ni concentration. For higher concentrations of transition metal atoms, it is difficult to predict how the contact will behave and what effects might occur. Naively, one would expect a single-atom contact to be randomly formed by either a noble (Cu, Au) or a transition metal (Ni, Pd) atom, with a probability depending on concentration. However, this does not take into account several effects that could influence the formation of the contact, such as segregation, diffusion and possibly electromigration. Indeed, here we show that this simple picture does not hold. In addition we report an unexpected transition in the conductance characteristics as a function of the bias voltage for random transition-noble metal alloys. Complementary to the experiments we have performed molecular dynamics simulations in order to obtain information on the structure and composition of the atomic point contacts that is not provided by a conductance measurement.

## II Experiment

The alloys that were used in this experiment, apart from copper-nickel, were silver-paladium and gold-platinum, all at a 1:1 concentration ratio. These materials were chosen for the following reasons: First, they are miscible at any concentration. Second, the conductance of the noble metal atoms (Cu, Ag, Au) is different from that of the transition metal atoms (Ni, Pd, Pt), (see above) which allows us, in principle, to distinguish single-atom contacts of noble metals and those of transition metal atoms.

Samples were made by arc-melting equal amounts of Cu (Ag, Au) and Ni (Pd, Pt) and quenching to room temperature. CuNi and AgPd form random alloys. Below $``$1260 C, AuPt may segregate, but we see little evidence for segregation in our data. Of these samples wires are made and annealed for several hours at 900 C. A 1.5 cm long piece of wire is then glued on a phosphor-bronze substrate by two drops of epoxy on both sides of a notch made in the middle of the wire. The wire is broken at low temperature (4.2 K) by applying a mechanical force on the substrate. The two broken ends are then brought back together and the bending is controlled by means of a piezo element allowing us to obtain a contact of any size, where the contact elongation can be controlled with an accuracy better than 0.01 Å. The conductance is measured in a four point configuration.

Individual conductance traces (plots of the conductance vs the elongation of the wire) differ from each other because the positions of the atoms are different each time a new contact is formed. For a more quantitative view of the conductance of the contact a histogram technique is used that will indicate preferred values of the conductance when the size of the contact decreases. A histogram is constructed as follows: the $`y`$-axis of an indivual conductance trace is divided into a number of bins and we count the number of data-points that fall in those bins. Peaks in the histogram represent preferred values of the conductance.

Histograms were created from $``$2500 individual conductance traces. Between these individual traces the broken ends of the wire are indented to a large contact size that is set by a predetermined value for the conductance. The usual indentation depth was to a conductance of 25–30 $`G_0`$, corresponding to a contact of approximately as many atoms.

For a further investigation of the noble-transition metal point contacts, classical molecular dynamics simulations have been performed on a CuNi point contact. The system was simulated using two-atom potentials and a Johnson-Oh embedding function formalism. The motion of the atoms in a point contact of a random alloy at low temperature was calculated while the contact was being stretched.

## III Results

In Fig. 1 histograms of AuPt recorded with four different voltages applied across the contact (bias voltage) are shown. All histograms are normalized by the area under the curves in the interval  $`G_0`$.The histogram at low bias voltage exhibits a strong peak positioned somewhat below 1 $`G_0`$. The amplitude of the peak initially decreases slightly as the bias voltage is increased but at $`V_{\mathrm{bias}}=`$ 400 mV it has disappeared and instead a peak positioned at G$``$1.9 $`G_0`$ appears. It is important to notice that this peak is located approximately at the position of the first peak of pure Pt.

Measurements of the conductance of copper-nickel and silver-paladium showed that also for these alloys the low bias voltage histogram resembles that of the noble metal having a dominant peak slightly below 1 $`G_0`$, while at high bias it resembles the conductance histogram of the transition metal and all the weight below 1.5 $`G_0`$ has disappeared.

After a high bias voltage measurement, decreasing the voltage across the contact did not result in a return to a noble-metal-like histogram, as can be seen in Fig. 2. For all alloys, only after full indentation of both ends of the wire (by putting the voltage on the piezo to zero), creating a contact of mesoscopic size ($`>`$ 1000 atoms), a noble-metal-like histogram was recovered.

From a conductance histogram it is impossible to tell whether the shift in conductance is gradual or sudden since the histogram is accumulated from many conductance scans. Therefore we decided to measure a number of subsequent histograms of a smaller number of scans (500–1000) to be able to observe possible time-evolution of the contact. Again, the results are similar for the three alloys. At intermediate bias voltage (between 300 and 400 mV) it is possible to observe the change from a noble to transition-metal-like conductance histogram over the course of a small number of histograms as can be seen in Fig. 3.

We have also studied the influence of the indentation depth of the broken ends of the wires on the transition in the contact. In this experiment, the largest contact size between individual traces of a histogram was varied for a number of histograms. This contact size between traces is a measure for the indentation depth, and is characterized by a certain conductance value. It is observed (Fig. 4) that at smaller indentation depth the change from noble to transition-metal-like appears sooner (both at lower bias voltage and in a shorter time period) than for larger indentation. The additional fine structure in the case of very small indentation (10 $`G_0`$) in Fig. 4 can be explained by the fact that at small indentation the number of possible atomic configurations is reduced compared to a large indentation, and consequently the number of conductance values is limited. At very low indentation depths this can result in a repetition of a limited number of evolution paths for the contact.

In order to determine if the transition is related to the direction of the current in the contact we have conducted the same experiment using an ac voltage instead of a dc one. The frequency of the bias voltage was $`\omega =`$10 kHz, which was limited by the frequency transfer characteristics of our setup. The current and voltage over the sample were detected using two lock-in amplifiers. We found that the results obtained from the ac-voltage experiment are similar to ones from the dc-experiments. We conclude that the transition is not related to the direction of the current, at least for time scales longer than 0.1 ms.

In the molecular dynamics simulations, a random alloy of copper and nickel atoms was created, shaped like an hourglass, which was subjected to about 30 cycles of indenting and breaking of the contact at a temperature of 4K. The main outcome of the simulations is that the randomness of the alloy is preserved; no segregation was observed for up to 18 cycles of indentation to a contact size of 10–50 atoms and subsequent breaking. Our model does not allow us to simulate an electrical potential across the contact or a current through it. Another significant drawback is the fact that even though the timescale of the simulation was $`10^9`$ s, which is 8 orders of magnitude smaller than in the experiment, the calculation time was one or two days for a single cycle, thus limiting the number of simulations we could perform.

## IV Discussion

As is illustrated for Au and Pt in Fig. 5, at low bias voltage the histograms of the alloys resemble the histograms of the noble metals (Cu, Ag, Au) rather than those of the transition metals (Ni, Pd, Pt), with a distinct peak slightly below 1 $`G_0`$, but no peak around $``$1.5–1.8 $`G_0`$. However, when a high bias voltage is applied to the contact the picture is reversed. In fact, for bias voltages above the threshold value the conductance histogram only exhibits a peak located at a conductance value of 1.5–1.8 $`G_0`$, and is therefore similar to a transition metal histogram.

This suggests that as a result of the bias voltage applied across the contact the noble metal atoms are expelled from it, and only transition metal atoms remain. The peak around 1 $`G_0`$ for low-bias voltage indicates the presence of noble metal atoms; the fact that there is no weight in the conductance histogram below 1.5 $`G_0`$ in case of a high-bias voltage indicates the absence of noble metal atoms.

From the time dependence of the current-induced changes we conclude that the composition of the contact is irreversibly modified, locally. Fig. 3 illustrates how this modification gradually builds up and Fig. 2 shows that the initial properties are only recovered after making a large indentation. These observations can be understood by assuming that Au (Cu, Ag) atoms are driven away from the contact at high bias, leaving a nearly pure Pt (Ni, Pd) contact. After indentation to large contacts fresh Au (Cu, Ag) atoms are anew mixed into the contact area. The importance of the irreversibility of the shift with bias voltage should not be underestimated. It indicates that the shift cannot be attributed to a straightforward effect in pure metals, since such an effect would be reversible.

This view is supported by the fact that the transition is also dependent on the indentation depth. At higher indentation depths the transition takes place at a higher bias voltage. This agrees with the larger region that has to be depleted of noble metal atoms for larger indentation depths in order to have a contact consisting purely of transition metal atoms. The number of atoms which have to be expelled therefore increases with the indentation.

There are two principle candidates for explaining the observed change in composition of the contacts when a bias voltage is applied across the contact: electromigration, and differential thermal diffusion. Electromigration is the motion of atoms in a conductor under the influence of an applied voltage. It is traditionally described in terms of a direct force and wind force. The direct force is believed to be due to the electric field and the effective charge of the atom, and the wind force is due to the scattering of the current carrying electrons resulting in a momentum transfer. More recent work suggests that it is not meaningful to distinguish these two components and that one should consider the total force to be due to an induced bias in activation barriers for migration .

Scattering of current-carrying electrons causes Joule heating of the contact  resulting in a temperature gradient leading to a diffusion of atoms away from the contact. This could lead to a change in the composition of the contact when the diffusion coefficients of the different types of atom are unequal. The type of atom with the largest diffusion coefficient would be expelled from the contact. We estimate that a bias voltage of 400 mV raises the lattice temperature locally to $``$400 K . The diffusion constants for the noble metals are usually much larger than those for the transition metals.

The experiments discussed above all support the view that the composition of the contact is changed, but unfortunately they do not shed any light on the mechanism behind it. The dependence of electromigration on the direction of the current could possibly allow us to distinguish between thermal diffusion and electromigration. In an ac-experiment the direction of the electrical force would change with the potential. When the polarity of the potential changes rapidly the average force on the atoms would be negligible, and only thermal diffusion would remain as a possible cause for the transition. The result of the ac-experiments were the same as those from the dc-experiments. Therefore we conclude that the transition in the contact is more likely due to the second effect. Electromigration is not fully excluded since the gradient in the current density would drive the atoms away from the contact area, making the return path much less effective.

The outcome of the simulations did not show a preference for the formation of noble-metal last atom contacts over transition-metal atom contacts. The fact that we could not simulate a voltage or current in the contact prevents us from investigating a possible preference for the formation of transition metal atom contacts at high bias. However, the simulations did shed some more light on the dynamics of the point contact. In the discussion we have previously assumed that the last contact would be formed by one atom between the two shoulders of the banks. In the simulations we observed that the last contact is often formed by two atoms, one in each shoulder, see Fig.6.

This changes the question of what the conductance of a noble- or transition-metal-atom contact is, to the question of what is the conductance of a contact between a noble and a transition-metal atom. Based on the work in Refs.  we speculate that the conductance of a contact between a copper and a nickel atom is determined by the copper atom, since it has the smallest number of valence orbitals or conductance channels, which would limit the conductance to a single channel with a conductance somewhat below 1 $`G_0`$. If, as the simulations suggest, most of the final contacts are formed by two atoms in series, then $``$75% of the final contacts would have approximately the conductance of copper. This may also explain why we cannot distinguish a nickel conductance peak in the CuNi low bias histogram, since one would expect only $``$25% of the contacts to be formed by two nickel atoms. Thus, the conductance histogram of CuNi would resemble the copper histogram with a peak below 1 $`G_0`$. The assumption that the conductance between two atoms is determined by the atom with the smallest number of conductance channels remains to be further investigated. However, recent work on a hydrogen molecule trapped between two Pt electrodes shows that this system has a single nearly fully transmitting conductance channel, presumably due to the single $`s`$ orbital of hydrogen .

In conclusion, we have investigated atomic point contacts made from the noble-transition metal alloys AuPt, CuNi, and AgPd at a concentration ratio of 1:1. At low bias voltage the histograms resemble the ones of the noble metals whereas this situation is reversed at high bias. Our interpretation is that the contact is depleted of noble metal atoms with high bias. The concentration can therefore be tuned with the bias voltage. The mechanism we propose relies on the strong local gradients in lattice temperture and/or current density, driving the noble-metal atoms away from the junction at high voltage bias.

We thank A. Lodder, R.H.M. Smit, C. Untiedt, and I.K. Yanson for helpful discussions, and R.W.A. Hendrikx and M.B.S. Hesselberth for assistance with the sample preparation. This research has been supported by a European Community Marie Curie Fellowship under contract number HPMF-CT-1999-00196.


text:
# Cylindrical surface plasmon mode coupling as a three-dimensional Kaluza- Klein theory

\[

## Abstract

Cylindrical surface plasmon (CSP) mode coupling in a nanowire or a nanochannel is described in terms of three-dimensional Kaluza-Klein theory in which the compactified third dimension is considered to be the angular coordinate of the metal cylinder. Higher ($`n0`$) CSP modes are shown to posses a quantized effective charge proportional to their angular momentum $`n`$. In a nanowire these slow moving charges exhibit strong long-range interaction via fast massless CSPs with zero angular momentum ($`n=0`$). Such a mode-coupling theory may be used in description of nonlinear optics of cylindrical nanowires and nanochannels (for example, in single-photon tunneling effect), and may be further extended to describe interaction of electrons with nonzero angular momenta in thin gold wires and carbon nanotubes.

PACS no.: 78.67.-n, 11.10.Kk, 41.20.Jb

\]

Quantum communication and information processing require unusual sources of light with strong quantum correlations between single photons. This necessity drives strong current interest in exploration and understanding of nonlinear optical phenomena at single photon levels. One example of such phenomena is photon blockade in a nonlinear optical cavity  introduced in close analogy with the well-known phenomenon of Coulomb blockade for quantum-well electrons. Another recent example is the proposed observation of sound waves and, possibly, superfluid state in a weakly interacting ”photon gas” in a nonlinear Fabry-Perot cavity . In the latter proposal photons in the nonlinear cavity are treated as massive particles weakly interacting with each other via some unspecified nonlinear optical mechanism.

Very recently strong evidences of nonlinear optical effects such as single- photon tunneling  and light-controlled photon tunneling  occurring at very low light intensities in natural nanometer-scale pinholes and artificial cylindrical nanochannels drilled using focused ion beam milling in thick gold films and membranes, and filled with nonlinear polymer have been observed. Transmission of nanoholes in gold and silver films is commonly accepted to happen via excitation of intermediary surface plasmons (SP) propagating along the metal film surfaces, and cylindrical surface plasmons (CSP) excited in the cylindrical nanochannels . In the case of long cylindrical nanochannels filled with nonlinear material theoretical description of nonlinear optical transmission of the nanochannels requires understanding of CSP behavior (such as mode propagation and interaction) at a level of single optical quanta in the presence of strong nonlinear optical interaction between them. So far to my knowledge, no such theory has been developed.

In this Letter I analyze CSP behavior in two closely related cylindrical geometries, such as a long cylindrical nanochannel drilled in a metal that supports surface plasmons (such as gold and silver) and filled with nonlinear material, and a long cylindrical metal nanowire embedded in a nonlinear material. In both geometries a specific case of chiral nonlinear media will be addressed (it is important to mention that metal itself becomes a chiral material in an external magnetic field , and all metal interfaces exhibit pronounced nonlinear behavior, so that chirality and nonlinear optical effects must be considered in cylindrical surface plasmon propagation). Since CSPs may be considered as if they live in a ”three-dimensional space-time” with two spatial dimensions being a very long z-dimension along the axis of the cylinder, and a small ”compactified” $`\varphi `$-dimension along the circumference of the cylinder, the theory of CSP mode propagation and interaction is formulated similar to the three-dimensional Kaluza-Klein theory . In the most interesting cylindrical nanowire case, higher ($`n>0`$) CSP modes are shown to posses a quantized effective charge (not of electric origin) proportional to their angular momentum. These slow moving effective charges exhibit strong long- range interaction via fast uncharged massless CSPs with zero angular momentum ($`n=0`$). The character of CSP mode interaction in a nanochannel is somewhat similar, but in this case the zero angular momentum modes acquire mass, so that interaction of higher mode quanta becomes short-ranged. Such a mode-coupling theory may be further extended to describe mode propagation and interaction in other chiral optical waveguides, as well as interaction of nonzero angular momentum electron states in thin gold wires and carbon nanotubes.

Let us recall the basic features of Kaluza-Klein theories. In modern Kaluza-Klein theories  the extra N-4 space-time dimensions are considered to be compact and small (with characteristic size on the order of the Planck length). The symmetries of this internal space are chosen to be the gauge symmetries of some gauge theory , so a unified theory would contain gravity together with the other observed fields. In the original form of the theory a five-dimensional space-time was introduced where the four dimensions $`x^1,\mathrm{},x^4`$ were identified with the observed space-time. The associated 10 components of the metric tensor $`g_{\alpha \beta }`$ were used to describe gravity. After a compactified fifth dimension $`x^5`$ with a small circumference $`L`$ was added, the extra four metric components $`g_{\alpha 5}`$ connecting $`x^5`$ to $`x^1,\mathrm{},x^4`$ gave four extra degrees of freedom, which were interpreted as the electromagnetic potential (here we use the following convention for greek and latin indices: $`\alpha =1,\mathrm{},4`$; $`i=1,\mathrm{},5`$). An additional scalar field $`g_{55}`$ or dilaton may be either set to a constant, or allowed to vary.

When a quantum field $`\psi `$ coupled to this metric via an equation

$$\mathrm{}_5\psi +a\psi =0$$
(1)
is considered, where $`\mathrm{}_5`$ is the covariant five-dimensional d’Alembert operator, the solutions for the field $`\psi `$ must be periodic in the $`x^5`$ coordinate. This leads to the appearance of an infinite ”tower” of solutions with quantized $`x^5`$-component of the momentum:

$$q_n^5=2\pi n/L$$
(2)
where $`n`$ is an integer. In our four-dimensional space-time on a large scale such solutions with $`n0`$ interact with the electromagnetic potential $`g_{\alpha 5}`$ as charged particles with an electric charge $`e_n`$ and mass $`m_n`$:

$$e_n=\mathrm{}q_n(16\pi G)^{1/2}/c$$
(3)
$$m_n=\mathrm{}(q_n^2a)^{1/2}/c$$
(4)
where $`G`$ is the gravitational constant (see for example the derivation in ). For the purposes of discussion below let us follow this derivation when an angular coordinate $`\varphi ^5`$ varying within an interval from 0 to $`2\pi `$ is introduced, so that

$$\varphi ^5=2\pi \frac{x^5}{L}$$
(5)
Now the metric can be written as

$$ds^2=g_{\alpha \beta }dx^\alpha dx^\beta +2g_{\alpha 5}dx^\alpha d\varphi ^5+g_{55}d\varphi ^5d\varphi ^5,$$
(6)
Here we are not interested in possible spatial dependence of $`g_{55}`$ and consider all $`g_{i5}`$ components to be independent of $`\varphi ^5`$. Equation (1) with $`a=0`$ for the quantum field $`\psi `$ in this metric should be written as

$`{\displaystyle \frac{}{x^\alpha }}(g^{\alpha \beta }{\displaystyle \frac{\psi }{x^\beta }})+{\displaystyle \frac{}{x^\alpha }}(g^{\alpha 5}{\displaystyle \frac{\psi }{\varphi ^5}})+`$ (7)
$`{\displaystyle \frac{}{\varphi ^5}}(g^{5\alpha }{\displaystyle \frac{\psi }{x^\alpha }})+{\displaystyle \frac{}{\varphi ^5}}(g^{55}{\displaystyle \frac{\psi }{\varphi ^5}})=0`$ (8)
Since we assume that the $`g_{i5}`$ do not depend explicitly on $`\varphi ^5`$, we should not address the terms like $`g^{i5}/\varphi ^5`$. Thus, we may search for the solutions in the usual form as $`\psi =\mathrm{\Psi }(x^\alpha )e^{iq\varphi ^5}`$, where periodicity in $`\varphi ^5`$ requires $`q_n=n`$. As a result, we obtain

$$\mathrm{}\psi q_n^2\frac{1g_{\alpha 5}g^{\alpha 5}}{g_{55}}\psi +2iq_ng^{\alpha 5}\frac{\psi }{x^\alpha }+iq_n\frac{g^{\alpha 5}}{x^\alpha }\psi =0$$
(9)
This is the same as the Klein-Gordon equation in the presence of an electromagnetic field: in four-dimensional space-time it describes a particle of mass

$$m=\frac{2\pi \mathrm{}q_n}{cg_{55}^{1/2}(q_n)},$$
(10)
which interacts with a vector field $`g^{\alpha 5}`$ through a quantized charge $`e_nq_n`$. In this theory the conservation of the electric charge is a simple consequence of the conservation of the $`x^5`$-component of the momentum. Similar Kaluza-Klein theory may be formulated in a three-dimensional space-time, which has one compactified spatial dimension. Again, quantized $`\varphi `$-component of the momentum will play a role of an effective charge, which interacts with a two-component vector field $`g^{\alpha 3}`$.

It is well-known that Maxwell equations in a general curved space-time background $`g_{ik}(x,t)`$ are equivalent to the macroscopic Maxwell equations in the presence of matter background with some nontrivial electric and magnetic permeability tensors $`ϵ_{ik}(x,t)`$ and $`\mu _{ik}(x,t)`$ . Thus, strong similarity between the results of the three-dimensional Kaluza-Klein theory and the solutions of Maxwell equations in some quasi-three-dimensional cylindrical waveguide geometries may be expected: in both cases we consider a massless field in a geometry which has a compactified cyclic $`\varphi `$-dimension in the presence of nontrivial space-time curvature or permeability tensors, respectively. This similarity is expected to be especially strong in the case of surface plasmon waveguides, since surface plasmons live in (almost) three- dimensional space-times on the metal interfaces.

Let us first consider general properties of Maxwell equations in a medium (media), which has cylindrically symmetric geometry. The insight gained from the previous discussion tells us that in order to look similar to the three- dimensional Kaluza-Klein theory, the medium should discriminate between left- and right- circular polarized waves (the waves which have opposite angular momenta, and thus expected to posses opposite effective Kaluza-Klein charges). This means that the medium should be chiral or optically active. There are different ways of introducing optical activity (gyration) tensor in the macroscopic Maxwell equations. It can be introduced in a symmetric form, which is sometimes called Condon relations :

$$\stackrel{}{D}=ϵ\stackrel{}{E}+\gamma \frac{\stackrel{}{B}}{t}$$
(11)
$$\stackrel{}{H}=\mu ^1\stackrel{}{B}+\gamma \frac{\stackrel{}{E}}{t}$$
(12)
Or it can be introduced only in an equation for $`\stackrel{}{D}`$ (see ). In our consideration I will follow Landau and Lifshitz , and for simplicity use only the following equation valid in isotropic or cubic-symmetry materials:

$$\stackrel{}{D}=ϵ\stackrel{}{E}+i\stackrel{}{E}\times \stackrel{}{g},$$
(13)
where $`\stackrel{}{g}`$ is called the gyration vector. If the medium exhibits magneto- optical effect and does not exhibit natural optical activity $`\stackrel{}{g}`$ is proportional to the magnetic field $`\stackrel{}{H}`$:

$$\stackrel{}{g}=f\stackrel{}{H},$$
(14)
where the constant $`f`$ may be either positive or negative. For metals in the Drude model at $`\omega >>eH/mc`$

$$f(\omega )=\frac{4\pi Ne^3}{cm^2\omega ^3}=\frac{e\omega _p^2}{mc\omega ^3},$$
(15)
where $`\omega _p`$ is the plasma frequency and $`m`$ is the electron mass .

Initially, let us consider a medium with $`\stackrel{}{g}=\stackrel{}{g}(r,z,t)`$ directed along the $`\varphi `$-coordinate. Such a distribution may be produced, for example, in a medium exhibiting magneto-optical effect around a cylindrical nanowire if a current is passed through the wire. After simple calculations we obtain a wave equation in the form:

$$\stackrel{}{}\times \stackrel{}{}\times \stackrel{}{B}=\mathrm{\Delta }\stackrel{}{B}=\frac{ϵ}{c^2}\frac{^2\stackrel{}{B}}{t^2}+\frac{i}{c}\frac{(\stackrel{}{}\times [\stackrel{}{E}\times \stackrel{}{g}])}{t}$$
(16)
The z-component of this wave equation for a solution $`e^{in\varphi }`$ is:

$`{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}(r{\displaystyle \frac{B_z}{r}})+{\displaystyle \frac{^2B_z}{z^2}}{\displaystyle \frac{ϵ^2B_z}{c^2t^2}}{\displaystyle \frac{n^2}{r^2}}B_z+`$ (17)
$`i{\displaystyle \frac{ng}{rc}}{\displaystyle \frac{(iE_z)}{t}}+{\displaystyle \frac{in}{cr}}({\displaystyle \frac{g}{t}})(iE_z)=0`$ (18)
It can be re-written in the form similar to equations (1) and (8) as follows:

$$\widehat{a_r}B_z+\mathrm{}_2B_z\frac{n^2}{g_{\varphi \varphi }}B_z+in(\frac{g}{r})(\frac{iE_z}{ct})+in\frac{(g/r)}{ct}E_z=0,$$
(19)
where $`\widehat{a_r}`$ plays the role of factor $`a`$ in equation (1), and $`g_{\varphi \varphi }=r^2`$. Similarity of equations (8) and (17) becomes more evident if we identify the Kaluza-Klein vector field $`g^{0\varphi }`$ as $`g/r`$, disregard terms higher than linear in $`g`$, and recall that for cylindrical surface plasmon modes $`iE_z=\alpha B_z`$, where the coefficient of proportionality $`\alpha `$ is real and is determined by the boundary conditions . The missing factor of 2 in front of the fourth term in equation (17) originates from the asymmetric way of introducing the gyration in the macroscopic Maxwell equations (see equations (10-12) and the relevant discussion). Thus, from the macroscopic Maxwell equations we arrive to a picture of $`n>0`$ waveguide modes interacting with the ”gyration potential” $`g/r`$ via quantized effective ”chiral charges” proportional to the angular momenta of the modes.

The dispersion law and electromagnetic field distribution of cylindrical surface plasmon modes of a cylindrical metal wire may be found in . The $`\omega (k)`$ of the $`n=0`$ CSP mode goes down to $`\omega =0`$ approaching the light line $`\omega =kc/ϵ^{1/2}`$ from the right, as $`k0`$. The field of this mode has only the following nonzero components: $`E_r`$, $`E_z`$, and $`H_\varphi `$, thus satisfying our initial requirements for the direction of $`\stackrel{}{g}`$, and if we recall equation (13), we notice that in the media exhibiting magneto- optical effect, such as the cylindrical metal wire itself, the higher $`n>0`$ CSP modes interact with the $`H_\varphi `$ field of the $`n=0`$ CSPs via their effective ”chiral charges” proportional to $`n`$. Thus, CSP quanta with $`n=0`$ may be considered as massless (in the $`k0`$ limit) quanta of the ”gyration field”. The dispersion laws of the CSP modes with $`n>0`$ start at some nonzero frequencies and intersect the light line (become nonradiative and de-couple from free-space photons) at some finite $`\omega `$. The modes of different $`n`$ are well separated from each other, and there is no crossing. In the $`k\mathrm{}`$ limit the $`\omega (k)`$ of all the CSP modes saturates at $`\omega =\omega _p/2^{1/2}`$. Thus, the group velocity $`d\omega /dk`$ of the higher modes goes to 0 in this limit: the ”charged” quanta are slow.

Let us now derive an analog of the Poisson equation for the ”gyration potential” and ”chiral charges” around the cylindrical metal wire. Let us search for the solutions of the wave equation (15) in the form $`\stackrel{}{B}=\stackrel{}{B_0}+\stackrel{}{B_n}`$ and $`\stackrel{}{E}=\stackrel{}{E_0}+\stackrel{}{E_n}`$, where $`\stackrel{}{B_0}`$ and $`\stackrel{}{B_n}`$ are the CSP fields with zero and nonzero ($`n`$) angular momenta, respectively, and the ”gyration field” in the magneto-active media is obtained in a self-consistent manner as $`\stackrel{}{g}=f(\stackrel{}{H}+\stackrel{}{B_0}+\stackrel{}{B_n})`$, where $`\stackrel{}{H}`$ is a constant external field. We are interested in the solution for the field $`\stackrel{}{B_0}`$ in the limit $`\omega _00`$ in the presence of the $`\stackrel{}{B_n}`$ field. The resulting nonlinear Maxwell equation may be simplified assuming that the field $`\stackrel{}{B_n}`$ is supposed to be the solution of linear Maxwell equation, and the terms proportional to $`f^2`$ and higher may be neglected. As a result, we obtain:

$`\mathrm{\Delta }\stackrel{}{B_0}={\displaystyle \frac{if}{c}}{\displaystyle \frac{(\stackrel{}{}\times [\stackrel{}{E_n}\times \stackrel{}{B_n}])}{t}}{\displaystyle \frac{if}{c}}{\displaystyle \frac{(\stackrel{}{}\times [\stackrel{}{E_0}\times \stackrel{}{B_n}])}{t}}`$ (20)
$`{\displaystyle \frac{if}{c}}{\displaystyle \frac{(\stackrel{}{}\times [\stackrel{}{E_n}\times \stackrel{}{B_0}])}{t}}`$ (21)
Since the fields $`\stackrel{}{B_0}`$ and $`\stackrel{}{B_n}`$ are not supposed to be coherent, their products disappear after time averaging, and we are left with the equivalent Poisson equation:

$$\mathrm{\Delta }\stackrel{}{B_0}=\frac{f\omega _n}{c}\stackrel{}{}\times [\stackrel{}{E_n}\times \stackrel{}{B_n}]=\frac{4\pi f\omega _n}{c^2}\stackrel{}{}\times \stackrel{}{S_n},$$
(22)
where $`\stackrel{}{S_n}`$ is the Pointing vector of the CSP field with the nonzero ($`n`$) angular momentum. As an interesting consequence of this equation, let us note that in an isotropic medium

$$\stackrel{}{}\times (\stackrel{}{}\times \stackrel{}{B_0}\frac{4\pi f\omega _n}{c^2}\stackrel{}{S_n})=0$$
(23)
This observation lets us to conclude what is the real physical nature of the ”chiral charge” in the normal language of four-dimensional physics. Due to magneto-optical effect, cylindrical surface plasmons with nonzero angular momenta behave as current loops bound to the metal wire. Unlike dipole-dipole interaction of the current loops in free space, the interaction of CSPs is mitigated by the presence of cylindrical surface plasmons with zero angular momentum. Thus, interaction of current loops becomes quasi one-dimensional and very long-ranged.

Using the same approximations as before, we can also derive an analog of the Gauss theorem for the ”chiral charges”. Let us consider a cylindrical volume $`V`$ around a cylindrical metal wire (see Fig.1), such that the side wall of the volume $`V`$ is located very far from the wire and the CSP fields are zero at this wall.

Using the vector calculus theorem we can write

$$_V\stackrel{}{}\times \stackrel{}{S}d^3x=_S\stackrel{}{n}\times \stackrel{}{S}𝑑a=_{S_2}\stackrel{}{N}\times \stackrel{}{S}𝑑a_{S_1}\stackrel{}{N}\times \stackrel{}{S}𝑑a,$$
(24)
where $`S`$ is the closed two-dimensional cylindrical surface bounding $`V`$, with area element $`da`$ and unit outward normal $`\stackrel{}{n}`$ at $`da`$, $`S_1`$ and $`S_2`$ are the front and the back surfaces of $`V`$, and $`\stackrel{}{N}`$ is the chosen direction of the wire. Using equation (20) we obtain

$`{\displaystyle _V}{\displaystyle \frac{4\pi f\omega _n}{c^2}}\stackrel{}{}\times \stackrel{}{S}d^3x={\displaystyle _{S_2}}\stackrel{}{N}\times [\stackrel{}{}\times \stackrel{}{B_0}]𝑑a`$ (25)
$`{\displaystyle _{S_1}}\stackrel{}{N}\times [\stackrel{}{}\times \stackrel{}{B_0}]𝑑a`$ (26)
Since $`\stackrel{}{N}\times [\stackrel{}{}\times \stackrel{}{B_0}]=B_{0\varphi }/z`$, we see that a ”chiral charge” produces a local step in the ”gyration field”. If this effective Gauss theorem would be precise at any distance, we would come up with an unphysical result that a localized ”chiral charge” would produce an infinite field $`B_\varphi `$ at infinity. In reality, one has to take into account all the effects that we neglected while deriving the effective Poisson equation and Gauss theorem, such as higher than quadratic terms in the nonlinear Maxwell equations, and also take into account final lifetime of the cylindrical surface plasmons. Nevertheless, the underlying Kaluza-Klein mechanism of the ”chiral charge” interaction lends credibility to the conclusion that CSPs of a cylindrical metal nanowire exhibit very strong long-range interaction. The case of cylindrical nanochannel is somewhat different, since in this case the zero angular momentum mode has nonzero frequency in the limit $`k0`$ (in other words it acquires some effective mass), so that interaction of higher mode quanta via exchange of zero- angular momentum CSPs becomes more short-ranged. Similar mode-coupling theory may be further developed to describe mode propagation and interaction in other chiral optical waveguides.

Let us estimate the strength of chiral charges interaction. Using equations (17) and (22) within the approximations described above the potential energy of two equal chiral charges $`n`$ separated by distance $`z`$ can be written as

$$W=z\frac{2\pi \mathrm{}\omega _nf^2nS_{n\varphi }}{ϵrc},$$
(27)
where $`S_{n\varphi }=n\omega _n(ϵE_z^2+B_z^2)/(4\pi \beta _n^2r)`$ is the $`\varphi `$-component of the Pointing vector, $`\beta _n^2=ϵ\omega _n^2/c^2k_n^2`$ is defined by the dispersion law of the $`n`$th mode, and $`r`$ is the radius of the cylindrical metal wire. Assuming $`f(\omega )`$ from (14) and $`\beta _nk_n`$ the order of magnitude estimate of $`W`$ may be written as

$$W=zn^2(\frac{e^2}{\mathrm{}c})\mathrm{}\omega _n(\frac{\omega _p^4}{\omega _n^4})(\frac{\lambda }{r})\frac{\mathrm{}^2}{r^3m^2c^2},$$
(28)
where $`\lambda `$ is the wavelength of $`\omega _n`$ light in vacuum. As may be expected, the nonlinear optical interaction of CSP modes grows inversely proportional to $`r^4`$, and achieves considerable strength in nanometer scale metallic nanowires and nanochannels even without addition of nonlinear dielectrics inside a nanochannel or around a nanowire. If we assume $`\mathrm{}\omega _n=1`$eV, $`\mathrm{}\omega _p=7`$eV, and $`r=1`$nm (disregarding nonzero skin depth) the potential energy of two unit chiral charges is of the order of $`W2\times 10^4z`$eV/cm, so even interaction of single CSP quanta separated by submicrometer distances is considerable.

It is also important to mention that electron states with nonzero angular momenta in metal nanowires (such as gold, silver, copper, aluminum and many other metals which support surface plasmons, as well as carbon nanotubes with metallic conductance) should also behave like chiral charges described above. Their long-range current-loop-like interaction via exchange of $`n=0`$ CSPs must be taken into account while considering mesoscopic conductance of thin wires. The strength of such interaction should be much stronger than the interaction of CSPs with each other, since the potential energy of two electron current-loops is proportional to $`f`$ (and not to $`f^2`$ as in equation (24)). Long range attractive interaction between two electron current loops may have important consequences for superconductive properties of metal nanowires and carbon nanotubes, with considerable strength of such interaction indicative of potential high $`T_c`$ superconductivity.


text:
# Constrained Opinion Dynamics: Freezing and Slow Evolution

## Abstract

We study opinion formation in a population of leftists, centrists, and rightist. In an interaction between neighboring agents, a centrist and a leftist can become both centrists or leftists (and similarly for a centrist and a rightist), while leftists and rightists do not affect each other. For any spatial dimension the final state is either consensus (of one of three possible opinions), or a frozen population of leftists and rightists. In one dimension, the opinion evolution is mapped onto a constrained spin-1 Ising model with zero-temperature Glauber kinetics. The approach to the final state is governed by a $`t^\psi `$ long-time tail, with $`\psi `$ a non-universal exponent that depends on the initial densities. In the frozen state, the length distribution of single-opinion domains has an algebraic small-size tail $`x^{2(1\psi )}`$ with average domain length $`L^{2\psi }`$, where $`L`$ is the length of the system.

One of the basic issues in opinion dynamics is to understand the conditions under which consensus or diversity is reached from an initial population of individuals (agents) with different opinions. Models for such evolution are typically based on each agent freely adopting a new state in response to opinions in a local neighborhood F74 . The attribute of incompatibility – in which agents with sufficiently disparate opinions do not interact – has recently been found to prevent ultimate consensus from being reached LN97 ; opinion . Related phenomenology arises in the Axelrod model axelrod ; CMV , a simple model for the formation and evolution of cultural domains. The goal of the present paper is to investigate the role of incompatibility within a minimal model for opinion dynamics. This constraint has a profound effect on the nature of the final state. Moreover, there is anomalously slow relaxation to the final state. While we primarily frame our discussion in terms of opinion dynamics, our results apply equally well to the coarsening of spin systems. In the latter context, we obtain a new non-universal kinetic exponent in one dimension that originates from topological constraints on the arrangement of spins.

We consider a ternary system in which each agent can adopt the opinions of leftist, centrist, and rightist. The agents populate a lattice and in a single microscopic event an agent adopts the opinion of a randomly-chosen neighbor, but with the crucial proviso that that leftists and rightists are considered to be so incompatible that they do not interact. While a leftist cannot directly become a rightist (and vice versa) in a single step, the indirect evolution leftist $``$ centrist $``$ rightist is possible. Our model is similar to the classical voter model voter  and also turns out to be isomorphic to the 2-trait 2-state Axelrod model axelrod ; CMV . Due to the incompatibility constraint in our model, the final opinion outcome can be either consensus or a frozen mixture of extremists with no centrists. Figure 1 shows a typical frozen state on the square lattice (with periodic boundary conditions). Notice the nested enclaves of opposite opinions and the clearly visible clustering.

We can exploit the connection between the voter model and our opinion dynamics model to infer the final state of the latter. If we temporarily disregard the difference between leftists and rightists, the resulting binary system of centrists and extremists reduces to the voter model, for which one of two absorbing states — either all centrists or all extremists — is eventually reached. In the context of the ternary opinion system, the latter event can mean either a consensus of extremists or a frozen mixed state of leftists and rightists, as depicted in Fig. 1.

Because of the underlying voter model dynamics, the average density of each species is globally conserved in any spatial dimension. Therefore $`\rho _i(t)=\rho _i(t=0)`$, where $`i`$ refers to one of the states $`(+,0,)`$ and the angle brackets denote an average over all realizations of the dynamics and over all initial states that are compatible with the specified densities. As a result of this conservation law, with probability $`P_0=\rho _0`$ the final state consists of all centrists and with probability $`1\rho _0`$ there are no centrists in the final state. In the latter case, there can be either a consensus of $`+`$ (this occurs with probability $`P_+`$), consensus of $``$ (probability $`P_{}`$), or a frozen mixed state (probability $`P_+`$). Figure 2 shows the dependence of these final state probabilities on $`\rho _0`$ in the mean-field limit (where all agents are interconnected) in the symmetric case $`\rho _+=\rho _{}=(1\rho _0)/2`$. In one and two dimensions, the final state probabilities are nearly identical to the mean-field predictions when $`\rho _+=\rho _{}`$, but differences become apparent in the strongly asymmetric cases of $`\rho _+\rho _{}`$. The final state probabilities also depend very weakly on the system size.

We now focus on the one dimensional case. Here our opinion dynamics model is equivalent to a constrained spin-1 Ising chain that is endowed with single-spin flip zero-temperature Glauber kinetics glauber , with leftist, centrist, and rightist opinions equivalent to the spin states $``$, 0, and $`+`$, respectively. The incompatibility constraint means that neighboring $`+`$ and $``$ spins do not interact.

This Ising model picture suggests that the best way to analyze the dynamics in one dimension is to reformulate the system in terms of domain walls. There are three types of domain walls: freely diffusing mobile domain walls between $`+0`$ and between $`0`$, denoted by $`M_+`$ and $`M_{}`$, respectively, and stationary domain walls $`S`$ between $`+`$. The mobile walls evolve by

$`M_\pm +M_\pm \mathrm{},M_\pm +M_{}S.`$ (1)
When a mobile wall hits a stationary wall, the former changes its sign while the latter is eliminated via the reaction

$`M_\pm +SM_{}.`$ (2)
Thus stationary domain walls are dynamically invisible; their only effect is that the sign of a mobile wall changes whenever it meets a stationary wall, after which the latter disappears (Fig. 3). The inertness of the stationary walls is reminiscent of kinetic constraints in models of glassy relaxation. These constraints typically lead to extremely slow kinetics JE ; SE ; MDG ; CRRS , as is also observed in our opinion dynamics model.

The rate equations corresponding to the processes in Eqs. (1) and (2) are

$$\dot{M}=2M^2\dot{S}=MS+M^2.$$
(3)
These give the asymptotic behaviors $`Mt^1`$ and $`St^{1/2}`$. Thus an approximate rate equation approach already predicts that stationary domain walls decay slower than mobile walls.

An important subtlety in the arrangement of domain walls is that an arbitrary initial opinion state necessarily leads to an even number of mobile walls between each pair of stationary walls. It is also easy to verify that domain wall sequences of the form $`\mathrm{}M_+M_{}M_+\mathrm{}`$ cannot arise from an underlying opinion state. These topological constraints play a crucial role in the kinetics.

The exact density of mobile walls can be obtained by again mapping the constrained spin-1 system onto a spin-1/2 system that is equivalent to the voter model. In this mapping, we consider both $`+`$ and $``$ spins as comprising the same (non-zero) spin state, while the zero spins comprise the other state. With this identification, the reduced model is just the spin-1/2 ferromagnetic Ising chain with zero-temperature Glauber dynamics and no kinetic constraint. In this reduced model, domain walls $`M_+`$ and $`M_{}`$ are indistinguishable and they diffuse and annihilate when upon colliding. The density of mobile walls $`M(t)=M_+(t)+M_{}(t)`$ is known exactly for arbitrary initial conditions from the original Glauber solution glauber . For initially uncorrelated opinions and if the magnetization of the spin-1/2 system – here the difference between the density of non-zero and zero spins – equals $`m_0`$, then amar

$$M(t)=\frac{1m_0^2}{2}e^{2t}\left[I_0(2t)+I_1(2t)\right],$$
(4)
where $`I_k`$ is the modified Bessel function of index $`k`$.

In the spin-1 system, $`m_0=\rho _++\rho _{}\rho _0`$, or $`m_0=12\rho _0`$ by normalization. If the initial densities of $`+`$ and $``$ opinions are equal, then $`M_+(t)=M_{}(t)`$ and their densities are

$`M_\pm (t)`$ $`=`$ $`\rho _0(1\rho _0)e^{2t}\left[I_0(2t)+I_1(2t)\right]`$ (5)
$``$ $`\rho _0(1\rho _0)(\pi t)^{1/2}.`$
As expected, the mobile wall density asymptotically decays as $`t^{1/2}`$ because of the underlying diffusive dynamics. However, we find numerically that the density of stationary domain walls $`S(t)`$ decays as

$$S(t)t^{\psi (\rho _0)},$$
(6)
with a non-universal exponent $`\psi (\rho _0)`$ that goes to zero as $`\rho _00`$ (Fig. 4).

To help understand the mechanism for the slow decay of the stationary domain wall density, we also simulated a test system with spatially uncorrelated domain walls. While such a domain wall state cannot arise from any initial set of opinions, we can prepare directly an uncorrelated arrangement of domain walls with prescribed densities. For any initial condition in this test system, the stationary wall density decays as $`t^{3/8}`$ (inset to Fig. 4), consistent with known results on persistence per ; BKR ; stauffer . Here persistence refers to the probability that a given lattice site is not hit by any diffusing domain wall. For the kinetic spin-1/2 Ising model, the persistence probability decays as $`t^\theta `$, with $`\theta =3/8`$ der , independent of the initial domain wall density, when the walls are initially uncorrelated. Thus the topological constraints imposed on the domain wall arrangement by the initial opinion state appear to control the dynamics.

These topological constraints lead to the initial-condition dependence of the amplitude in the mobile wall density (Eq. (5)). This arises because for $`\rho _00`$ the system initially consists of long strings of stationary walls that are interspersed by pairs of mobile walls, and their survival probability is proportional to their initial (unit) separation redner , leading to the asymptotic density for mobile walls is $`M2\rho _0/\sqrt{\pi t}`$ (Eq. (5)). We now exploit this observation to estimate the density of stationary walls as $`\rho _00`$. Within a rate-equation approximation, the density of stationary domain walls decays according to

$$\dot{S}=kMS.$$
(7)
While such an equation is generally inapplicable in low spatial dimension, we can adapt it to one dimension by employing an effective time-dependent reaction rate $`k\sqrt{2/\pi t}`$ BKR ; redner . This is just the time-dependent flux to an absorbing point due to a uniform initial background of diffusing particles; such a rate phenomenologically accounts for effects of spatial fluctuations in one dimension. Substituting the asymptotic expression for $`M(t)`$ from Eq. (4) and the reaction rate $`k\sqrt{2/\pi t}`$ into this rate equation, we find that the density of stationary walls decays as $`t^\psi `$ with $`\psi (\rho _0)=\sqrt{8}\rho _0/\pi `$ as $`\rho _00`$. It is the amplitude in the density of mobile domain walls that ultimately causes the slow decay in the stationary wall density.

A more compelling way to determine $`\psi (\rho _0)`$ is via persistence in the $`q`$-state Potts model. Because the initial magnetization in the reduced spin-1/2 system is $`m_0=12\rho _0`$, it has been argued (see e.g. sire ) that this system should be identified with the $`q`$-state Potts model with $`m_0=2/q1`$, or $`q=(1\rho _0)^1`$. Using the exact persistence exponent for the $`q`$-state Potts model with Glauber kinetics der , and identifying $`\psi `$ with this persistence exponent, we obtain

$$\psi (\rho _0)=\frac{1}{8}+\frac{2}{\pi ^2}\left[\mathrm{cos}^1\left(\frac{12\rho _0}{\sqrt{2}}\right)\right]^2,$$
(8)
with the limiting behavior $`\psi (\rho _0)2\rho _0/\pi `$ as $`\rho _00`$. This asymptotics agrees with our numerical results for $`\rho _00.4`$ (Fig. 5) but then deviates for larger $`\rho _0`$, where $`\psi (\rho _0)`$ must monotonically approach 1/2 as $`\rho _01`$. It should be noted that in this identification with persistence, we have ignored the creation of stationary interfaces due to the meeting of mobile domain walls. This creation process occurs with a rate $`(dS/dt)_{\mathrm{gain}}t^{3/2}`$ and is subdominant for $`\psi <1/2`$ FKB .

An important characteristic of the system is the spatial arrangement of domain walls. From our simulations, the mean distances between nearest-neighbor $`MM`$ and $`MS`$ walls both appear to grow as $`t^{1/2}`$ due to the diffusive motion of mobile domain walls. The distributions of these two distances both obey scaling, with scaling function of the form $`ze^{z^2}`$, where $`z=x(t)/x(t)`$ is the scaled separation between walls. In contrast, the distances between neighboring stationary walls $`x_{SS}`$ appear to be characterized by two length scales. There are large gaps of length of the order of $`t^{1/2}`$ that are cleared out by mobile walls before they annihilate, but there are also many very short distances remaining from the initial state (Fig. 3). The corresponding moments $`x_{SS}^k(t)^{1/k}`$ reflect this multiplicity of scales, with $`x_{SS}^k(t)^{1/k}`$ approaching a $`t^{1/2}`$ growth law as $`k\mathrm{}`$, and growing extremely slowly in time for $`k0`$.

Finally, we quantify the frozen final state by the magnetization distribution $`P(m)`$, namely, the density difference between $`+`$ and $``$ spins. This distribution becomes broader as $`\rho _0`$ increases (Fig. 6), reflecting the fact that there is progressively more evolution before the system ultimately freezes. For small $`\rho _0`$, $`P(m)`$ has a $`m^2`$ tail. We may explain this result by considering the evolution of a single pair of mobile walls. This pair annihilates at time $`t`$ with probability density $`\mathrm{\Pi }(t)t^{3/2}`$. The total magnetization of the resulting frozen state scales as $`t^{1/2}`$ since the domain wall pair annihilates at a distance $`xt^{1/2}`$ from its starting point. Then from $`P(m)dm=\mathrm{\Pi }(t)dt`$, together with $`\mathrm{\Pi }(t)t^{3/2}`$ and $`m\rho _0t^{1/2}`$, we obtain $`P(m)\rho _0m^2`$. While the argument has been formulated in one dimension, we expect this result to apply in all spatial dimension.

The frozen state is reached when $`t=T_fL^2`$; this is the time needed for mobile domain walls to diffuse throughout the system and thus be eliminated. At this time, the density of stationary domain walls is of the order of $`ST_f^\psi L^{2\psi }`$. Thus the average length of single-opinion domains is $`xL^{2\psi }`$. Numerically we find that the domain length distribution has a power law tail, $`F(x)x^\mu `$, with $`1<\mu <2`$. The lower bound ensures normalizability while the upper bound implies that $`x`$ diverges as $`L\mathrm{}`$. From the above power law form, $`x=𝑑xxF(x)L^{2\mu }`$. This matches our previous estimate of $`xL^{2\psi }`$ when $`\mu =2(1\psi )`$. Figure 7 shows the length distribution of single-opinion domains in the frozen state. Direct estimates of the exponent $`\mu `$ from this plot are in good agreement with the exponent relation $`\mu =2(1\psi )`$, with $`\psi `$ obtained from the time dependence of the mobile wall density in Fig. 4.

In summary, the constraint that extremists with opposite opinions cannot influence each other substantially slows down opinion dynamics. In one dimension, the density of stationary interfaces between neighboring $`+`$ and $``$ spins decays as $`t^\psi `$, with $`\psi (\rho _0)2\rho _0/\pi `$ as $`\rho _00`$. This slow dynamics is a consequence of the subtle topological constraints on the arrangement of domain walls, together with the incompatibility constraint that neighboring $`+`$ and $``$ spins do not interact.

The final opinion outcome depends non-trivially on the initial densities of the leftist, rightist and centrist opinion states. With probability $`\rho _0`$ the final population consists of only centrists, while with probability $`1\rho _0`$ the final population does not contain any centrists. Analytically, we only know that consensus of centrists is reached with probability $`P_0=\rho _0`$; the determination of the complementary final state probabilities for general initial conditions $`P_+`$, $`P_{}`$, and $`P_+`$ remains an open challenge. Finally, a frozen mixture final state would not exist if there is a non-zero rate (however small) for opposite extremists to influence each other directly. Perhaps this lack of direct influence can account for the sad phenomenon of the proximity of incompatible populations in too many locations around the world.

We thank S. Majumdar for helpful correspondence. We are also grateful to NSF grant DMR9978902 for partial support of this work.