text:
# Numerical Analysis of Jets Produced by Intense Laser

## 1 Introduction

Astrophysical jets found in supersonic bipolar flows from black holes, active galactic nuclei(AGN), protostars and similar phenomena are subject to active research. Several problems concerning jets are remaining to be solved, for example, why do the jets propagate such a long distance so stably. Many numerical simulations of astrophysical jets have been performed so far. 2D non-relativistic calculations have been done for AGN jets, by Norman et al. (1982). They found some dimensionless parameters relevant to the propagation of jets over a long distance. 2D and 3D relativistic calculations have been carried out by Martí et al. (1997); Aloy et al. (1999). Protostellar jets got more attention, detailed numerical simulations have been performed by several authors (see Reipurth, Bally & Devine (1997); Reipurth et al. (1997); Reipurth & Hearthcote (1997); Bally & Reipurth (2001) for observational data). Since the radiation cooling and/or external magnetic field are supposed to have very important effect on the plasma jets propagation, radiation hydrodynamic/magnetohydrodynamic codes are essential. For example, Blondin et al. (1990); Stone & Norman (1993) showed with their radiation hydrodynamic codes that the radiation cooling plays an important role for the collimation, leading to the formation of a dense cold shell at the head of the jet. Todo et al. (1993) studied, by means of with MHD code the propagation of jets from young stellar objects into ambient matter, and showed that the azimuthal magnetic field of $``$ 70 $`\mu `$G is strong enough to drive a helical kink instability. Since then, many efforts have been put to simulations in 2D and 3D with both radiation cooling and magnetic fields included (Cerqueira & Gouveia Dal Pino, 1997; Frank et al., 1998; Gardiner et al., 2000; Frank et al., 2000). The parametric researches that vary the form of the injection, steady or pulsed jets, and the configuration of magnetic field (helical or longitudal) helped to clarify the systematics of jets. In particular, it was shown that a helical magnetic field can affect the morphology of the head of jets. It was also demonstrated that Kelvin-Helmholtz instability is crucial for the morphology of jets, and linear as well as nonlinear stability analyses were carried out with or without magnetic field to characterize the influence of these factors on the jet shape (Xu et al., 2000; Hardee & Stone, 1997; Stone & Hardee, 2000; Ryu et al., 2000).

So far these studies are purely theoretical. Recently, however, with the advance of femtosecond laser technology, an experimental investigation of these astrophysical phenomena in the laboratory seems to be possible (Remington et al., 1999, 2000b). This new scientific field is called laboratory astrophysics or laser astrophysics. This progress enable the investigation of radiation hydrodynamics, hydrodynamic instabilities, and similar astrophysical phenomena in model experiments with intense lasers. Such experiments might be able to give us a new and improved insight into those phenomena which could not be obtained by studies of extraterrestrial observations or numerical simulations. Some experiments along these lines have already been carried out, in which hydrodynamic instabilities and the interactions of a strong shock with matter have been measured. Other experiments have been proposed (Remington et al., 2000a). Backup by numerical simulations are indispensable, because they can connect between the processes occurring in laser produced plasmas to astrophysical phenomena that are the aim of our studies. This paper is one of those attempts along this line. In this paper, we will propose a different experiment to understand astrophysical jets.

Some experiments relevant to astrophysical jets have been performed with intense lasers (Farley et al., 1999; Shigemori et al., 2000; Stone & Turner, 2000). The so-called ’cone’ target (Fig.1) is employed in those experiments. Another type of experiment has been recently done by Lebedev et al. (2001), using a conical array of fine metric wires. In general, intense laser matter interaction generates high temperature and high density plasma which starts to ablate immediately. If the intensity of laser is $`10^{14}\text{W cm}^2`$, the expansion velocity of the plasma reaches a few hundreds of kilometers per second. This is of the same order as the velocity of protostellar jets. This indicates that jets which realistically mimic astrophysical ones can be produced in laboratory. In these experiments, however, the plasma expands into vacuum, and one can not the investigate influence of ambient gas, which does exist for the astrophysical jets. As a matter of fact, the ratio of the mass density of the jet to that of the interstellar gas ($`\rho _{jet}/\rho _{ambient}`$) is $`110`$ for protostellar jets.

In this paper, we investigate the propagation of laser-produced jets not only in the vacuum but also in the ambient gas using the 2D hydrodynamics code developed recently (Mizuta, 2001). Although the code has two versions, relativistic and non-relativistic ones, in the present paper we use only the non-relativistic code.

This paper is organized as follows. In sections 2 and 3, the basic formulas, the numerical method and the initial conditions are presented. The main results are shown in section 4. We discuss some implications to astrophysical jets in section 5 and give a short summary in section 6.

## 2 Basic Equations and Numerical Method

Assuming a plasma having axial symmetry 2D Euler equation for perfect fluid is solved:

$`{\displaystyle \frac{𝒖}{t}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{(r𝒇)}{r}}+{\displaystyle \frac{𝒈}{z}}=𝒔,`$ (1)
where the conserved variable vector $`𝒖`$, the flux vectors $`𝒇`$ and $`𝒈`$, and the source vector $`𝒔`$ are defined as,

$`𝒖=(\rho ,\rho v_r,\rho v_z,E)^T,`$ (2)
$`𝒇=(\rho v_r,\rho v_r^2+p,\rho v_rv_z,(E+p)v_r)^T,`$ (3)
$`𝒈=(\rho v_z,\rho v_rv_z,\rho v_z^2+p,(E+p)v_z)^T,`$ (4)
$`𝒔=(0,p/r,0,0)^T,`$ (5)
$`E=\rho (ϵ+(v_r^2+v_z^2)/2).`$ (6)
In the above equations, $`\rho ,v_i,p,`$ and $`ϵ`$are mass density, $`i`$-component of velocity, pressure, and specific internal energy, respectively. We ignore viscosity and heat conduction of the fluid. These assumptions are justified as follows. As for the collisions, the parameter $`\delta =\lambda _{mfp}/L_{sys}`$ defined as the ratio of the collisional mean free path $`\lambda _{mfp}`$ of electrons or ions and, the characteristic length of the system $`L_{sys}`$, is $`10^410^6`$ for our models described below. Hence a hydrodynamic description of the jets is appropriate. Other dimensionless parameters characterizing the dissipative processes are the Reynolds number $`Rer_jv_j/\nu `$ which is a measure of the local viscosity and the Peclet number $`Per_jv_j/\nu _h`$ which describes the heat conduction inside the jet. Here $`r_j`$, $`v_j`$, $`\nu `$ and $`\nu _h`$ are the radius and velocity of the jet, the viscosity and thermal diffusivity coefficients respectively. Criteria for the influence of these dissipative phenomena on the flow are $`Re1`$ and $`Pe1`$. In our models we find $`Re10^9`$ and $`Pe10^5`$. Thus the neglection of these dissipative processes is well justified. These numbers are evaluated at $`t=2.0`$ ns at the center of the inflow region, see below.

An ideal gas equation is used in the computations,

$`p=\rho ϵ(\gamma 1),`$ (7)
where $`\gamma `$ is the adiabatic exponent (=5/3 in this paper).

We adopt the Marquina’s flux formula (Donat & Marquina, 1996) which is based on an approximate Riemann solver derived from the spectral decomposition of the Jacobian matrix of Euler equations. We have developed a computational code which solves Eq. (1) and (7) The code uses cylindrical coordinates with $`300(r)\times 1500(z)`$ grid points. The accuracy of the code is the second order in space due to a MUSCL method(van Leer, 1977, 1979) and is the first order in time.

It is assumed that the plasma is optically thin, so the radiation cooling effect is important. In the cooling term $`|J_{rad}|`$, only for the bremsstrahlung is included (Zel’dovich & Raizer, 1966),

$`|J_{rad}|=1.4210^{27}Z_{}^{}{}_{}{}^{2}T^{1/2}n_en_{ion}\text{erg cm}^3\text{s}^1,`$ (8)
where $`Z^{}`$,$`T`$,$`n_e`$ and $`n_{ion}`$ are the average ionic charge, temperature and number densities of electron and ion, respectively. The temperature and average ionic charge are related to the specific internal energy $`ϵ`$ by,

$`ϵ={\displaystyle \frac{3}{2}}{\displaystyle \frac{(Z^{}+1)T}{m_i}},`$ (9)
where $`m_i`$ is the mass of ion. The average ionic charge state $`Z^{}`$ and temperature $`T`$ were calculated iteratively from local instantaneous density and specific internal energy, as obtained from the hydrodynamic calculation by using Eq.(9) and approximate expression based on the Thomas Fermi model for $`Z^{}`$ (Salzmann, 1998). The formulas for the line emission processes were computed in the usual manner and are not included here for simplicity. As already mentioned the radiation transport is also neglected in this paper. Since we are interested mainly in the qualitative effects of radiation cooling, these approximations are adequate.

We have to distinguish the target matter from the ambient matter in order to make the radiation cooling effective only for the former. For this purpose we introduced another continuity equation for the target matter:

$`{\displaystyle \frac{(\rho f)}{t}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{(r\rho v_rf)}{r}}+{\displaystyle \frac{(\rho v_zf)}{z}}=0,`$ (10)
where $`f`$ is the fraction function. This equation is solved simultaneously with the hydrodynamic equations. With $`Z^{},T,`$ and $`f`$ thus obtained, the cooling term is calculated in a separate step. For numerical convenience, the radiation cooling term is turned off for very low temperatures, less than 70 eV. The radiation cooling term is neglected also in the vicinity of the target, where the radiative processes are very complicated and the basic picture of our model is not valid anyway.

## 3 Initial Condition

In order to simulate the experiments, we put a ’cone’ target at an end of the computational region. At $`t=0`$, the laser deposits all its energy on the surface. The depth of the surface is $`10\mu \text{m}`$. A Gaussian shaped laser pulse of duration of $`100\text{ps}`$ was assumed. This is much shorter than the dynamical time scale, which is more than a few nanoseconds. When the target made of gold is irradiated by intense laser of $`10^{14}\text{W cm}^2`$, the temperature of the target rises to a few hundred eV up to 1 keV, and the ablation takes place. Because the velocity of the ablated plasma is at most its sound velocity, $`v_{ab}c_s10^7\text{cm s}^1`$, the ablated plasma expands only up to $`10\mu \text{m}`$ during the laser irradiation $`100\text{ps}`$. This is much smaller than the scale of the system we consider in this paper. The laser energy deposited into thermal energy is $`E_l=526\text{J}`$ for all simulations. Energy of ionization is not included, therefore in a real experiment the incident laser energy, of course, has to be significantly larger.

Four simulations have been done in this paper. In cases 1 and 2, the ambient gas is very dilute ($`\rho _a=10^6\text{g cm}^3`$). In the first case, we neglect the radiation cooling, while in the second case the radiation cooling is turned on. These intend to simulate the experiments. In the other two cases, 3 and 4, we let the plasma expansion in the presence of an ambient gas. Case 3 neglects the radiation cooling, and case 4 is done with the radiation cooling term. No corresponding experiments have been carried out as yet. These models serve to clarify the effect of the ambient gas.

The temperature of the target surface, where the laser energy is deposited, increase to $`400`$ eV. The angle of the cone is $`126^0`$, which is similar to the experiments (see Fig.2). The target matter is chosen to be gold ($`{}_{}{}^{79}\text{Au}`$) so that we can obtain large effect of the radiation cooling due to a high-Z plasma. The density of the target $`\rho _t`$ is initially uniform, $`\rho _t=19.2\text{g cm}^3`$, and equals to the solid density of the gold. Table 1 summarizes the parameters of the initial conditions.

## 4 Results

### 4.1 Vacuum case

First, we discuss cases 1 and 2 which simulate plasma expansion into vacuum. This condition is similar to the experiments of Farley et al. (1999); Shigemori et al. (2000). For numerical reasons, we assumed in these cases as well a cold ($`T_a=0.026\text{eV}`$) and very low density ($`\rho _a=10^6\text{g cm}^3`$) ambient gas. Due to this very low density this ambient gas had no influence on any of our results. Figure 3 shows the density contours at different times with the radiation cooling off (Fig.3a) and on (Fig.3b), respectively.

In the case of no radiation cooling (Fig.3a), we can see two regions, that is, the inflow region in which matter flows to the axis, and the outflow region. This is schematically shown in Fig.4. Most of the laser-heated matter expands from the target, forming plasma with velocity of a few hundreds of kilometers per second, farther away from the target it converges to the symmetry axis, making a nozzle-like structure in the inflow region. Then, the plasma in the nozzle spouts out from the tip of the nozzle, forming the outflow region. Because the ambient matter is very dilute, there is no pressure support to sustain this nozzle structure. As a result, the outflow spreads out in all directions. This flow structure is suitable for low Z targets, because the effect of the radiation cooling is indeed negligible.

An analysis of the nozzle structure in the inflow region is in order. The converging flow from the cone target to the central axis is very similar with the structure considered by Cantó et al. (1988). From their theoretical model we can estimate the beam radius $`r_j`$ and velocity $`v_j`$ from the converging flow angle $`\theta `$, the reciprocal of the compression ratio $`\xi `$ between the density of the nozzle to that of the inflow and the velocity of the converging flow $`v_0`$ (see Fig.5). These are $`\theta =26.5^0`$ and $`\xi =1/4`$ in our case. Then $`r_j`$ and $`v_j`$ are determined as

$`r_j={\displaystyle \frac{\mathrm{tan}\alpha }{\mathrm{tan}\theta +\mathrm{tan}\alpha }}y_0,`$ (11)
$`v_j={\displaystyle \frac{\mathrm{cos}(\theta +\alpha )}{\mathrm{cos}\alpha }}v_0,`$ (12)
where $`y_0`$ is the width of the inflow and $`\alpha `$ is the angle between the conical shock and z-axis and given as

$`\mathrm{tan}\alpha ={\displaystyle \frac{(1\xi )+[(1\xi )^24\xi \mathrm{tan}^2\theta ]^{1/2}}{2\mathrm{tan}\theta }}.`$ (13)
Applying this theory to our case with $`y_0=250\mu \text{m}`$ and $`v_03\times 10^7\text{km s}^1`$, we obtain $`r_j=180\mu \text{m}`$ and $`v_j1\times 10^7\text{km s}^1`$, in good agreement with our numerical results.

For high Z targets, we have to consider radiation cooling effects. When radiation cooling is taken into account, the flow structure is changed dramatically (Fig.3b). The plasma in the inflow region is collimated strongly and shows a very thin nozzle at a later time (the radius of the structure is $`40\mu \text{m}`$). Although the collimation of the plasma is sustained in the inflow region by converging inflow from the target, the matter again spreads out in the outflow region, resulting in almost the same structure as that for the case without radiation cooling. These experiments are not suitable for the study of the propagation of collimated jets.

Figures 6 a, b and c are the density, pressure and temperature profiles along $`r`$-axis, respectively, at $`z=1000\mu \text{m},t=2.0\text{ns}`$. In the case with radiation cooling, the plasma is cooled efficiently around the symmetry axis, leading to the increase of the pressure gradient. As the density increases, the plasma is further cooled, and as a result a jet-like nozzle is formed, which was actually observed in the experiments. Altogether, our simulations support the experimental results that the radiation cooling is a crucial ingredient for the formation of this nozzle.

### 4.2 Dense gas case

In cases 3 and 4, the ambient matter density was much higher, $`\rho _a=10^3\text{g cm}^3`$. The time evolution of hydrodynamics leading to the formation of the inflow and outflow regions is essentially unchanged. Figure 7 shows the density contours at different times with the radiation cooling off (Fig.7a) and on (Fig.7b), respectively. The bow shock, which accelerates the ambient matter, can be seen clearly.

In the case with radiation cooling, the flow structure is dramatically changed as is obvious in Fig.7b. The effect of radiation cooling manifests itself in both the inflow and the outflow regions. In the inflow region, the nozzle structure becomes very thin just like in the vacuum case in the later time. The main difference appears in the outflow region. Unlike all the other cases, the width of the outflow region is very small. Opposed to intuition, this narrow structure is not a direct outcome of the thin nozzle in the inflow region. This is understood from the second panel of Fig.7b, in which there is already a very thin outflow region while there is no collimated structure in the inflow region. When the matter enters the outflow region, it tends to spread out as in vacuum case. However, the shocked ambient matter slows the expansion in the direction perpendicular to the symmetry axis, and the density and temperature of the flow remain high enough for the radiation cooling to be effective. As a result, a collimated plasma is produced in the outflow region.

This difference of structure in the outer region is remarkable, as is further evident in Fig.8, in which the density, pressure and Mach number profiles along $`r`$-axis, respectively, at $`z=2000\mu \text{m},t=4.0\text{ns}`$ are shown. The radius of jet is about 40$`\mu \text{m}`$, and the density ratio is $`\eta =\rho _{jet}/\rho _a100`$, namely, this is a so-called dense jet. From Fig.8b, we see that the narrow structure in the outflow region is sustained by the shocked ambient matter, whose pressure is about 1000 times larger than that of the unshocked ambient matter. The jet is supersonic $`M_{jet}70`$. The velocity of jet is $`400\text{km s}^1`$. This is of the same order as the velocity of protostellar jets. In order to study the propagation of astrophysical jets in the laboratory, we think it is necessary to do experiments with ambient gas.

## 5 Discussion

Firstly we consider the scaling law between a laboratory jet and a protostellar jet. The velocity of these two jets is almost the same-a few hundred kilometers per second. The ratio of the time scales is, however, widely different $`10^{24}\text{yr}`$ for protostellar jets and $`10^9\text{s}`$ for laboratory jets, a scale up of a factor of $`10^{1820}`$. For the experiment to simulate efficiently the astrophysical jet, the ratio of the lengths of jets has to be $`10^{1820}`$. The typical ratio of the radii of the jets ($`r_{prot}/r_{lab}`$) is $`10^{15}\text{cm}/10^3\text{cm}=10^{18}`$. Thus, the aspect ratios of the jets are close to each other. The temperature is a few eV for both jets, which is consistent with the scaling of velocity mentioned above. Thus, a good scaling law holds between the protstellar jet and the laboratory jet considered in this paper.

