text:
# Coupled map gas: structure formation and dynamics of interacting motile elements with internal dynamics

## 1 Introduction

High-dimensional dynamical systems with many interacting nonlinear elements have been studied extensively. Here we call such a system “coupled dynamical system”. Among such systems, both those with local coupling and global coupling have been studied at length. For instance, in the study of pattern formation and spatiotemporal chaos, partial differential equations and coupled map lattices (CML)  have been employed extensively. In a CML, dynamical elements that are spatially distributed interact with each other through local interactions such as diffusive coupling. The pattern dynamics that arise through the interplay between the internal dynamics of individual elements and their interaction is studied by changing the coupling strength. Globally coupled systems with all-to-all interactions among dynamical elements also represent a typical type of coupled dynamical system . They provide the simplest type of model for complex dynamic networks. Using globally coupled systems, several important concepts in the study of high-dimensional dynamical systems have been developed. Such globally coupled dynamical systems can be considered idealizations of many types of physical, chemical and biological systems.

Coupled dynamical systems cannot always be categorized as either locally coupled or globally coupled. In particular, living systems consist of a huge number of dynamical elements with a variety of scales and with a variety of types of interactions, ranging from local to global. They exhibit interesting collective behavior at the macroscopic level, whose spatiotemporal scales can be much larger than that of individual elements. In such complex systems, even if the properties of each element and the mechanism of its interaction with other elements are well understood, the collective properties of the system are generally much too complicated to be understood without extensive investigation.

In most models of coupled dynamical systems that have been employed to this time, the interaction structure itself is fixed in the sense that the elements with which any given element interacts as well as the nature of the interaction are predetermined. In particular, interacting elements are never decoupled, and non-interacting pairs are never coupled.

In some biological and other types of complex systems, however, the interactions among the dynamical elements can also change in time. For instance, consider a system of motile organisms (e.g. bacteria). Each organism in this system has internal dynamics and interacts with the others. The interaction structure here changes in time according to the motion of the elements. Even within a cell, interactions among molecules change in time. There, enzymatic reactions consist of a finite number of molecules, each of which is a dynamical element having internal states . By sharing a reactive chemical resource, these molecules interact so as to constitute a reaction network. The interactions among the molecules should be time dependent, due to the motion of the molecules. Another example is provided by neural networks, where the coupling between elements changes in time, through the formation and destruction of synaptic connections between neurons.

In order to describe systems of these kinds, we need to construct model coupled dynamical systems whose couplings are dynamic, changing in a manner that depends on the states of individual elements and the interactions between them. In such models, we may find a novel class of collective behavior and self-organization. We thus propose to study the characteristic dynamics of systems in which the interactions between dynamical elements change in time.

There are many possible types of models with internal dynamics, interactions, and motile elements that exhibit interaction structure. Instead of constructing a realistic model specific to some phenomena, we introduce a simple abstract model possessing the above three features. In the present paper, we consider a coupled map whose elements are characterized by spatial positions that change in time. By setting a finite range for the interaction of these motile elements, the existence of an interaction between two given elements will change in time. In this system, we classify several phases with distinct pattern dynamics that depend on the parameters governing the coupling strength, the coupling range, and the internal dynamics. The characteristic dynamics of each phase are elucidated.

This paper is organized as follows. In Section 2, we introduce a simple model of interacting motile dynamical elements constituting a combination and extension of a couple map lattice and a globally coupled map. The model exhibits a variety of phenomena, such as spatial clustering, fusion and fission of the clusters, and intermittent diffusion of the elements. In Section 3, these phenomena are studied for models in one- and two-dimensional space. The global phase diagram is presented in Section 4. Spatial clustering is one of the characteristic dynamical phenomena in the present system. In Section 5, we study this clustering in detail for the case in which the clusters that form are not completely isolated in space but rather interact with each other. Then, in Section 6, we study the case in which each cluster is completely isolated after formation. In Section 7, we study the case in which clusters form structures in space with a scale larger than that of a single cluster. These structures are maintained by the interplay between the internal dynamics of elements and the interactions among elements. The dynamics of the formation of these structures are studied. The paper is concluded in Section 8 with the summary and discussion.

## 2 Couple map gas as a simple model of interacting motile dynamical elements

To this time, two types of coupled dynamical systems have been investigated, one with local coupling and the other global coupling. Such systems consist of elements with internal dynamics and interactions among them. In the former case, the interactions among elements are restricted locally in space, as in the case of reaction-diffusion equations and coupled map lattices (CML). In particular CMLs with local chaotic dynamics and nearest neighbor interactions have been extensively studied as models of spatiotemporal chaos. The simplest CML model of this kinds is given by

$$x_{n+1}(i)=(1\epsilon )f(x_n(i))+\frac{\epsilon }{2}\left(f(x_n(i+1))+f(x_n(i1))\right),$$
(1)
where $`n`$ is a discrete time step, $`i`$ is the site index, and $`f(x)=1ax^2`$. In the globally coupled case, each element interacts with all other elements. The simplest such model is a globally coupled map (GCM), given by

$$x_{n+1}(i)=(1\epsilon )f(x_n(i))+\frac{\epsilon }{N}\underset{j=1}{\overset{N}{}}f(x_n(j)),$$
(2)
where $`n`$ is a discrete time step, $`i`$ is the element index, and $`N`$ is the total number of elements. This model can be regarded as a mean-field version of a CML, or as an infinite dimensional CML. GCMs have also been studied as the simplest models of complex dynamical networks.

In most studies of coupled systems carried out to this time, elements with internal (oscillatory) dynamics are fixed in space, and the structure of the connections between elements is fixed in time; i.e., two elements that are coupled initially are always coupled. The coupling strengths between elements are also constant in most previously studied models, although time dependent coupling strengths have been introduced in some models , with connection topologies fixed.

In the present paper, we study a system in which the connections between elements are time dependent in accordance with the motivation expressed in Section 1. For this purpose, we consider the introduction of “motile” elements into CML and GCM models as a natural extension of these systems. Two such elements interact when they are located within a given range in space, and thus the structure of the couplings between elements changes in time due to the motion of the elements. Because the position of elements changes in time, two elements coupled at some time can be decoupled at some other time.

Here, we are interested in general aspects of coupled dynamical systems of motile elements, and for this reason we wish to consider as simple a model as possible. We impose the following conditions on this model:

1. Each element has time-dependent internal state and spatial position.
2. Each element is active in the sense that, even in isolation, its internal state exhibits oscillatory (chaotic) dynamics.
3. The dynamics of the internal state of a given element are affected by local interactions between this and other elements.
4. The elements move through the influence of the force produced by the local interaction, which depends on the internal states of the interacting pair.

In order to construct a model satisfying the above conditions, we stipulate that each element has an internal state represented by a scalar value and is characterized by a time-dependent position in real space. For the dynamics of local internal state, a one-dimensional map that can exhibit chaotic oscillation is employed. The dynamics of the internal state of a given element are influenced by other elements within a distance $`R`$ from this element. This part of the model is basically given in the form of coupled maps. Furthermore, there is a force between elements within the same range, which depends on the internal states of the pair of interacting elements. The motion of the elements in space is governed by these forces. We call the class of models defined in this manner “coupled map gases (CMGs)”.

A simple example of a CMG is given by

$$\begin{array}{ccc}x_{n+1}(i)\hfill & =\hfill & (1\epsilon )f(x_n(i))+\frac{\epsilon }{N_n(i)}\underset{j𝒩(i)}{}f(x_n(j)),\hfill \\ \stackrel{}{r}_{n+1}(i)\hfill & =\hfill & \stackrel{}{r}_n(i)+\underset{j𝒩(i)}{}\frac{\stackrel{}{r}_n(j)\stackrel{}{r}_n(i)}{|\stackrel{}{r}_n(j)\stackrel{}{r}_n(i)|}(x_{n+1}(i),x_{n+1}(j)),\hfill \end{array}$$
(3)
where $`x_n(i)`$ is the internal state of the $`i`$-th element at time step $`n`$, and $`\stackrel{}{r}_n(i)`$ is its position in $`d`$-dimensional space. As the map governing the internal dynamics, we adopt the logistic map $`f(x)=1ax^2`$, because CMLs and GCMs employing this map have been extensively studied as prototype models of spatiotemporal chaos and network dynamics. In the above equations, $`|\stackrel{}{r}(j)\stackrel{}{r}(i)|`$ is the distance between the $`i`$-th and $`j`$-th elements, the set of elements interacting with the $`i`$-th element is denoted by $`𝒩(i)=\{j:|\stackrel{}{r}_n(j)\stackrel{}{r}_n(i)|R\}`$, and the number of elements in $`𝒩(i)`$ is denoted by $`N_n(i)`$. Here, $`(x_{n+1}(i),x_{n+1}(j))`$ is a “force”, which depends on the internal states $`x`$ of the $`i`$-th and $`j`$-th elements. In the present paper, we adopt $`(x(i),x(j))=Fx(i)x(j)`$ as the form of the force, using a fixed parameter $`F`$. The elements, whose number is fixed at $`N`$, are located within the fixed spatial domain $`[0,L]^d`$, for which we use periodic boundary conditions.

