text:
# First-principles calculations of the self-trapped exciton in crystalline NaCl

## Abstract

The atomic and electronic structure of the lowest triplet state of the off-center (C$`_{\text{2v}}`$ symmetry) self-trapped exciton (STE) in crystalline NaCl is calculated using the local-spin-density (LSDA) approximation. In addition, the Franck-Condon broadening of the luminescence peak and the a$`{}_{1g}{}^{}\text{b}_{3u}`$ absorption peak are calculated and compared to experiment. LSDA accurately predicts transition energies if the initial and final states are both localized or delocalized, but 1 eV discrepancies with experiment occur if one state is localized and the other is delocalized.

Unlike a molecule, an extended system such as a solid can support both spatially localized and delocalized single particle states and excitations. Physical properties such as luminescence can differ dramatically depending on which type of state occurs. Deciding theoretically whether a localized or delocalized solution exists is a challenging problem . Here we examine NaCl, a classic example  where electronic excitations self-localize by coupling to the lattice, creating local lattice distortions. Because the degree of localization will affect the Coulomb energies, approaches that incompletely cancel the self-interaction contribution to the exchange energy (e.g., the local density approximation) sometimes fail to predict the actual localized solution.

The ground state of alkali halides with one electron removed is the VK center : the resulting hole does not delocalize at the top of the valence band, but rather (symmetrically) attracts two Cl- ions into a tightly bound molecule , effectively becoming a Cl$`{}_{}{}^{}{}_{2}{}^{}`$ molecular ion. The local symmetry of this atomic configuration is D$`_{\text{2h}}`$. An excess electron in bulk NaCl forms a large mobile (Fröhlich) polaron , but in the presence of a self-trapped hole forms a self-trapped exciton (STE). It has been suggested  that the self-trapped exciton state breaks symmetry and sits off-center with C$`_{\text{2v}}`$ symmetry.

Although there have been many theoretical studies  of the self-trapped exciton and VK center problems, no density-functional calculations have been reported. In this letter we report local-spin-density approximation (LSDA) calculations of the STE and VK center. The advantage of LSDA calculations is that they provide one of the simplest tools capable of providing a realistic model of this competition. In both cases, our LSDA calculations give undistorted, delocalized solutions with lower energy than the self-trapped solution, contrary to experiment. For the VK center, no metastable local minimum trapped solution was found; however, for the (neutral) STE we find locally metastable solutions, with the on-center STE 0.14 eV higher in energy than the off-center STE solution, which in turn is higher by 0.20 eV than a free electron-hole pair. Even with this discrepancy in the total energy, the atomic positions for the STE solution (see Table I) are reasonable and agree well with Hartree-Fock second-order Møller-Plesset (MP2) perturbation theory . We focus on the properties of the local minimum solution for the off-center STE, which provides a test of the ability of density functional methods to treat the strong coupling of electronic and lattice degrees of freedom, and calculate the spectral properties of excited states.

To solve the LSDA equations we use a plane wave pseudopotential method  with a spin-dependent exchange-correlation potential , and full structural relaxation in a supercell approach. For most calculations, we used a 32 atom supercell with translation vectors (2,0,0), (0,2,0) and (1,1,1), giving a nearest neighbor distance between STEs of 9.4 Å, and used four special k-points in the irreducible wedge for the Brillouin zone integrations. Tests varying the number of k-points and supercell size suggest that these parameters are adequate . One of the in $`xy`$ plane nearest Cl and two Na atoms are most displaced while the rest of the atoms move by a much smaller amount. Some calculated structural parameters for the STE are given in Table I.

The calculated atomic displacements from the ideal NaCl structure for the off-center STE are shown in Fig. 1. The b3u hole state in the triplet STE is localized on the Cl$`{}_{}{}^{}{}_{2}{}^{}`$ molecule (Fig. 1a), with nearly equal weight on the two Cl ions. The last spin up electron (a1g) is mostly localized on the (1/2,1/2,0) vacant halogen site (Fig. 1b), as in the case of the F-center. The formation of the Cl$`{}_{2}{}^{}{}_{}{}^{}`$ “molecule” in the STE is mainly due to the shift of a single Cl. This asymmetric shift can be rationalized by noting that the Madelung energy (with canonical charges of $`\pm `$1 for Na and Cl) of the D$`_{\text{2h}}`$ configuration is about 0.17 eV higher than the C$`_{\text{2v}}`$ one, i.e., the ionic Madelung terms that favor the rock salt structure in the first place favor keeping one of the Cl ions on a lattice site; in addition, the extra electron in the STE (compared to the $`V_K`$ center) can lower its energy by this distortion.

To obtain vibrational properties of the off-center STE, we made finite displacements from the equilibrium geometry. The Cl$`{}_{}{}^{}{}_{2}{}^{}`$ stretching mode $`\omega _{\text{str}}=242\text{cm}^1`$ is smaller than the experimental Raman frequency 361 cm-1 . ¿From the force matrix associated with the Cl$`{}_{}{}^{}{}_{2}{}^{}`$ stretching mode, we found that only the two neighboring Na ions at ($`\frac{1}{2}`$,1,0) and (1,$`\frac{1}{2}`$,0) couple significantly. Unlike H center calculations , coupling to the Na atoms yields only a small shift in the frequency of the Cl$`{}_{}{}^{}{}_{2}{}^{}`$ mode to $`\omega _{\text{str}}=234\text{cm}^1`$. Since the LSDA places the STE too high in energy, it is not too surprising that the curvature of the STE local minimum is underestimated,

Figure 2 shows the computed density of states (DOS) of the perfect NaCl and the inverse participation ratio (IPR=$`\text{a}_0^3/8|\mathrm{\Psi }_i(\stackrel{}{r})|^4d^3r`$) for states of majority and minority spin for the off-center STE. (The IPR is a measure of the localization of a state.) The localized hole b3u and electron a1g states lie in a gap of about 6 eV between the conduction and valence bands of the perfect crystal, where as usual the LSDA underestimates the gap. Rather than using the single-particle eigenvalues, we obtain estimates of the excitation energies as the difference between total energies of different electronic configurations. To calculate the energy to create a free electron-hole pair (the gap energy), we occupy spin up states with one extra electron, while spin down states have one empty state. The energy difference between the two solutions (6.44 eV) for the same atomic positions should correspond to the free electron-hole pair (experimental value 7.96 eV ).

In the distorted STE solution, the energy of electron-hole pair recombination was found by comparing energies of the STE solution with the energy of NaCl having the same atomic displacements as that of the STE. This energy of 4.25 eV is roughly the same as the energy difference between a1g and b3u states given by the DOS, and compares to a value of 3.35 eV obtained from a luminescence experiment  at 11 K. The distorted NaCl, with a lattice distortion energy of 2.5 eV relative to the ideal NaCl positions, is in a highly excited vibrational state. Thus, significant Franck-Condon effects in the spectral properties of the STE are expected. To qualitatively describe the luminescence, a ground state potential curve was calculated for the configuration coordinate $`\alpha `$. It was assumed that all atoms move back to the perfect crystal positions proportionally to their distortions in the STE solution. The result is shown in Fig. 3, along with a quadratic fit $`k_a\alpha ^2/2`$ with one adjustable parameter $`k_a`$. The effective one-dimensional Schrödinger equation for $`\mathrm{\Psi }(\alpha )`$ describes a harmonic oscillator with frequency $`\omega _a=\sqrt{k_a/I}=109\text{cm}^1`$. The moment of inertia $`I=M_n\delta \stackrel{}{R}_n^2`$ was chosen so that $`I\dot{\alpha }^2/2`$ equals the kinetic energy of the atoms with mass $`M_n`$, when they move from the initial displacements $`\delta \stackrel{}{R_n}`$.

The same type of the potential energy curve was calculated for the STE. The quadratic fit works only in the close vicinity of the exciton metastable minimum. The resulting frequency is $`\omega _b=123\text{cm}^1`$. The experimental  temperature T=11 K justifies a zero temperature approximation (the STE initial state is the vibrational ground state). Since the luminescence peak position corresponds to the vibrational level $`n`$$``$170 of the electronic ground state, quasiclassical wavefunctions were used in the numerical integral evaluation. The sequence of vibrational sidebands should be replaced by a sequence of convolved densities of phonon states $`D(\omega )`$. We approximate this by a Gaussian, $`D(\omega )\mathrm{exp}(\omega ^2/2\gamma ^2)/\sqrt{2\pi }\gamma `$ with the width $`\gamma `$=43 cm-1 chosen such that the first three moments coincide with the experimental  phonon DOS. This gives a luminescence width of 0.43 eV, while the experimental width is 0.63 eV.

The optical response $`\sigma (\omega )`$ of the long-lived triplet STE also has been measured . The diagonal part of the optical conductivity tensor is:

$`\sigma _\alpha (\omega )={\displaystyle \frac{\pi e^2N}{m^2\omega \mathrm{\Omega }}}{\displaystyle \underset{\mathrm{}\mathrm{}^{}}{}}{\displaystyle \underset{\stackrel{}{k}}{}}(f_\mathrm{}\stackrel{}{k}f_\mathrm{}^{}\stackrel{}{k})`$ (1)
$`\mathrm{}\stackrel{}{k}|p_\alpha |\mathrm{}^{}\stackrel{}{k}^2\delta (\mathrm{}\omega E_\mathrm{}^{}\stackrel{}{k}+E_\mathrm{}\stackrel{}{k})`$ (2)
where $`f_\mathrm{}\stackrel{}{k}`$ is the occupancy of the state $`|\mathrm{}\stackrel{}{k}`$. The spin state index is included in the band index $`\mathrm{}`$. Integration over the zone has been performed using 26 k-points in the irreducible zone. In Fig. 4a, absolute optical conductivity curves are shown for the three polarizations, $`\stackrel{}{E}`$(1,-1,0), (0,0,1), and (1,1,0) (parallel to the the Cl$`{}_{}{}^{}{}_{2}{}^{}`$ molecular axis). Fig. 4b shows the average conductivity $`\sigma (\omega )=_i\sigma _i(\omega )/3`$, which is compared with experiment , rescaled so that the total weights under the both curves are the same.

The first peak in Fig. 4b is centered at 0.95 eV; the splitting between the peaks for the three different polarizations is not resolved in the calculations. Most of the weight in these peaks comes from transitions of the last localized spin-up a1g electron into empty conduction band delocalized states (see Fig. 2b).

The second peak centered at 3.58 eV is the a$`{}_{1g}{}^{}\text{b}_{3u}`$ transition for the spin down electron (see Fig. 2c). The energy difference between the ground state and electronically excited exciton state with the same atomic configuration turns out to be the same as the eigenvalue difference of b3u and a1g states of the ground STE. The excited exciton will lower its energy by moving atoms back to the undistorted positions of perfect NaCl. To apply the Franck-Condon principle, we repeat the same type of calculations for the excited exciton as we did for NaCl and the ground state STE (Fig. 3). When two Cl atoms move away from each other, the a1g empty state merges with the valence Cl 3$`p`$ band, which makes it very difficult to choose which state to depopulate during the iterations. Instead we used the results for the ground state STE to obtain the energy of the electronically excited state by adding eigenvalue difference $`\lambda _\mathrm{}1\lambda _\mathrm{}2`$ between two states for which dipole matrix element $`\mathrm{}1|p_3|\mathrm{}2^2`$ is the largest. Results are shown on Fig. 3 along with a quadratic fit ($`\omega _c`$=116 cm-1) which works for the entire range of the parameter $`\alpha `$. The delta-function of Eq. (2) corresponding to the a1g$``$b3u transition was replaced by a sequence of convolved Gaussian peaks.

The transition energy (3.58 eV) between the two localized states a1g$``$b3u (see Fig. 2c) agrees well with the experimental peak at 3.8 eV. But when initial and final states have different degrees of localization, errors of an eV in the luminescence (LSDA: 4.25 eV, expt. : 3.35 eV) and in the optical excitation of the bound electron into the conduction band (LSDA: 0.95 eV, expt. : 1.95, 2.13 and 2.00 eV for three different field polarizations) occur. A possible explanation for this error is that in LSDA the incomplete cancellation of the large repulsive self-interaction energy $`𝑑\stackrel{}{r}𝑑\stackrel{}{r^{}}n_i(\stackrel{}{r})n_i(\stackrel{}{r^{}})/|\stackrel{}{r}\stackrel{}{r^{}}|`$ of a localized state $`n_i=|\mathrm{\Psi }_i|^2`$ and the corresponding exchange term may destabilize a localized solution in favor of a delocalized one. Corrections such as self-interaction corrections (SIC)  or LDA+U  may reduce the error; in SIC, shifts in the energy of a state $`i`$ on the order of $`d^3rn_i^{4/3}(\stackrel{}{r})`$  are expected. If the degrees of localization of the two states are similar, however, then the LSDA transition energies are reasonable.

In summary, we have presented LSDA calculations for the STE in NaCl, including the coupling between the lattice and electronic states. The off-center STE is found to be more stable than the on-center STE, but both are metastable compared to free electron-hole pairs. Both luminescence and optical conductivity, including vibrational Franck-Condon effects, were also calculated. The hole state in the off-center STE is found to be rather evenly split between the two Cl atoms, but the electron state is localized to the vacant site left by the shifted Cl ion. The density functional description of electronic transitions between localized states, such as the a$`{}_{1g}{}^{}\text{b}_{3u}`$ absorption peak, agree well with experiment. For transitions between the localized and delocalized states, discrepancies of order 1 eV with experiment arise. Although the LSDA can capture many features of the STE states, when a localized solution competes with a delocalized solution, the incomplete cancellation of the self-interaction may destabilize the localized solution. Given the usefulness of the LSDA method for unraveling complex materials, it is important to test and develop approaches that can treat localized and delocalized states on the same footing.

###### Acknowledgements.

We thank M. L. Cohen, G. W. Fernando, P. M. Johnson, S.G. Louie, W. E. Pickett for discussions. This work was supported in part by NSF Grant No. DMR-9725037 and by DOE Grant No. DE-AC-02-98CH10886.


format: text

text:
# Auditory sensitivity provided by self-tuned critical oscillations of hair cells

## I Generic aspects

Amplification and frequency filtering of a Hopf bifurcation. We discuss the behavior of a dynamical system which is controlled by a parameter $`C`$. Above a critical value, $`C>C_c`$, the system is stable; for $`C<C_c`$ it oscillates spontaneously. At the critical point (or Hopf bifurcation) $`C=C_c`$ the system shows remarkable response and amplification properties which do not depend on the physical mechanism at the basis of the bifurcation. These generic properties can be described as follows. Since we are interested in the response to a periodic stimulus with frequency $`\nu =\omega /2\pi `$, we express the hair-bundle deflection $`x(t)`$ by a Fourier series $`x(t)=x_ne^{in\omega t}`$ with complex amplitudes $`x_n=x_n^{}`$. In the vicinity of the bifurcation, the mode $`n=\pm 1`$ is dominant and the response to an externally-applied sinusoidal stimulus force $`f(t)=f_1e^{i\omega t}+f_1e^{i\omega t}`$ can be expressed in terms of a systematic expansion in $`x_1`$. Symmetry arguments (see Appendix A) imply that the first nonlinear term is cubic:

$$f_1=𝒜x_1+|x_1|^2x_1+\mathrm{}$$
(1)
where $`𝒜(\omega ,C)`$ and $`(\omega ,C)`$ are two complex functions. The bifurcation point is characterized by the fact that $`𝒜`$ vanishes for the critical frequency, $`𝒜(\omega _c,C_c)=0`$. For $`C<C_c`$ and no external force, the system oscillates with $`|x_1|^2\mathrm{\Delta }^2(C_cC)/C_c`$, where $`\mathrm{\Delta }`$ is a characteristic amplitude. For $`C=C_c`$ the response to a stimulus at the critical frequency has amplitude

$$|x_1|||^{1/3}|f_1|^{1/3}.$$
(2)
This represents an amplified response at the critical frequency with a gain

$$r=\frac{|x_1|}{|f_1|}|f_1|^{2/3}$$
(3)
that becomes arbitrarily large for small forces.

If the stimulus frequency differs from the critical frequency, the linear term in Eq. (1) is non-zero and can be expressed to first order as $`𝒜(\omega ,C_c)A_1(\omega \omega _c)`$. The dramatic amplification of weak signals, implied by Eq. (3), is maintained as long as this term does not exceed the cubic term in Eq. (1). If the frequency mismatch increases such that $`|\omega \omega _c||f_1|^{2/3}||^{1/3}/|A_1|`$, the response becomes linear

$$|x_1|\frac{|f_1|}{|(\omega \omega _c)A_1|}$$
(4)
i.e. the gain is independent of the strength of the stimulus.

Thus the Hopf resonance acts as a sharply tuned high-gain amplifier for weak stimuli, and as a low-gain filter for strong stimuli. This generic behavior is illustrated in Fig. 1 with data obtained by numerical simulation. Laser interferometry measurements of the motion of the basilar membrane, when the live cochlea is stimulated by pure tones, display strikingly similar features . In particular, the peak response as a function of force amplitude has been demonstrated to obey a power law $`|x_1||f_1|^{0.4\pm 0.2}`$ . This strongly suggests that the membrane is being driven by a dynamical system which is poised at a Hopf bifurcation.

Self-tuned critical oscillations. How does the system come to be so precisely balanced at the critical point?***Self-tuning to a coexistence point has been discussed previously in certain dynamic first-order transitions . The control parameter must be tuned to $`CC_c`$, otherwise the nonlinear amplification is lost. Moreover, the value of $`C_c`$ differs for hair cells with different characteristic frequency. We propose a feedback mechanism which allows the dynamical system to operate automatically close to the bifurcation point, whatever its characteristic frequency. Without loss of generality, we assume that the control parameter decreases as long as the system does not oscillate. After some time, critical conditions are reached and spontaneous oscillations ensue. The onset of oscillations triggers an increase of the control parameter which tends to restore stability. Hence the system converges to an operating point close to the bifurcation point. To illustrate this general idea, we consider the following simple feedback which changes $`C`$ in response to deflections $`x`$:

$$\frac{1}{C}\frac{C}{t}=\frac{1}{\tau }\left(\frac{x^2}{\delta ^2}1\right)$$
(5)
where $`\delta `$ is a typical amplitude. If no external force is applied, this feedback, after a relaxation time $`\tau `$, tunes the control parameter to a value $`C_\delta `$ for which spontaneous oscillations with $`|x_1|\delta `$ occur. If $`\delta `$ is small compared to the characteristic amplitude $`\mathrm{\Delta }`$, this is on the oscillating side close to the bifurcation, $`|C_\delta C_c|/C_c(\delta /\mathrm{\Delta })^2`$. Two modes of signal detection are possible: (i) For transient stimuli short compared to $`\tau `$ the system operates at $`C_\delta `$. The amplitude $`|x_1|`$ shows the characteristic nonlinear response discussed above. (ii) For stimuli sustained over longer times, self-tuning maintains $`|x_1|\delta `$ constant for different stimulus amplitudes. This effect of the feedback represents a perfect adaptation mechanism. However, in the presence of noise, phase-locking of the response (to be discussed later) occurs as soon as an external stimulus is applied, and this can be detected.

## II Model

Mechanosensitivity and self-tuning mechanism provided by transduction channels. We demonstrate the general principles introduced above by devising a specific model for the amplification of acoustic stimuli by hair bundles in non-mammalian vertebrates. A schematic representation of a hair bundle is shown in Fig. 2. It consists of several stereocilia and a single kinocilium . Transduction of hair bundle deflection to a chemical signal occurs via channels located near the tip of each stereocilium. Tip links which connect neighboring stereocilia are believed to be the gating springs of the transduction channels. If the hair bundle is deflected, tension in the tip links triggers the opening of the channels. The subsequent influx of ions (principally $`\mathrm{K}^+`$, but also $`\mathrm{Ca}^{2+}`$) causes a corresponding change of the membrane potential of the hair cell which, in turn, generates a nervous signal.

The mechanosensor can be used for self-tuning, as well as signal detection. Many physiological processes are regulated by ionic concentrations, so it is natural to identify the $`\mathrm{Ca}^{2+}`$ concentration with the control parameter $`C`$. We assume that $`C`$ decreases if the transduction channels are closed, owing to the action of pumps in the cell membrane. When the hair bundle is deflected by $`x`$, the transduction channels open with probability $`P_o(x)`$. We therefore characterize the mechanosensor by the equation

$$\frac{C}{t}=\frac{C}{\tau }+J_oP_o(x)$$
(6)
where $`J_o`$ is the $`\mathrm{Ca}^{2+}`$ flux through open transduction channels. Note that this equation provides self-tuning. In this case it replaces the more simple but less realistic Eq. (5). For our numerical examples, we use a two-state model for the channels with

$$P_o(x)=\frac{1}{1+Ae^{x/\delta }}$$
(7)
where $`(1+A)^11`$ is the probability that a channel is open when the hair cell is quiescent and $`\delta `$ is the characteristic amplitude of motion to which the system is sensitive. For a sufficiently long relaxation time $`\tau `$, the slow variation $`C_0`$ of $`CC_0(t)+C_1e^{i\omega t}`$ can be separated from the small-amplitude oscillations, giving

$$_tC_0\frac{C_0}{\tau }+J_o\stackrel{~}{P}_o(|x_1|^2)$$
(8)
where $`\stackrel{~}{P}_o=_0^{1/\nu }𝑑tP_o(x_1e^{i\omega t}+x_1e^{i\omega t})`$ is the averaged probability of channel opening in the presence of oscillations, which increases monotonically with their amplitude $`x_1`$. For physically relevant parameter values, the system reaches a steady state close to the bifurcation point, independent of the initial value of $`C`$.

Oscillations generated by molecular motors. We still have to specify the nature and location of the oscillator within the hair bundle. It has been suggested that the transduction channels might be the source of the instability ; or, that myosin motors in the stereocilia might generate the force necessary to move the bundle . We propose a third possibility: the kinocilium could vibrate using its internal dynein motors. Recently, a simple physical mechanism has been proposed which allows motor proteins operating in collections to generate spontaneous oscillations by traversing a Hopf bifurcation . Typically, motors move along cytoskeletal filaments and elastic elements oppose this motion. In this case, two possibilities exist: The system either reaches a stable balance between opposing forces, or it oscillates around the balanced state. Three time scales characterize this behavior: The relaxation time $`\lambda /K`$ of passive relaxation, where $`\lambda `$ is the total friction and $`K`$ is the elastic modulus; and the times $`\omega _1^1`$ and $`\omega _2^1`$, where $`\omega _1`$ is the kinetic rate at which a motor detaches from a filament, and $`\omega _2`$ is the attachment rate. An explicit solution of a simple model is derived in Appendix B. We find that for an appropriate value of a control parameter $`C`$, which is related to the ratio $`\omega _1/\omega _2`$, a Hopf bifurcation occurs with critical frequency

$$\omega _c\left(\frac{K\alpha }{\lambda }\right)^{1/2},$$
(9)
which is the geometric mean of the passive relaxation rate $`K/\lambda `$ and the typical ATP hydrolysis rate $`\alpha =\omega _1+\omega _2`$. The above identification of $`C`$ with the $`\mathrm{Ca}^{2+}`$ concentration is consistent with the fact that $`\mathrm{Ca}^{2+}`$ regulates motor protein activity . The system oscillates if the elastic modulus does not exceed a maximal value $`K_{\mathrm{max}}k_0N`$, where $`k_0`$ is the crossbridge elasticity of a motor and $`N`$ is the total number of motors. The maximal frequency, obtained when $`K=K_{\mathrm{max}}`$, can be significantly higher than the ATP hydrolysis rate $`\alpha `$.

Characteristic frequency of a vibrating kinocilium. The kinocilium is a true cilium containing a cylindrical arrangement of microtubule doublets and dynein motors. Since cilia have the well-established tendency to beat and vibrate with frequencies from tens of Hz up to at least 1 kHz , the kinocilium is a natural candidate to be responsible for the Hopf bifurcation. A simple two-dimensional model can be used to discuss the main physical properties of a vibrating cilium near a Hopf bifurcation . In this model, motors induce the bending of a pair of elastic filaments separated by a distance $`a`$ (corresponding to the distance between neighboring microtubule doublets in the axoneme). An isolated kinocilium of length $`L`$ and bending rigidity $`\kappa `$, fixed at the basal end but free at its tip, will vibrate in a wave-like mode with wavelength $`\mathrm{\Lambda }4L`$. A typical displacement $`z`$ of the motors along the filaments leads to bending of the filament pair and a deflection of the tip by a distance $`xzL/a`$. The elastic bending energy is of order $`(\kappa /L^3)x^2(\kappa /a^2L)z^2`$. Therefore, the total elastic modulus experienced by the motors is given by

$$K\frac{\kappa }{La^2}.$$
(10)
The viscous energy dissipation per unit time due to motion of the kinocilium is of order $`\eta L(_tx)^2(\eta L^3/a^2)(_tz)^2`$. Therefore, the total friction experienced by the motors can be written as

$$\lambda \frac{\eta L^3}{a^2}+\lambda _0\rho L,$$
(11)
where $`\lambda _0`$ describes the dissipation within the kinocilium per motor and $`\rho =N/L`$ denotes the number of dynein motors per unit length along the axoneme. If internal friction can be neglected, i.e. if $`LL_0=(\lambda _0a^2\rho /\eta )^{1/2}`$, the frequency of a vibrating cilium at the bifurcation point is given by

$$\omega _c\left(\frac{\kappa \alpha }{\eta }\right)^{1/2}\frac{1}{L^2}.$$
(12)
Using typical values $`\lambda _010^9\mathrm{kg}/\mathrm{s}`$, $`a20\mathrm{nm}`$, $`\rho \mathrm{5\hspace{0.33em}10}^8\mathrm{m}^1`$ and $`\eta 10^3\mathrm{kg}/\mathrm{ms}`$, we find $`L_0200\mathrm{nm}`$, shorter than typical kinocilia. Using $`\alpha 10^3\mathrm{s}^1`$, $`k_010^3\mathrm{N}/\mathrm{m}`$ and $`\kappa \mathrm{4\hspace{0.33em}10}^{22}\mathrm{Nm}^2`$ (the bending rigidity of 20 microtubules), the frequency range between 100$`\mathrm{Hz}`$ and 10$`\mathrm{kHz}`$ can naturally be spanned by changing the length of the kinocilium between $`1\mu \mathrm{m}`$ and $`10\mu \mathrm{m}`$.

The above argument neglects the contribution of the stereocilia to elasticity and motion. The elastic response of stereocilia to hair bundle deflections has been measured . It can be well described by an angular spring at the base of each stereocilium which contributes an elastic energy per stereocilium of the order of $`k_s(x/L)^2`$, where $`k_s`$ is an angular elastic modulus. The kinocilium length $`L`$ in a hair bundle is approximately inversely proportional to the number $`N_s`$ of stereocilia , and we write $`N_sl_s/L`$, where $`l_s10^4\mathrm{m}`$ is the total length of stereocilia. With this assumption, we find an additional elastic modulus

