text:
# Thermal Resonance in Signal Transmission

## I Introduction

In the past few years it has become abundantly clear that the presence of noise in nonlinear systems may lead to an enhancement of a number of often desirable features such as energy localization and mobility and the detection and propagation of weak signals. The interplay of stochasticity and nonlinearity that amplifies the system response is a cooperative phenomenon whose detailed nature depends on the particular structure of the system and the forces acting upon it. One manifestation of the interplay is found in the phenomenon called stochastic resonance, which has been invoked in a wide range of physical , chemical , geological , and biological  systems. Recent literature, including our own work,  has focused on spatially extended systems  including noise-enhanced propagation in coupled arrays of bistable units , excitable media , reaction-diffusion systems , and dynamics and signal propagation in cardiac tissue . It has been repeatedly noted that discrete extended systems pose particular mathematical challenges that have barely been explored in spite of the fact that many physical systems are intrinsically discrete .

The ubiquitous picture of stochastic resonance involves a particle moving in a double-well potential subject to a weak external signal that periodically changes the potential by alternately raising and lowering the wells . The signal is “weak” if the periodic force is too small to cause the particle to scale the barrier between the wells. Nevertheless, an appropriate random force is sufficient to cause the particle to cross over the barrier even in the absence of a deterministic signal. In the simultaneous presence of a weak signal and a sufficiently weak noise, the transitions over the barrier occur rarely and at a rate determined by the noise intensity. These transitions are slow compared to the frequency of the deterministic signal; the transition rate then carries little information about the signal. At the other extreme, when the noise is strong it induces rapid transitions that are again essentially uninfluenced by the frequency of the signal. At an optimal noise intensity, however, the mean first passage time associated with the noise and the frequency of the signal are in synchrony (stochastic resonance), and the passage from one well to the other carries maximal information about the signal frequency.

A less ubiquitous but nonetheless important occurrence of stochastic resonance (that has been called “nonconventional” by its discoverers) arises for particles moving in nonlinear monostable potentials . It is argued that stochastic resonance can be expected to occur in any single-well underdamped system for which the spectral density of the fluctuations of the system in the absence of a periodic signal exhibits a well-resolved narrow peak that grows faster than quadratically with temperature. The effect is confirmed via analog simulations of a single-well Duffing oscillator  and of a SQUID loop . More recently, stochastic resonance at higher harmonics in monostable systems was ascertained for an overdamped system when the nonlinearity is not concentrated at the equilibrium position .

Recent developments in the field have generalized these ideas to linearly coupled arrays of bistable oscillators . A signal with the help of the noise in these arrays can cause a “phase jump.” If the noise is sufficiently weak, the phase jump travels in the form of a moving kink (strong noise causes random phase jumps that make it difficult to separately identify a phase jump associated with the signal). The creation or destruction of a phase kink is an activated process, i.e., the signal and/or noise must be sufficiently strong to cause a transition from one well of one of the bistable potentials to the other. The presence of such kinks is associated with an “energy gap:” it takes a finite amount of energy to destroy a kink. The language used in this description was originally borrowed from the kink soliton context.

Reported instances of stochastic enhancement and stochastic resonance in extended arrays involve coupled (overdamped) bistable units. Herein we show that enhanced propagation can be achieved through thermal tuning of even simpler discrete arrays of masses connected by monostable anharmonic springs (with no local potentials). The signal is identified with an amplitude that exceeds (by a predetermined amount) that due to the thermal background. Here there is no activation process and no energy gap, and the signal can simply disperse or dissipate. The language appropriate to this case is akin to that originally associated with envelope solitons. We focus on the propagation distance and amplitude along the chain of a signal continuously applied at one site of a 1D array, and compare results for harmonic and hard anharmonic chains. In particular, we show that when the anharmonic chain is immersed in a thermal bath, it is possible to maximize the distance of propagation and the amplitude of the signal at a given site by tuning the temperature to particular optimal values. We call this phenomenon thermal resonance. Our systems are in general not overdamped and thus include inertial contributions to the motions of the masses. Noise and damping represent a realistic thermal environment with a tunable temperature and dissipation that obey an appropriate fluctuation-dissipation relation. We have found and reported elsewhere  that the propagation of an energy pulse in a hard anharmonic array can be enhanced by immersion in a thermal bath, and that hard anharmonicity in the springs causes a tight and persistent packing of the energy. Those results suggest the possibility of a thermal resonance in the transport of a sustained external signal in these simple arrays.

In Section II we present our model and some details of the numerical integration of the equations of motion. Our characterization of a thermal resonance is presented in Section III, and our main results are shown in Section IV. In Section V we discuss the dependence of our results on different parameter models. Section VI contains our summary and conclusions.

## II The Model and Numerical Procedure

Our model consists of a one-dimensional chain of $`N`$ unit-mass sites, each connected to its nearest neighbors by either harmonic (quadratic) or hard anharmonic (quartic) springs. Accordingly, the Hamiltonian of the array is

$$H=\underset{n=0}{\overset{N1}{}}\left(\frac{p_n^2}{2}+\frac{k}{2}(x_nx_{n1})^2+\frac{k^{}}{4}(x_nx_{n1})^4\right),$$
(1)
where $`k`$ and $`k^{}`$ are the harmonic and anharmonic coupling constants respectively. Thermalization of the chains is achieved through a Langevin prescription for coupling the system to a heat bath. The stochastic equations of motion for sites $`n=1,\mathrm{},N1`$ are obtained from the Hamiltonian augmented by the usual Langevin forces:

$`\ddot{x}_n`$ $`=`$ $`k(x_{n1}x_n)k(x_nx_{n+1})+k^{}(x_{n1}x_n)^3`$ (3)
$`k^{}(x_nx_{n+1})^3\gamma \dot{x}_n+f_n(t),`$
where a dot represents a derivative with respect to time. The $`f_n(t)`$ are zero-centered, Gaussian, $`\delta `$-correlated fluctuations that satisfy the fluctuation-dissipation relation at temperature $`T`$,

$$f_n(t)f_n^{}(t^{})=2\gamma k_BT\delta _{n,n^{}}\delta (tt^{})$$
(4)
($`k_B`$ is Boltzmann’s constant). We impose periodic boundary conditions, so that $`x_Nx_0`$. A sustained signal is applied to the site $`n=0`$ that determines its velocity at all times:

$$\dot{x}_0=A\mathrm{sin}(\omega _0t).$$
(5)
The positions and momenta of all the other sites are otherwise “free” and determined by the equations of motion. We study the propagation of this signal along the chain as a function of the temperature. This particular way of applying a signal is of course not unique (e.g. one might apply an oscillatory force instead), but we have ascertained that the results are insensitive to the detailed choice.

An analytic solution of this problem is not available for an anharmonic chain, so we must rely on numerical integration, which is performed using the second order Heun’s method (equivalent to a second order Runge-Kutta integration) . The time step is determined by the period of oscillation of the velocity of the first site, $`\tau =2\pi /\omega _0`$, through the relation $`\mathrm{\Delta }t=\tau /2^{12}`$. For each simulation, the system is initially allowed to relax to thermal equilibrium. For all the simulations presented, this is achieved in less than $`20`$ units of the dimensionless time. Typically, after a transient that is longer the further the site is from site 0 (and thus a measure of the velocity of propagation), each site settles into stationary behavior that is a mixture of thermal motion and response to the signal. Sites that are far from the signal never exhibit this transient (thus indicating a finite distance of propagation) and simply continue their thermal motion. At any given site that is reached by the signal, one can observe the amplitude of the motion associated with the signal over and above the thermal motion.

