text:
# Simulated Radio Images and Light Curves of Young Supernovae

## 1 Introduction

Although supernovae (SNe) have been observed over many centuries, it is only in the past few decades that it has been possible to detect them at radio wavelengths soon after outburst. Since the first radio SN (RSN) detection by Gottesman et al. (1972), many more RSNe have been observed, and a few have been discovered by virtue of their radio emission alone (e.g. SN1981K and SN1986J). One of the earliest RSNe to be discovered, SN1979C, has been extensively monitored at radio wavelengths for almost two decades. Many others, such as SN1986J, SN1980K, SN1993J and the well-known SN1987A, have also been the subject of extended studies.

Weiler et al. (1986) and Weiler, Sramek & Pangia (1996) identified several common characteristics of most RSNe: a) the emission is non-thermal with very high brightness temperatures; b) they appear first at high frequencies and then progressively at lower frequencies; c) at each frequency there is an initial rapid increase and then, after reaching a maximum, the radio flux density declines as a power-law in time; d) the spectral index $`\alpha `$ (where the flux density $`S\nu ^\alpha `$) between any two frequencies is initially negative, increases rapidly, and ultimately approaches a constant positive value. It is generally accepted that the radio emission from SNe is synchrotron radiation from ultra-relativistic electrons accelerated at the expanding supernova shock (or shocks). The declining power-law frequency spectrum seen after the SN radio emission peaks is consistent with optically-thin synchrotron emission from a power-law energy distribution of ultrarelativistic electrons. The observed characteristics (a)–(d) are thus often loosely summarized by the statement that radio SNe exhibit an optically-thick rise followed by an optically-thin power-law decline. Recent observations of deviations of RSNe light curves from long-term trends have been attributed to changes in the ambient density (Montes et al., 1998, 2000; Panagia et al., 1999).

The radio emission from supernovae is thought to arise in the hydrodynamically unstable interaction region between the SN ejecta and a circumstellar medium (CSM) formed by mass loss from the progenitor star. This region contains both shocked supernova ejecta and shocked swept-up material, separated by a contact discontinuity, and is bounded on the outside by a forward shock and on the inside by a reverse shock that travels back into the ejecta with respect to the contact discontinuity. The recent observation of a hard X-ray tail in the spectrum of SN1006 has shown that electrons can be accelerated up to TeV energies at a supernova shock (Koyama et al., 1995). The electrons may be subject to a magnetic field which is essentially that of the progenitor star wind, or which is amplified as a result of turbulence behind the supernova shock. However, there are as yet no detailed models that satisfactorily treat the hydrodynamics of the supernova expansion, the electron acceleration, and the origin and evolution of the magnetic field, that together conspire to give the observed radio properties.

An analytic treatment of the evolution of an idealized RSN by Chevalier (1982a, 1984) is often referred to as the “mini-shell model”. This model assumes that the expanding supernova is self-similar, with the shock radius expanding as a power-law in time, $`r_st^m`$, where $`m`$ is a constant (Chevalier, 1982a). The energy distribution of the relativistic electrons is assumed to be a power law of constant spectral index, and the energy density of the relativistic electrons and that of the magnetic field are assumed to each be a constant fraction of the thermal energy density behind the expanding supernova shock. Free-free absorption in the ionized CSM (Chevalier, 1982b), which decreases as the SN shock expands through the overlying CSM, can produce an initial optically-thick increase in flux density. Decreasing synchrotron self-absorption (SSA) in the expanding radio source may also be important (Chevalier & Dwarkadas, 1995; Chevalier, 1998).

The mini-shell model, and a range of empirical variants of it (cf. Van Dyk et al. 1994a; Weiler et al. 1996), have been used to parameterize the radio emission from a number of young supernovae. The mini-shell model works well within its limitations (Chevalier, 1986) and can successfully reproduce the evolution of some RSNe. However, some of the empirical extensions of the mini-shell model which have been used to fit recent RSNe no longer reflect the physical basis of the original model. Furthermore, the mini-shell model assumes self-similarity and is therefore not applicable to RSNe whose ejecta are not represented by a power-law density profile, or whose expansion is not self-similar for other reasons.

Attempts to compute RSN models that incorporate a more detailed treatment of the physical processes involved have generally focused on trying to describe in more detail the acceleration of electrons to relativistic energies. Duffy, Ball, & Kirk (1995) considered a time-dependent treatment of diffusive shock acceleration in their study of SN1987A. They used a 2-fluid model in which the accelerated ions provided an additional component to the total pressure. Fransson & Björnsson (1998) adopted a different approach in their study of SN1993J, opting to fit the individual radio spectra at the various observed dates, rather than concentrating on the light curves at different wavelengths. Their calculated light curves are consistent with observations up to an age of about 4 years. Jun & Jones (1999) have included relativistic electron transport in their magnetohydrodynamic simulations of a young supernova remnant interacting with an interstellar cloud.

We have developed a robust and general method to relate the hydrodynamics of a supernova explosion to the evolving radio emission from these sources. We first compute the evolution of the SNe using spherically-symmetric, high-resolution numerical hydrodynamic simulations. Our method allows for the use of arbitrary profiles for the SN ejecta and circumstellar medium, a significant generalization from analytic models which are restricted to the self-similar expansion which results from power-law ejecta and CSM profiles. The light curves are produced by a detailed radiative transfer calculation using ray tracing methods. Synchrotron self-absorption is implicitly included in the calculation, free-free absorption in the surrounding medium is explicitly calculated, and in one case we include internal free-free absorption. The relativistic particle energy density and the magnetic energy density can be independently specified, and in general are a function of radius. In the simplest implementation, used for the results presented in this paper, these are taken to be a fixed fraction of the local thermal energy density, similar to the assumptions used in the mini-shell model. A further generalization of the method, computing the local relativistic particle distribution using advanced particle acceleration codes, is already underway.

The results of calculations for a model for SN1980K with power-law ejecta and CSM density profiles – giving a self-similar expansion – are presented as a test of our method. We then model the emission from a model of SN1993J with a complicated ejecta profile which is motivated by detailed modeling of the optical and x-ray emission, and which results in an expansion that is far from self similar.

In this paper we first describe our method, give the results of calculations of the radio emission from model RSNe, and compare the model emission with the observational data for SN1980K and SN1993J. The radiation transfer calculation and the manner of computation are outlined in §2. In §3 we describe the process of calculating the model radio emission, including the hydrodynamic simulations that are used as an input to the radiative transfer code. As a test of the method and the code we calculate the radio evolution from a hydrodynamical simulation of a self-similar shock expansion, and compare the results with the behavior expected from an analytic treatment of this special case. In §4 we discuss the process of fitting observed radio data, and present the results of model calculations for SN1980K, which can be fitted using the minishell model, and SN1993J which cannot. Finally, in §5 we evaluate the results, discuss the implications of this work, and outline future prospects.

## 2 Radiation Transfer

The basic radiation transfer equation (e.g. Rybicki & Lightman 1979) is

$$\frac{\mathrm{d}I_\nu (x)}{\mathrm{d}x}=j_\nu (x)\kappa _\nu (x)I_\nu (x),$$
(1)
where $`I_\nu (x)`$ is the flux density, $`j_\nu (x)`$ is the emissivity and $`\kappa _\nu (x)`$ is the opacity, all at position $`x`$. The optical depth between position $`x_0`$ and $`x`$ is defined as

$$\tau _\nu (x_0,x)=_{x_0}^x\kappa _\nu (x^{})dx^{}.$$
(2)
The solution to Equation (1) between $`x_0`$ and $`x`$ can then be written as

$$I_\nu (x)=I_0(x_0)e^{\tau _\nu (x_0,x)}+_{x_0}^xe^{\tau _\nu (x^{},x)}j_\nu (x^{})dx^{}.$$
(3)
We assume that there is no background radiation (i.e., $`I_0=0`$), and denote the diameter of the SN by $`r_s`$ so that the equation we use for our radiation transfer calculations is:

$$I_\nu (x>r_s\mathrm{sin}\theta )=_0^{r_s\mathrm{sin}\theta }e^{\tau _\nu (x^{},r_s\mathrm{sin}\theta )}j_\nu (x^{})dx^{}$$
(4)
where $`\theta `$ is the angle between the line of sight and the radial line to the intersection of the line of sight and the shock.

The hydrodynamical simulations discussed here are 1-dimensional. The code converts these data into a 3-dimensional sphere embedded in a Cartesian grid. The hydrodynamical quantities at each grid point are used to calculate $`\kappa _\nu `$ and $`j_\nu `$ (as described in the next two sections). These are then integrated along each line of sight to calculate $`\tau _\nu `$ and finally the flux density, $`I_\nu (r)`$ for that line of sight. The integration is performed using the trapezoidal rule. The numerical grid is sufficiently fine that the method is adequate, and the trapezoidal rule is very robust. A 2-dimensional array of surface brightness (as a function of position across the source) is built, and this array can be used to produce a synthetic image or summed to obtain a point on a light curve.

### 2.1 Synchrotron Emission & Absorption

The synchrotron emissivity ($`j_\nu `$) and opacity ($`\kappa _{\nu ,SSA}`$) depend on the distribution of the radiating electrons, the magnetic field ($`B`$), and the frequency ($`\nu `$). We assume that the energy distribution of ultra-relativistic electrons is a power law $`n(E)=n_0E^\gamma `$ where $`E`$ is the particle energy, $`\gamma `$ and $`n_0`$ are constants, and the total number density is $`n(E)dE`$. The hydrodynamical simulations do not provide direct information on the magnetic field or the high energy particle distribution. Here we follow Chevalier (1982b) and assume that the magnetic energy density, $`u_B`$, and the relativistic particle energy density, $`u_{rel}`$, are proportional to the thermal particle internal energy density, $`u_e`$, such that

$`u_B=`$ $`B^2/(2\mu _0)=\zeta _Bu_e,`$ (5)
$`u_{rel}=`$ $`{\displaystyle n(E)EdE}=\zeta _{rel}u_e,`$ (6)
where $`\zeta _B`$ and $`\zeta _{rel}`$ are constants. In the Chevalier (1982b) mini-shell model these constants are set to around 1%. The common assumption of energy equipartition corresponds to $`\zeta _B=\zeta _{rel}`$. There is some justification for such a relation between the energy densities: $`u_e`$ is high between the forward and reverse shock, which is where the magnetic field is likely enhanced and the particles accelerated to relativistic speeds. However such arguments are weak and ultimately need to be replaced by physical processes, such as particle acceleration.