$$K_s\frac{l_sk_s}{La^2}$$
(13)
contributed by the stereocilia to $`K`$. The contribution of stereocilia to the friction coefficient $`\lambda `$ can be estimated as

$$\lambda _s\eta l_sL^2/a^2.$$
(14)
The measured value of $`k_s`$  indicates that $`K_s`$ dominates the contribution to $`K`$ given by Eq. (10). Similarly, since $`l_s>L`$, $`\lambda _s`$ should dominate friction. In this case, we expect

$$\omega _c\left(\frac{k_s\alpha }{\eta L^3}\right)^{1/2}$$
(15)
and the range of frequencies is somewhat reduced.

If the stiffness of the ensemble of stereocilia greatly exceeds that of the kinocilium, a new situation arises. Since the kinocilium is attached to the stereocilia at its tip , movement of the tip is strongly reduced and the kinocilium preferentially vibrates in a mode with wavelength $`\mathrm{\Lambda }2L`$ for which the relation given by Eq. (12) is again valid.

Adaptation motors. Finally, we note that in order to obtain a robust self-tuning, the feedback mechanism must be sensitive only to the oscillation amplitude $`x_1`$ and not to a stationary displacement $`x_0`$; we assumed in Eq. (8) that $`\stackrel{~}{P}_o`$ is independent of $`x_0`$. The transduction mechanism of stereocilia and their tension-gated channels does indeed have this property: It is well known that an ATP-dependent adaptation mechanism  exists which removes the dependence of the channel current on a constant displacement $`x_0`$. It is widely believed that this adaptation involves the motion of myosin motors, which maintains constant the steady-state tension in the tip links that control the transduction channels . Therefore, the stereociliar transduction mechanism has precisely the required properties to be used as a feedback signal for self-tuning the bifurcation.

## III Simulation

The above analysis does not include the effects of noise. Brownian motion is one source of fluctuations in the movement of the hair bundle . Another source of noise is caused by stochastic fluctuations in motor protein force as dynein molecules bind to, and detach from the microtubules in the kinocilium. In order to investigate the consequences of this randomness, we performed a Monte Carlo simulation of the two-state model described in , using a realistic number of motor proteins. The motors, when attached, experienced a potential of amplitude $`U`$ and period $`l`$ (the corresponding crossbridge elasticity is $`k_0=U/l^2`$). The attachment rate was constant and independent of position, $`\omega _2=\alpha `$. Detachment was localized to a region of width $`0.1l`$, centred on the potential minimum and the detachment rate $`\omega _1`$ was regulated by the $`\mathrm{Ca}^{2+}`$ concentration $`C`$. We chose $`\omega _1/\omega _2=(C_q/C)^3`$, where $`C_q`$ is the steady-state $`\mathrm{Ca}^{2+}`$ concentration in a quiescent hair cell; the precise functional dependence is unimportant as long as $`\omega _1/\omega _2`$ decreases monotonically with increasing $`C`$ in a fairly sensitive way. We simulated systems with $`N`$ = 1000-4000 motors, representing hair cells with different kinocilium length and stereociliary number.

Self-tuning and the characteristic frequency of spontaneous oscillations. The self-tuning of a hair cell to the vicinity of the bifurcation, where small-amplitude spontaneous oscillations occur, is demonstrated in Fig. 3. When a change of the internal $`\mathrm{Ca}^{2+}`$ concentration is imposed, the system is transiently perturbed; but after an interval of time of order $`\tau `$, it returns to the same steady state. Hair cells with different numbers of motors acquire different internal concentrations of $`\mathrm{Ca}^{2+}`$, in order to adjust the motor detachment rate in such a way that the system approaches the critical point. The spontaneous oscillations of three different hair bundles are shown in Fig. 4. Note that the characteristic frequency is approximately proportional to the inverse-square of the number of motors, as expected from Eq. (12), and that it can exceed the typical ATP cycle rate $`\alpha `$ when the total number of motors is small (short kinocilium). All three hair cells execute spontaneous oscillations with a similar amplitude, as expected from our arguments above. However, the noise introduces a significant new effect: The oscillations are irregular. The incoherence of the phase of the oscillation is evident in the Fourier transform of the displacement, which exhibits a broad peak centred on the characteristic frequency.

Dynamic response to a tone at the characteristic frequency. The response of a self-tuned hair bundle to a sinusoidal force with a frequency approximately equal to the bundle’s characteristic frequency is illustrated in Fig. 5. For weak stimuli,the amplitude of the oscillation does not increase with the amplitude $`|f_1|`$ of the applied force; this is because the small response to the stimulus is masked by the noisy, spontaneous motion. Instead, the phase of the hair-bundle oscillation becomes more regular; as it does so, a peak emerges from the Fourier spectrum at the driving frequency. The height of the peak grows approximately as $`|f_1|^{1/3}`$ for intermediate values of $`|f_1|`$, and approximately linearly for very weak stimuli, and also for very strong stimuli (for which the response is essentially passive). Thus the Fourier component of the hair-bundle displacement at the driving frequency responds to the stimulus in the generic manner discussed above.

Phase-locking of the nervous signal. The flow of ions through the transduction channels depolarizes the cell membrane which, in turn, opens voltage-gated channels at the base of the hair cell and generates a synaptic current . This sequence of events happens fast enough for variations at the synapse faithfully to reflect the hair-bundle motion at frequencies below 1 kHz. Information about the auditory stimulus is subsequently passed along the auditory nerve in the form of a spike train. Simplifying this transduction process, we assume that a spike is elicited whenever the transduction channel current $`J`$ passes a threshold value of $`0.5J_o`$. The resulting nervous response is shown in Fig. 6. In the absence of a stimulus, the self-tuned critical oscillations of the bundle cause the nerve to fire stochastically, at a low rate. When a weak sinusoidal stimulus is applied at the characteristic frequency, the firing rate does not increase above the spontaneous rate, but the spike train becomes detectably phase-locked to the stimulus. The degree of phase-locking increases rapidly as the amplitude of the stimulus increases, reflecting the growing regularity of the hair-bundle motion. It is only when the neural response is almost completely phase-locked that the firing rate begins to rise. Eventually, for strong stimuli, the spike rate saturates at the stimulus frequency. This behavior is strikingly similar to that which is observed experimentally. In particular, it is well known that the threshold for phase-locking in the auditory nerve fiber is 10-20 dB lower than the threshold at which the firing rate begins to rise .

## IV Discussion

The self-tuned critical oscillator which we have introduced as a system for signal detection has characteristic amplification properties which are generic and do not depend on the choice of a specific model. The oscillation frequency, however, depends on the physical mechanism involved. In our model of non-mammalian hair cells, the frequency is determined by the geometry of the hair bundle. A simple morphological gradient along the basilar membrane would endow the ear with the ability to analyze a wide range of frequencies. Experimentally, it is well established that the height of hair bundles progressively increases along the cochlea, and that concurrently the characteristic frequency of the hair cells declines . Our proposition that the kinocilium is likely to play an active role in non-mammalian vertebrate hair cells suggests experiments which study the motility of the kinocilium and its potential for generating oscillations. If the kinocilium is the only source of oscillations, its removal should suppress hair-bundle vibrations. If the self-tuning mechanism is removed instead (by cutting the tip links, for example) our model predicts that the kinocilium should exhibit stronger spontaneous oscillations. Careful control of the extracellular ion concentrations, such as $`\mathrm{Ca}^{2+}`$ would be essential in such experiments.

The simulations of our model reveal how a hair bundle can achieve its remarkable sensitivity to weak stimuli. By profiting from the periodicity of a sinusoidal input, and measuring phase-locking rather than the amplitude of response, the mechanosensor can detect forces considerably weaker than those exerted by a single molecular motor (if the bundle were a simple, passive structure, its response to such forces would be smaller than its Brownian motion). An important implication of this detection mechanism is that even though the hair cell selects a certain frequency, the signal must still be encoded by the interval between spikes elicited in the auditory nerve. Paradoxically, the stochastic noise caused by the motor proteins serves a useful purpose. It ensures that the self-tuned critical oscillations of the hair bundle are incoherent, so that the pattern of spontaneous firing in the nerve is irregular. Against this background, the regular response to a periodic stimulus can easily be detected.The benefits of noise have been discussed in a variety of situations such as those involving stochastic resonance .

Another beneficial feature of noise arises from the fact that weak stimuli do not increase the amplitude of oscillation above the spontaneous amplitude. Thus the $`\mathrm{Ca}^{2+}`$ concentration remains constant, the hair bundle stays in the critical regime, and active amplification can be sustained indefinitely. Stronger stimuli cause the system to drift away from the critical point, so that the degree of amplification diminishes over time. It is well known that both the perception of loudness  and the firing rate of the auditory nerve  decrease over a period of a few seconds, when a stimulus of moderate intensity is maintained. This phenomenon, which is usually referred to as ‘adaptation’, is consistent with our self-tuning mechanism. Following the sustained presentation of a loud stimulus, the spontaneous firing rate diminishes and the threshold to weak stimuli is augmented . Such ‘fatigue’ is also accounted for by our self-tuning mechanism.

