text:
# On the Equation of State for Scalar Field

( )

## Abstract

We consider Friedmann cosmologies with minimally coupled scalar field. Exact solutions are found, many of them elementary, for which the scalar field energy density, $`\rho _f`$, and pressure, $`p_f`$, obey the equation of state (EOS) $`p_f=w_f\rho _f`$. For any constant $`|w_f|<1`$ there exists a two-parameter family of potentials allowing for such solutions; the range includes, in particular, the quintessence ($`1<w_f<0`$) and ‘dust’ ($`w_f=0`$). The potentials are monotonic and behave either as a power or as an exponent for large values of the field. For a class of potentials satisfying certain inequalities involving their first and second logarithmic derivatives, the EOS holds in which $`w_f=w_f(\phi )`$ varies with the field slowly, as compared to the potential.

1. Introduction Quintessence is a dynamic energy with negative pressure-to-density ratio . The idea of quintessence provides new ‘degrees of freedom’ in cosmology , extends the variety of modern field models to include extreme forms of energy , and may stimulate a better understanding of the fundamental problem of interplay between gravity and field theories . Quintessence is also shown to give rise to spacetimes with new interesting properties .

However, one needs to find out what kind of physical substance obeys the quintessence EOS, that is, to model quintessence, in some way. It was common understanding from the very beginning that a scalar field could act as quintessence; nevertheless, it should be clearly demosntrated.

The EOS for the scalar field has been intensively studied in  for the so called tracking cosmological solutions introduced in , and some classes of potentials allowing for the field EOS were described. In this letter we examine this question for any cosmological solutions. We first give an explicit description of all potentials allowing for solutions with scalar field satifying the linear EOS with constant parameter $`w_f`$. Then, by way of generalization, we derive the conditions on the potential which provide the EOS where the parameter varies slowly with the field, $`w_f=w_f(\phi )`$. 2. Exact solutions with scalar field satisfying linear EOS

We consider Friedmann universe described by the Friedmann-Robertson-Walker metric,

$$ds^2=dt^2+a^2(t)\left(\frac{dr^2}{1kr^2}+r^2d\mathrm{\Omega }^2\right),$$
(1)
where $`k=1,0`$, or $`1`$, according to whether the universe is open, flat, or closed ($`G=c=1`$). In the presence of matter of the density, $`\rho `$, and pressure, $`p`$, and a minimally coupled scalar field, $`\phi `$, with the potential $`V(\phi )0`$, the evolution equations can be written in the form:

$`3{\displaystyle \frac{\dot{a}^2}{a^2}}=8\pi \rho +8\pi \rho _f3{\displaystyle \frac{k}{a^2}};`$ (2)
$`\dot{\rho }=\mathrm{\hspace{0.17em}3}{\displaystyle \frac{\dot{a}}{a}}(\rho +p);`$ (3)
$`\dot{\rho _f}=\mathrm{\hspace{0.17em}3}{\displaystyle \frac{\dot{a}}{a}}(\rho _f+p_f).`$ (4)
The proper definitions of the scalar field energy density, $`\rho _f`$, and pressure, $`p_f`$, are:

$$8\pi \rho _f\dot{\phi }^2+V,8\pi p_f\dot{\phi }^2V.$$
(5)
The first of them is implied by equation (2), the second one is validated by the very equation (4), written as the field energy conservation equation, identical in its form with the conservation law (3) for matter. Also, the expression for the acceleration, following from eqs. (2)–(3) and (5), reads:

$`3\ddot{a}/a=4\pi (\rho +3p)+V2\dot{\phi }^24\pi [(\rho +3p)+(\rho _f+3p_f))];`$ (6)
the density and pressure of both matter and field are involved here uniformly.

We assume the linear equation of state for the matter, $`p=w\rho ,w=\mathrm{const}`$. Our main interest is concerned with the ‘normal’ matter with $`w0`$; however, all the following results hold for any $`w1`$. Using the EOS, equation (3) immediately integrates to express the density via the scale factor ($`C>0`$ is an arbitrary constant):

$$8\pi \rho =C/a^{3(1+w)}.$$
(7)
Global dynamics of cosmological expansion governed by the system (2)–(4) was described in . In this paper we are looking for special solutions subject to an additional constraint, namely, a linear EOS for the scalar field, $`p_f=w_f\rho _f`$. First, we assume $`w_f`$ constant, and look for potentials which would allow solutions satisfying such EOS. Introducing expressions (5) converts the EOS into an equivalent relation,

$$\dot{\phi }^2=\frac{1+w_f}{1w_f}V(\phi ),$$
(8)
showing that the acceptable range of the EOS parameter is $`|w_f|<1`$. Evidently, the EOS of the Zel’dovich superstiff fluid, $`w_f=1`$, holds only for the free field, $`V(\phi )0`$; the vacuum EOS ($`w_f=1`$) is only valid for the constant field $`\phi _c`$ such that, by (4), $`V^{^{}}(\phi _c)=0`$. This is the classical case of a universe with the cosmological constant, studied recently in detail in .

For any given $`V(\phi )`$ and $`w_f`$, equation (8) completely specifies the time evolution of the scalar field. On the other hand, the EOS of the scalar field allows for an immediate integration of the field energy conservation equation (4), in exactly the same fashion as with the equation (3) for the matter, so that

$$8\pi \rho _f=C_f/a^{3(1+w_f)};$$
(9)
here $`C_f>0`$ is another arbitrary constant of motion. Combining this and (8), one relates the potential (as a function of time for the solution under investigation) to the scale factor,

$$V=(1w_f)C_f/2a^{3(1+w_f)}.$$
(10)
At this point, only the equation (2), of the whole system (2)–(4), remains to be satisfied. By (7) and (9), it transforms into the first order equation for the scale factor (the chosen positive sign corresponds to expansion):

$$\dot{a}=\left(C/3a^{1+3w}+C_f/3a^{1+3w_f}k\right)^{1/2}.$$
(11)
This, in fact, completes the solution: integration of equation (11) gives the scale factor as a function of time, next, relation (10) provides the potential, also as a function of time. Scalar field $`\phi (t)`$ is then determined by integrating a known function of time:

$$\dot{\phi }(t)=\pm \sqrt{\frac{1+w_f}{1w_f}V}=\pm \sqrt{\frac{(1+w_f)C_f}{\mathrm{\hspace{0.17em}\hspace{0.17em}2}[a(t)]^{3(1+w_f)}}}.$$
(12)
Formula (12) shows that $`\phi =\phi (t)`$ is a monotonic function of time, so that the inverse to function $`t=t(\phi )`$ is well defined. Hence, the potential can be finally determined as the function of the field, $`V(\phi )=(1w_f)C_f/2\left[a(t(\phi ))\right]^{3(1+w_f)}`$, according to the expression (10).

Using the EOS for both matter and scalar field, and then (7) and (9), the acceleration equation (6) for the found solution becomes

$$3\ddot{a}/a=(1/2)\left[(1+3w)C/a^{3(1+w)}+(1+3w_f)C_f/a^{3(1+w_f)}\right].$$
(13)
Even if the matter is ‘normal’, $`w0`$, but the field acts as quintessence with $`1<w_f<1/3`$, the expansion clearly accelerates, at least at large enoguh values of $`a`$. Also, if $`w_f0`$, then the scalar field can play the role of dark matter in the universe (whether the scalar field can clamp to ‘normal’ matter, as the dark matter seemingly does, requires separate investigation).

The found solution contains, along with $`C,C_f`$, two more arbitrary constants coming from integration of (11) and (12). Both of them are not significant, since they just shift the values of time and scalar field. Thus we have demonstrated that, given any $`|w_f|<1`$ (and $`w>1`$), a two–parameter family of potentials $`V=V(\phi ,C,C_f)`$ and the corresponding cosmological solution exists, for which the EOS of the scalar field holds throughout the whole evolution. If $`1<w_f<0`$, then the scalar field acts exactly as quintessence, and, consequently, for $`1<w_f<1/3`$ the expansion of the universe accelerates, at least at its late stages.

3. Explicit description of the potential

The above derivation lacks a direct formula relating the potential with the scalar field, $`V=V(\phi )`$. The formula is not difficult to establish: one just writes the evolution equation (2) not in terms of the scale factor, as in (11), but rather in terms of the potential. Indeed, from (10), one has $`a=[(1w_f)C_f/2V]^{1/3(1+w_f)}`$; differentiating this in time using (8), one obtains also a formula for $`\dot{a}/a`$ through $`V`$ and $`dV/d\phi `$. Replacing all the terms in the evolution equation (11) with their ready expressions via $`V`$ results in the first order differential equation for $`V(\phi )`$:

$$Q(V,w_f)dV/V=\sqrt{6(1+w_f)}d\phi ,$$
(14)
where

$`Q(V,w_f)=\left[1+AV^{\frac{ww_f}{1+w_f}}kBV^{\frac{1+3w_f}{3(1+w_f)}}\right]^{1/2}`$ (15)
$`A=\left(C/C_f\right)\left[(1w_f)C_f/2\right]^{\frac{w_fw}{1+w_f}},B=\left(3/C_f\right)\left[(1w_f)C_f/2\right]^{\frac{1+3w_f}{3(1+w_f)}},`$ (16)
and the signs correspond to the signs in (12).

After a straightforward integration of (14) (that leads, in many cases, to an elementary result, see the next section), we arrive at the desired expression for the potential as function of the scalar field, $`V=V(\phi )`$. The dynamics of the field, i. e., the function $`\phi =\phi (t)`$, is then determined by integrating the first expression in (12), and the scale factor evolution is found by its expression through $`V(\phi (t))`$ above. This expression shows that always $`V=\mathrm{}`$ at the Big Bang ($`a=0`$), and $`V0`$ when $`a\mathrm{}`$; it also demonstrates that the potential is a monotonic function of scalar field. Moreover, equation (14) allows us to specify the behavior of the potential in both limits; the most interesting are the late stages of the cosmological expansion.

If the EOS parameter satisfies

$$w_f\mathrm{min}\{w,1/3\},$$
(17)
then the potential drops exponentially, $`V(\phi )\mathrm{exp}(\alpha \phi )`$. The constant $`\alpha >0`$ equals to $`\sqrt{6(1+w_f)}`$ when no equality holds in (17), and is expressed through all the parameters otherwise. So, under the condition (17) the found solutions implement the runaway scenario , with $`\phi \pm \mathrm{},V(\phi )0`$ at large times. Remarkably, all this happens exactly in the case when the scalar field behaves as a ‘strong’ quintessence with the pressuere equally or more negative than the one of the Einstein quitessence, for which $`w_f=1/3`$. When condition (17) is invalid, the potential behaves as a power of the scalar field when it is small, $`V(\phi )\phi ^\nu `$; here $`\nu >0`$ depends only on $`w`$ and $`w_f`$.

Note, however, that the large–time behavior of the scale factor $`a(t)`$ and the potential $`V(\phi (t))`$ \- as fucntions of time – is qualitatively the same for all possible parameter values: these functions are, respectively, a positive and negative power of time related by (10). It is only the scalar field, $`\phi (t)`$, that evolves differently at large times, depending on whether condition (17) is true or not: in the first case, the field is logarithmic, in the second case it has a power asymptotics, same as the other two dynamic variables.

4. Examples of closed-form solutions

Several classes of physically meaningful closed–form solutions are described by the same formula for the potential, differing just by the values of the parameters involved. The reason for this is the structure of the key equation (14): there are two terms added to unity under the square root in $`Q(V,w_f)`$ \[see (15)\]. In all the cases when either of the terms is a constant, the integral of (14) is elementary and, of course, of the same functional form. We give the expression only for the potential \[if $`V(\phi )`$ is known, $`\phi (t)`$ and $`a(t)`$ are found as described in the previous section\]. This potential, allowing for solutions with scalar field satisfying a linear EOS, is

$$V(\phi )=V_0\left[\mathrm{sinh}^2m(\phi \phi _{})\right]^\mu ,\phi _{}=\mathrm{const}.$$
(18)
The expression is valid in the following five cases : 1) flat universe, $`k=0`$; 2) universe with scalar field only, $`C=0`$; 3) same EOS for both matter and field, $`w=w_f`$; 4) scalar field as Einstein quintessence, $`w_f=1/3`$; 5) Einstein quintessence as matter, $`w=1/3`$. The cases differ only by the values of $`\mu `$, $`V_0`$, and $`m`$, whose expressions are cumbersome. Thus we show just one example, for the first case of flat universe:

$$\mu =\frac{1+w_f}{w_fw},V_0=\frac{(1w_f)C_f}{2}\left(\frac{C}{C_f}\right)^\mu ,m=|w_fw|\left[\frac{3}{2(1+w_f)}\right]^{1/2}.$$

In each of the above cases it is taken that only one of the five assumptions designating the cases is fulfilled (for instance, in the first case $`C0,w_fw,w_f1/3`$, $`w1/3`$, etc.). On the other hand, many combinations of two of those assumptions (such as $`k=0`$ and $`C=0`$, $`k=0`$ and $`w=w_f`$, etc.) imply that $`Q(V,w_f)`$ in the equation (14) becomes a constant. Therefore, the potential proves to be exactly an exponent in every such case ,

$$V(\phi )=V_0\mathrm{exp}(\alpha \phi ),$$
(19)
with arbitrary $`V_0`$ and $`\alpha >0`$ depending on the parameters. The scalar field then is a logarithmic function of time, which results in the potential proportional to the inverse square of time, $`V(\phi (t))t^2`$; the scale factor is also a power, $`a(t)t^{2/3(1+w_f)}`$.

Many other closed–form solutions are not related to the degeneracy of $`Q(V,w_f)`$. To give one example, for pressureless matter ($`w=0`$) and scalar field obeying the EOS of quintessence with $`w_f=2/3`$ in the open universe ($`k=1`$), the potential is:

$$V(\phi )=4\overline{e}(\phi )\left\{\left[\overline{e}(\phi )B\right]^24A\right\}^1,\overline{e}(\phi )\mathrm{exp}\left[\pm \sqrt{6(1+w_f)}(\phi \phi _{})\right],$$
(20)
where $`A=6C/5C_f^3,B=2/5`$, according to (16) for this particular case. The stability of cosmological solutions driven by quitenssence with the pressure to density ratio $`(2/3)`$ has been recently studied in .

5. Generalizations. Linear EOS with slowly varying parameter, $`w_f(\phi )`$

Generalization of the above results to the case of arbitrary number, $`N`$, of non-interacting matter species and arbitrary number, $`N_f`$, of non-interacting scalar fields $`\phi _n`$ satisfying (different) linear EOS, is straightforward. Such solutions can, in particular, fit the concordant cosmological data, if $`N=N_f=2`$ and $`w^{(1)}=0`$ for baryonic matter, $`w^{(2)}=1/3`$ for radiation, scalar field $`\phi _1`$ plays the role of the dark matter with $`w_f^{(1)}=0`$, and scalar field $`\phi _2`$ acts as quintessence with $`w_f^{(2)}=1÷0.6`$ .

Also, the EOS of the field holds approximately for long periods of time for the potentials which are small perturbations of those found in the previous sections. However, the scalar field can satisfy the EOS approximately in a much broader range of situations. Suppose the parameter in the EOS for the field is not a true constant, but a function of the field, $`w_f=w_f(\phi )`$, so that the EOS, or the equivalent equation (8), is no longer a constraint, rather, a definition of a new dynamic variable. The only constraint we now impose is that this variable, $`w_f(\phi )`$, varies slowly with the field as compared to its density $`\rho _f`$, so that the integration of the field energy conservation equatuion (4) resulting in (9) is still approximately valid. By (8), the condition of slow variation proves to be

$$\frac{2}{1w_f}\left|\frac{dw_f}{d\phi }\right|\frac{1}{V}\left|\frac{dV}{d\phi }\right|,\mathrm{or}\mathrm{just}\left|\frac{dw_f}{d\phi }\right|\left|\frac{d\mathrm{ln}V}{d\phi }\right|$$
(21)
provided that $`w_f(\phi )`$ is not too close to unity; particularly, the second form of the condition is true for all the negative values of $`w_f(\phi )`$.

Under the condition (21), all calculations of sections 2 nad 3 hold to lowest order in the slow variation. The key equality (14), which we now write as

$$w_f=\frac{Q^2(V,w_f)}{6}\left(\frac{d\mathrm{ln}V}{d\phi }\right)^2\mathrm{\hspace{0.17em}1},$$
(22)
is no longer a differential equation for $`V(\phi )`$, but, for a given potential, a transcendetal equation for $`w_f(\phi )`$, instead. If its proper solution exists, then condition (21) becomes an inequality which limits the class of potentials in question.

Since this is an important general result, we formulate it accurately. Suppose that for some values of $`w`$, $`k=0,\pm 1`$ and constants $`C,C_f`$ involved in $`Q`$, equation (22) has a solution $`w_f=W(V,d\mathrm{ln}V/d\phi )`$. Then the inequalities

$$W(V,\frac{d\mathrm{ln}V}{d\phi })<1,\left|\frac{d}{d\phi }W(V,\frac{d\mathrm{ln}V}{d\phi })\right|\left|\frac{d\mathrm{ln}V}{d\phi }\right|$$
(23)
describe all the potentials for which cosmological scalar field satisfies a linear EOS with the parameter varying slowly with the field. To lowest order in this slow variation, the dynamics of cosmological expansion, i. e., the functions $`\phi (t)`$ and $`a(t)`$, are found successively from (12) and (10). If the solution $`w_f=W(V,d\mathrm{ln}V/d\phi )`$ and conditions (23) are valid not for all the values of the field, but only for some range of it, then, accordingly, the EOS holds only for the corresponding period of cosmological evolution, with the values of $`\phi `$ within this range.

Equation (22) is generally very complicated \[recall that $`w_f(\phi )`$ is involved in the powers of $`V`$, as well as in the ‘constants’ $`A`$ and $`B`$, see (15), (16)\] and should be studied in detail separately. However, it dramatically simplifies for later stages of cosmological expansion under the condition (17), which includes the most interesting case $`w0,w_f<1/3`$ (‘normal’ matter and ‘strong’ quintessence). Indeed, at the later stages the scale factor is large, $`a1`$, and, by (10), the potential is small, $`V/C_f1`$. Due to this and condition (17), the factor $`Q`$ gets close to unity, so that (22) turns into an approximate explicit expression for $`w_f(\phi )`$ through the potential,

$$w_f=\frac{1}{6}\left(\frac{d\mathrm{ln}V}{d\phi }\right)^2\mathrm{\hspace{0.17em}1},$$
(24)
and conditions (23) also become explicit and rather simple:

$$\left|\frac{d\mathrm{ln}V}{d\phi }\right|<2,\left|\frac{d^2\mathrm{ln}V}{d\phi ^2}\right|3.$$
(25)
The set of potentials satisfying (25) is rather wide. For instance, $`V(\phi )=V_0(\phi ^2+b^2)\mathrm{exp}(\alpha \phi )`$ with $`b^21`$, $`|\alpha |<2`$,and $`\alpha 0`$ (to make the potential small at large $`|\phi |`$ of the proper sign), satisfies (25) for all values of $`\phi `$, with the varying part of $`w_f`$ of the order of $`1/b^2`$,

$$w_f=1+\alpha ^2/6+O(1/b^2).$$

The factor $`(\phi ^2+b^2)`$ in this example can evidently be replaced with any polynomial of an even degree and with no real roots, under a single additional constraint on its coefficients guaranteeing the smallness of its second logarithmic derivative. Next, a ratio of two such polynomials can be taken, in which case the parameter $`\alpha `$ in the exponent could even be zero, if the degree of the polynomial in the denominator is larger than in the numerator, etc.

All the above examples correspond to the runaway scenario $`\phi (t)\mathrm{},V(\phi (t))0`$ at $`t+\mathrm{}`$. This is true in general for the solution (22), (25) at $`V/C_f1`$. Indeed, if the scalar field goes to a finite value $`\phi _0`$ at large times, then, by (10), $`V(\phi _0)=0`$. Since the potential is nonnegative, $`\phi _0`$ is the point of its minimum, and $`V^{^{}}(\phi _0)=0`$. But then the logarithmic derivative $`(\mathrm{ln}V)^{^{}}`$ tends to infinity when $`\phi (t)\phi _0`$, and equation (23) cannot be satisfied in the vicinity of $`\phi _0`$. In fact, the EOS of the scalar field proves to be nonlinear in this vicinity.

Thus we have shown that if matter satifies the EOS with $`w>1/3`$, and the potential drops at large enough values of the scalar field in such way that inequalities (25) hold, then the runaway scenario is possible with the field obeying the EOS of a ‘strong’ quintessence, $`1<w_f(\phi )<1/3`$, at later stages of expansion. This regime takes on the earlier in the expansion, the larger the abundance of the field (the ratio $`C_f/C`$), and the closer its EOS to the EOS of vacuum \[$`w_f(\phi )`$ to $`1`$\] is. So, under conditions (25) on the potential going to zero at large values of $`\phi `$, scalar field quintessence dominates the expansion of the universe at its later stages in all the runaway solutions, and the expansion is then accelerating.

Acknowledgments This work was supported by NASA grant NAS 8-39225 to Gravity Probe B. The author is grateful to A.D.Chernin, D.I.Santiago, and R.V.Wagoner for their valuable remarks.


format: text

text:
# 115In-NQR evidence for unconventional superconductivity in CeIn3 under pressure

## I Introduction

Since the discovery of cerium (Ce) based heavy-fermion (HF) superconductivity in CeCu2Si2 , many experimental works on CeCu2Ge2, CeIn3  and CePd2Si2  have suggested that the antiferromagnetism and superconductivity are closely related each other. The discovery of pressure ($`P`$) -induced HF superconductivity in Ce-based HF antiferromagnetic (AF) compounds has stimulated further experimental works under $`P`$. Recently, a new HF material CeRhIn5 consisting of alternating layers of CeIn3 and RhIn2 was found to reveal an AF to superconducting (SC) transition at a relatively lower critical pressure $`P_c=`$ 1.63 GPa than in previous examples. The SC transition temperature $`T_c`$ = 2.2 K is the highest one to date among $`P`$-induced superconductors. This finding has opened a way to investigate the $`P`$-induced evolution of both magnetic and SC properties over a wide $`P`$ range. In the previous paper , the 115In-NQR study of CeRhIn5 has clarified the $`P`$-induced anomalous magnetism and unconventional superconductivity. At $`P`$ = 2.1 GPa, the nuclear spin-lattice relaxation rate $`1/T_1`$ reveals a $`T^3`$ dependence below the SC transition temperature $`T_c`$, which shows the existence of line-node in the gap function. Note that the unconventional superconductivity in CeMIn5($`M`$ = Co, Rh, Ir) occurs at the border of the AF phase.

It is suggested that $`P`$-induced superconductivity in CeIn3, which is the host crystal of CeMIn5, is also mediated by AF magnetic fluctuations from resistivity measurements under $`P`$. However, evidence for a bulk superconductivity and whether it is mediated by AF magnetic fluctuations or not in CeIn3 has not been established because its superconductivity under $`P`$ is observed only in the resistivity measurements. It is still unclear whether any problem concerning a sample surface or its quality that could prevent the onset of bulk $`P`$-induced superconductivity is ruled out or not. This is in contrast to the case of CeRhIn5 where the bulk nature of the SC phase was fully established  over a wide $`P`$ range 2 $``$ 5 GPa . Note that the $`P`$-induced superconductivity in CeRhIn5 occurs in the non-Fermi-liquid state in the whole $`P`$ range.

CeIn3 crystallizes in the cubic AuCu3 structure and orders antiferromagnetically ($`T_N=10`$ K) at $`P=0`$ with an ordering vector Q = (1/2,1/2,1/2). $`T_N`$ monotonically decreases with $`P`$ and superconductivity appears below $`T_c`$ 0.2 K at a critical pressure $`P_c=2.55`$ GPa and the onset of superconductivity has been observed by the resistivity measurements and is limited in a narrow $`P`$ range of about 0.5 GPa. Also it is reported that non-Fermi-liquid regime in CeIn3 under $`P`$ is limited in a quite narrow $`P`$ range around the peak of resistive $`T_c`$ in the SC phase. The previous 115In-NQR result on CeIn3 revealed that the AF order disappears around $`P^{}`$ = 2.44 GPa close to $`P_c`$ and the Fermi-liquid behavior was observed below 10 K at $`P=2.74`$ GPa . In the previous paper , we reported that the localized magnetic character is robust against the application of $`P`$ up to 1.9 GPa, beyond which $`T^{}`$, which marks an evolution into an itinerant magnetic regime, increases rapidly. It is of great interest that the superconductivity in CeIn3 emerges in such an itinerant regime. The window for the 3D AF fluctuation regime with respect to an external pressure is much narrower in CeIn3 than in CeRhIn5, which may also be partly responsible for the narrower SC region in CeIn3.