The thermal particle energy density is related to the thermal pressure ($`p`$), which is output by the hydrodynamic simulations, through the adiabatic index ($`\mathrm{\Gamma }`$):

$$p(r,t)=(\mathrm{\Gamma }1)u_e(r,t)$$
(7)
Following Mioduszewski, Hughes, & Duncan (1997) equations (5), (6) and (7) can be used to write the synchrotron emissivity and opacity as

$`j_\nu (r,t)(\gamma 2)\mathrm{\Delta }\zeta _B^{\frac{1}{4}}\left({\displaystyle \frac{p(r,t)}{\mathrm{\Gamma }1}}\right)^{\frac{\gamma +5}{4}}\nu ^{\frac{\gamma 1}{2}},`$ (8)
$`\kappa _{\nu ,SSA}(r,t)(\gamma 2)\mathrm{\Delta }\left({\displaystyle \frac{p(r,t)}{\mathrm{\Gamma }1}}\right)^{\frac{\gamma +6}{4}}\nu ^{\frac{\gamma +4}{2}},`$ (9)
where $`\mathrm{\Gamma }`$, $`\gamma `$, $`\zeta _B`$, $`\zeta _{rel}`$ and $`\nu `$ are input parameters, $`\mathrm{\Delta }=\zeta _B^{\frac{\gamma +2}{4}}\zeta _{rel}`$, and $`p(r,t)`$ is obtained from the hydrodynamic simulations.

### 2.2 Free-Free Absorption

In contrast to the synchrotron emission, free-free absorption (FFA) depends upon the thermal density and temperature of the plasma and therefore can be calculated with a minimum of assumptions. The opacity ($`\kappa _{\nu ,FFA}`$) goes as

$$\kappa _{\nu ,FFA}\frac{n_en_i}{n_\nu \nu ^2T^{3/2}}\mathrm{ln}C(T,\nu )$$
(10)
where $`n_\nu `$ is the refractive index and $`C(T,\nu )`$ is a constant (Lang, 1980). For most of the models for RSNe considered here the free-free absorption arises from the overlying circumstellar medium (CSM) and $`n_\nu 1`$. In the hydrodynamical simulations the CSM density $`\rho _{\mathrm{CSM}}`$ is assumed to fall off as a power law, $`\rho _{\mathrm{CSM}}r^s`$ where $`r`$ is the radial distance from the progenitor and $`s`$ is a constant. We also assume that a constant fraction of $`\rho _{\mathrm{CSM}}`$ is ionized. The temperature ($`T`$) of the CSM is assumed to be constant with radius and is an input parameter. The free-free optical depth, $`\tau _{\nu ,FFA}`$ is obtained by integrating from the outer edge of the SN, determined by the radius of the forward shock ($`r_s`$), to infinity, whence:

$$\tau _{\nu ,FFA}=\frac{r_s^{12s}}{2s1}K_{FFA}$$
(11)
where $`K_{FFA}=\kappa _{\nu ,FFA}r_s^{2s}`$ is independent of $`r_s`$.

Although the external temperature is a parameter in the hydrodynamical simulations, it does not have a strong effect on the pressure between the forward and reverse shocks (because the pressure of the CSM is much less than the pressure in the post-shocked region), where most of the radiation is produced. It is therefore reasonable to use a single hydrodynamical simulation to calculate the resulting synchrotron emission for a range of assumed CSM temperatures.

For internal free-free absorption ccombined with synchrotron self-absorption (only considered in §4.2.2), $`\kappa _{\nu ,FFA}`$ is added to $`\kappa _{\nu ,SSA}`$ and the total $`\kappa _\nu `$ is used in the radiation transfer calculation outlined at the beginning of §2.

## 3 Computing Light Curves

We calculate light curves in two steps. First, a hydrodynamical simulation is run for a given SN ejecta profile and a given CSM profile. The results of a simulation comprise the spatial profile of hydrodynamical quantities at different epochs in the evolution of the SN. The particular outputs we require are the pressure profiles at each epoch, which are used as inputs to the radiation transfer code to calculate the radio flux densities of the SN at each epoch.

### 3.1 Hydrodynamic Simulations

We model the hydrodynamics of the interaction of the expanding supernova ejecta with the circumstellar medium using high-resolution simulations. All the simulations described here were carried out with the VH-1 code, a 3D finite-difference code (here used as 1D) based on the Piecewise Parabolic Method of Colella & Woodward (1984). The calculations are carried out in Lagrangian co-ordinates and re-mapped onto an Eulerian grid.

The VH-1 code uses an expanding grid which is essential for simulations where the dimensions change by orders of magnitude over a single run. The grid tracks the outer shock and the outer grid radius increases accordingly. The rate of expansion of the grid is adjusted such that most of the grid is taken up by the interaction region, bounded by the inner (or reverse) shock, and the outer (or forward) shock. Resolving this high-pressure region in detail is the highest priority since it is the source of most of the radio emission. Over a large part of the evolution almost a half to two-thirds of the grid is occupied by the interaction region, ensuring that it is sufficiently well sampled.

Setting up the initial conditions for a simulation requires a model for the SN ejecta and a model for the CSM that the ejecta are to expand into. For Type II supernovae the most commonly used ejecta models have a density profile which is uniform inside some reference radius and decreases as a power law with radius in the outer parts. The interaction of power-law ejecta with a power-law circumstellar medium can be described analytically by the self-similar solutions of Chevalier (1982a) and Nadyozhin (1985).

If the CSM is formed via a constant wind from the progenitor star before it exploded as a SN, then the wind density of this medium would vary as $`r^2`$ where $`r`$ is the radial distance from the star. Ultimately, as the wind merges with the interstellar medium, the density profile must approach a uniform value. Many factors, such as evolution of the wind or structure in the pre-existing ambient medium, can lead to different density distributions at smaller radii.

### 3.2 Testing the code: Self-Similar Solutions

We first test the code by computing the emission from a self-similarly evolving SN which can be compared with analytic results. We consider SN ejecta with a power-law density profile $`\rho v^n`$ where $`\rho `$ is the density, $`v=r/t`$ is the velocity of the freely-expanding ejecta at radius $`r`$, and $`n>5`$ is a constant. When such ejecta expand into a circumstellar medium with density $`\rho _{CSM}r^s`$, the radius of the contact discontinuity ($`r_1`$) evolves as $`r_1t^{(n3)/(ns)}=t^m`$ (Chevalier, 1982a). The intrinsic (unabsorbed) flux density, given the assumed scalings of the magnetic field and relativistic electrons, then scales with frequency and time as $`S\nu ^\alpha t^\beta `$ where $`\beta =\alpha 3(1m)`$ (Chevalier, 1982b) and $`\alpha =(\gamma 1)/2`$.

Most observed SN light curves are characterized by a rapid, optically-thick rise and a more gradual optically-thin decline (for examples see Weiler, Sramek & Pangia (1996)). The optically-thick rise in these models is caused by decreasing absorption as the SN expands. This absorption can be due to free-free absorption and/or synchrotron self-absorption. Since the free-free absorption is directly related to the density of the overlying CSM, the decreasing CSM density with radius causes the free-free absorption to decline. Less directly, the synchrotron self-absorption decreases because the declining density in the CSM causes the the pressure between the forward and reverse shocks to decrease and thus the magnetic field and density of relativistic particles decrease. We consider three limits in which the time dependence of the radio flux density can be expressed analytically (Chevalier, 1982a, 1984):

1. optically thick due to internal synchrotron self-absorption

$$Ft^{(2\frac{2s}{4})m+\frac{1}{2}}$$
(12)
2. optically thin (i.e. with negligible optical depth due to either synchrotron self-absorption or free-free absorption)

$$Ft^{3m+[(2s)m2]\frac{\alpha +3}{2}}$$
(13)
3. optically thick due to external free-free absorption.

$$\mathrm{ln}\left(\frac{F_{\mathrm{thin}}}{F}\right)t^{(12s)m}$$
(14)
where $`F_{\mathrm{thin}}`$ is the flux density in the absence of absorption.

For the specific self-similar case we consider, $`n=12`$ and $`s=2`$ and we set $`\gamma =2.5`$. Thus $`m=(n3)/(ns)=0.9`$. Table 1 shows the analytic values of the time indices for these three cases together with the values obtained from fitting the light curves produced by the radiation transfer code for the corresponding hydrodynamical simulation. The agreement is remarkably good.

## 4 Fitting Light Curves

Given the results of a numerical hydrodynamical simulation, 4 further parameters must be specified to calculate the evolving radio emission: the spectral index ($`\alpha `$ or $`\gamma `$), the temperature of the CSM ($`T_{CSM}`$), and the energy density constants, $`\zeta _B`$ and $`\zeta _{rel}`$. The different effects that these 4 parameters have on the light curves, and their implications for fits to observed data, can be summarized as follows:

influences only the external free-free absorption. Most importantly, this parameter determines the time when the optical depth is small enough that the source can first be observed (i.e. when the SN turns on in the radio) if free-free absorption is the dominant absorption mechanism at that time. Increasing the ionization fraction of the CSM can mimic decreasing the temperature (see equation 10). A model in which the CSM temperature decreases with radius could also produce a low optical depth early in the SN evolution. The fitted model parameter $`T_{CSM}`$ therefore does not necessarily reflect the true temperature of the CSM.

determines the frequency dependence and the slope of the optically-thin decline, as can be seen from inspection of equation (13). Its other effects can be balanced by the constants, $`\zeta _B`$ and $`\zeta _{rel}`$.

affect the relative levels of the optically-thin and thick parts of the light curve, but do not change the slope of the rise or fall of the light curve. When the late optically-thick part of the light curve is dominated by synchrotron self-absorption these parameters determine the position of the peak of the light curve.