Since self-tuning positions the system slightly on the oscillating side of the critical point, self-tuned criticality provides a natural explanation for otoacoustic emissions . In its normal working state, the inner ear would generate faint sounds with a broad range of frequencies. If the feedback mechanism were to fail in certain cells, the spontaneous oscillations could become large enough for distinct tones to be emitted.

A self-tuned Hopf bifurcation is ideal for sound detection because it provides sharp frequency selectivity and a nonlinear gain which compresses a wide range of stimulus intensities into a narrow range of response. We therefore believe that the concept applies to all vertebrate hearing systems, and potentially to other mechanoreceptor systems. Kinocilia are typically absent in mammalian cochlea, and we suggest that their force-generating role has been assumed by the outer hair cells. Self-tuning of these motile cells might be realized by a mechanism similar to that presented here, using transduction channels in their hair bundles. It could, however, work very differently; for example, it might involve feedback from the inner hair cells, via the efferent nervous system.

Tuning to the proximity of a critical point is likely to be a general strategy adopted by sensory systems. Simple molecular receptors , as well as the physiological sensors of higher organisms, can enhance their response to weak stimuli in this way. We propose that the physics of self-tuned criticality is the ‘central science of transducer physiology’ spoken of by Delbrück .

Acknowledgements: We thank P. Martin and C. Petit for useful discussions, and A.F. Huxley for referring us to the work of T. Gold. T. Duke is grateful for the hospitality of Institut Curie and the Niels Bohr Institute, and acknowledges the support of the Royal Society. After submission of our manuscript we learned that Eguiluz, Ospeck, Choe, Hudspeth and Magnasco have independently described the generic response of a system near a Hopf bifurcation; we thank them for communicating their unpublished results.

## A Generic behavior at a Hopf-bifurcation

### 1 Nonlinear relation between periodic stimulus and displacements

We are interested in the response $`x(t)`$ of a nonlinear system to a periodic stimulus force $`f(t)`$. If only one frequency $`\nu =\omega /2\pi `$ is present we use the Fourier expansions

$`f(t)`$ $`=`$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}f_ne^{in\omega t}`$ (A1)
$`x(t)`$ $`=`$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}x_ne^{in\omega t},`$ (A2)
where the complex coefficients $`x_n`$ and $`f_n`$ obey $`x_n=x_n^{}`$ and $`f_n=f_n^{}`$. This representation implies that we focus on the limit cycle solution and ignore all transient relaxation phenomena. We consider the class of systems for which the force at a given time depends in a nonlinear way on the history of the displacements $`x(t)`$ alone; as we will discuss in section D more complex cases do not change the basic properties. In this situation, the relation between $`x`$ and $`f`$ can be expressed as a systematic expansion of the force amplitudes $`f_n`$ in the amplitudes $`x_n`$:

$`f_k=F_{kl}^{(1)}x_l`$ $`+`$ $`F_{klm}^{(2)}x_lx_m`$ (A3)
$`+`$ $`F_{klmn}^{(3)}x_lx_mx_n+O(x^4),`$ (A4)
where the expansion coefficients $`F_{k,k_1,..,k_n}^{(n)}`$ are symmetric with respect to permutations of the indices $`k_1..k_n`$. The limit cycle solutions are invariant with respect to translations in time $`tt+\mathrm{\Delta }t`$. Under these transformations the amplitudes change as $`x_nx_ne^{in\omega \mathrm{\Delta }t}`$ and $`f_nf_ne^{in\omega \mathrm{\Delta }t}`$. Inspection of Eq. (A4) shows that the time translation symmetry allows only for those terms $`F_{k,k_1,..,k_n}^{(n)}x_{k_1}..x_{k_n}`$ for which $`k=k_1+\mathrm{}+k_n`$. For all other cases $`F_{k,k_1,..,k_n}^{(n)}`$ must vanish which significantly restricts the number of terms.

### 2 Hopf bifurcation

The nonlinear system exhibits spontaneous oscillations and a Hopf-bifurcation if nontrivial solutions to Eq. (A4) with $`x_n0`$ exist in the case where all $`f_k=0`$, i.e. if no stimulus force is applied. Without loss of generality, we consider here an instability of the mode $`x_1`$. In this case, the dominant terms allowed by symmetry read ($`f_k=0`$)

$`0`$ $``$ $`F_{11}^{(1)}x_1+2F_{1,2,1}^{(2)}x_1x_2`$ (A5)
$`+`$ $`6F_{1,1,1,1}^{(3)}x_1^2x_1+6F_{1,1,2,2}^{(3)}x_2x_2x_1`$ (A6)
$`0`$ $``$ $`F_{22}^{(2)}x_2+2F_{211}^{(2)}x_1^2.`$ (A7)
Eq. (A7) determines $`x_22(F_{211}^{(2)}/F_{22}^{(2)})x_1^2`$. Inserting this relation in Eq. (A6), we obtain to lowest order

$$0𝒜x_1+|x_1|^2x_1,$$
(A8)
where $`𝒜F_{11}^{(1)}`$ and $`3F_{1,1,1,1}^{(3)}4F_{211}^{(2)}F_{1,2,1}^{(2)}/F_{22}^{(2)}`$.

The coefficients $`𝒜(\omega ,C)`$ and $`(\omega ,C)`$ are complex and in general depend on frequency $`\omega `$ and a control parameter which we denote by $`C`$. A Hopf bifurcation occurs at a critical point $`C=C_c`$ at which $`𝒜`$ vanishes for a frequency $`\omega _c`$, i.e. $`𝒜(\omega _c,C_c)=0`$. This can be demonstrated as follows: A spontaneously oscillating solution satisfies

$$|x_1|^2=\frac{𝒜}{}$$
(A9)
Note, that such a solution can only exist if $`𝒜/`$ is real and negative. At the bifurcation point, $`𝒜=0`$ and $`𝒜/`$ is therefore real for $`\omega =\omega _c`$, however the corresponding amplitude $`|x_1|^2`$ vanishes. In the vicinity of this point we expect to find solutions with finite amplitude. We use the expansion

$$𝒜(\omega ,C)(\omega \omega _c)A_1+(CC_c)A_2$$
(A10)
where $`A_1`$ and $`A_2`$ are complex coefficients and we neglect higher order terms. Spontaneous oscillating solutions exist only if $`𝒜/`$ is real. This condition is satisfied for a particular frequency $`\omega =\omega _s`$ with

$$\omega _s=\omega _c+\frac{Im(A_2/)}{Im(A_1/)}(C_cC).$$
(A11)
The ratio $`𝒜/`$ at this frequency $`\omega _s`$ changes sign for $`C=C_c`$; here we assume without loss of generality that it is positive for $`C<C_c`$. In this case, the system oscillates spontaneously with an amplitude which according to Eq. (A9) behaves as $`|x_1|^2=\mathrm{\Delta }^2(C_cC)/C_c`$, where

$$\mathrm{\Delta }^2=C_c\left(Re(A_2/)Re(A_1/)\frac{Im(A_2/)}{Im(A_1/)}\right)$$
(A12)
is a typical amplitude. We have thus demonstrated that Eq. (A8) characterizes a Hopf-bifurcation if the complex coefficient $`𝒜`$ vanishes at a critical point $`C_c`$ for a critical frequency $`\omega _c`$.

### 3 Amplified response to sinusoidal stimuli

If a sinusoidal stimulus $`f(t)=f_1e^{i\omega t}+f_1e^{i\omega t}`$, for which all $`f_n`$ with $`n\pm 1`$ vanish, Eq. (A8) becomes

$$f_1𝒜x_1+|x_1|^2x_1.$$
(A13)
We consider a system that is tuned exactly to the bifurcation, $`C=C_c`$. In this situation spontaneous oscillations do not occur and $`𝒜=(\omega \omega _c)A_1`$. If the imposed frequency is equal to the critical frequency $`\omega =\omega _c`$, the coefficient $`𝒜`$ vanishes and we can solve Eq. (A13) for $`|x_1|`$ to find the nonlinear response

$$|x_1|||^{1/3}|f_1|^{1/3},$$
(A14)
as a function of the force amplitude $`|f_1|`$. This behavior represents an amplified response with a gain

$$r=\frac{|x_1|}{|f_1|}|f_1|^{2/3}$$
(A15)
that becomes arbitrarily large for small forces. If the frequency $`\omega `$ is different from $`\omega _c`$, this nonlinear response still holds as long as the linear term in Eq. (A13) is small compared to the cubic term and can be neglected. This is the case if $`|x_1|^2|𝒜/|=|\omega \omega _c||A_1/|`$. Therefore, the nonlinear regime characterized by Eq. (A14) holds for sufficiently large force amplitudes, $`|f_1||(\omega \omega _c)A_1|^{3/2}/||^{1/2}`$, or if the frequency is sufficiently close to the critical frequency, $`|\omega \omega _c||f_1|^{2/3}||^{1/3}/|A_1|`$.