Our interest here lies in demonstrating resonances in the dependence on temperature of the propagation distance and velocity and of the response amplitudes once stationarity has been achieved (in all our simulations, each site has settled into its long-time behavior after 100 first-site oscillations). Possible energy return effects around the periodic chain are prevented by making the chain sufficiently long and/or sufficiently increasing the dissipation parameter of distant sites. Our chains typically consist of 70 sites with a large dissipation at sites 27–32 (these numbers can easily be varied). Our “measurements” are then taken over 80 oscillation periods and non-zero-temperature results are averaged over 500 realizations. All of these choices (equilibration time before applying the signal, integration time step, transients, length of chain, and number of realizations) have been carefully tested.

## III Characterization of Thermal Resonance

We must choose a sensible response variable to characterize the behavior of our array. When stochastic resonance is studied in arrays of bistable potentials, the system response is usually analyzed in terms of a crossover time series that characterizes the transitions of each bistable element from one well to the other. In excitable media a reasonable response variable involves firing times of the individual elements. In our system the most convenient choice is the velocity of each site because it corresponds directly to the applied signal, and because the time average of the velocity in the stationary state vanishes at any temperature.

Fig. 1 is a dramatic but typical demonstration of a realization of thermal resonance. It shows the mean velocity $`v`$ (averaged over realizations) as a function of time for the $`12^{th}`$ site of a hard anharmonic chain at different temperatures. The first panel shows the results for zero temperature and for a very low temperature; the temperature increases in subsequent panels. At zero temperature the $`12^{th}`$ site hardly moves because the signal has been dissipated to the bath before reaching this site (in the corresponding harmonic chain the signal reaches the $`12^{th}`$ site quite vigorously at zero temperature, a confirmation of the fact that a given dissipation is much more effective in a hard anharmonic potential than in a harmonic one  – see Fig. 2). A very small temperature increase (still first panel) causes a large enhancement of the signal, which clearly now reaches the site. This is apparent in the oscillatory behavior of the velocity over and above the noisy background. The temperature in the second panel is close to its optimal value, that is, the value that most enhances the signal at this particular site relative to the thermal background. Hence the motion of the $`12^{th}`$ site at this temperature is mostly driven by the periodic forcing of the first site. A further increase in the temperature (third panel) causes the average velocity to become increasingly noisy because ever larger fluctuations dominate the dynamics. Finally, at a sufficiently high temperature (fourth panel), the signal is essentially buried in the fluctuations, and the motion is simply that imposed by the thermal bath. In contrast, in a harmonic array ($`k=3`$) the signal at any site simply degrades with increasing temperature. This is illustrated in Fig. 2.

To provide a quantitative measure of the thermal resonance, we define the power spectral density $`S_j(\omega )`$ at each site $`j`$ as

$$S_j(\omega )=_{\mathrm{}}^+\mathrm{}e^{i\omega \tau }v_j(t)v_j(t+\tau )𝑑\tau ,$$
(6)
where the brackets denote an ensemble average over realizations and an average over time. In addition to the thermal fluctuations, this function contains the spectral information about that part of the signal that has reached site $`j`$. An example of a portion of the spectrum for the case that we will call our “standard case” ($`k^{}=5`$, $`\gamma =0.2`$, $`\omega _0=1.0`$, $`A=0.5`$) is shown in Fig. 3. The signal extraction from background noise that characterizes stochastic resonance is traditionally performed via a signal-to-noise ratio ($`SNR`$:

$$SNR(j)\mathrm{log}_{10}\left(\frac{\mathrm{signal}\mathrm{power}(j)\times \mathrm{\Delta }\omega }{\mathrm{thermal}\mathrm{power}(j)}\right),$$
(7)
where the signal power is the value $`S_j(\omega _0)S_{j,noise}(\omega _0)`$, the thermal power $`S_{j,noise}(\omega _0)`$ is estimated by performing a $`4^{th}`$-order polynomial fit to $`S_j(\omega )`$ around – but not including – the forcing frequency $`\omega _0`$, and $`\mathrm{\Delta }\omega `$ denotes the frequency integration step and is equal to 0.0125 (the inverse of the 80 oscillation periods used as our measurement time) throughout this paper. This definition of the $`SNR`$ is not unique, but our results are robust with respect to variations in this definition.

## IV Thermal Resonances

We first present $`SNR`$ results for the harmonic chain, so as to clarify later the ways in which the anharmonic chain behaves differently. The analytic calculations associated with the harmonic chain are presented in Appendix A. These results serve as a test for our numerical simulations. The first panel in Fig. 4 shows $`SNR`$ curves as a function of temperature for different sites, and the second panel shows the same results as a function of distance from the forced site for different temperatures. The results are exactly as shown in the Appendix and as one would expect: the $`SNR`$ decreases monotonically with increasing temperature and with increasing distance from the applied signal. Note that $`SNR(j)`$ decreases with increasing temperature because the numerator in Eq. (7) is essentially independent of temperature while the denominator increases (see Appendix A). A point to note is that the decay of $`SNR(j)`$ with $`j`$ at a given temperature provides a measure of the shape of the stationary front of the signal at that temperature.

We now turn to the anharmonic chain. The first panel in Fig. 5 shows $`SNR`$ curves as a function of temperature for different sites. The thermal resonance is identified with the $`j`$-dependent maximum of $`SNR(j)`$ as a function of temperature. The signal at the first few sites is monotonically weakened as a function of the temperature (as it is at all sites in the harmonic chain), but the $`SNR`$ is enhanced with increasing temperature for sites farther away until it reaches a maximum (thermal resonance), beyond which it decreases. Note that the optimal temperature increases with distance from the first site. For the particular parameters used in this simulation the thermal resonance occurs most clearly at temperatures in the range $`k_BT<4`$ for sites between the $`8^{th}`$ and the $`15^{th}`$; these details can of course be modified by changing the parameters (see Section V).

Complementary results for $`SNR`$ curves as a function of the site for different temperatures are shown in the second panel in Fig. 5. The resonance effects are evidenced by the crossings of the different curves. The crossings reflect the rise and subsequent drop in the $`SNR`$ at a given site, and the fact that the curves cross at different sites confirms that the thermal resonance temperature varies from site to site.

The non-monotonic behavior of the $`SNR`$ with temperature is due to the fact that while the thermal power increases with $`T`$, the signal power at first also increases, and more rapidly than the thermal power. This is exactly the condition for the existence of stochastic resonance formulated by Dykman et al.  for single potentials, but that theory has not yet been generalized to extended systems. The fast increase of the signal power is consistent with trends found earlier for transmission of pulses along discrete arrays . We found that a pulse in a hard anharmonic chain travels more rapidly with increasing temperature (except for the lowest temperatures, see below), while in a harmonic array a change in temperature has no effect on pulse speed. While a pulse broadens with increasing temperature in both arrays, the pulse in the anharmonic array remains relatively more compact. At sufficiently high temperatures the signal power at a given site becomes independent of $`T`$ (as in the harmonic chain) because the signal response has reached its maximum value; a further increase in temperature only affects the signal further down the chain. The “crossover” temperature of the signal power from the increasing to the saturated behavior depends on the chain and signal parameters, a dependence explored in the next section.

The fact that the $`SNR`$ for the sites shown is very low (essentially zero) for the lowest temperatures at sites beyond the first few shows that for a purely anharmonic chain the signal essentially stops beyond the first few sites, whereas in the harmonic example the signal reaches all the sites shown even at the lowest temperatures shown. This is also consistent with the behavior shown earlier for the transmission of a pulse : the pulse velocity at sufficiently low temperatures is actually lower in the hard chain than in the harmonic case. Had we included a harmonic potential contribution equal to that of the harmonic chain, the $`SNR`$ vs $`k_BT`$ curves would start at the same values as in Fig. 4, but for sufficiently distant sites from the first they would still be non-monotonic. We have omitted a harmonic contribution to present the thermal resonance effect in its purest form.

The apparently monotonic behavior of the first few sites is due to the fact that for the temperatures shown the signal reaches these sites in any case. A resonance at these sites would be seen for different parameter values and/or at even lower temperatures.

Our preceding descriptions point to another interesting measure of a thermal resonance, namely, the propagation length $`\mathrm{\Lambda }`$, defined as the number of sites (i.e., distance along the chain) for which the $`SNR`$ exceeds a certain threshold value. Fig. 6 presents the temperature dependence of $`\mathrm{\Lambda }`$ for harmonic and anharmonic arrays and an arbitrarily chosen $`SNR`$-threshold value of $`1.2`$. Since the harmonic array does not exhibit thermal resonance, a monotonic decay of $`\mathrm{\Lambda }`$ with increasing temperature is observed. On the other hand, the hard chain shows a maximum for a moderate temperature. Again, the particular values of optimal temperature ($`k_BT3`$) and optimal distance ($`\mathrm{\Lambda }17`$) can be modified by choosing different parameters and/or different $`SNR`$ thresholds, but the qualitative behavior persists as seen in Fig. 6.

## V Parameter Dependences

The parameters that can be varied in our model are the amplitude $`A`$ and frequency $`\omega _0`$ of the velocity of the first site, the coupling parameter $`k^{}`$, and the damping coefficient $`\gamma `$.

The first panel of Fig. 7 shows the $`SNR`$ curves as a function of temperature for a lower amplitude than in the standard case. Comparison with the first panel of Fig. 5 shows that the overall $`SNR`$ is now (of course) lower and that the resonance temperature at each site has increased. The latter behavior indicates that for lower amplitudes the crossover temperature from an increasing signal to a saturated signal increases as the signal weakens. The second panel shows results for weaker coupling. Again the overall $`SNR`$ is lower and the resonance temperature higher at each site. A weaker coupling thus has the effect of weakening the signal.

The first panel of Fig. 8 shows the $`SNR`$ curves as a function of temperature for higher damping. Comparison with the first panel of Fig. 5 shows that the overall $`SNR`$ is again lower than in the standard case and the resonance temperature higher. The same trends are observed with a higher driving frequency as shown in the second panel of Fig. 8. Again, each of these changes leads to an effectively weakened signal.

Figures 7 and 8 illustrate clear trends. Decreasing $`A`$ or $`k^{}`$, or increasing $`\gamma `$ or $`\omega _0`$, lead to the following consequences: a) the $`SNR`$ at any given site and temperature decreases; b) the $`SNR`$ resonance temperature increases at any given site; c) at a given temperature the resonance occurs at a site closer to the first. These trends are consistent with those found earlier for transmission of pulses along similar arrays . We found that a pulse in a hard anharmonic chain travels more rapidly and relatively more compactly with increasing amplitude, decreasing damping, increasing coupling, or increasing driving frequency. By contrast, in a harmonic array changes in temperature or in the signal and damping parameters have no effect on pulse speed.