We calculate the flux density in arbitrary units, and finally scale the light curves to the observed data. The values of $`\zeta _B`$ and $`\zeta _{rel}`$ must satisfy the physical restriction $`\zeta _{rel}+\zeta _B1`$; but cannot be determined independently. When synchrotron self-absorption is unimportant $`\zeta _B`$ and $`\zeta _{rel}`$ can be absorbed into the overall scaling parameter and so are not determined by the model fit, and are not quoted. When synchrotron self-absorption is important, any two models with the the same value of $`\mathrm{\Delta }=\zeta _{rel}\zeta _B^{\frac{\gamma +2}{4}}`$ produce the same radio emission apart from an arbitrary scale factor.

We follow an iterative procedure to fit the observed evolution of a radio supernova. Having chosen and run an appropriate hydrodynamical simulation, we choose a set of the four additional model parameters and calculate the radiative transfer to obtain the radio evolution. The four parameters are adjusted on the basis of a comparison of the computed radio emission with the observed data, and the radiative transfer calculation is repeated for the same hydrodynamical simulation results.

### 4.1 SN1980K

SN1980K is a Type II supernova that was searched for at radio wavelengths soon after optical maximum, but not detected until some 35 days later (Weiler et al., 1986). Weiler et al. (1992) fit a Chevalier mini-shell model, with free-free absorption in an $`r^2`$ CSM, to the data and obtain a reduced $`\chi ^2`$ of 2.7. A fit to a more extensive data set by Montes et al. (1998) provides a similar value of reduced $`\chi ^2`$, but the model used is not directly comparable to our calculations.

Our method is a generalization of the Chevalier mini-shell model, and it provides a similarly satisfactory fit to the data for a simple power-law ejecta density profile expanding into a CSM with $`\rho _{CSM}r^2`$. Figure 1 shows the observed radio evolution of SN1980K at two frequencies together with two such model fits. The left panel shows a model with only free-free absorption in the CSM. The implied $`T_{CSM}`$ is quite high ($`3\times 10^5\mathrm{K}`$), but is consistent with the value $`T_{CSM}10^510^6\mathrm{K}`$ found by Lundqvist & Fransson (1988) at early times. In any event, $`T_{CSM}`$ may not indicate the true temperature of the CSM, as explained at the beginning of this section. A fit using the same hydrodynamical model when synchrotron self absorption is important, but there is negligible free-free absorption in the CSM, is shown in the right panel of Figure 1. Qualitatively, the two alternative model fits to the 5 GHz light curves are very similar; the free-free absorption model fits the optically thick rise at 1.4 GHz better, but the synchrotron self-absorption model better matches the peak. Quantitatively, the fits give comparable values of $`\chi ^2`$ and we cannot distinguish between the two absorption mechanisms. The $`\chi ^2`$ values of both of these fits (calculated using the data up to day 3380) are very similar to that of Weiler et al. (1992).

Chevalier (1998) estimated the minimum radius of the SN at the peak of the light curve if synchrotron self-absorption was the dominant absorption mechanism, and argued synchrotron self-absorption was probably not significant in SN1980K because the implied expansion velocity is much smaller than observed. However, this argument assumes equipartition of energy between the relativistic particle energy density and the magnetic energy density (i.e., $`\zeta _{rel}=\zeta _B`$). The model that we fit allows $`\zeta _B`$ and $`\zeta _{rel}`$ to take values which are far from equipartition, in which case the arguments of Chevalier (1998) imply an expansion velocity for SN1980K which is similar to the measured velocity and synchrotron self-absorption is not ruled out.

Although neither the fits to the light curves, nor the expansion velocity arguments, determine whether the observed optically-thick rise of SN1980K is dominated by free-free absorption or by synchrotron self-absorption, the spatial distribution of surface brightness in these two cases is quite different. Our method of calculation of the radio emission allows us to generate synthetic images, and two examples are shown in Figure 2. These synthetic radio images and the corresponding radial profiles of the surface brightness shown below them are for the same simulations used to produce the light curves in Figure 1. They show the emission at 1.48 GHz on day 172 after explosion, when the source is still optically thick at this frequency. The free-free absorption dominated image is a distinct ring which is brightest at a radius roughly half the total radius of the source. In contrast, the synchrotron self-absorption dominated image is a uniformly-bright disk. These differences were predicted analytically by Marscher (1985) and would be easily distinguishable (in principle) via high resolution images at epochs very soon after the explosion of a radio supernova. Unfortunately, SN1980K could not have been imaged at sufficient resolution to reveal these differences, and even the two closest supernovae detected to date, SN1993J and SN1987A, were not imaged with sufficient resolution at early times to make this distinction.

### 4.2 SN1993J

SN1993J, in the nearby galaxy M81 at a distance of 3.6 Mpc (Freedman et al., 1994), is one of the closest SNe to explode in modern history. It is believed to be a Type II SN that may have been a member of a binary system (Höflich, Langer & Duschinger, 1993; Podsiadlowski et al., 1993; Woosley et al., 1994). The radio evolution of SN1993J has been followed in detail throughout most of its lifetime, and thanks to its proximity it has been possible to image the evolving RSN extensively with VLBI (Marcaide et al., 1997; Bartel et al., 2000). Marcaide et al. (1997) found that the expansion of the radio source was consistent with a self-similar $`rt^m`$ dependence. However, Bartel et al. (2000), using careful astrometry on a larger sample of observations obtained over a longer time period, have reported that the expansion is not self-similar. The simple mini-shell model cannot reproduce the rapid turn-on and subsequent relatively slow optically-thick rise of the radio emission from SN1993J (Van Dyk et al., 1994b). Furthermore, since the evolution of the SN1993J radio source is not self-similar, a direct application of the mini-shell model is not appropriate. It is therefore an ideal test case for our method of calculation of the radio emission.

Fransson & Björnsson (1998) assumed a power-law ejecta density profile for SN1993J and produced simulated radio light curves which compare well with the early data. However, Nomoto, Suzuki, Iwamoto and other coworkers have suggested that the optical and X-ray evolution of SN1993J implies that the ejecta distribution is more complicated. Iwamoto et al. (1997) computed the 2D hydrodynamical evolution of SN1993J and calculated the optical light curve. They found a series of ejecta models that gave an acceptable fit to the observed light curve if the extent of the mixing for nickel was assumed to be larger than seen in their 2D simulations. Similar models were used by Suzuki & Nomoto (1995) to fit the X-ray emission with reasonable results.

We use a specific ejecta model supplied by this group – model 4H47 (courtesy of Ken Nomoto, private communication) – and calculate the hydrodynamics and radio emission from the model supernova. The model ejecta mass is 3.12 M, the pre-SN progenitor star has a radius of 350 R and an envelope mass of 0.47 M. The fraction of helium in the envelope mass was about 0.79. The explosion energy is 10<sup>51</sup> ergs.

The density distribution of the CSM surrounding SN1993J is poorly constrained by observations. Many of the models that reproduce the optical, X-ray and radio emission require a density profile that falls off more slowly than $`r^2`$ (eg. Suzuki & Nomoto 1995; Fransson, Lundqvist & Chevalier 1996; Nomoto and Suzuki 1998; Van Dyk et al. 1998). Fransson & Björnsson (1998) however find the standard $`r^2`$ distribution to be adequate in their self similar modeling of the radio light curve. We therefore consider three CSM distributions, with density varying as $`r^{1.5},r^{1.7}`$, and $`r^2`$.

Figure 3 shows the density profiles at four timesteps during the hydrodynamic evolution for a CSM density profile of $`\rho _{CSM}r^{1.7}`$. The interaction of the ejecta with the CSM leads to the formation of an outer shock which expands through the CSM, and a reverse shock that propagates back into the ejecta, separated by a contact discontinuity. The 4H47 ejecta model contains local maxima, two of which are indicated by arrows in the top left panel of Figure 3. The inclusion of a cooling function in our hydrodynamical simulation means that these peaks remain well defined as the ejecta expand. There is a significant density jump just ahead of the outermost local density maximum (labelled 4). The interactions of this density jump and the local peak with the reverse shock occur at ages of around 1300 days and 2300 days respectively, the latter of which is close to the present epoch. These interactions have a dramatic effect on the pressure profiles, and hence on the expansion of the supernova and the radio emission.

The structure in the ejecta profile leads to a distinctly non-self-similar evolution. We define the expansion parameter $`M(t)=\mathrm{d}\mathrm{ln}r_s/\mathrm{d}\mathrm{ln}t`$ so that for self-similar expansion $`M(t)=m`$, a constant. Figure 4 shows the variation of the expansion parameter with time during the early stages of the simulation of the interaction of the 4H47 ejecta with the $`r^{1.7}`$ CSM. The decrease of the expansion parameter with time over the first 2300 days, representing an increasing rate of deceleration of the supernova shock, is consistent with the observations of Bartel et al. (2000) who fit an expansion parameter $`M(t)0.937`$ up to 306 days and $`M(t)0.77`$ between 582 and 1893 days (shown as dashed lines in figure 4). At around day 2300 the interaction of the outermost local maximum in the ejecta density profile with the reverse shock leads to a significant increase in the pressure between the reverse and forward shocks. This has the effect of accelerating the forward shock and leads to the rapid increase in the expansion parameter seen in Figure 4 around days 2300–4000.

#### 4.2.1 Radio Supernova 1993J

The radio evolution of SN1993J is most extensively followed at 8.4 GHz. We investigate the radio emission calculated from the hydrodynamical simulation of the 4H47 ejecta model by first trying to fit the observed 8.4 GHz light curve, and then considering the data at other frequencies.

None of the fits obtained using only free-free absorption, or only synchrotron self-absorption, provides a good match to the early data for any of the CSM profiles considered. A much better fit is possible when both free-free absorption and synchrotron self-absorption are included. Figure 5 shows the best fits obtained for $`\rho _{CSM}r^2`$ (solid line), $`\rho _{CSM}r^{1.7}`$ (dashed line), and $`\rho _{CSM}r^{1.5}`$ (dash-dotted line). All provide an adequate fit to the rising phase, an indication of the flexibility provided by the two different absorption processes. The differences between the models become apparent near the peak of the light curves, and increase as the flux density decreases. In the declining phase neither absorption process is important, and the only model parameters that affect the light curve for a given hydrodynamical simulation are the spectral index $`\gamma `$ (or $`\alpha `$), $`\mathrm{\Delta }`$ and the overall scaling. The fit for a CSM density profile of $`\rho _{CSM}r^2`$ is very good at times up to about day 200, soon after the peak flux density, but at later times the model flux density drops off more rapidly than the observations. The CSM density profile of $`\rho _{CSM}r^{1.7}`$ provides the best fit up to the present epoch. The decline of the model for $`\rho _{CSM}r^{1.5}`$ is considerably slower than the observed dependence, and the peak is much broader than is observed.