The relevant parameters which characterize the jet are the density ratio $`\eta `$, the Mach number and velocity of jets $`M_{jet},v_{jet}`$, the pressure ratio $`K`$ and the cooling parameter $`\chi `$. As the magnetic fields were not included our code, no comparison can be made with magnetohydrodynamic simulations (Todo et al., 1993; Cerqueira & Gouveia Dal Pino, 1997; Frank et al., 1998; Gardiner et al., 2000; Frank et al., 2000). As is obvious from the results obtained both for the vacuum and the ambient gas cases, the radiation cooling is the most crucial ingredient in the collimation of the jets.

The cooling parameter is defined from the cooling term which has already been employed in the hydrodynanmic simulations, $`|J_{rad}|`$ in Eq.(8). It is the ratio

$`\chi ={\displaystyle \frac{\text{cooling length}}{\text{jet radius}}}{\displaystyle \frac{v_{jet}\tau _{rad}}{r_{jet}}},`$ (14)
where $`v_{jet}`$ and $`r_{jet}`$ are the velocity and radius of jet, respectively, and $`\tau _{rad}`$ is a characteristic radiation cooling time. $`\tau _{rad}`$ is defined from the ratio of thermal energy density to emission power as

$`\tau _{rad}={\displaystyle \frac{e_{thr}}{|J_{rad}|}}.`$ (15)
The smaller $`\chi `$ is, the more effective the cooling is. The cooling term adopted in Eq.(8) is different from those used in other papers (Blondin et al., 1990; Stone & Norman, 1993). We assumed in this paper that the jet matter is Au, rather than hydrogen in most other papers. Our choice was influenced by the need to enhance the radiation cooling so that the hydrodynamics is affected within a few nanoseconds. In fact, the average ionic charge state $`Z^{}`$ is about $`2040`$ for $`T=50100\text{eV}`$. Although the detailed dependence on $`T`$ is different between our cooling term and others, it is emphasized that the dynamics should be similar as long as the cooling parameter $`\chi `$ has a similar value during the greater part of the time evolution.

In case 4 (with dense ambient gas and radiation cooling), the typical radiation cooling time is $`\tau _{rad}10\text{ns}`$ at $`t=1.0\text{ns}`$, and $`\tau _{rad}100\text{ns}`$ for $`t=4.0\text{ns}`$. These were evaluated at the points with the largest emission power. Accordingly, the radiation cooling parameters for each time are $`\chi 10`$ and $`\chi 100`$, respectively. Thus, the radiation cooling is effective for the earlier phase, while the jet is cool enough, while the radiation cooling is ineffective for the later phase. However, it is noted that we underestimate the radiation cooling, since we neglected the line emissions which could be substantial in our cases, and take into account only the bremsstrahlung. This was done mainly to keep numerical simplicity.

The aim of our calculations and future experiments is to understand better the physics of protostellar jets. Although the observational values of physical quantities are rather uncertain, they are typically, $`n_e10^3\text{cm}^3`$, $`T_j10^4\text{K}`$, $`v_j100\text{km s}^1`$ and $`r_j=10^{15}\text{cm}`$ (Bally & Reipurth, 2001). These values roughly correspond to $`\eta 110^3,M_{jet}1`$ with the main uncertainty coming from the lack of data on the density and temperature of the ambient matter. From these parameters, it is estimated that the total emission power is $`|J_{rad}|10^{18}\text{erg cm}^3\text{s}^1`$ for collisional by excited hydrogen atoms (Dalgarno & McRay, 1972). The emission power due to the bremsstrahlung is $`|J_{rad}|10^{19}\text{erg cm}^3\text{s}^1`$. Thus the bremsstrahlung is substantially weaker, and in general the bound-bound emission is dominant in such plasmas. The emission power $`|J_{rad}|10^{18}\text{erg cm}^3\text{s}^1`$ gives a cooling time of $`\tau _{rad}10^9\text{s}`$, which then leads to a cooling parameter of $`\chi 10`$. It is evident from these values that the experimental values used in this paper are appropriate for studying astrophysical jets. As mentioned in the introduction, Blondin et al. (1990); Stone & Norman (1993) studied the properties of propagating jets by varying $`\eta `$ between 1 and 10. Their parameters correspond to $`\chi `$ between 0.1 and 10. We think that our results have overlapping region to the initial conditions of those papers. This means that the experiments proposed here might give some new insights into the astrophysical jets studied theoretically in their papers.

## 6 Conclusion

In this paper results of four hydrodynamic simulations have been shown, in which jet were generated from ’cone’ target irradiated by intense laser. In all the simulations, we have found a common feature, that is, the existence of an inflow region and an outflow region. Depending on whether the radiation cooling is turned on or off, and whether there is some ambient gas or not, the flow structure inside these regions is very different from each other.

When the plasma flows into a very low density ambient gas, the collimation of the plasma jet occurs only due to radiation cooling. With radiation cooling, the collimated plasma shows up in the inflow region. However, since the structure is sustained by converging inflow from the target, it disappears in the outflow region, in which there is no pressure support to maintain it. This feature is very similar to that found in the experiments. We have also found that the nozzle structure with no radiation cooling is well described with the analytical model by Cantó et al. (1988)

In the case with ambient gas, a bow shock appears very clearly. The main difference between the vacuum and dense gas cases is that in the latter case the collimated flow appears not only in the inflow region but also in the outflow region. Particularly, the collimated plasma in the outflow region is sustained by the pressure of the shocked ambient matter.

We have shown the possibility to carry out a laboratory experiment, using high intensity laser-plasma interaction, to generate collimated plasma jets propagating into an ambient gas. The parameters of such plasma jet are similar to those of protostellar jets. If such an experiment is performed, we can get important information about the physics of astrophysical jets.

External magnetic fields were not included in our calculations. Experiments with magnetic fields have already proposed and will be carried out in the future. We intend to adjust our simulations to the needs of future experiments.

This work was carried out on NEC SX4, Cybermedia Center Osaka Univ. and HITACHI SR8000, Institute of Laser Engineering, Osaka Univ.

One of authors (A. Mizuta) would like to thank K. Sawada, H. Nagatomo and N. Ohnishi for good suggestions for numerical methods.


format: text

text:
# Flux flow noise and braided rivers of superconducting vortices

\[

## Abstract

Current-voltage measurements of type-II superconductors are described by a coarse-grained model of superconducting vortex dynamics. We find that the power spectra of the voltage fluctuations, and the noise power, are related to the large scale morphology of the plastic flux flow. At currents corresponding to the peak in differential resistance, the flux flow forms a braided river, the noise power is maximal, and the power spectra has a $`1/f^\alpha `$ form with $`\alpha 1.8`$ over a wide frequency range. This agrees with recent experiments on $`NbSe_2`$. The observed variation of $`\alpha `$ with applied current is a crossover phenomenon.

\]

Temporally correlated fluctuations with a power spectrum characterized by a frequency distribution of the form $`1/f^\alpha `$ appear ubiquitously in nature. These correlations are often observed in time series of nonequilibrium systems such as fluvial rivers , nerve cells , quasars , and Barkhausen emission , to name a few. Despite attempts at a general explanation , the concrete sources of these fractal temporal fluctuations are usually unknown. Here we show, using a cellular model of vortex dynamics, that the $`1/f^\alpha `$ spectrum observed in voltage noise measurements of type II superconductors is related to the complex spatio-temporal pattern of the plasticly flowing magnetic flux through the system. This suggests a general link between long range spatial correlations observed in driven dissipative systems and long range temporal correlations obtained in time series measurements of them. This link was first proposed in the context of sandpile models , but the original models did not display nontrivial $`1/f^\alpha `$ noise . The sandpile type model presented here to describe plastic flux flow, however, clearly demonstrates this link.

Marley, Higgins, and Bhattacharya (MHB) found in current-voltage $`(IV)`$ measurements on the type-II superconductor $`NbSe_2`$ that the voltage fluctuations show a $`1/f^\alpha `$ decay in the power spectrum with an exponent $`\alpha `$ that varied with the applied transport current. Holding the applied magnetic field constant and increasing the transport current, they found that $`\alpha `$ varied, increasing from its initial value near 1 to about 1.8, and then falling again. The low frequency noise power $`fS(f)`$ was also measured. At a constant applied field, $`fS(f)`$ has a pronounced maximum at the transport current for which the largest value of $`\alpha `$ was measured, suggesting that both phenomena are related.

Those experiments were performed in the so-called “peak-effect regime” where the soft flux line lattice formed by the superconducting vortices is thought to deform by tearing, giving plastic flow . At magnetic fields $`H_{pl}<H<H_p`$, the critical current, $`I_c`$, is an increasing function of magnetic field, and the $`IV`$ relation is highly nonlinear. As $`H`$ is increased, the size of the peak in the differential resistance, $`dI/dV`$, grows, reaching a maximum value at $`H_c1.88T`$ which is significantly different from $`H_{pl}1.81T`$ and $`H_p1.95T`$. The magnitude of the differential resistance peak shrinks and it occurs at increasingly larger $`I`$, with increasing $`H>H_c`$. At $`H_c`$ itself, the peak in the differential resistance occurs at $`I_c`$, which marks the onset of the flow of superconducting vortices where the voltage, $`V`$, becomes non-zero.