There are four basic control parameters in our model, $`a`$, $`\epsilon `$, $`N(\frac{R}{L})^d`$ and $`\frac{F}{R}`$. The role of each control parameter is as follows: $`a`$ determines the nonlinearity of the local dynamics, $`\epsilon `$ is the coupling strength between elements, $`N(\frac{R}{L})^d`$ is the density of elements within the interaction range, and $`\frac{F}{R}`$ is the effective motility of the elements. Throughout this paper, we fix $`\frac{NR}{L}10`$, and $`\frac{F}{R}0.01`$. The effective motility is chosen to be relatively small so that the motion of each element in space is rather smooth.

Note that the elements in the present model are “active” in the sense that, even in isolation, their internal dynamics can consist of periodic or chaotic oscillation. This strongly contrasts with the situation in recently studied models of the collective motion of flocks  and ants , in which passive elements are employed.

In the present CMG, the motion in real space is “overdamped”, i.e., in the equation describing this motion, the force term is essentially equated with the velocity, without an acceleration term included. It is rather straightforward to introduce a momentum variable for each element and with it include an acceleration term in the equation of motion. Indeed we have carried out several simulations, with such equations. However, we have found that as long as there is friction force, the important results are unchanged. For this reason, we consider only the overdamped case here.

The choice of the force in Eq.(3) is rather arbitrary. We have imposed the condition that the force between two elements is attractive or repulsive, depending on the signs of the internal states of the interacting elements. With this stipulation, the results we obtained seem to be rather universal for any choice of $`(x(i),x(j))`$. As an alternative type of force, we have also studied the case $`(x(i),x(j))=\mathrm{cos}[2\pi (x(i)x(j))]`$, and we found that the qualitative features are the same.

When the parameter $`F`$ is positive, two elements interact attractively if the signs of their internal states are the same, and repulsively if they are opposite. Since the variable $`x(i)`$ oscillates within the range $`[1,1]`$, the sign of the force is correlated with the degree of synchronization of the interacting elements’ oscillations. It should be noted, however, that interaction between two arbitrary elements is attractive on the average if $`F`$ is positive, because the average of $`x`$ for a single logistic map (and in the present CMG) is positive for $`a<2.0`$.

As a type of coupled map system, CMG is a natural extension of a CML in which the condition of the confinement of an element to a single lattice point is removed. If the elements are spontaneously located at positions separated by a distance about $`R`$, the dynamics of their internal states are those of a CML. On the other hand, for a set of elements located within a distance $`R`$and separated from other elements, the dynamics of their internal states are those of a GCM. Thus the CMG spontaneously chooses the coupling structure of their internal states ranging from those of a CML to those of a GCM.

In the case of CML or GCM, the variable $`x(i)`$ of each element can be regarded as a field variable. In the present CMG, in which elements move in space, internal states cannot be regarded as field variables on fixed space. The motion of elements in space can be regarded as dynamics external to the dynamics of the internal state of each element. We are concerned with the “interplay” between such external dynamics and the internal dynamics.

It is also interesting to consider our model from a different viewpoint. Our model (in the 1-dimensional case) can be rewritten as

$$r_{n+1}(i)=r_n(i)+fx_{n+1}(i)\left(\underset{j𝒩_l(i)}{}x_{n+1}(j)\underset{j𝒩_r(i)}{}x_{n+1}(j)\right)$$
(4)
where $`𝒩_l(i)=\{j;j𝒩(i),r(j)<r(i)\}`$ and $`𝒩_r(i)=\{j;j𝒩(i),r(j)>r(i)\}`$. In this form, the second term on the right hand side is regarded as the gradient of the field at the position $`r(i)`$, which is determined by the internal state variables $`x(i)`$ of the elements located within a distance $`R`$ of $`r(i)`$. Particles move in this chaotic field, which they themselves form. This interpretation of our model is possible for any spatial dimension.

## 3 Characteristic phenomena observed in a coupled map gas

In this section, we give a rough survey of the phenomena observed in the present model, leaving detailed analysis to later sections. Spatiotemporal diagrams of the 1-dimensional system are given in Figure 1, where the trajectories of the elements in real space are plotted by drawing lines between the successive positions of the elements.

To describe the dynamics, it is useful to introduce the notion of a “cluster” in space. Here, a “cluster” is defined as a set of elements that are located in close proximity in such a manner that, as a whole, they are distinguishable from the other elements.111The term “cluster” is used with a different meaning in the case of GCMs. In such models, a “cluster” is defined as a set of synchronized elements With regard to the dynamics of clusters, we note the following three distinct types of behavior:

1. Elements gather to form a cluster whose members are not fixed; i.e., elements enter and leave the cluster over time. Such clusters often split into two, and two clusters often merge into one. This formation and division of clusters occurs repeatedly in time, as shown in Figures 1(a), (b) and (c).
2. Elements form clusters that remain separated by distances approximately equal to $`R`$ from their neighboring clusters. Accordingly elements in each cluster can interact with those in neighboring clusters. In this case, neighboring clusters exchange elements intermittently as shown in Figure 1(d).
3. Clusters are all separated by distances larger than $`R`$. Hence, there are no interactions between elements in different clusters. In this case, the members of each cluster are fixed. Here, the size of each cluster is smaller than $`R`$, so that the dynamics within each cluster are described by a single GCM, as shown in Figures 1(e) and (f).

This classification of the spatiotemporal dynamics displayed by our model is also valid in the 2-dimensional case. Examples of successive instantaneous patterns in the 2-dimensional case are shown in Figure 2. In a 2-dimensional system, there appears a variety of cluster shapes, in addition to those described above. In this situation described by Figure 2, elements start to gather to form a circular cluster. Then, the cluster divides into two new circular clusters. Note that the time evolution begins from random initial conditions. The figure displays only patterns appearing after transient behavior has died away, and after the cluster has first formed. The formation and division of circular clusters seen here is generally observed in the 2-dimensional system.

Another example of the 2-dimensional case is depicted in Figure 3. In this example, the value of $`F`$ is negative. Nevertheless, a lattice structure of clusters is formed. Within a cluster, the forces between elements are repulsive, while the forces between elements of neighboring clusters are attractive. As far as we have studied, the formation of clusters is generic in this class of systems, not restricted to the present CMG.

Although the pattern formation in the 2-dimensional case is interesting, the behavior we have found for the 2-dimensional system to this time can conceptually be understood from the results of the 1-dimensional case. For this reason, we focus on the 1-dimensional case in the remainder of the paper.

## 4 Global phase diagram

In this section we characterize the behavior discussed in Section 3 by introducing some statistical measures, and also present the phase diagram of our model. As the measures, we introduce the following three quantifiers, corresponding to the motion in real space, the dynamics of the interaction structure, and the coherence of the oscillatory internal dynamics of the elements.

1. Diffusion of elements: In the simulations we carried out, elements can either move throughout the space or be localized in space, depending on the value of the parameters. In the former case, diffusive, rather than ballistic, motion is observed. Therefore it is useful to introduce a quantity to measure the diffusion of elements. For this purpose we first define the mean square displacement (MSD) of the positions of elements as

$$\delta r_t^2=\frac{1}{N}\underset{i=1}{\overset{N}{}}\left(r_{t+t_0}(i)r_{t_0}(i)\right)^2,$$
(5)
where $``$ denotes the average over an ensemble of samples. When the MSD increases linearly with $`t`$, it is possible to define the diffusion coefficient of the elements, characterizing their diffusion in real space, as

$$D=\underset{t\mathrm{}}{lim}\frac{1}{t}\delta r_t^2.$$
(6)
2. Change of coupling: Due to the motion of the elements, a pair of elements coupled at one time step can be decoupled at the next time step. The frequency that this change in coupling takes place is a good measure for classifying phases with regard to the fluidity of the interaction. As this measure, we numerically study the conditional probability that a pair of coupled elements remains coupled at the next step (denoted by $`P_{CC}`$) and the conditional probability that a pair of decoupled elements remains decoupled at the next step (denoted by $`P_{DD}`$). In the present model, these probabilities are given by

$`P_{CC}=𝒫(|r_n(i)r_n(j)|R||r_{n1}(i)r_{n1}(j)|R)`$ (7)
$`P_{DD}=𝒫(|r_n(i)r_n(j)|>R||r_{n1}(i)r_{n1}(j)|>R)`$ (8)
for arbitrary $`i`$ and $`j`$, where $`𝒫(𝖠|𝖡)`$ is the conditional probability that $`𝖠`$ happens given that $`𝖡`$ happens. Obviously, the probabilities of the change from coupled to decoupled ($`P_{CD}`$) and vice versa ($`P_{DC}`$) are given by$`P_{CD}=1P_{CC}`$ and $`P_{DC}=1P_{DD}`$. If there are no coupling changes, $`P_{CC}`$ and $`P_{DD}`$ are unity. When there are frequent changes, $`P_{CC}`$ and $`P_{DD}`$ take small values.
3. Coherence among the dynamics of internal states: In globally coupled maps, (chaotic) oscillations of all elements are synchronized in the strong coupling regime. Similarly, in the present model, the internal states of elements located within the interaction range $`R`$ can be synchronized when tightly coupled. When the coupling of these elements changes, such synchronization is partially or completely destroyed. Also, weak coherence can be sustained over all elements, even when they are separated by distance greater $`R`$. Therefore, to characterize the system, it is useful to introduce a measure of the degree of coherence.