## VI Conclusion

We have demonstrated new thermal resonances in simple one-dimensional arrays of masses connected by hard springs. Our chains are not overdamped and hence include inertial effects. We have shown that the distance and amplitude of propagation of a signal imposed at one end of the chain can each be optimized by tuning the temperature of the system. The resonance behavior reflects the temperature dependence of the distance traveled by the front of a signal and the dispersion of this front once steady state has been reached. At a fixed temperature, any parameter change that leads to an increase in velocity of propagation tends to increase the range of sites where the signal can be detected as well as the $`SNR`$ at a given site. Increased velocity of propagation is also associated with a lowering of the resonance temperature for a given site.

A complete understanding of the behavior of these simple anharmonic arrays in a thermal environment requires and relies on a number of other inquiries, some of which we have undertaken. One concerns the distribution of energy, and the persistence and mobility of energy fluctuations, in thermal equilibrium . We have found that soft anharmonic chains (that is, chains with soft interaction potentials) experience greater energy fluctuations than harmonic chains, which in turn experience greater fluctuations than hard chains. This is a direct consequence of the virial theorem. We have also established that fluctuations are mobile in harmonic and hard chains but not in the soft. Most interestingly, thermal fluctuations travel most rapidly and remain localized over considerably greater distances in the hard chain. These results in turn can lead to very different transition rate statistics and effective damping coefficients for a bistable impurity embedded in each of these arrays at a given temperature .

Another interesting line of inquiry concerns the propagation of an energy pulse along such arrays. As mentioned earlier, we have found that the propagation of an energy pulse in a hard anharmonic array can be enhanced by immersing the array in a thermal bath, and that hard anharmonicity in the springs causes a tight and persistent packing of the energy .

Almost all of our results are numerical. Analytic results for the harmonic array provide some insights and point toward the possibility that an approach based on linearization methods  together with a linear response theory  might lead to analytic insights for nonlinear arrays. We are currently exploring these possibilities.

The thermal resonances and other interesting behavior that we have found for simple hard anharmonic chains are likely to be prototypical and therefore applicable to many other discrete systems with anharmonic interactions. While many other systems have been investigated in which a careful balance of parameters at zero temperature or a manipulation of external noise lead to interesting localization, resonance, and synchronization phenomena, our systems are among the simplest generic systems in which the temperature can be used as the tuning parameter to achieve similar effects.

## Acknowledgments

R. R. gratefully acknowledges the support of this research by the Ministerio de Educación y Cultura through Postdoctoral Grant No. PF-98-46573147. A. S. acknowledges sabbatical support from DGAPA-UNAM. This work was supported in part by the Engineering Research Program of the Office of Basic Energy Sciences at the U.S. Department of Energy under Grant No. DE-FG03-86ER13606.

## A Harmonic Chain

In this appendix we consider a harmonic chain and calculate analytically some of the numerical results presented or invoked in the body of the paper. Although the behavior of linear chains is in general well understood, it is useful to present results in the particular context of this work. These are not readily available in the literature.

### 1 Zero Temperature, With Signal

Consider first an infinite linear chain ($`\mathrm{}<n<\mathrm{}`$) at zero temperature, that is, in the absence of fluctuations (but with the dissipative contribution). The equations of motion then are

$$\ddot{x}_n=k(x_{n+1}+x_{n1}2x_n)\gamma \dot{x}_n,$$
(A1)
with the signal $`\dot{x}_0=A\mathrm{sin}(\omega _0t)`$ applied at site $`n=0`$. We conjecture the quasi-stationary solution

$$x_n^s(t)=\frac{A}{\omega _0}e^{\mu |n|}\mathrm{cos}(\omega _0t+|n|b)+c_n$$
(A2)
where the superscript $`s`$ denotes the presence of the signal. The constants $`\mu `$ and $`b`$ are to be determined. The additive constant $`c_n`$ must be independent of $`n`$ – it simply represents an overall translation of the chain (because it is not anchored) but otherwise does not contribute to the velocity analysis. The velocities associated with Eq. (A2) are