In this paper, we report measurements of $`1/T_1`$ and ac-susceptibility at $`P`$ = 2.65 GPa down to $`T`$ = 50 mK which provide evidence for unconventional and bulk nature of superconductivity in CeIn3. At $`P`$ = 2.65 GPa, from the relation of $`T_1T`$ = const. below $`T_{FL}`$ 5K characteristic for the Fermi-liquid state, we conclude that the $`P`$-induced superconductivity occurs well below $`T_{FL}`$ without any signature for a non-Fermi-liquid behavior governed by AF spin-fluctuations. Our result is consistent with the recent resistivity measurements. This is in contrast to the case of CeRhIn5 where superconductivity occurs $`\mathrm{𝑜𝑛𝑙𝑦}`$ in the non-Fermi-liquid regime.

## II Experimental

The single crystals of CeIn3 were grown by Czochralski method as reported elsewhere  and moderately crushed into grains in order to make rf pulses penetrate into samples easily. To avoid some crystal distortions, however, the size of the grains is kept with diameters larger than 100 $`\mu `$m. In a small piece of CeIn3 from the same batch used in this work, the zero resistance transition was confirmed in the $`P`$ range of 2.2 $``$ 2.8 GPa. The phase diagram inferred from the resistivity data is shown in the inset of Fig. 4b, which is in good agreement with the previous reports. At $`P=`$ 2.5 GPa that is close to a critical pressure $`P_c`$, the $`T_c`$ determined by the resistivity measurement reaches a maximum value of $`T_c`$ = 0.19 K and its transition width becomes the sharpest. Measurement of $`1/T_1`$ was made by the conventional saturation-recovery method in the temperature ($`T`$) range of 0.05 $``$ 100 K at $`P`$ = 2.65 GPa. The 115In-NQR $`1/T_1`$ was measured at the transition of 3$`\nu _Q`$ ($`\pm 5/2\pm 7/2`$) above $`T`$ = 1.4 K, but at 1$`\nu _Q`$ ($`\pm 1/2\pm 3/2`$) below $`T`$ = 1.4 K. The ac-susceptibility (ac$`\chi `$) was measured by using an in-situ NQR coil. Hydrostatic pressure was applied by utilizing a NiCrAl-BeCu piston-cylinder cell, filled with Si-based organic liquid as a pressure-transmitting medium. To calibrate the pressure at low temperatures, the shift in $`T_\mathrm{c}`$ of Sn metal under pressure was monitored by the resistivity measurement.  To reach the lowest temperature of 50 mK, a 3He-4He dilution refrigerator is used.

## III experimental result

Figures 1 and 2 show the respective $`T`$ dependencies of $`1/T_1`$ and $`1/T_1T`$ at $`P=`$ 0, 2.35, and 2.65 GPa. At $`P`$ = 0, the AF order occurs at $`T_N`$ = 10 K without any clear signature for an itinerant HF regime. Rather, in the paramagnetic state, the 4$`f`$-electron derived moments behave as if localized, since $`1/T_1`$ stays a constant. It was reported that $`T_N`$ monotonically decreases with pressure and eventually is reduced to $`T_N`$ = 5 K at $`P`$ = 2.35 GPa. At this pressure, it is noteworthy that $`1/T_1`$ starts to decrease below $`T^{}`$ 15 K. This is defined as a characteristic temperature below which the system crosses over to the itinerant HF regime. At pressures exceeding $`P`$ = 2.35 GPa, the internal field at the 115In site due to the AF ordering decreases, revealing a same $`P`$ dependence as $`T_N`$.  At $`P`$ = 2.65 GPa, it is remarkable, as presented in Figs. 1 and 2, that $`1/T_1`$ shows a significant decrease below $`T^{}`$ 30 K, followed by a behavior of $`T_1T`$ = const. that is characteristic for the Fermi-liquid state to be realized below $`T_{FL}=5`$ K. This result proves that the Fermi-liquid description is valid over a wide $`T`$ range 0.1 $``$ 5 K at $`P`$ = 2.65 GPa. Any trace of AF spin-fluctuations at $`P`$ = 2.65 GPa is not appreciable due to an increase of the hybridization of Ce-4$`f`$ electrons and conduction electrons as discussed in the previous paper. This is in good agreement with the recent resistivity measurement.

Next, we discuss about the superconductivity in CeIn3 at $`P`$ = 2.65 GPa. The SC transition at $`P`$ = 2.65 GPa is evidenced by the $`1/T_1`$ measurement. As seen in Figs. 1 and 2, below $`T_c`$ = 95 mK, $`1/T_1`$ decreases suddenly, signaling the onset of bulk SC transition. In addition, the absence of the coherence peak in $`1/T_1`$ just below $`T_c`$ is consistent with an anisotropic SC gap. The $`T_1`$ measurement below 50 mK was too hard to determine the nodal structure because of a heat-up problem arising from the exciting rf pulses.

The SC transition is also corroborated by the ac$`\chi `$ measurement. The $`T`$ dependence of the ac$`\chi `$ at $`P`$ = 2.65 GPa is shown in Fig. 3. The absolute value of susceptibility was determined by using a single crystals of CeIrIn5 as a reference in which bulk superconductivity is fully established. The ac$`\chi `$ decreases below $`T_c`$ 95 mK due to the SC diamagnetism, indicating a clear signature of a SC transition. The diamagnetic signal of CeIn3 reaches 80$`\%`$ of the value of CeIrIn5. It indicates the $`P`$-induced superconductivity in CeIn3 is the bulk one. Note that the $`T_c`$ probed by the $`T_1`$ and ac$`\chi `$ measurements is lower than $`T_c^\rho `$ $``$ 0.15 K at zero resistance on the same sample. This suggests that the surface of the sample makes a short cut for the electrical current flow, exhibiting a rather higher $`T_c`$ value than that for the bulk one. This phenomenon is similar to the SC characteristics in CeIrIn5, in which the bulk $`T_c`$ 0.4 K confirmed by the measurements of the specific heat, the susceptibility, and the NQR measurements  is much lower than $`T_c`$ 1 K determined by the zero resistance. A reason why this inconsistency occurs is still unresolved.

Note that, at $`P`$ = 2.65 GPa, the normal state is dominated by the Fermi-liquid excitations below $`T_{FL}=5`$ K that is far above $`T_c=95`$ mK, but not by AF spin-fluctuations. The superconductivity occurs thus in the Fermi-liquid regime, contrasting with the unconventional superconductivity in CeCu2Si2 and CeMIn5 both of which reveal a higher $`T_c`$ near the border of AF and SC transitions in the non-Fermi-liquid state. Fig. 4a indicates the $`T`$ dependence of $`1/T_1T`$ in CeRhIn5 at $`P`$ = 2.1 GPa where the superconductivity is fully established. This $`T`$ dependence was consistent with the AF spin-fluctuation model near a magnetic criticality. Most remarkably, the $`P`$-induced superconductivity in CeRhIn5 occurs under the dominance of strong AF spin-fluctuations near the border to the AF phase. As seen in Figs. 4a and 4b, the behavior of $`1/T_1T`$ = const. in CeIn3 strongly contrasts with that in CeRhIn5. This difference in magnetic fluctuations is suggested to be intimately related to the fact that the maximum $`T_c^{max}=2.1`$ K in CeRhIn5 is an order magnitude larger than $`T_c^{max}`$ = 0.2 K in CeIn3. Note that $`T_c`$ is much smaller for CeIn3 because non-Fermi-liquid and AF spin-fluctuations regime is limited in quite narrow $`P`$ range than that in CeCu2Si2 and CeRhIn5. For the latter, the superconductivity occurs at the border of the AF phase where the non-Fermi liquid takes place and even coexists with it. This may be related to the robustness of SC region against increasing $`P`$ for the latter.  It is an issue to study in future whether or not the superconductivity and the antiferromagnetism in CeIn3 coexist and whether or not the AF magnetic fluctuations in the non-Fermi-liquid state raise $`T_c`$ near the border of AF and SC phases.

## IV conclusion

In conclusion, we have provided evidences for the bulk and unconventional superconductivity in CeIn3 below 95 mK at $`P`$ = 2.65 GPa on the basis of the measurements of $`T_1`$ and ac-susceptibility. We have found that the system changes the magnetic character from the localized to itinerant regime below $`T^{}`$ 30 K and the Fermi-liquid description becomes valid below $`T_{FL}`$ 5K without any appreciable trace of AF spin-fluctuations developing. $`1/T_1`$ reveals no coherence peak just below $`T_c`$ = 95 mK, indicative of unconventional superconductivity. We remark that the $`P`$-induced SC phase in CeIn3 is characterized by a lower $`T_c`$ and a narrower SC region than in CeRhIn5 where SC phase coexists with the AF phase and $`T_c=2.1`$ K is obtained in virtue of the non-Fermi-liquid state dominated by strong AF spin-fluctuations.

## Acknowledgements

We thank K. Ishida and C. Thessieu for useful discussions. This work was supported by the COE Research (10CE2004) in Grant-in-Aid for Scientific Research from the Ministry of Education, Sport, Science and Culture of Japan. One of the author ( Yu Kawasaki ) has been supported by JSPS Research Fellowships for Young Scientists.


text:
# Response of parametrically-driven nonlinear coupled oscillators with application to micro- and nanomechanical resonator arrays

## I  Motivation: Non-linearity of MEMS and NEMS resonators

Recent technological advances have enabled the fabrication of mechanical resonators down to micrometer and even nanometer scales, with frequencies almost in reach of the GHz range.Craighead (2000); Roukes (2001a) Even though the main thrust in this field of research comes from the need to produce smaller, lighter, faster, and more efficient electromechanical systems, there is new basic physics to be learned along the way.Roukes (2001b) One particularly interesting aspect of the physical behavior of micro- and nanoelectromechanical systems (MEMS and NEMS) is their nonlinear mechanical response at relatively small deviations from equilibrium. This nonlinear behavior has been observed experimentally,Craighead (2000); Buks and Roukes (2001) and also exploited to achieve mechanical signal amplification and mechanical noise squeezing Rugar and Grütter (1991); Carr et al. (2000); Harrington and Roukes (2002) in single resonators. In addition, MEMS and NEMS facilitates the fabrication of large arrays of resonators, for which the coherent response might be useful for signal enhancement and noise reduction. It is important to understand the nonlinear behavior of MEMS and NEMS resonators in order to improve their future designs. At the same time, one can take advantage of these systems for the experimental study of nonlinear dynamics.

This paper is motivated by a recent experiment by Buks and Roukes (Buks and Roukes, 2002, henceforth BR) who fabricated an array of 67 fully-suspended doubly-clamped micromechanical resonating beams. Each beam was $`270\mu m\times 1\mu m\times 0.25\mu m`$ in size, and the distance between neighboring beams was $`4\mu m`$. The substrate beneath the array was completely etched away, forming a suspended diffraction grating with optical access from both sides. All even-numbered beams were electrically connected to one electrode and all odd-numbered beams to a second electrode. This allowed the application of electrostatic forces to induce coupling between the beams. The system was driven parametrically by introducing an AC component to the potential difference between the even-numbered and odd-numbered beams. The collective response of the array, as a function of the driving frequency and the DC component of the potential difference, was measured using optical diffraction. The response that BR inferred from their measurement was surprising in that (i) instead of showing a band consisting of a sequence of resonance peaks at the 67 normal frequencies of the array, the typical response as the frequency was swept up showed a small number of wide peaks where the response gradually increased and very abruptly decreased; and (ii) the array responded at frequencies beyond the expected top edge of the band.