Four typical types of evolution of the MSD with time are plotted in Figure 4. The diffusion coefficient of the elements given in Eq.(6) calculated from the MSD is plotted as a function of the coupling strength $`\epsilon `$ in Figure 5, with $`a=1.8`$ fixed. From Figure 5, it can be seen that there are three regimes in which the diffusion coefficient is very small. These three regimes are separated by regimes with larger values of the diffusion coefficient. Since the diffusion coefficient plotted in Figure 5 is obtained from the MSD over a finite time ($`10^5`$ steps), it is not possible to conclude that the diffusion coefficient in these three regimes is zero, but we can reach a somewhat certain conclusion by investigating whether the MSD increases with time or eventually stops increasing. Examples of the time evolution of the MSD are displayed in Figure 4.

From the results given in Figures 4 and 5, the following three types of regimes are distinguished.

1. Strong diffusion: There are two such regimes ($`0.1\epsilon 0.13`$ and $`0.2\epsilon 0.24`$ in Figure 5), in which the time evolution of the MSD is clearly proportional to $`t`$ and the diffusion constant is large.
2. Weak diffusion: In the regime with intermediate coupling strength ($`\epsilon 0.17`$) among the three regimes with small diffusion coefficient, the MSD increases almost linearly with time. Here, the motion of the elements exhibits a diffusive behavior, whose diffusion constant is distinctly smaller than that in the regimes with strong diffusion.
3. No diffusion: For the other two regimes with small values of the diffusion coefficient in Figure 5, i.e., for $`\epsilon 0.25`$ and $`\epsilon 0.1`$, the MSD eventually stops increasing in time. In this case, the elements are localized in space.

Using the above classification based on the diffusion behavior, we can classify the types of cluster dynamics discussed in the last section. These three types of cluster dynamics are characterized as follows.

1. The strong diffusion regimes (a) correspond to case 1 in the last section (see Figures 1(a), (b) and (c)), where the formation and division of clusters continues indefinitely.
2. The weak diffusion regime (b) corresponds to case 2 of the last section (see Figure 1(d)), where clusters form a lattice, and elements are exchanged intermittently between clusters.
3. The no diffusion regime (c) corresponds to case 3 of the last section (see Figures 1(f) and (e)) where each element is localized in a cluster, and each cluster is isolated from the other clusters.

The probabilities $`P_{CC}`$ and $`P_{DD}`$ for the changes of couplings are plotted in Figure 6 as a functions of the coupling strength. (There, all other parameters are fixed with the same values as in Figures 4 and 5.) For $`\epsilon 0.1`$ and $`\epsilon 0.25`$, $`P_{CC}=P_{DD}=1`$, since in these regimes the coupling never changes in time. In the intermediate coupling regime, with a small diffusion coefficient, the probabilities are much smaller than 1. This indicates that changes of couplings take place frequently, even though the diffusion is much slower.

With respect to the three types of quantities discussed above, the MSD, the probabilities $`P_{CD}`$ and $`P_{DC}`$ and the coherence, we have succeeded in classifying the behavior of the present CMG into the following five phases: 1) coherent phase, 2) fluid phase, 3) intermittent phase, 4) desynchronized phase, 5) checkered phase. The phase diagram is given in Figure 7 and the characteristic features of each phase are summarized in Table 1

We now describe these phases in detail.

When the coupling strength is sufficiently large, elements do not diffuse in space, as shown in Figures 4(a) and 5. In this case, the coupling among elements does not change, as shown in Figure 6. The spatiotemporal diagram corresponding to this phase is given by Figure 1(f). Here, it is seen that elements split into disconnected clusters separated by distances larger than the interaction range $`R`$. Therefore, there are no interactions between clusters. Within each cluster, all elements are coupled to each other. Since all clusters are completely decoupled, the dynamics of the elements in each cluster are reduced to those of a globally coupled map whose total number of elements is the number of elements in the cluster. Within each cluster, the internal oscillations of elements are synchronized. Thus the coherence in this phase is that of the complete synchronization case of the GCM.

When the coupling strength is sufficiently smaller, elements diffuse rapidly, as shown in Figures 4(b) and 5. In this case, the couplings among elements change frequently, as shown in Figure 6. In this phase, the formation and division of clusters is repeated indefinitely, as shown in the spatiotemporal diagram given by Figures 1(a), (b) and (c). In this phase, the cluster structure is dynamic, and its evolution is determined by the interplay between the internal dynamics of interacting elements and the motion of the elements. We discuss these dynamics in more detail in Section 5.

In the midst of the fluid phase, there exists a region in which the diffusion coefficient is distinctively small but non-zero, as shown in Figures 4(c) and 5. Even though the diffusion here is quite weak, the coupling among elements change quite frequently as shown in Figure 6. In this case, the clusters form a lattice structure depicted in Figure 1(d), which is almost static with only the intermittent exchanges of elements.

When the coupling strength is even smaller, beyond a certain value, the elements no longer diffuse, as shown in Figures 4(c) and 5. In this case, the coupling among elements also does not change, as shown in Figure 6. Here, clusters are separated by distances greater than $`R`$ (as in Figure 1(e)) and therefore do not interact, while within each cluster all elements interact with each other, and the oscillations of elements are desynchronized. Because all clusters are completely separated, the internal dynamics here are those of the desynchronized state of the GCM.

There exists another phase in a region with a smaller value of the nonlinearity parameter $`a`$, in which the couplings among elements change frequently, but elements do not diffuse in space. In this case, two neighboring clusters are separated by a distance approximately equal to $`R`$. Some elements between neighboring clusters interact with each other at some time steps, but they can be separated by distances greater than $`R`$ at some other times. This is the reason that here the probabilities $`P_{CC}`$ and $`P_{DD}`$ are smaller than 1. Within each cluster, oscillations are strongly correlated, while oscillations of neighboring cluster are out of phase. In this phase, elements exhibit period-2 band oscillations. Denoting the state of $`x(i)`$ only by its sign, all elements in a given cluster oscillate as $`+++\mathrm{}`$, while those in a neighboring cluster oscillate perfectly out of phase as $`++\mathrm{}`$. Thus the states $`x`$, when represented in this way, form a checkered spatial pattern among cluster.

## 5 Formation and collapse of clusters in the fluid phase

In a coupled dynamical system, the dynamics of elements often become synchronized so that they come to form groups in phase space. Such groups are often also referred to as “clusters” . In the present case, such synchronization of oscillations influences the formation and collapse of clusters in real space. Conversely, the motion of elements in real space influences the internal dynamics of elements. As a result, formation, fusion and fission of the clusters occur repeatedly in the fluid phase. Studying the interrelation between the internal dynamics and the motion in real space, we have found that the process of cluster formation, fission and fusion in space can be described as follows.

The interactions among elements tend to synchronize the internal states of elements. However, chaotic nature of the internal dynamics tends to destroy this synchronization. When synchronization dominates, forces between elements on the average become attractive, and clusters are formed.

Within each cluster, the dynamics are of a GCM type (to the extent that the influence of neighboring clusters can be ignored). A number of synchronized groups are formed among the elements in a cluster.

In some of such groups the sign of the internal states evolve as $`++\mathrm{}`$, while that of in other groups evolve as $`++\mathrm{}`$. Then, the force between two elements from such two types of groups is repulsive, while the force between the elements within the same type of group is attractive.

Due to the repulsive forces between groups, the cluster divides into two new clusters, if the cumulative force is sufficiently strong. This repulsion continues until the distance between the two new clusters is larger than $`R`$. Consequently, the two new clusters come to be located at a distance of about $`R`$.

When the distance between the two new clusters becomes larger than $`R`$, they no longer interact. Initially, almost all elements within one of these clusters (i.e., elements within the range $`R`$) are nearly synchronized, but this state is unstable, because synchronized chaos in the GCM is unstable in this parameter region. Thus the process returns to step 2. Thus each cluster has the potentiality of repeated division.

If a cluster interacts with other clusters, depending on the internal dynamics of these clusters, these clusters can fuse together. Then, the evolution of this newly merged cluster returns to step 2.

In this way, the fusion and fission processes of clusters continues indefinitely.

In order to understand the above described process in a more quantifiable manner, we now characterize the behavior of the spatial motion and the internal degrees of freedom quantitatively.

In the fluid phase, one of the most important processes in the internal dynamics is the synchronization. In the synchronization process, the difference between the internal variables of two elements decreases with time. This process can be effectively seen by measuring the information creation per bit $`I_n^{(m)}`$, introduced in Ref., where $`n`$ is a discreate time step and $`m`$ is a particular integer value. For this, we partiton the space of the internal state into intervals of a size $`2^m`$. Then, we define the distribution function $`P_n^{(m)}(j),(j=0,\mathrm{},2^{m+1}1)`$ as the fraction of elements whose internal states take values in the interval $`[j\times 2^m1,(j+1)\times 2^m1](j=0,\mathrm{},2^{m+1}1)`$ at time step $`n`$. Here, we disregard the position of each element. Then, the information creation per bit $`I_n^{(m)}`$ at time step $`n`$ is defined by