If the frequency mismatch $`|\omega \omega _c|`$ becomes large, or if forces $`|f_1|`$ are small, a new regime occurs for which the linear term in (A13) dominates. In this regime, the response is linear,

$$|x_1|\frac{|f_1|}{|(\omega \omega _c)A_1|},$$
(A16)
and the gain is constant. This is a passive response if the stimulus frequency is too far from the critical frequency.

### 4 Additional remarks

The above derivation is based on an expansion (A4) in the displacements $`x_n`$. This excludes some nonlinearities in the force which can lead to additional nonlinear terms in Eq. (A13). The most general form of Eq. (A13) is

$`f_1𝒜x_1`$ $`+`$ $`|x_1|^2x_1+𝒞x_1|f_1|^2+𝒟x_1f_1^2`$ (A17)
$`+`$ $`|x_1|^2f_1+x_1^2f_1+𝒢|f_1|^2f_1.`$ (A18)
However, for small forces $`f_1`$ and small amplitudes $`x_1`$, the results derived above are not affected. The regime of nonlinear reponse $`|f_1||x_1|^3`$, as well as the linear response regime $`|f_1||x_1|`$ still exist. If $`|f_1||x_1|`$, the nonlinear terms in $`f_1`$ renormalize the third order term, which in this regime is negligable. If $`|f_1||x_1|^3`$, the nonlinear terms in $`f_1`$ are of even higher order and can be neglected.

## B Oscillations generated by molecular motors

### 1 Two state model

The two state model describes force-generation as a result of transitions between two states, a bound state and a detached state of a motor and its track filament. The interaction between a motor at position $`z`$ along the filament in states $`1`$ and $`2`$ is characterized by two periodic potentials $`W_1(z)=W_1(z+l)`$ and $`W_2(z)=W_2(z+l)`$ where $`l`$ is the period. We introduce the relative position $`\xi =z\mathrm{mod}l`$ with respect to the potential period. Detachment and attachment rates are denoted $`\omega _1(\xi )`$ and $`\omega _2(\xi )`$, respectively. Oscillations can occur in this model if a large number $`N`$ of motors move collectively against an external elastic element of modulus $`K`$.

We introduce the probability $`P_1(\xi )`$ and $`P_2(\xi )`$ of finding a motor bound at position $`\xi `$ in state $`1`$ or $`2`$, which satisfy the normalization condition

$$_0^l𝑑\xi (P_1+P_2)=1$$
(B1)
For a large number of motors collectively moving with the same velocity $`v`$ the dynamic equations read

$`_tP_1+v_\xi P_1`$ $`=`$ $`\omega _1P_1+\omega _2P_2`$ (B2)
$`_tP_2+v_\xi P_2`$ $`=`$ $`\omega _1P_1\omega _2P_2`$ (B3)
The velocity $`v`$ is determined by the force-balance condition

$$f=\lambda v+Kz+N_0^l𝑑\xi (P_1_\xi W_1+P_2_\xi W_2)$$
(B4)
where $`\lambda `$ is a friction coefficient describing the total friction and $`z`$ is the displacement of the motors, $`_tz=v`$. For an incommensurate arrangement of motors with respect to the track filament and a large number $`N`$ of motors, $`P_1(\xi )+P_2(\xi )=1/l`$ and the equations of motions simplify:

$$_tP+v_\xi P=(\omega _1+\omega _2)P+\omega _2/l,$$
(B5)
where we denote for simplicity $`P(\xi )=P_1(\xi )`$.

We discuss a simple choice for the potentials and transition rates for which the Hopf bifurcation is easy to determine analytically. We consider the potential

$$W_1(\xi )=U\mathrm{cos}(2\pi \xi /l)$$
(B6)
with amplitude $`U`$, and the potential $`W_2`$ to be constant. The transition rates are chosen to be periodic functions

$`\omega _1(\xi )`$ $`=`$ $`\beta \beta \mathrm{cos}(2\pi \xi /l)`$ (B7)
$`\omega _2(\xi )`$ $`=`$ $`\alpha \beta +\beta \mathrm{cos}(2\pi \xi /l)`$ (B8)
parameterized by two coefficients $`\alpha `$ and $`\beta `$. With this choice,

$$\omega _1(\xi )+\omega _2(\xi )=\alpha $$
(B9)
is constant and the fact that $`\omega _1`$ and $`\omega _2`$ are positive restricts $`\beta `$ to the interval $`0\beta \alpha /2`$.

### 2 Linear response function

In order to determine the linear coefficient $`𝒜`$ which determines the stability of the system, we look for small amplitude oscillations close to the resting state with $`v=0`$. We write

$`P`$ $``$ $`p_0+p_1e^{i\omega t}`$ (B10)
$`f`$ $``$ $`f_1e^{i\omega t}`$ (B11)
$`z`$ $``$ $`z_1e^{i\omega t}`$ (B12)
where $`p_0=\omega _2/\alpha l`$. To linear order in $`z_1`$, we find from Eq. (B5)

$$p_1=\frac{i\omega z_1}{i\omega +\alpha }_xp_0$$
(B13)
The corresponding force is given by

$$f_1𝒜z_1$$
(B14)
with

$$𝒜=i\omega \lambda +K+\chi ,$$
(B15)
where the active response $`\chi `$ of the motors is given by

$$\chi =N_0^l𝑑\xi \frac{i\omega }{i\omega +\alpha }_\xi p_0_\xi W_1$$
(B16)
For the choice of Eq. (B6) and (B8) the integral can be calculated and we obtain

$$𝒜(C,\omega )=i\omega \lambda +KNk_0C\frac{i\omega /\alpha +(\omega /\alpha )^2}{1+(\omega /\alpha )^2}.$$
(B17)
Here, we have introduced the dimensionless control parameter $`C2\pi ^2\beta /\alpha `$ with $`0<C<\pi ^2`$ and the cross-bridge elasticity $`k_0U/l^2`$ of the motors.

### 3 Hopf bifurcation

A Hopf bifurcation occurs if there is a pair of values $`(C,\omega )`$ for which $`𝒜`$ as given by Eq. (B17) vanishes. Such a point indeed exists. For the critical value

$$C_c=\frac{\lambda \alpha +K}{Nk_0}$$
(B18)
the bifurcation occurs for the critical frequency

$$\omega _c=\left(\frac{K\alpha }{\lambda }\right)^{1/2}$$
(B19)
The critical frequency is bounded by the fact that $`C_c<\pi ^2`$. The maximal frequency occurs for the maximal possible value of $`K`$

$$K_{\mathrm{max}}=Nk_0\pi ^2\lambda \alpha $$
(B20)
for which $`C_c=\pi ^2`$. This frequency is given by

$$\omega _{\mathrm{max}}=\alpha \left(N\pi ^2\frac{k_0}{\lambda \alpha }1\right)^{1/2}$$
(B21)
Note, that the maximal frequency can be significantly higher than the typical rate $`\alpha `$ of the chemical cycle.


text:
# The Enigmatic Radio Afterglow of GRB 991216

## 1 Introduction

The intense gamma-ray burst GRB 991216 was detected on 1999 December 16.67 UT by the Burst and Transient Experiment (BATSE) on board the Compton Gamma Ray Observatory satellite (Kippen, Preece & Giblin (1999)). Follow-up observations with the PCA instrument on the Rossi X-ray Timing Explorer (RXTE) satellite resulted in the detection of a previously uncataloged X-ray source, which was subsequently seen to fade by a factor of five, seven hours later (Takeshima et al. (1999)). Uglesich et al. (1999) identified a fading optical source, at a position consistent with the RXTE transient, and shortly thereafter the radio counterpart was discovered (Taylor & Berger (1999)).

Here we present radio measurements of this burst from 1 GHz to 350 GHz. While the emission from X-ray and optical afterglow was fairly typical (Halpern et al. (2000)), the radio afterglow of GRB 991216 was unusual in two respects. First, the onset of the decay began much earlier than that in most radio afterglows. Second, the temporal decay indices in the radio, optical and X-ray bands are markedly different from each other. We explore a number of possible explanations for these behaviors.

## 2 Observations

Very Large Array (VLA111The NRAO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. NRAO operates both the VLA and the VLBA.): A log of the observations and flux density measurements are summarized in Table 1. We used J0509+1011 (at 8.46 GHz and 4.86 GHz) and J0530+135 (at 1.43 GHz) for phase calibration. J0542+498 was used for flux calibration at all frequencies.

Very Long Baseline Array (VLBA6): A single 2-hr observation was carried out at 8.42 GHz and 2-bit samples of a 64 MHz bandwidth signal in one hand of polarization were recorded. The nearby ($`<1.1^{}`$) calibrator J0509+1011, a core-jet source, was observed every 3 minutes for delay, fringe rate and fringe phase calibration. The total flux density of the calibrator was found to be 9.5% less than was measured by the VLA on the same day. Given that the jet of J0509+1011 is likely to have some extended emission that is not detectable by the VLBA, it is likely that the absolute flux calibration of the VLBA is well within its nominal value of 5%.