We show below that both of these effects are a direct result of the fact that the restoring forces acting on the resonators as well as the damping that they undergo are both nonlinear. In Section II we describe the simplest equations of motion that are required to model the nonlinear resonator array. In Section III we solve the response of a single nonlinear resonator to parametric excitation at twice its resonance frequency using secular perturbation theory (for comparison, we solve in Appendix A the response to parametric excitation at the resonance frequency). In Section IV we use the same method to calculate the response of the coupled resonator array and obtain exact results for a few (two or three) resonators. Understanding the analytical results of these two sections allows us to interpret the results of Section V where we numerically integrate the equations of motion for an array of 67 resonators. Our results agree qualitatively with the observations of BR, explaining the two points mentioned above, but we suggest that further experiments be performed in order to test our theoretical calculations in a more quantitative fashion.

## II  Equations of motion

We seek the simplest set of equations of motion (EOM) that capture the important physical aspects of the array of coupled micromechanical beams. We first note that the normal frequencies of an individual beam are sufficiently separated that the frequency bands, formed by the coupling of the beams in the array, are well separated by gaps in which the system cannot respond. We therefore assume that we can treat the lowest band separately from all the others, so that each individual beam is oscillating strictly in its fundamental mode of vibration. Each beam can therefore be described by a single degree of freedom $`x_n`$, giving its displacement from equilibrium. We neglect any inhomogeneities in the fabrication of the beams and assume that all beams are identical. BR have actually examined each beam individually and report that their beams have a fairly uniform distribution of resonance frequencies, with an average of $`\omega _B=179.3`$kHz, and a standard deviation of 0.53kHz. There is a much larger variation in the quality factors of the beams prior to the application of electrostatic interaction between them, but this variation disappears when a small potential difference is introduced between the beams.Buks

The coordinates $`x_n`$ are all assumed small so that only terms to lowest order in $`x_n`$, necessary to capture the physical behavior of the system, will be kept in the EOM. Two types of forces act on the beams. The elastic restoring force of each beam and the electrostatic forces between the beams. Experiments done by Buks and Roukes Buks and Roukes (2001) on single beams of the type used in the array show that their response is like that of a Duffing oscillator—an oscillator whose restoring force contains a term proportional to the cube of the displacement and makes the oscillator stiffer than it would be within the harmonic approximation. Assuming a symmetric restoring force, and therefore no term proportional to an even power of $`x_n`$, and neglecting higher than cubic-order nonlinear corrections, the elastic force acting on the $`n^{th}`$ individual beam is

$$F_{elastic}^{(n)}=m\omega _B^2x_nm\alpha x_n^3,$$
(1)
where $`m`$ is the effective mass of a beam oscillating in its fundamental mode, whose frequency is $`\omega _B`$.

Even though the electrostatic force between two parallel charged wires decays only as $`1/r`$, for simplicity we consider only the attractive interactions between nearest-neighbor beams. Within this approximation each term in the EOM depends either on the variables $`x_n`$, describing the displacement of an individual beam from its equilibrium position, or on the difference variables $`x_{n+1}x_n`$, describing the relative displacements of a pair of neighboring beams. To keep the equations as simple as possible we restrict each type of nonlinear term in the EOM to depend either on $`x_n`$ or on $`x_{n+1}x_n`$, depending on whether it is mostly influenced by the elastic forces of the beams or the electrostatic interaction between them, respectively.

The cubic term in the expansion of the nearest-neighbor electrostatic interaction tends to pull the beams away from equilibrium, acting against the cubic term in the expansion of the elastic force in Eq. (1). Because, as we shall confirm later, the response of the array is consistent with having a cubic term which stiffens the beams, the elastic contribution to the cubic term is stronger than the electrostatic one. We therefore ignore the cubic as well as all higher terms in the electrostatic interaction, which we write as

$$F_{electric}^{(n)}=\frac{1}{2}m\mathrm{\Delta }^2[1+H\mathrm{cos}\omega _pt](x_{n+1}2x_n+x_{n1}).$$
(2)
Note that the linear electrostatic force constant $`\frac{1}{2}m\mathrm{\Delta }^2`$, which is modulated with a relative amplitude $`H1`$, representing the DC and the AC components of the applied voltage, is positive, acting to soften the elastic restoring force. The factor of $`1/2`$ is used with the difference variable for convenience.

Parametric excitation, as it appears in the bare Mathieu equation for a single oscillator of frequency $`\omega _0`$ without all the additional terms that we have here, is an instability of the system that occurs whenever the drive frequency is around one of the special values $`\omega _p=2\omega _0/n`$, where $`n`$ is an integer that labels the so-called instability tongues of the system (named after the tongue-shaped instability curves in the frequency-amplitude plane).Landau and Lifshitz (1976) We choose the parametric driving frequency $`\omega _p`$ to be around twice some value $`\omega _0`$ within the array’s band of normal frequencies. We are therefore exciting the system in its first ($`n=1`$) instability tongue. Thus,

$$\omega _p=2\omega _0+ϵ\mathrm{\Omega },$$
(3)
where $`ϵ`$ is a small parameter. In the BR experiment the system was actually excited in its second instability tongue, i.e. $`\omega _p`$ was chosen around some frequency in the band. It turns out that the response at the second tongue, apart from a few differences, is quite similar to that of the first tongue. We therefore prefer to carry out full calculations only for the first tongue which is somewhat easier to handle, and just for comparison, we calculate in Appendix A the response of a single nonlinear oscillator, excited at its second tongue.

There is good reason to believe that most of the dissipation in the coupled system is a result of the electrostatic interaction which causes currents to flow through the beams. This assumption is based on the observation of Buks and Roukes Buks  that the quality factors greatly diminish as the DC component of the electrostatic potential is increased. We therefore make the simplifying approximation that dissipation occurs predominantly as a result of currents, all other dissipation mechanisms being relatively negligible. This approximation relieves us of the worry about the variation in the quality factors of the individual beams before application of the electrostatic potential. The dissipative forces in the EOM are therefore written with respect to the difference variable,

$`F_{diss}^{(n)}`$ $`={\displaystyle \frac{1}{2}}m\omega _B\mathrm{\Gamma }(\dot{x}_{n+1}2\dot{x}_n+\dot{x}_{n1})`$ (4)
$`+{\displaystyle \frac{1}{2}}m\omega _B\alpha \eta [(x_{n+1}x_n)^2(\dot{x}_{n+1}\dot{x}_n)(x_nx_{n1})^2(\dot{x}_n\dot{x}_{n1})],`$
where we have included a nonlinear contribution to the dissipation, of the same order as the nonlinear elastic force (1). When putting all the pieces together, we (a) divide out the effective mass $`m`$ of a beam; (b) scale time $`tt/\omega _B`$ so that all frequencies (including $`\mathrm{\Delta }`$) are measured in units of $`\omega _B`$; and (c) scale length $`xx/\sqrt{\alpha }`$ to get rid of the dependence on $`\alpha `$. The equation of motion for the $`n^{th}`$ beam becomes

$`\ddot{x}_n`$ $`+`$ $`x_n+x_n^3+{\displaystyle \frac{1}{2}}\mathrm{\Delta }^2\left[1+H\mathrm{cos}(2\omega _0+ϵ\mathrm{\Omega })t\right](x_{n+1}2x_n+x_{n1})`$ (5)
$``$ $`{\displaystyle \frac{1}{2}}\mathrm{\Gamma }(\dot{x}_{n+1}2\dot{x}_n+\dot{x}_{n1})`$
$``$ $`{\displaystyle \frac{1}{2}}\eta \left[(x_{n+1}x_n)^2(\dot{x}_{n+1}\dot{x}_n)(x_nx_{n1})^2(\dot{x}_n\dot{x}_{n1})\right]=0.`$
In the following sections we shall solve these equations using secular perturbation theory. The physical parameter that allows us to use this approach is the linear damping coefficient which is assumed to be small. We therefore express it as $`\mathrm{\Gamma }=ϵ\gamma `$, taking $`ϵ`$ to be our small expansion parameter. The parametric instability of the system then occurs for small driving amplitude near resonance, and if, in addition, we consider the system near the onset of the instability, we can assume that the effects of nonlinearity are small as well. Thus, for small displacements $`x_n`$ all the nontrivial physical effects, namely, the parametric excitation, the cubic elastic restoring force, and both the linear and the amplitude-dependent dissipation, all enter the EOM as perturbative corrections to the simple linear equations. All these perturbative terms can be chosen to enter the EOM in the same order of the small parameter $`ϵ`$ by taking the leading order in $`x_n`$ to be of order $`ϵ^{1/2}`$, and expressing $`\mathrm{\Delta }^2H=ϵh`$. This ensures, as we shall confirm later on, that all the terms will contribute to the lowest-order solution we are seeking. The final form of the EOM is therefore

$`\ddot{x}_n`$ $`+`$ $`x_n+x_n^3+{\displaystyle \frac{1}{2}}\left[\mathrm{\Delta }^2+ϵh\mathrm{cos}(2\omega _0+ϵ\mathrm{\Omega })t\right](x_{n+1}2x_n+x_{n1})`$ (6)
$``$ $`{\displaystyle \frac{1}{2}}ϵ\gamma (\dot{x}_{n+1}2\dot{x}_n+\dot{x}_{n1})`$
$``$ $`{\displaystyle \frac{1}{2}}\eta \left[(x_{n+1}x_n)^2(\dot{x}_{n+1}\dot{x}_n)(x_nx_{n1})^2(\dot{x}_n\dot{x}_{n1})\right]=0.`$
As for boundary conditions, we follow the experiment of BR who had two additional fixed beams, identical to all the rest, at both ends of the array. This means that we define two extra variables and set them to zero, $`x_0=x_{N+1}=0`$.

## III  Response of a single parametrically-driven nonlinear oscillator

We begin by calculating the response of a single nonlinear oscillator to parametric excitation. Previous calculations of this problem exist in the literature(Litak et al., 1999, and references therein), nevertheless, we solve it here as a precursor to the many-oscillator case, treated in the next section. The equation of motion (6) for the single-oscillator case becomes

$`\ddot{x}`$ $`+`$ $`\left[\omega ^2ϵh\mathrm{cos}(2\omega +ϵ\mathrm{\Omega })t\right]x`$ (7)
$`+`$ $`ϵ\gamma \dot{x}+x^3+\eta x^2\dot{x}=0,`$
where we choose $`\omega _0`$ to be $`\omega =\sqrt{1\mathrm{\Delta }^2}`$, the resonance frequency of the beam in the harmonic approximation. The parametric excitation is performed around twice the actual resonance frequency of the oscillator. (In Appendix A we treat the case where the excitation is performed around the resonance frequency of the resonator.)

We calculate the correction to linear response by using secular perturbation theory.Hand and Finch (1998); Strogatz (1994) Recalling that the motion of the oscillator away from equilibrium is on the order of $`ϵ^{1/2}`$ we try a solution of the form

$$x(t)=ϵ^{1/2}(A(T)e^{i\omega t}+c.c.)+ϵ^{3/2}x_1(t)+\mathrm{}$$
(8)
where $`T=ϵt`$ is a slow time variable, allowing the complex amplitude $`A(T)`$ to vary slowly in time. The variation of $`A(T)`$ gives us the extra freedom to eliminate secular terms and ensure that the perturbative correction $`x_1(t)`$, as well as all higher-order corrections to the linear response, do not diverge (as they do if one uses naive perturbation theory). Using the relation

$$\dot{A}=\frac{dA}{dt}=ϵ\frac{dA}{dT}ϵA^{},$$
(9)
we calculate the time derivatives of the trial solution (8)

$`\dot{x}`$ $`=`$ $`ϵ^{1/2}([i\omega A+ϵA^{}]e^{i\omega t}+c.c.)+ϵ^{3/2}\dot{x}_1(t)+\mathrm{}`$ (10a)
$`\ddot{x}`$ $`=`$ $`ϵ^{1/2}([\omega ^2A+2i\omega ϵA^{}+ϵ^2A^{\prime \prime }]e^{i\omega t}+c.c.)+ϵ^{3/2}\ddot{x}_1(t)+\mathrm{}`$ (10b)
Substituting these expressions back into the equation of motion (7), and picking out all terms of order $`ϵ^{3/2}`$, we get the following equation for the first perturbative correction