In this letter, we describe the $`IV`$ noise experiments of MHB using a coarse-grained cellular model of superconducting vortex dynamics , reminiscent of a sandpile model. We show that a fractal power spectrum exists only at the maximum of the differential resistance while for other values of the applied transport current, $`I`$, significant deviations from the power-law behavior occur. This can be easily understood in terms of the underlying morphology of the vortex flow, which has been shown to form a braided river at the current corresponding to the peak in the differential resistance . The braided structure is also responsible for the maximum of the low frequency noise power. In particular, our results support claims  that complex patterns of the vortex flow (which were not experimentally observed themselves) are responsible for the observed noise in the $`IV`$ measurements.

The cellular vortex model  used here is defined as follows. Consider a square lattice of size $`L\times L`$ in which each site $`i`$ is occupied by an integer number of vortices, $`m_i`$. The total energy of the vortex system includes the short-ranged repulsive pairwise interaction between the vortices and the attractive interaction of vortices with the pinning potential, $`\widehat{V}`$. For a given configuration of vortex number $`\{m_i\}`$, the total energy of the system is:

$$H\left(\{m_i\}\right)=\underset{i,j}{}J_{ij}m_im_j\underset{i}{}\widehat{V}_im_i,$$
(1)
where $`J_{ij}`$ includes an on-site interaction, and a weaker nearest-neighbor interaction. The model describes a vortex system coarse-grained to the scale of the London length.

The motion of the vortices is assumed to be highly overdamped, described by the equation of motion $`\stackrel{}{v}\stackrel{}{f}`$. The force to move a unit vortex from $`x`$ to $`y`$ is calculated by taking a discrete gradient of Eqn. (1),

$`F_{xy}`$ $`=`$ $`\widehat{V}_y\widehat{V}_x+[m_xm_y1]`$ (3)
$`+r[m_{x1}+m_{x2}+m_{x3}m_{y1}m_{y2}m_{y3}],`$
where the nearest neighbors of $`x`$ are $`y`$, $`x1`$, $`x2`$, and $`x3`$, and the nearest neighbors cells of $`y`$ are $`x`$, $`y1`$, $`y2`$, and $`y3`$. The onsite vortex interactions have unit strength, while $`r`$ describes the strength of the nearest neighbor vortex interactions. The normalized pinning potential $`\widehat{V}_x`$ is a random number taken from a uniform distribution in the interval between zero and $`V_{max}`$. In each time step, all cells of the lattice are updated in parallel. A single vortex moves from a cell to a neighboring cell if the force in that direction is positive, or equivalently if the total energy of the system is lowered.

Many alternatives exist to handle the situation when more than one unstable direction appears for a vortex to move. In this work, we use two of them. In the first version of the model, the unstable direction that has the largest force is chosen and the vortex moves in that direction. The fluctuation spectrum in this case is generated solely by a deterministic dynamics. The second version of the model is stochastic. If more than one unstable direction exists, one of them is chosen randomly and the vortex moves in that direction.

Since the parameter $`r`$ controls the strength of the vortex repulsions, increasing $`r`$ decreases the steady state slope of the vortex pile. Thus, $`r`$ controls the critical current, which in the simulations is the slope of the vortex pile where vortices begin to continuously move. In the experiments, the critical current can be controlled by the applied magnetic field. Therefore, increasing the applied field $`H_{pl}<H<H_p`$ in the experiments is represented by decreasing $`r`$ in the simulations. Periodic boundary conditions apply at the top and bottom boundaries.

Initially, the system was prepared with a V-shaped flux density profile characteristic of the Bean state. This describes the application of a magnetic field, and was accomplished by starting with an empty lattice and then raising the left and right boundaries to equal heights, letting vortices enter and penetrate into the system until a stable state was reached. Then, a transport current $`I`$ was applied. According to Ampere’s law, a transport current will induce a magnetic field that raises the magnetic field on one boundary and lowers it on the opposite side. Thus, a transport current is modeled by shifting the height of, say, the left boundary verses the right boundary. Those changes were made in half integer steps. At small $`I`$, the shift between the left and right boundary conditions is small, and the system reaches a pinned state with no vortex flow. However, at larger $`I`$, the average slope is steep enough that vortices continue to flow from the high (left) boundary to the low (right) boundary, where they are removed from the system. The $`IV`$ characteristic then is the relation between the magnitude of the shift (representing the applied transport current) and the average flow rate of vortices (representing the voltage) when the critical current is exceeded.

The differential resistance ($`dI/dV`$) curves for a number of different values of $`r`$ are shown in Fig. 1  for the deterministic model. For $`r=0.13`$, the $`IV`$ curve is entirely convex, so the corresponding differential resistance curve has a large peak at the onset of vortex motion. This value of $`r`$ thus corresponds to $`H_c`$ in the experiments of MHB. As $`r`$ decreases, the $`IV`$ curve becomes S-shaped, and the peak in the differential resistance gets smaller and occurs at increasing $`I`$. These results are entirely consistent with the experimental $`IV`$ measurements in the range $`H_c<H<H_p`$. As shown in , the shape of the $`IV`$ curve is related to the large scale morphology of the vortex flow. In particular, the peak of the differential resistance, or equivalently the inflection point of the $`IV`$ curve, occurs where the vortex flow has a braided river morphology which is characterized by a self-affine (multi)fractal structure. At small $`r`$ the curve is quite noisy, especially near onset. This phenomena has been seen in experiments where it has been called “IV fingerprints”. As will be discussed below, this noisiness is due to glassy channel switching of the vortex flow.

The power spectra for a few values of the applied current at $`r=0.13`$ are shown in Fig. 2a. The current $`I=0.185`$ is just above onset, where for this value of $`r`$ the vortex flow is a braided river. The power spectrum for this value of $`I`$ shows $`1/f^\alpha `$ type behavior with $`\alpha 1.8`$ over a range of frequencies from about $`5.0\times 10^4`$ to about $`5.0\times 10^2`$. At higher frequencies, there is a crossover to behavior that decays less rapidly with $`f`$, presumably due to effects associated with the discrete lattice used in the simulations. At lower frequencies, there is a crossover to white noise which is related to the longest time scale over which correlations are retained in the system.

As $`I`$ increases, the range of frequency where $`1/f^\alpha `$ noise occurs shrinks and the white noise roll off at low frequencies sets in at a higher frequency. This indicates that the long range correlations in the system fade away. That happens presumably because, as the braided river floods, which occurs as $`I`$ increases, the spatial correlation length shrinks. The frequency at which the high frequency crossover occurs is independent of $`I`$, but the magnitude of the power spectrum at high frequency grows with $`I`$ because of the larger number of vortices moving. If one were to fit a value for the exponent $`\alpha `$ as a function of $`I`$, a varying exponent would be measured due to crossover. The exponent $`\alpha `$ would range from about 1.8 near the peak in the differential resistance to near 0 at large $`I`$. This range of $`\alpha `$ is precisely what was reported by MHB close to $`H_c`$. Agreement is also found for the low frequency noise power which is obviously largest for $`I=0.185`$, and decreases with increasing $`I`$. In particular, the crossover picture is consistent with the very limited range ($`1`$ decade) of power-law behavior found in experiments . Note that for $`r=0.13`$, both the deterministic and stochastic versions of the model give essentially the same results.

Mostly similar results are also found for the power spectra at $`r=0.01`$, which, for the deterministic version of the model, are shown in Fig. 2b. However, there are also some differences. At this value of $`r`$, the peak in the differential resistance does not occur at onset, which is at $`I_c=0.310`$, but rather at $`I=0.375`$. Near $`I=0.375`$, the power spectrum has a $`1/f^\alpha `$ form with $`\alpha 1.8`$ over a broad range of frequencies. There are again crossovers at both large and small frequencies to regions where the power spectrum decays less rapidly with $`f`$. The frequency at which the low frequency crossover occurs gets larger with increasing $`I`$, as the braided river floods and the spatial correlation length decreases.

The low frequency noise power as a function of $`I`$ for $`r=0.01`$ is shown in Fig. 3. It is calculated by adding all of the points of $`S(f)`$, with $`f`$ between $`10^5`$ and $`10^4`$. The result rises many orders of magnitude and peaks where the river flow is braided, at the peak in the differential resistance. This result is also seen in the experiments.

At transport currents $`I`$ closer to $`I_c`$, where the S-shaped IV curve is concave, and the vortices flow through isolated filamentary channels , the voltage fluctuations are periodic in the deterministic model. Also, in this regime, long-time scale transient and glassy behavior occurs. This is due to channel switching. At $`I`$ close to $`I_c`$, the system typically switches between a number of different sets of channels before settling into the steady-state. The system can continue switching for $`10^6`$ timesteps, or more, before settling. It does finally settle into a stationary configuration, but only after many different sets of channels are probed. As the current increases, the transient time becomes longer than we can measure, possibly infinite, near the peak in the differential resistance.

For the stochastic model, the results at $`r=0.01`$ differ somewhat from those of the deterministic model. First, the critical current $`I_c`$ is larger in the stochastic case. Apparently, the stochastic lattice updating can frustrate flow through isolated filamentary channels. Second, in the stochastic model for $`I>I_c`$ the system never settles into a stationary river configuration. Instead, channel switching continues to occur for at least $`10^7`$ lattice updates, and presumably continues indefinitely. However, $`1/f^\alpha `$ noise with $`\alpha 1.8`$ also occurs in this version of the model near the peak in the differential resistance.