$$\dot{x}_n^s(t)=Ae^{|n|\mu }\mathrm{sin}(\omega _0t+|n|b).$$
(A3)
This result is clearly consistent with the imposed signal at $`n=0`$. The solution (A3) can be used to calculate the zero-temperature signal according to Eq. (6):

$$S_n^s(\omega _0)=_{\mathrm{}}^+\mathrm{}e^{i\omega _0\tau }\dot{x}_n(t)\dot{x}_n(t+\tau )𝑑\tau .$$
(A4)
The average (which in this deterministic case is only over time) eliminates rapidly oscillating contributions, leaving only a $`\delta `$-function-type of contribution. The result is the power spectral density

$$S_n^s(\omega _0)=\frac{A^2}{2}e^{2|n|\mu }.$$
(A5)
To find the constant $`\mu `$ we substitute (A2) into (A1) and set the coefficients of $`e^{i\omega _0t}`$ and those of $`e^{i\omega _0t}`$ equal to zero (one resulting equation is simply the complex conjugate of the other). This immediately leads to the relation

$$\omega _0^2=k(e^{\mu +ib}+e^{\mu ib}2)i\gamma \omega _0.$$
(A6)
We make the substitution $`ue^{\mu +ib}`$ and note the symmetry in $`u`$ and $`1/u`$:

$$k\left(u+\frac{1}{u}2\right)+\omega _0^2i\gamma \omega _0=0.$$
(A7)
Next we multiply through by $`u`$ and solve the resulting quadratic equation:

$$u_\pm =\frac{(\omega _0^22ki\gamma \omega _0)\pm \sqrt{(\omega _0^22ki\gamma \omega _0)^24k^2}}{2k}.$$
(A8)
Note that $`u_+u_{}=1`$, so if one solution is $`e^{\mu +ib}`$, then the other is $`e^{\mu ib}`$. Which is which is not clear at this point. $`\mu `$ must be positive for a physically acceptable solution, and this will be used below to sort out the choice.

To extract $`\mu `$ we calculate the ratio $`u_+/u_{}^{}`$, which is either $`e^{2\mu }`$ or $`e^{2\mu }`$. If $`u_+=e^{\mu +ib}`$ then $`u_+/u_{}=e^{2\mu }`$. If, on the other hand, $`u_+=e^{\mu ib}`$ then $`u_+/u_{}=e^{2\mu }`$. Which it is will be seen at the end when we inspect the magnitude of the result.

In calculating the ratio $`u_+/u_{}^{}`$ one must exercise caution in taking complex conjugates because of the square root in Eq. (A8); it is not appropriate to simply change every $`i`$ to a $`i`$ in such an expression if one does not know the signs of the terms inside the square root. It is helpful to write each complex number in terms of an amplitude and a phase:

$$(\omega _0^22ki\gamma \omega _0)2k\alpha e^{i\beta }$$
(A9)
so that

$`2k\alpha `$ $`=`$ $`\left[(\omega _0^22k)^2+\gamma ^2\omega _0^2\right]^{1/2},`$ (A10)
$`\beta `$ $`=`$ $`\mathrm{tan}^1\left({\displaystyle \frac{\gamma \omega _0}{\omega _0^22k}}\right).`$ (A11)
Similarly, with

$$(\omega _0^22ki\gamma \omega _0)^24k^24k^2\nu ^2e^{2iϵ}$$
(A12)
we find

$`4k^2\nu ^2`$ $`=`$ $`\omega _0\left[\omega _0^2(\omega _0^24k\gamma ^2)^2+4\gamma ^2(\omega _0^22k)^2\right]^{1/2},`$ (A13)
$`2ϵ`$ $`=`$ $`\mathrm{tan}^1{\displaystyle \frac{2\gamma (\omega _0^22k)}{\omega _0(\omega _0^24k\gamma ^2)}}.`$ (A14)
Then

$$u_+=\alpha e^{i\beta }+\nu e^{iϵ},u_{}=\alpha e^{i\beta }\nu e^{iϵ},$$
(A15)
and the ratio $`u_+/u_{}^{}`$ is thus real. Multiplying top and bottom of this ratio by $`u_{}`$ gives

$$\frac{u_+}{u_{}^{}}=\frac{u_+u_{}}{u_{}^{}u_{}}=\frac{1}{|u_{}|^2}.$$
(A16)
Further and explicitly

$$|u_{}|^2=(\alpha e^{i\beta }\nu e^{iϵ})\times (\alpha e^{i\beta }\nu e^{iϵ})=\alpha ^2+\nu ^2+2\alpha \nu \mathrm{cos}(\beta ϵ).$$
(A17)
All that remains is the evaluation of the resulting expressions. In summary

$$\alpha ^2+\nu ^2+2\alpha \nu \mathrm{cos}(\beta ϵ)=e^{2\mu }$$
(A18)
or

$$\alpha ^2+\nu ^2+2\alpha \nu \mathrm{cos}(\beta ϵ)=e^{+2\mu }.$$
(A19)
The choice is determined by whether the result is $`<1`$ (in which case the first equality holds) or $`>1`$ (in which case it is the second). The quantities $`\alpha `$, $`\nu `$, $`\beta `$, and $`ϵ`$ are thus completely defined in terms of the system parameters.

We have tested our numerical simulations against the prediction (A5) for large ranges of parameter values and have found agreement to five significant figures.

The solutions presented in this section do not obey $`N`$-periodic boundary conditions. It is fairly straightforward using an image-like method to construct explicitly periodic solution by applying the signal at $`0,\pm N,\pm 2N,\mathrm{}`$ but the changes would be exponentially small in $`N`$ and unimportant for sufficiently long chains. The quasi-stationary behavior assumed above for sites not too distant from $`n=0`$ sets in long before the signal at $`n=0`$ reaches sites $`\pm N`$.

### 2 Finite Temperature, No Signal

Now we solve the equations of motion

$$\ddot{x}_n=k(x_{n+1}+x_{n1}2x_n)\gamma \dot{x}_n+f_n(t)$$
(A20)
without a signal. We are again interested in the stationary behavior.

We define the Fourier transform and its inverse,

$$y_q=\underset{n=0}{\overset{N1}{}}x_ne^{2\pi iqn},x_n=\frac{1}{N}\underset{q=0}{\overset{N1}{}}y_qe^{2\pi inq}.$$
(A21)
Transforming Eq. (A20) immediately leads to

$$\ddot{y}_q+4k\mathrm{sin}^2\left(\frac{\pi q}{N}\right)y_q+\gamma \dot{y}_q=F_q(t)$$
(A22)
where the inhomogeneous term is the transform of the noise:

$$F_q(t)\underset{n=0}{\overset{N1}{}}f_n(t)e^{2\pi iqn}.$$
(A23)
The solution of this second order inhomogeneous differential equation is of standard form. The initial conditions (which we will take to be zero) are unimportant since we seek the long-time behavior:

$$y_q^0(t)=\underset{t\mathrm{}}{lim}\frac{1}{[r_2(q)r_1(q)]}_0^tF_q(\tau )\left(e^{r_2(q)(t\tau )}e^{r_1(q)(t\tau )}\right)𝑑\tau $$
(A24)
and consequently