$`\ddot{x}_1`$ $`+`$ $`\omega ^2x_1=(2i\omega A^{}e^{i\omega t}+c.c.)+h\mathrm{cos}[(2\omega +ϵ\mathrm{\Omega })t](Ae^{i\omega t}+c.c.)`$ (11)
$``$ $`\gamma (i\omega Ae^{i\omega t}+c.c.)(Ae^{i\omega t}+c.c.)^3\eta (Ae^{i\omega t}+c.c.)^2(i\omega Ae^{i\omega t}+c.c.).`$
The collection of terms proportional to $`e^{i\omega t}`$ on the right-hand side of Eq. (11), called the secular terms, act like a force, driving the simple harmonic oscillator on the left-hand side at its resonance frequency. The sum of all the secular terms must vanish so that the perturbative correction $`x_1(t)`$ will not diverge. This gives us an equation for determining the slowly varying amplitude $`A(T)`$. After expressing the cosine as a sum of exponentials we get

$$2i\omega \frac{dA}{dT}\frac{h}{2}A^{}e^{i\mathrm{\Omega }T}+i\omega \gamma A+3|A|^2A+i\omega \eta |A|^2A=0.$$
(12)
Ignoring initial transients, and assuming that the nonlinear terms in the equation are sufficient to saturate the growth of the instability, we try a steady-state solution of the form

$$A(T)=ae^{i\frac{\mathrm{\Omega }}{2}T}.$$
(13)
The solution to the equation of motion (7) is therefore

$$x(t)=ϵ^{1/2}(ae^{i(\omega +ϵ\mathrm{\Omega }/2)t}+c.c.)+O(ϵ^{3/2}),$$
(14)
where we are not interested in the correction $`x_1(t)`$ of order $`ϵ^{3/2}`$, but rather in the fixed amplitude $`a`$ of the lowest order term. This amplitude $`a`$ can be any solution of the equation

$$\left[(3|a|^2\omega \mathrm{\Omega })+i\omega (\gamma +\eta |a|^2)\right]a=\frac{h}{2}a^{},$$
(15)
obtained by substituting the steady-state solution (13) into the equation (12) of the secular terms. We immediately see that having no response ($`a=0`$) is always a possible solution regardless of the excitation frequency $`\mathrm{\Omega }`$. We divide both sides of the last equation by $`\gamma \omega `$ and define the rescaled variables $`\overline{a}=a/\sqrt{\gamma \omega }`$, $`\overline{\mathrm{\Omega }}=\mathrm{\Omega }/\gamma `$, $`\overline{\eta }=\omega \eta `$, and $`\overline{h}=h/2\gamma \omega `$, in terms of which the equation for the fixed complex amplitude $`a`$ becomes

$$\left[(3|\overline{a}|^2\overline{\mathrm{\Omega }})+i(1+\overline{\eta }|\overline{a}|^2)\right]\overline{a}=\overline{h}\overline{a}^{}.$$
(16)
Expressing $`\overline{a}=|\overline{a}|e^{i\varphi }`$ we obtain, after taking the magnitude squared of both sides, the intensity $`|\overline{a}|^2`$ of the non-trivial response as all positive roots of the equation

$$\left(3|\overline{a}|^2\overline{\mathrm{\Omega }}\right)^2+\left(1+\overline{\eta }|\overline{a}|^2\right)^2=\overline{h}^2.$$
(17)
This has the form of a distorted ellipse in the $`(\overline{\mathrm{\Omega }},|\overline{a}|^2)`$ plane, and a parabola in the $`(|\overline{a}|^2,\overline{h})`$ plane. In addition, we obtain for the relative phase of the response

$$\varphi =\frac{i}{2}\mathrm{ln}\frac{\overline{a}^{}}{\overline{a}}=\frac{1}{2}\mathrm{arctan}\frac{1+\overline{\eta }|\overline{a}|^2}{3|\overline{a}|^2\overline{\mathrm{\Omega }}}.$$
(18)
In Figure 1 we plot the response intensity $`|\overline{a}|^2`$ of a single oscillator to parametric excitation as a function of frequency $`\overline{\mathrm{\Omega }}`$, for fixed amplitude $`\overline{h}=1.5`$, in terms of the rescaled variables. Solid curves indicate stable solutions, and dashed curves are solutions that are unstable to small perturbations. Thin curves show the response without nonlinear damping ($`\overline{\eta }=0`$) which grows indefinitely with frequency $`\overline{\mathrm{\Omega }}`$ and is therefore incompatible with the experimental observations of BR and the assumptions of our calculation. Thick curves show the response with finite nonlinear damping ($`\overline{\eta }=1`$). With finite $`\overline{\eta }`$ there is a maximum value for the response $`|\overline{a}|_{max}^2=(\overline{h}1)/\overline{\eta }`$, and a maximum frequency,

$$\overline{\mathrm{\Omega }}_{SN}=\overline{h}\sqrt{1+\left(\frac{3}{\overline{\eta }}\right)^2}\frac{3}{\overline{\eta }},$$
(19)
at which the stability of the solution changes (known as a saddle-node bifurcation). For frequencies above $`\overline{\mathrm{\Omega }}_{SN}`$ the only solution is the trivial one $`\overline{a}=0`$. These values are indicated by horizontal and vertical dotted lines in Figure 1.

The threshold for the instability of the trivial solution is easily calculated by setting $`\overline{a}=0`$ in the expression (17) for the nontrivial solution. As seen in Figure 1, for a given $`\overline{h}`$ the threshold is situated at $`\overline{\mathrm{\Omega }}=\pm \sqrt{\overline{h}^21}`$. The threshold is plotted in Figure 2, in the $`(\mathrm{\Omega },h)`$ plane. Note that the minimal amplitude needed for instability is obtained on resonance ($`\overline{\mathrm{\Omega }}=0`$) and its value is $`\overline{h}=1`$, or $`h=2\gamma \omega `$, so that it scales as the linear damping coefficient $`\gamma `$.

Finally, in Figure 3 we plot the response intensity $`|\overline{a}|^2`$ of the oscillator as a function of amplitude $`\overline{h}`$, for fixed frequency $`\overline{\mathrm{\Omega }}`$ and finite nonlinear damping $`\overline{\eta }=1`$. Again, solid curves indicate stable solutions, and dashed curves unstable solutions. Thick curves show the response for $`\overline{\mathrm{\Omega }}=1`$, and thin curves show the response for $`\overline{\mathrm{\Omega }}=\overline{\eta }/3`$ and $`\overline{\mathrm{\Omega }}=1`$. The intersection of the trivial and the nontrivial solutions, which corresponds to the instability threshold, occurs at $`\overline{h}=\sqrt{\overline{\mathrm{\Omega }}^2+1}`$. For $`\overline{\mathrm{\Omega }}<\overline{\eta }/3`$ the nontrivial solution for $`|\overline{a}|^2`$ grows continuously for $`\overline{h}`$ above threshold and is stable. This is a supercritical bifurcation. On the other hand, for $`\overline{\mathrm{\Omega }}>\overline{\eta }/3`$ the bifurcation is subcritical—the nontrivial solution grows for $`\overline{h}`$ below threshold. This solution is unstable until the curve of $`|\overline{a}|^2`$ as a function of $`\overline{h}`$ bends around at a saddle-node bifurcation at

$$\overline{h}_{SN}=\frac{1+\frac{\overline{\eta }}{3}\overline{\mathrm{\Omega }}}{\sqrt{1+\left(\frac{\overline{\eta }}{3}\right)^2}},$$
(20)
where the solution becomes stable and $`|\overline{a}|^2`$ is once more an increasing function of $`\overline{h}`$. For amplitudes $`\overline{h}<\overline{h}_{SN}`$ the only solution is the trivial one $`\overline{a}=0`$.

Like the response of a forced Duffing oscillator, the response of a parametrically excited Duffing oscillator also exhibits hysteresis in a frequency scan. If the frequency $`\overline{\mathrm{\Omega }}`$ starts out at negative values and is increased gradually with a fixed amplitude $`\overline{h}`$, the response will gradually increase along the thick solid curve in Figure 1, until $`\overline{\mathrm{\Omega }}`$ reaches $`\overline{\mathrm{\Omega }}_{SN}`$ and the response drops abruptly to zero. If the frequency is then decreased gradually, the response will remain zero until $`\overline{\mathrm{\Omega }}`$ reaches the upper instability threshold $`\sqrt{\overline{h}^21}`$, and the response will jump abruptly to the thick solid curve above, and then gradually decrease to zero along this curve. A similar hysteretic behavior will be observed if the amplitude $`\overline{h}`$ is varied with a fixed frequency $`\overline{\mathrm{\Omega }}>\overline{\eta }/3`$, as can be inferred from Figure 3.

## IV  Response of a parametrically-driven array of nonlinear coupled oscillators—Secular perturbation theory

Consider now the coupled array of nonlinear oscillators as described by the general EOM (6). We calculate its response to parametric excitation, again using secular perturbation theory. We expand $`x_n(t)`$ as a sum of standing wave modes with slowly varying amplitudes

$$x_n(t)=ϵ^{1/2}\underset{m=1}{\overset{N}{}}(A_m(T)\mathrm{sin}(nq_m)e^{i\omega _mt}+c.c.)+ϵ^{3/2}x_n^{(1)}(t)+\mathrm{},n=1\mathrm{}N.$$
(21)
Recall that the boundary conditions are such that there are two additional fixed beams, labeled 0 and $`N+1`$, exerting electrostatic forces on the first and the last beams of the array. With these boundary conditions ($`x_0=x_{N+1}=0`$), the possible wave vectors $`q_m`$ are given by

$$q_m=\frac{m\pi }{N+1},m=1\mathrm{}N.$$
(22)
We substitute the trial solution (21) into the EOM term by term. Up to order $`ϵ^{3/2}`$ we have:

$$\ddot{x}_n=ϵ^{1/2}\underset{m}{}\mathrm{sin}(nq_m)([\omega ^2A_m+2i\omega ϵA_m^{}]e^{i\omega _mt}+c.c.)+ϵ^{3/2}\ddot{x}_n^{(1)}(t),$$
(23a)
$`x_{n+1}2x_n+x_{n1}`$ $`=`$ $`4ϵ^{1/2}{\displaystyle \underset{m}{}}\mathrm{sin}^2\left({\displaystyle \frac{q_m}{2}}\right)\mathrm{sin}(nq_m)(A_me^{i\omega _mt}+c.c.)`$ (23b)
$`+`$ $`ϵ^{3/2}\left(x_{n+1}^{(1)}2x_n^{(1)}+x_{n1}^{(1)}\right),`$
$$\frac{1}{2}ϵ\gamma (\dot{x}_{n+1}2\dot{x}_n+\dot{x}_{n1})=2ϵ^{3/2}\gamma \underset{m}{}\omega _m\mathrm{sin}^2\left(\frac{q_m}{2}\right)\mathrm{sin}(nq_m)(iA_me^{i\omega _mt}+c.c.),$$
(23c)
$`x_n^3`$ $`=`$ $`ϵ^{3/2}{\displaystyle \underset{j,k,l}{}}\mathrm{sin}(nq_j)\mathrm{sin}(nq_k)\mathrm{sin}(nq_l)`$ (23d)
$`\times `$ $`(A_je^{i\omega _jt}+c.c.)(A_ke^{i\omega _kt}+c.c.)(A_le^{i\omega _lt}+c.c.)`$
$`=`$ $`{\displaystyle \frac{ϵ^{3/2}}{4}}{\displaystyle \underset{j,k,l}{}}\{\mathrm{sin}[n(q_j+q_k+q_l)]+\mathrm{sin}[n(q_jq_k+q_l)]`$
$`+\mathrm{sin}[n(q_j+q_kq_l)]\mathrm{sin}[n(q_j+q_k+q_l)]\}`$
$`\times `$ $`\{A_jA_kA_le^{i(\omega _j+\omega _k+\omega _l)t}+3A_jA_kA_l^{}e^{i(\omega _j+\omega _k\omega _l)t}+c.c.\},`$
and
$`{\displaystyle \frac{1}{2}}\eta [(x_{n+1}x_n)^2(\dot{x}_{n+1}\dot{x}_n)`$ $``$ $`(x_nx_{n1})^2(\dot{x}_n\dot{x}_{n1})]`$ (23e)
$`=`$ $`2\eta ϵ^{3/2}{\displaystyle \underset{j,k,l}{}}\mathrm{sin}{\displaystyle \frac{q_j}{2}}\mathrm{sin}{\displaystyle \frac{q_k}{2}}\mathrm{sin}{\displaystyle \frac{q_l}{2}}`$
$`\times `$ $`\{\mathrm{sin}\left[{\displaystyle \frac{q_j+q_k+q_l}{2}}\right]\mathrm{sin}[n(q_j+q_k+q_l)]`$
$`+\mathrm{sin}[[{\displaystyle \frac{q_jq_k+q_l}{2}}]\mathrm{sin}[n(q_jq_k+q_l)]`$
$`+\mathrm{sin}[[{\displaystyle \frac{q_j+q_kq_l}{2}}]\mathrm{sin}[n(q_j+q_kq_l)]`$
$`+\mathrm{sin}\left[[{\displaystyle \frac{q_j+q_k+q_l}{2}}]\mathrm{sin}[n(q_j+q_k+q_l)]\right\}`$
$`\times `$ $`(A_je^{i\omega _jt}+c.c.)(A_ke^{i\omega _kt}+c.c.)(i\omega _lA_le^{i\omega _lt}+c.c.).`$
At order $`ϵ^{1/2}`$ we simply get the linear dispersion relation, given by