Additionally, it is interesting to note that, in the braided regime, the change in the total number of vortices contained within the system also show $`1/f^\alpha `$ fluctuations with $`\alpha 1.15`$. This agrees quantitatively with the fall-off power spectrum measured in the field ramping experiments of Field and co-workers .

Other theoretical studies of vortex dynamics have also found $`1/f^\alpha `$ behavior in the power spectrum of the voltage noise. Driving the system in a different fashion designed to maintain a critical state, Mohler and Stroud  using the cellular model  and Olson, et al. , using molecular dynamics simulations found fractal behavior in the noise spectrum. Those results are consistent with the fact that in both cases the vortex flow has a braided river morphology . However, the value of $`\alpha `$ they measured deviates from that measured here due to the different driving. Another approach was taken by Domínguez . Integrating a set of equations that results from a discretization of the London model, he found two regimes of vortex flow, filamentary motion near onset, and homogeneous motion at higher currents. In the filamentary channel regime, he discovered $`1/f`$ behavior of the voltage noise power spectrum. Using molecular dynamics simulations, Kolton et al. , also found the same two regimes and a peak in the low frequency noise power at the transition between them. They also suggested that the transition is from plastic flow in the filamentary regime to smectic flow in the homogeneous regime, thereby associating the transition with a reordering of the flux line lattice. The cellular model used here cannot describe any possible structural transitions in the flux line lattice. Nevertheless, we find remarkable quantitative similarity with experimental voltage noise measurements involving the transition from filamentary, through braided, to homogeneous motion of the vortices.

In summary, voltage noise in $`IV`$ measurements has been studied in a simple, cellular model of vortex dynamics. The model captures the essential features of the plastic flow of the flux line lattice formed by the superconducting vortices: $`1/f^\alpha `$ behavior of the power spectrum was observed where the differential resistance has a peak, and the vortex flow has been shown to form a braided river. This result is robust, occurring over all values of the vortex interaction strength $`r`$, and for two different versions of the model, one completely deterministic and the other stochastic. The results model experiments conducted on $`NbSe_2`$, in the peak-regime where the FLL moves plasticly. Most notably, quantitative agreement is found in the behavior of the exponent $`\alpha `$ between simulations and experiments, both of which find $`\alpha 1.8`$ at the peak of the differential resistance. Both also find that the low frequency noise power is maximal at the peak in the differential resistance. Therefore, we conclude that is the self-organized spatio-temporal structure of the braided river flow that is responsible for the observed $`1/f^\alpha `$ behavior of the $`IV`$ power spectra in superconductors.

This work was supported by the NSF through grant #DMR-0074613. KEB also thanks the Alfred P. Sloan Foundation and the Texas Center for Superconductivity for support. JD would like to thank the DAAD for financial support and the Imperial College for support and hospitality.


text:
# Nonlinear Analysis of the Eckhaus Instability: Modulated Amplitude Waves and Phase Chaos with Non-zero Average Phase Gradient

## 1 Introduction

The emergence of chaotic behaviour from ordered states in spatially extended systems has been the subject of many recent experimental and theoretical investigations . Nonetheless, the mechanisms leading from stationary regimes to chaotic (or spatially irregular) phases still pose many challenging questions. One of the most studied instabilities in extended oscillatory systems is the Eckhaus instability of plane waves .

The occurrence of this instability has been experimentally observed in many quasi one-dimensional systems like the oscillatory instability of a Rayleigh-Bénard convection pattern , hydrothermal waves , heated wire convection , sidewall convection , the Taylor-Dean system  and internal waves excited by the Marangoni effect . The Eckhaus instability also plays an important role in the radial dynamics of spiral waves in the Belousov-Zhabotinsky reaction .

The complex Ginzburg-Landau equation (CGLE)  is the appropriate amplitude equation to describe the slow dynamics near a supercritical transition to unidirectional traveling waves. In one spatial dimension, the CGLE reads:

$$_tA=A+(1+ic_1)_x^2A(1ic_3)|A|^2A,$$
(1)
where $`c_1`$ and $`c_3`$ are real coefficients and the field $`A=A(x,t)=|A(x,t)|\mathrm{e}^{i\phi (x,t)}`$ has complex values. As exact solutions the CGLE admits plane waves of the form $`A_q(x,t)=a_q\text{e}^{i(qx\omega _qt)}`$, where $`q`$ indicates the wavenumber, $`a_q=\sqrt{1q^2}`$ and $`\omega _q=c_3+q^2(c_1+c_3)`$.

A linear stability analysis  of these solutions can be performed by considering the perturbed solution $`\stackrel{~}{A}_q(x,t)=(a_q+\delta a)\text{e}^{i(qx\omega _qt)}`$, where $`\delta a\text{e}^{ikx}\text{e}^{\sigma (k)t}`$. The growth rates associated to the complex perturbation $`\delta a`$ is

$`\sigma (k)`$ $`=`$ $`k^22iqc_1k(1q^2)`$ (2)
$`\pm \sqrt{(1+c_3^2)(1q^2)^2[c_1k^22iqkc_3(1q^2)]^2}.`$
The plane waves become linearly unstable to long wavelength perturbations ($`k0`$) for $`q=q_E\sqrt{(1c_1c_3)/(2(1+c_3^2)+1c_1c_3)}`$. This limit is called Eckhaus instability . Above the Benjamin-Feir-Newell (BFN) line $`1c_1c_3=0`$, all plane waves are unstable to homogeneous perturbations. For a given $`qq_E`$, the corresponding plane wave is linearly unstable against perturbations with wavenumbers $`k`$ inside the interval $`0<|k|<k_c`$. $`k_c`$ increases for increasing values of the parameters $`c_1`$ and $`c_3`$ .

As noticed in , the Eckhaus instability is a convective instability. Thus, it is relevant in the systems with periodic boundary conditions considered here, while it would be suppressed in fixed boundary conditions , e. g. zero-flux or Dirichlet boundary conditions. In the latter geometries, the absolute instability of the plane waves has to be computed. It occurs for sufficiently large $`q`$ or/and $`c_1`$,$`c_3`$ values inside the Eckhaus unstable range.

Another interesting aspect of the Eckhaus instability is found when the nonlinearities of the CGLE are taken into account. A weakly nonlinear analysis  revealed that for

$$c_1^2(16c_3^2)+c_1(2c_3^3+16c_3)(8+c_3^2)>0$$
(3)
the Eckhaus instability becomes supercritical, i.e. the instabilities are saturated and the emerging quasi-periodic solutions (resp. modulated amplitude waves) coexist with the unstable plane waves.

Numerical simulations  provided examples of such stable modulated amplitude waves (MAWs). Stable MAWs have also been observed in experiments on surface-tension-driven hydrothermal waves  as well as on the Taylor-Dean system  and on internal waves excited by the Marangoni effect.

In addition, the CGLE exhibits two qualitatively different spatiotemporal chaotic states known as phase chaos (when the modulus of the field $`|A|`$ is bounded away from zero) and defect chaos (when the phase of $`A`$ displays singularities where $`|A|=0`$) . The subclass of MAWs with zero average phase gradient is important for understanding the transition from phase to defect chaos (see full curve in Fig. 1) . In the phase chaos regime states with nonzero average phase gradient $`\nu `$ have a dynamics quite different from that at $`\nu 0`$. In particular, these states can be either chaotic or regular depending on the initial conditions and on the parameters $`c_1,c_3`$ and $`\nu `$. In this paper we will focus on MAWs with $`\nu 0`$ and on the dynamical regime associated to them, that is referred to as “wound-up” phase chaos . It will be shown that MAWs and wound-up phase chaos exist between the dashed and the full curve in Fig. 1.

In Section 2 the analysis of MAWs introduced in  is extended to arbitrary values of the average phase gradients of the field. The two parameter family of MAW solutions is parametrized by the spatial period $`P`$ of the modulation and by the average phase gradient

$$\nu :=\frac{1}{P}_0^P𝑑x\phi _x.$$
(4)
For plane wave solutions, $`\nu `$ equals the wavenumber $`q`$. In analogy, the phase gradient $`\phi _x`$ is often called “local wavenumber”. A linear stability analysis will show that MAWs with $`\nu 0`$ can be stable even in infinitely large systems. In contrast, MAWs with $`\nu =0`$ are always unstable in systems of large length $`LP`$. For systems with periodic boundary conditions the average phase gradient of the whole system can only be changed, if a space-time defect occurs : $`|A(x,t)|`$ drops to zero and $`\phi _x`$ locally diverges at a defect. Persistent phase chaos with conserved $`\nu \nu _M0`$ has been observed in numerical simulations of the CGLE (1) . The maximum conserved average phase gradient $`\nu _M`$ decreases as function of the coefficients $`c_1,c_3`$  and vanishes at the apparent transition from phase to defect chaos. $`\nu _M`$ was therefore suggested  as an order parameter for this transition. We extended the numerical determination of $`\nu _M`$ towards smaller $`c_3`$ and report the corresponding data in Fig. 2.

