text:
# Comparison of Molecular Dynamics with Hybrid Continuum-Molecular Dynamics for a Single Tethered Polymer in a Solvent

## I Introduction

Molecular dynamics (MD) simulations have long been used to model complex fluids both in and out of equilibrium. As computers get more powerful there has been an increasing desire for more chemically accurate models of these fluids. This means that simulations are becoming larger and more accurate, but also that much simulation time is being devoted to model in detail parts of the computational system of little direct scientific interest. Hybrid methods combine regions of relatively high degree of chemical accuracy in a specific region of interest, and a more coarse-grained model where the dynamics can be solved in a less computationally intensive way, further away from the specific region of interest. We focus, in particular, on hybrid methods that combine Lennard-Jones type specificity with larger scale continuum methods. Such hybrid methods have been applied in a number of fields, including Lennard-Jones fluids O’Connell and Thompson (1995), biophysics Ayton et al. (2001) and MD/CFD coupling Li et al. (1998, 1999). This type of simulation technique is particularly useful in studying interface problems, where the region of interest is a localized part of the entire system.

Typical hybrid methods consist of three regions: a traditional region where dynamics are simulated using well-established techniques such as molecular dynamics Allen and Tildesley (1987), a continuum region where CFD or elasticity differential equations are solved using classical techniques, and an overlap region where the necessary transport information of the MD and continuum regions are exchanged. The primary motivation for using a hybrid scheme is to reduce computer time devoted to simulating bulk regions of little direct interest. As such, a hybrid scheme is ideally suited to studying interfacial systems.

In this paper we apply a hybrid technique to a single polymer tethered to a wall with explicit solvent. The complex dynamics arising from this system have attracted a degree of interest from experimentalists, who used fluorescence microscopy and videomicroscopy to investigate the dynamic properties of individual DNA chains in a shear flow, either tethered to a wall Doyle et al. (2000) or free LeDuc et al. (1999). These experiments reveal that the structural quantities, such as the mean elongation of the polymer, are very sensitive to flow environment and that the dynamical properties depend strongly on the initial conformation. Moreover, care needs to be taken to control the finite size effects, such as those due to long-ranged hydrodynamic interaction between the polymer and the walls LeDuc et al. (1999). This large “sensitivity” of the tethered polymer dynamics is in fact a valuable test for the hybrid model. First, the hybrid model reduces the size of the MD simulation box while avoiding finite size effects and, second, the coupling has to be able to perfectly reproduce flows at very small shear rates. As shown in recent work Delgado-Buscalioni et al. (2003), this second task is non-trivial because the signal-to-noise ratio of the stress that one needs to communicate from the particle to the continuum system is very small.

The problem of tethered polymers under flow has a geometry which is ideal for a hybrid scheme. The scientific interest lies around the polymer although, in a standard MD simulation, the solvent particles within the bulk flow require most of the computational time. Single polymers in a bath of explicit molecular solvent have been the focus of a great deal of attention in the last decade Ahlrichs and Dünweg (1999); Dünweg and Kremer (1993); Aust et al. (1999); Pierleoni and Ryckaert (1995); B. Dünweg and Kremer (1998); Jendrejack et al. (2002). Many of these studies are devoted to examining a free chain in solution in order to make comparisons with theoretical predictions, explore the dynamics regime beyond the short-time dynamics or extract scaling laws as a function of polymer length or shear flow. In these studies, the solvent is explicity simulated. For example, the study by Dünweg and Kremer Dünweg and Kremer (1993) uses a polymer of length $`L=60`$ beads in a bath of $`7940`$ Lennard-Jones spheres. Aust, Kröger and Hess Aust et al. (1999) simulate polymers of length $`L=10`$ to $`60`$ in systems where the total number of particles including solvent ranges from $`1000`$ to $`5832`$. It is clear in these cases that most of the computational effort is devoted to solving the equations of motion of the solvent particles when the real scientific interest lies in the polymer behavior.

The single polymer we study is tethered to a wall, and a variety of shear rates is imposed as a model interfacial problem to compare classical MD techniques to the hybrid simulation. In classical MD, we sandwich the polymer and solution between two explicit walls, and impose periodic boundary conditions in the remaining two directions. The polymer is tethered to the bottom wall, and shear is created by moving the top wall at constant velocity in a direction parallel to the wall. In the hybrid case, we model one wall, the polymer and some of the solvent explicitly using MD, and impose shear by a boundary condition in the CFD regime of the calculation. The shear is translated to the MD regime via energy and momentum flux transfers in the overlap region. We compare various conformational properties of the polymer for the two techniques.

Our paper is organized as follows. In the following section we briefly outline both the classical MD simulation and the hybrid simulation techniques. In Sec. III we compare the conformation of the polymer as calculated by each simulation method. The computational costs and benefits of the hybrid scheme are compared to classical MD. We also compare our results to experimental data of tethered DNA under shear flow. We conclude with a discussion in Sec. IV.

## II  Method

We describe in this section both the molecular dynamics and hybrid dynamics models used in our simulations. The MD part of the hybrid scheme was the same as the classical MD used in the pure molecular dynamics simulations.

### II.0.1 Molecular Dynamics

The polymer model and simulation techniques are similar to those used in previous work Thompson and Robbins (1990); Barsky and Robbins (2001). The polymer potential is based on the bead-spring model developed by Kremer and Grest Kremer and Grest (1990). Linear polymers containing $`N=60`$ beads each are created by linking nearest neighbors on a chain with the potential