$$\omega _m^2=12\mathrm{\Delta }^2\mathrm{sin}^2(\frac{q_m}{2}),m=1\mathrm{}N.$$
(24)
At order $`ϵ^{3/2}`$ we get $`N`$ equations of the form

$$\ddot{x}_n^{(1)}+x_n^{(1)}+\frac{1}{2}\mathrm{\Delta }^2\left(x_{n+1}^{(1)}2x_n^{(1)}+x_{n1}^{(1)}\right)=\underset{m}{}\left(m^{th}\mathrm{secular}\mathrm{term}\right)e^{i\omega _mt}+\mathrm{other}\mathrm{terms},$$
(25)
where the left-hand sides are, again, coupled linear harmonic oscillators, with a dispersion relation given by (24). On the right-hand sides we have $`N`$ secular terms which act to drive the coupled oscillators $`x_n^{(1)}`$ at their resonance frequencies. As we did for a single oscillator in section III, here too we require that all the secular terms vanish so that the $`x_n^{(1)}`$ remain finite, and thus obtain equations for the slowly varying amplitudes $`A_m(T)`$. To extract the equation for the $`m^{th}`$ amplitude $`A_m(T)`$ we make use of the orthogonality of the modes, multiplying all the terms by $`\mathrm{sin}(nq_m)`$ and summing over $`n`$. We also express all normal frequencies relative to the same reference frequency $`\omega _0`$, used to define the excitation frequency $`\omega _p`$ in Eq. (3), so that

$`\omega _m=\omega _0+ϵ\mathrm{\Omega }_m.`$ (26)
We find that the coefficient of the $`m^{th}`$ secular term, which is required to vanish, is given by

$``$ $`2i\omega _m{\displaystyle \frac{dA_m}{dT}}2i\gamma \omega _m\mathrm{sin}^2\left({\displaystyle \frac{q_m}{2}}\right)+hA_m^{}\mathrm{sin}^2\left({\displaystyle \frac{q_m}{2}}\right)e^{i\left(\mathrm{\Omega }2\mathrm{\Omega }_m\right)T}`$ (27)
$``$ $`{\displaystyle \frac{3}{4}}{\displaystyle \underset{j,k,l}{}}A_jA_kA_l^{}e^{i\left(\mathrm{\Omega }_j+\mathrm{\Omega }_k\mathrm{\Omega }_l\mathrm{\Omega }_m\right)T}\mathrm{\Delta }_{jkl;m}^{(1)}`$
$``$ $`2\eta {\displaystyle \underset{j,k,l}{}}\{[2i\omega _lA_j^{}A_kA_le^{i\left(\mathrm{\Omega }_j+\mathrm{\Omega }_k+\mathrm{\Omega }_l\mathrm{\Omega }_m\right)T}i\omega _lA_jA_kA_l^{}e^{i\left(\mathrm{\Omega }_j+\mathrm{\Omega }_k\mathrm{\Omega }_l\mathrm{\Omega }_m\right)T}]`$
$`\times \mathrm{\Delta }_{jkl;m}^{(2)}\mathrm{sin}{\displaystyle \frac{q_j}{2}}\mathrm{sin}{\displaystyle \frac{q_k}{2}}\mathrm{sin}{\displaystyle \frac{q_l}{2}}\mathrm{sin}{\displaystyle \frac{q_m}{2}}\}=0,`$
where we have introduced two $`\mathrm{\Delta }`$ functions, defined in terms of Kronecker deltas as:

$`\mathrm{\Delta }_{jkl;m}^{(1)}`$ $`=`$ $`\delta _{j+k+l,m}\delta _{j+k+l,m}\delta _{j+k+l,2(N+1)m}`$ (28a)
$`+`$ $`\delta _{jk+l,m}\delta _{jk+l,m}\delta _{jk+l,2(N+1)m}`$
$`+`$ $`\delta _{j+kl,m}\delta _{j+kl,m}\delta _{j+kl,2(N+1)m}`$
$``$ $`\delta _{j+k+l,m}+\delta _{j+k+l,2(N+1)m}\delta _{j+k+l,2(N+1)+m},`$
and
$`\mathrm{\Delta }_{jkl;m}^{(2)}`$ $`=`$ $`\delta _{j+k+l,m}+\delta _{j+k+l,m}\delta _{j+k+l,2(N+1)m}`$ (28b)
$`+`$ $`\delta _{jk+l,m}+\delta _{jk+l,m}\delta _{jk+l,2(N+1)m}`$
$`+`$ $`\delta _{j+kl,m}+\delta _{j+kl,m}\delta _{j+kl,2(N+1)m}`$
$`+`$ $`\delta _{j+k+l,m}\delta _{j+k+l,2(N+1)m}\delta _{j+k+l,2(N+1)+m}.`$
These $`\mathrm{\Delta }`$ functions ensure the conservation of lattice momentum—the conservation of momentum to within the non-uniqueness of the specification of the normal modes due to the fact that $`\mathrm{sin}(nq_m)=\mathrm{sin}(nq_{2k(N+1)\pm m})`$ for any integer $`k`$. The first Kronecker delta in each line is a condition of direct momentum conservation, and the other two are the so-called umklapp conditions where only lattice momentum is conserved.

As for the single oscillator, we again try a steady-state solution, this time of the form

$$A_m(T)=a_me^{i\left(\frac{\mathrm{\Omega }}{2}\mathrm{\Omega }_m\right)T},$$
(29)
so that the solutions to the EOM, after substitution of (29) into (21), become

$$x_n(t)=ϵ^{1/2}\underset{m}{}(a_m\mathrm{sin}(nq_m)e^{i\left(\omega _0+\frac{ϵ\mathrm{\Omega }}{2}\right)t}+c.c.)+O(ϵ^{3/2}),$$
(30)
where all modes are oscillating at half the parametric excitation frequency $`\omega _p/2`$.

As before, we are not interested in the corrections of order $`ϵ^{3/2}`$ but only in the values of the fixed amplitudes $`a_m`$ as functions of all the parameters of the original EOM. Substituting the steady state solution (29) into the equations (27) for the time-varying amplitudes $`A_m(T)`$, we obtain the required equations for the fixed complex amplitudes $`a_m`$

$`(\mathrm{\Omega }2\mathrm{\Omega }_m)\omega _ma_m2i\gamma \omega _ma_m\mathrm{sin}^2\left({\displaystyle \frac{q_m}{2}}\right)+ha_m^{}\mathrm{sin}^2\left({\displaystyle \frac{q_m}{2}}\right){\displaystyle \frac{3}{4}}{\displaystyle \underset{j,k,l}{}}a_ja_ka_l^{}\mathrm{\Delta }_{jkl;m}^{(1)}`$
$`2i\eta \mathrm{sin}{\displaystyle \frac{q_m}{2}}{\displaystyle \underset{j,k,l}{}}\omega _l\left[2a_j^{}a_ka_la_ja_ka_l^{}\right]\mathrm{sin}{\displaystyle \frac{q_j}{2}}\mathrm{sin}{\displaystyle \frac{q_k}{2}}\mathrm{sin}{\displaystyle \frac{q_l}{2}}\mathrm{\Delta }_{jkl;m}^{(2)}=0.`$ (31)
We can change to rescaled variables as we did in the case of a single oscillator by dividing the equations for the amplitudes (IV) by $`(\gamma \omega _0)^{3/2}`$ and defining as before $`\overline{a}_j=a_j/\sqrt{\gamma \omega _0}`$, $`\overline{\mathrm{\Omega }}=\mathrm{\Omega }/\gamma `$, $`\overline{\eta }=\omega _0\eta `$, and $`\overline{h}=h/2\gamma \omega _0`$, and in addition $`r_m=\omega _m/\omega _0`$ and $`\delta _m=2\mathrm{\Omega }_m/\gamma `$. After doing so we obtain the rescaled equations

$`(\overline{\mathrm{\Omega }}\delta _m)r_m\overline{a}_m2ir_m\mathrm{sin}^2\left({\displaystyle \frac{q_m}{2}}\right)\overline{a}_m+2\overline{h}\mathrm{sin}^2\left({\displaystyle \frac{q_m}{2}}\right)\overline{a}_m^{}{\displaystyle \frac{3}{4}}{\displaystyle \underset{j,k,l}{}}\overline{a}_j\overline{a}_k\overline{a}_l^{}\mathrm{\Delta }_{jkl;m}^{(1)}`$
$`2i\overline{\eta }\mathrm{sin}{\displaystyle \frac{q_m}{2}}{\displaystyle \underset{j,k,l}{}}r_l\left[2\overline{a}_j^{}\overline{a}_k\overline{a}_l\overline{a}_j\overline{a}_k\overline{a}_l^{}\right]\mathrm{sin}{\displaystyle \frac{q_j}{2}}\mathrm{sin}{\displaystyle \frac{q_k}{2}}\mathrm{sin}{\displaystyle \frac{q_l}{2}}\mathrm{\Delta }_{jkl;m}^{(2)}=0.`$ (32)
This is the main result of the perturbative calculation. We have managed to replace $`N`$ coupled differential equations (6) for the oscillator coordinates $`x_n(t)`$ by $`N`$ coupled algebraic equations (IV) for the time-independent mode amplitudes $`a_m`$. All that remains, in order to obtain the overall collective response of the array as a function of the parameters of the original EOM, is to solve these coupled algebraic equations.

Before doing so we should note the following general statements. First, one can easily verify that for a single oscillator ($`N=j=k=l=m=1`$), the general equation (IV) reduces to the single-oscillator equation (15), we derived in section III. Next, one can also see that the trivial solution, $`a_m=0`$ for all $`m`$, always satisfies the equations, though, as we have seen in the case of a single oscillator, it is not always a stable solution. Finally, one can also verify that whenever for a given $`m`$, $`\mathrm{\Delta }_{mmm;j}^{(1)}=\mathrm{\Delta }_{mmm;j}^{(2)}=0`$ for all $`jm`$, then a single-mode solution exists with $`a_m0`$ and $`a_j=0`$ for all $`jm`$. These single-mode solutions have the elliptical shape of the single-oscillator solution given in Eq. (17), and satisfy the equation

$$\frac{1}{4\mathrm{sin}^4(q_m/2)}\left(\frac{3}{4}\mathrm{\Delta }_{mmm;m}^{(1)}|\overline{a}_m|^2\overline{\mathrm{\Omega }}\right)^2+\left(1+\mathrm{sin}^2\frac{q_m}{2}\mathrm{\Delta }_{mmm;m}^{(2)}\overline{\eta }|\overline{a}_m|^2\right)^2=\overline{h}^2,$$
(33)
where for each solution we have set $`\omega _0=\omega _m`$, so that $`\delta _m=0`$ and $`r_m=1`$. Note that generically $`\mathrm{\Delta }_{mmm;m}^{(1)}=\mathrm{\Delta }_{mmm;m}^{(2)}=3`$, except when umklapp conditions are satisfied.