$$I_n^{(m)}=s_n^{(m)}s_n^{(m1)},$$
(9)
where $`s_n^{(m)}`$ is the $`m`$-bit entropy at time step $`n`$ given by

$$s_n^{(m)}=\underset{j=0}{\overset{2\times 2^m1}{}}P_n^{(m)}(j)\mathrm{log}P_n^{(m)}(j).$$
(10)
Since the $`m`$-bit entropy is a nondecreasing function of $`m`$, $`I^{(m)}`$ cannot take a negative value. If two elements are synchronized to a precision of $`2^{\mathrm{}+1}`$, the value of $`s^{(m)}`$ increases at $`m=\mathrm{}`$. Therefore $`I^{(m)}`$ takes a positive value when $`m=\mathrm{}`$, while $`I^{(m)}`$ equals zero for values of $`m`$ less than $`\mathrm{}`$. If the two elements are becoming synchronized, the largest value $`m=\mathrm{}`$ at which $`s^{(m)}`$ increases and $`I^{(m)}`$ is non-zero increases with time. Thus, the synchronization process can be viewed as the increase of the largest value of $`m`$ for which $`I^{(m)}`$ is greater than zero.

At the bottom of Figure 8, $`I_n^{(m)}`$ for $`m=2,\mathrm{},30`$ is plotted as a function of time. Several synchronization processes can be seen in this figure, as the increase of the largest value of $`m`$ for which $`I_n^{(m)}>0`$. In Figure 8, the time evolution of the $`m`$-bit entropy $`s_n^{(m)}`$ is also plotted for $`m=10`$. This entropy measures the degree to which elements are synchronized. The synchronization processes can also be identified by the decrease of $`s_n^{(10)}`$, while once synchronization has been established, its breakdown can be identified by the increase of $`s_n^{(10)}`$.

Change in the spatial configuration of elements is characterized by the spatial entropy. First, we partition the entire space into intervals of a given size. Then we define the the fraction of elements in each interval $`j`$ as $`Q_n(j)`$. Here, the width of an interval is chosen to be approximately the same as the effective width of the cluster (see the caption of Figure 9). In our simulations, we chose the value $`L/64`$. The spatial entropy is computed as

$$𝒮_n=\underset{j=1}{}Q_n(j)\mathrm{log}Q_n(j).$$
(11)
From the time evolution of $`𝒮_n`$ displaying in Figure 8, one can see that $`𝒮_n`$ decreases in the process of cluster formation and increases in the process of cluster collapse.

In Figure 8, $`𝒮_n`$ is plotted as a function of time. As the synchronization proceeds, as is shown by the increase of the largest $`m`$ for which $`I^{(m)}>0`$, the decrease of the spatial entropy $`𝒮_n`$ is observed. This indicates cluster formation processes. After the synchronization process becomes completed, the spatial entropy $`𝒮_n`$ stops decreasing. Then, $`𝒮_n`$ begins to increase, indicating that the process of cluster collapse has begun. The collapse of cluster accompanies the breakdown of the synchronization among the internal dynamics of its element. As a result, the $`m`$-bit entropy $`s_n^{(10)}`$ also begins to increase. In Figure 8, the conditional probability $`P_{CD}+P_{DC}`$ is plotted as a function of time, indicating the frequency of coupling change (see Section 4). From this, it is seen that the collapse of a cluster leads to an increase in the frequency of coupling change.

In this way, the formation and collapse of clusters take place repeatedly due to the interplay between the internal dynamics and the spatial motion of the elements. Such processes are expected to take place at any point in space. Although $`I_n^{(m)}`$, $`s_n^{(m)}`$, $`P_{CD}+P_{DC}`$ and $`𝒮_n`$ are computed by treating all elements identically without taking account of spatial structures, we can clearly see the processes of cluster formation and collapse reflected in the behavior of these averaged quantities.

## 6 Absence of fusion and fission of clusters in the coherent and desynchronized phases

In this section, we study cases in which neither fusion nor fission of clusters takes place once clusters are formed. One such case is in the coherent phase with the coupling strength sufficiently large, and the other is in the desynchronized phase with the coupling strength sufficiently small.

In the case of the coherent phase, elements in a cluster interact with each other strongly, so that the internal dynamics of the elements in the same cluster are highly synchronized. It is the attraction among elements resulting from synchronization that causes a cluster to form. For a system in the synchronized state, beginning from random initial conditions, after transient behavior dies away, clusters are formed and exist isolated in space, with no interaction among them. Hence, the internal dynamics of the elements in any given cluster are those of a GCM in its coherent state. If coherent oscillation of these isolated GCM systems is stable, this stable cluster formation is guaranteed, and the resulting clusters will be separated by a distance greater than $`R`$.

In the case of the desynchronized phase, elements are coupled weakly, so that there is no synchronization, and the coherence among the internal dynamics of elements is weak. Even in such a case, the elements form clusters that are also separated by distances greater than $`R`$, and again there exists no interaction among clusters. In this case, within a cluster, the internal dynamics of the elements are those of the desynchronized state of a GCM. Then, the forces between elements change sign frequently, and the interactions between elements can often be repulsive. Thus the reason that no elements escape from a cluster is not immediately evident.

A possible reason for this fact may be that the sum of forces among the elements gives an attraction towards the center of the cluster, even if each two-body force is repulsive. In order to check if this could be the case, the direction and strength of the force acting on each element was calculated as a function of the distance from the center of mass of the cluster. The strength of the force was calculated as follows. The strength of the force acting on an element is given by the second term on the right hand side of Eq.(3). Then, the position of an element in the cluster is defined by its distance from the center of the cluster. Here, the center of a cluster is defined as the center of mass of the elements that are located within a distance $`d`$ of a given element at a given time step. Since even elements that directly interact with each other can be members of different clusters, $`d`$ should be smaller than $`R`$. The strength of the force on an element and its fluctuation as functions of the position can depend on the distance $`d`$. However, if a well-defined cluster exists, there should be a range of values of $`d`$ for which the value of the force calculated in this manner does not strongly depend on $`d`$. We take the value of $`d`$ within such an interval. The position of the center of a cluster computed in this way, by choosing an arbitrary element, is approximately the same for any chosen element within the cluster.

In Figures 9(a) and (b), the strength of the force experienced by an element and its fluctuation are plotted as functions of the position, and in Figure 9(c) and (d), the position of an element at one time step as a result of the force experienced at the previous time step is plotted as a function of the positions. In both the cases of the fluid phase, in (a) and (c), and the desynchronized phase in (b) and (d), elements around the center of the cluster are attracted to the center, on average. However, by taking the fluctuations into account, these two phases are clearly distinguished.

For the fluid phase, the fluctuation (indicated by the error bars in the figure) always crosses over the line $`y=x`$ in Figure 9(c). It is thus seen that elements can occasionally escape from the cluster. For the desynchronized phase, by contrast, the fluctuation does not cross over the line $`y=x`$ for positions within a particular range near when attractive force disappears. This indicates that in the desynchronized phase, the elements around the edge of the cluster are with certainly attracted to the center, and therefore elements cannot escape from the cluster. We thus understand that the collective attraction in the desynchronized phase is the source of cluster formation.

## 7 Structure formation and weak diffusion in the intermittent phase

In this section we study the intermittent phase by considering structure formation and intermittent behavior.

In the intermittent phase, after transient behavior has died out, through the repeated fusion and fission of clusters, clusters eventually form a lattice structure with a nearly fixed spacing $`R`$ (see Figure 1(d)). This lattice structure can be clearly characterized by the power spectrum defined as

$$P(k)=\frac{1}{N}\underset{j=1}{\overset{N}{}}x(j)e^{2\pi \mathrm{i}kr(j)},$$
(12)
in which the internal degrees of freedom are also taken into account. Figure 10 displays typical power spectra for the intermittent phase and the fluid phase. The sharp peaks in Figure 10(a) indicate the formation of lattice structure over a long spatial distance. The first peak is located at the wave number $`\frac{1}{2R}`$ not $`\frac{1}{R}`$, implying some structure with wavelength $`2R`$. This is because the internal degrees of freedom of the elements in a cluster are approximately in phase, whereas those of the elements of neighboring clusters are of nearly opposite phase. Hence the clusters form a checkered pattern, and the wavelength is $`2R`$. Even in the case of the fluid phase, depicted in Figure 10(b), there is a peak at around the wave number $`\frac{1}{2R}`$. This peak is not as sharp here as it is in the case of the intermittent phase, due to the successive formation and collapse of clusters.

In the fluid phase, synchronization among elements leads to the formation of clusters, as shown in Section 5, where the synchronization process was characterized by the information creation per bit, $`I_n^{(m)}`$. In the intermittent phase, the internal dynamics of elements within a cluster are roughly in phase. In contrast to the fluid phase, these dynamics do not exhibit any symptoms of synchronization among elements. In Figure 11, $`I_n^{(m)}`$ is plotted as a function of time. No increase of the largest value of $`m`$ for which $`I_n^{(m)}>0`$ is observed. This indicates that no synchronization process takes place.