$$\dot{y}_q^0(t)=\underset{t\mathrm{}}{lim}\frac{1}{[r_2(q)r_1(q)]}_0^tF_q(\tau )\left(r_2(q)e^{r_2(q)(t\tau )}r_1(q)e^{r_1(q)(t\tau )}\right)𝑑\tau ,$$
(A25)
where

$$r_{1,2}(q)=\frac{\gamma }{2}\pm \sqrt{\left(\frac{\gamma }{2}\right)^24k\mathrm{sin}^2\left(\frac{\pi q}{N}\right)}.$$
(A26)
The superscript $`0`$ is used to stress the absence of a signal.

With this result and the correlation function for the thermal fluctuations that follows immediately from Eq. (4),

$$F_q(t_1)F_q^{}(t_2)=2\gamma k_BT\delta (t_1t_2)N\delta _{q,q^{}},$$
(A27)
we obtain upon integration and Fourier inversion the velocity correlation function

$`C_n^0(\tau )`$ $``$ $`\dot{x}_n(t)\dot{x}_n(t+\tau )`$ (A28)
$`=`$ $`{\displaystyle \frac{\gamma k_BT}{N}}{\displaystyle \underset{q=0}{\overset{N1}{}}}{\displaystyle \frac{1}{[r_2^2(q)r_1^2(q)]}}\left(r_1(q)e^{r_1(q)t}r_2(q)e^{r_2(q)t}\right).`$ (A29)
We note in passing that this reduces to the standard correlation function for a single harmonic oscillator in a heat bath, $`C_n^0(\tau )k_BTe^{\gamma \tau }`$, when the coupling coefficient $`k0`$. The associated power spectral density is

$`S_n^0(\omega )`$ $`=`$ $`2{\displaystyle _0^{\mathrm{}}}𝑑\tau C_n^0(\tau )\mathrm{cos}(\omega \tau )`$ (A30)
$`=`$ $`{\displaystyle \frac{2\gamma k_BT\omega ^2}{N}}{\displaystyle \underset{q=0}{\overset{N1}{}}}{\displaystyle \frac{1}{[r_1^2(q)+\omega ^2][r_2^2(q)+\omega ^2]}}.`$ (A31)
Figure 9 shows the predicted spectrum (smooth curves) and the corresponding simulation results for two different temperatures. The agreement is typical of broad parameter ranges.

### 3 Finite Temperature and Signal: SNR

Finally, we deal with the full equation of motion

$$\ddot{x}_n=k(x_{n+1}+x_{n1}2x_n)\gamma \dot{x}_n+f_n(t)$$
(A32)
with the signal $`\dot{x}_0=A\mathrm{sin}(\omega _0t)`$ applied at site $`n=0`$. The fluctuation $`f_0`$ is set to zero since the velocity of that site is fixed by the signal.

The Fourier transform is defined as before, but in transforming the equations of motion we must take into account that (A32) does not hold for $`n=0`$. We thus multiply (A32) by $`e^{2\pi iqn}`$ and sum over $`n`$ but only from $`n=1`$ to $`N1`$:

$$\underset{n=1}{\overset{N1}{}}x_ne^{2\pi iqn}=y_qx_0.$$
(A33)
Proceeding in this manner we obtain in place of Eq. (A22)

$$\ddot{y}_q+4k\mathrm{sin}^2\left(\frac{\pi q}{N}\right)y_q+\gamma \dot{y}_q=U_q(t)$$
(A34)
where the inhomogeneous term now is

$$U_q(t)=F_q(t)+k\left(2x_0x_1x_{N1}\right)+\gamma \dot{x}_0+\ddot{x}_0.$$
(A35)
The terms containing $`x_0`$ and its derivatives are known; $`x_1`$ and $`x_{N1}`$ (they are equal by symmetry) would have to be found “self-consistently” ex-post (see below).

Equation (A34) is again a linear inhomogeneous differential equation with constant coefficients. Its solution has two additive contributions:

$$y_q(t)=y_q^{hom}(t)+y_q^{inh}(t)$$
(A36)
The portion $`y_q^{hom}(t)`$ is the solution of the homogeneous equation, that is, with $`U_q(t)=0`$, and the contribution $`y_q^{inh}(t)`$ is due to the inhomogeneity.

Suppose that $`F_q(t)`$ were equal to zero. The solution in this case, which we denote by $`y_q^s(t)`$ is exactly the solution of Sec. A 1; that is, the inverse transform of

$$y_q^s(t)=y_q^{s,hom}(t)+y_q^{s,inh}(t)$$
(A37)
at long times must be precisely Eq. (A2). Now consider the consequences of again including $`F_q(t)`$. The homogeneous part of the solution has to be exactly as before, i.e., $`y_q^{hom}=y_q^{s,hom}`$ because the homogeneous part of the equation has not changed, the constraint at $`x_0`$ also has not changed, and initial conditions are in any case immaterial at long times. However, the inhomogeneous part of the solution does change, and it does so in two ways. First and most important is the direct contribution of the thermal fluctuations to the solution. Second, and less important, is the change in $`x_1`$ and $`x_1`$ (which is the same as that of $`x_1`$) caused by the addition of the noise, and the change that this in turn imposes on the other displacements. This contribution is more difficult to calculate and affects only sites near $`n=0`$. Ignoring the latter we thus assume the solution

$$y_q(t)=y_q^s(t)+y_q^0(t)$$
(A38)
where $`y_q^s(t)`$ is given in Eq. (A37) and $`y_q^0(t)`$ in Eq. (A24).

The correlation function and spectrum are then simply the sums of the correlation functions and spectra for the system with the signal but no fluctuations (zero temperature) and those for the thermalized system in the absence of the signal. Thus, the spectrum is the sum of (A5) (appropriately weighted by a delta function that places it at $`\omega _0`$) and (A31). In particular, the SNR for the harmonic system then finally is

$$SNR_n(\omega _0)=\mathrm{log}_{10}\left(\frac{\frac{A^2}{2}e^{2n\mu }}{\frac{2\gamma k_BT\omega ^2}{N}_{q=0}^{N1}\frac{1}{\left[r_1^2\left(q\right)+\omega ^2\right]\left[r_2^2\left(q\right)+\omega ^2\right]}}\right)$$
(A39)
which decreases monotonically with $`T`$ as well as with $`n`$. This result is consistent with the linear response theory approach introduced and widely applied by Dykman et al. . Comparisons of this analytic result with numerical simulations are shown in Fig. 10 (cf. Fig. 4). Equation (A39) clearly captures the correct behavior; the small discrepancies between analytic and numerical results (which always lie above the analytic curves) reflect the omitted terms discussed above.


format: text

text:
# Solitons and vortices in ultracold fermionic gases

## ACKNOWLEDGMENTS

K.R. acknowledges support of the Subsidy by Foundation for Polish Science.


text:
# AGN and Cooling Flows

## 1. Introduction

Andy Fabian’s talk was packed with images from the new generation of X-ray telescopes. Images and spectra from Chandra and XMM will enliven the debate between the proponents of the steady-state cooling-flow paradigm and those who argue that unsteady heating by AGN is important. I want to mention some things that need to be borne in mind if the debate is to be fruitful. First let us recall the essentials of the cooling-flow model.

## 2. The Cooling-Flow Paradigm

The cooling-flow model was first outlined by Cowie & Binney (1977) and was elaborated by Fabian & Nulsen (1977), Nulsen (1986) and others. Extensive reviews can be found in Sarazin (1988) and Fabian (1994).