Additional solutions, involving more than a single mode, exist in general but are hard to obtain analytically. We calculate these multi-mode solutions below for the case of two and three oscillators by finding the roots of the coupled algebraic equations numerically. In Appendix B we show the explicit sets of coupled mode-amplitude equations for these cases.

In Fig. 4 we show the solutions for the response intensity of two oscillators as a function of frequency, for a particular choice of the equation parameters. The top graph shows the square of the amplitude of the antisymmetric mode $`\overline{a}_2`$, whereas the bottom graph shows the square of the amplitude of the symmetric mode $`\overline{a}_1`$. Solid curves indicate stable solutions and dashed curves indicate unstable solutions. The two elliptical single-mode solution branches, mentioned in the previous paragraph are easily spotted. These branches are labeled by $`S_1`$ and $`S_2`$ (In Appendix B, Equations (49), we give the analytical expressions for these two solution branches). In addition, we find two double-mode solution branches, labeled $`D_1`$ and $`D_2`$, involving the excitation of both modes simultaneously. Note that the two branches of double-mode solutions intersect at a point where they switch their stability.

With two oscillators we obtain regions in frequency where three stable solutions can exist. If all the stable solution branches are accessible experimentally then the observed effects of hysteresis might be more complex than in the simple case of a single oscillator. This is demonstrated in Fig. 5 where we compare our analytical solutions with a numerical integration of the differential equations of motion (6) for two oscillators. The response intensity, plotted here, is the time and space averages of the square of the oscillator displacements

$$I=\frac{1}{N}\underset{n=1}{\overset{N}{}}x_n^2,$$
(34)
where the angular brackets denote time average, and here $`N=2`$. A solid curve shows the response intensity for frequency swept upwards, and a dashed curve shows the response intensity for frequency swept downwards. Small circles show the analytical response intensity, using the fact that $`I=3(|\overline{a}_1|^2+|\overline{a}_2|^2)/2`$, for the stable regions of the four solution branches shown in Fig. 4. With the analytical solution in the background, one can easily understand all the discontinuous jumps, as well as the hysteresis effects, that are obtained in the numerical solution of the equations of motion. Note the the $`S_1`$ branch is missed in the upwards frequency sweep and is only accessed by the system in the downwards frequency sweep. One could trace the whole stable region of the $`S_1`$ branch by changing the sweep direction after jumping onto the branch, thereby climbing all the way up to the end of the $`S_1`$ branch and then falling onto the tip of the $`D_1`$ branch or to zero.

In Fig. 6 we show the solutions for the response intensity of three oscillators as a function of frequency, for a particular choice of the equation parameters. The graphs show the squares of the amplitudes of the three different modes. Solid curves indicate stable solutions and dashed curves indicate unstable ones. For three oscillators there is only one elliptical single-mode solution branch, of the form of Eq. (33), whose exact analytical expression is given in Eq. (53). This branch is labeled by $`S_2`$. In addition, we find a host of nontrivial multi-mode solution branches, including one that is disconnected from all other branches. We show these plots, not only to demonstrate that it is possible to obtain such solutions exactly, but also to emphasize the large number and nontrivial structure of the solution branches one finds, even for such a small number of oscillators. This can only serve as a hint for the multi-mode solutions one can expect to find when the number of oscillators is large, as in the BR experiment.

## V  Response of parametrically-driven nonlinear coupled oscillators—Numerical integration of the equations

The equations of motion (6) were integrated numerically for an array of $`67`$ oscillators, as in the BR experiment. The results for the response intensity (34) as a function of parametric drive frequency $`\omega _p`$ (measured in units of the top band edge frequency $`\omega _0`$) are shown in Figure 7. These results must be considered illustrative only, since many of the parameters of the experimental system are not known. The parameters used to construct the figure, $`\mathrm{\Delta }^2=0.02,h=0.016,\gamma =0.004,\eta =6.0`$, were chosen using the insights gained from the two and three beam results. We should emphasize that the structure of the response branches depends very strongly on the equation parameters.

A number of the important features of the experimental data are reproduced by these calculations. We concentrate on the solid curve in the figure, which is for frequency swept upwards, since this is the protocol that was used in the experiment. In particular, the response intensity shows features that span a range of frequencies that is large compared with the mode spacing (which is about 0.0006 for the parameter used). The lowest frequency feature, from about $`\omega _p/\omega _0=1.94`$ to $`\omega _p/\omega _0=1.97`$, can be identified as the response to the parametric drive of a single mode at or very near the band edge at $`\omega /\omega _0=0.98`$, analogous to the one mode response shown in Figure 1. Furthermore, the variation of the response with frequency shows abrupt jumps, particularly on the high frequency side of the features as the frequency is raised. Finally, the response extends to frequencies higher than the band edge for the linear modes, which would give a response only up to $`\omega _p/\omega _0=2.0`$. All of these features are understood now that we have seen the analytical solutions for small numbers of oscillators. In particular, the wide features compared with the mode spacing are explained by the simple fact the as the frequency is swept upwards a particular solution branch is followed as long as it remains stable. In the meantime many other stable solutions, that may be as close to each other as the mode spacing, are simply skipped over.

Comparing the two traces in Figure 7 shows that the response for a downward frequency sweep is significantly different, with a less dramatic variation of the response. In particular, note that the downwards sweep was able to access additional stable solution branches that were missed in the upwards sweep. There is also no response above $`\omega _p/\omega _0=2.0`$ in this case. This is because the zero displacement state is stable for $`\omega _p/\omega _0>2.0`$, and the system will remain in this state as the frequency is lowered, unless a large enough disturbance kicks it onto another of the solution branches. The hysteresis on reversing the frequency sweep was not looked at in the first experiments, and it would be interesting to test this prediction in further experiments.

## VI  conclusions

We have calculated the response of nonlinear coupled oscillators to parametric excitation. Our calculations agree qualitatively with the experimental measurements of Buks and Roukes Buks and Roukes (2002) and explain the main features observed in the experiment. The abrupt drops in the response as the frequency is swept upwards, the response continuing beyond the upper edge of the frequency band, and the large size of the response features compared with the mode spacing, are all qualitatively explained.

Nevertheless, we propose that a more systematic study be conducted on systems of coupled nonlinear resonators so that our theoretical predictions could be tested more quantitatively. For example, successive measurements on systems containing only one, two, and three coupled resonators, that are made as identically as possible, could be used to extract the nonlinear parameters of the resonators. These could then be used to predict and explain the response of a large resonator array more quantitatively.

Furthermore, we have demonstrated that as the number of oscillators is increased the number of the solution branches for the response of the system increases and the effects of hysteresis become more and more complicated. This suggests that the appropriate experimental protocol for studying a system with many oscillators should be—in addition to the standard up-sweep and down-sweep in frequency—to change the direction of the sweep after every abrupt change in the response intensity. Such a protocol may provide more information about the response curve by accessing additional branches of the solution and fully tracing them out.

## Appendix A  Parametric excitation of a single oscillator at its second instability tongue

For a single nonlinear oscillator, like the one studied in section III, which is parametrically-excited at its second instability tongue, Eq. (5) becomes

$`\ddot{x}`$ $`+`$ $`\left[\omega ^2\mathrm{\Delta }^2H\mathrm{cos}(\omega +ϵ\mathrm{\Omega })t\right]x`$ (35)
$`+`$ $`ϵ\gamma \dot{x}+x^3+\eta x^2\dot{x}=0,`$
where again $`\omega =\sqrt{1\mathrm{\Delta }^2}`$ is the resonance frequency of the oscillator in the harmonic approximation, but the parametric excitation is performed around $`\omega `$ and not around $`2\omega `$. In this case the scaling of $`\mathrm{\Delta }^2H`$ with respect to $`ϵ`$ needs to be redetermined. The technical reason for this is that if we naively take $`\mathrm{\Delta }^2H=ϵh`$, as before, then the parametric driving term does not contribute to the order $`ϵ^{3/2}`$ secular term which we use to find the response, and the order $`ϵ^{1/2}`$ term in $`x`$ becomes identically zero.

The remedy for this situation is to let $`\mathrm{\Delta }^2H`$ scale like $`ϵ^p`$, with $`p<1`$, so that there will be a non-secular correction to $`x`$ at a lower order than $`ϵ^{3/2}`$. The value of $`p`$ will be chosen such that this correction will contribute to the order $`ϵ^{3/2}`$ secular term and will give us the required response. The equation of motion (35) becomes

$$\ddot{x}+\omega ^2x=\frac{hϵ^p}{2}(e^{i(\omega t+\mathrm{\Omega }T)}+c.c.)xϵ\gamma \dot{x}x^3\eta x^2\dot{x},$$
(36)
and we try an expansion of the solution of the form

$$x(t)=ϵ^{1/2}(A(T)e^{i\omega t}+c.c.)+ϵ^{p+1/2}x_p(t)+ϵ^{3/2}x_1(t)+\mathrm{}$$
(37)
Substituting this expansion into the equation of motion (36) we obtain at order $`ϵ^{1/2}`$ the linear equation as usual, and at order $`ϵ^{p+1/2}`$

$$\ddot{x}_p+\omega ^2x_p=\frac{h}{2}(Ae^{i(2\omega t+\mathrm{\Omega }T)}+A^{}e^{i\mathrm{\Omega }T}+c.c.).$$
(38)
As expected, there is no secular term on the right-hand side, so we can immediately solve for $`x_p`$,

$$x_p(t)=\frac{h}{2}(\frac{A}{3\omega ^2}e^{i(2\omega t+\mathrm{\Omega }T)}+\frac{A^{}}{\omega ^2}e^{i\mathrm{\Omega }T}+c.c.)+O(ϵ).$$
(39)
Substituting the solution for $`x_p`$ into the expansion (37), and the expansion back into the equation of motion (36), contributes an additional term from the parametric driving which has the form

$`ϵ^{2p+1/2}{\displaystyle \frac{h^2}{4}}({\displaystyle \frac{A}{3\omega ^2}}e^{i(2\omega t+\mathrm{\Omega }T)}+{\displaystyle \frac{A^{}}{\omega ^2}}e^{i\mathrm{\Omega }T}+c.c.)(e^{i(\omega t+\mathrm{\Omega }T)}+c.c.)`$
$`=ϵ^{2p+1/2}{\displaystyle \frac{h^2}{4\omega ^2}}\left({\displaystyle \frac{2}{3}}A+A^{}e^{i2\mathrm{\Omega }T}\right)e^{i\omega t}+c.c.+\mathrm{non}\mathrm{secular}\mathrm{terms}.`$ (40)
To contribute to the order $`ϵ^{3/2}`$ secular term, we see that we must set $`p=1/2`$. This gives us the required contribution to the equation for the vanishing secular terms. All other terms remain as they were in Eq. (12), so that the new equation for determining $`A(T)`$ becomes

$$2i\omega \frac{dA}{dT}\frac{h^2}{4\omega ^2}\left(\frac{2}{3}A+A^{}e^{i2\mathrm{\Omega }T}\right)+i\omega \gamma A+3|A|^2A+i\omega \eta |A|^2A=0.$$
(41)
Again, ignoring initial transients, and assuming that the nonlinear terms in the equation are sufficient to saturate the growth of the instability, we try a steady-state solution, this time of the form

$$A(T)=ae^{i\mathrm{\Omega }T}.$$
(42)
The solution to the equation of motion (7) is therefore

$$x(t)=ϵ^{1/2}(ae^{i(\omega +ϵ\mathrm{\Omega })t}+c.c.)+O(ϵ),$$
(43)
where the correction $`x_{1/2}`$ of order $`ϵ`$ is given in Eq. (39) and, as before, we are not interested in the correction $`x_1(t)`$ of order $`ϵ^{3/2}`$, but rather in the fixed amplitude $`a`$ of the lowest order term. We substitute the steady-state solution (42) into the equation (41) of the secular terms and obtain