$$U_{nn}(r_{ij})=\{\begin{array}{cc}\frac{1}{2}kR_0^2\mathrm{ln}\left[1\left(r_{ij}/R_0\right)^2\right]\hfill & r_{ij}<R_0\hfill \\ \mathrm{}\hfill & r_{ij}R_0,\hfill \end{array}$$
(1)
where $`r_{ij}`$ is the distance between beads $`i`$ and $`j`$, $`R_0=1.5\sigma `$, $`k=30ϵ/\sigma ^2`$, and $`\sigma `$ and $`ϵ`$ set the length and energy scales, respectively. The monomers in the solvent and in the polymer interact through a truncated Lennard-Jones (LJ) potential

$$U_{LJ}(r_{ij})=\{\begin{array}{cc}4ϵ\left[\left(\sigma /r_{ij}\right)^{12}\left(\sigma /r_{ij}\right)^6\right]\hfill & r_{ij}<r_c\hfill \\ 0\hfill & r_{ij}r_c.\hfill \end{array}$$
(2)
The cutoff is set at $`r_c=2^{1/6}\sigma `$ for all fluid particles to produce a purely repulsive interaction between beads.

The bounds of the simulation cell are periodic in the $`x`$ and $`y`$ directions, with periods $`L_x38.5\sigma `$ and $`L_y33.4\sigma `$, respectively. In the $`z`$ direction the cell is bounded by top and bottom walls. Each wall contains two layers of $`1600`$ spheres strongly tethered to the sites of a $`(1,1,1)`$ plane of an fcc lattice by harmonic springs of stiffness $`\kappa =1320ϵ\sigma ^2`$. The wall atoms do not interact with each other, and the wall-fluid interaction is LJ with an increased cutoff of $`r_c=1.25\sigma `$ and increased energy scale of $`ϵ_{wf}=\sqrt{1.7}ϵ`$. The increased cutoff and energy ensure sufficient adhesion of the fluid to the wall so that the slip at the wall is minimized for the shear rates considered here. The polymer is anchored to the wall by enforcing the tethering potential, Eq. (1) between the end of the polymer and one wall atom.

The walls are $`48\sigma `$ apart for the pure MD simulation; in the hybrid simulation the molecular dynamics region persists for $`19\sigma `$. There are sufficient solvent monomers to yield a mean fluid density of approximately $`\rho =0.8\sigma ^3`$ in the center of the simulation cell, although density oscillations are induced within a few $`\sigma `$ of the walls Khare et al. (1996).

The equations of motion are integrated using a velocity Verlet algorithm Allen and Tildesley (1987), with a time step $`\delta t=0.0075\tau `$, where $`\tau =\sigma \sqrt{m/ϵ}`$ is the basic unit of time, and $`m`$ is the mass of a monomer. A constant temperature of $`k_BT=1.0ϵ`$ is maintained with a Langevin thermostat Kremer and Grest (1990). To ensure that this thermostat does not bias the shear profile, the Gaussian white noise and damping terms are only added to the equations of motion for the velocity components normal to the mean flow, that is the $`y`$ and $`z`$ directions Thompson and Robbins (1990).

The shear flow in our pure MD simulation is induced by moving the atomistic top wall at a constant speed $`v_x`$ in the $`x`$ direction. In the hybrid simulation, a shear boundary condition is used for the continuum regime, and this resulted in shear in the MD regime by the exchange of momentum in the overlap domain. The starting configuration was that of a single polymer tethered to the wall, in an equilibriated solvent. We repeat the simulation for two different starting configurations, i.e. each configuration has a polymer tethered to a different wall atom, and the initial conformation of the polymer is different. The initial polymer configurations were either taken from previous simulations on melts Barsky and Robbins (2001), or generated from a random walk. Although over long periods of time we expect that different starting configurations will give the same configurational averages, previous work Dünweg and Kremer (1993) has shown that hundreds of different initial configurations are required to arrive at reasonable ensemble averages. In view of this, we used two initial configurations for the pure MD simulations; although this falls short of the number of initial configurations required to achieve ensemble averages, it does give us a window over which to compare the hybrid simulation.

The local shear rate, $`\dot{\gamma }`$, of the fluid is calculated by computing the local change in the $`x`$ component of velocity, $`v_x`$, as a function of $`z`$, i.e. $`\dot{\gamma }=v_x/z`$. The upper wall velocity was chosen so that the shear rate $`\dot{\gamma }`$ assumed the values $`0.0,0.0005,0.001,0.002,0.005,0.01\tau ^1`$. Simulations at higher shear rates were created from lower shear rates by increasing the wall velocity or boundary condition and allowing the system to achieve steady state. The simulations were done for at least one million time steps, and the runs of $`\dot{\gamma }=0.001\tau ^1`$ and $`0.005\tau ^1`$ were simulated for at least ten million time steps, corresponding to a total run time of $`75,000\tau `$. In the analysis in the following section the first $`250,000`$ time steps of data for each given shear rate were ignored, to allow the system to reach steady state; this length of time was determined to be the longest time necessary for the system to reach steady state once a new shear rate was imposed.

### II.0.2 Hybrid Dynamics

Our hybrid dynamics code Delgado-Buscalioni and Coveney (2003a) consists of three domains: the particle domain (P) which was by the same molecular dynamics method described above, the continuum domain (C) treated by standard continuum fluid dynamics and, an overlap region where information from the other two domains is exchanged. The hybrid scheme is a protocol to exchange fluxes of conserved quantities, specifically mass, momentum and energy between both classically treated regimes. To implement the two-way flux exchange, the overlap region consists of two different subdomains: the P$``$C and the C$``$P cells. Within the C$``$P region, the fluxes from the continuum domain are imposed on the particle domain, whereas within the P$``$C cell the microscopic fluxes are coarse-grained in time and space Delgado-Buscalioni and Coveney (2003a) to supply boundary conditions for the continuum domain.

The spatial decomposition used for the present set-up is shown in Fig. 1. The molecular dynamics domain ranges from the atomistic wall at $`z0`$ and extends to $`z=l_{CP}=19\sigma `$. The continuum fluid dynamics domain comprises $`l_{PC}zL_z`$, where $`l_{PC}14.5\sigma `$ is the $`z`$coordinate of the P$``$C interface and $`L_z=50\sigma `$ is the extent of the whole simulation domain. The center of the P$``$C cell is located at $`z=13.4\sigma `$; it has a volume $`V_{PC}=\mathrm{\Delta }z_{PC}A`$ where $`A=L_x\times L_y`$ and $`\mathrm{\Delta }z_{PC}2.2\sigma `$ is the extension along the $`z`$ direction. The C$``$P cell is placed at a distance $`2.2\sigma `$ from the end of the P$``$C cell and covers a region of $`\mathrm{\Delta }z_{CP}2.2\sigma `$, from $`z16.8\sigma `$ to $`z=l_{CP}19\sigma `$.

In what follows we outline the coupling protocol and provide the numerical details used in the present implementation. The C$``$P coupling represents the most complicated part of the hybrid scheme; a more detailed explanation of the method in the frame of the general case of unsteady flows with mass, momentum and energy exchanges can be found in reference Delgado-Buscalioni and Coveney (2003b). The steady flow considered here only carries momentum along the $`x`$ direction. Although the mean flux of mass and energy across the C and P interfaces is zero, fluctuations in the particle system produce perturbative mass currents along the $`z`$ direction which need to be taken into account. This part of the C$``$P scheme is presented in Appendix A.

We now focus on how the momentum flux is exchanged between the C and P domains, starting with a discussion of the C$``$P coupling. For the pure Couette shear flow considered here, the momentum flux due to the C flow across any $`z=`$constant surface is given by

$$𝚷=P𝐤\eta \dot{\gamma }𝐢,$$
(3)
where $`P=P(\rho ,T)`$ is the hydrostatic pressure, $`\eta `$ is the dynamic viscosity and $`\dot{\gamma }v_x/z`$ is the shear rate. The value of the dynamic viscosity was measured in a previous pure MD simulation via the standard non-equilibrium procedure Aust et al. (1999); Barsky and Robbins (2001); Evans and Morris (1990); for $`\rho =0.8`$ and $`T=1.0`$ we obtained $`\eta =1.75\pm 0.04`$. The stress induced by the C-flow in the P domain is given by the local momentum flux at the C$``$P interface, $`𝚷_{\mathrm{𝐂𝐏}}`$. In order to introduce this stress into the molecular dynamics domain we add an external force $`𝐅_{\mathrm{𝐞𝐱𝐭}}=𝚷_{\mathrm{𝐂𝐏}}A`$ to those molecules within the C$``$P cell. At any instant of time, $`t`$, this force is equally distributed among the $`N_{CP}(t)`$ particles inside the C$``$P cell, so the external force per particle is $`𝐅_{\mathrm{𝐞𝐱𝐭}}/N_{PC}=𝚷_{\mathrm{𝐂𝐏}}A/N_{CP}`$. Note that this external force has a component normal to the C$``$P interface, which provides the hydrostatic pressure, and a tangential component providing the shear stress. The molecules are free to enter or leave the C$``$P region, so the number of molecules within this region, $`N_{CP}(t)`$, and the value of the overall external force fluctuate in time. The average “pressure force” per particle is $`P(\rho ,T)/(A\overline{N}_{CP})`$, where $`\overline{N}_{CP}2000`$ is the mean number of particles within the C$``$P cell (see Appendix A), $`A=L_x\times L_y=1286\sigma ^2`$ and the pressure $`P(\rho ,T)`$ is given from the equation of state provided by Hess et al. Hess et al. (1998): $`P(0.8,1)6.5ϵ/\sigma ^3`$. Such a force prevents the escape of particles and, although it induces some ripples on the density profile over the C$``$P cell, it maintains the correct value of the density along the inner part of the MD domain, as seen in Fig. 4b and discussed further in the Appendix A.

The shear force is distributed over the particles in the same way as described above for the pressure force. In this case, the flux of $`x`$-momentum to be injected in the particle system is $`\eta \dot{\gamma }_{CP}`$, where $`\dot{\gamma }_{CP}`$ is the local shear rate of the C-flow measured at the C$``$P interface.

We next discuss the P$``$C coupling. The continuum domain is a coarse-grained description of the fluid, therefore any information transferred from the molecular to the continuum system needs to be averaged in space and time. These averages need to be local in the continuum space and time coordinates. To that end, the particle quantities are averaged within the P$``$C cell and over a time interval $`\mathrm{\Delta }t_{av}`$. It is important to stress that within the P$``$C cell each particle’s dynamics are not directly modified by any external artifact, in other words the motion of each particle is uniquely determined by the usual molecular dynamics scheme. To ensure consistency within the hybrid scheme, $`\mathrm{\Delta }t_{av}`$ and the volume of the P$``$C cell are restricted 111As explained in references Delgado-Buscalioni and Coveney (2003b, c), these restrictions come from the C flow (stability of the scheme and spatio-temporal resolution) and from the P system (local thermodynamic equilibrium, signal-to-noise ratio).. For the steady flow employed in this study, the most compelling condition is to guarantee that the signal-to-noise ratio of the momentum flux is larger than one and for that reason $`\mathrm{\Delta }t_{av}`$ needs to increase as $`\dot{\gamma }`$ decreases Delgado-Buscalioni et al. (2003). We used $`\mathrm{\Delta }t_{av}=100\tau `$ for $`\dot{\gamma }0.001`$ and reduced it gradually to $`10\tau `$ for the fastest flows considered.

To solve the equations of motion in the continuum domain we used the finite volume formulation Patankar (1980) because it matches by construction the fluxes across cells. Since the solvent flow is isothermal, incompressible and there is a uniform pressure, the mean $`x`$-velocity is governed by $`v_x/t=\nu ^2v_x/z^2`$, where $`\nu =\eta /\rho `$ is the kinematic viscosity and $`v_x`$ is the velocity in the $`x`$ direction. At the top of the simulation cell we impose a smooth wall in the CFD sense. This wall moves at a constant velocity $`v_x(L_z,t)=u_{wall}`$ which creates the shear flow in the simulation. The protocol for the P$``$C coupling establishes the boundary condition for the continuum domain at the P$``$C interface, $`z=l_{PC}`$. The coarse-grained microscopic flux of $`x`$-momentum across the P$``$C interface, whose expression is given in Ref. Delgado-Buscalioni and Coveney (2003b), is set equal to the corresponding value for the C flow at the $`z=l_{PC}`$ boundary, $`\eta \dot{\gamma }_{PC}`$, where $`\dot{\gamma }_{PC}(v_x(l_{PC}+\mathrm{\Delta }z)v_x(l_{PC}))/\mathrm{\Delta }z`$. This condition gives the desired velocity to be imposed at the boundary $`v_x(l_{PC})`$. The continuity of velocity is ensured by adding a relaxing term in the flux equation which drives the C-velocity at the interface towards the corresponding averaged P-velocity (see references Delgado-Buscalioni and Coveney (2003c); Delgado-Buscalioni et al. (2003) for details).

## III  Results and Discussion

In this section we compare the conformational behavior of the polymers from the MD and hybrid simulations; we use two independent MD simulations for comparison. We study two MD systems because it is well known Jendrejack et al. (2002) that two simulations, or experiments, on a tethered polymer may exhibit considerable variation in conformational behavior, even at rather high shear rates. We conclude this section with a discussion of the computational costs and benefits of the hybrid and classical MD techniques.

In Figure 2 we show the mean-square end-to-end distance $`R^2`$ of the polymer in each of the $`x,y`$ and $`z`$ directions. This is a standard measure of polymer conformation Doi and Edwards (1986), and its values are related to the values of the radius of gyration. At low shear rates the polymer conformation calculated from the hybrid simulation is well within the measured conformations of the MD simulations. At the highest shear rates, Fig. 2(a) shows the conformation of the hybrid polymer to be about $`10\%`$ larger than the polymer in the MD simulation; this difference is within the standard deviation of $`R_{xx}^2`$, which is about $`15\%`$. Figs.2(b) and (c) show the $`y`$ and $`z`$ components of $`R^2`$ as a function of shear rate $`\dot{\gamma }`$. In both cases the two MD simulations serve as good indicators of the variability of the conformational behavior of a single polymer; the conformation of the polymer in the hybrid simulation is well within the variations found in the two MD simulations, at all shear rates.

Figure 3 shows the probability of the maximum extension from as a function of distance along $`x,y`$ and $`z`$ directions, for a shear rate of $`\dot{\gamma }=0.001\tau ^1`$. It is clear that the variation of the distributions obtained with the hybrid simulation is well within the distribution of the two MD simulations. This indicates that not only is the average conformation comparable between the two simulation techniques, but that the probability distributions also compare favorably.

In Figure 4 we show the density, of both the solution monomers and polymer, $`\rho =N/V`$ as a function of distance from the wall for both the MD and hybrid simulations. The density is calculated in slices of approximately $`0.01\sigma `$ perpendicular to the wall. The regular spacing of the wall monomers, as two monolayers of a $`(1,1,1)`$ face of an fcc crystal, induces an ordering in the fluid; this ordering is well established Thompson and Robbins (1990) and persists for approximately $`5\sigma `$. At the wall the monomer density variations are identical for both the MD and hybrid simulations as seen in Fig 4(a). In Fig. 4(b) we see that the density of both simulations remains the same until the monomers in the hybrid system feel the effects of the constant pressure condition imposed on the overlap region. The constant pressure is implemented as a simple normal force per particle on all monomers in the C$``$P regime, as discussed in the previous section. This force induces a local ordering in the monomers, which in turn creates density fluctuations. It is noteworthy, however, that these density oscillations are much lower than at the atomistic wall, shown for comparison. More recent work on the hybrid scheme has established that we can reduce these density fluctations even further, as discussed in the Appendix A. In the MD simulation, the upper wall is identical to the lower one, and thus the density fluctuations near the former are the same as those at the latter.

Fig. 5 shows the probability of finding any polymer bead in a plane, where the plane slices are $`0.2\sigma `$ in thickness. The two-dimensional probabilities were calculated in an analogous way to the one-dimensional probabilities discussed above. The shear rate shown is $`\dot{\gamma }=0.001`$. Inspection of the two-dimensional bead distributions indicates that below a distance of $`5\sigma `$ to the wall, the beads tend to be ordered in layers parallel to the wall plane. This result is not only a consequence of the polymer-wall interaction but also an effect of the interaction with the solvent. Near the wall the solvent is ordered in layers, as in Fig. 4(a), and the polymer minimizes the monomer-solvent potential energy by adapting its distribution to match the locations of the solvent layers. The order induced by the wall in the polymer structure can be noticed even in the isovalues of the probability distribution along the wall plane $`xy`$, shown in Fig. 5(b), and along the $`zy`$ plane in Fig. 5(c). Over a distance of $`6\sigma `$ around the attachment position the isovalues of the probability distribution in the $`xy`$ plane delineate the minimum energy lines of the wall atoms LJ potential. In this model, the size of the wall atoms was chosen to be the same as those of the monomers and solvent particles $`1\sigma `$. In view of Fig. 5(b), one should expect that the structure of the polymer is quite sensitive to any modification in the details of the wall-fluid interaction, owing to either changes in the size of the wall atoms or details of the interaction potential.

In Fig. 6 we present a comparison of the radius of gyration $`R_g`$, as calculated from MD simulations in this work and that of Aust, Kröger and Hess Aust et al. (1999) (AKR) who studied a single free polymer in a bath of solvent molecules, at a variety of imposed shear rates. The potential used to describe the polymer and solvent were the same in both AKR’s work and ours; however AKR used a slightly higher density, $`\rho =0.85\sigma ^3`$, compared to our value of $`\rho =0.8\sigma ^3`$. The simulation of AKR used no walls, so the polymer was free to respond to the imposed shear so as to best lower the free energy of the system. Hence the usefulness of the comparison lies primarily in exploring the effect of the wall on the polymer. We see that at extremely low shear rates the values of the radii of gyration are quite comparable. As the shear rate increases the value of $`R_g`$ that we calculate becomes much larger than for the equivalent free polymer. This is due entirely to the attraction the polymer has with the wall.

Our results are in agreement with the experimental findings of Doyle et al. Doyle et al. (2000) for individual tethered DNA chains under shear flow. For a quantitative comparison with these experimental data we evaluated the Weissemberg number, Wi, defined as the product of the shear rate and the longest relaxation time of the polymer, i.e. the relaxation time at zero shear rate $`\tau _0`$. We calculate $`\tau _0`$ from the autocorrelation of the polymer extension at $`\dot{\gamma }=0`$ and obtain $`\tau _02000\tau `$. Also, the fractional extension is calculated by normalizing the polymer extension with its contour length: $`0.965\times (N1)`$, where $`N=60`$ is the the number of monomers and $`0.965\sigma `$ is the mean separation between two consecutive beads 222We checked that the monomer separation remains unchanged for the range of $`\dot{\gamma }`$ considered (it varies in less than $`0.005\sigma `$.). Using this value of $`\tau _0`$ we plot in Fig. 7 the normalized mean fractional extension along the flow direction versus the Weissemberg number, along with the the expermimental results of Doyle et al. The results obtained with both, the MD and hybrid simulations are in very good agreement with the experimental data for the range of shear rates considered here.

Figure 8 shows the end-to-end volume of the polymer, measured as the product of the three components of the end-to-end vector. This quantity gives an estimate of the space that the polymer explores during its motion. This volume increases for increasing shear rate and reaches a maximum value around $`\dot{\gamma }0.002\tau ^1`$. As the shear rate is further increased the volume accessible to the polymer decreases monotonically. This behavior of the end-to-end volume is quite similar to the findings of Doyle et al. Doyle et al. (2000) concerning the amplitude of the fluctuation of the chain extension. As the shear rate was varied, they found that fluctuations reached a maximum size at $`\mathrm{Wi}5.1`$. Using the estimate $`\tau _02000\tau `$, we find that the maximum end-to-end volume occurs at about $`\mathrm{Wi}4`$. In fact, the size of the fluctuations is determined by the magnitude of the volume made available by the polymer motion; or, in other words, larger fluctuations increase the explored volume. We shall present a more detailed comparison with the results of Doyle et al. Doyle et al. (2000) in future work.

There are computational costs to the hybrid method that are not present for classical MD. These include simulating the continuum regime and the calculations arising from the coupling procedure within the overlap region, e.g. particle insertion and deletion and the evaluation of the particle stress tensor. For the flow considered here, the solution of the continuum flow required around $`0.01\%`$ the time needed for a LJ force calculation. In general, the computational time spent in simulating the continuum region depends on the problem considered, but in any case it will always be much smaller than the MD force evaluation for the solvent. Furthermore, the calculation of the C-flow occurs once for every $`20`$ LJ force calculations, which ensures extra savings in computational time. As shown in Appendix B, the coupling protocol, within the overlap region, is very efficient: only $`0.01\%`$ of the total computational time was spent in particle insertion and deletion while around $`5\%`$ in the evaluation of the particle stress tensor. The hybrid code as tested here needs less than half the solvent particles, thus the overall savings in computational time is considerable.

## IV Conclusion

In this paper we have compared a newly developed hybrid MD/CFD code to a traditional MD simulation for a single polymer tethered to a wall undergoing shear flow in Couette geometry. We find that the two methods give comparable results for the conformation of the polymer within measured uncertainty.

Our results indicate that the coupling protocol of the hybrid code requires around $`5\%`$ of the computer time compared to the Lennard-Jones part of the code. Most of the CPU time devoted in the “coupling” protocol is spent in the evaluation of the particle momentum flux; while insertion and extraction of particles is rather fast, taking less than $`1\%`$ of the overall CPU time.

This implies that, compared with a traditional MD simulation, the amount of computational time saved by the hybrid scheme is proportional to the volume of the simulation that is described by the coarse-grained model (CFD). In traditional MD simulations of interfacial phenomena finite size effects significantly alter the local interfacial dynamics, and they can only be reduced by increasing the volume of the simulation box that surrounds the interfacial region of interest. This means that most of the computational cost is likely to be spent in the resolution of the bulk flow. In this paper we have shown that this drawback disappears when using a proper hybrid MD-CFD scheme. To that end, we considered a problem which is very sensitive to small changes in the surrounding fluid environment: the motion of a single tethered polymer under shear flow. The excellent agreement found in the comparisons with the full MD results indicates that the hybrid scheme indeed eliminates finite size effects even in relatively small systems. This means that hybrid simulations can be expected to significantly reduce the computational cost of appropriate interfacial problems.

Apart from the savings in CPU time, the hybrid scheme enables to us gather information from all relevant time and length scales, so it is well suited to treat multiscale problems where bulk fluid flow plays an important rôle; other examples include crystal growth from fluid phases, wetting phenomena and membrane dynamics under flow.

We regard the results of the present work as an encouraging sign for future simulations, and we plan to explore various selected interfacial systems in forthcoming research.

## V Acknowledgments

We acknowledge fruitful discussions on the hybrid formalism with Eirik Flekkøy. This research was supported by the European Commission through a Marie Curie Fellowship (HPMF-CT-2001-01210) and by the EPSRC RealityGrid project grant GR/R67699. R. D-B also acknowledges support from project BFM2001-0290.

## Appendix A Mass, longitudinal momentum and energy fluctuations

In this work the mean solvent flow carries no longitudinal momentum along the $`z`$-direction and has a constant mean density. However, we observe that the polymer motion induces density and longitudinal velocity fluctuations within the particle region that induce currents of mass and longitudinal momentum travelling along the simulation box. These perturbative currents have to be controlled at the C$``$P interface. We need to ensure that the mean mass flux across the $`z=l_{CP}`$ interface is zero, but in such a way that any pressure waves leave the simulation box once they reach the C$``$P interface. In other words, we need to prevent any pressure waves from bouncing back at the C$``$P interface in the MD region.

The average number of particles crossing the C$``$P interface per unit time is given by $`\dot{N}_{CP}=A\rho v_z_{CP}`$. Zero mass flux is ensured by equating this rate $`\dot{N}_{PC}`$ to the rate of insertion of molecules into the particle system Delgado-Buscalioni et al. (2003); Delgado-Buscalioni and Coveney (2003a). In the calculations presented here we used another control equation which provides a finer control on the particle density near the C$``$P interface. This approach is based on relaxing the local density at the C$``$P buffer to a prespecified value $`\rho _O`$,

$$\dot{N}_{CP}=\frac{V_{PC}}{\tau _m}\left(\rho _{CP}\rho _O\right)$$
(4)
where $`V_{CP}`$ is the volume of the C$``$P cell, $`\rho _{CP}`$ is its local the particle density averaged over $`\mathrm{\Delta }t_{av}`$ and $`\tau _m`$ is a relaxation time which controls the rate at which the density fluctuations within C$``$P cell are smoothed out. We set the value of $`\tau _m`$ slightly smaller than the time needed by a sound wave to cross the C$``$P cell ($`O(1)\tau `$). This procedure ensures that fluctuations carrying mass and longitudinal currents are damped at the C$``$P cell and do not bounce back to the inner part of the MD domain. According to Eq. (4), particles are extracted if $`\dot{N}_{PC}<0`$ and, as explained in Delgado-Buscalioni and Coveney (2003b), the first particles to be extracted are those closer to the C$``$P interface. If $`\dot{N}_{PC}>0`$, new particles are inserted with a velocity extracted from a Maxwellian distribution with mean velocity $`v_y=v_z=0`$ and $`v_x=\dot{\gamma }z`$ and temperature $`T=1.0`$. The insertion of particles in liquids is not a trivial task, however, and it is addressed by the usher algorithm for particle insertion Delgado-Buscalioni and Coveney (2003a). The value of $`\rho _O`$ in Eq.(4) was set to a slightly smaller value, $`\rho _O=0.65`$, than the mean density $`0.8`$. This reason for this choice is, first, to alleviate the computational cost of insertion (see Appendix B) and, second to reduce the amplitude of the ripples of the density profile at the C$``$P buffer, as shown below. For a liquid with $`\rho 0.8`$ the usher algorithm needs around 30 iterations to insert a LJ atom at a location where the potential energy equals the mean specific potential energy of the system Delgado-Buscalioni and Coveney (2003a), where each iteration corresponds to the evaluation of a single-particle force. If the density is decreased to $`0.65`$, it only needs about 15 iterations. Figure 6 compares the density profile resulting from using $`\rho =0.65`$ in Eq. (4) with that arising from a pure MD simulation. The “hybrid” density profile presents some ripples whose amplitude is damped after around $`3\sigma `$, whereas inside the P$``$C cell the hybrid density profile perfectly matches the density within the bulk.

As long as the fluid is isothermal and there are no mean pressure gradients, the mean energy flux across the C$``$P interface is zero. We therefore only need to guarantee that the specific energy of the newly inserted particles matches that of the ensemble. The kinetic energy is matched by inserting new particles with a Maxwellian distribution, as stated above. In order to match the potential energy, new particles are inserted at sites where the inter-particle potential energy equals the chemical potential of the system, thereby ensuring the Widom insertion criterion.

## Appendix B Computational cost of the hybrid method

We compare the computational cost of the coupling subroutines with those pertaining to the MD part of the hybrid scheme. This comparison was made using the gprof command available in the package of the f77 compiler. One of the parts of the hybrid scheme for which one may expect a certain cost in computational time is particle insertion. Table 1 presents some results obtained for different shear rates and values of the density $`\rho _O`$ in Eq. (4). Typically, Eq. (4) requires around 5 insertions per time interval $`\tau `$ and around 15 iterations per particle (each interation involving a single-force evaluation). Therefore, for a time step of $`\mathrm{\Delta }t_P=0.0075\tau `$, the insertion of new particles needs typically about one extra force evaluation per time step. This number is very small when compared with the number of force evaluations needed in the MD system, which is on the order of the number of particles $`N_p10^4`$. This estimate is consistent with our findings concerning the computational cost. As shown in table 1, in hybrid calculations using $`\rho _O=0.65`$, the time spent in the insertion/extraction subroutines was about $`1.5\times 10^4`$ times the time spent in the force evaluation and around $`0.9\times 10^5`$ if one includes the Verlet list evaluation. This performance confirms the extremely high efficiency of the usher algorithm for particle insertion.

As a matter of fact, the dominant cost of the hybrid scheme resides in the evaluation of the particle momentum flux. Its cost in CPU time was about 0.06 times the cost of the force plus Verlet list subroutines. We note that the implementation of this part of our code was not constructed in an efficient way because we evaluated the particle momentum flux at each MD time step. Considering that for the evaluation of $`j_p`$ we used measurements of $`j_p`$ separated by its decorrelation time, about $`0.06\tau `$ Delgado-Buscalioni et al. (2003), we could have measured the particle flux roughly every $`10`$ time steps and further reduced that ratio by a factor 10. Finally, the time needed to solve the diffusion equation in the continuum domain was very small compared with the MD force subroutine, by a factor of less than $`10^4`$. In general, the computational time required to solve the continuum system will, of course, depend on the specific problem solved.


format: text

text:
# Spatial patterns of desynchronization bursts in networks

## I Introduction

In this paper we study the synchronization of networks of coupled chaotic units that are nearly, but not exactly, identical. In particular, we will be concerned with the spatial patterns of desynchronization bursts that appear when this synchronization is present but intermittent.

When two or more identical dynamical systems are coupled, they can synchronize under appropriate circumstances. The synchronization of chaotic units has been studied extensively pecora1 ; pikovsky  and is of significance in biology elson -mosekilde , laser physics roy -uchida , and other areas hudson ; cuomo . At the same time, the importance of complex networks has been recently appreciated, and progress has been made towards their understanding, including characteristics that might help distinguish qualitatively different networks newman1 -newman2 . The dynamics of a network of coupled oscillators, and, in particular, its synchronization, has therefore emerged as a subject of great interest.

Pecora and Carroll pecora2  have proposed a model and analysis method (the master stability function) for the study of the stability of the synchronous state of networks of identical coupled chaotic units, and this technique has recently been extensively applied pecora3 ; takashi  to study the synchronization properties of different kinds of networks of identical noiseless chaotic units. These networks include small world smallworld  and scale-free networks barabasi .

The analysis of network synchronization by use of the master stability function technique has so far assumed all the units to be identical and noise-free, so that an exact synchronized state is possible. In practice, however, even if one strives to make the oscillators the same, they are still expected to have a small amount of parameter mismatch, and a small amount of noise is also expected to be present. Under such circumstances, it is known that the synchronization can be interrupted by sporadic periods of desynchronization (bursts). The bursts are typically caused by a periodic orbit that is embedded in the synchronized chaotic attractor and is unstable in a direction transverse to the synchronization manifold. This phenomenon is commonly referred to as bubbling ashwin -bhott2 , and has been studied extensively for two coupled oscillators rulkov ; pecora4 .

Our purpose in this paper is to study desynchronization bursts in networks of coupled chaotic nonidentical units. (Noise has a similar effect but will not be treated in this paper.) We will use the master stability function approach and, in order to account for the possibility of bubbling, we will also extend this approach to include the stability of embedded periodic orbits. In this case, the bursts have the added feature of having spatial patterns on the network, and we find that these patterns can be predicted from the network connectivity matrix. We will show how these bursts affect different parts of the network in different ways. In particular, we will see how adding connections in a ring can destabilize precisely those nodes that are the most connected, leaving other parts of the network substantially synchronized. (This a somewhat counterintuitive effect related to the fact that, in some cases, increasing the coupling strength destabilizes the synchronous state pecora2 ,pecora5 .)

Arbitrarily small amounts of mismatch will eventually, through the bubbling mechanism, induce desynchronization bursts. We will show that some of the spatial patterns of this possibly microscopic mismatch might get amplified to a macroscopic size in the bursts. We will discuss how one can use knowledge of the parameter mismatch of the dynamical units in the network to decrease the effective size of the mismatch driving the bursts, thereby improving the robustness of the synchronization.

If synchronization is desired, the network and the parameters should be constructed so that the synchronous state for the identical oscillator system is robustly stable (this implies the absence of noise or mismatch induced desynchronization bursts). Even then, the synchronization will not be perfect if the oscillators have parameter mismatch. We will describe the characteristics of the deviations from exact synchronization in terms of the mismatch and the master stability function.

This paper is organized as follows. In Sec. II we review the master stability function approach and apply it to the case of coupled Rössler units. We also discuss the bubbling mechanism by including the embedded periodic orbits in the master stability function analysis. In Sec. III we numerically consider particular networks as examples and show the resulting bursts and their spatial patterns. The patterns we obtain are long and short wavelength modes in a ring and localized bursts produced by strengthening of a single connection in a ring. In Sec. IV we study the effects of the spatial patterns of the mismatch in the development of the bursts. In Sec. V we study the deviations from the synchronous state caused by the mismatch when the synchronous state of the identical oscillator system is stable. In Sec. VI we summarize our conclusions.

## II Master stability function and bubbling

We now briefly review the master stability function approach introduced in pecora2 . Consider a system of $`N`$ dynamical units, each one of which, when isolated, satisfies $`\dot{X}_i=F(X_i,\mu _i)`$, where $`i=1,2,\mathrm{}N`$, and $`X_i`$ is the $`d`$-dimensional state vector for unit $`i`$. In pecora2  the parameter vectors $`\mu _i`$ are taken to be the same, $`\mu _i=\mu `$. Here, however, the parameter vectors $`\mu _i`$ are in general different for each unit, but we assume the difference, or mismatch, between them to be small. Generalizing the situation treated in Ref. pecora2  to the case where the individual units are not identical (i.e., the $`\mu _i`$ are not all equal), the system of coupled dynamical units is taken to be of the form

$$\dot{X}_i=F(X_i,\mu _i)g\underset{j=1}{\overset{N}{}}G_{ij}H(X_j),$$
(1)
where the coupling function $`H`$ is independent of $`i`$ and $`j`$, and the matrix $`G`$ is a Laplacian matrix ($`_jG_{ij}=0`$) describing the topology of network connections. For $`ij`$, the entry $`G_{ij}`$ is zero if oscillator $`i`$ is not connected to oscillator $`j`$ and nonzero otherwise. The nondiagonal entries of $`G`$ are determined by the connections, and the diagonal elements are the negative of the sum of the nondiagonal matrix elements in their row. The coupling constant $`g`$ determines the global strength of the coupling.

Assume first that all the dynamical units are identical, that is, $`\mu _i=\mu `$. We will refer to this situation as the idealized case. In this case there is an exactly synchronized solution $`X_1=X_2=\mathrm{}=X_N=s(t)`$ whose time evolution is the same as the uncoupled dynamics of a single unit, $`\dot{s}=F(s,\mu )`$. This convenient result arises because the Pecora-Carroll model uses the particular choice of coupling in (1) that ensures that the summation is identically zero when all of the $`X_j`$ are equal. We will denote this synchronization manifold, $`X_1=X_2=\mathrm{}=X_N`$, by $`M`$. This manifold is a $`d`$ \- dimensional surface within the $`Nd`$ \- dimensional phase space of Eq. (1).

The stability of the synchronized state can be determined from the variational equations obtained by considering an infinitesimal perturbation $`ϵ_i`$ from the synchronous state, $`X_i(t)=s(t)+ϵ_i(t)`$,

$$\dot{ϵ}_i=DF(s)ϵ_ig\underset{j=1}{\overset{N}{}}G_{ij}DH(s)ϵ_j.$$
(2)
Let $`ϵ=[ϵ_1,ϵ_2,\mathrm{},ϵ_N]`$ be the $`d\times N`$ matrix representing the deviation of the entire network from the synchronized state. In matrix notation, Eq. (2) becomes

$$\dot{ϵ}=DF(s)ϵgDH(s)ϵG^T.$$
(3)
While (3) allows for nonsymmetric coupling, we henceforth assume the coupling matrix $`G`$ to be symmetric, $`G=G^T`$. We write the symmetric matrix $`G`$ as $`G=L\mathrm{\Lambda }L^T`$, where $`\mathrm{\Lambda }`$ is the diagonal matrix of real eigenvalues $`\lambda _1,\lambda _2,\mathrm{},\lambda _N`$ of $`G`$ and $`L`$ is the orthogonal matrix whose columns are the corresponding real orthonormal eigenvectors of $`G`$ ($`L^TL=I`$). Define the $`d\times N`$ matrix $`\eta =[\eta _1,\eta _2,\mathrm{},\eta _N]`$ by $`ϵ=\eta L^T`$. Then Eq. (3) is equivalent to

$$\dot{\eta }=DF(s)\eta gDH(s)\eta \mathrm{\Lambda }.$$
(4)
Componentwise,

$$\dot{\eta }_k=\left(DF(s)g\lambda _kDH(s)\right)\eta _k.$$
(5)
The quantity $`\eta _k`$ is the weight of the $`k^{th}`$ eigenvector of $`G`$ in the perturbation $`ϵ`$. The linear stability of each ‘spatial’ mode $`k`$ is determined by the stability of the zero solution of (5). As a consequence of the condition $`_jG_{ij}=0`$, there is a special eigenvalue, $`\lambda =0`$, whose eigenvector is $`ϵ_N=[1,1,1,\mathrm{},1]`$, corresponding to perturbations in the synchronization manifold $`M`$. Since these are not perturbations from the synchronous state, the analysis is focused on the perturbations corresponding to nonzero eigenvalues.

By introducing a scalar variable $`\alpha =g\lambda _k`$, the set of equations given by (5) can be encapsulated in the single equation,

$$\dot{\eta }=\left(DF(s)\alpha DH(s)\right)\eta .$$
(6)
The master stability function $`\mathrm{\Psi }(\alpha )`$ pecora2  is the largest Lyapunov exponent for this equation for a typical trajectory in the attractor. This function depends only on the coupling function $`H`$ and the chaotic dynamics of an individual uncoupled element, but not on the network connectivity. The network connectivity determines the eigenvalues $`\lambda _k`$ (independent of details of the dynamics of the chaotic units). In the sense of typical Lyapunov exponents, the stability of the synchronized state of the network is determined by $`\mathrm{\Psi }_{}=sup_k\mathrm{\Psi }(g\lambda _k)`$, where $`\mathrm{\Psi }_{}>0`$ indicates instability. Thus the Pecora-Carroll model cleanly breaks the stability problem into two components, one from the dynamics (obtaining $`\mathrm{\Psi }(\alpha )`$) and one from the network (determining the eigenvalues $`\lambda _k`$).

In contrast to previous work using the master stability function technique, in this paper we are interested in the dynamics of systems in which a small parameter mismatch is present. (Even though in this paper our examples are restricted to the case of mismatch, we emphasize that the same type of bursting phenomenon is expected for identical oscillators if noise is present ashwin -rulkov .) Although the synchronization manifold $`M`$ present in the dynamics of the idealized system is, in general, not invariant for the system with mismatch, it still may provide a useful approximation to the dynamics in systems with small mismatch. If $`M`$ is stable for the idealized system, and the mismatch is small enough, then trajectories near $`M`$ will tend to stay near $`M`$, and we regard the vicinity of $`M`$ to be the “synchronized” state. However, stability of $`M`$ in the idealized case of identical oscillators is not sufficient to guarantee robust synchronization in a real system where the oscillators are not identicalashwin -rulkov . While in the vicinity of the synchronization manifold $`M`$, a typical trajectory will eventually follow very closely a periodic orbit embedded in the attractor of the idealized system. Some of these periodic orbits may be unstable in a direction transverse to $`M`$. When in the vicinity of a transversally unstable periodic orbit, mismatch (or noise) will cause the trajectory have a component in the direction transverse to $`M`$ and hence to leave the vicinity of the synchronization manifold $`M`$. If there are no other attractors, the trajectory will eventually return to the vicinity of $`M`$, and the process will repeat, the result being bursts of desynchronization sporadically interrupting long intervals of near synchronization. This type of dymanics is called bubbling ashwin .

Thus, in the presence of mismatch (or noise), to determine the robustness of synchronization, it is necessary to determine the transverse stability of the embedded periodic orbits for the noiseless system of identical oscillators. For coupling as in (1), this analysis is independent of the network, and such analyses have been carried out before, e.g., for the analysis of two coupled oscillators in Ref. pecora4 . Equation (6) can be used as before to construct the master stability function for each periodic orbit, if the appropriate periodic trajectories are inserted for $`s(t)`$ in (3).

As an example, in this paper we work with the Rössler system rossler :

$`\dot{x}=(y+z),`$ (7)
$`\dot{y}=x+ay,`$
$`\dot{z}=b+z(xc).`$
In terms of our previous notation, $`d=3`$, $`\mu =[a,b,c]^T`$, and $`X=[x,y,z]^T`$. We choose the parameters of the idealized system to be $`a=b=0.2`$, $`c=7`$. For these parameters, the system has a chaotic attractor (see Fig. 1).

We found the periodic orbits embedded in this attractor up to period five, and performed the analysis described above on them. We found these orbits by looking at the Poincare surface of section $`\{y=0,x<0\}`$. To a good approximation, in this surface of section the dynamics is well described by a one dimensional map $`x_{n+1}=f(x_n)`$, which we approximated using a polynomial fit. From this approximation to $`f`$, we determined periodic orbits of period $`p`$ by using Newton’s method to find the roots of $`x=f^p(x)`$, where $`f^p`$ denotes the $`p`$ times composition of $`f`$. We found one period $`1`$ orbit, one period $`2`$ orbit, two period $`3`$ orbits, three period $`4`$ orbits, and four period $`5`$ orbits. Using coupling through the $`x`$ coordinate,

$$H([x,y,z]^T)=[x,0,0]^T,$$
(8)
we obtained a stability function $`\mathrm{\Psi }(\alpha )`$ for each orbit, the largest of which will determine if the synchronization is robust. Results are shown in Fig. 2. For all values of $`\alpha `$, we found that the master stability function corresponding to the period $`1`$ orbit (thick dashed curve) is larger than that for a typical chaotic orbit (thick continuous curve), as well as those for the other periodic orbits we have found (several of which are shown as thin curves).

Based on the discussion above, bubbling induced bursting should occur whenever the master stability function for a typical chaotic orbit in the attractor is negative for $`\alpha =g\lambda _k`$ and all $`k`$, while the period one orbit has positive master stability function for $`\alpha =g\lambda _k`$ for some value of $`k`$. Denoting the master stability function for a typical chaotic orbit by $`\mathrm{\Psi }_0(\alpha )`$ (Thick continuous curve in Fig. 2) and for the period one orbit by $`\mathrm{\Psi }_1(\alpha )`$ (Thick dashed curve in Fig. 2), the bubbling region of $`\alpha `$ corresponds to $`\mathrm{\Psi }_0(\alpha )<0`$, $`\mathrm{\Psi }_1(\alpha )>0`$. In our example, this region corresponds to $`0.16<\alpha <0.48`$ or $`3.8<\alpha <4.5`$. The range $`0.48<\alpha <3.8`$ will be referred to as the stable region, and the remaining zone will be called the unstable region.

If a network of slightly mismatched chaotic systems coupled according to 1is to be robustly synchronizable without bursts of desynchronization, $`g\lambda _k`$ must lie in the stable region for all $`k`$, where $`\lambda _k`$ is the $`k`$th eigenvalue of $`G`$. If $`g\lambda _k`$ lies in the stable region for some $`k`$ and in the bubbling region for other $`k`$, then bubbling will typically occur.

## III Examples

In this Section we provide examples of spatially patterned bursting by considering different configurations of the chaotic units. We will first work with the units connected in a ring with each connection of equal strength. The Laplacian matrix $`G`$ for this arrangement is

$$G=\left(\begin{array}{ccccccc}2& 1& 0& 0& \mathrm{}& 0& 1\\ 1& 2& 1& 0& \mathrm{}& 0& 0\\ 0& 1& 2& 1& \mathrm{}& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 1& 0& \mathrm{}& 0& 0& 1& 2\end{array}\right),$$
(9)
and its eigenvalues are given by $`\lambda _k=4\mathrm{sin}^2(\frac{\pi k}{N})`$. Since $`\lambda _k=\lambda _{Nk}`$, each eigenvalue has multiplicity two, with the exception of $`\lambda _N=0`$, and, if $`N`$ is even, $`\lambda _{\frac{N}{2}}=4`$. The matrix $`G`$ is shift invariant, that is, its entries satisfy, modulo $`N`$, $`G_{i,j}=G_{0,ij}`$. Under these conditions, the diagonalization procedure described above corresponds to a discrete Fourier transform pecora5 . For the eigenvalue $`\lambda _k`$ we choose the eigenvector $`w_k`$ given by $`w_k[\mathrm{sin}(\frac{2\pi jk}{N})]_{j=1}^N`$ for $`1k<\frac{N}{2}`$, and by $`w_k[\mathrm{cos}(\frac{2\pi jk}{N})]_{j=1}^N`$ for $`\frac{N}{2}kN`$. (Due to the degeneracy of the eigenvalues in this case, there is some arbitrariness in choosing the eigenvectors.) Thus, the longest wavelength modes have the smallest eigenvalues, and viceversa.

### III.1 Long wavelength burst

First we consider a case in which bursting of the longest wavelength mode occurs. We consider $`N=12`$ and $`g=0.71`$. With these values, the longest wavelength mode corresponds to $`\alpha =g\lambda _10.19`$. This value is in the bubbling region, and all other modes are in the stable region.

To introduce heterogeneity in the dynamical units, we imagine that we have mismatch predominantly in one of the parameters, say $`a`$. We simulate this mismatch by adding random perturbations to the parameter $`a`$ of each oscillator. These perturbations are uniformly distributed within a $`\pm 0.5\%`$ range; i.e., $`a_i`$ is chosen randomly in the interval $`[0.995a,1.005a]`$, where $`a`$ is the parameter value of the unperturbed system ($`a=0.2`$). The parameters $`b`$ and $`c`$ were taken to be the same for each oscillator, $`b_i=b=0.2`$, $`c_i=c=7`$. How a particular choice of the mismatch affects the bubbling process will be discussed in Section IV.

We solved the $`12`$ coupled differential equations (Eq. (1)) with the initial conditions chosen near the attractor in the synchronization manifold. In Fig. 3 we plot the quantity $`x_1x_6`$ for $`1000t1600`$.

Most of the time, this variable is close to zero, as expected if the oscillators are synchronized. Approximately at the time $`t=1380`$, this difference grows, reaching magnitudes close to $`3`$. By time $`t=1500`$, the difference has decreased and is again close to zero.

To confirm the mediating role of the embedded unstable periodic orbits in the development of the desynchronization burst, we show in Fig. 4 a plot of $`x_1`$ versus $`y_1`$ from $`t=1372`$ to $`t=1392`$, which is near the start of the burst. During this time, the trajectory closely follows the period $`1`$ orbit, which is the most transversally unstable of the periodic orbits. Similar observations have been previously reported for two coupled chaotic systems pecora4 .

Finally, in Fig. 5 we plot $`x_jx_{j1}`$ as a function of $`j`$, the oscillator index, for $`t=1360`$ (open triangles), $`t=1385`$ (open circles), and $`t=1410`$ (open squares). The desynchronization burst can be observed developing mainly at the longest possible wavelength.

When subsequent bursts were studied in the same way, it was found that the phase of the long wavelength burst assumed only one value. This is due to the fact that the mismatch is ‘frozen’, that is, each oscillator has a given set of parameters which differs by a given amount from the mean values. This fixed spatial heterogeneity favors certain spatial patterns over others. We will discuss this in more detail in Section IV.

### III.2 Short wavelength burst

Short wavelength bursting can be expected, for example, when $`N=8`$ and $`g=1.09`$. In this case the value of $`\lambda _k`$ corresponding to the shortest wavelength mode yields $`g\lambda _k=4.36`$, which is in the bubbling region, while all the other modes are in the stable region. In this case the observation of the bursts is more difficult, as the transversal instability of the orbits and the transversal stability of the attractor are less pronounced (compare $`\mathrm{\Psi }(4.36)`$ for this case vs. $`\mathrm{\Psi }(0.19)`$ for the previous example in Fig. 2). Accordingly, the perturbations of the parameter $`a`$ were made larger, with perturbations randomly chosen with uniform density within a $`\pm 6\%`$ range of the ideal values of the parameter ($`a=0.2`$). In principle this is not necessary, as a burst will eventually occur after long enough time. In practice, however, it is necessary to reduce the waiting time to a reasonable value. As before, the coupled equations were solved with an initial condition on the synchronization manifold. In Fig. 6 we show $`y_2y_1`$ as a function of time for one choice of initial conditions.

The difference $`y_2y_1`$ is usually positive and of magnitude close to $`1`$. This asymmetry is not a surprise since the oscillators are slightly different. For the relatively large value of the mismatch used, this is the “synchronized state”. It is seen in Fig. 6 that the difference $`y_2y_1`$ increases rapidly at around $`t18400`$, and soon reaches values higher than $`7`$. It remains large for a longer time than in the case of the long wavelength burst (see Fig. 3) and decays more slowly as well. This is in qualitative agreement with the smaller absolute values of the master stability functions for the short wavelength mode, both for typical orbits on the attractor and for the periodic orbits.

In Fig. 7 we plot $`y_jy_{j1}`$ as a function of $`j`$, the oscillator index, for $`t=17970`$, $`t=18270`$ and $`t=18570`$. As expected, the burst mainly affects the shortest wavelength mode. In order to perform a more quantitative assessment of these observations, we define $`\xi _k=\{([\eta _k]_x)^2+([\eta _{Nk}]_x)^2\}^{\frac{1}{2}}`$ for $`1k<\frac{N}{2}`$ and $`\xi _{\frac{N}{2}}=\left|[\eta _{\frac{N}{2}}]_x\right|`$, where $`[\eta _k]_x`$ is the $`x`$ component of the three dimensional vector $`\eta _k`$. The quantity $`\eta _k`$ is proportional to the magnitude of the vector $`[\eta _{1,k}\mathrm{sin}\left(\frac{2\pi jk}{N}\right)+\eta _{1,Nk}\mathrm{cos}\left(\frac{2\pi jk}{N}\right)]`$, the component in the Fourier decomposition of the perturbation $`ϵ`$ associated with modes $`k`$ and $`Nk`$. Thus, the quantity $`\xi _k`$ represents the weigth of the modes associated to the eigenvalue $`\lambda _k`$, which have, for $`1k\frac{N}{2}`$, a wavelength of $`\frac{1}{k}`$’th of a full wavelength. In Fig. 8, we plot as a function of time the quantities $`\xi _k^2`$ for $`k=1,2,3,4`$. The short wavelength mode ($`k=4`$, upper curve) is dominant during the burst.

### III.3 Localized burst

In the above examples all links had equal weights. As an example of a case with unequal link weights we consider the case where the previous network is modified by doubling the strength of one of the links. Let the link whose strength is doubled be the link that connects nodes $`p`$ and $`p+1`$. For example, for $`p=4`$, $`N=8`$, this yields the Laplacian matrix

$$G=\left(\begin{array}{cccccccc}2& 1& 0& 0& 0& 0& 0& 1\\ 1& 2& 1& 0& 0& 0& 0& 0\\ 0& 1& 2& 1& 0& 0& 0& 0\\ 0& 0& 1& 3& 2& 0& 0& 0\\ 0& 0& 0& 2& 3& 1& 0& 0\\ 0& 0& 0& 0& 1& 2& 1& 0\\ 0& 0& 0& 0& 0& 1& 2& 1\\ 1& 0& 0& 0& 0& 0& 1& 2\end{array}\right),$$
(10)
Adopting the analysis technique of Ref. sott , we can show that such an enhanced connection has the consequence that the largest eigenvalue of $`G`$ corresponds to an eigenfunction that is exponentially localized to the region near the strong connection. That is, for large $`N`$, the components of this eigenfunction decay exponentially as the distance between the localized region and the node corresponding to a component increases. Using the ideas of Ref. sott , we now provide this analysis. The equations for the eigenvector $`w`$ and eigenvalue $`\lambda `$ are

$`2w_{p+1}w_{p1}+3w_p=\lambda w_p,`$ (11)
$`w_{p+2}2w_p+3w_{p+1}=\lambda w_{p+1},`$
$`w_{j1}w_{j+1}+2w_j=\lambda w_j,`$
for, respectively, nodes $`p`$, $`p+1`$ and $`j`$ different from $`p`$ or $`p+1`$. Since nodes $`p`$ and $`p+1`$ are identical, we can assume, for $`k>0`$, $`w_{p+1+k}=\pm w_{pk}`$. Furthermore, we propose the ansatz $`w_{p+1+k}t^k`$, for $`k>0`$. This will be a good approximation if the mode is localized (i.e., $`|t|<1`$), and the network is big enough that $`|t|^{\frac{N}{2}}1`$. In the antisymmetric case, $`w_{p+1+k}=w_{pk}`$, Eqs. (11) yield,

$`5t=\lambda ,`$ (12)
$`tt^1+2=\lambda `$
which gives

$$t=\frac{1}{3}\text{ , }\lambda =\frac{16}{3}.$$
(13)
Compare this eigenvalue with the largest eigenvalue for the network in which all links have equal strength, which has a value of $`4`$. The symmetric solution, $`w_{p+1+k}=w_{pk}`$, yields $`t=1`$ and $`\lambda =0`$, corresponding to the eigenvector $`[1,1,\mathrm{}1]`$ of perturbations in the synchronization manifold. The smallest nonzero eigenvalue remains unchanged.

As an example, we show the localized desynchronization bursts produced by one of these strengthened connections for the case $`N=8`$, corresponding to $`G`$ given by (10) and the illustration in Fig. 9. The parameters of the idealized system are again $`a=b=0.2`$, and $`c=7`$, with a coupling strength of $`g=0.79`$.

It is remarkable that despite the small number of nodes, the actual localized eigenvector and eigenvalue agree well with (13) ($`\lambda =5.334\mathrm{}`$ and $`\frac{w_6}{w_5}=0.334\mathrm{}`$).

In Fig. 10 we show $`x_5x_4`$ as a function of time.

As in the short wavelength case, the burst is not very sharp due to the small magnitude of the transversal Lyapunov exponents. Nevertheless, it can be seen that the difference $`x_5x_4`$ increases approximately at $`t=9000`$ and returns to a relatively small value after reaching values considerably above the average.

In Fig. 11a we plot the difference between the $`x`$ coordinate of node $`j`$ and its mean over all nodes, $`x_j\overline{x}`$, where $`\overline{x}=\frac{1}{N}_{j=1}^Nx_j`$, as a function of the oscillator index $`j`$, for $`t=8750`$ (open triangles), $`t=9000`$ (open circles), and $`t=9250`$ (open squares). In Fig. 11b we show the localized eigenvector of the Laplacian $`G`$ found numerically.

As discussed before, the desynchronization burst affects mainly nodes $`4`$ and $`5`$ (those which share the strengthened connection) and the ones adjacent to them. Nodes $`1`$,$`2`$,$`7`$ and $`8`$, however, maintain approximate synchronization during the burst.

In Fig. 12 we show the mode weights corresponding to the $`x`$ coordinate as a function of time. The top curve corresponds to $`[\eta _4]_x^2`$ (for the localized mode), and the curves close to the horizontal axis to $`[\eta _k]_x^2`$, $`k4`$, for the other modes. (The degeneracy of the eigenvalues is broken by the strengthened connection, so we do not combine $`[\eta _k]_x`$ and $`[\eta _{Nk}]_x`$ as before.) Confirming the qualitative similarity between the eigenvector and the spatial pattern of the desynchronization burst observed in Fig. 11, the weight corresponding to the localized eigenvector is seen to be dominant during the period of time in which the burst occurs.

## IV Effects of the mismatch spatial patterns

In this section we will discuss the effects that the mismatch spatial patterns have on the development of the desynchronization bursts. For these purposes, it will be convenient to rewrite Eq. (1) in the form

$$\dot{X}_i=\overline{F}(X_i)g\underset{j=1}{\overset{N}{}}G_{ij}H(X_j)+Q_i(X_i),$$
(14)
where $`\overline{F}(X_i)=F(X_i,\overline{\mu })`$ with $`\overline{\mu }=\frac{1}{N}_{j=1}^N\mu _j`$, and $`Q_i(X_i)=F(X_i,\mu _i)\overline{F}(X_i)`$. The term $`Q_i`$ represents the effect of the mismatch and is assumed to be small. As before, we linearize around the synchronous state to get

$$\dot{ϵ}_i=D\overline{F}(s)ϵ_ig\underset{j=1}{\overset{N}{}}G_{ij}DH(s)ϵ_j+Q_i(s),$$
(15)
where we have discarded terms of order $`Qϵ`$. With the previous notation and $`Q=[Q_1,Q_2,\mathrm{}Q_N]`$, we obtain after the diagonalization

$$\dot{\eta }_k=\left(D\overline{F}(s)g\lambda _kDH(s)\right)\eta _k+(QL)_k,$$
(16)
where $`(QL)_k`$ is the $`k`$’th column of the $`d\times N`$ matrix $`QL`$. In the ring with equal coupling along each link, the diagonalization procedure corresponds to a Fourier transform. In this case, we see that the mismatch affects the different modes according to the weigth, $`(QL)_k`$, of this particular mode in its Fourier expansion. In other cases, for example in the localized eigenvector, the strength of the mismatch affecting the localized mode is proportional to the weigth of the localized eigenvector in the eigenvector decomposition of the mismatch. We will now discuss two applications of these results.

### IV.1 Amplification of mismatch patterns when modes with the same eigenvalue burst

We have shown that the modes of the mismatch force the corresponding modes of the deviations from the synchronous state. When bubbling induced bursting is expected, the size of the mismatch determines the average time between bursts bhott2 . Thus, the size of the mismatch component in mode $`k`$ determines the average interburst time when that mode is in the bubbling regime.

When the spectrum of the matrix $`G`$ is degenerate, the spatial modes of the mismatch play an extra role. All the modes sharing the same eigenvalue $`\lambda `$ have the same stability properties, and thus, when the corresponding value $`g\lambda `$ is in the bubbling zone, all eigenvectors with this eigenvalue are equally likely to appear. The only difference between these modes is the strength with which they are forced, which is determined by the mismatch component in that mode as shown in Eq. (16) (or, if noise is present, by the noise component in that mode).

An example of this situation is the ring with connections of equal strength in the long wavelength bursting scenario. Since the ring is invariant with respect to rotations, the phase of the long wavelength oscillations can not be determined only from the network and dynamics part of the problem. The two modes with the longest wavelength (corresponding to sinusoidal and cosinusoidal oscillations) have the same eigenvalue. It is the mismatch that in this case determines the phase of the long wavelength burst.

We will show how one can determine the phase of the long wavelength desynchronization burst in the case of coupled Rössler systems in a ring with equal coupling along each link. For this system, the mismatch vector $`Q_j(X_j)`$ is given by

$$Q_j([x_j,y_j,z_j]^T)=\left(\begin{array}{c}0\\ y_j\delta a_j\\ \delta b_jz_j\delta c_j\end{array}\right),$$
(17)
where $`\delta a_j=a_j\overline{a}`$ und similarly for $`\delta b_j`$ and $`\delta c_j`$. We define $`_k(u)=_{j=1}^Nu_j\widehat{w}_j^k`$, where $`\widehat{w}_j^k`$ is the normalized $`j`$’th component of the $`k`$ eigenvector described at the beginning of Section III. With this convention, the term $`(QL)_k`$ in equation (16) is given by

$$(QL)_k=\left(\begin{array}{c}0\\ y_k(\delta a)\\ _k(\delta b)z_k(\delta c)\end{array}\right).$$
(18)
Here $`\delta a=[\delta a_1,\delta _2,\mathrm{},\delta _N]`$ and similarly for $`\delta b`$, $`\delta c`$, and $`y,z`$ are the trajectories around which the linearization was made.

We consider the case in which mismatch in one parameter is dominant, for example $`a`$. The mismatch in the parameters $`b`$ and $`c`$ will be assumed negligible compared with that in $`a`$, so that $`\delta b`$, $`\delta c\delta a`$. In this case, only the second component of (18) is of relevance. Thus modes $`\eta _1`$ and $`\eta _{N1}`$ are excited with a strength proportional, respectively, to $`_1(\delta a)`$ and $`_{N1}(\delta a)`$; see (16). The magnitude of $`\eta _k`$ will be proportional to $`_k(\delta a)`$, and thus the excitation of the long wavelength mode (which is the only one for which perturbations grow) is proportional to

$$_1(\delta a)\mathrm{sin}\left(\frac{2\pi j}{N}\right)+_{N1}(\delta a)\mathrm{cos}\left(\frac{2\pi j}{N}\right)$$
(19)
$$\mathrm{sin}\left(\frac{2\pi j}{N}+\varphi \right),$$
(20)
where $`\mathrm{tan}\varphi =_{N1}(\delta a)/_1(\delta a)`$.

We now show results of numerical simulations illustrating the above. The parameters $`N`$ and $`g`$ will be as in the long wavelength example in the previous section. We use the same random set of perturbations used in that example. As described above, we obtained the phase $`\varphi `$ of the long wavelength component of the vector $`\delta a_i`$. In Fig. 13 we plot $`y_jy_{j1}`$ for different times during a burst (filled symbols). In the same Figure, we plot a scaled version of $`\mathrm{sin}\left(\frac{2\pi j}{12}+\varphi \right)\mathrm{sin}\left(\frac{2\pi (j1)}{12}+\varphi \right)`$ (open circles). The phase of the desynchronization burst is in agreement with that of the long wavelength component of the mismatch.

When the mismatch affects predominantly one parameter as in this case, the phase of the bursts can be predicted as described above. When mismatch in different parameters is comparable, the phases of the long wavelength modes of the different parameter mismatches compete and the bursts develop with one of these phases or with a combination of them.

It must be emphasized that this analysis is possible only when there is a degeneracy of the eigenvalues. For example, the location of the localized burst can not be determined in this way, as it is fixed in the position of the strengthened link. In this case, the mismatch component in the localized mode would only affect the average time between bursts.

### IV.2 Artificial supression of unstable modes using knowledge of the mismatch

We will now discuss another consequence of Eq. (16). We imagine a situation where we are given a number of nearly identical oscillators that we are to connect in a network which we desire to be in synchronism as much as possible. Furthermore, we imagine that, through measurements made individually on each oscillator, we are aware of the amount of mismatch in each oscillator. The question we address is this: Using our knowledge of the individual mismatches, how should we arrange the oscillators in the network so as to best supress the frequency of desynchronism bursts? To answer this question, we note that, according to the previous discussion, we should reduce the mismatch component in the mode which is in the bubbling region. Since the size of the mismatch affects the average interburst time bhott2 , reducing this component is desirable if one wants to improve the quality of the synchronization. This can be accomplished by judiciously arranging the dynamical units so that the $`k`$’th mode of the mismatch is minimized when the corresponding value $`g\lambda _k`$ is in the bubbling region. For example, to supress long wavelength bursts, one may arrange the units so that the parameter errors alternate above and below the mean. To supress the localized bursting described in the previous section, one could arrange the units so that those with the more similar parameters are the ones in the region of the strengthened connection.

As a concrete example, we test this idea using simulations for the case of short wavelength bursting presented in the previous Section. We again assume for simplicity that mismatch in the parameter $`a`$ is dominant. We generate random perturbations in the parameter $`a`$ within a $`\pm 6\%`$ range of the value $`a=0.2`$, as explained in the previous section. With this set of parameters given, we set up the dynamical units in the ring using two different permutations of their positions. One of them ($`a_s`$) has a smaller and the other ($`a_l`$) a larger short wavelength component $`_4(a)`$ than the original random sequence. The ratio $`_4(a_l)/_4(a_s)`$ is approximately $`15`$. In Fig 14 we plot $`x_1x_2`$ as a function of time for configuration $`a_{large}`$ (top curve) and for configuration $`a_{small}`$ (bottom curve).

The difference $`x_1x_2`$ is much smaller in the former case than in the latter, roughly by a factor of $`15`$, as can be expected from the ratio $`_4(a_{large})/_4(a_{small})`$. This qualitative example illustrates how one can use knowledge of the mismatch to supress undesired instabilities.

## V Spatial patterns of deviations from the stable synchronous state

So far, we have concentrated in the case in which the value of $`g\lambda _k`$ is in the bubbling regime for one mode $`k`$ and in the stable regime for the other modes, so that desynchronization bursts occur sporadically. As we have seen, these bursts present spatial patterns on the network.

If synchronization is desired, one would might try to avoid the bubbling regime by designing the network and adjusting the coupling strength so that all the modes lie in the stable zone. One would also strive to reduce the mismatch, but as mentioned before, there are practical limitations on how much one can make the oscillators exactly the same.

If $`\mathrm{\Psi }(g\lambda _k)`$ is negative for all modes (indicating transversal stability of the synchronous state) one can have, depending on the degree of transversal stability, fair synchronization even with relatively large amounts of mismatch. If one is to operate under such conditions, it is important to know the characteristics of the deviations from the synchronous state.

Thus we ask in this scenario: How large are the spatial patterns of the deviations from the synchronous state, and how does this depend on the mismatch and on the degree of transversal stability?

The spatial modes of these deviations obey Eq. (16). In the absence of the term $`(QL)_k`$, the zero solution is stable, and typical perturbations from it decay, having a negative Lyapunov exponent given by $`h_k\mathrm{\Psi }(g\lambda _k)`$. The first term in the right hand side of Eq. (16) can be thougt of as a damping term with a damping rate given by $`h_k`$, and the second term, $`(QL)_k`$, as a forcing term. Since we are considering the stable case, these two factors, on average, cancel each other. By definition, the Lyapunov exponent for the system without mismatch is given by $`h_k=\frac{\eta _k^T(D\overline{F}g\lambda _kDH)\eta _k}{\left|\eta _k\right|^2}`$, where the angle brackets indicate time average. Assuming a solution $`\eta _k`$ of the system with mismatch to yield the same value of this time average, we left multiply Eq. (16) by $`\eta _k^T\left|\eta _k\right|^2`$ and average to obtain

$$\left|h_k\right|\frac{\eta _k^T(QL)_k}{\left|\eta _k\right|^2}\frac{\left|(QL)_k\right|}{\left|\eta _k\right|},$$
(21)
where the angle brackets indicate time average. This leads to the following rough estimate,

$$\left|\eta _k\right|\frac{\left|(QL)_k\right|}{\left|h_k\right|}$$
(22)
As an example we consider Rössler units in a ring with all connections of equal strength. We choose $`N=8`$, $`g=0.6`$ ($`\mathrm{\Psi }(g\lambda _k)<0`$ for all values of $`k`$). Furthermore, we add a random perturbation to the parameter $`a`$ of each oscillator chosen uniformly from within a $`\pm 0.1\%`$ range of $`a=0.2`$.

In Fig. 15 we show, for $`k=1,\mathrm{},7`$, the quantities $`\left|\eta _k\right|`$ (squares), $`\left|(QL)_k\right|`$ (triangles), and $`\frac{\left|(QL)_k\right|}{\left|h_k\right|}`$ (stars).

The magnitudes of the forcing term for the different modes ($`\left|(QL)_k\right|`$) span roughly two orders of magnitude, and the magnitude of the response ($`\left|\eta _k\right|`$) looks roughly proportional to the latter. When the forcing term is corrected by dividing it by the magnitude of the corresponding Lyapunov vector $`\left|h_k\right|`$, the resulting quantity ($`\frac{\left|(QL)_k\right|}{\left|h_k\right|}`$) matches very well the observed response.

## VI Conclusions

We have studied the stability properties of the synchronized state in a network of coupled chaotic dynamical units when these have a small heterogeneity. We have shown that when the dynamical units that are coupled in a network are sligthly different, the synchronized state can be interrupted by large infrequent desynchronization bursts for some values of the parameters. The range of the parameters for which this phenomenon is expected can be obtained by performing a master stability function analysis of the chaotic attractor and of the periodic orbits embedded in it.

The desynchronization bursts are induced by the bubbling phenomenon, and have spatial patterns on the network. These spatial patterns can be predicted from the eigenvectors of the Laplacian matrix $`G`$ and the master stability functions mentioned above. We showed examples illustrating the development of bursts with spatial patterns. One of our examples showed that the strengthening of a single connection might destabilize the nodes near this connection, while leaving the rest of the network approximately synchronized.

Direct measurement of the parameter mismatch in the elements of a network might prove useful. We discussed how this knowledge could be used to reduce the frequency of bursts and to predict the relative weights of different spatial patterns in a burst. We also discussed how one could, from knowledge of the mismatch and of the master stability function, describe the spatial patterns and magnitude of the deviations from the synchronized state when the synchronization of the corresponding identical unit system is robust.

We emphasize that although we did not discuss the effects of noise, the phenomenon described in this paper also occurs for noisy identical oscillators. Desynchronization bursts with spatial patterns are expected for noisy, identical oscillators if one has them for noiseless, nonidentical oscillators. The difference is that the parameter mismatch is always ‘frozen’, in the sense that the mismatch is always the same for each oscillator, whereas for noise this is not the case.

Acknowledgements: This work was sponsored by ONR (Physics) and by NSF (contracts PHYS 0098632 and DMS 0104087).


text:
# Note on the Robustness of the Neutrino Mass Bounds from Cosmology

## I Introduction

Recently, there has been a lot of progress in neutrino physics. Measurements of solar and atmospheric neutrino fluxes have yielded tight bounds on the mass differences between the different neutrino mass eigenstates (see e.g. rev  for recent reviews). The bounds on the absolute masses from direct measurements, however, are much weaker. Currently, the bounds on the absolute neutrino masses coming from indirect cosmological considerations are much stronger. In particular, given the recent high precision data on cosmic microwave background (CMB) anisotropies and the increased information about the galaxy power spectrum based on the recently completed large-scale galaxy redshift surveys, it is possible to deduce tight bounds on the sum of the neutrino masses list . These analyses, however, make use of certain assumptions about the evolution of the early Universe. Specifically, they make use of theoretical predictions for the power spectra of the galaxy distribution, derived under the assumption that the primordial spectrum of cosmological perturbations was purely adiabatic and coherent, as predicted in simple single field models of inflation. In this paper we show that the bounds derived in list  are robust against the addition of a contribution to the primordial power spectrum which is nonadiabatic and incoherent, namely a contribution coming from a distribution of topological defects.

The cosmological limits on the sum of the neutrino masses are derived by combining the high precision data of CMB anisotropies with the measurements of the galaxy power spectrum which have recently improved in accuracy upon completion of the 2dF (two degree field) large scale structure survey. Given the assumption that the primordial power spectrum is purely adiabatic and coherent (as predicted in simple single scalar field models of inflation), the CMB angular power spectrum determines to good accuracy the slope of the primordial power spectrum and the cosmological parameters (see e.g. parsrev  for recent reviews). The most sensitive dependence on the neutrino masses comes from the smallest scales, those which are not (yet) probed by CMB anisotropies but rather by the galaxy power spectrum. The reason for this dependence lies in the phenomenon of neutrino free streaming nfs . Because of their small mass, the neutrinos have large velocities at $`t_{eq}`$, the time of equal matter and radiation. Hence, neutrino density perturbations on comoving scales smaller than the neutrino free streaming length (the distance the neutrinos travel in one Hubble expansion time) at $`t_{eq}`$ are suppressed (the suppression factor increases exponentially as a function of the wavenumber). Hence, if the contribution of the neutrinos to $`\mathrm{\Omega }`$ (the energy density in units of the critical density) is too large, there will be insufficient power to explain the observed magnitude of the galaxy power spectrum. A lower bound on the galaxy power spectrum thus translates into an upper bound on the contribution of neutrinos to $`\mathrm{\Omega }`$, and thus to an upper bound on the sum of the neutrino masses. Note that the most stringent constraints come from the smallest length scales for which the galaxy power spectrum can be reliably determined.

However, the assumption that the primordial spectrum of fluctuations is purely adiabatic and coherent is a very restrictive one. Even in the context of scalar field driven inflationary cosmology, as soon as one considers models with more than one scalar field, it is possible to obtain a contribution of entropy fluctuations LindeMukh . Many inflationary models based on grand unified theories (see e.g. Rachel ) or on the brane inflation scenario (see e.g. Tye ) predict a contribution to the power spectrum from cosmic strings, yielding a contribution of isocurvature and incoherent fluctuations. It has been known for a long time BKST ; BKT  that the transfer function which relates the primordial power spectrum to the present power spectrum changes dramatically on small distance scales if a contribution of seed perturbations such as those created by cosmic strings is added. The basic point is that the cosmic strings constitute density fluctuations seeds which are not erased by neutrino free streaming. Thus, in a seed model, the accretion of neutrinos on small scales is delayed (and thus reduced in amplitude) but not prevented. Hence, it appears at first sight that even a small addition of cosmic string seed fluctuations to the primordial power spectrum might dramatically loosen the cosmological bounds on the sum of the neutrino masses. Here, we shall demonstrate that this is in fact not the case.

The outline of our analysis is as follows. We will first consider the largest contribution of cosmic strings to the angular power spectrum on large angular scales which is consistent with the current data on CMB anisotropies, making use of the recent Wilkinson Microwave Anisotropy Probe (WMAP) results WMAP . This result will then determine the contribution of the cosmic strings to the mass power spectrum. In this context, it is important to take into account the fact that the perturbations produced by the strings are isocurvature in nature (as emphasized e.g. in TTB ). This will effect the ratio of their contributions to the angular CMB power spectrum and the matter power spectrum. We then use the (approximately) known transfer function for strings and cold dark matter (which is larger on small scales than the standard transfer function for adiabatic fluctuations with cold dark matter) to estimate the total matter power spectrum on small scales, and analyze whether the loss of power on small scales when increasing the neutrino contribution to $`\mathrm{\Omega }`$ can be compensated by the increase in power coming from the presence of cosmic string seed perturbations.

## II Analysis

The recent WMAP data on the angular power spectrum of CMB anisotropies has mapped out with high precision the region of the power spectrum corresponding to the first acoustic peak. This peak is narrowly centered at a value $`l=220\pm 1`$ WMAP , in good agreement with a cosmology in which the Universe is spatially flat and the primordial spectrum of fluctuations is coherent and adiabatic. Cosmic strings, on the other hand, give rise to isocurvature fluctuations which are incoherent. As a consequence, there are no marked acoustic oscillations in the angular power spectrum PST ; Albrecht ; Battye ; Allen , but only a fairly broad Doppler peak Periv . Thus, the present data tightly constrain the maximal contribution of strings to the angular power spectrum of CMB anisotropies on large angular scales.

Since the theoretical predictions for the spectrum of CMB anisotropies resulting from cosmic strings in the acoustic peak region are quite uncertain, the bounds on the contribution of cosmic strings to the CMB angular power spectrum on large scales are also not certain. Different investigations have yielded bounds between $`1\%`$ Avelino  and $`10\%`$ Pogosian  (see also Endo ). We will denote the upper bound as $`f`$.

We will assume that the fraction $`f_{CDM}`$ of the dark matter is cold, and the rest hot. There are two sources of fluctuations, the conventional adiabatic perturbations (e.g. produced by quantum fluctuations during inflation) and the cosmic strings. We will assume that the cosmic strings are not correlated with the inflationary perturbations. In this case, the cross terms (between the string fluctuations and the adiabatic perturbations) to the angular CMB power spectrum and to the mass power spectrum vanish, and the mass power spectrum $`P(k)`$ calculated in linear cosmological perturbation theory can be written as

$$P(k)=(1f)P_a(k)+fP_{CS}(k),$$
(2.1)
where $`P_a`$ and $`P_{CS}`$ are the power spectra for a pure adiabatic model and a pure cosmic string model, respectively.

Since the tightest constraints on the fraction of hot dark matter come from the observed power spectrum of matter on the smallest scales for which linear perturbation theory is adequate and for which observational results are robust, we will in the following focus on these scales. For these scales, the contribution of hot dark matter to the matter power spectrum for adiabatic fluctuations is negligible, and thus

$$P_a(k)f_{CDM}^2P_0(k),$$
(2.2)
where $`P_0`$ is the power spectrum of the current concordance model (scale invariant adiabatic fluctuations in the $`\mathrm{\Lambda }`$CDM model 111We will consider the same background cosmological parameters $`\mathrm{\Lambda }`$, $`h`$, etc., as in this concordance model. As emphasized in Endo , the addition of cosmic strings to the model may yield a slightly different best fit value for the parameters. We will return to this issue in future work.). The result (2.2) is not exact for two reasons: first, since cold dark matter does cluster about the primordial perturbations, neutrino free streaming delays but does not totally prevent the clustering of neutrinos. This effect would lead to a spectrum larger than given in Eq. (2.2). On the other hand, the local clustering of the cold dark matter occurs as if the Universe were slightly open. This effect tends to reduce the spectrum from the result (2.2). However, both of these effects are small if the fraction of hot dark matter is not too large. Even if we use the exact transfer function of adiabatic fluctuations WDM , the final result does not change significantly.

On small scales, the transfer function for a model with cosmic strings and cold dark matter is different from the transfer function in the concordance model AS1 . The transfer function is, in fact, much larger. Thus, as noted in AS1 , a pure cosmic string model which is $`\sigma _8`$-normalized (i.e. normalized such that the power spectrum on a length scale of $`8h^1`$Mpc agrees with observations) produces too much small scale structure. Hot dark matter clustering onto cosmic string loops, on the other hand, is delayed but not prevented. Thus, the hot dark matter power spectrum in a cosmic string model on small scales is not exponentially suppressed. However, the hot dark matter cosmic string power spectrum on small scales is much smaller than the cold dark matter power spectrum, so that to a good approximation we can use

$$P_{CS}(k)f_{CDM}^2P_{(0,CS)},$$
(2.3)
where $`P_{(0,CS)}`$ denotes the power spectrum in a pure cosmic string model with cold dark matter and the current concordance parameters.

The question we would like to ask is whether the contribution of cosmic strings to the power spectrum, given by Eqs. (2.1) and (2.3) is large enough to compensate for the loss in small scale power from the adiabatic fluctuations \[given by Eq. (2.2)\] if $`f_{CDM}<1`$ and to render the theory compatible with observations. More specifically, we are interested in finding the smallest value of the cold dark matter fraction $`f_{CDM}`$ for which the theoretical power spectrum $`P(k)`$ exceeds the observational result $`P_{obs}(k)`$, i.e.

$$P(k)>P_{obs}(k)$$
(2.4)
on the smallest scales for which high quality observational results exist and for which we can trust the results of linearized perturbation analysis, namely $`k=0.2h\text{Mpc}^1`$ Hannestad ..

The results of AS2  indicate that on scales of $`k=0.2h\text{Mpc}^1`$ the cosmic string cold dark matter power spectrum $`P_{(0,CS)}`$ does not exceed the adiabatic cold dark matter spectrum $`P_0(k)`$ by more than a factor of 2. In fact, even on scales as low as $`1h^1`$Mpc, the difference is not more than a factor of 10. The results of AS2  are for models without a cosmological constant. However, as demonstrated in Wu , on the smaller scales of relevance here, the matter power spectrum does not depend on the value of the cosmological constant in a significant way, and thus we can make use of the results of AS2 . On the other hand, it is well known that the perturbations seeded by strings are isocurvature in nature (see e.g. TTB ). For isocurvature fluctuations, the CMB temperature anisotropies on large scales are given by

$$\frac{\delta T}{T}=\mathrm{\hspace{0.17em}2}\mathrm{\Phi },$$
(2.5)
where $`\mathrm{\Phi }`$ is the relativistic gravitational potential (see e.g. MFB ) which determines the magnitude of the mass perturbations. This is to be compared with

$$\frac{\delta T}{T}=\frac{1}{3}\mathrm{\Phi },$$
(2.6)
in the case of adiabatic fluctuations. Thus, if we normalize the power spectrum by the CMB anisotropies on large angular scales, the cosmic string power spectrum on large scales is a factor of 36 smaller than the power spectrum for adiabatic fluctuations. Hence, we see that on scales relevant to the large-scale structure, the decrease in the matter power spectrum which results from the addition of strings as a secondary source of cosmological fluctuations far exceeds the increase which is obtained by the change in the transfer function. Using the optimistic value 2 for the increase in the power due to the change in the transfer function on the relevant scale of $`k=0.2h\text{Mpc}^1`$, we conclude that in a model with both cosmic string and adiabatic scale-invariant fluctuations

$`P(k)`$ $`=`$ $`(1f)P_a(k)+fP_{CS}(k)`$ (2.7)
$``$ $`(1f)f_{CDM}^2P_0(k)+ff_{CDM}^2{\displaystyle \frac{2}{36}}P_0(k)`$
$`<`$ $`f_{CDM}^2P_0(k).`$
In a model with only adiabatic perturbations, the upper bound on the sum of the neutrino masses obtained by combining the large-scale structure and CMB data comes from demanding that

$$P(k)P_l(k),$$
(2.8)
where $`P_l(k)`$ is the observational lower bound on the mass spectrum. Since in this case

$$P(k)f_{CDM}^2P_0(k),$$
(2.9)
and knowing that the $`P_0(k)`$ fits the date well, this gives a lower bound on $`f_{CDM}`$ which is equivalent to an upper limit on the sum of the neutrino masses.

Combining these last three equations, we conclude that the inclusion of a contribution of cosmic strings to the primordial power spectrum does not lead to a relaxation of the cosmological bounds on the sum of the neutrino masses.

## III Discussion

In this paper we have shown that the current cosmological neutrino mass limits, which were derived under the assumption that the spectrum of cosmological perturbations is purely adiabatic, are robust against the addition of the maximal contribution which cosmic strings could make to the primordial power spectrum. A priori, this is not obvious, since cosmic string seeds survive neutrino free streaming and, therefore, clustering on small scales (scales which yield the tightest constraints on the neutrino masses) is not prevented but only delayed. However, if the power spectrum is normalized by the large-scale CMB anisotropies, then the contribution of strings to the matter power spectrum is suppressed by a factor of $`36`$ since the primordial string perturbations are isocurvature instead of adiabatic.

To obtain the above result, we have assumed that the string-induced fluctuations and the adiabatic perturbations are statistically uncorrelated. If they were correlated, a cross term in the power spectrum would appear which would only be suppressed by a factor of $`6`$. We have also assumed that bias is unimportant, i.e. that the observed galaxy correlation function equals the calculated matter power spectrum. However, in a cosmic string model we expect significant biasing (in particular in a model in which the dark matter has a significant component of hot matter). However, in the absence of a better understanding of the dynamics of cosmic strings (in particular whether string loops TuB  or string wakes SiVi ; Va ; StVe  are dominant), no firm conclusions can be drawn.

Even if it turns out that the addition of a cosmic string component to the primordial power spectrum does not change the cosmological neutrino bounds, this addition will be important for many cosmological issues, in particular for early structure formation Moessner  and reionization Avelino .

Acknowledgments

We are grateful to Guy Moore for interesting conversations at the beginning of this project. R.B. acknowledges the hospitality of the McGill Physics Department and the Fields Institute during the course of this work. At Brown, this work was supported in part by the U.S. Department of Energy under Contract DE-FG02-91ER40688, TASK A. M.Y. is partially supported by the Japanese Grant-in-Aid for Scientific Research from the Ministry of Culture, Sports, Science, and Technology. A.M. acknowledges support from CITA and NSERC.