In Section 3, a nonlinear analysis of the Eckhaus instability allows estimates for which parameter values defects occur. Lower and upper bounds for the limit $`\nu _M`$ of wound-up phase chaos are derived from the existence and stability properties of the MAWs. For increasing $`\nu `$ the degree of chaoticity associated with the wound-up phase chaos decreases . For large enough $`\nu `$ the dynamics can even become regular and stable quasi-periodic MAWs appear. The diamonds in Fig. 2 indicate this stability limit for numerical simulations with fixed $`c_1=3.5`$. The analysis in Subsection 2.3 will clarify this observation.

The large number of parameters ($`c_1,c_3,\nu ,P`$) calls for restrictions. We limit our analysis to fixed $`c_1=3.5`$ since most previous numerical work has been done at this value . The results will be presented as projections of the $`P`$ direction onto the ($`c_3,\nu `$) plane as well as in cuts through the parameter space spanned by $`c_3,\nu `$ and $`P`$. Additional investigations of the existence domains of MAWs revealed qualitatively similar results for fixed $`c_1=0.4,1.2,2.1`$ and $`5`$ and variable $`c_3`$ as well as for fixed $`c_3=0.83`$ and variable $`c_1`$. Two of these choices were studied by numerical simulations in . A similarity transformation maps coherent structures onto each other along curves $`(c_1+c_3)/(1c_1c_3)=const`$ in coefficient space . The parameters $`\nu ,P,\omega ,v`$ are transformed accordingly. One can thereby extend the results presented here to other values of the coefficient $`c_1`$. Section 4 discusses possible observations of MAWs in experimental systems. Finally, Section 5 summarizes the main results.

## 2 Existence and stability of modulated amplitude waves

### 2.1 Coherent structure approach

In analogy with the linear analysis of the Eckhaus instability of plane waves we make the following ansatz for saturated modulations

$$A(x,t)=a(z)\text{e}^{i\stackrel{~}{\varphi }\left(z\right)}\text{e}^{i\left(qx\stackrel{~}{\omega }t\right)}$$
(5)
and rewrite it as

$$A(x,t)=a(z)\text{e}^{i\varphi \left(z\right)}\text{e}^{i\omega t},$$
(6)
where $`a`$ and $`\varphi `$ are real-valued functions of $`z:=xvt`$ and $`\varphi (z)=\stackrel{~}{\varphi }(z)+qz,\omega =qv\stackrel{~}{\omega }`$. Here $`a(z)`$ and $`\varphi (z)`$ represent coherent structures . Coherent structures have been studied extensively  and play an important role in various regimes of the CGLE .

Substitution of ansatz (6) into the CGLE (1) yields the set of three coupled nonlinear ordinary differential equations (ODEs)

$`a_z`$ $`=`$ $`b`$
$`b_z`$ $`=`$ $`\psi ^2a\gamma ^1[(1+c_1\omega )a+v(b+c_1\psi a)(1c_1c_3)a^3]`$ (7)
$`\psi _z`$ $`=`$ $`2b\psi /a+\gamma ^1[c_1\omega +v(c_1b/a\psi )(c_1+c_3)a^2]`$
where $`b:=a_z,\psi :=\varphi _z`$ and $`\gamma :=1+c_1^2`$ . The continuation software AUTO97  is used to compute the periodic orbits of the ODEs (7) that correspond to spatially periodic functions $`a(z),\varphi (z)`$. In order to choose a unique solution from the continuous two-parameter family of periodic orbits we set the system size $`L`$ equal to the period $`P`$ of the periodic orbit and fix its average phase gradient by $`\nu =\frac{1}{L}_0^L\psi 𝑑z`$ .

The continuation procedure starts from a fixed point $`(a,b,\psi )=(\sqrt{1q^2},0,q)`$ that corresponds to a plane wave solution. Varying $`c_3`$, a Hopf (HB) bifurcation (filled square in Fig. 3) is detected in the ODEs where the mode with the smallest possible wave number $`k_{HB}=2\pi /P`$ destabilizes the plane wave. Continuing the resulting branch of MAWs the free parameters $`\omega `$ and $`v`$ are adjusted by the continuation algorithm. The continuation follows a unique branch of MAWs with $`\nu =q`$ and $`P=L`$. Fig. 3 shows examples of resulting bifurcation diagrams.

### 2.2 Existence limits of MAWs

Upon increasing of $`c_3`$ amplitude modulations grow and develop a localized depression $`|A|_{min}`$ where $`\varphi _x`$ has a maximum (see Fig. 3b,c). As for $`\nu =0`$, these MAWs are called the lower branch in contrast to the coexisting upper branch MAWs. The upper branch MAWs are always unstable, while the MAWs of the lower branch can be stable in appropriate parameter regions. Examples of these lower branch MAWs have been obtained by numerical simulations earlier . They have been analyzed in detail in . Numerical simulations can neither uncover unstable upper branch MAWs nor elucidate the limits of existence of MAWs. The bifurcation analysis presented here reveals that upper and lower branch meet and terminate in a saddle-node (SN) bifurcation (filled triangle in Fig. 3 a,d,e ). Due to the SN bifurcation the upper branch MAWs always have at least one unstable eigenmode, see also . The upper branch continues to negative $`c_3`$ and there connects to another instability of the plane wave with identical $`\nu `$ and $`P`$. In the following we will concentrate on the lower branch MAWs.

For large $`\nu `$ and small $`P`$, the Hopf bifurcation is no longer supercritical and an unstable branch emerges directly from the plane wave. This is in agreement with analytical predictions . Fig. 4 shows an example which also includes the second HB at negative $`c_3`$. For $`\nu =0`$ the MAWs emerge stationary  and acquire $`v0`$ above a subsequent drift pitchfork bifurcation . In the present case $`\nu 0`$ the plane wave already breaks the reflection symmetry, the initial MAW has a nonzero velocity and the drift pitchfork (DP) bifurcation (filled diamond) is unfolded. See Fig. 5 for an example at fixed $`c_3=2`$. The branch emerging at the HB in Fig. 5b represents the MAWs as discussed above. The second branch in Fig. 5b emerges at the period doubling (PD) bifurcation (open square) of MAWs with half the period. It always has unstable eigenmodes that drive the dynamics away from it to the coexisting MAWs of shorter period. Therefore this (upper) branch plays no essential role and is not treated further.

#### 2.2.1 Infinite system size

We have analyzed the existence of lower branch MAWs in the entire parameter space ($`c_3,\nu ,P`$) at fixed $`c_1=3.5`$. The system size is assumed infinitely large in order to allow for arbitrary periods $`P`$ of MAWs. Fig. 6 shows examples of existence domains for $`P=15,P=30`$ and $`P\mathrm{}`$. We find that both HB and SN shift to larger $`c_3`$ as the period $`P`$ is decreased. The same behavior has already been observed in the special case $`\nu =0`$ . The dotted curve in Fig. 6 indicates the “envelope” of all SN bifurcations for MAWs of any period and therefore is the upper boundary for the existence domain of the MAWs.

#### 2.2.2 Medium system size

Experimental setups and numerical simulations are restricted to finite system size $`L`$. Often periodic boundary conditions (corresponding to an annular geometry) are used in order to study bulk effects of extended systems and to minimize boundary effects. The periodic boundary conditions also restrict possible modes of perturbations. As described by Eq. (2) the instability threshold of plane wave solutions depends on the wavenumber $`k`$ of the perturbation. Since in the studied range of coefficients the Eckhaus instability is a long-wavelength instability the plane waves will be stabilized in small systems. The instability threshold is shifted to larger values of the coefficients $`c_1,c_3`$ and can be computed from Eq. (2) setting $`k=2\pi /L`$.

Clearly the selection of perturbations by periodic boundary conditions also restricts possible MAWs. Their average phase gradient $`\nu `$ and the period $`P`$ have to be consistent with the system size and this renders the two-parameter family of MAWs discrete. It is thus convenient to parametrize MAWs by the average phase gradient $`\nu `$ and the ratio $`n`$ of wavelength

$$n:=\frac{P}{2\pi /\nu }.$$
(8)
The ratio $`n`$ takes values of integer fractions where the nominator counts the number of underlying wave length and the denominator the number of humps of the modulation. Hence this quantity is easily accessible in experiments. The existence domains of MAWs with respective $`n`$ are presented in the ($`c_3,\nu `$) parameter plane in Fig. 7.

#### 2.2.3 Small system size

Here we focus on the extreme case. The shortest possible system with periodic boundary conditions only contains one wavelength of the plane wave, consequently its length $`L`$ is given by $`L=2\pi /\nu `$. In  the quantity $`\nu _U`$ was determined in analogy to $`\nu _M`$ for large systems. $`\nu _U`$ denotes the largest $`\nu `$ for which none of the random initial conditions (different realizations of noise added to a plane wave) produced a defect. In the following these data (symbols in Fig. 8) are compared to the existence domains of MAWs.

Within the light shaded area in Fig. 8 plane wave solutions with wavenumber $`\nu `$ are stable in the short system. The stability area extends over the phase chaos and into the defect chaos region. We stress this result because from the experimental observation of stable plane waves one cannot necessarily infer that the dynamics of the system may be reproduced by the CGLE with coefficients $`c_1,c_3`$ in the Benjamin-Feir stable region. The dashed curve denotes a subcritical instability. Only unstable upper branch MAWs exist to the left of this curve. For smaller $`\nu `$ the instability again turns supercritical and stable lower branch MAWs exist inside the dark shaded region. The thick full curve gives the SN bifurcation for MAWs with $`P=L`$. The thin curves show the respective limits of MAWs with shorter period. Defects are expected beyond the SN  and the subcritical instability which well reproduces the data from numerical simulations  except at small $`\nu `$. Simulations with $`\nu \nu _U`$ resulted in modulations with a single hump (squares) or with two (triangles) or three (stars) humps of different size . The latter two are observed above the SN of MAWs with $`P=L`$. Here the initial conditions select MAWs with shorter period which only coexist at small $`\nu `$. The SN with $`P=L`$ nevertheless gives a lower bound for the formation of defects.