For our best fit model, $`T_{CSM}=7\times 10^4\mathrm{K}`$, which is smaller than the value $`2\times 10^5\mathrm{K}`$ quoted by Fransson & Björnsson (1998), and than the values suggested by Fransson et al. (1996) who considered the CSM temperature in detail. For $`T_{CSM}<7\times 10^4\mathrm{K}`$ the free-free opacity in our model is too high at early times.

All the model curves exhibit a temporary flattening in the decline, at times varying from day $`1000`$ for $`\rho _{CSM}r^2`$, to day $`2000`$ for $`\rho _{CSM}r^{1.5}`$. This is directly attributable to the interaction of the outer density jump and local maximum in the 4H47 model ejecta with the reverse supernova shock, which temporarily halts the decline in the pressure between the forward and reverse supernova shocks. The observed 8.4 GHz data show a similar flattening, starting around day 1500 and extending to at least day 2500. Such features in observed RSN light curves have previously been attributed to interactions between the SN shock with density structures in the CSM. Our simulations suggest that light curve features may also arise from interactions between the SN shock and density structures in the ejecta, but it will be difficult to distinguish between these two possibilities solely on the basis of the light curves.

Figure 6 shows synthetic images and surface brightness profiles of SN1993J at 8.4 GHz at 20, 86, 390 and 1301 days after explosion, for the the same $`\rho _{CSM}r^{1.7}`$ simulation as for Figures 4 and 5. On day 20 – panel (a) – the source is optically thick and synchrotron self-absorption dominates. As expected, this image is similar in character to the right panel in Figure 2. On day 86 – panel (b) – SN1993J is becoming optically thin. On day 390 – panel (c) – the source and the CSM are optically thin and the flux density at all frequencies is declining with time. On day 1301 – panel (d) – the model 8.4 GHz light curve has flattened as a result of the interaction of the density jump in the ejecta with the reverse shock. As expected, the ring seen in Figure 6 a, b and c is bound on the inside by the reverse shock, and on the outside by the forward shock. The VLBI images of SN1993J of Bartel et al. (2000) reveal a clumpy ring of emission that when averaged azimuthally is qualitatively similar to the simulated images shown in Figure 6. The radius and thickness of the bright ring of the synthetic images are similar to those of the observed data, but the radial profile of the surface brightness in the synthetic images is quantitatively different from that observed, a point we discuss in more detail in §4.2.2. Panel (d) shows a distinct brightening on the inner edge of the ring as a result of the pressure increase associated with the interaction of a density discontinuity in the ejecta with the reverse SN shock, which also causes the flattening seen in the model light curve. If the observed flattening were due to an interaction between the forward SN shock and a density structure in the CSM, the result would be a brightening on the outer edge of the emission ring. Such differences may ultimately provide a means of distinguishing between SN shock interactions with structures in the ejecta and those in the CSM. However images of SN1993J obtained to date are not of sufficiently high resolution to show these kinds of differences, and the clumpy nature of the images complicates their interpretation.

Although Figure 5 shows the model flux densities for only a single frequency, the model parameters determine the flux density at all frequencies and at all epochs. There are useful observations of the evolving radio flux density from SN1993J at four frequencies, 4.9, 8.4, 15.0 and 22 GHz. Figure 7 shows the observed data and the model flux densities for the parameters determined by the best fit in Figure 5 with $`\rho _{CSM}r^{1.7}`$. The limitations of the model are immediately apparent. Despite the relatively good fit of the model to the data at 8.4 GHz, the fits at the other frequencies are relatively poor at all times after the flux density has increased to a significant fraction of the maximum. This may be attributed to the model itself, not the particular choice of parameters used. The fundamental model assumptions that the magnetic energy density is proportional to the thermal energy density, and that the relativistic electron energy distribution is a (fixed) power law with an energy density which is also proportional to the thermal energy density, provide an unbreakable link between the frequency and time dependence of the model flux density. At late times when the source and the CSM are optically thin both the frequency and the time dependence of the model flux density are completely determined by the model parameter $`\alpha `$ (or $`\gamma `$) for a given hydrodynamical simulation. The other model parameters, $`T_{CSM}`$, $`\mathrm{\Delta }`$, and the overall scaling, have no effect on the frequency or temporal dependence at late times. In particular, the model fit to the 8.4 GHz data requires the choice $`\gamma =2.1`$ to fit the observed temporal evolution. This value of $`\gamma `$ implies $`\alpha =0.55`$, which determines the frequency dependence of the model. If the simulation was self similar the model flux density in the optically-thin regime would therefore have a frequency dependence of $`S\nu ^{0.55}`$. A simple fit to the observed spectra at times after about 150 days indicates $`S\nu ^{0.8}`$, a difference which is much too large to be adequately fitted by the simulation.

#### 4.2.2 Internal Free-Free Absorption

Figure 8 shows the azimuthally averaged profile of the observed 8.4 GHz emission from SN1993J (solid line) plotted with the profile of the best fit model simulation at approximately the same age (dashed line). The diameter and width of the bright ring of emission in the observed and synthetic images are very similar. However the ratio of the surface brightness at the center of the image and at the brightest point on the ring is significantly smaller in the observed data than in the simulated images. The most likely reason for this difference is that internal free-free absorption, i.e. free-free absorption occurring within the SN, has not been included in our simulations. For the first tens of years the temperatures and densities internal to the reverse shock are such that most, if not all, the emission that passes through the interior is absorbed, significantly decreasing the surface brightness seen from the center of the SN. The dotted line in Figure 8 shows the surface brightness profile of the simulated image when internal free-free absorption is included. It matches the observed profile well, lying within the RMS scatter at all radii, and we conclude that internal free-free absorption is important in SN1993J. Figure 9 shows the comparison of the simulated 8.4 GHz light curves for SN1993J with and without internal free-free absorption. For the the first 8 or so years the light curves have a similar shape, not changing our conclusions in the previous sections of this paper. After year 8 there is a slight flattening in the light curve until year 55 as the internal free-free absorption becomes unimportant due to the decrease in the density inside the reverse shock as the SN expands. During this time (starting at around the present epoch) the ratio of the surface brightness at the center of the SN to that at the peak of the bright ring of the simulated SN increases steadily, and we therefore expect the same to occur in the VLBI images of SN1993J. At higher frequencies the internal free-free absorption becomes unimportant at earlier epochs. These conclusions and predictions regarding the internal free-free absorption are independent of the assumptions used to calculate the synchrotron radiation. If the surface brightness ratio does not increase as predicted, this has implications for the internal density and temperature of the SN.

If other RSNe have a similar level of internal free-free absorption to that implied for SN1993J, radio emission from an embedded pulsar would be completely absorbed until the SN ejecta have expanded sufficiently for the internal free-free optical depth to become negligible. Our model for SN1993J suggests that the optical depth at 1 GHz falls to unity around 150 years (50 years for 8.4 GHz) after the explosion. Detection of a radio pulsar embedded in a similar nearby SN would be impossible before that time. For most extragalactic SN this is a moot point since radio pulsars at extragalactic distances are unlikely to be detectable even if not absorbed. Modeling the radio plerion created by the pulsar, Reynolds & Chevalier (1984) also predict that one could not observe a pulsar for 200 years at radio frequencies due to internal free-free absorption.

## 5 Discussion and conclusions

Our results indicate the viability of the calculation of supernova radio emission from high-resolution hydrodynamical simulations. The calculations presented in this paper retain much of the essential physics of the mini-shell model, particularly the connection between the model parameters and physical parameters describing the shock, the magnetic field, and the radiating electron distribution. When the ejecta and the CSM have simple power-law density profiles our generalized treatment reproduces the results of the analytic self-similar model. The advantage of our new method is that it provides a means to calculate the radio emission in cases where the expansion in not self-similar, which we demonstrate with calculations for an ejecta model for SN1993J which has a complicated profile motivated by detailed modeling of the optical and x-ray emission.

Attempts to fit the observed time and frequency dependence of the radio emission from SN1993J, which is not well modeled by the analytic mini-shell, are encouraging but not yet completely successful. Fits to the evolving radio flux density of SN1993J at a single frequency indicate that neither synchrotron self-absorption nor free-free absorption alone can produce the observed increase, but that when both effects are included the rise in the flux density can be well modeled for a range of parameters. This conclusion reinforces the results of Fransson & Björnsson (1998) and Chevalier (1998) who argued that both synchrotron self-absorption and free-free absorption were important for SN1993J. Most other RSNe first become detectable at much later epochs than did SN1993J, and their subsequent radio evolution can be fitted adequately using only free-free absorption. This may be because in other SNe the conditions are such that the optical depth due to synchrotron self-absorption is negligible by the time the free-free optical depth through the CSM becomes small enough for the SN to be observed. An alternative is that the sparse sampling of the optically-thick increase of most RSNe simply precludes a determination of the relative importance of synchrotron self-absorption and free-free absorption, as shown in the SN1980K model fits of §4.1.

Our simulations show that the different absorption mechanisms produce quite different spatial surface brightness distributions at very early epochs. They also show that features in the radio light curves of SNe may result from the interaction of density structures in the supernova ejecta with the reverse supernova shock. Such interactions may be detectable through visible signatures in the simulated images, such as a brightening at one edge of the ring of radio emission. VLBI imaging to date does not have sufficient resolution to reveal such effects, but they may prove fruitful goals for future studies. However comparison of the profiles of VLBI images of SN1993J at 8.4 GHz (Bartel et al., 2000) with our simulated images implies that internal free-free absorption is important now and will remain so for a few tens of years. We predict that the ratio of the surface brightness at the center of the SN to that at the bright ring should increase, starting at about the present epoch, due to the gradual decrease of the internal free-free absorption as the SN expands. A similar level of internal free-free absorption in a nearby SN would absorb the emission from an embedded radio pulsar, rendering it undetectable at 1 GHz until around 150 years after explosion.