$$\left[\left(3|a|^22\omega \mathrm{\Omega }\frac{2}{3}\frac{h^2}{4\omega ^2}\right)+i\omega \left(\gamma +\eta |a|^2\right)\right]a=\frac{h^2}{4\omega ^2}a^{}.$$
(44)
We divide both sides of the last equation by $`\gamma \omega `$ and define the rescaled variables: $`\overline{a}=a/\sqrt{\gamma \omega }`$, $`\overline{\mathrm{\Omega }}=\mathrm{\Omega }/\gamma `$, $`\overline{\eta }=\omega \eta `$, and $`\overline{h}=h/2\sqrt{\gamma \omega ^3}`$, in terms of which we obtain after taking the magnitude squared of both sides, in addition to the trivial solution $`a=0`$, the non-trivial response

$$\left(3|\overline{a}|^22\overline{\mathrm{\Omega }}\frac{2}{3}\overline{h}^2\right)^2+\left(1+\overline{\eta }|\overline{a}|^2\right)^2=\overline{h}^4.$$
(45)
Fig. 8 shows the response intensity $`|\overline{a}|^2`$ as a function of the frequency $`\overline{\mathrm{\Omega }}`$, for fixed amplitude $`\overline{h}=1.5`$ in the second instability tongue. The solution looks very similar to the response shown in Fig. 1 for the first instability tongue, though we should point out two important differences. The first is that the orientation of the ellipse, indicated by the slope of the curves for $`\overline{\eta }=0`$, is different. The slope here is $`2/3`$ whereas for the first instability tongue the slope is $`1/3`$. The second is the change in the definition of $`\overline{h}`$. The lowest amplitude required for having an instability, is again on resonance ($`\overline{\mathrm{\Omega }}=0`$), and its value is again $`\overline{h}=1`$, but now this implies that $`h=2\sqrt{\gamma \omega ^3}`$, or that $`h`$ it scales as $`\sqrt{\gamma }`$. This is consistent with the well known result (see, for example, Ref. Landau and Lifshitz, 1976) that the minimal amplitude for the instability of the $`n^{th}`$ tongue scales as $`\gamma ^{1/n}`$.

## Appendix B  Explicit equations for two and three coupled oscillators

### B.1 Two coupled nonlinear oscillators

For two coupled oscillators ($`N=2`$) we have

$`q_1={\displaystyle \frac{\pi }{3}},q_2={\displaystyle \frac{2\pi }{3}},`$ (46)
$`\omega _1^2=1{\displaystyle \frac{1}{2}}\mathrm{\Delta }^2,\omega _2^2=1{\displaystyle \frac{3}{2}}\mathrm{\Delta }^2,`$ (47)
and we choose the reference frequency $`\omega _0`$ to be $`\omega _2`$, so that $`\delta _2=\mathrm{\Omega }_2=0`$, $`r_2=1`$, $`\delta _1=2(\omega _1\omega _2)/\gamma ϵ\delta >0`$, and $`r_1=\omega _1/\omega _2r`$. For $`\mathrm{\Delta }1`$, $`\delta \mathrm{\Delta }^2/ϵ\gamma `$ and $`r1+\mathrm{\Delta }^2/2`$. The first mode is the symmetric one with $`x_1(t)=x_2(t)`$ and the second mode is antisymmetric with $`x_1(t)=x_2(t)`$. The equations (IV) for the rescaled complex amplitudes $`\overline{a}_1`$ and $`\overline{a}_2`$ are

$`(\overline{\mathrm{\Omega }}\delta )r\overline{a}_1i{\displaystyle \frac{r}{2}}\overline{a}_1+{\displaystyle \frac{\overline{h}}{2}}\overline{a}_1^{}{\displaystyle \frac{9}{4}}\left(|\overline{a}_1|^2\overline{a}_1+2|\overline{a}_2|^2\overline{a}_1+\overline{a}_2^2\overline{a}_1^{}\right)`$
$`{\displaystyle \frac{3}{8}}i\overline{\eta }\left[r|\overline{a}_1|^2\overline{a}_1+2r|\overline{a}_2|^2\overline{a}_1+(2r)\overline{a}_2^2\overline{a}_1^{}\right]=0,`$ (48a)
$`\overline{\mathrm{\Omega }}\overline{a}_2i{\displaystyle \frac{3}{2}}\overline{a}_2+{\displaystyle \frac{3}{2}}\overline{h}\overline{a}_2^{}{\displaystyle \frac{9}{4}}\left(|\overline{a}_2|^2\overline{a}_2+2|\overline{a}_1|^2\overline{a}_2+\overline{a}_1^2\overline{a}_2^{}\right)`$
$`{\displaystyle \frac{3}{8}}i\overline{\eta }\left[9|\overline{a}_2|^2\overline{a}_2+2|\overline{a}_1|^2\overline{a}_2+(2r1)\overline{a}_1^2\overline{a}_2^{}\right]=0.`$ (48b)
The two single-mode solution branches, having the general form of Eq. (33) and labeled $`S_1`$ and $`S_2`$ in Fig. 4, are easily obtained by setting $`\overline{a}_2`$ or $`\overline{a}_1`$ to zero in the coupled equations above, respectively. This yields the analytical forms of these solutions which are

$`S_1:`$ $`\left({\displaystyle \frac{9}{2}}|\overline{a}_1|^22r(\overline{\mathrm{\Omega }}\delta )\right)^2+r^2\left(1+{\displaystyle \frac{3}{4}}\overline{\eta }|\overline{a}_1|^2\right)^2=\overline{h}^2,`$ (49a)
$`S_2:`$ $`\left({\displaystyle \frac{3}{2}}|\overline{a}_2|^2{\displaystyle \frac{2}{3}}\overline{\mathrm{\Omega }}\right)^2+\left(1+{\displaystyle \frac{9}{4}}\overline{\eta }|\overline{a}_2|^2\right)^2=\overline{h}^2.`$ (49b)
### B.2 Three coupled nonlinear oscillators

For three coupled oscillators ($`N=3`$) we have

$`q_1={\displaystyle \frac{\pi }{4}},q_2={\displaystyle \frac{\pi }{2}},q_3={\displaystyle \frac{3\pi }{4}}`$ (50)
$`\omega _1^2=1\mathrm{\Delta }^2+\mathrm{\Delta }^2/\sqrt{2},\omega _2^2=1\mathrm{\Delta }^2,\omega _3^2=1\mathrm{\Delta }^2\mathrm{\Delta }^2/\sqrt{2},`$ (51)
and we choose the reference frequency $`\omega _0`$ to be $`\omega _2`$, so that $`\delta _2=0`$, $`r_2=1`$, and $`\delta _1=2(\omega _1\omega _2)/\gamma ϵ=\delta _3\delta >0`$. For $`\mathrm{\Delta }1`$, $`\delta \mathrm{\Delta }^2/\sqrt{2}ϵ\gamma `$ and $`r_{1,3}1\pm \mathrm{\Delta }^2/2\sqrt{2}`$. The equations (IV) for the rescaled complex amplitudes $`\overline{a}_1`$, $`\overline{a}_2`$, and $`\overline{a}_3`$ are

$`(\overline{\mathrm{\Omega }}\delta )r_1\overline{a}_1i{\displaystyle \frac{2\sqrt{2}}{2}}r_1\overline{a}_1+{\displaystyle \frac{2\sqrt{2}}{2}}\overline{h}\overline{a}_1^{}`$
$`{\displaystyle \frac{3}{4}}\left(3|\overline{a}_1|^2\overline{a}_1|\overline{a}_3|^2\overline{a}_3+4|\overline{a}_2|^2\overline{a}_3+2\overline{a}_2^2\overline{a}_3^{}+4|\overline{a}_2|^2\overline{a}_1+2\overline{a}_2^2\overline{a}_1^{}+6|\overline{a}_3|^2\overline{a}_1+3\overline{a}_3^2\overline{a}_1^{}2|\overline{a}_1|^2\overline{a}_3\overline{a}_1^2\overline{a}_3^{}\right)`$
$`i\overline{\eta }[{\displaystyle \frac{32\sqrt{2}}{4}}3r_1|\overline{a}_1|^2\overline{a}_1{\displaystyle \frac{\sqrt{2}+1}{4}}r_3|\overline{a}_3|^2\overline{a}_3+{\displaystyle \frac{2\sqrt{2}}{2}}(2r_1|\overline{a}_2|^2\overline{a}_1+(2r_1)\overline{a}_2^2\overline{a}_1^{})`$
$`+{\displaystyle \frac{\sqrt{2}1}{4}}(2r_3|\overline{a}_1|^2\overline{a}_3+(2r_1r_3)\overline{a}_1^2\overline{a}_3^{})+{\displaystyle \frac{1}{4}}(2r_1|\overline{a}_3|^2\overline{a}_1+(2r_3r_1)\overline{a}_3^2\overline{a}_1^{})]=0,`$ (52a)
$`\overline{\mathrm{\Omega }}\overline{a}_2i\overline{a}_2+\overline{h}\overline{a}_2^{}{\displaystyle \frac{3}{2}}\left(2|\overline{a}_2|^2\overline{a}_2+2|\overline{a}_1|^2\overline{a}_2+\overline{a}_1^2\overline{a}_2^{}+2|\overline{a}_3|^2\overline{a}_2+\overline{a}_3^2\overline{a}_2^{}+\overline{a}_1^{}\overline{a}_2\overline{a}_3+\overline{a}_1\overline{a}_2^{}\overline{a}_3+\overline{a}_1\overline{a}_2\overline{a}_3^{}\right)`$
$`i\overline{\eta }[|\overline{a}_2|^2\overline{a}_2+{\displaystyle \frac{2\sqrt{2}}{2}}(2|\overline{a}_1|^2\overline{a}_2+(2r_11)\overline{a}_1^2\overline{a}_2^{})`$
$`+{\displaystyle \frac{2+\sqrt{2}}{2}}(2|\overline{a}_3|^2\overline{a}_2+(2r_31)\overline{a}_3^2\overline{a}_2^{})]=0,`$ (52b)
$`(\overline{\mathrm{\Omega }}+\delta )r_3\overline{a}_3i{\displaystyle \frac{2+\sqrt{2}}{2}}r_3\overline{a}_3+{\displaystyle \frac{2+\sqrt{2}}{2}}\overline{h}\overline{a}_3^{}`$
$`{\displaystyle \frac{3}{4}}\left(3|\overline{a}_3|^2\overline{a}_3|\overline{a}_1|^2\overline{a}_1+4|\overline{a}_2|^2\overline{a}_3+2\overline{a}_2^2\overline{a}_3^{}+4|\overline{a}_2|^2\overline{a}_1+2\overline{a}_2^2\overline{a}_1^{}+6|\overline{a}_1|^2\overline{a}_3+3\overline{a}_1^2\overline{a}_3^{}2|\overline{a}_3|^2\overline{a}_1\overline{a}_3^2\overline{a}_1^{}\right)`$
$`i\overline{\eta }[{\displaystyle \frac{3+2\sqrt{2}}{4}}3r_3|\overline{a}_3|^2\overline{a}_3+{\displaystyle \frac{\sqrt{2}1}{4}}r_1|\overline{a}_1|^2\overline{a}_1+{\displaystyle \frac{2+\sqrt{2}}{2}}(2r_3|\overline{a}_2|^2\overline{a}_3+(2r_3)\overline{a}_2^2\overline{a}_3^{})`$
$`{\displaystyle \frac{\sqrt{2}+1}{4}}(2r_1|\overline{a}_3|^2\overline{a}_1+(2r_3r_1)\overline{a}_3^2\overline{a}_1^{})+{\displaystyle \frac{1}{4}}(2r_3|\overline{a}_1|^2\overline{a}_3+(2r_1r_3)\overline{a}_1^2\overline{a}_3^{})]=0.`$ (52c)
Only one single-mode solution of the form of Eq. (33) exists in the case of three oscillators, and involves the second mode. It is obtained by setting $`\overline{a}_1=\overline{a}_3=0`$ in the coupled equations above. The analytical expression for this solution is

$$S_2:(3|\overline{a}_2|^2\overline{\mathrm{\Omega }})^2+(1+\overline{\eta }|\overline{a}_2|^2)^2=\overline{h}^2.$$
(53)
###### Acknowledgements.

This work was supported by Grant No. 1999458 from the U.S.-Israel Binational Science Foundation (BSF), and NSF Grant No. DMR-9873573. We thank Eyal Buks, Michael Roukes, and Inna Kozinsky for many useful discussions and for sharing their experimental results prior to publication.