If one cluster is isolated from the lattice structure, in phase motion within the cluster cannot be maintained and as a result, the cluster divides into two. Thus, in the intermittent phase, the interactions among clusters are essential for the formation of clusters. As shown in Figure 6, the high activity of interaction between clusters also suggests the importance of interactions.

In spite of the absence of synchronization among elements and the frequent change of the interactions among clusters, the diffusion constant in the intermittent phase is extremely small but not zero, as shown in Figure 5. This indicates that the structure is stable most of the time. Then, we are faced with the question of why this structure is so stable.

To answer this question, we recall the collective attraction phenomena in a cluster investigated in the previous section. In Section 6, it was seen that in the desynchronized phase, when we plot the direction and the strength of force acting on an element as a function of the position of the elements in the cluster, we find that there is a collective attraction of elements that acts to form a cluster. Hence, in this case, there exists no interaction among clusters and no diffusion of elements. We also calculated the force and its fluctuation in the present case. In Figure 12(a), the average force experienced by an element is plotted as a function of its position, and the position at one time step as a result of the force experienced at the previous time step is plotted in Figure 12(b). As shown, there is an average attractive force exerted an element around the center of the cluster, and there is a region in which the fluctuation of the force is not large enough to change its sign. The position at one time step, as a function of that of at the previous time step does not cross over the diagonal line $`y=x`$, except the center. This indicates the elements remain in a cluster for a very long time. We should note that in the present case, it is the interactions among the neighboring clusters that causes the structure to be stable. On the other hand, there are some neutral zones, where the force is nearly zero (and the next position is around the $`y=x`$ line) at positions larger than about 1.0. This leads to the intermittent switch of elements among neighboring clusters.

## 8 Summary and Discussion

In this paper we have proposed an abstract model, a coupled map gas, in order to study coupled dynamical systems whose couplings change in time in a manner that depends on the states of the system elements and their interactions. For this purpose, we introduced “motility” in space for the elements in our model, which represents a combination and extension of CML and GCM systems. In the present model, active elements with internal dynamics move through space and as determined by their interactions with the other elements. The internal state of each element is determined by these interactions as well as by its own intrinsic dynamics. As a result of their motion in space, the couplings among elements can change in time.

The present model exhibits a variety of phenomena. One of the most characteristic types of behavior of this system is the formation of clusters in space. These clusters are not isolated, but rather interact with each other. As a result of such interactions, clusters exchange elements. The formation and collapse of clusters occur repeatedly. Depending on the parameter values, clusters can form a lattice structure. In such a case, the exchange of elements between clusters is only intermittent.

The above mentioned phenomena were studied here from the viewpoint of the interplay between the internal and external dynamics. As studied in Section 7, the clusters form a stable lattice structure, in which the internal oscillations of the elements in neighboring clusters are out of phase with each other. Thus, there are two kinds of elements, whose internal states have one of the two different phases. However, this formation of lattice structure is not a simple pattern formation of the ordering of two kind of elements, as is indicated by the intermittent diffusion of elements. On the one hand, the internal dynamics of an element depend on its position in the structure. As a result, two kinds of elements with different types of internal dynamics emerge. On the other hand, the existence of two kinds of elements leads to the formation of lattice structure. Some balance between the effects of the dynamics of spatial degrees of freedom and the dynamics of internal degrees of freedom is essential to the formation of stable structure here.

As the coupling strength is increased or the nonlinear parameter of the elements is decreased, there is an increasing tendency toward the synchronization. In this case, elements in a cluster can form synchronized groups with oscillations of different phases. As a result, the lattice structure becomes destabilized, and the formation and collapse of clusters take place repeatedly. Conversely, as the coupling strength is decreased or the nonlinear parameter is increased, the dynamics of each element become too chaotic to allow the formation of groups with correlated phases of oscillation. In this case too, lattice structure is destabilized, and again clusters form and collapse repeatedly.

In this paper we have discussed a novel type of structure formation and dynamics in a system of interacting motile elements with internal dynamics. Often, this structure is not rigid and is flexible, due to the interplay between the internal dynamics and the interactions. Structure formation and collective dynamics of active elements are widely seen in a biological systems as well as in a biological-type physico-chemical systems. The presently studied coupled map gas system may be too simple and abstract to capture all the complexity of such biological systems. Still, due to the universality in the class of systems possessing interaction motile elements of the behavior we have investigated, we expect that the present study will sheds a new light on the understanding of collective dynamics in biological systems.

The authors would like to thank S. Sasa, T. Ikegami, M. Sano and A. Mikhailov for useful discussions. This research was supported by Grants-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan (11CE2006). T. S. gratefully acknowledges the support from the Alexander von Humboldt Foundation (Germany).


format: text

text:
# The Broadband Afterglow of GRB 980329

## 1 Introduction

The discovery of long-lived afterglow emission from gamma-ray bursts (GRBs) has revolutionized the field by enabling redshift determinations and host galaxy identifications. In addition, by interpreting the emission in the theoretical framework of a relativistic shock wave expanding into a circumburst medium, broadband afterglow measurements can constrain the explosion geometry and energetics, as well as the properties of the surrounding medium (Wijers & Galama, 1999; Granot, Piran & Sari, 1999; Harrison et al., 2001; Panaitescu & Kumar, 2001b). The hydrodynamic evolution of the blast-wave is strongly influenced by the the total energy in the shock, the geometry of the outflow and the density structure of the medium into which it is expanding. The time-dependence of the radiated emission from the shock depends on the hydrodynamic evolution, as well as on the partition of energy between the radiating electrons and the magnetic field. With data of sufficient quality in conjunction with a theoretical model, broadband measurements of the spectral evolution of the afterglow allow us to deduce fundamental physical parameters of the explosion.

Bright gamma-ray bursts such as GRB 980329 have been targets of extensive broadband followup. GRB 980329 was well-localized in the gamma-rays, and the position quickly refined as a result of the X-ray detection by In’t Zand et al. (1998). However, initial searches for an optical afterglow were unsuccessful until variable emission was identified at radio wavelengths (Taylor et al., 1998). Subsequent observations of the radio position uncovered a faint optical counterpart (Djorgovski et al., 1998), as well as a relatively bright near-infrared (NIR) transient (Klose, 1998). Because of the delay in the identification of the optical afterglow, the early optical monitoring was somewhat sparse. In spite of the eventual detection of optical and NIR counterparts, no redshift has been determined, due to the faintness of the host emission and its lack of prominent emission lines.

In this paper we present broadband observations of the afterglow of GRB 980329. We monitored the emission at multiple radio frequencies over times extending beyond the first month (early time data were reported by Taylor et al. (1998)). In addition, we present late-time optical and NIR observations in the R, I, H and K bands, as well as a re-analysis of submillimeter data reported by Smith et al. (1999). We fit the broadband emission to a fireball model, allowing us to derive the physical parameters of the afterglow, and measure the properties of the host galaxy.

## 2 Observations and Data Reduction

### 2.1 Optical/NIR Data

Optical and NIR observations were made using the Keck111The W. M. Keck Observatory is operated by the California Association for Research in Astronomy, a scientific partnership among California Institute of Technology, the University of California and the National Aeronautics and Space Administration. 10-m telescopes on Mauna Kea, Hawaii in the R, I, and K bands. Except for a re-calibration of the K-band points from Larkin et al. (1998a,b), all the photometric measurements presented in Table 1 have not been previously published. Deep R- and I-band images were obtained on several epochs using the Low Resolution Imaging Spectrometer (LRIS; Oke et al. (1995)), and on one epoch (2001 January) with the Echelle Spectrograph and Imager (ESI; Sheinis et al. (2000)). The data consist of multiple CCD exposures with typical durations of 300 s taken in good photometric conditions at an airmass between 1.1 to 1.3. We flat fielded and combined the data following standard practice. For calibration purposes, we used an image of the GRB 980329 field taken at the Palomar 60-inch in January 1999 under photometric conditions. Its photometric zeropoint was determined using four field stars which were in common with Reichart et al. (1999).

The K-band images of the GRB 980329 field were taken with the Near Infrared Camera (NIRC; Matthews & Soifer 1994) on the Keck I telescope. We used IRAF to reduce the data, and the DIMSUM package222http://iraf.noao.edu/iraf/ftp/contrib/dimsumV2/ to combine images and subtract sky background. We used observations of Persson et al. (1998) standards on each photometric night to calibrate a sequence of stars in the field, against which we performed relative photometry of the afterglow. We estimate that the calibration is accurate to approximately 5%.

We also made use of an H-band observation from the Hubble Space Telescope (HST) archive. HST observed GRB 980329 on 1998 October 16 through the F160W filter with NICMOS Camera 2 as part of GO7863 (PI: A. Fruchter). We used the standard NICMOS pipeline developed at STScI with the best available reference and calibration files to process the images. We then shifted and added the exposures to form a final, stacked image. We used the IRAF aperture photometry tools to measure the afterglow flux and the NICMOS photometric calibration published by Stephens et al. (2000) to convert this measurement to the standard H-band magnitude.

### 2.2 Very Large Array (VLA)