The relatively slow decline of the SN1993J flux density at a given frequency may indicate that the magnetic field in the emitting region declines more slowly than the scaling with thermal energy density implies (e.g. Chevalier 1998), or that the density of the CSM encountered by the expanding supernova shock to date declines more slowly than the $`r^2`$ dependence expected for a constant wind. We find a good fit for the latter case with $`\rho _{CSM}r^{1.7}`$. A density distribution less steep than $`r^2`$ has also been proposed by a number of other researchers (Suzuki & Nomoto, 1995; Fransson et al., 1996; Van Dyk et al., 1998). Such a dependence could result from a time variation in the mass-loss rate of the wind emitted by the pre-SN progenitor star, or a variation in the wind velocity, or both. If the wind was steady during an earlier phase of evolution of the progenitor, then at some distance further away from the star one would expect the CSM density distribution to revert to an $`r^2`$ profile. Suzuki & Nomoto (1995) and Fransson et al. (1996) suggest that, in order to match the X-ray and optical light curves and spectra, the density distribution must change close to the star, within a radius of about 10<sup>16</sup> cm. However, there is no break in the radio light curves at an early stage that would indicate a sudden change in the density distribution, and in our simulations we assume that the initial CSM density distribution persists throughout the evolution.

Finally, we find that the generalization of the mini-shell model which results from a realistic treatment of the hydrodynamics is insufficient to account for the observed evolution and flux density dependence of SN1993J. The link between the frequency and time dependencies of the model flux density due to the assumed energy density scalings remains implicit in the implementation of the generalized model used here, which therefore doesn’t have sufficient flexibility to fit the data. The inclusion of any one of a number of aspects of the physics which is not addressed by the model as presented here may break this nexus. An obvious next step is to investigate alternative scalings for the magnetic field and/or the relativistic electron distribution (see for example Chevalier 1996 and Chevalier 1998). An alternative we are presently pursuing is to include a physically-justifiable calculation of the time-dependent acceleration of the radiating electrons resulting from the calculated evolution of the expanding supernova shock, along the lines of the treatment by Duffy, Ball, & Kirk (1995).

We thank the VLBI collaboration led by Norbert Bartel, and in particular Michael Bietenholz and Michael Rupen, for providing us with the flux density measurements of SN1993J prior to publication and with access to their VLBI images, and Ken-Ichi Nomoto for providing the 4H47 ejecta model. AJM thanks Michael Rupen and Mark Wardle for suggestions which greatly improved this paper. VVD thanks Roger Chevalier for helpful discussions. We are grateful to Kurt Weiler and colleagues for providing RSN data sets on the web page http://rsd-www.nrl.navy.mil/7214/weiler/sne-home.html. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. VVD’s research was partially supported by a grant from NASA administered by the American Astronomical Society.


format: text

text:
# Dynamics of Coupling Functions in Globally Coupled Maps: Size, Periodicity and Stability of Clusters

## I Introduction

There has been much interest in the study of the collective dynamics of chaotic systems subjected to global interactions. Such systems arise naturally in the description of arrays of Josephson junctions, charge density waves, multimode lasers, neural dynamics, evolutionary, chemical and social networks . The globally coupled map (GCM) lattice  constitutes a prototype model for such global-coupling dynamics. It has recently been argued that GCM systems yield universal classes of collective phenomena . Specifically, a GCM system can exhibit a variety of collective behaviors such as clustering (i.e., the formation of differentiated subsets of synchronized elements in the network) ; non-statistical properties in the fluctuations of the mean field of the ensemble ; global quasiperiodic motion ; and different collective phases depending on the parameters of the system . It has been shown that a GCM system is closely related to a single map subjected to an external drive and that this analogy may be used to describe the emergence of clusters in GCM systems in geometrical terms .

In particular, the phenomenon of clustering is relevant as it can provide a simple mechanism for segregation, ordering and onset of differentiation of elements in many physical and biological systems. In addition to GCM systems, dynamical clustering has also been found in a globally coupled Rössler oscillators , neural networks , and coupled biochemical reactions . The interest in this phenomenon has recently grown, since dynamical clusters have been observed experimentally in an array of electrochemical oscillators interacting through a global coupling .

In this paper, we investigate the process of cluster formation in general globally coupled map systems by focusing on the dynamics of their global coupling functions. In most studies on GCM systems, the mean field of the network has been used as the global coupling function. Here, we study GCM systems subjected to different global coupling functions and show how they can be analyzed under a common framework. We investigate how the distribution of elements among a few clusters and their periodicities depend on the functional form of the global coupling.

Section II contains a description of the dynamics of different global coupling functions in GCM systems and a calculation of the possible periodicities and cluster sizes when two clusters emerge in these systems. The driven map analogy is employed in Sec. III to interpret the clustering behavior of GCM systems. In Sec. IV the dynamical properties of periodic clusters in systems exhibiting a constant global coupling are predicted; and for a particular family of global coupling functions, the stability condition for these clustered states is derived in an Appendix. Conclusions are presented in Sec. V.

## II Dynamics of global coupling functions

Consider a general globally coupled map system

$$x_{t+1}(i)=(1ϵ)f\left(x_t(i)\right)+ϵH(x_t(1),x_t(2),\mathrm{},x_t(N)),$$
(1)
where $`x_t(i)`$ gives the state of the element $`i`$ $`(i=1,2,\mathrm{},N)`$ at discrete time $`t`$; $`N`$ is the size of the system; $`ϵ`$ is the coupling parameter; $`f(x)`$ describes the (nonlinear) local dynamics, which in the present article is chosen to be the quadratic map $`f(x)=1rx^2`$; and $`H_t(x_t(i),x_t(2),\mathrm{},x_t(N))`$ is the global coupling function. We shall consider a general class of global coupling functions of $`N`$ variables such that $`H(\mathrm{},x_t(i),\mathrm{},x_t(j),\mathrm{})=H(\mathrm{},x_t(j),\mathrm{},x_t(i),\mathrm{})`$, $`i,j`$; that is, $`H`$ is assumed to be invariant to argument permutations. This property of the coupling function ensures that, at any time, each element of the globally coupled system is subjected to the same influence of the coupling term. Some examples of coupling functions belonging to this class are

$`H=`$ $`x`$ $`={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}x_t(i);`$ (2)
$`H=`$ $`f`$ $`={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}f(x_t(i));`$ (3)
$`H=`$ $`\mathrm{\Delta }x`$ $`=\sqrt{{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(x_t(i)x\right)^2};`$ (4)
$`H=`$ $`\overline{x}`$ $`={\displaystyle \underset{i=1}{\overset{N}{}}}|x_t(i)|^{\frac{1}{N}}.`$ (5)
The first two examples correspond to forward and backward mean field coupling, respectively, and they have been widely used in GCM studies. The third global coupling function is the usual dispersion or mean square deviation of $`N`$ variables, and it may describe systems whose elements do not interact when they are synchronized. This kind of global interaction might be relevant in some biological or social systems where the members of a community are driven by their deviations from the mean behavior. The last example is the geometric mean. This type of multiplicative coupling occurs, for instance, in a system of $`N`$ sequential amplifiers where the gain of element $`i`$ is a function of the magnitude of its state $`x_t(i)`$, and $`H`$ is proportional to the total gain of the system. Many statistical functions of $`N`$ variables share the property of invariance under argument permutations and they could as well be taken as global coupling functions in GCM systems given by Eq.(1).

For some range of its parameters the GCM system in Eq.(1) reaches an asymptotic collective behavior characterized by the segregation of the elements into $`K`$ clusters, each exhibiting a period $`P`$, where the $`k\text{th}`$ cluster has a number $`N_k`$ of elements, with $`_{k=1}^KN_k=N`$. The fraction $`p_k`$ of elements in the $`k\text{th}`$ cluster is $`p_k=N_k/N`$. The evolution of the $`k\text{th}`$ cluster may be described by a variable $`\chi _t(k)`$ which gives the common state of the $`N_k`$ elements belonging to this cluster at time $`t`$. The periodic orbit adopted by the state $`\chi _t(k)`$ of the $`k\text{th}`$ cluster can be expressed as a sequence of $`P`$ values $`[\chi _1(k),\chi _2(k),\mathrm{},\chi _P(k)]`$. The specific partition $`\{p_1,p_2,\mathrm{},p_K\}`$ into $`K`$ clusters and the specific values taken by the periodic orbit of each cluster depend on initial conditions and parameters of the system.

When a GCM system falls into $`K`$ periodic clusters, the coupling function $`H`$ also shows a periodic motion. As an illustration of this behavior, Fig. 1(a) shows a typical situation in which a GCM system, with global coupling function $`H=\overline{x}`$, displays two clusters, each in period two. In this case, the coupling function follows a period-two motion.

Collective states consisting of two clusters have recently been observed in an experimental array of globally interacting chemical elements . This clustered collective state is interesting to analyze in globally coupled maps. In this simple situation, there is a fraction $`p`$ of elements in one cluster and a fraction $`(1p)`$ in the other cluster. Thus the global coupling functions from Eqs.(2)-(5) simplify to:

$`H=`$ $`x`$ $`=p\chi _t(1)+(1p)\chi _t(2),`$ (6)
$`H=`$ $`f`$ $`=pf(\chi _t(2))+(1p)f(\chi _t(2)),`$ (7)
$`H=`$ $`\mathrm{\Delta }x`$ $`=\sqrt{p(1p)(\chi _t(2)\chi _t(1))^2},`$ (8)
$`H=`$ $`\overline{x}`$ $`=|\chi _t(1)|^p|\chi _t(2)|^{(1p)}.`$ (9)
When the GCM system reaches a two-cluster state, the dynamics of the system reduces to the two coupled maps

