text:
# On the Aging Dynamics in an Immune Network Model

## 1 Introduction

In this paper we discuss a model for the evolution of the immune repertoire of B cells, which are responsible for the humoral immune responses. B cells belong to one of the main classes of white blood cells: the lymphocytes. These cells carry on their surface the order of $`10^5`$ molecular receptors (proteins) and once activated they produce antibodies, which are copies of their molecular receptor. During the life time of an individual the immune system is able to produce the order of $`10^{11}`$ different antibodies or different populations of B cells. The antigen (virus, bacteria, poison, cellular residue, etc) is not recognized as a whole but by its epitopes, which are patches on its structure that may be recognized by specific sites of the antibody molecules. By pattern recognition different antibodies will mark the epitopes of a given antigen, therefore forming a complex that will be eliminated by macrophages (another class of white blood cells) .

According to clonal selection theory , elements that challenge the immune system will determine the populations (clones) of B cells that will proliferate: those populations will produce antibodies which will be able to recognize different epitopes of the specific antigen. The immune network theory , however, is based on the fact that the antibodies (and molecular receptors) are able to recognize and to be recognized, and therefore during the immune response there are different types of interaction: antigen-antibodies and antibodies-B cells. In other words, when a given population of B cells is activated by the presence of a given antigen the produced antibodies will not only mark the specific antigens but also activate new B cell populations with complementary molecular receptors, which in turn will recognize them. The increase on the concentration of these complementary populations, on their turn, will maintain the proliferation of the antigen recognizing population, installing a feedback mechanism that will keep several populations activated. This kind of dynamics will generate a functional multi-connected network among different populations of B cells that will be dynamically regulated by mechanisms of activation and suppression. The network will then play an important role on the regulation of the immune responses. Although the immune network theory is part of the current immunological thinking, there are only few experiments supporting the interaction among clones with complementary receptors and the existence of such a network . According to these experimental findings, if the network exists only $`1020\%`$ of the populations will belong to it, the rest of the immunocompetent populations remaining at rest.

Recently we have successfully used a mathematical model  (inspired in a previous one proposed by de Boer, Segel and Perelson ) which takes into account the main features of Jerne’s immune network theory, to simulate experiments on immunization and aging performed with mice  that could not be explained by the clonal selection theory . The simulations allowed to interpret the experimental results from the point of view of the immune network theory.

The model allows to follow the evolution of the concentrations of the different populations of B cells in discrete shape space, a formalism which maps all possible molecular receptors of a given organism into points of a $`d`$-dimensional space. To each point (receptor) we associate a clone that corresponds to the population of B cells and antibodies characterized by this molecular receptor. The concentration of each clone will be described by a three-state automaton representing low, medium and high concentrations, and the interaction among different clones is based on complementarity. As far as we know, the model  corresponds to the first successful attempt in describing the dynamics of the immune network as proposed by Jerne  and recently it could be used to reproduce experimental results performed with mice . Besides the biological implications of the results obtained in Ref. , the dynamical behavior exhibited by this model in the biologically relevant parameter region is quite interesting by itself and should be better investigated. This region has been shown  to exist for dimensions $`d2`$ and comprises a broad stripe near the transition between stable and chaotic behaviors, in which the model describes a multi-connected functional network . In this region we found the majority of the populations in the resting state (low concentration) while the activated ones may reach $`1020\%`$ of the total number of populations. The activated populations are aggregated in clusters of different sizes, which fuse and split as time passes following an aggregation and disaggregation dynamics. Therefore the largest cluster at each time step is found in a different region .

In the experiments described in ref. , 6 mice of the same linage are subjected to the immunization protocol as follows: the researchers inject ova by means of an intra-venal injection, wait for 14 days and measure the number of specific antibodies. Then they inject again the same amount of antigen, wait for 7 days and measure the amount of antibodies, and continue by repeating the same protocol every 7 days. In order to simulate the immunization protocol, the CA is subjected to specific perturbations by flipping chosen resting populations (low concentration) to their activated state (high concentration). Depending on the perturbation (damage) size, it may disappear after a few time steps or part of the damage (activated populations) may be incorporated to the network. In other words, the memory about the perturbations is due to the ability of the system to adapt (plasticity) and incorporate information about them. We have used two types of perturbations: random small ones, which correspond to the noise that mice are subjected to during the experiment, caused by the environment, and periodic large ones that will simulate the immunization protocol of multiple antigen perturbations under the same conditions . After few presentations there is a saturation of the response of the system to the perturbation. During this process, the system incorporates new information and some of the previous information is lost, keeping the number of activated sites in the network almost unchanged. This kind of behavior was also observed experimentally in mice , where saturation is related to a refractory behavior of the immune system.

We have also observed aging effects on the dynamics of this system . An older system is more rigid: the network loses plasticity to incorporate new information . A recent study  has shown that the distribution of cluster sizes during the time evolution of the system has a characteristic cluster size (exponential behavior), but the distribution of persistence times (the period during which a given population remains activated and belongs to the network) exhibits a power law behavior. While the existence of a characteristic cluster size may be related to the loss of plasticity, the power law behavior of the persistence time may be associated to the memory generated by the dynamics of the system.

The question we address here refers to how the learning process takes place dynamically and what is the cause of the loss of plasticity. The slow dynamics observed in this system presents analogies with the physical aging effects observed and reported on glass studies : as the system gets older (ages) there is a loss of plasticity for structural or molecular relaxation and less changes are observed during the relaxation time. The mechanisms underlying the slow dynamics of glassy systems are the spatial (or geometric) disorder related with the difficulty to satisfy simultaneously all microscopic interactions, a characteristic called frustration. Glasses and spin glasses exhibit a rough energy landscape with many local minima which is responsible for slowing down the relaxation towards equilibrium. Their dynamics depends on the past history of the sample and under small perturbations they relax slowly. In our biologically motivated model the non-local interactions are based on complementarity (both perfect and slightly defective matches) and reflect the activation and suppression mechanisms in the complementary regions of the shape space. The analogue to frustration, in the network, is generated by the inability of the system to incorporate information and to satisfy the constraints of the activation and suppression mechanisms. We show that the CA dynamics gives rise to a large family of periodic attractors which are very robust, and could be regarded as the analogues of the minima in the glassy systems. The lost of plasticity and aging effects will therefore be related to the non-ergodicity in the phase space.

In this paper the dynamics of the immunological responses in the network model are investigated using the common tools adopted in the study of glassy systems . In particular we will focus on the relaxation of two-time autocorrelations of the B-cell populations, the structure of the attractors and effects of perturbations and immunizations. The paper is organized as follows: in section 2 we describe the model and the procedures adopted in order to measure the correlations among different configurations of the system; in section 3 we present and discuss the results obtained for small and large perturbations and in section 4 we present a summary, prospectives for future work and some conclusions.

## 2 The Model

The model under study here is a modified version  of the model proposed by Stauffer and Weisbuch , which in turn was inspired in a previous model proposed by de Boer, Segel and Perelson  to describe the time evolution of the immune repertoire. It is a deterministic window cellular automaton model based on the shape space formalism , which describes the interactions of B-lymphocytes and antibodies, and the main mechanism underlying these interactions, which is pattern recognition (lock-key interaction). The dynamics of the system is regulated by specific interactions involving complementary molecular receptors of the different clones. The memory about the relevant antigens, presented to the system during its past history, emerges from the dynamics of the system, rather than being stored in a static registry.

To each point of a $`d`$-dimensional discrete lattice we associate a given receptor, which in turn will be described by $`d`$-coordinates representing important physical-chemical features of the receptor that differentiate one from the other . Since clones are classified according to their molecular receptor, to each point $`\stackrel{}{r}`$ of the discrete shape space we associate a three-state automaton $`B(\stackrel{}{r},t)`$ that will describe the concentration of its population over the time: low ($`B(\stackrel{}{r},t)=0`$), intermediate ($`B(\stackrel{}{r},t)=1`$) and high ($`B(\stackrel{}{r},t)=2`$).

The time evolution of the cellular automaton is based in a non-local rule: population $`B(\stackrel{}{r},t)`$ at site $`\stackrel{}{r}`$ is influenced by the populations at site $`\stackrel{}{r}`$ (its mirror image or complementary shape) and its nearest-neighbors ($`\stackrel{}{r}+\delta \stackrel{}{r}`$) (representing defective lock-key interactions). The influence on the population at site $`\stackrel{}{r}`$ due to its complementary populations is described by the field $`h(\stackrel{}{r},t)`$:

$$h(\stackrel{}{r},t)=\underset{\stackrel{}{r}^{}(\stackrel{}{r}+\delta \stackrel{}{r})}{}B(\stackrel{}{r}^{},t)$$
(1)
where for a given $`\stackrel{}{r}`$ the sum runs over the complementary shape $`\stackrel{}{r}^{}=\stackrel{}{r}`$ and its nearest neighbors. Due to the finite number of states of the population $`B`$, the maximum value of the field $`h(\stackrel{}{r})`$ is $`h_{max}=2(2d+1)`$. The updating rule is based on a window of activation, which describes the dose response function involved in B-cell activation . There is a minimum field necessary to activate the proliferation of the receptor populations ($`\theta _1`$), but for high doses of activation (greater than $`\theta _2`$) the proliferation is suppressed. The updating rule may be summarized as:

$$B(\stackrel{}{r},t+1)=\{\begin{array}{cc}B(\stackrel{}{r},t)+1\hfill & \mathrm{if}\theta _1h(\stackrel{}{r},t)\theta _2\hfill \\ B(\stackrel{}{r},t)1\hfill & \mathrm{otherwise}\hfill \end{array}$$
(2)
but no change is made if it would lead to $`B=1`$ or $`B=3`$. We define the densities of sites in state $`i`$ at time $`t`$ as $`B_i(t)`$ ($`i=1`$, 2, 3).

The initial configurations are randomly generated according to the following concentrations: $`B_1(0)=B_2(0)=x/2`$, while the remaining $`L^d(1x)`$ sites are initiated with $`B(\stackrel{}{r},0)=0`$. This model may exhibit stable or chaotic behaviors depending on the values of $`x`$, $`\theta _1`$ and $`\theta _2\theta _1`$. However it is on the transition region between the two behaviors that the model behaves like a multi-connected network .

In order to simulate the immunization protocol performed in the mice experiments we have followed the procedure which is described in Ref.  and summarized below. We have adopted the scale of 5 time steps corresponding to 1 day . While the system evolves according to the deterministic dynamics (eq. 2), small and large perturbations can be produced, by setting the state of the chosen sites at $`B(\stackrel{}{r},t)=2`$.

#### Small Perturbations

The small perturbations account for the immunological stimuli (noise) coming from the environment. The time interval between two consecutive small perturbations is a random number uniformly distributed between 1 and 100 time steps. Each perturbation corresponds to a random number of damages (from 1 to 3) introduced at regions of resting populations ($`B=0`$) which are randomly drawn (at every perturbation). The size of each damage may also vary randomly from 1 to 3 (onion-like) concentric layers around a central site (containing 7, 25 or 63 populations, respectively, in 3 dimensions).

#### Large Perturbations

The large perturbations correspond to the immunization protocol which starts at a predetermined age of the mice. Like in the experiments, we stimulate the system periodically (every 35 steps $``$ 1 week), and always in the same region (which is initially chosen at random but kept unchanged along the simulation). The damage size in this case corresponds to six layers (377 populations) around an specific site.

Previous results on this model have shown  that the response to the immunizations presents a strong dependence on the initial time (“age”) at which the periodic protocol starts (fitting experimental data extremely well for mice whose immunization protocol started at different ages). This has led us to the conjecture that the dynamical behavior of the system might be at least qualitatively similar to that of some glassy systems, despite the non-Hamiltonian nature of the CA dynamics.

One of the quantities commonly used in the study of glassy systems is the two-time autocorrelation function between the system configurations at two given times $`t`$ and $`t_w`$. A common experiment in glassy systems consists in preparing the system at a high temperature and suddenly making a quench to a low temperature. Then the system is allowed to relax up to a waiting time $`t_w`$, whose configuration is recorded. As the system continues to relax the autocorrelations between the instantaneous configurations at time $`t>t_w`$ and that at time $`t_w`$ are computed. The waiting time $`t_w`$ is called the age of the system in context of glassy systems. The monitoring of the two-time autocorrelations gives important insights on the relaxation process. The dynamics, whether stationary or not, can be readily recognized on the $`t_w`$ dependence of the autocorrelations, since in a stationary process two-time quantities depend only on time differences. Consequently by plotting correlations as a function of time difference $`tt_w`$ for different waiting times it is possible to distinguish between an essential out of equilibrium process from a stationary one. The aging processes observed in glassy systems are then related to the lack of temporal invariance .

Inspired on this approach, we will analyze the multi-connected network dynamics by defining and analyzing quantities analogous to the two-time autocorrelations for the CA:

$`C_{tot}(t,t_w)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{r}}{\overset{N}{}}}\delta (B(\stackrel{}{r},t_w),B(\stackrel{}{r},t))`$ (3)
$`C_{22}(t,t_w)`$ $`=`$ $`{\displaystyle \frac{1}{B_2(t_w)N}}{\displaystyle \underset{\stackrel{}{r}|B(\stackrel{}{r},t_w)=2}{\overset{N}{}}}\delta (B(\stackrel{}{r},t),2),`$ (4)
where $`N=L^d`$ is the total number of populations in the system. $`C_{tot}`$ and $`C_{22}`$ amount to normalized proximities (using Hamming distance as a measure) and from now on will be referred to simply as correlation functions. These quantities will be analyzed for different protocols: without any perturbation, with small perturbations (noise) and with large perturbations (immunization protocol) in order to differentiate the effects of the different kinds of perturbations.

A complementary view of the long-term behavior of the system can be obtained by looking at the attractors to which the system evolves. In the present case it will consist mainly of cycles, which will be reflected by the periodicity of the correlations. For this purpose, we have obtained the return maps of the consecutive maxima of the correlation functions. Note that the period of the maxima thus obtained does not correspond to the period of the real system, which is at least twice as long .

## 3 Results

### 3.1 Without perturbations

In previous works  it has been shown that depending on the initial concentration $`x`$ of active populations the system may exhibit periodic or chaotic behavior. Systems with low initial concentrations of active sites ($`x<x_c`$) evolve to either fixed points or orbits with short periods, while for $`x>x_c`$ chaotic attractors appear. However, the biologically relevant region is in the transition region between these two behaviors, where the system reaches one of several periodic orbits (as will be seen below) with a very long period and after a long transient. From now on, all the results have been obtained using the same parameters adopted in Ref. : $`d=3`$, $`\theta _1=h_{max}/3`$, $`\theta _2=2h_{max}/3`$ and $`x=0.26`$ (on the transition region).

Without any perturbation, the system evolves after a transient time towards a cycle with a large period, as shown in the return map of Fig. 1. We have also varied the waiting time from $`10`$ to $`100000`$ (not shown). The greater the waiting time, the closer to the attractor the system is, which is revealed by larger values of the autocorrelation functions. Once $`t_w`$ is larger than the transient, the time series for $`C_{tot}(t,t_w)`$ (and also $`C_{22}(t,t_w)`$ will include unity (see Fig. 1) and will not change for larger values of $`t_w`$.

A typical result for the time evolution of the correlation functions for $`t_w=100`$ is shown in Fig. 2. The upper panel corresponds to the time evolution of the densities $`B_1`$ (intermediate concentrations) and $`B_2`$ (high concentrations), while the lower panel corresponds to $`C_{tot}(t,t_w)`$ (open circles) and $`C_{22}(t,t_w)`$ (filled circles). Notice that while the concentrations relax to approximately constant values in a short time, $`C_{tot}`$ and $`C_{22}`$ take much longer to reach their attractors (note the logarithmic time scale). In order to study how the system relaxes towards the attractor it is more convenient to make use of $`C_{22}`$, since it measures the changes in the (more relevant) network of activated populations. In the case shown in Fig. 2 the transient time needed to attain the attractor is $`𝒪(10^4)`$ steps. It is necessary to point out, however, that this typical relaxation time is important only for the physical aspects of the dynamics. When mapped into the biological problem, it would correspond to $`5.5`$ years, which is much longer than the average life time of the mice used in the kind of experiment we simulated. Therefore the relevant behavior, from the biological point of view, happens to be in the transient of the model and not in its stationary state, a result which is interesting on its own.

### 3.2 Random Small Perturbations

How does the behavior of the system change when random small perturbations are produced on the parameter region used to simulate the real experiments performed with mice ? In order to investigate this issue, only the small perturbations described in Section 2 were produced, starting at time zero.

In Fig. 3 we compare the time evolution of the densities and correlation functions obtained for single runs for the system without perturbations and with small random perturbations. Note that in the case with small perturbations the correlations initially follow those of the purely deterministic system, showing only small differences. After about $`10^310^4`$ steps, however, they start decreasing faster, indicating some sort of cumulative effect that drives the system away from the region in phase space that it had approached until $`t_w`$. These effects are more easily noted for $`C_{22}`$. Moreover, the perturbations do not alter the behavior of the densities, as expected, since the number of activated populations is kept approximately constant by the self-regulatory mechanisms embedded in the dynamical rules. The changes are observed only in the populations that belong to the active network: in order to incorporate new information (new populations), older ones are deactivated. From the biological point of view, the realizations of the perturbations (small and large) differ from one individual to another, building the identity of each individual. At the end of life each individual will have a different history in terms of perturbations (antigen presentations) translated into the configuration of the populations belonging to the active network. This is nicely illustrated by Fig. 4, where two initially identical copies of a system undergo different realizations of the small perturbations according to the protocol described in section 2. The Hamming distance between them grows on a long time scale, revealing the mechanisms behind the behavior of the correlations in Fig. 3.

Returning to Fig. 3, it is important to stress that the changes observed in the concentrations for very long times (upper panel) are due to finite size effects. They are caused by the fact that all small perturbations are produced on regions of resting populations. For a finite system, after a long time all the possibilities will have already been explored. Increasing the size of the system, the changes on the densities disappear. This is shown in Fig. 5, where we repeat the simulations, under the same conditions of Fig. 3, for a larger system ($`L=100`$).

In order to test out the robustness of these results, we have varied the parameters controlling the protocol of the random small perturbations. For instance, by increasing the maximum number of different perturbations from 3 to 6 at each presentation, and/or by changing the maximum time interval between consecutive perturbations from $`100`$ to $`10`$, we observe the same qualitative behavior, with a faster decrease of the correlations — as expected, since in both cases the noise has been increased.

What happens when the system’s autocorrelations decrease due to the small random perturbations after it has reached a periodic attractor? Is the system driven to another cycle? To answer this question we have studied the stability of the cycles using the following procedure: we let the system evolve without perturbation towards its attractor for $`t_w=10000`$ time steps, after which we perturb the system with random small perturbations during a time interval $`\mathrm{\Delta }t_1`$. Then we turn off the perturbations and allow the system to relax during another time interval $`\mathrm{\Delta }t_2`$ after which we obtain the return map of the correlations in the following 200 time steps. If we produce only one perturbation at $`t_w`$ ($`\mathrm{\Delta }t_1=1`$) the system remains in its original cycle, yielding a return map similar to that of Fig. 1. Then, under the same conditions we repeat the simulations adopting now $`\mathrm{\Delta }t_1=1000`$ and $`\mathrm{\Delta }t_2=10000`$. We observe in Fig. 6 that the return map has changed slightly when compared to Fig. 1. In particular, it no longer shows $`C_{tot}=1`$ in the time series, but remains periodic, signaling that the system has shifted to a different cycle due to the perturbations. Apparently the periods observed in Figs 1 and 6 are the same. The distribution of periods and transient times are currently under investigation, results will be published elsewhere. Fig. 6 remains the same by increasing $`\mathrm{\Delta }t_2`$, which guarantees that the differences with respect to Fig. 1 are not a transient effect.

According to these results, the role the noise plays, if allowed enough time to perturb the system significantly, is to drive the system from one attractor to a nearby one, which suggests that there is a family of periodic attractors which can be very close to one another in phase space. These effects, however, take place on a time scale which is longer \[$`𝒪(10^310^4)`$ time steps\] than the lifetime of the mice \[$`𝒪(10^210^3)`$ time steps\]. The small perturbations are therefore of little importance to the CA dynamics in the biologically relevant time scale, as suggested in previous work .

The behavior of the autocorrelation functions in Figs. 3 and 5 is reminiscent of what is observed in glassy systems. Here the CA approaches a periodic attractor, being thereafter deflected by the small perturbations. In glassy systems, temperature drives the system from one local minimum to another, preventing it from getting stuck in local minima of the potential energy surface. As the physical age of the system (characterized by the value of the waiting time $`t_w`$) grows, it finds itself exploring deeper and deeper regions of a rough potential energy surface, diffusing towards equilibrium . Due to the roughness of the potential energy, equilibrium is only attained on very long time scales and relaxation is very slow. The older the system is the more it gets confined to a restricted region of phase space and the time scales for relaxation get longer and longer. Eventually, if we wait long enough, the system equilibrates and the dynamics becomes stationary, losing sensitivity to the waiting time. This picture of an aging physical system is reminiscent to the loss of plasticity for adaptation in a living organism as it gets older. In either case, one observes a strong dependence on the waiting time $`t_w`$, evident when measuring two-time quantities like autocorrelation functions and responses. Results for the CA model are shown in Fig. 7, where we see the decay of autocorrelations for three systems with different ages or waiting times. Note that the horizontal axis is the time difference between the total time and the waiting time. The three curves should collapse in the case that the dynamical evolution is stationary.

### 3.3 Large immunizations

What is the role of the large perturbations on the dynamics of the system? From previous studies we would expect the large perturbations to accelerate the aging process: while the random small perturbations would change the route to the natural attractor of the system, the large ones would reduce the transient time, a conjecture that would explain the loss of plasticity of the older mice . The protocol adopted is the one described in Section 2, with six-layer perturbations every 35 time steps, always at the same sites.

In Fig. 8 we compare the results obtained for the unperturbed system and those of the system subjected only to large immunizations starting at $`t_w`$ (note that the densities are now also plotted as functions of $`tt_w`$). Somewhat surprisingly, the decrease of the correlations for large perturbations is small, when compared to the case of the small perturbations (compare with Fig. 3). In hindsight, however, this can be understood because the large perturbations are always produced in the same region. For very long times the correlations will attain a stationary regime whose active populations contain at least part of the region involved in the immunizations . In other words, the driving produced by the large periodic perturbations is of a completely different nature than that of the small perturbations. The large perturbations seem to play a selective role: the cycles that the system can reach are restricted to those that contain at least part of the populations incorporated during the immunization. According to this picture, the aging effects observed in Ref.  could be simply related to the exploration of the phase space: the older the system the less possibilities of choosing new cycles it would have. In Section 4 we discuss this issue further.

In Fig. 9 we compare the time evolution of the densities $`B_1`$ and $`B_2`$ and the autocorrelation functions for the system subjected only to large perturbations and for the system with both small perturbations and large immunizations. The results, as expected, confirm the dominance of the small perturbations, over the large ones, in driving the system faster to a different attractor. The increase of the densities around $`t10^4`$ for $`L=50`$ corresponds to the same finite size effects occurring in Fig. 3, being associated to the active populations which have been incorporated into the network by means of the immunization protocol.

## 4 Concluding Remarks

We have analyzed the dynamics of a cellular automata model for the B-cell repertoire which has reproduced experimental results of immunizations on mice. Our analyses provide a broader context in which memory and plasticity take place, as discussed in Ref. .

Defining quantities analogous to the autocorrelation functions used in glassy systems, we have shown that the biologically relevant phenomena takes place in the transient regime of the model. The dependence of these functions on both $`t`$ and $`t_w`$ (as opposed to the difference $`tt_w`$ only) is reminiscent of the “aging” behavior observed in glassy systems, despite the fact that the underlying dynamics of both systems are controlled by completely different mechanisms.

Starting from different initial conditions the deterministic dynamics takes the system towards a family of long-period attractors. When subjected to random small perturbations (Fig 4), the system is driven towards a new attractor of the family, revealing that most of the noise is assimilated. However, only part of the large perturbations is incorporated, due to the mechanisms of activation and suppression, leading to a saturation of the learning process .

From the biological point of view, the history of the mouse (sample) will be written by the different antigen presentations (random perturbations), starting from its “birth” (initial condition). Since the system is large but finite there is a maximum amount of information it may incorporate. The closer it is to its “destiny” attractor, the less information it is able to learn, since the deeper it already is in a given basin of attraction. Therefore the biological aging corroborated by the experimental results may be simply a consequence of this dynamical feature.

From the results obtained up to now there are evidences that the purely deterministic dynamics is therefore non-ergodic. Despite their large number, however, the family of periodic attractors occupy only a fraction of the phase space. Evidences of this compression in phase space has been obtained in the following computer experiment: selecting randomly two initial configurations (with the same initial concentration), we measured the Hamming distance between them as a function of time as both systems evolve without any perturbation. In Fig. 10 we show the evolution of the average HD for $`100`$ pairs (error bars are standard deviations). During the first 10 time steps the HD decreases very rapidly, but for $`t>100`$ we observe a quasi-stationary regime, reflecting the slow driving to the attractors (see inset). Note that this behavior is somewhat similar to that of the autocorrelation functions for small $`t_w`$.

Still focusing on the evidence of phase space compression, in Fig. 11 we present the distribution of HD sampled from $`500(5001)/2`$ pairs, for two consecutive time steps ($`t=4000`$ and $`t=4001`$). The distributions can be well described by a Gaussian, with a width that remains approximately constant (even for long times). Notice that the average value for $`t=4001`$ is slightly larger than for $`t=4000`$, reflecting the oscillatory behavior of the average HD as the pair of samples reaches their periodic attractors. It should be noted that, for large $`N`$, the Central Limit Theorem assures that randomly chosen initial configurations (with the same $`x`$) naturally give rise to a Gaussian distribution of the HD between any two of them, at time zero. Interestingly, the CA dynamics does not change the shape of the distribution, its only effect being to essentially shift the Gaussian towards lower mean values, on a long time scale. The spatio-temporal structure of the cycles, as well as their transients and basins of attraction, would be the object of further study and will be published elsewhere.

This work was partially supported by CNPq, Faperj, Finep, Capes and Projeto Enxoval-UFPE. MC and RMZS acknowledge IF-UFF, where part of this work was done during their stay there.


format: text

text:
# Optical Trapping and Manipulation of Nano-objects with an Apertureless Probe

## Abstract

We propose a novel way to trap and manipulate nano-objects above a dielectric substrate using an apertureless near-field probe. A combination of evanescent illumination and light scattering at the probe apex is used to shape the optical field into a localized, three dimensional optical trap. We use the coupled-dipole method and the Maxwell stress tensor to provide a self-consistent description of the optical force, including retardation and the influence of the substrate. We show that small objects can be selectively captured and manipulated under realistic conditions.

Since the seminal work of Ashkin on radiation pressure, the possibility to exploit the mechanical action of optical fields to trap and manipulate neutral particles has spawned a wide range of applications. From atomic and nonlinear physics to biology, optical forces have provided a convenient way to control the dynamics of small particles (see Ref.\[REFERENCES\] for a review). Optical tweezers, for example, have proved useful not only for trapping particles, but also for assembling objects ranging from microspheres to biological cells . However, most of these manipulations involve particles whose size is between one and several micrometers. While for much smaller particles, such as atoms or molecules, the scanning tunneling microscope provides a powerful tool for manipulation and engineering , dealing with neutral particles of a few nanometers requires new experimental approaches.

The idea of using a metallic probe to trap small particles was reported by Novotny et al . Their calculations showed that strong field enhancement from light scattering at a gold tip could generate a trapping potential deep enough to overcome Brownian motion and capture a nanometric particle in water (a related work by Okamoto and Kawata demonstrates theoretically the trapping of a nanometric sphere in water near an aperture probe ). This technique should be delicate to implement in practice for three reasons. First, before the particle can be captured, Brownian fluctuations will have a disruptive effect. Second, radiation pressure from the illuminating laser will impart momentum to the particle. Therefore, one would have to capture a moving object. Finally, it will be rather difficult to use a near-field probe to find in water a particle a few nanometers in size. One might wait for a particle to wander in the trapping region, but such an operating mode does not allow for selective capture.

In this letter, we propose an experimental scheme to selectively capture and manipulate nanoparticles in vacuum or air above a substrate, using the tungsten probe of an apertureless near-field microscope. The particles are not in a liquid environment, hence there is no Brownian motion and the apertureless probe can be used as a near-field optical probe to localize and select the particles . This is an important asset when different particles have to be placed according to a specific pattern or when interactions between particles are investigated. For example, with our configuration quantum-dot nanocrystals (CdSe/ZnS nanospheres a few nanometers in diameter) could be placed in a specific geometry to study dot-dot interactions, or one nanocrystal could be isolated to study single-dot properties .

We first consider a spherical particle with radius 10 nm placed in air on a substrate (Fig. 1). Unless it is stated otherwise the permittivity of both the particle and the substrate is $`ϵ=2.25`$. The particle is illuminated (wavelength 500 nm) by two counter-propagating evanescent waves created by the total internal reflection of plane waves at the substrate/air interface (angle $`\theta >\theta _c=41.8^{}`$ with $`\sqrt{ϵ}\mathrm{sin}\theta _c=1`$). These two waves have the same polarization and a random phase relation. This symmetric illumination ensures the lateral force vanishes when the sphere is just below the tip, thus avoiding that the sphere be pushed away from the tip.

We study the interactions between the particle and a tungsten probe (commonly used in apertureless microscopy) with a radius at the apex of 10 nm. The theory used to compute the optical forces has been presented in detail elsewhere . Here we give a succinct summary. We use the coupled-dipole method  to model the light scattering and find the electromagnetic field inside the tip and the particle. Note that this procedure takes into account the multiple scattering between the sphere, the tungsten tip, and the substrate. We then use the Maxwell stress tensor technique to derive the optical forces . We emphasize that the stress tensor approach is exact and does not rely on any assumption regarding the nature of the field (whether evanescent or propagating) or of the objects.

All forces are computed for an irradiance of 0.05W/$`\mu `$m2 (this corresponds, for example, to a 5W laser beam focused over an area of 100$`\mu `$m2). Figure 2 shows the $`z`$ component of the force experienced by the sphere versus the vertical position of the tip above the sphere, for both TE and TM polarizations and for two angles of illumination ($`\theta =43^{}`$ and $`\theta =50^{}`$). As the tip moves closer to the sphere, the evolution of the force for the two polarizations is radically different. For TM illumination there is a large enhancement of the field near the apex of the probe because of the discontinuity of the $`z`$ component of the electric field across the air/tungsten boundary . This enhancement is responsible for the force being positive at short distances. Note that for the sphere, the force of gravity is on the order of $`10^7`$pN. The $`z`$ component of the force experienced by the sphere (when the tip is in contact with the sphere) is about $`3`$pN which is $`10^8`$ times its weight. Hence, gravity can be neglected. Fig. 2a shows that for TM polarization the force is larger for $`\theta =43^{}`$ than for $`\theta =50^{}`$. This is related to the slower decay of the evanescent wave for the smaller angle, which results in a weaker coupling between the sphere and the substrate. Moreover, the slower the decay of the evanescent field, the larger the field that reaches the tip, and the larger the (positive) gradient force caused by the field enhancement at the apex of the probe. As a consequence, when the tip approaches the sphere, the sign reversal (negative to positive) of the $`z`$ component of the force occurs at a larger distance for $`\theta =43^{}`$ ($`z`$=21 nm) than for $`\theta =50^{}`$ ($`z`$=13 nm). This is shown in the inset of Fig. 2a. On the other hand, for TE polarization (Fig. 2b), the magnitude of the $`z`$ component of the force increases while the force remains negative (directed toward the substrate and away from the tip) as the tip gets closer to the particle. This prevents any trapping. Our calculations show that this behavior is caused by the decrease of the field inside the upper part of the sphere which, in turn, causes a decrease of the gradient force. Another way of explaining this is to note that because the apex of the tip and the sphere are small compared to the wavelength, they can be approximated by two dipoles. For TE polarization these dipoles are essentially parallel to the substrate and as shown in Ref.\[REFERENCES\], two parallel dipoles tend to repel each other. Since the magnitude of the dipole is stronger for $`\theta =43^{}`$ (due to the larger intensity of the field), the repulsive force is also stronger. Note that for a silver tip at the plasmon resonance frequency, the force acting on the sphere is positive for both polarizations. At the frequency considered here, tungsten behaves like an absorbing dielectric and the force is positive only for TM polarization.

To fully assess the probe-particle coupling we study the evolution of the force experienced by the particle as the probe is moved laterally, the sphere remaining fixed. The components of the force acting on the sphere are plotted in Figs. 3 and 4 for two distances between the tip and the substrate (20 nm and 31 nm), and for two angles of illumination ($`43^{}`$ and $`50^{}`$). For TM polarization (Fig. 3a), the $`z`$ component of the force is negative when the tip is far from the particle. When the tip gets closer, the particle starts experiencing a positive force along $`z`$. The change of sign occurs between $`|x|=`$30 nm for $`\theta =43^{}`$ and $`z=20`$ nm, and $`|x|=`$7 nm for $`\theta =50^{}`$ and $`z=`$31 nm. For TE polarization the force is negative (Fig. 3b). Similarly, we plot in Fig. 4 the lateral force. The symmetry of the force plot is a consequence of the symmetric illumination. For TM polarization (Fig. 4a) this force tends to push the particle toward the tip. For example, if the tip is located at $`x=`$10 nm, we can see in Fig. 4a that the $`x`$ component of the force experienced by the sphere is positive, hence the lateral force pushes the sphere toward the tip. If the tip is located at $`x=`$-10 nm, the $`x`$ component of the force is negative and, again, it pushes the sphere toward the tip. Therefore, when the tip and the particle are close enough for the $`z`$ component of the force to be positive and to lift the particle off the surface, the lateral force actually helps bringing the particle in the trap. Again TE polarization gives a different result. The lateral force always pushes the particle away from the tip. However, since the magnitude of the (downward) $`z`$ component of the force is larger than the $`x`$ component by a factor of 5, we expect that the sphere is not displaced when the tip is scanned over it under TE illumination. Note that apertureless probes are often used in tapping mode when imaging a surface. This mode minimizes the lateral motion imparted to the object by the optical force.

We have shown that a tungsten probe can be used to trap efficiently a nanometric object above a surface using TM illumination. For nanomanipulation it is important to assess the stability of the trap as the probe lifts the particle off the substrate. Figure 5 shows the $`z`$ component of the force when the sphere is fixed at the apex of the tip and the tip is moved vertically (solid and dashed curves). The optical force remains larger than the weight of the particle over a range of several tens of nanometers. The particle can therefore be manipulated vertically as well as horizontally. Note that the evolution of the force with the distance to the substrate is linear rather than exponential. The particle experience a negative gradient force due the exponential decay of the intensity of the illumination. At the same time, the particle suffers a positive gradient force due to the field enhancement at the tip apex, which also decreases exponentially with $`z`$ because this enhancement depends on the intensity of the evanescent illuminating light. The competition between these two contributions results in an weak decrease of the trapping force as the particle is moved away from the substrate. If we change the nature, size or shape of the object the magnitude of the force changes but the conclusions are qualitatively the same (Fig. 5, curves with symbols).

It is fundamental to know whether our procedure also works if several particles are clustered together. We consider a set of three spheres (radius 10 nm, permittivity 2.25) aligned along $`x`$. The probe is placed above the middle sphere. We account for the multiple scattering between the three spheres, the substrate and the tip, therefore the optical binding experienced by the spheres  is included in our description. For this configuration, TE illumination again does not lead to trapping. For TM illumination, we plot in Fig. 6 the $`z`$ component of the force experienced by the middle sphere and by those on the sides as a function of the vertical distance between the probe and the middle sphere. For an angle of incidence $`\theta =43^{}`$, the $`z`$ component of the force, although largest for the middle sphere, is positive for all the spheres. This could be a problem if one wanted to manipulate only one particle. The central particle can be selectively trapped by increasing the angle of incidence of the illuminating beams to tighten the trap in the $`x`$ direction (Fig. 4). In Fig. 6 we see that for $`\theta =50^{}`$ the optical force induced by the probe is positive only for the middle sphere. This remains true for three spheres aligned along $`y`$. Figure 7 shows the extraction of the middle sphere by the tip. Our calculation shows that the vertical force experienced by the two side spheres remains negative when the probe moves away from the substrate. Therefore, the spheres on the sides do not hinder the capture of the middle sphere.

In Fig. 7, the curve with circles shows that the selective trapping works as well when the middle sphere is twice as big as the others. Once the chosen particle has been trapped, it can be moved above the substrate. Our calculations show that the presence of other particles on the substrate does not destroy the trap, provided that during the manipulation the trapped sphere is kept at least 5 nm above the spheres that are on the substrate. We have checked that if the optical binding causes the side spheres to move laterally toward the middle sphere, increasing the angle of illumination still creates a negative (downward) force on the two side spheres, and the selective capture of the middle sphere is not hampered.

The procedure to trap a small object with a tungsten tip is therefore the following: first TE illumination is used while the tip scans the surface in tapping mode or in constant-height mode if the area under investigation is small enough. Such modes avoid the displacement of the particle by the tip. Once an object has been selected, the probe is placed above the object and the polarization of the illumination is rotated to TM. The probe is then brought down over the particle and captures it. The probe can then move the particle above the substrate, both horizontally and vertically, to a new position where it can be released by switching back to TE polarization (note that if, for some reason, one wishes to move the particle over distances larger than the size of the illumination spot, one could move *the sample* with piezoelectric device once the sphere is trapped at the apex of the tip). The lack of trapping under TE illumination is actually an important advantage during the imaging (selection) and release phases of the manipulation. Indeed, under TE illumination, when the tip is above a particle it actually increases the downward optical force, which contributes to prevent the tip from sweeping the particle away.

In conclusion, we propose a new method to trap and manipulate nanometric particles in air above a substrate, using an apertureless tungsten probe. The probe is used to scatter two counter-propagating evanescent waves, generating an optical trap at the apex. We showed that an object of a few nanometers, can be selectively trapped and manipulated. An interesting extension of this work will be a study of the influence of different illuminations (e.g. focused beam) and a systematic study of the influence of the nature of the particle and the tip on the trapping mechanism. For example, the strong spectral dependence of the electromagnetic response of metal particles (or resonance excitation of both dielectric and metallic particles) could lead to a material selective trapping.

P. C. C. and A. R. would like to thank Stéphanie Emonin for many fruitful discussions.


text:
# Variable polarization in the optical afterglow of GRB 021004Based on observations made with the Nordic Optical Telescope; based on observations collected at the European Southern Observatory, Chile, by GRACE (Gamma-Ray Burst Afterglow Collaboration at ESO), under programme 70.D-0523(A)

## 1 Introduction

The generally accepted source of gamma-ray burst (GRB) afterglow emission is synchrotron radiation, produced when the initial relativistic blast wave hits the circumburst matter and starts radiating (Rees & Mészáros 1992; Paczyński & Rhoads 1993; Mészáros & Rees 1997; Wijers et al. 1997; Wijers & Galama 1999) . Synchrotron radiation is highly polarized, up to 75% (Rybicki & Lightman 1979), and polarization has indeed been measured for 6 GRB afterglows (Wijers et al. 1999; Covino et al. 1999; Rol et al. 2000; Björnsson et al. 2002; Covino et al. 2002a, e; Bersier et al. 2003a; Masetti et al. 2003). See also the reviews by (Björnsson 2003) and (Covino et al. 2003a). These measurements and obtained upper limits (Hjorth et al. 1999; Covino et al. 2002d) show that the level of polarization is generally small, presumably because the intrinsically high polarization is averaged out to the few percent observed (Gruzinov & Waxman 1999; Gruzinov 1999; Medvedev & Loeb 1999).

If the outflow of the blast wave is collimated into a jet (Rhoads 1997, 1999; Sari et al. 1999), several models predict changes in the degree of linear polarization from a few up to 30% (Sari 1999; Ghisellini & Lazzati 1999; Rossi et al. 2002). So far, only hints for these variations have been seen (Rol et al. 2000), mainly because of the low polarization values that are measured and the difficulties involved in obtaining a time series of accurate polarization measurements of GRB afterglows. One exception is possibly GRB 020405, for which Bersier et al. (2003a) find a large variation in the degree of linear polarization within a short time interval (see also Covino et al. 2003b), which cannot be reconciled with any current model.

GRB 021004 was localized with the wide-field X-ray Monitor (WXM) on board the High-Energy Transient Explorer-II (HETE-II) with an initial positional error of 10$`\mathrm{}`$. The position was immediately issued to the community, which allowed the rapid discovery of its afterglow with the Oschin/NEAT robotic telescope (Fox 2002).

The afterglow light curve is well covered and shows some deviations from a standard power-law decay, for which various explanations have been offered, such as variations in the burst energy or variations in the density of the surrounding medium (see for example Lazzati et al. 2002; Heyl & Perna 2003; Nakar et al. 2003; Dado et al. 2003). Holland et al. (2003) measure a break in the light curve between 3.5 and 7 days after the burst. The redshift of the afterglow plus host galaxy was determined to be $`z=2.33`$ (see for example Møller et al. 2002), while the spectrum consists of a complex of absorption systems (Salamanca et al. 2002; Mirabal et al. 2002; Møller et al. 2002) at the redshift of the host.

Polarization measurements were obtained by various groups (Covino et al. 2002b; Wang et al. 2003; Rol et al. 2002). Here, we report on early polarimetric observations obtained by our group with the Nordic Optical Telescope (NOT) at the Observatorio del Roque de los Muchachos on the Canary Islands, and later with the Very Large Telescope (VLT) at the European Southern Observatory (ESO) in Chile.

## 2 Data reduction and analysis

Polarimetric observations of the afterglow of GRB 021004 at the NOT were performed from October 4.859 UT until 4.908 UT ($``$8 to 10 hours after the burst) with the Andalucia Faint Object Spectrograph and Camera (ALFOSC), using two calcite plates at different orientations and a Bessel R filter. The calcite plate yields two overlapping images of the field-of-view (FOV) of the telescope, separated by about 15 arcsec. One image allows for the measurement of the ordinary ray, the other image gives the extra-ordinary ray, which together allow for the measurement of one of the linear Stokes parameters. For each observation, the orientation of the calcite plate was either 45 degrees or 90 degrees. We obtained 3 pairs of exposures, each exposure with an integration time of 600 seconds, which allows the determination of Stokes $`Q`$ (0/90 degrees polarization orientation) and $`U`$ (45/135 degrees) for each pair.

The polarimetry observations with the FOcal Reducer/low dispersion Spectrograph (FORS 1) at the VLT Antu were performed a few hours later, from October 5.151 UT until 5.196 UT. We used a Wollaston prism with a rotatable half-wave plate at four different angles. The images do not overlap but are separated by a mask covering half the FOV. Each angle allows a measurement of both the ordinary and extra-ordinary ray, which in turn allows the determination of the Stokes $`Q`$ and $`U`$ parameters. The broad-band filter applied here was Bessel V. The observations consisted of three sets of four exposures, with an integration time of 120 seconds for each of the first four exposures, and an integration time of 300 seconds for each exposure in the last two sets.

All data were reduced using the IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. software suite. The images were first bias subtracted. The flatfielding of the ALFOSC images was performed with the orientation of the calcite plate for the flatfield identical to that of the target image. The FORS 1 images were flatfielded using flatfields without the Wollaston prism and half-wave plate. Artefacts introduced by the prism or half-wave plate can be corrected for using two observations, one with the half-wave plate oriented at either 0 or 22.5 degrees, and one with the half-wave plate oriented at 45 degrees difference.

The ALFOSC images show filter scratches which could not be completely corrected for by flatfielding. To allow for the detection of errors introduced by such artefacts, the position of the afterglow on the CCD was slightly offset from the centre on the third set of observations, while it was centered for the first two sets. Observations of standard stars verified that a source positioned in the centre gives the correct polarization.

To measure the flux, we used aperture photometry on both the ALFOSC and FORS 1 images, where the apertures sizes were adapted to the measured seeing. Due to the small FOV of ALFOSC in polarimetric mode, only a few stars could be measured, and it was not possible to derive an accurate point-spread function (PSF) for the two images and perform PSF photometry.

For the determination of the Stokes parameters, only the relative flux of the ordinary and extra-ordinary ray is required. In the ALFOSC case, however, the calcite plate projects images for the ordinary and extra-ordinary ray differently, which results in the PSFs being different for the two images, and, for small apertures, the incorrect ratio of light to be measured. To obtain the most accurate as well as precise result, we chose the smallest aperture size for which the resultant Stokes parameters are still consistent with those measured using larger apertures. An aperture with a radius equal to $`10`$ pixels, or 1.5 times the seeing, provided the best results. This aperture is somewhat larger than the one used by Björnsson et al. (2002), who performed polarimetric observations of GRB 010222 with ALFOSC, but their observations were performed with better seeing ($`0\stackrel{}{.}9`$ compared to $`1\stackrel{}{.}2`$ for our observations).

Calibration and verification of the procedures was done with both zero polarization standard stars, to check for any instrumental polarization, and high polarization stars. For ALFOSC we used BD +28°4211 and BD +32°3739 as zero polarization stars, and HD 204827 as high polarization standard star. For FORS 1, BD $`12`$°5133 was used as high polarization standard star.

We performed aperture photometry on the FORS 1 data and verified the results by PSF photometry. We used an aperture of 1.5 times the seeing. As there is no difference for FORS 1 in the PSF between the images of the ordinary and the extra-ordinary ray, this resulted in modest errors. We used the other stars in the FORS 1 FOV to calculate the polarization introduced by interstellar matter, assuming that the net polarization of the field stars is zero and that the largest fraction of interstellar matter is between these field stars and the observer. Any resultant polarization is then caused by the interstellar matter, and was found to be $`P=(0.58\pm 0.33)\%,\theta =(105\pm 16)\mathrm{°}`$ Since the spread in the ($`Q`$$`U`$) values of the field stars is rather large, as can be seen in Figure 1, we used this spread as the error in the above values, instead of the error in the weighted mean of the field stars, which is much smaller.

After calculation of the Stokes parameters, we corrected for polarization induced by the ISM. The degree of linear polarization and its position angle were calculated from both ISM-corrected and uncorrected Stokes parameters, after which the polarization degree was corrected following Wardle & Kronberg (1974) and Simmons & Stewart (1985) for bias resulting from the fact that $`P`$ is a definite positive quantity.

## 3 Results

For both ALFOSC and FORS 1, we obtained 3 sets of polarization measurements. The results from FORS 1 are constant from one set to another, showing no variations on this short time scale (15 to 16 hours after the burst). We therefore used the summed image to obtain a higher signal to noise ratio and one final result. For ALFOSC, the three sets do not show evidence for such constancy, and we calculated the results separately (see Table 1). Since the polarization during the ALFOSC observation could have been variable, the assumption of constant polarization, needed to calculate $`P`$ and $`\theta `$ from the separate measurements with the calcite plate at 45 and 90 degrees, is not valid. To see whether the results in Table 1 are still representative of the polarization at the times of observation, we paired the observations differently, obtaining two more sets. The resulting values for $`P`$ and $`\theta `$ are consistent with the first two sets of the original values, that is, they are intermediate values. The results for the third set are not consistent with the previous two results and either show a very rapid change in the polarization, or are due to an artefact in the data. The latter could result from to the aforementioned defects in the filter, caused by the afterglow positioned at a faulty position on the CCD (the first two measurements do not suffer from this, as outlined above).

When comparing the ALFOSC data with the FORS 1 data, we assume there is no significant difference in the polarization due to the different filters (V and R) we have used. Wang et al. (2003) mention a wavelength dependent change in the polarization for their spectropolarimetric measurements, but this only occurs below $`405`$ nm and would not affect our comparison.

We have plotted the resultant degree of linear polarization and position angle as a function of time in Figure 2, where we have also included the results by Covino et al. (2002b, c) and Wang et al. (2003). From the data uncorrected for ISM polarization, we see rapid changes in the degree of polarization between 8 and 10 hours after the burst. The polarization measured by FORS 1 could be entirely due to the ISM polarization, since the spread in the ISM polarization is rather large, as also remarked by Covino et al. (2002c) and Rol et al. (2002). However, the change in the polarization from our FORS 1 data point to the one measured by Covino et al. (2002c) shows that the polarization is at least partially intrinsic to the afterglow.

A change in the position angle and degree of linear polarization is entirely due to the afterglow, assuming that the polarization of the field stars and ISM is constant in time. The first two ALFOSC measurements are consistent with having a constant position angle; the later FORS 1 point shows a change by about 45° in the position angle at a 5 sigma level. Inclusion of the measurement by Covino et al. (2002c) even indicates a change of about 90° from 9 to 89 hours after the burst. For clarity, we have plotted the $`Q`$ and $`U`$ Stokes parameters for all the data together with the Stokes parameters for several field stars in Figure 1.

To estimate the significance of this change, we calculated the probability that the measured values originated from one constant value, by calculating the $`\chi ^2`$ values for $`Q`$ and $`U`$, $`\chi _x^2=\mathrm{\Sigma }_i\left(\frac{x_i\overline{x}}{\sigma _{x_i}}\right)^2`$, where $`x`$ = $`Q`$ or $`U`$, and $`\overline{x}`$ is the weighted mean of the 7 available measurements for the corresponding Stokes parameter. Applying an F-test to the resultant $`\chi ^2`$ gives the requested probability. We also applied this procedure with $`\overline{x}`$ being the average ISM polarization value, and $`\sigma _x`$ the corresponding spread therein. All probabilities are small. There is a 1.2% chance that $`U`$ can entirely be attributed to ISM polarization, but this probability for $`Q`$ is almost zero, as is then the probability for $`P`$. The probability that either $`U`$ or $`Q`$ belongs to one average value is less than $`10^6`$.

## 4 Discussion

Several models explain changes in the polarization angle in the context of jetted outflow of the gamma-ray burst ejecta. The break seen in the light curve of GRB 021004 (Holland et al. 2003; Bersier et al. 2003b) is indicative of such a jetted outflow, though the many bumps in the light curve hamper the detection of such a break. Holland et al. (2003) also exclude a spectral change due to the passing of the cooling frequency $`\nu _c`$ through the optical as a cause of this break.

We now proceed under the assumption that the change in the polarization is caused by a jetted outflow.

The model proposed by Rossi et al. (2002) explains changes in the polarization using a structured jet: a jet with brighter (and possibly faster) core surrounded by dimmer (and slower) wings, with a standard energy reservoir. However, this model does not predict a change in the position angle and is thus ruled out by our measurements for this burst. The model suggested by Sari (1999) and Ghisellini & Lazzati (1999) predicts a 90 degrees change in the position angle, roughly around the time of the break. The precise moment of this change, as well as the maximum observable polarization, depends on the ratio of $`\theta /\theta _0`$, where $`\theta `$ is the angle between the jet axis and the line of sight, and $`\theta _0`$ is the initial opening angle of the jet. From our measurements, we deduce that the break time, where $`\theta _01/\mathrm{\Gamma }`$ ($`\mathrm{\Gamma }`$ being the bulk Lorentz factor), is somewhere between 10 hours and 1 day after the burst. The dependence of the change in polarization on $`t`$ is via $`\mathrm{\Gamma }`$. If $`\mathrm{\Gamma }(t)`$ is different than assumed in Sari (1999) and Ghisellini & Lazzati (1999), which is not unlikely in view of the complex behaviour of the light curve and its explanations, the time of the break will be different than estimated above. Color changes have also been seen in the optical part of the energy distribution of the burst (Bersier et al. 2003b; Matheson et al. 2003). Such behaviour requires detailed models for the polarization, whereas we have here used the general model for a smooth afterglow behaviour. It is likely, however, that these color changes are of little influence to the models described above, since they were seen past one day, while the largest change in our data is before one day.

Holland et al. (2003) find a jet break time of $`t_b=6`$ days after the burst, from fitting a broken power law to the data . These estimates are in stark contrast with our findings. However, Holland et al. (2003) find that the break is gradual and occurred over a period of 3.5 to 7 days after the burst; they also note that their estimate of the break time might be too high, possibly putting the break around 2 days. With this latter value, our data would agree more with the jet model for polarization. The various bumps in the light curve might further obscure the detection of an early time jet break.

The smooth and gradual break in the light curve would mean that the ratio $`\theta /\theta _0`$ is high and we are viewing the jet close to the edge (see Sari 1999; Ghisellini & Lazzati 1999). However, this should then give rise to significantly higher values in the degree of linear polarization than measured, unless we observed close to the moment where the position angle changed by 90 degrees. That could then also explain the third set of ALFOSC observations, and would show that those observations where taken very close to the jet break time.

## 5 Conclusions

Our polarimetric observations clearly show a change in the polarization of GRB 021004, most distinct in the position angle. The latter changes by 45° between 9 and 14 hours, and the inclusion of a later data point by Covino et al. (2002c) indicates even a 90° change over a 3.5 day period. Within the currently proposed GRB jet models, this would mean that we are looking at a uniform jet, with a break time of $`t_b1`$ days after the burst. This is in contrast with the result obtained by Holland et al. (2003), who obtained $`t_b6`$ days, but with a large spread in this value ($`3.5`$ to $`7`$ days), which could still be reconciled with our findings. The structured jet model as proposed by Rossi et al. (2002) is ruled out by the fact that this model does not predict a change in the polarization angle.

###### Acknowledgements.

ER acknowledges support from NWO grant nr. 614-51-003. JPUF gratefully acknowledges support from the Carlsberg Foundation. This work was supported by the Danish Natural Science Research Council (SNF). JMCC acknowledges the receipt of a FPI doctoral fellowship from Spain’s Ministerio de Ciencia y Tecnología. The data presented here have been taken using ALFOSC, which is owned by the Instituto de Astrofísica de Andalucía (IAA) and operated at the NOT under agreement between the IAA and the NBIfAFG of the Astronomical Observatory of Copenhagen. The NOT is operated on the island of La Palma jointly by Denmark, Finland, Iceland, Norway, and Sweden, in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. The FORS 1 data were obtained as part of an ESO Service Mode run for ToO programme 70.D-0523(A). The authors acknowledge benefits from collaboration within the Research Training Network ”Gamma-Ray Bursts: An Enigma and a Tool”, funded by the EU under contract HPRN-CT-2002-00294.


text:
# Non-exponential Dissipation in a Lossy Elastodynamic Billiard, Comparison with Porter-Thomas and Random Matrix Predictions

## Abstract

We study the dissipation of diffuse ultrasonic energy in a reverberant body coupled to a waveguide, an analog for a mesoscopic electron in a quantum dot. A simple model predicts a Porter-Thomas like distribution of level widths and corresponding nonexponential dissipation, a behavior largely confirmed by measurements. For the case of fully open channels, however, measurements deviate from this model to a statistically significant degree. A random matrix supersymmetric calculation is found to accurately model the observed behaviors at all coupling strengths.

A random matrix model of wave scattering from an ideal lead, off a chaotic region, has been employed in nuclear theory for decades \[1-6\]. It originated in the theory of nuclear reactions and is widely used in the analysis of compound nuclei, mesoscopic quantum dots and microwave cavities \[1-6\]. It has provided a model of intrinsic loss mechanisms in classical wave chaotic systems, with implications for microwave cavities and reverberant ultrasonics. In particular it has been shown that dissipation is not necessarily exponential in time, and made specific predictions for how that decay should behave.

A simple argument that relates such behavior to the distribution of resonance widths leads to a non-exponential decay law under the assumption of a Porter-Thomas (chi-square) width distribution key-5 ; key-15 . That distribution follows from first order perturbation theory on a system with Gaussian mode shape statistics and a finite number of loss channels. A more general argument may be constructed by the information-theoretic approach key-2  or the random Hamiltonian approach key-4 ; key-8 . These also predict non-exponential decay and corresponding delay time distributions, but with different details key-4 ; key-5 ; key-11 . Acoustic key-15  and microwave key-14  measurements have confirmed non-exponential decays. Critical comparisons with the models have not been undertaken.

In this Letter we address diffuse energy decay in an open acoustical system. While non-exponential decays have long been observed, and modeled by means of Porter-Thomas distributions of resonance widths, deviations from those predictions have not yet been detected. Indeed, it has not yet been clear whether measurements can distinguish between the predictions of the simple argument and those of random matrix theory (RMT) \[3–9\]. It is therefore interesting to compare both the RMT predictions key-8  and the simple argument predictions key-15  with measurements conducted on an experimental realization of a chaotic billiard attached to a waveguide. In most cases we find that a Porter-Thomas model does an adequate job of fitting observed decay profiles. In the case of strongly coupled channels, however, a full RMT prediction is superior. The difference is small, but statistically significant.

Consider a reverberant body as in Fig. 1 coupled to a waveguide. The body has mean diffuse energy density (per mode) of $`\epsilon `$. The total energy in a band of width $`\mathrm{\Delta }\omega `$ is $`\epsilon D\mathrm{\Delta }\omega `$, where the body’s modal density $`D=N_{Body}/\omega `$. An attached waveguide has mean energy per outgoing mode no greater than $`\epsilon `$. Incoming modes have no energy. Each guided mode (i.e. each channel) therefore carries power at a rate no greater than $`(d\mathrm{\Delta }\omega )(\epsilon /2)v_g`$, where $`d`$ is the lineal modal density per length in the channel, $`d=N_{channel}/\omega L=1/\pi v_g`$ and $`v_g`$ is the group velocity of that mode. The factor ($`1/2`$) is due to neglect of the incoming modes. Thus each channel conducts outgoing power $`\mathrm{\Pi }(1/2\pi )\epsilon \mathrm{\Delta }\omega `$, independent of dispersion in the channel. At maximal coupling each channel contributes the same mean partial width, an energy decay rate of $`\mathrm{\Pi }/\epsilon D\mathrm{\Delta }\omega =1/[2\pi N_{Body}/\omega ]1/t_{Heisenberg}`$. This picture does not apply to individual normal modes of the body, but only to the mean. Those normal modes of the body which overlap well with the waveguide will dissipate rapidly; those which overlap poorly will do so slowly. While the average level width (and early time decay rate) may well correspond with the above estimate, the less strongly dissipated modes will eventually dominate a transient decaying field; the apparent dissipation rate will appear to diminish.

Modeling the diffuse field as a superposition of independent real normal modes with Gaussian statistics (this is not correct, lossy systems generically have complex eigenmodes), each with a proper decay rate given by first order perturbation theory, leads to a chisquare-like distribution of modal decay rates, and a net transient energy decay

$$E(t)=E_0\underset{i=1}{\overset{M}{}}\left(1+2\sigma _it\right)^{1/2}$$
(1)
where $`M`$ is the number of open outgoing channels, and $`\sigma _i`$ is the decay rate through the $`i`$th channel, $`\sigma _i1/t_H`$. In certain limits, it reduces $`E(t)E_0\mathrm{exp}\left(_i\sigma _it\right)`$.

In the special case of equipotent channels, $`\sigma _i=\sigma /M`$, one recovers the simplest Porter-Thomas model key-15 ,

$$E(t)=E_0\left(1+2\sigma t/M\right)^{M/2}$$
(2)
This form for the dissipation of ultrasonic energy density has been confirmed phenomenologically by taking $`E_0`$, $`M`$, and $`\sigma `$ to be adjustable parameters. Observed profiles $`E(t)`$ have been found to fit remarkably well key-15 . A particularly noteworthy case is that of key-14  in which the case of three weakly coupled ($`\sigma _i1/t_H`$) channels was successfully fit to Eq. (1). In many of these fits it has not been clear whether the recovered parameter $`M`$ is meaningful, i.e. whether it does indeed correspond to a discrete number of equipotent effective loss channels.

A better theory, beyond first order perturbation, is provided by considering a random matrix model key-3 ; key-4  for which level width distributions are known to be non-chisquare key-4 . Here the dynamics is governed by a (GOE) Gaussian random Hamiltonian, consistent with the assumed chaotic ray trajectories of the body, plus an anti-Hermitian part corresponding to a discrete number $`M`$ of outgoing channels. Supersymmetric techniques key-4 ; key-8  allow construction of various averages. Amongst the more easily constructed is $`E(t)=\left|G_{ij}\left(t\right)\right|^2`$, the mean square response at a site $`j`$ distinct from the site $`i`$ of the source, each distinct from any dissipative sites. $`G_{ij}\left(t\right)`$ is the time-domain Green’s function. On inverse Fourier transforming the result of a supersymmetric calculation key-8  one finds:

$`E(\tau ){\displaystyle _1^{\mathrm{}}}{\displaystyle _1^{\mathrm{}}}𝑑\lambda _1𝑑\lambda _2\mathrm{\Pi }(\tau ,\lambda _1,\lambda _2)f(\tau ,\lambda _1,\lambda _2)`$
$`\times {\displaystyle \frac{\theta \left(\lambda _1\lambda _22\tau +1\right)\theta \left(2\tau \lambda _1\lambda _2+1\right)\left(1\left(2\tau +\lambda _1\lambda _2\right)^2\right)}{\left(\lambda _1^2+\lambda _2^2+\left(2\tau +\lambda _1\lambda _2\right)^22\lambda _1\lambda _2\left(2\tau +\lambda _1\lambda _2\right)1\right)^2}}`$ (3)
where, $`\tau =t/t_H`$, $`M`$ is the number of channels, each characterized by a coupling parameter $`g_i1`$, $`f=\left(\lambda _1^21\right)\lambda _2^2+\left(\lambda _2^21\right)\lambda _1^2+1\left(2\tau +\lambda _1\lambda _2\right)^2`$, $`\mathrm{\Pi }=_{i=1}^M\left(g_i+2\tau +\lambda _1\lambda _2\right)\left(g_i^2+2g_i\lambda _1\lambda _2+\lambda _1^2+\lambda _2^21\right)^{1/2}`$. For a large number of weak ($`g_i1`$) channels, this reduces to: $`E(t)=E_0\mathrm{exp}\left(\tau 2/\left(g_i+1\right)\right)`$, in agreement with the naive model. At finite $`M`$ the two models differ slightly. A comparison at $`M=4`$, $`g_i=1.5`$, $`\sigma _i=2/\left(g_i+1\right)t_H=4/5t_H`$ for all $`i`$, is shown in Fig. 2. The naive model overestimates curvature.

We study the aluminum body (volume $`561cm^3`$, free surface $`451cm^2`$) pictured in Fig. 3. Non parallel faces and defocusing surfaces enhance ray-chaos. After preliminary baseline measurements, it was welded to an aluminum wire ($`1100`$ alloy, $`3.18mm`$ diameter, $`3m`$ length). Tests were carried out with the spiral part submerged in a water bath. Attenuation in the water assures negligible reflected energy and thus the presence of only outgoing waves. This was confirmed by separate measurements. The guided elastic waves of a circular waveguide are described by the Pochammer dispersion relation key-16 . All modes with azimuthal number $`n>0`$ are two-fold degenerate. At low frequency, below the first cutoff at $`580kHz`$, there are $`M=4`$ propagating guided waves. Two are flexural ($`n=1`$), one is extensional ($`n=0`$) and one torsional ($`n=0`$). A calculation of cutoff frequencies gives the number of open channels. The strength with which these channels are coupled is not known apriori. Measurement of reflections of waves incident from the wire onto the block indicate that the coupling is good; reflections are generally weak. Dispersion and the possibility of mode conversion complicate any attempt to be more quantitative.

Profiles $`E(t)`$ were constructed both before and after the waveguide was attached. For each case, piezoelectric pulses of negligible duration were applied to pin-transducers in light oil contact with the body, as in Fig. 3. Waveforms of durations of up to $`100msec`$, and bandwidths to $`2MHz`$ were recorded at a digitization rate of $`5MSa/sec`$. A low pass filter with a cutoff of $`2.25MHz`$ prevented aliasing. Repetition averaging improved signal to noise ratios and extended the system’s dynamic range. The resulting waveforms were time-windowed into $`62`$ successive $`1.64msec`$ sections with tapered edges. Each windowed waveform was Fourier transformed and squared and integrated over rectangular bins of width $`25`$ or $`50kHz`$. The result was an array of spectral energy densities versus time for each of several narrow frequency bands. On repeating this for $`16`$ distinct source and receiver positions an average $`E(t)`$ was constructed for each band. A typical profile is seen in Fig. 4. As expected, the waveguide has augmented the decay rate.

The upper curve is the reference case without the waveguide; it shows a decay which is very nearly exponential, i.e. consistent with intrinsic decay mechanisms being widely distributed and corresponding to a large number of weakly coupled dissipative channels. As in key-15 , the behavior fits well to Eq. (2); chisquares are excellent. While the parameter M extracted from that fit is perhaps not meaningful, we do take Eq.(2) as a valid way to smooth the reference data. The smoothed reference $`\mathrm{ln}E(t)`$ is then subtracted from the measured $`\mathrm{ln}E(t)`$ in the waveguide-attached case to give a difference $`\mathrm{ln}E`$, what we would have measured in the attached case if our reference block had been non-dissipative. This profile shows substantial curvature, consistent with a hypothesis of a small number of outgoing channels.

The observed difference $`\mathrm{ln}E`$’s were then themselves fit to Eq. (2). Figs. 5 and 6 show the parameters $`\sigma `$ and $`M`$ extracted from those fits, and compares them with expectations based on assuming perfect coupling: $`M_{Pochammer}=`$ number of propagating modes; $`\sigma =M_{Pochammer}/t_H`$. That the $`\sigma `$ of the difference $`\mathrm{ln}E`$ is not in excess of this theory, even at high frequencies, is an indication that the welding process has not significantly increased intrinsic absorption. The correspondence is remarkably good for such a simple theory; $`M`$ and $`\sigma `$ show the predicted features, in particular those associated with the onsets of new guided modes at $`580`$, $`800`$ and $`1000kHz`$. At low frequencies, the correspondence is poor; the model assumption of perfect coupling is incorrect.

The reduced chisquares of these fits are shown in the inset. At high frequencies they are within one standard deviation of the expected value of unity. Plots of the residuals show that fluctuations are not systematic. We conclude that the higher-frequency data is consistent with Eq. (2); small curvatures are well modeled by a single phenomenological parameter $`M_{effective}`$; the high frequency data does not permit conclusions in regard to the relative virtues of the various theories. At low frequency the chisquares exceed unity. This is in part due to the inadequacy of our spatial averaging there; different transducer positions are often within a wavelength and are correlated. It is in part also due to some systematic deviations that indicate a need for a better theory; unphysical values of $`M`$ (Fig. 6) are a further indication of that need.

At low frequencies ($`<580kHz`$), where there are only four open channels, we attempt a fit to the richer theories ($`1`$) and ($`3`$) by fixing $`M`$ at four and adjusting the coupling strengths. We take the two flexural waves to have identical coupling (the weld is axisymmetric); set $`g_3=g_4`$, $`\sigma _3=\sigma _4`$, and adjust only $`E_0`$ and three values of $`g_i1`$ (or $`\sigma _i1/t_H`$). The fits’ values for total loss coefficient ($`t_H\sigma _i4`$ and $`2/\left(1+g_i\right)4`$) are shown in Fig. 7. Representative plots of the data and fits are shown in Fig. 8. In cases with weak coupling, both models do well. At higher frequencies where coupling is efficient (all $`g`$’s close to unity), and first order perturbation theory is invalid, there are significant discrepancies. At several of the frequencies only the full RMT calculation Eq. (3) is statistically acceptable. At a few other frequencies where Eq. (1) gives an acceptably low chisquare, its call for all four channels to be perfectly coupled, $`\sigma _i=1/t_H`$, is implausible. The fits at these several frequencies may be the first demonstration of non-Porter-Thomas like level width distributions; it is also evidence for the superior applicability of a full RMT calculation over that of simpler theories. This has implications for RMT and wave chaos, for mesoscopics, nuclear physics, microwave physics, diffuse field ultrasonics, and also for structural acoustics where modeling of losses, both intrinsic and radiative key-17 ; key-18 , is gathering increased attention.

This work was supported by NSF grant CMS-0201346.