Observations were made with the VLA333The NRAO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. NRAO operates the VLA. at three frequencies, 8.46 GHz, 4.86 GHz and 1.43 GHz. All observations employed the standard VLA continuum mode, with data being recorded in two adjacent 50 MHz bandpasses. Calibration of the array phase was derived from observations of the nearby calibrators J0653+370 and J0713+438. Calibration of the flux density scale was done using J1331+305 or J0542+498, or by extrapolating past measurements of J0713+438, which was has a very stable flux density. Table 2 contains a log of the observations and a list of the measured flux densities. The VLA data for the first month after the burst were published by Taylor et al. (1998).

### 2.3 Owens Valley Radio Observatory (OVRO)

A continuum observation with OVRO was made at a central frequency of 99.98 GHz on 1998 Dec 13.42 and 14.40 UT, several months after the burst. The observations totalled approximately 14 hours on source in good 3 mm weather conditions. The total bandwidth was 2 GHz, resulting in an rms of $`0.7`$ mJy. The flux was calibrated using the extragalactic source 3C273, and the phase from the nearby quasar 0552+398.

A further observation was made at a central frequency of 99.5 GHz on 2001 November 26.48 UT, nearly 1400 days after the burst. The observation consisted of a single, 11-hr long track (9 hours on source) taken under excellent 3 mm conditions with four antennas. The total bandwidth was 4 GHz, resulting in an rms of $`0.5`$ mJy. We set the flux density scale using the extragalactic sources 3C84 as well as 3C273, and derived phase calibration from J0646+448.

### 2.4 James Clark Maxwell Telescope (JCMT)

We re-analyzed archival JCMT444The JCMT is operated by The Joint Astronomy Centre on behalf of the Particle Physics and Astronomy Research Council of the UK, the Netherlands Organization for Scientific Research, and the National Research Council of Canada. observations of this burst taken at frequencies of 650 GHz, 350 GHz and 220 GHz. Details of the observing procedure can be found in Smith et al. (1998). The data were reduced using the SCUBA User Reduction Facility (Jenness & Lightfoot, 1998) in the same manner as that for the recent GRB 010222 (Frail et al., 2002). Raw signals were flat-fielded to account for the small differences in bolometer response, extinction corrected, and de-spiked to remove anomalous signals above the 3-sigma level. Short time-scale sky variations were also removed using pixels around the edge of the array containing no source emission (Jenness, Lightfoot & Holland, 1998). A flux calibration factor was then applied to convert to Jy. Flux calibration factors (FCF) of 240 $`\pm `$ 15 Jy/V, 197 $`\pm `$ 13 Jy/V and 384 $`\pm `$ 82 Jy/V were applied to the 220 GHz, 350 GHz and 660 GHz data, respectively (Coulson, 2000). Table 3 contains a log of the observations and a list of the measured flux densities.

## 3 Afterglow Model

We interpret the data in the framework of the cosmological fireball model (e.g., Piran (1999)), in which an energetic explosion accelerates a small amount of matter to ultrarelativistic velocities. Internal shocks within this flow produce the burst event itself, while the relativistic shock propagating into the surrounding medium produces the afterglow.

In the model, we assume the ultrarelativistic shock transfers a constant fraction of its total energy to the magnetic field ($`ϵ_B`$) as well as a constant fraction to shocked, thermalized electrons ($`ϵ_e`$). As in strong subrelativistic shocks, the electrons are assumed to be accelerated into a power–law distribution of energies (P$`(\gamma _e)\gamma ^p`$) and they radiate via synchrotron emission under the influence of magnetic fields. This produces a broken power–law spectrum whose peak and spectral breaks evolve in time according to the shock’s behaviour (Sari, Piran & Narayan, 1998), set by the total energy in the shock, the geometry of the ejecta and the density of matter surrounding it.

The specific model we used to fit the data is outlined in some detail in Harrison et al. (2001). It includes the effects of inverse Compton (IC) scattering on the shock evolution and emitted radiation spectrum as prescribed by Sari & Esin (2001). It allows for a conical (jet-like) outflow geometry with half-opening angle $`\theta _{jet}`$ (using the treatment given by Sari, Piran & Halpern (1999)), and for expansion into a medium of constant density $`n`$, or a medium with a density gradient, i.e., $`nr^2`$. The latter density profile would be typical of a medium altered by the wind of a massive star (in this case the GRB progenitor). In addition, we calculate and include radiative corrections to the shock energy. This is a refined method as compared to the adiabatic evolution typically assumed in our previous work (Harrison et al., 2001; Berger et al., 2001a, 2000). We approximate the radiative evolution by treating the shock at each moment as though it were instantaneously adiabatic, with an energy calculated at that particular time. This is an appropriate treatment as long as losses are moderate, with the change in energy being slow.

In addition to the basic input physics describing the evolution of the fireball outlined above, the broadband model for GRB afterglows incorporates several additional effects resulting from propagation of the radiation between the fireball and the observer. These include interstellar scintillation in the centimeter radio regime, dust extinction in the optical/NIR regime, and a contribution to the emission from the host galaxy of the GRB. These features are evident in the data, and must be included to derive accurate model parameters.

Interstellar scintillation (ISS), due to the turbulent ionized gas of our Galaxy distorting wavefronts propagating to the observer, can be important at radio wavelengths (Goodman, 1997; Frail et al., 1997). We account for ISS by first estimating the fractional variation in the flux density expected by these distortions (Walker, 1998). This uncertainty in the model flux is added in quadrature to the statistical uncertainties in the measured values when estimating $`\chi ^2`$.

Dust within GRB host galaxies (either in the circumburst environment or along the line of sight) is a source of extinction for optical/NIR afterglows (Klose et al., 2000; Sokolov et al., 2001; Djorgovski et al., 2001) which must be accounted for in modeling the optical spectrum. The frequency dependence of the extinction curve is uncertain, since we know little about the ISM properties of the GRB 980329 host galaxy. For GRB 980329, the observed red optical to near-IR afterglow spectrum (Palazzi et al., 1998) suggests a steep extinction curve. The best characterized (relatively) steep extinction law is for the Small Magellanic Cloud (SMC). We therefore adopt the SMC bar extinction curve from Weingartner & Draine (2001) , which is parameterized by the extinction level $`A_V`$ in the rest frame of the host galaxy. Other extinction curves were used in model attempts and did not substantially change the results.

Depending on its luminosity, the host galaxy of the GRB may contribute a background level which dominates the total brightness at late times. This is observed as a flattening of the light curves, chiefly in the optical and IR. For GRB 980329, we include a constant term for the R, I, H and K to represent the host emission, and the values of these terms are fit by the model. The J-band only includes a single data point, so we interpolate the host term between the I and H fits. Evidence for host emission has also been observed at centimeter wavelengths (Berger, Kulkarni & Frail, 2001), again deduced from late-time flattening of the lightcurve. The radiation at these wavelengths would most likely result from synchrotron and thermal bremsstrahlung in galaxies undergoing substantial star formation (Condon, 1992). We include the possibility of centimeter host emission in our model, where we fit for the normalization at 1.43 GHz, and scale other bands as $`\nu ^{0.8}`$. In light of recent claimed detections of GRB hosts at submillimeter wavelengths (Hanlon et al., 2000; Berger, Kulkarni & Frail, 2001; Frail et al., 2002; Berger et al., 2001b, 2002), we consider the possibility of a submillimeter host, and allow for such a component in the fit, with a modified blackbody as in Frail et al. (2001), parameterized by its flux at 350 GHz. This last host component was found not to be required for our best model.

## 4 Best-Fit Broadband Model

The afterglow model described in §3 was fit to the data summarized in Tables 1-3. In addition, we included all previously published data in the X-ray (In ’t Zand et al., 1998), optical (Palazzi et al., 1998; Gorosabel et al., 1999; Reichart et al., 1999), and radio (Taylor et al., 1998) bands. We converted the X-ray measurements to flux values using a photon index of 2.4. We corrected the optical data for absorption in our Galaxy (Schlegel, Finkbeiner & Davis, 1998) before converting to flux densities using the factors in Bessell (1979) for the optical and Bessell & Brett (1988) for the near-IR bands. To account for any cross-calibration uncertainty, we have added in quadrature 5% uncertainties to all the measured fluxes.

The model which best describes this broadband dataset is a collimated outflow expanding into a constant density medium. The best model parameters, derived from least-squares maximizing the fit probability, are summarized in Table 4 and Figs. 1-6. The $`\chi ^2`$ for the fit is 116.4 for 92 degrees of freedom. Although we derived fits for three representative redshifts ($`z=1,2,3`$), the results of which are all shown in Table 4, we confine our detailed comments below to the $`z=2`$ solution. This choice of representative redshifts was made based upon the range of likely $`z`$ for this burst. Very high redshifts $`z>5`$ are not considered as they are not compatible with the underlying host’s colors, as detailed in §4.1. A redshift $`z<1`$ is considered implausible due to the lack of lines expected to be detected (if $`z<1`$) in several spectra taken of the host. The host is visible at optical wavelengths, and thus not completely obscured, so it is quite unlikely that the prominent star-formation-related oxygen line \[OII\]$`\lambda `$3727 or the H$`\alpha `$ line would not have been observable if the host is a faint galaxy at $`z<1`$. Many of the basic conclusions do not depend on the redshift, or can be easily scaled given the information below. In §5 we discuss some of the limitations of our best fit, as well as some alternate models which also fit the data, but only with unphysical parameters.