$`\chi _{t+1}(1)`$ $`=`$ $`(1ϵ)f(\chi _t(1))+ϵH,`$ (10)
$`\chi _{t+1}(2)`$ $`=`$ $`(1ϵ)f(\chi _t(2))+ϵH;`$ (11)
where $`H=H(\chi _t(1),\chi _t(2),p)`$ is the reduced, two-cluster coupling function, as in Eqs.(6)-(9). If both cluster states $`\chi _t(1)`$ and $`\chi _t(2)`$ fall in period-two orbits, the coupling function $`H`$ follows, in general, a period-two motion, as shown in Fig. 1(a), although $`H`$ may become constant in some circumstances (see Fig. 1(b) and Sec. IV). Let $`H_1`$ and $`H_2`$ be the values adopted alternatively in time by $`H`$ in its period-two orbit for a given partition $`\{p,1p\}`$, as indicated on Fig. 1(a). The values $`H_1`$ and $`H_2`$ depend on i) the functional form of $`H`$; ii) the parameters of the GCM system Eq.(1); and iii) the fraction $`p`$.

Consider then several GCM systems with the same parameters $`r=1.7`$ and $`ϵ=0.2`$ but with different coupling functions such as those in Eqs. (2)-(5). Two clusters in period two can emerge in each of these systems for some range of the fraction $`p`$. The resulting asymptotic orbits $`[H_1,H_2]`$ of the respective coupling functions are shown in Fig. 2(a) as $`p`$ varies, giving rise to a curve in each case. Note that each function $`H`$ possesses period-two orbits only for a limited range of the fraction $`p`$. Figure 2(b) is a magnification of Fig. 2(a) which shows that the dynamics of the backward and the mean field coupling functions become equal for $`p=0.5`$; i.e., when the two clusters have equal sizes. In this case, both coupling functions reach the constant value $`H_1=H_2=0.365`$. Notice that the dispersion coupling function, $`H=\mathrm{\Delta }x`$, only displays states with $`H_1=H_2`$; that is, even when the two clusters in period two may have different sizes, this particular global coupling always reaches a stationary value. The curves for the other coupling functions are symmetrical with respect to the diagonal in Fig. 2(b), which they cross for $`p=0.5`$. On the diagonal, the coupling functions are constant and the two clusters evolve out of phase with respect to each other (Sec. IV).

Note also that the different global coupling functions perform a period-two motion only on a restricted region of the plane $`(H_1,H_2)`$. It will be shown in Sec. III that period-two orbits of any permutable coupling function will fall within the dashed contour in Fig. 2.

In general, a coupling function $`H`$ of a GCM system in a collective state of two clusters can reach various asymptotic periodic orbits for appropriate initial conditions. Each Fig. 3(a) to 3(d) shows the regions on the space of parameters $`(p,ϵ)`$ for which a coupling function $`H`$ of a GCM in a two-cluster state displays different periodic motions. The local parameter is fixed at $`r=2`$. Figures 3(a) and 3(b) correspond to the backward and forward mean field coupling, respectively. Note the very different distributions of periodic regions for the coupling functions in Figs. 3(a) and 3(b). It should be noticed that, besides the collective periodic states for two clusters shown in Figs. 3(a)-3(d), there can exist other states in a GCM system consisting of more than two periodic clusters for the same values of the parameters $`r`$ and $`ϵ`$, but corresponding to different initial conditions.

The inverse problem of determining the global coupling function in experimental systems is relevant since in general the specific functional form of the acting coupling is not known. This can be a complicated problem because, in addition, the exact form of the local dynamics may not be extracted in most situations. However, if the local dynamics is known some insight on the function $`H`$ of a globally coupled system can be gained within the framework presented here. For example, in the case of a dynamical system showing two period-two clusters with partition $`p`$, the resulting asymptotic orbit $`[H_1,H_2]`$ can be obtained by measuring the cluster orbits and using Eqs. (10) and (11). For different realizations of partition $`p`$, the curve $`[H_1,H_2]`$ can be drawn as a function of $`p`$ on the plane $`(H_1,H_2)`$, and compared with curves $`[H_1,H_2]`$ corresponding to known functional forms $`H`$ such as those in Fig. 2.

## III Driven map analogy

As shown in Ref. , the clustering behavior of GCM systems can be analyzed through its analogy with the dynamics of a single map subjected to an external drive. In order to interpret the results of the preceding section, let us consider an associated driven map

$$s_{t+1}=(1ϵ)f(s_t)+ϵL_t,$$
(12)
where $`s_t`$ is the state of the map at discrete time $`t`$, $`f(s_t)`$ is the same local dynamics as in Eq.(1), and $`L_t`$ is an external driving term assumed to be periodic with period $`P`$. We denote the sequence of $`P`$ values adopted by the periodic drive $`L_t`$ by $`[L_1,L_2,\mathrm{},L_P]`$. The analogy between a GCM system and a driven map arises because in the former system (Eq. (1)) all the elements are affected by the global coupling function in exactly the same way at all times, and therefore the behavior of any element $`x_t(i)`$ in the GCM is equivalent to the behavior of a single driven map (Eq. (12)) with $`L_t=H`$ and initial condition $`s_o=x_o(i)`$. Additionally, if a GCM system reaches a clustered, periodic collective state, its corresponding coupling function $`H`$ follows in general a periodic motion. Thus the associated driven map (Eq. (12)) with a periodic drive $`L_t`$ should display a behavior similar to that of an element belonging to a periodic cluster in the GCM system. In particular, periodic drives resulting in periodic orbits of $`s_t`$ in Eq. (12) may be employed to predict the emergence of clustered, periodic states in a GCM (Eq. (1)), regardless of the specific functional form of the global coupling $`H`$ and without doing direct simulations on the entire GCM system.

The driven map is multistable; i.e., there can exist several attractors for the same parameter values $`r`$ and $`ϵ`$. Specifically, for a given periodic drive $`[L_1,L_2,\mathrm{},L_P]`$, the map $`s_t`$ may reach a number of distinct asymptotic responses $`\overline{s}_t(j)`$; $`(j=1,2,\mathrm{},J)`$, all with the same period, and depending on initial conditions. The orbits of $`s_t`$ with period $`M`$ can be expressed as a sequence of values $`[\overline{s}_1(j),\overline{s}_2(j),\mathrm{},\overline{s}_M(j)]`$. The correspondence between a GCM system (Eq. (1)) in a state of $`K`$ clusters with period $`P`$ and its associated driven map (Eq. (12)) can be established when $`M=P`$ and $`J=K`$.

Using this analogy, the main features in Fig. 2 can now be explained. In terms of a driven map subjected to a period-two drive $`L_t=[L_1,L_2]`$ and the same parameters $`r=1.7`$ and $`ϵ=0.2`$ as in the GCM systems considered in Fig. 2, the bounded region on that figure contains the values of $`L_1`$ and $`L_2`$ for which the driven map $`s_t`$ just possesses two distinct asymptotic orbits $`(J=2)`$ of period two $`(M=2)`$, denoted by $`[\overline{s}_1(1),\overline{s}_2(1)]`$ and $`[\overline{s}_1(2),\overline{s}_2(2)]`$. For values of $`L_1`$ and $`L_2`$ outside this diamond shaped region, the map $`s_t`$ may also reach a number of asymptotic periodic orbits but none with both $`J=2`$ or $`M=2`$, at least on the interval $`1L_1,L_21`$. Because of the analogy drawn above, all period-two motions $`[H_1,H_2]`$ of permutable coupling functions in GCM systems given by Eq. (1) with $`r=1.7,ϵ=0.2`$ and displaying two clusters in period two will fall on this bounded region of the plane $`(H_1,H_2)`$. Equivalently, a collective state of two clusters in period two can emerge in a GCM system only if its global coupling function has an orbit $`[H_1,H_2]`$ with values $`H_1=L_1`$ and $`H_2=L_2`$ contained within the bounded region of Fig. 2. The boundaries of that region vary as the parameters $`r`$ and $`ϵ`$ are changed. It is out of the scope of the present work to establish how the bounded region for two clusters in period two, as well as other regions for different clustered, periodic states, depend on parameters. However, it is worth noticing that this dependence can serve to characterize GCM systems with permutable coupling functions, and that this characterization can be obtained by the sole use of an associated driven map.

The analogy between a GCM system with a given coupling function $`H`$ in a two-cluster state and an associated driven map can be carried further by defining an associated coupling function

$$\mathrm{\Theta }_H=H(\overline{s}_t(1),\overline{s}_t(2),p);$$
(13)
that is, $`\mathrm{\Theta }_H`$ is similar to the reduced two-cluster coupling function, such as Eqs.(6)-(9) but with the arguments $`\chi _t(1)`$ and $`\chi _t(2)`$ from the cluster trajectories replaced by the driven map orbits $`\overline{s}_t(1)`$ and $`\overline{s}_t(2)`$. The associated coupling function links the dynamics of clusters in a GCM system to the dynamics of a single associated driven map. The function $`\mathrm{\Theta }_H`$ depends on the functional form of the coupling function $`H`$, on the partition $`\{p,1p\}`$ among the two clusters in the GCM, and on the orbits $`\overline{s}_t(1)`$ and $`\overline{s}_t(2)`$, which themselves are function of the period-two drive $`[L_1,L_2]`$ and the parameters $`r`$ and $`ϵ`$. Thus for fixed $`r`$ and $`ϵ`$, and a given $`H`$, we have $`\mathrm{\Theta }_H=\mathrm{\Theta }_H(L_1,L_2,p)`$. The main point is that, an equivalence between a GCM system Eq. (1) in a two-cluster, period-two state, and an associated driven map, Eq. (12), with a period-two drive occurs when the following conditions are fulfilled