The radio afterglow was detected at a position of (epoch J2000) $`\alpha =5^h9^m31.2983^s`$, $`\delta =+11^{}17^{}7.262^{\prime \prime }`$, with (conservative 1-$`\sigma `$ error of 0.001′′ in each coordinate). The source is unresolved with a size of less than 0.001′′.

Ryle Telescope: Observations at 15 GHz with the Ryle Telescope at Cambridge (UK) were made by interleaving 15 minute scans of GRB 991216 with short scans of the phase calibrator J0509+1011. The flux density scale was tied to observations of 3C 48 and 3C 286.

Owens Valley Radio Observatory Interferometer (OVRO): The source was observed for a single 13 hr track in two continuum 1 GHz bands (central frequencies 98.481 GHz and 101.481 GHz) under good 3-mm weather conditions. Gain calibration used the quasar 0528$`+`$134, while observations of Uranus and 3C 454.3 provided the flux density calibration scale with an estimated uncertainty of $`20`$%. See Shepherd et al. (1998), for details of the calibration and imaging. No source was detected (see Table 1).

James Clark Maxwell Telescope (JCMT222The JCMT is operated by The Joint Astronomy Centre on behalf of the Particle Physics and Astronomy Research Council of the UK, the Netherlands Organization for Scientific Research, and the National Research Council of Canada.): Observations in the 350 GHz band were made using the Sub-millimeter Common-User Bolometer Array (Holland et al. (1999)). The data were taken under good sky conditions on both nights. For flux calibration we used the source CRL618, and assumed its flux density to be 4.57 $`\pm `$ 0.21 Jy. The pointing was monitored and found to vary by less than 2″. See Kulkarni et al. (1999) for details of data reduction. The source was not detected at either epoch (see Table 1). At the position of GRB 991216 we derive an average flux of $`0.28\pm 1.1`$ mJy.

## 3 Results

In Figure 1 we display the 8.46 GHz light curve, as well as the X-ray and optical (R-band) light curves obtained from measurements reported in the GRB Coordinates Network (GCN)333$`\mathrm{http}://\mathrm{lheawww}.\mathrm{gsfc}.\mathrm{nasa}.\mathrm{gov}/\mathrm{docs}/\mathrm{gamcosray}/\mathrm{legr}/\mathrm{bacodine}/\mathrm{gcn}\mathrm{\_}\mathrm{main}.\mathrm{html}.`$ and Halpern et al. (2000). A noise-weighted least squares fit of the form $`F_\nu t^{\alpha _\nu }`$ was made to each of these light curves. Using all the 8.46 GHz data, including the upper limits, we derive $`\alpha _r=0.82\pm 0.02`$ ($`\chi _r^2=26.5/15`$; here $`\chi _r^2`$ is the reduced $`\chi ^2`$).

A similar least squares fit of the optical and X-ray data over the first three days (Figure 1), yields $`\alpha _o=1.33\pm 0.01`$ ($`\chi _r^2=11/28`$) and $`\alpha _x=1.61\pm 0.06`$ ($`\chi _r^2=7.7/3`$). From a more extensive dataset, Halpern et al. (2000) fit $`\alpha _o=1.07_{0.08}^{+0.17}`$ over the same time range. We will use their value of $`\alpha _o`$ in the discussion to follow. The relatively large value of $`\chi _r^2`$ for fit to the X-ray data presumably reflects the uncertainties inherent in converting the counts measured by three different instruments (RXTE-ASM, RXTE-PCA and Chandra ACIS) into Jansky flux units. Using the RXTE-PCA data alone avoids this cross-calibration issue and yields $`\alpha _x=1.61\pm 0.05`$.

## 4 The Unusual Nature of the Radio Afterglow: The Failure of the Basic Afterglow Model

The radio afterglow from GRB 991216 is unusual on two counts. First, the radio afterglow in the centimeter band does not show the usual rise to a peak value $`f_m`$ (at epoch $`t_m`$) before undergoing a power law decay. The radio flux appears to decline continuously starting from the epoch of the first observation. Thus $`t_m<1.49`$ d as compared to the 10–100 d seen in other bursts (e.g. Frail, Waxman & Kulkarni (2000), Frail et al. (1999), Frail et al. (1999)). Second, the temporal decay indices ($`\alpha _\nu `$) in the radio, optical and X-ray bands are markedly different from each other. Proceeding from radio to higher frequencies, $`\alpha _\nu `$ steepens by $``$0.4 every four decades in frequency.

In contrast, the optical and X-ray afterglow appears to find a straightforward explanation in the standard afterglow model in which a jet geometry is invoked (Halpern et al. (2000)). Below we show that the radio observations cannot be reconciled with a standard jet (or sphere) afterglow model. We then explore possible modifications to the standard model.

The simplest afterglow model is one in which the broad-band afterglow emission arises from the forward shock of a relativistic blast wave propagating into a constant density medium (Sari, Piran & Narayan (1998)). It is assumed that the electrons in the forward shock region are accelerated to a power law distribution for $`\gamma _e>\gamma _m`$, $`dN/d\gamma _e\gamma _e^p`$; here $`\gamma _e`$ is the Lorentz factor of the electrons, $`p`$ is the power law index and $`\gamma _m`$ is the minimum Lorentz factor. Gyration of these electrons in strong post-shocked magnetic fields gives rise to broad-band afterglow emission. Two modifications to this picture are routinely considered. (1) An inhomogeneous circumburst medium (specifically, $`\rho (r)r^2`$; here $`\rho `$ is the density at distance $`r`$ from the source). Such a circumburst medium is expected should GRBs originate from massive stars (Chevalier & Li (1999)). (2) A jet-like geometry for the blast wave (Sari, Piran & Halpern (1999)). This modification is motivated by the propensity of jets in astrophysical sources.

Regardless of these modifications, the broad band spectrum is composed of three characteristic frequencies: $`\nu _a`$, the synchrotron self-absorption frequency; $`\nu _m`$, the frequency at which the emission peaks (and attributed to the electrons with Lorentz factor $`\gamma _m`$), and $`\nu _c`$, the cooling frequency. Electrons radiating photons with frequency $`>\nu _c`$ cool on timescales faster than the age of the blast wave. The evolution of these frequencies is determined by the dynamics of the blast wave. The usual ordering of these frequencies at epochs relevant to the discussion here is (going from low to high frequencies) $`\nu _a`$, $`\nu _m`$ and $`\nu _c`$.

For GRB 991216 the early radio decay implies that $`\nu _m`$ is already below the centimeter radio band at 1.49 days. The steepening of the afterglow emission from optical to X-ray can be explained by placing $`\nu _c`$ between the optical and X-ray bands. The expected steeping $`\mathrm{\Delta }\alpha `$ is 1/4 which is marginally consistent with $`\alpha _o\alpha _x=0.54_{0.10}^{+0.18}`$. However, even if we ignore this, we are simply unable to explain the decay in the radio band, since no additional steepening is expected between $`\nu _m`$ and $`\nu _c`$.

The standard afterglow model can be made to agree with the light curves by postulating an energy slope $`p`$ which gradually steepens with increasing electron energy $`\gamma _e`$. We use the spherical, constant density afterglow model (Sari, Piran & Narayan (1998)) to convert, in each band, the observed decay index to $`p`$ and obtain: $`p=2.09\pm 0.03`$ (radio), $`p=2.43_{0.11}^{+0.23}`$ (optical), and $`p=2.81\pm 0.08`$ (if $`\nu _c`$ is below the X-ray band) or $`p=3.15`$ (if $`\nu _c`$ above the X-ray band). We are justified in applying the spherical model for early times ($`t<t_J`$ 2–5 d) since the jet geometry is manifested only for $`t>t_J`$ (Halpern et al. (2000)).

Curvature is both observed and modeled in the synchrotron spectra of the non-relativistic shocks from supernova remnants which are accelerating cosmic rays (e.g. Baring et al. 1999). To date, models of ultra-relativistic shocks favor a universal value of $`p`$, independent of energy (Vietri (2000), Gallant et al. (2000)), but non-linear effects have yet to be treated.

Nonetheless, the invocation of curvature in the energy distribution of the electrons cannot explain the observed broad-band spectrum (Figure 2) of the afterglow on December 18 (corresponding to 1.33 days after the burst). A plausible fit to the entire data is obtained with $`\nu _a=1.3`$ GHz, $`\nu _m=270`$ GHz and $`\nu _c=7\times 10^{16}`$ Hz and $`f_m=3.4`$ mJy; this fit is displayed by the dashed line in Figure 2. As the blast wave slows down, $`\nu _m`$ moves to lower values while preserving $`f_m`$ and thus we expect the flux in the centimeter band to rise, whereas the observed flux falls. If the afterglow has a jet-like geometry then the radio afterglow is expected to rise until the epoch $`t_J`$, and subsequently decay very slowly ($`f_\nu t^{1/3}`$) until $`\nu _m`$ passes through the centimeter band, after which we expect to see a decline similar to that seen in the optical ($`f_\nu t^{2.2}`$) (Harrison et al. (1999)). As can be seen from Figure 2, the radio observations are grossly inconsistent with these expectations, particularly the decay is much faster than $`t^{1/3}`$.