In the best model for $`z=2`$, the isotropic-equivalent fireball energy at the time when the fireball evolution becomes nearly adiabatic ($`E_{\mathrm{iso}}(t_{\nu _c=\nu _m})`$) is approximately $`10^{54}`$ ergs. The measured gamma-ray fluence for this GRB is $`F_\gamma =5.5\times 10^5`$ erg cm-2 (In ’t Zand et al. (1998)), so the isotropic gamma-ray energy is $`E_{\mathrm{iso}}(\gamma )=4\pi F_\gamma d_L^2(1+z)^16\times 10^{53}`$ erg. The large energy budgets inferred for both the shock $`E_{\mathrm{iso}}(t_{\nu _c=\nu _m})`$ and the emitted gamma-ray radiation $`E_{\mathrm{iso}}(\gamma )`$ derived assuming isotropy are greatly reduced in this model by the relatively large degree of collimation (jet opening angle of $`\theta _{jet}2^{}`$). A similar degree of collimation has been inferred previously in GRB afterglows: GRB 990510 has $`\theta _{jet}=3^{}`$ (Harrison et al., 1999) and GRB 000911 has $`\theta _{jet}=2^{}`$ (Price et al., 2002). For a two-sided jet this implies a total energy in the fireball shock of $`8.3\times 10^{50}`$ ergs, similar to the energy released in supernovae. Likewise, for $`z=2`$ the geometry-corrected gamma-ray energy is reduced to $`4.0\times 10^{50}`$ erg, a value that is in good agreement with the mean of $`5\times 10^{50}`$ erg derived from a larger sample (Frail et al., 2001). We note the total energy is similar ($`5.4\times 10^{50}`$) for $`z=1`$, and a factor of about three higher for $`z=3`$.

The energy quoted in Table 4 is a lower limit on the true initial energy of the blastwave since it is derived at a time $`t_{\nu _\mathrm{c}=\nu _\mathrm{m}}`$=6.1 d, when the lowest energy electrons can cool within the dynamical timescale of the system. Observationally, this corresponds to the time when the cooling break $`\nu _c`$ crosses the spectral break which results from the low-energy cutoff in the input electron energy spectrum, $`\nu _m`$. This criterion separates the two regimes of radiative losses; early times when radiant energy results in a decrease in the blastwave energy with time, and late times when the blast-wave evolution is adiabatic. For our best model, radiative losses are important, but not extreme even at early times, since the fraction of energy in radiating electrons is not dominant ($`ϵ_e`$=12%). We estimate that from the time the GRB ends ($``$ 10-100s post-trigger), to when the blastwave is nearly adiabatic, the energy drops by a factor of five. If we restrict the interval to begin when the first afterglow data were measured, ($`t=0.25`$ d for the first data point to $`t_{\nu _\mathrm{c}=\nu _\mathrm{m}}`$), the energy drops by a factor of 1.6. From the time $`t=t_{\nu _\mathrm{c}=\nu _\mathrm{m}}`$ until late times the energy drops by only 15% .

The energy in the fireball derived from our model significantly exceeds the emitted gamma-ray energy, i.e., $`E_{\mathrm{iso}}(t_{\nu _c=\nu _m})>E_{\mathrm{iso}}(\gamma )`$. As discussed above, the initial fireball energy will be even larger if radiative losses are taken into account. This is the case for the majority of GRB afterglows with energies derived from model fits (see e.g. Panaitescu & Kumar (2001a)). In the fireball model, it is likely that the energy remaining in the shock during the afterglow exceeds that of the prompt gamma-ray emission of the GRB event itself; this will be the case for a radiative efficiency of $`<`$ 50% during the GRB. The radiative efficiency of internal shocks driving the prompt GRB are expected to be $``$ 10%, not $``$ 50%, leaving most of the initial shock energy in the fireball (Guetta, Spada & Waxman, 2001).

The ratio of the energy fraction in magnetic field, $`ϵ_B=`$17%, to that in the electrons, $`ϵ_e`$=12%, determines the relative importance of Comptonization. Compton scattering can contribute significantly to the total cooling rate if the ratio of electron to magnetic energy fractions is greater than unity (Sari & Esin, 2001). In our best model for GRB 980329, this ratio is of order unity, so Comptonization does not dominate the electron cooling, but it is a non-negligible effect. Flux from Compton scattering can in fact be seen peaking in the X-rays in Figs. 4 and 6, largely as a result of the steep electron spectral index. For an electron energy spectral index of 2, with equal energies in each logarithmic frequency interval (the infinite-energy limit), the IC luminosity would be lower than the synchrotron luminosity at all frequencies, and no IC peak would be observable in the spectrum. The index of $`p=2.9`$ derived for the best model puts less energy in each successive decade above the peak frequency. With a significant circumburst electron density providing a non-negligeable opacity to Compton scattering, the peak of the IC flux density, above the synchrotron peak, dominates the total flux near the X-rays.

The circumburst medium density derived from the model, $`n=20`$ cm-3, is comparable to that of GRB 000926 (Harrison et al., 2001) and that of several other bursts (Panaitescu & Kumar (2001a)). This relatively high density, in reference to an average galactic ISM, can be inferred from the measured value of the self-absorption break $`\nu _a`$ (Granot, Piran & Sari, 1999), which is shown in Fig.6. The frequency $`\nu _a`$ depends upon other fundamental parameters besides the ambient density ($`\nu _an^{3/5}ϵ_e^1ϵ_B^{1/5}E_{iso}^{1/5}`$), and in this particular case the high $`\nu _a`$ results from the relatively high density, combined with a moderate electron energy fraction. Models with low circumburst densities (i.e., $`n10`$ cm-3) cannot be fit to the data by adjusting the energy and electron fraction. Highly radiative models, with large electron energy fractions, can provide reasonable fits to the data, but the densities (depending on the redshift of the fit) vary from approximately the same as to an order of magnitude greater than those for the best model (see §5).

### 4.1 Host Galaxy

We infer the presence of measurable flux from the GRB host galaxy in the optical from the late-time flattening of the lightcurves (Fig. 3). In the final epoch of our optical observations the measured brightness of R=26.53$`\pm `$0.22 (see Table 1) is essentially entirely due to the host, and is only 1.5 mag fainter than the median R magnitude for known GRB host galaxies (Djorgovski et al., 2001). Our measured magnitude differs from the preliminary value of R=27.7$`\pm `$0.3 given by Holland et al. (2000). This discrepancy seems to be due to a mis-identification of the host in their HST images which they identify with a source 0.5′′ southwest of the GRB position. Further work by that group (Jaunsen et al., 2002) has identified the host at a position consistent with that of Bloom et al. (2002), which established with improved astrometry the host whose R-band magnitude is given here. Bloom et al.’s (2002) measured offset for the GRB from the host center is only 37$`\pm `$48 mas, corresponding to a host-normalized offset of 0.215$`\pm `$0.291, placing GRB 980329 within the half light radius of its compact host. Our I band host flux values from Table 1 are in good agreement with Jaunsen et al.’s (2002) late-time I band host measurement.

In Fig. 5 we plot the spectral energy distribution of the afterglow at 0.7 days after the burst along with the late-time measurements from Table 1, assumed to be due to the host galaxy. We corrected all points for extinction in our Galaxy. Palazzi et al. (1998) first noted the steep spectral slope between the R and I bands seen in the early time afterglow, and Fruchter (1999) suggested that this “dropout” of the R band could be produced by absorption from the Ly$`\alpha `$ forest if the redshift $`z`$ of this burst was greater than five. This steep slope, however, is not reflected in the host spectrum. The afterglow at 0.7 days after the burst is significantly redder (R$``$I=2.7$`\pm `$0.4) than the host galaxy (R$``$I=0.2$`\pm `$0.3) itself (note that the quoted R$``$I above are corrected for Galactic extinction). From the absence of a strong $`RI`$ break in the host spectrum, we can rule out a redshift of $`z>5`$ for GRB 980329.

Another result that emerges from the modeling is the presence of significant dust extinction in the host galaxy. Early attempts to model the optical data for this burst (Palazzi et al., 1998; Lamb, Castander & Reichart, 1999) also found that the spectrum was substantially reddened by dust. Our fitted host $`A_\mathrm{V}`$ corresponds to hydrogen column density of $`N_H2\times 10^{21}`$ cm-2, assuming a gas-to-dust ratio similar to that of the Milky Way (Predehl & Schmitt, 1995; Reichart, 2001). In’t Zand et al. (1998) used the X-ray spectrum of the afterglow to derive a column density $`N_H=1.0\pm 0.4\times 10^{22}`$ cm-2, with a 99% confidence range of $`1.3\times 10^{21}`$-$`1.5\times 10^{22}`$ cm-2, after subtracting a Galactic contribution of approximately 0.8$`\times `$1021cm-2 (Dickey & Lockman, 1990). To translate this measured column into the restframe of the host galaxy requires multiplying by a factor $`(1+z)^{8/3}`$. For a redshift of 2, this exceeds the value derived from optical extinction. Such discrepancies have been noted before (Vreeswijk et al., 1999; Galama & Wijers, 2001), and are taken as evidence of significant dust destruction in the circumburst medium out to a radius of order 10-20 pc (Waxman & Draine, 2000; Fruchter, Krolik & Rhoads, 2001; Reichart, 2001). In the case of GRB 980329, however, the redshift is not known, and the uncertainties in the dust extinction law make it difficult to claim evidence for dust destruction.