$`\mathrm{\Theta }_H(\overline{s}_1(1),\overline{s}_1(2),p)`$ $`=`$ $`L_1,`$ (14)
$`\mathrm{\Theta }_H(\overline{s}_2(1),\overline{s}_2(2),p)`$ $`=`$ $`L_2.`$ (15)
Eqs. (14)-(15) constitute a set of two nonlinear equations for $`L_1`$ and $`L_2`$, for a given $`p`$. The solution $`[L_1^{},L_2^{}]`$ of Eqs. (14)-(15) predicts that the GCM possesses a state characterized by the coupling function motion $`[H_1=L_1^{},H_2=L_2^{}]`$ and cluster orbits $`[\chi _1(1),\chi _2(1)]=[\overline{s}_1(1),\overline{s}_2(1)]`$, $`[\chi _1(2),\chi _2(2)]=[\overline{s}_1(2),\overline{s}_2(2)]`$. The succession of solutions $`[L_1^{},L_2^{}]`$, as $`p`$ varies, yields the curve corresponding to a given $`H`$ in Fig. 2. Thus, each curve on the plane $`(H_1,H_2)`$ is parameterized by the fraction $`p`$. Moreover, there exist solutions $`[L_1^{},L_2^{}]`$ to Eqs. (14)-(15) only for an interval of $`p`$. Therefore, the curves in Fig. 2 can, in principle, be calculated a priori by using an associated driven map and just the specific functional form of $`H`$ in each case. The range of possible cluster sizes, described by the values of the fraction $`p`$ for which exist solutions to Eqs. (14)-(15), can also be predicted by this method. Similarly, the regions of period two in Figs. 3(a)-3(d) can be obtained by varying the parameter $`ϵ`$ and calculating the interval of $`p`$ for which Eqs. (14)-(15) have solutions.

## IV Prediction and stability of clusters in systems with constant global coupling

Another simple clustered collective state in GCM systems occurs when the coupling function $`H`$ remains constant in time, i.e., $`H=C`$. This behavior may take place in a GCM system with a permutable coupling function when $`K`$ clusters, each having $`N/K`$ elements and period $`K`$, are evolving with shifted phases in order to yield a constant value for $`H`$. That is, if the periodic orbits of $`K`$ identical-size clusters are cyclically permuting in time, the resulting $`H`$ becomes constant. For those collective states, the behavior of any of such clusters in the GCM system can be emulated by an associated driven map subjected to a constant forcing $`L_t=C`$ . In the case of a GCM displaying two equal size clusters in period two, this situation corresponds to the intersection of $`H`$ with the diagonal in Fig. 2. The cluster orbits are then related as $`\chi _t(1)=[\chi _1(1),\chi _2(1)]=[a,b]`$ and $`\chi _t(2)=[\chi _1(2),\chi _2(2)]=[b,a]`$. On the other hand, the associated driven map with $`L_t=C`$ has a unique asymptotic period-two orbit $`\overline{s}_t=[\alpha ,\beta ]`$ on a range of $`C`$, where $`\alpha `$ and $`\beta `$ are functions of $`C`$. The associated coupling function $`\mathrm{\Theta }_H`$ also simplifies in such case. For the reduced, two-cluster couplings in Eqs (6)-(9), the corresponding associated coupling functions become