### 2.3 Instabilities of MAWs

In contrast to the case of MAWs with zero average phase gradients some MAWs with non-zero $`\nu `$ are stable even in very large systems. A linear stability analysis of MAWs as in  yields the spectrum of eigenvalues as shown in Fig. 9 for a typical example. From Fig. 9 we conclude that for this example the entire spectrum in the infinite system will be confined to the left half-plane. Thus MAWs should be found in experiments, that can be well described by the CGLE for appropriate control parameters. In this section we present a detailed study of the stability properties of MAWs. MAWs with a single hump per period $`P`$ will be called “single MAWs”. Their existence domains were studied in the previous section. However, the effective interaction between adjacent periods of a single MAW can be repulsive or attractive (see Figs. 11 and 12). Period doubling (PD) bifurcations (open squares) occur at the transitions from repulsive to attractive interaction . There, new branches of MAWs with longer period but many humps per period emerge from the primary branch of single MAWs. We will call these solutions “multi-hump MAWs”. In their profile some humps gain more space and others are compressed in an alternating fashion. The new branches extend to larger $`c_3`$ than the corresponding single MAWs. Fig. 10 shows how these branches arrange in a system with 4 interacting humps ($`L=4P`$). As long as the PD bifurcations are supercritical, the multi-hump MAWs are stable. They represent the saturated solution for attractive interaction between subsequent modulations. For large systems a whole sequence of period doubling bifurcations will lead to multi-hump MAWs with an overall period equal to the system size. Hence they appear as an erratic spatial sequence of humps and depressions. This spatial sequence propagates in a coherent fashion. We named these patterns multi-hump MAWs to emphasize the connection among the coherent structures.

Examples of these stable aperiodic patterns were already observed in numerical simulations. R. Montagne et al.  denote this behavior as “frozen phase turbulence” while A. Torcini et al.  use the term “solutions of type $`\beta `$”.

The observed coexistence of a large number of stable multi-hump MAWs results in a strong dependence of the final state on the initial conditions of the numerical simulation. Although each regular final configuration must be consistent with a particular single or multi-hump MAW it is difficult to predict how the selected final patterns depend on the initial conditions.

MAWs with large period $`P`$ undergo a “splitting” instability as in the limit case $`\nu =0`$ . Roughly, the spatial profiles of these MAWs consist of a localized hump and a plane wave part. Since the extended plane wave is linearly unstable the splitting instability is reminiscent of the Eckhaus instability. It creates more humps on the plateau of the unstable MAW and reduces the period $`P`$ of MAWs on average.

Figs. 11 and 12 represent cuts through the parameter space at fixed $`\nu =0.25`$ and $`c_3=0.5`$, respectively. They show the typical arrangement of stable and unstable parameter regions for single MAWs. Other examined cuts for $`c_3=0.1,0.2,0.3,0.4,0.6`$ and $`0.7`$ qualitatively show the same order.

The cut through parameter space $`c_3,\nu ,P`$ at $`\nu =0.25`$ is shown in Fig. 11. The HB (dashed curve in the figure) approaches the Eckhaus instability for $`P\mathrm{}`$ as the lower bound of the existence domain. From above the domain is limited by the SN (solid curve). For small $`P`$ (large $`c_3`$) the HB is subcritical and only unstable upper branch MAWs exist to the left. In the infinite system MAWs are found to be linearly stable for a broad range of parameters (dark shaded area). At low $`P`$ the interaction instability occurs (white area) whereas at large $`P`$ the long plateau of the MAW is unstable to splitting (light shaded area). For $`c_3<0.45`$ most random initial conditions will evolve to stable MAWs.

A cut perpendicular to the previous one is shown in Fig. 12. Curves and shadings have the same meaning as discussed above. Starting from random initial conditions at $`\nu >0.1`$ a transient may again lead to a stable MAW with local periods $`P`$ inside the stable windows. At lower $`\nu <0.1`$ the probability of approaching a stable configuration decreases since only a third of the previous stable $`P`$ intervals remains. Below $`\nu =0.02`$ no stable state can be prepared at all. Instead one observes wound-up phase chaos with an associated maximal Lyapunov exponent that increases for decreasing $`\nu `$. As in the limit case $`\nu =0`$ (phase chaos), the dynamics is driven by the attractive interaction and annihilation of localized modulations in competition with the splitting instability that produces new peaks in the modulations. In particular for decreasing $`\nu `$, the splitting instability extends to shorter periods $`P`$ and significantly overlaps with the interaction instability. With the above arguments many results obtained by numerical simulations of the CGLE can be well interpreted. In particular it has been observed in  for the same choice of parameters ($`c_1=3.5`$ and $`c_3=0.5`$) that the maximal Lyapunov exponent (averaged over many different initial conditions) is positive for $`\nu =0`$ and decreases monotonously towards zero for increasing $`\nu `$. Above $`\nu =0.09`$ no chaotic solutions have been observed.

## 3 Defect formation in wound-up phase chaos

In this section the formation of defects and the resulting change of the average phase gradient are studied. For $`\nu 0`$, the scenario of defect formation past the saddle-node bifurcation of the relevant MAW is analogous to the previously studied case $`\nu =0`$ , see Subsection 3.1. In particular, the dependence of the final selected average phase gradient $`\nu _f`$ on the initial value $`\nu _i`$ found in numerical simulations  can be interpreted. Subsection 3.2 is then devoted to the limit $`\nu _M`$ of wound-up phase chaos. For a certain range of parameters the limit $`\nu _M`$ is reproduced by means of the stability properties of MAWs at the saddle-node bifurcation. These arguments work well for defect creation with $`\nu >0.1`$. At smaller values of $`\nu `$ various instabilities (splitting and interaction) of MAWs compete and a general statement is more difficult, compare also .

### 3.1 Beyond the saddle-node bifurcation

The role of the SN bifurcation for the dynamics has been studied in  for the limit case $`\nu =0`$. For $`\nu 0`$ we find similar behavior. Fig. 13 gives examples for $`\nu =0.25,P=2\pi /\nu `$. Perturbations of a plane wave lead to defects only above the SN, whereas below the SN such perturbations have to be very large to overcome the saddle-type upper branch MAW.

There are no SNs for parameters below the SN corresponding to $`P\mathrm{}`$. Thus, starting from random initial conditions defects may only form at parameters above the SN of $`P\mathrm{}`$. The SN of $`P\mathrm{}`$ represents a lower bound for defect formation.

For large systems the formation of defects depends on the local period of initial perturbations in a similar way as for $`\nu =0`$ . The peak to peak distances of $`\phi _x(x,t)`$ are used to determine local periods $`p`$. In that context, the following explanation of defect formation has been proposed. Defects are observed in the phase turbulent regime whenever local structures, similar to MAWs, with spatial periods $`p`$ larger then $`P_{SN}`$ occur in the system. Where $`P_{SN}`$ denotes the period of the MAWs at the SN (that coincides with the maximal MAW-period) for the considered choice of parameters. For larger values of $`\nu `$ or $`c_3`$ the SN occurs for smaller $`P_{SN}`$ as shown in Figs. 11 and 12. Therefore, at larger $`\nu ,c_3`$ local periods $`p`$ beyond the SN and subsequent defect formation are more probable. In contrast to the case $`\nu =0`$, there is only a short transient of phase chaos in the simulations with nonzero initial $`\nu _i>\nu _M`$. The distribution of local periods $`p`$ of the perturbations is given by the realization of the noise in the initial condition. For local periods slightly above the relevant SN (as in Fig. 13e), the perturbation first increases to a modulation similar to MAWs and appears almost saturated for some transient time until finally a defect appears. This transient of defect formation becomes shorter as the distance to the SN grows. If initial conditions are prepared with $`\nu _i\nu _M`$, then some local periods will be far beyond the corresponding SNs. The transients of defect formation are shorter in this case.

These two observations suggest an interpretation of the curve $`\nu (t)`$ representing the temporal evolution of $`\nu `$ during transients with $`\nu _i>\nu _M`$ . The larger the initial $`\nu _i`$ is chosen the smaller is the final value $`\nu _f`$ . The time scales of competing processes have to be considered. Local defect formation will not instantaneously effect distant spatial locations along the system. Instead the local change of the average phase gradient $`\nu `$ via a defect will take a transient time to relax over the entire system. For $`\nu _i`$ much larger than $`\nu _M`$, defect formation happens on a short time scale and independently leads to defects at many different spatial locations before the relaxation of the decreased average phase gradient $`\nu _f\nu _M`$ can stop defect formation. For $`\nu _i`$ just above $`\nu _M`$ defects form slowly and the reduced $`\nu _f`$ can relax the phase gradient at distant locations before other defects occur.