To summarize, while the optical and X-ray observations can be accounted for by a jet model, the radio observations are inconsistent with the standard model. This forces us to consider afterglow models in which the radio emission (at least in bulk) arises from a source other than the usual forward shock.

## 5 A Forward and Reverse Shock Model

The most natural explanation for two components would be an early contribution from a reverse shock followed by a forward shock element at later times. This is the explanation invoked to account for the early (1-2 day) radio emission from the afterglow of GRB 990123 (Sari & Piran (1999), Kulkarni et al. (1999) but see Galama et al. (1999)). The two bursts share several common features: in both cases, a jet has been deduced with $`t_J`$ few days, both were quite bright at gamma-ray energies and finally both had a seemingly small value of $`t_m`$ (as measured in the centimeter band). However, in the case of GRB 990123, the peak flux of the forward shock was $`f_m<260\mu `$Jy (Kulkarni et al. (1999)) and the radio light curve was dominated by the reverse shock. In contrast, the forward shock for GRB 991216 appears to be quite strong. This difference then explains the seemingly different radio light curves.

At late times (i.e. timescales greater than the duration of the burst) the flux from the reverse shock is expected to fall as $`t^{1.8}`$ (Kobayashi & Sari (2000)). In contrast, the forward shock emission rises as $`t^{1/2}`$ for $`t<t_J`$ and then slowly decays, $`t^{1/3}`$ until the $`\nu _m`$ moves into the centimeter band. Since $`t_J`$ is known from optical observations (Halpern et al. (2000)), the remaining unknowns are the strength of the reverse and forward shock emission.

In this picture, the reverse shock dominates the radio emission for the first few days and the model fit consists of mainly fitting a power law with $`f_\nu t^{1.8}`$. We note that at day 1.5, the VLA 8.46 GHz flux and the Ryle 15 GHz flux are comparable. This suggests that the reverse shock is already optically thin at 8.46 GHz at this epoch – similar to the situation for GRB 990123 (Kulkarni et al. (1999)). We deduce the parameters of the forward shock by fitting the radio to optical spectrum around $`t_J`$=5 days to the forward shock model (the contribution of the reverse shock is expected to be negligible thanks to the steep decay and since $`t`$ is comparable to $`t_J`$, the spherical fireball model is still applicable); we find $`\nu _m1.4\times 10^{12}`$ Hz and $`f_m=1`$ mJy. As can be seen from Figure 3 this reverse-forward model provides a reasonable fit to the observations.

There are two predictions of this model. First, we expect $`\nu _m`$ to cross the centimeter band at $`t_b=t_J(\nu _m/8.46\mathrm{GHz})^{1/2}64`$ d. For $`t>t_b`$, we expect the radio flux to decline as steeply as the optical flux does for $`t>t_J`$. The low flux values as measured at the VLA around this epoch are in agreement with this model.

A second prediction is that for $`t<t_b`$, we expect, the spectrum to rise as $`\nu ^{1/3}`$ for $`\nu <8.46`$ GHz. Unfortunately, the data are too sparse to rigorously test this expectation. Nonetheless, we note that at day 17.44, the spectrum between 1.43 and 8.46 GHz can be described by a simple power-law with slope $`\beta _r=0.45`$, steeper than a $`\nu ^{1/3}`$ slope by 3.6$`\sigma `$. We consider this to be the weakest point of the model but do not consider it fatal since the quoted uncertainties include only instrumental errors and do not include external effects such as interstellar scintillation.

The strongest confirmation of this model would have been the detection of an optical flash, as in the case of GRB 990123 (Akerlof et al. (1999)). The strong radio emission from the reverse shock allows us to predict (by scaling from the optical and radio observations of GRB 990123, Sari & Piran (1999)) that the flash would have been 8th magnitude or even brighter. Unfortunately, this event occurred during daytime and therefore was not observed by existing prompt optical counterpart experiments (LOTIS – H. S. Park; ROSTE – C. Akerlof; pers. comm.)

We end this section by noting a worrisome and puzzling issue: we are unable to provide a consistent explanation for the near-IR, optical and X-ray observations with a standard fireball afterglow spectrum. As noted in Figure 2, there is a broad maximum around $`2\times 10^{14}`$ Hz – suggesting that this is the peak frequency ($`\nu _m`$) of the fireball. Fitting a template afterglow spectrum we obtain the following: $`\nu _{m1}=2.1\pm 0.6\times 10^{14}`$ Hz, $`f_{m1}=150\pm 10`$ $`\mu `$Jy and $`\nu _{c1}=2\times 10^{16}`$ Hz. We note that a similar broad peak in the near IR (and attributed to $`\nu _m3\times 10^{14}`$ Hz at $`\mathrm{\Delta }t=0.5`$ d) was observed for GRB 971214 (Ramaprakash et al. (1998)). However, if we evolve this $`\nu _m`$ back in time (with $`\nu _mt^{3.2}`$) we predict a rising R-band light curve, inconsistent with the observations (Figure 1). Moving $`\nu _m`$ to lower frequencies solves this problem but we are left with no explanation for the “near-IR” bump.

## 6 A Two-Component Forward Shock Model

We now consider a model in which much of the radio emission arises as the forward shock of an additional fireball (hereafter the second fireball). The principal attraction of the second fireball is that we no longer need to relate the radio decay rate to those at optical and radio frequencies. We clarify that the optical and X-ray observations are explained by the forward shock of the fireball (the first fireball) discussed in the previous section. As noted earlier, there is good evidence suggesting that the first fireball is a jet. Thus the second fireball must be a more isotropic fireball and move at a smaller Lorentz factor. Indeed, in some GRB models, the central engine is expected to inject two fireballs: a high $`\mathrm{\Gamma }`$ jet and a low $`\mathrm{\Gamma }`$ spherical wind.

A reasonable fit to the radio data of this second fireball (FS 2; see Figure 2) on day 1.33 is provided by $`f_{m2}1.2`$ mJy, $`\nu _{m2}`$=7 GHz, and $`\nu _{a2}`$=2 GHz. The location of the cooling frequency $`\nu _{c2}`$ is unconstrained. As a test, we evolved the afterglow spectrum forward in time. The model does an excellent job reproducing the declining flux density from 1.43 and 8.46 GHz at 17.44 days (an observation which the reverse-forward shock model fails to explain), but at day 60.40 it predicts a 1.43 GHz flux of $`100`$ $`\mu `$Jy, where only an upper limit of $`57\pm 44`$ $`\mu `$Jy is measured. Again we consider this 3-$`\sigma `$ discrepancy as a major, but not fatal, weakness of this model.

The three inferred parameters ($`\nu _{m2}`$, $`f_{m2}`$, $`\nu _{a2}`$) allow us to obtain the energy of the blast wave and the density of the ambient medium (Wijers & Galama (1999)): $`E_{52}10^2`$ erg and $`n10^4`$ cm-3; these values are relatively insensitive to the value of the unknown $`\nu _c`$ (which is however constrained to lie above the optical band). The large $`E`$ and small $`n`$ are primarily due to the small value of $`t_m`$.

If this interpretation is correct then we have uncovered the first example of a GRB exploding in a very low density medium – perhaps the halo of a host galaxy. The dynamics of explosions is governed by the ratio $`E/n`$, and as noted above, this ratio is perhaps $`10^5`$ larger than that typically derived in other afterglow. For this reason, both fireballs, the high $`\mathrm{\Gamma }`$ and the low $`\mathrm{\Gamma }`$ fireballs, would then be expanding at high Lorentz factors days after the burst. The $`\mathrm{\Gamma }`$ for the low $`\mathrm{\Gamma }`$ fireball would be an impressive 20 one day after the burst, and the fireball would have had a size of 100 microarcseconds three weeks after the burst – almost within reach of measurable with VLBI techniques (cf. Taylor et al. (1997)). The jet fireball would be expanding even faster in which case the opening angle of the jet is not 6 (Halpern et al. (2000)) but only 1.

To summarize, the radio afterglow of GRB 991216 is unusual and cannot be explained by the standard forward shock model. A conventional reverse-forward shock model or an exotic two-component forward shock model can account for the observations, but each has one major (but not necessarily fatal) weakness. Finally, we have no explanation for the near-IR bump seen on day 1.33. GRB 991216 shows that there may be yet new surprises in GRB afterglows.

DAF thanks Chris Fassnacht, Steve Myers, Lin Yan and Jim Ulvestad for generously giving up portions of their VLA time so that GRB 991216 could be observed. DAF thanks J. Halpern for making his paper available prior to publication and Y. Gallant and M. Vietri for useful discussions. We would like to thank S. Jogee for preparing the first observations of this burst at the Owens Valley Observatory and A. Sargent for generously allocating the time on short notice. Research at the Owens Valley Radio Observatory is supported by the National Science Foundation through NSF grant number AST 96-13717. SRK’s research is supported by grants from NSF and NASA. RS and TJG are supported by Sherman Fairchild Fellowships.