There is a suggestion of a 1.43 GHz radio host in the data, with the model requiring a flat, positive component on average. Although the addition of this component improves the fit, the significance of a nonzero radio host parameter is only 3$`\sigma `$. If real, a radio host at $`25\mu `$Jy would be about 1/3 the host flux density found for GRB 980703 (Berger, Kulkarni & Frail, 2001). We note that the submillimeter data are in good agreement with the afterglow model (see Fig.2), and we do not require any host contribution in this band.

## 5 Alternate Models and Limitations

The best fit model given in §4, while not a unique interpretation of the data, is a self-consistent solution that derives reasonable values for the blastwave energy $`E_{\mathrm{iso}}(t_{\nu _c=\nu _m})`$, the opening angle $`\theta _{jet}`$, the ambient density $`n`$ and the parameters of the shock ($`p`$, $`ϵ_e`$, $`ϵ_B`$). A collimated outflow expanding into a constant density medium describes all the data well, addressing puzzling features described in previous attempts to fit this afterglow (Fruchter, 1999; Lamb, Castander & Reichart, 1999; Smith et al., 1999) without invoking a very high redshift or other additional components. The unusual features of this data set include the submillimeter excess, the very red $`RI`$ afterglow colors and their relation to the host, and the observed decline of the peak flux density $`F_m`$ with time (or equivalently with decreasing frequency $`\nu _m`$, since $`\nu _m`$ decreases with time).

This last feature warrants further explanation, since it is an effect that can only result from a finite number of physical causes. The “peak flux cascade” can be readily seen in Figs. 1 and 2, where the peak flux density is 2.5 mJy at 350 GHz but declines to 1.5 mJy at 90 GHz and further falls to 0.35, 0.2 and $`<`$0.1 mJy at 8.46 GHz, 4.86 GHz, and 1.43 GHz, respectively. A fit to the data near maximum between 4.86 GHz and 350 GHz gives a power-law slope of 0.59$`\pm `$0.07. A similar behavior was also observed for the afterglow of GRB 970508 (Galama et al., 1998; Frail, Waxman & Kulkarni, 2000). Within the context of the standard fireball model there are three ways to produce this behavior. First, if the flux evolution is observed after collimation of the ejecta becomes evident in the light curve decay (i.e., post-jet), then $`F_m\nu _m^{1/2}`$ (Sari, Piran & Halpern (1999)). For our best model, this is what produces the observed peak flux cascade. Alternately, a density gradient $`nr^2`$ in the surrounding medium (a stellar wind model) will give $`F_m\nu _m^{1/3}`$ (Chevalier & Li, 1999). Finally, radiative losses can produce a peak flux cascade. However, unless these are severe (namely most of the shock energy in electrons), the effect is quite weak. For example, for $`ϵ_e0.1`$, $`F_m\nu _m^{0.08}`$, while for $`ϵ_e1`$, $`F_m\nu _m^{0.37}`$ (Cohen, Piran & Sari, 1998).

The other models we derived that fit the primary characteristics of the data all required unusual physical assumptions. For example, we found a solution with extreme radiative corrections (100% of the shock energy going into electrons) that could reproduce the observed peak flux cascade. Formally, this model fits the data better than the “best model” presented in the previous section, however, it reaches the unphysical edge of parameter space, and with extreme radiative losses our near-adiabatic treatment breaks down and cannot be fully trusted. The highly radiative model is isotropic, and for $`z=2`$, has the following parameters: $`E_{\mathrm{iso}}(t_{\nu _c=\nu _m})2\times 10^{52}`$ ergs, $`n20`$ cm-3, $`p2.02`$, $`ϵ_B0.17`$, $`ϵ_e1`$, A(V) $``$ 1.1 and a centimeter host of $``$ 17 $`\mu `$Jy at 1.43 GHz. In addition to the unphysical assumption about the electron energy partition, this model also has an electron spectral index approaching two, and hence a diverging total energy. This radiative model only accounts for $``$ 1/2 of the 350 GHz flux, suggesting an underlying submillimeter host of $``$ 0.7 mJy (this component improves the fit at $``$ 3 $`\sigma `$ level). This submillimeter host flux level is just below the sensitivity limit of current instruments and would likely not be detectable at late times, if it indeed exists. We consider the collimated solution presented in §4 to be the best model as it is the best fit of the models with realistic parameters. This best model reproduces the flux cascade with a relatively narrow collimation angle (early jet-break). We note that the $`p=2.88`$ we derive for the best $`z=2`$ model is somewhat larger than found for other afterglows, which generally fall in the range $`p=2.22.4`$. It is, however, physically reasonable, and we regard all the parameters associated with the best model as acceptable. The highly radiative solution is plotted along with the best model in the lightcurves presented in Figs. 14. The fit is visibly somewhat better, but at the cost of unphysical assumptions concerning the underlying parameters of the fireball.

The most serious limitation to modeling the afterglow of GRB 980329 is the lack of a good redshift estimate. Even a fairly comprehensive dataset such as this cannot constrain fundamental parameters without knowing the distance. This is because the synchrotron emission can be reproduced at different $`z`$ simply by rescaling the physical parameters by appropriate powers of $`(1+z)`$. Only “second order” effects such as host extinction, the IC component and radiative corrections do not directly re-scale with $`(1+z)`$. In principle, we could include $`z`$ as a free parameter in the model, and fit for the best value. In practice, however, the combination of the sparseness of the real dataset, and the uncertainties in the model prevent any unique redshift determination. This is evident from Table 4, where we show a good fit with reasonable physical parameters for all three redshifts.

The absence of a strong break in the host galaxy spectral colors (Fig. 5) allows us to place an upper limit on the host redshift. Spectral energy distribution (SED) fittings to host colors may place a stronger constraint (e.g. Castander & Lamb (1999), concerning the 970228 event), especially if the factor of $``$ four decline from K-band to I-band can be modelled by a Balmer break. Further HST imaging in optical bands might prove fruitful in this regard.

The absence of a redshift is responsible in part for the relatively large and uncertain estimate of the electron index $`p`$ in our best model. The post-jet evolution of the optical light curves is determined by $`p`$ (Sari, Piran & Halpern, 1999). Unfortunately, this data is sparsely sampled at early times (Fig. 3), prior to when the host galaxy dominates the light. Likewise, the index $`p`$ determines the shape of the synchrotron spectrum, and should therefore be derivable from the measurements. At optical wavelengths, however, there is a degeneracy between $`p`$ and the dust extinction law — the latter of which can depend sensitively on $`z`$. Lacking any knowledge about the extinction properties of dust in high redshift galaxies we adopted the SMC extinction law. Neither the X-ray spectral slope, nor the X-ray to optical flux ratio can break this degeneracy, since the contribution from inverse Compton scattering alters the X-ray flux normalization as well as the spectrum.

Finally, we note that both our best solution and the highly radiative solution fail to predict the early radio emission at 8.46 GHz ($`t<3`$ d). This level of fluctuation is too great to be accounted for by the estimated ISS effects. Prompt, short-lived radio emission in excess of the normal afterglow component has been detected toward other GRBs (Kulkarni et al., 1999; Frail et al., 2000; Harrison et al., 2001). This is usually attributed to radiation from a reverse shock (Sari & Piran, 1999), and we suggest that this may explain the bright early-time radio point for GRB 980329.

## 6 Conclusions

The new observations presented here, along with the data in the literature, form a complete collection of broadband measurements for the afterglow from the bright burst GRB 980329 beginning from early to late times. We have used these data together with a comprehensive fireball model to deduce fundamental physical parameters of the event (energy, geometry, density) and to measure the properties of the host galaxy.

The late-time optical/NIR data shows that the host is significantly bluer than the afterglow. Thus we reject the hypothesis that the very red afterglow colors were due to Lyman-$`\alpha `$ absorption in the intergalactic medium to a very high redshift source.

All of the afterglow’s features can be explained over a wide range of $`z`$ by a model in which the ejecta are collimated in a jet. Significant dust extinction is inferred within the host galaxy and a moderately high circumburst density $`n20`$ cm-3 is required. Although this collimated model is not a unique solution to the data, it explains the red optical/NIR color and the cascade in the peak flux from submillimeter to centimeter wavelengths (see §5) without resorting to an extreme redshift $`z>5`$ or requiring additional complications in the host galaxy’s properties. Models which invoke isotropic geometry require such complications, circumburst densities up to 10 times higher, as well as large, unphysical radiative corrections. Correcting the isotropic-equivalent gamma-ray energy release for the collimation, for $`z=2`$ we obtain $`E(\gamma )=E_{\mathrm{iso}}(\gamma )\theta ^2/2=4.0\times 10^{50}`$ ergs. This value is typical of other events to date and is easily accounted for by current progenitor models.

RS acknowledges support from the Fairchild Foundation and from a NASA ATP grant. JSB acknowledges a grant from the Hertz foundation. Research with the Owens Valley Radio Telescope, operated by Caltech, is supported by NSF grant AST96-13717