Let us now verify if the mechanisms proposed to explain defect formation in the phase chaos regime for solutions with $`\nu 0`$ still hold for the wound-up phase chaos regime. For $`\nu 0`$, defects form if and only if the period $`p`$ of local structures is larger than $`P_{SN}`$, where $`P_{SN}`$ is the period for which a SN occurs at the chosen coefficients $`c_1`$ and $`c_3`$ . We have considered two cuts in the parameter space at $`c_3=0.5`$ and $`c_3=0.65`$ and investigated the distributions of periods $`p`$ for solutions with average phase gradient $`\nu \nu _M`$. In particular, 50 realizations of a system of length $`L=512`$ initialized as a plane wave with wavenumber $`\nu `$ plus noise in the amplitude and in the phase have been considered. Then these different initial conditions have been followed for an integration time $`t=500,0001,000,000`$. The last part of the run has been examined in order to extract the length of the coherent structure with the maximal period $`p_{max}`$ occurring during the evolution.

From these simulations we obtain the following: if defects occur then $`p_{max}>P_{SN}`$ in all observed cases. However, it is not true that a defect is formed any time we observe $`p_{max}>P_{SN}`$. If we let the system relax for a long time ($`t=500,000`$) and measure $`p_{max}`$, then the number of initial conditions leading to a $`p_{max}>P_{SN}`$ without defect formation is noticeably reduced. For $`c_3=0.65`$ the maximal conserved phase gradient is $`\nu _M0.1`$ . In the simulations we do not observe a defect for $`\nu =0.086`$ and $`\nu =0.098`$ but in the first case only 2 % of all runs show $`p_{max}>P_{SN}`$, while in the latter case this percentage increases to 8 %. Increasing $`\nu `$ the maximal period $`p_{max}`$ increases and more and more situations with $`p_{max}>P_{SN}`$ are found upon approach to $`\nu _M`$.

The differences to the $`\nu =0`$ case may be explained by the coexistence of chaotic and stable not-chaotic attractors. Depending on the initial conditions, the solution of the CGLE can evolve towards one or the other. Therefore the system may exhibit local structures similar to multi-hump MAWs that possess SN bifurcations at parameter values larger than those of single MAWs (compare Fig. 10). In that case, some periods $`p`$ may even exceed $`P_{SN}`$.

### 3.2 Limit of wound-up phase chaos

For random initial conditions with $`\nu `$ in the narrow range between the SN of $`P\mathrm{}`$ and the existence limit of MAWs (see Fig. 6), it depends on the specific realization of the noise whether a defect can form or a stable MAW results. No defects form below the line $`\nu _M(c_3)`$. In order to understand this observed limit $`\nu _M(c_3)`$ of wound-up phase chaos, it is sufficient to consider the SNs of single MAWs.

Although initial conditions with large $`P`$ beyond a SN could lead to defects, this is often prevented by the action of the splitting instability. Then the period is decreased before a defect can form. Following the SN curve in Figs. 11 and 12 one encounters a transition from SNs with a splitting instability at large $`P`$ to SNs without this instability at short $`P`$. Defect formation in wound-up phase chaos with $`\nu >0.1`$ occurs for parameters where the splitting instability is not present near the SN, i.e. above the dotted line in Fig. 14.

Fig. 14 summarizes the bounds found so far for the limit of wound-up phase chaos. The domain of stable plane waves at low $`c_3`$ is limited by the Eckhaus instability (thin curve). Within the shaded area only supercritical HBs of different period $`P`$ occur but no SNs. This area is limited by the lowest SN curve of $`P\mathrm{}`$ (thick solid curve). No defects do form from random initial conditions within the shaded area. The dashed curve denotes the upper limit of the existence domain of MAWs with any period. Saddle-node bifurcations exist in the window between this dashed curve and the thick solid curve for the SN with $`P\mathrm{}`$. Defects will always form for any choice of initial conditions at parameters above the dashed curve. The dotted curve gives the transition from active (below) to inactive (above) splitting modes at the SN. This transition is computed by linear stability analysis along cuts like Figs. 11 and 12. Splitting inhibits defect formation below this dotted curve. Filled circles correspond to $`\nu _M`$ obtained from numerical simulations as in . Diamonds refer to the transition from chaotic (below) to non-chaotic (above) asymptotic states.

The SN for $`P\mathrm{}`$ (thick solid) is a lower bound for defect formation which also holds in the limit $`\nu =0`$. The point L<sub>1</sub> marks the transition from phase to defect chaos studied earlier .

As long as the dynamics is regular ($`\nu >0.1`$), the upper bound for the limit (onset of splitting at SN, dotted curve) of wound-up phase chaos reproduces well the numerical observations. For chaotic states with $`\nu <0.1`$ defect formation eventually becomes possible even despite the presence of the splitting instability. This coincides with the increasing weight of the instability of single MAWs to attractive interaction of subsequent amplitude peaks.

## 4 Experimental observations

A variety of experimental observations in quasi one-dimensional geometries can be well interpreted by MAWs. These systems shall also serve for testing further properties of MAWs .

(i) For the oscillatory instability of Rayleigh-Bénard convection patterns in an annular cell, B. Janiaud et al. report long living transients of modulated waves that eventually cause defects . The underlying Eckhaus instability was found to be subcritical. Then, we expect the modulation to grow roughly exponentially whereas the long transient of an almost saturated modulation is similar to the dynamics near the saddle-node bifurcation as shown in Fig. 13e. In the latter case the Eckhaus instability is supercritical and stable MAWs may exist for nearby values of the experimental parameters, see also the discussion in . In all cases a single pulse-shaped modulation with the period equal to the cell length was present in the system and led to a defect. The excitation of several modulations per cell and thereby smaller period can stabilize the modulated pattern and provide more examples of the dynamics near the saddle-node bifurcation.

(ii) Hydrothermal waves have been studied in ring-shaped cells  as well as in linear cells . N. Mukolobwiez et al. report a supercritical Eckhaus instability, stable MAWs with the period equal to the cell length and defect formation after a parameter change . N. Garnier et al. observe modulated waves with both the wavenumber and the period of modulations being selected by one of the longitudinal boundaries . Thermal or mechanical forcing at the boundary may yield more insight into the multistability of the two-parameter family of MAWs.

(iii) In rotating Rayleigh-Bénard convection Y. Liu et al. observe the subcritical Eckhaus instability of a traveling wave sidewall mode . The authors suggest higher order corrections to the CGLE in order to explain the observed discrepancy between the linear group velocity and the observed velocity of finite wavelength perturbations. However, this difference already follows from the linear analysis of the Eckhaus instability  and can be used as a test for the assumed coefficients of the CGLE.

(iv) For the Taylor-Dean system I. Mutabazi et al. report a stable MAW that they called “triplet state” because of the length scale ratio of modulation and underlying wave . They also observe the formation of defects. Clearly the triplet state is just one realization of the two-parameter family of MAWs.

(v) Finally we mention the heated wire convection , internal waves  and the oscillatory variant of the Belousov-Zhabotinsky reaction  where further investigations of the observed Eckhaus instabilities may also reveal MAWs.

## 5 Discussion

The bifurcation analysis of modulated amplitude waves (MAWs) has been extended to nonzero average phase gradient ($`\nu 0`$). Small amplitude MAWs (“lower branch”-MAWs) of specific spatial period $`P`$ exist between a supercritical Hopf bifurcation (HB) and a saddle-node (SN) bifurcation. The HB asymptotically reaches the Eckhaus instability from above as $`P`$ goes to infinity. MAWs are a direct consequence of the Eckhaus instability of plane waves; they are obtained by a computer-assisted nonlinear analysis of this instability. We encounter SNs with decreasing values of $`P`$ as $`c_1,c_3`$ and $`\nu `$ are increased. These SNs govern the formation of defects from random initial conditions as well as many aspects of the evolution of wound-up phase chaos. The SNs bound the existence region of MAWs in the Eckhaus unstable regime.

A linear stability analysis of MAWs revealed that they can be linearly stable even in systems of infinite size. These domains are limited by the interaction instability at low and the splitting instability at high values of the spatial period $`P`$ of the MAW. The competition of the two instabilities drives wound-up phase chaos and determines the degree of chaoticity of the dynamics. For fixed coefficients $`c_1`$ and $`c_3`$, the SN associated to $`P\mathrm{}`$ occurs at the lowest value of $`\nu `$. This SN establishes a “lower” bound for defect formation and thereby for the limit $`\nu _M`$ of wound-up phase chaos. The splitting instability can inhibit defect formation if the SN occurs at large $`P`$ hence defects are created more frequently for parameters above a second curve where the splitting instability vanishes at the SN. Thereby an “upper” bound for the limit $`\nu _M`$ of wound-up phase chaos is obtained. Earlier numerical observations on $`\nu _M(c_3)`$ are well reproduced for $`\nu >0.1`$, respectively small $`c_3`$. For $`\nu <0.1`$, the description of phase chaos relies on considerations similar to those already discussed in the limit case $`\nu =0`$ . Finally, several experimental observations were interpreted in terms of MAWs.

We acknowledge fruitful discussions with M. van Hecke, N. Garnier, F. Daviaud, W. van de Water and L. Kramer. AT was partially supported by the MURST-COFIN00 grant “Chaos and localization in classical and quantum mechanics”.