$`\mathrm{\Theta }_x(\alpha ,\beta )`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\alpha +\beta )`$ (16)
$`\mathrm{\Theta }_f(\alpha ,\beta )`$ $`=`$ $`{\displaystyle \frac{1}{2}}[f(\alpha )+f(\beta )]`$ (17)
$`\mathrm{\Theta }_{\mathrm{\Delta }x}(\alpha ,\beta )`$ $`=`$ $`{\displaystyle \frac{1}{2}}|\alpha \beta |`$ (18)
$`\mathrm{\Theta }_{\overline{x}}(\alpha ,\beta )`$ $`=`$ $`|\alpha |^{1/2}|\beta |^{1/2}.`$ (19)
Then Eqs. (14)-(15) with $`L_1=L_2=C`$ reduce to the single equation

$$\mathrm{\Theta }_H(\alpha ,\beta )=C,$$
(20)
which can be seen as an equation for $`C`$, for given values of the parameters $`ϵ`$ and $`r`$. The solution $`C=C_{}`$ of Eq. (20) provides a complete description of the GCM state since then $`a=\alpha (C_{}),b=\beta (C_{})`$ and $`H=C_{}`$. Figure 4 shows the bifurcation diagram of $`s_t`$, Eq. (12), as a function of the constant drive $`L_t=C`$ up to period two, with fixed parameters $`r=2`$ and $`ϵ=0.4`$. The fixed point region in this diagram corresponds to one stationary cluster (i.e., a synchronized collective state) in the GCM, Eq. (1), with constant $`H`$. The period-two window corresponds to the values $`\alpha `$ and $`\beta `$ adopted by the driven map on this range of $`C`$. Once $`\alpha (C)`$ and $`\beta (C)`$ are known from the bifurcation diagram, the function $`\mathrm{\Theta }_H(C)`$ associated to any global coupling function in a GCM can be readily obtained, assuming that the GCM is in a state of two equal size clusters $`(p=0.5)`$, evolving out of phase with respect to each other. In Figure 4, the $`\mathrm{\Theta }_H`$ functions associated to the four global couplings in Eqs.(2)-(5) with fixed $`ϵ`$ are shown as function of $`C`$. As stated above, the solutions $`C_{}`$ to Eq. (20) correspond to states in GCM systems with a coupling function reaching a stationary value $`H=C_{}`$. Thus, the intersections of the $`\mathrm{\Theta }_H`$ curves with the diagonal in Fig. 4 give all the possible states of GCM systems that maintain a constant $`H`$, either with one stationary cluster (if the intersection occurs on the fixed point window of the bifurcation diagram of $`s_t`$) or with two clusters in period two (if the intersection occurs on the period two window of the diagram). Note that both the backward and the forward mean field couplings have the same two-cluster, period-two solution $`x=f=C_{}=0.416`$, but these couplings have different functional dependence on the constant drive $`L_t=C`$. The coincidence of the couplings $`x`$ and $`f`$ for $`p=0.5`$ was already seen in Fig. 2(b). Similarly, the geometric mean coupling $`H=\overline{x}`$ gives only one two-cluster, period-two solution at $`\overline{x}=C_{}=0.305`$. In contrast, the dispersion global coupling, $`H=\mathrm{\Delta }x`$, has three intersections with the diagonal: one corresponds to the synchronized stationary state in the associated GCM, with $`\mathrm{\Delta }x=C_{}=0`$, and the other two correspond to different clustered states of the GCM, each consisting of two equal size clusters in period two, with $`\mathrm{\Delta }x=C_{}=0.083`$ and $`\mathrm{\Delta }x=C_{}=0.25`$, respectively. All the states predicted by the intersections of the different $`\mathrm{\Theta }_H`$ with the diagonal in Fig. 4, except one, are readily found in simulations on the corresponding GCM systems for appropriated initial conditions in each case. Actually, for a GCM with the coupling $`H=\mathrm{\Delta }x`$, the predicted two-cluster, period-two state with $`\mathrm{\Delta }x=C_{}=0.083`$ is unstable: it is never achieved in simulations on the GCM system, even when the initial conditions are chosen very close to that state. What is observed, instead, is the evolution of the GCM system towards either the stationary one cluster state with $`\mathrm{\Delta }x=0`$ or the two-cluster, period-two state with $`\mathrm{\Delta }x=0.25`$. Thus, in addition to being predicted by the solutions of the equation $`\mathrm{\Theta }_H=C`$, the observed clustered states of a GCM displaying constant coupling must be stable, which implies some stability condition on the solutions. It can be shown (see Appendix) that for coupling functions satisfying $`_{i=1}^NH/x_i=0`$, the condition $`d\mathrm{\Theta }_H/dC>1`$ at the intersection with the diagonal implies that the corresponding solution is unstable. This is the case of the global coupling $`H=\mathrm{\Delta }x`$. Note that the solution at $`C_{}=0.083`$ is the only one for which $`\frac{d\mathrm{\Theta }}{dC}|_C_{}>1`$ in Fig. 4 and therefore it is unstable, independently of the cluster fraction $`p`$. For the coupling $`H=f`$, a stability analysis of states consisting of two or three clusters in period three has been performed by Shimada and Kikuchi . However, the simple criterium for instability $`\frac{d\mathrm{\Theta }}{dC}|_C_{}>1`$ does not apply for $`H=f`$.

Constant coupling functions may also occur in GCM systems with different cluster sizes; that is the case of a GCM possessing dispersion global coupling $`H=\mathrm{\Delta }x`$ and displaying two clusters with any partition $`\{p,1p\}`$, as seen in Fig. 1(b). Figure 5 shows the associated function $`\mathrm{\Theta }_{\mathrm{\Delta }x}`$ with fixed $`ϵ`$ as a function of the constant drive $`C`$ for several values of the fraction $`p`$. There exist a critical fraction $`p_c=0.25`$ bellow which only one solution corresponding to one stationary cluster, i.e., synchronization, can appear in the GCM system. Above this critical fraction, two states, each consisting of two clusters in period two, are additionally predicted by the solutions $`C=C_{}`$ of $`\mathrm{\Theta }_{\mathrm{\Delta }x}=C`$. These solutions emerge as a pair: one solution is always unstable since $`d\mathrm{\Theta }_H/dC|_C_{}>1`$ there, and the other is the one two-cluster, period-two collective state with constant $`\mathrm{\Delta }x`$ that is actually observed in simulations on the corresponding GCM system. For the fraction $`p_c=0.25`$, there is a two-cluster, period-two solution $`\mathrm{\Theta }_{\mathrm{\Delta }x}=C_{}=0.125`$ that is marginally stable.

## V Conclusions

Most studies on GCM lattices and other globally coupled systems have assumed mean field coupling. However, other forms of global coupling may be relevant in some situations. We have analyzed, in a general framework, the clustering behavior in GCM systems subjected to permutable global coupling functions by considering the dynamics of the coupling functions.

We have shown that different GCM systems can be represented by the orbits of their coupling functions on a common space. For simplicity, only collective states in GCM systems consisting of two clusters in period two were considered. We have shown that the functional form of the global coupling in a GCM system determines the periodicity of its motion and the possible distributions of elements among the clusters. The existence of a well defined interval of possible partitions among two clusters, out of which no clusters emerge in the system, has been observed experimentally . In experimental or natural situations where clustering occurs, the specific functional form of the coupling is in general unknown. The present study may be useful to obtain insight into the acting global coupling function in practical situations.

We have employed a previously introduced analogy between a GCM system and a single externally driven map  in order to give a unified interpretation of the observed clustering behavior of the GCM systems considered in this article. A periodically driven map with local periodic windows can display multiple asymptotic periodic responses which are similar to cluster orbits in a GCM system with permutable $`H`$. This analogy implies that dynamical clustering can occur in any GCM system with a permutable coupling function and periodic windows in the local dynamics. The presence of windows of stable periodic orbits in the local map is essential for the emergence of clusters. In fact, no clustering is observed in a GCM system if the local maps do not have periodic windows ; what is observed instead is synchronization or nontrivial collective behavior, i.e., an ordered temporal evolution of statistical quantities coexisting with local chaos.

The associated coupling function derived from the driven map analogy is particularly simple to use in the prediction of clustered states in GCM systems with two equal size clusters and exhibiting constant global coupling $`H=C`$. The associated coupling function $`\mathrm{\Theta }_H`$ can be directly constructed from the bifurcation diagram of the steadily driven map. The cluster states are obtained from the solutions of Eq. (20) and can be represented graphically in a simple way. Although Eq. (20) has been used for the case of two clusters in period two, it can also be applied to find GCM states consisting of $`K`$ equal size clusters in period $`K`$. In addition, the associated coupling function $`\mathrm{\Theta }_H`$ carries information about the stability of the predicted two-cluster states. In particular, for the family of coupling function satisfying property (29), the stability condition of clustered states in a GCM with constant $`H=C_{}`$ is directly given by the slope $`d\mathrm{\Theta }_H/dC|_C_{}`$. The example of a GCM system with dispersion coupling function $`H=\mathrm{\Delta }x`$ reveals that a constant coupling can also be maintained by clusters of different sizes. Our method based on Eq. (20) also predicts successfully the cluster states in these situations.

The driven map analogy suggests that the emergence of clusters should be a common phenomenon which can be expected in various dynamical systems formed by globally interacting elements possessing stable periodic orbits on some parameter range of their individual dynamics. The examples presented here show that progress in the understanding of the collective behavior of globally coupled systems can be achieved by investigating their relation to a driven oscillator.

## Acknowledgment

This work has been supported by Consejo de Desarrollo Científico, Humanístico y Tecnológico of the Universidad de Los Andes, Mérida, Venezuela.

## On the relation between Stability and the Associated Coupled Function

Consider a general GCM system with any global permutable coupling function. Suppose that the system reaches a state consisting of two clusters. Then the dynamics of the system reduces to two coupled maps Eqs. (10)-(11), i.e.,

$$\begin{array}{ccc}\chi _{t+1}(1)=\hfill & (1ϵ)f(\chi _t(1))+ϵH(\chi _t(1),\chi _t(2),p)=\hfill & F(\chi _t(1),\chi _t(2)),\hfill \\ \chi _{t+1}(2)=\hfill & (1ϵ)f(\chi _t(2))+ϵH(\chi _t(1),\chi _t(2),p)=\hfill & G(\chi _t(1),\chi _t(2)).\hfill \end{array}$$
(21)
If the two clusters are in period-two orbits, the stability of this collective state in the GCM is given by the eigenvalues of the product of Jacobian matrices

$$𝐉=𝐉_1𝐉_2=\underset{i=1}{\overset{2}{}}\left(\begin{array}{cc}\frac{F}{\chi _i(1)}& \frac{F}{\chi _i(2)}\\ & \\ \frac{G}{\chi _i(1)}& \frac{G}{\chi _i(2)}\end{array}\right).$$
(22)
If the two clusters move out of phase, the asymptotic state of the GCM can be described by the two orbits $`\chi _t(1)[\chi _1(1),\chi _2(1)]=[a,b]`$ and $`\chi _t(2)[\chi _1(2),\chi _2(2)]=[b,a]`$, which satisfy

$$\begin{array}{cc}b=\hfill & (1ϵ)f(a)+ϵH(a,b,p),\hfill \\ a=\hfill & (1ϵ)f(b)+ϵH(a,b,p).\hfill \end{array}$$
(23)
For the local dynamics $`f(x)=1rx^2`$, one gets

$$𝐉_1=\left(\begin{array}{cc}2(1ϵ)ra+ϵH_a& ϵH_b\\ & \\ ϵH_a& 2(1ϵ)rb+ϵH_b\end{array}\right),$$
(24)
$$𝐉_2=\left(\begin{array}{cc}2(1ϵ)rb+ϵH_b& ϵH_a\\ & \\ ϵH_b& 2(1ϵ)ra+ϵH_a\end{array}\right),$$
(25)
where

$$H_a\frac{H(\chi _t(1),\chi _t(2))}{\chi _t(1)}|_{\chi _t(1)=a,\chi _t(2)=b}=\frac{H(\chi _t(1),\chi _t(2))}{\chi _t(2)}|_{\chi _t(1)=b,\chi _t(2)=a},$$
(26)
and

$$H_b\frac{H(\chi _t(1),\chi _t(2))}{\chi _t(2)}|_{\chi _t(1)=a,\chi _t(2)=b}=\frac{H(\chi _t(1),\chi _t(2))}{\chi _t(1)}|_{\chi _t(1)=b,\chi _t(2)=a}.$$
(27)
Consider now the dispersion coupling function, Eq. (4). This coupling belongs to the family of functions of $`N`$ variables

$$H(x_1,x_2,\mathrm{},x_N)=\underset{i=1}{\overset{N}{}}h_e(x_ix),$$
(28)
where $`h_e`$ is any even function of its argument. It can be straightforwardly shown that this family of functions possesses the property

$$\underset{i=1}{\overset{N}{}}\frac{H}{x_i}=0.$$
(29)
Therefore, in a two cluster state, any $`H`$ in this family of global coupling functions satisfies

$$\frac{H(\chi _t(1),\chi _t(2))}{\chi _t(1)}+\frac{H(\chi _t(1),\chi _t(2))}{\chi _t(2)}=0.$$
(30)
If the two clusters evolve out of phase with respect to each other, and additionally the GCM has a coupling $`H`$ with property (29), then the two eigenvalues of the matrix $`𝐉`$ in Eq. (22) become identical and their value is

$$\lambda =2rϵ(ϵ1)(aH_b+bH_a)+4r^2ab(ϵ1)^2.$$
(31)
The stability criterion of this state is given by the modulus of the eigenvalue $`\lambda `$; that is, $`|\lambda |>1`$ ($`|\lambda |<1`$) implies that the state is unstable (stable). The values $`a`$ and $`b`$ are, respectively, the values of $`\alpha `$ and $`\beta `$ at the intersection of the function $`\mathrm{\Theta }_{\mathrm{\Delta }x}`$ with the diagonal in Fig. 4.

Let us analyze the relationship between the eigenvalue $`\lambda `$ and the derivative $`d\mathrm{\Theta }_{\mathrm{\Delta }x}/dC`$ at the intersection points with the diagonal in Fig. 4 or Fig. 5. In general,

$$\frac{d\mathrm{\Theta }}{dC}=\frac{\mathrm{\Theta }}{\alpha }\frac{\alpha }{C}+\frac{\mathrm{\Theta }}{\beta }\frac{\beta }{C},$$
(32)
where

$$\alpha =\frac{1+R}{2r(ϵ1)},$$
(33)
$$\beta =\frac{1R}{2r(ϵ1)},$$
(34)
and

$$R=\left(34rϵ^2C+4rϵC+4r8rϵ+4rϵ^2\right)^{1/2}.$$
(35)
Since

$$\frac{\alpha }{C}=\frac{\beta }{C}=\frac{1}{2r(ϵ1)}\frac{R}{C}=\frac{ϵ}{R},$$
(36)
then

$$\frac{d\mathrm{\Theta }}{dC}=\frac{ϵ}{R}\left(\frac{\mathrm{\Theta }}{\alpha }\frac{\mathrm{\Theta }}{\beta }\right).$$
(37)
Since $`\mathrm{\Theta }_{\mathrm{\Delta }x}`$ has the same functional form as $`H=\mathrm{\Delta }x`$, then $`\mathrm{\Theta }_{\mathrm{\Delta }x}`$ also satisfies property (30), that is

$$\frac{\mathrm{\Theta }_{\mathrm{\Delta }x}}{\alpha }=\frac{\mathrm{\Theta }_{\mathrm{\Delta }x}}{\beta },$$
(38)
and therefore

$$\frac{d\mathrm{\Theta }_{\mathrm{\Delta }x}}{dC}=\frac{2ϵ}{R}\frac{\mathrm{\Theta }_{\mathrm{\Delta }x}}{\alpha }.$$
(39)
Let $`C=C_{}`$ be a value of $`C`$ corresponding to the intersection of $`\mathrm{\Theta }_{\mathrm{\Delta }x}`$ with the diagonal in Fig. 4. Then $`\alpha (C_{})=a`$ and $`\beta (C_{})=b`$; and $`\mathrm{\Theta }_{\mathrm{\Delta }x}(C_{})=H(a,b)`$, which gives

$$\frac{d\mathrm{\Theta }_{\mathrm{\Delta }x}}{dC}|_C_{}=\frac{2ϵ}{R}H_a.$$
(40)
Using the fact that $`H_a=H_b`$ from Eq. (30), the eigenvalue $`\lambda `$ becomes

$$\lambda =2rϵ(ϵ1)H_a(ba)+4r^2ab(ϵ1).$$
(41)
Eqs.(33) and (34) with $`C=C_{}`$ give the values $`a`$ and $`b`$, respectively. Then, substitution of these values and $`H_a`$ from Eq.(40) in Eq. (41) yields

$$\lambda =R^2\left(\frac{d\mathrm{\Theta }_{\mathrm{\Delta }x}}{dC}|_C_{}1\right)+1.$$
(42)
Therefore, the condition $`d\mathrm{\Theta }/dC|_C_{}>1`$ implies that $`|\lambda |>1`$, and thus the two-cluster, period-two solution with $`C_{}=0.07`$ given by the intersection of $`\mathrm{\Theta }_{\mathrm{\Delta }x}`$ with the diagonal in Fig. 4 is unstable. Similarly, the solutions $`C_{}`$ of $`\mathrm{\Theta }_{\mathrm{\Delta }x}=C`$ for the different curves in Fig. 5 for which $`d\mathrm{\Theta }/dC|_C_{}>1`$, are unstable.

Note that the above stability result for $`H=\mathrm{\Delta }x`$ is also valid for any global coupling function satisfying property (29).