The X-ray emitting gas is considered to be in a quasi-steady state. It is extremely inhomogeneous in that at any radius there is hot diffuse gas and cooler, denser gas in approximate pressure equilibrium. The coolest gas rapidly cools to very low temperatures and disappears from the X-ray sky. Its place is taken by initially warmer gas, that cools to form cool gas. The place of the warmer gas is taken by slightly warmer gas still, which cools to form just warm gas, and so on right up to the hottest gas. Overall we have a steady inward flow of material of different densities, which gains internal energy as the gravitational field compresses it. Since in any comoving volume of the matter there is a steady trickle of material ‘dropping out’ at extremely low temperatures and high densities, the density of the X-ray observable gas does not rise as the gas moves inwards as rapidly as it would if the gas were homogeneous.

By adjusting the radial profile of the mass dropout, or, equivalently, the rate $`\dot{M}(r)`$ at which X-ray observable gas flows in across a sphere of radius $`r`$, it is possible to obtain a perfect match to the azimuthally-averaged X-ray surface-brightness profile. One finds that $`\dot{M}r`$ approximately. Once the brightness profile has been fitted to the data and a metallicity chosen, the X-ray spectrum at each radius can be predicted.

## 3. Problems of the Cooling-Flow Model

Prior to 1987 it was thought that Field’s (1965) thermal instability would naturally give rise to the multiphase medium that permits mass dropout. Malagoli, Rosner & Bodo (1987) showed that Field’s analysis was entirely misleading in a case such as this, in which the cooling medium has a clear gradient in specific entropy and is in approximate hydrostatic equilibrium in a gravitational field. In this case a cool, over-dense parcel of gas (in which the specific entropy will be abnormally low) simply sinks until it reaches the radius at which the bulk of the gas shares its specific entropy. Subsequently, the parcel executes vertical oscillations, that will in practice be strongly damped by turbulence.

While some of us abandoned the cooling-flow model on reading Malagoli et al, others went in with rubber bands to fix the model up. The rubber bands were called magnetic field lines, and they pinned over-dense parcels in place. Clearly, these field lines have to connect cooler gas to the hot embedding matrix, and be fairly straight because they are under significant tension. Consequently, heat will be efficiently conducted down them. Fields of the expected strength can only pin clouds smaller in radius than $`20`$pc, and these will be evaporated away within a Hubble time unless the thermal conductivity along field lines is several orders of magnitude smaller than expected.

For some years the spectra predicted by the cooling-flow model have conflicted with observation in that the data show fewer photons at low energies, $`1`$keV, than are predicted. This has been attributed to ‘excess absorption’ by cool gas, probably the very same gas that has dropped out of the flow. The new generation of X-ray telescopes have shown that any absorption is caused by gas at $`10^6`$K, no cooler gas being detectable (e.g., Böhringer et al., 2000). Searches at radio frequencies for emission by very cold gas have repeatedly failed to detect significant reservoirs of gas, and very strong limits have been set on populations of young stars, from the highest masses down to significantly below one solar mass (Prestwich et al., 1997).

Before the work of Malagoli et al., mass dropout seemed plausible in that the inhomogeneities upon which it depended would be generated by thermal instability from tiny seed inhomogeneities. Once one had read Malagoli et al, it was apparent that one not only had to posit the existence of strongly non-linear inhomogeneities in the initial conditions, but one also had to take steps to prevent their prompt erasure. In these circumstances, it seemed to some of us advisable to stand back and recall that cooling is fastest where the density is highest, and this is known to be at the centre of the cooling flow. Surely cooling to very low temperatures will occur first at the centre? And will such cooling not feed the central black hole and provoke it to an outburst? Might that outburst not restructure the X-ray emitting gas, and thus regulate the cooling process?

## 4. Heating by AGN

This proposal that cooling flows are in a dynamical equilibrium with a central AGN has the attraction of a fresh start. It sweeps away cobwebs in the form of magnetically pinned over-densities and strange absorbing screens. But one has to admit that at this stage it is little more than a start: the detail is complex and formidably difficult to work out. Of course we must seek critical confrontation with the data wherever possible, but we must avoid damning the overall picture if we find that the data do not support a naive expectation of what the model might predict.

Working out the observational consequences of AGN feedback on cooling flows is hard because

* the system is constantly far from a steady state;
* the geometry is inherently complex;
* to predict the radio data one needs to know the distribution of both the ultra-relativistic electrons and the magnetic field.

The cooling-flow model is very much simpler and has been rather completely worked out. Hence the collision between the cooling-flow model and the heating picture is asymmetrical: on the one side we have a model with well defined predictions, and on the other a general picture from which it is very hard to wring firm predictions. We rightly have a preference for the simple over the complex, for the definite over the vague. Consequently, the majority of people working in the field have clung tenaciously to the cooling-flow model in the face of ever-increasing difficulties. Never the less, I submit that the heating model enjoys a decisive advantage in another important criterion by which we judge competing theories: a priori plausibility in the light of our knowledge of the whole of physics and astronomy. We must make a determined effort to work the picture up into a well-explored model and judge it sternly only when we have extracted secure predictions from it.

### 4.1. Time-dependence

The fundamental difference between the heating picture and the cooling-flow model is that in the former partial derivatives with respect to time will always be important. One expects an individual source to cycle through a series of stages. In the first stage the system cools towards a catastrophe. In the second, the AGN is rejuvenated and blasts its immediate surroundings. In the third turbulence dies down, and the entropy stratification of the gas is restored before the first phase of the cycle resumes.

At each stage of the cycle one expects the source to look quite different. The X-ray gas is expected to respond more slowly than the synchrotron-emitting plasma, which will respond less quickly than any central continuum source. In principle it should be possible to make predictions about how many sources should be seen in each stage, but we have neither such predictions nor the demographic knowledge that would be required to test them. Proponents of the cooling-flow model are accustomed to imagining the objects to be in steady states and one has to watch for them implicitly carrying this assumption over to the heating picture.

Some important stages, such as that in which material cools to below $`1`$keV, may be very short-lived and rarely observed. The phase of violent energy output by the central black hole is likely to be subdivided into extremely short-lived episodes, since the relevant dynamical time is a matter of months only.

### 4.2. Complex geometry

Even if both the gravitational potential and the central heat source were spherical, the resulting X-ray source would not be because heating from below is inherently unstable: hot plasma will rise in some directions, while cool plasma falls in others. Much of the energy released at the centre is likely to emerge as ordered kinetic energy in a jet or other strongly collimated outflow. One cannot credibly model the system in spherical symmetry (though I have tried! – see Binney & Tabor, 1995). I am skeptical that much useful can be achieved even in cylindrical symmetry (but see D’Ercole & Ciotti, 1998). Fortunately, fully three-dimensional simulations are now just about feasible, and we can hope for fairly realistic models to become available. This will remain a computationally challenging problem for some time, however, because the range of spatial scales and densities involved is extreme – from the parsec scale of the energy source to the megaparsec scale of the X-ray halo; from a proton density in excess of $`10^8`$m-3 near the centre of the cooling flow to densities many orders of magnitude lower in the plasma (e±?) of a relativistic jet; from inflow speeds of $`10`$km$`\text{s}^1`$ in the bulk of the cooling flow to the speed of light in the jet. A further complication is the possibility that the IGM becomes multiphase – during the cooling phase this will not occur, but gas heated from below is liable to generate distinct phases, as has long been observed in the disk of the Milky Way.

One should beware simplistic assumptions about how a source’s observable properties will be affected by a jet. Displacement of thermal plasma by relativistic plasma will suppress X-ray emission. On the other hand, X-ray emission will be enhanced when thermal plasma is shocked or pushed up from lower, cooler layers. Correspondingly, X-ray emission will be suppressed when higher-entropy plasma sinks downwards to take up space left by the upward movement of material. From the case of M87 we know too that the synchrotron and inverse-compton processes can sometimes enhance X-ray emission.

### 4.3. Radio emission

There is a great temptation to feel that one is seeing a true picture of the relativistic plasma when one looks at a map of radio-continuum emission. In fact, as Katherine Blundell (this volume) has emphasized, the synchrotron emissivity at a given wavelength and for a given electron spectrum is a sensitive function of the magnetic field strength. A simple argument suggests that the latter is probably very non-uniform.

This argument starts from the observation that in a perfectly conducting medium the magnetic field obeys the same equation as the vorticity (e.g. Kulsrud et al., 1997). Observations of rivers in spate and wind patterns teach us that in high Reynolds number flows, vorticity is strongly concentrated into sheets and vortex lines; in the language of non-linear systems, vorticity is ‘intermittent’. We must expect magnetic fields in turbulent conducting media to be similarly intermittent, and this expectation is confirmed by recent observations of the Sun (Hagenaar et al., 1999). Hence, when we look at the delicate filamentary pattern of synchrotron emission in an X-ray halo such as that of M87, we are probably seeing merely the regions of enhanced $`B^2`$, while the relativistic particles are more uniformly distributed. In particular, it is likely that both X-ray emission and synchrotron emission derive from the same plasma; the synchrotron emission comes from the high-energy tail of the particle distribution, and is localized merely because $`B^2`$ is.

If it is true that ultra-relativistic electrons are space-filling out to the outer boundary of the radio halo (which in the case of M87 lies at more than half the cooling radius), substantial $`P\mathrm{d}V`$ work is likely to have been done by the radio source on the surrounding thermal plasma.

The radiating electrons are ultra-relativistic, and it is dangerous to assume that these exotic particles dominate the total energy in suprathermal particles; certainly in the solar neighbourhood, the cosmic-ray energy is dominated by the least energetic particles. Suprathermal particles with Lorentz factors $`\gamma 1`$ have extremely long lifetimes because their synchrotron radiation is negligible and they are too fast to be significantly Coulomb scattered. Consequently, they can persist and be dynamically important long after synchrotron emission has become unobservably faint.

Observations of jets suggest that the ultra-relativistic particles that dominate radio maps are generated by jets with bulk Lorentz factors of a few. These jets are unlikely to be produced by the accretion torus around the AGN; they are probably produced by an exotic process, such as the Blandford-Znajek (1977) effect, that taps the rotational energy of the black hole. The fate of the comparable or larger quantity of energy that is released within the accretion disk has long been a puzzle. The ADAF model (Narayan & Yi, 1995) posits that it is swallowed by the black hole. Blandford & Begelman (1999) have strongly criticised this model \[which has recently encountered difficulty matching observed spectral energy distributions (Di Matteo et al., 2000)\] and argue that the energy is probably carried away by a wind off the accretion torus, rather than swallowed by the black hole. This conjecture chimes with observations of accreting stars, from SS 433 to the sources of Herbig-Haro objects, which imply that an accreting system is invariably associated with a bipolar flow at velocities that range up to the largest Kepler velocity of accreting material.

The implications for the heating picture of a significant fraction of the accretion energy from an AGN emerging as a sub-relativistic bi-polar flow would be considerable. The flow would probably be more steady than the observationally much more conspicuous high-$`\gamma `$ jet that would occasionally flicker to life at its centre. It would involve very much more momentum than the high-$`\gamma `$ jet, and it would be hard to detect through synchrotron radiation because in the sub-relativistic shocks that would arrest it, very little energy would go into ultra-relativistic particles. It is likely to generate structure in X-ray maps that is hard to trace in radio maps because it is not associated with strong synchrotron emission. Hence one should treat with caution claims that a weak correlation between the radio and X-ray structures of an object imply that the AGN is not heating the gas.

## References

Binney, J.J. & Tabor, G.R., 1995, MNRAS, 276, 663

Blandford, R.D. & Begelman, M.C., 1999, MNRAS, 303 L1

Blandford, R.D. & Znajek, R.L., 1977, MNRAS, 179, 433

Böhringer et al., 2000, A&A in press (astro-ph/0011459)

Cowie, L.L. & Binney, J.J., 1977, ApJ, 215, 723

D’Ercole, A. & Ciotti, L., 1998, ApJ, 494, 535

Di Matteo, T., Quataert, E., Allen, S.W., Narayan, R. & Fabian, A.C., 2000, MNRAS, 311, 507

Fabian, A.C., 1994, ARAA, 32, 277

Fabian, A.C. & Nulsen, P.E.J., 1977, MNRAS, 180, 479

Field, G.B., 1965, ApJ, 142, 531

Hagenaar, H.J., Schrijver, C.J. Title, A.M. & Shine, 1999

Kulsrud, R.M., Cen, R., Ostriker, J.P. & Ryu, D., 1997, ApJ, 480, 481

Malagoli, A., Rosner, R. & Bodo, G., 1987, ApJ, 319, 632

Narayan, R. & Yi, 1995, ApJ, 452, 710

Nulsen, P.E.J., 1986, MNRAS, 221, 377

Prestwich, A.H., Joy, J., Luginbuhl, C.B., Sulkanen, M. & Newberry, M., 1997, ApJ, 477, 144

Sarazin, G., 1988, ‘X-ray Emission from Clusters of Galaxies’, Cambridge University Press, Cambridge.


text:
# Diffusion of Point Defects in Two-Dimensional Colloidal Crystals

## Abstract

We report the first study of the dynamics of point defects, mono and di-vacancies, in a confined 2-D colloidal crystal in real space and time using digital video microscopy. The defects are introduced by manipulating individual particles with optical tweezers. The diffusion rates are measured to be $`D_{mono}/a^23.27\pm 0.03`$Hz for mono-vacancies and $`D_{di}/a^23.71\pm 0.03`$Hz for di-vacancies. The elementary diffusion processes are identified and it is found that the diffusion of di-vacancies is enhanced by a dislocation dissociation-recombination mechanism. Furthermore, the defects do not follow a simple random walk but their hopping exhibits memory effects, due to the reduced symmetry (compared to the triangular lattice) of their stable configurations, and the slow relaxation rates of the lattice modes.

Uniform colloidal spheres suspended in a solvent, under appropriate conditions can self-assemble into ordered crystalline structures . Using these crystals as a model system, a great variety of problems in materials science, physical chemistry and condensed matter physics have been investigated during the past two decades. However, the properties of simple structural defects, such as vacancies and interstitials, dislocations, etc., ubiquitous in real solids, have received little, if any, attention. Considering the unique possibility of real space and real time imaging, colloidal systems provide a great potential for research in this direction.

Recently, it has been demonstrated that point defects can be created in two-dimensional colloidal crystals by manipulating individual particles with optical tweezers. Direct imaging of these defects verified that their stable configurations have lower symmetry than the underlying triangular lattice, as predicted by numerical simulations for a number of two-dimensional systems. The same study revealed that point defects can dissociate into pairs of well-separated dislocations, a topological excitation permitted only in two but not three dimensions. In this Letter we report the first study of the dynamics of mono and di-vacancies in two-dimensional colloidal crystals. We discover the first evidence that topological excitations into dislocation pairs enhance the diffusion of di-vacancies. Moreover, the hopping of the defects does not follow a pure random walk, but exhibits surprising memory effects.

The colloidal crystals are prepared by 0.360 $`\mu `$m diameter negatively charged polystyrene-sulfate micro spheres (Duke Scientific, #5036), suspended in highly deionized water. The particles crystallize due to strong electrostatic interaction. A single layer of particles is confined between two fused silica substrates separated by $``$2 $`\mu `$m, resulting in a single layer colloidal crystal with a lattice constant $`a1.1\mu `$m. The experimental details can be found in. Point defects are created by trapping a particle with optical tweezers and dragging it away from its lattice site. The dynamics of such single, isolated point defects are recorded on videotape, using a CCD video-camera. Each individual frame (1/60th of a second) was acquired on a PC and processed using a particle-tracking algorithm. A triangulation algorithm is used to identify the mis-coordinated particles and characterize the configuration of the defect.

We tracked a number of mono and di-vacancies for long enough times (up to 40 seconds) to obtain reliable measurements of their diffusion constants. Typical trajectories are shown in Fig.1. The position of the defect $`\stackrel{}{x}(t)`$ is defined as the average of the positions of the 5-fold coordinated particles around the core of the defect. $`\stackrel{}{x}(t)`$ is determined with a resolution of $`0.1a`$. Thus one can identify the defect in a particular lattice site or in-between two neighboring sites. To estimate the diffusion constants, we take the average of the squared displacements $`<\mathrm{\Delta }\stackrel{}{x}^2(\delta t)>`$ as a function of the time separation $`\delta t`$, running through all the available trajectories. A linear relationship $`<\mathrm{\Delta }\stackrel{}{x}^2>`$ vs. $`\delta t`$ is obtained for the range 1-5 seconds (Fig.2A). The data in this range are weighted by their errors and fitted to a zero y-intercept linear relation, whose slope gives the respective diffusion constants. The results for mono and di-vacancies are: $`D_{mono}/a^23.27\pm 0.03`$Hz, $`D_{di}/a^23.71\pm 0.03`$Hz.

The diffusion constants measured above give the combined contribution of all the possible diffusion pathways. In this Letter, we reveal the role of topological excitations on the diffusion of point defects. As was shown in , the possible stable configurations of mono and di-vacancies are not just the ”well-identified” point defect configurations but there are also configurations where the point defects appear as separated dislocation pairs. The defects are observed to dissociate into two dislocations that glide from site to site and finally recombine at a new position that can be several lattice constants away from the original one. The individual dislocations move faster than the hopping of mono or di-vacancies and this topological excitation has a lifetime $`\tau _o`$ of a fraction of a second. Such a process is a distinct diffusion pathway, giving a contribution to the diffusion constant $`D_{pair}\tau _o^1l^2`$, where $`l^2`$ is the mean square separation of the dislocation pair.

Mono-vacancies spend considerably less time as a dislocation pair than di-vacancies and the aforementioned diffusion mechanism is not important in their case. One can clearly recognize the different behavior exhibited by mono and di-vacancies, as illustrated in Fig.3, where the separation of the dislocations comprising the defect is shown as a function of time. The mono-vacancy is observed primarily to fluctuate between configurations in which the dislocations are bound, with very rare excursions into configurations where the dislocations dissociate. In contrast, the di-vacancy frequently spends long intervals (a few seconds) as a dislocation pair, with the two dislocations dissociating and recombining several times during this period.

The combined diffusion constant for di-vacancies is $`D_{di}=pD_{pair}+(1p)D_{di}^o`$, where $`p0.25`$ is the fraction of time excited into a dislocation pair, and $`D_{pair}`$,$`D_{di}^o`$ are the contributions from the dissociation-recombination mechanism and from the ”bare” di-vacancy hopping. In order to estimate the diffusion constants for the two different mechanisms, we link into two separate trajectories the intervals spent as a dislocation pair and as a ”bare” di-vacancy. To avoid systematically picking up hoping events, the shortest interval was 1 second (larger than the time between jumps). Carrying out the same procedure as for the original trajectories, the two diffusion constants are found to be $`D_{pair}/a^2=5.90\pm 0.09`$Hz and $`D_{di}^o/a^2=2.97\pm 0.03`$Hz (Fig.2B). These numbers combine to give the measured value $`D_{di}/a^2=3.70`$Hz.

We now turn our attention to the details of the underlying stochastic hopping processes. Naively, one would expect the trajectories of the defects to be pure random walks on a triangular lattice, i.e. when hopping, the system chooses with equal probability between the six neighbor sites and consecutive events are uncorrelated. However, a careful analysis reveals a peculiar feature: their motion looks nothing like a 2D random walk, since typically it is restrained along particular directions. The existence of easy directions for diffusion is related to the low symmetry of the stable configurations of the relaxed lattice around the defect. Therefore, the probability for hopping towards the six nearest neighbors is not the same for all of them, but is higher along particular lines of symmetry, Fig.4. For di-vacancies excited into a pair of separated dislocations, the easy direction coincides with the gliding direction of the two dislocations. The easy direction for diffusion remains unchanged for up to a few seconds.

The existence of easy directions for hopping makes the trajectories of the defects quasi-one-dimensional, at least for a few consecutive steps. This effect is demonstrated by measuring the probability $`P_o(n)`$ for returning to the original position after the defect has executed $`n`$ steps. It is found that the measured $`P_o(n)`$ for small $`n`$ resembles the one for random walk on a 1-D chain, rather that for random walk on a 2-D triangular lattice. This is illustrated in Fig.5: $`P_o(n)`$ clearly deviates from the expected 2-D behavior, following the characteristic oscillations of the 1-D random walk of period=2 steps.

A second surprising peculiarity involves memory during site-to-site hopping. The defects appear to spend some time around a lattice site, attempting to hop towards neighboring sites by changing into a split intermediate configuration between the two lattice points. Depending on the symmetry of the relaxed lattice in the original configuration (2 or 3-fold), the defect is seen to attempt hopping towards either two or three of the six neighbors. Oddly, the defect typically hops back to the original site, finally getting to the next site only after several attempts. The back-hopping to forward-hopping ratio is found to be $`3:1`$.

This memory effect originates from the slow relaxation rates of distotrions in the overdamped colloidal crystal: the core region of the defect can rearrange rapidly ($`10Hz`$), however the rest of the lattice around the defect cannot follow at the same rate (e.g. distortions with wavelength $`10a`$ relax at $`1Hz`$). As a vacancy changes from a configuration at one site to a split intermediate configuration between two sites, the rest of the lattice does not immediately relax, and seems to ”push” the vacancy back to its initial site.

In summary, the dynamics of point defects in colloidal crystals are studied in real space and time. Di-vacancies in particular are observed to dissociate into two dislocations that can separate several lattice spacings before recombination. This distinct mechanism enhances the diffusion rate of di-vacancies compared to their ”bare” one. The reduced symmetry of the defects’ stable configurations results in easy directions for diffusion. We have also found that the site-to-site hopping seems to exhibit memory effects accounted for by the slow relaxation of the short wavelength modes of the crystal. We acknowledge helpful discussions with Prof. S.C. Ying and Prof. D.A. Weitz. This work was supported by NSF Grant No. DMR-9804083, the Petroleum Research Fund, the Research Corporation, and the A.P. Sloan Foundation.