text:
# Driving Bose-Einstein condensate vorticity with a rotating normal cloud

## Abstract

We have developed an evaporative cooling technique that accelerates the circulation of an ultra-cold 87Rb gas, confined in a static harmonic potential. As a normal gas is evaporatively spun up and cooled below quantum degeneracy, it is found to nucleate vorticity in a Bose-Einstein condensate. Measurements of the condensate’s aspect ratio and surface-wave excitations are consistent with effective rigid-body rotation. Rotation rates of up to 94$`\%`$ of the centrifugal limit are inferred. A threshold in the normal cloud’s rotation is observed for the intrinsic nucleation of the first vortex. The threshold value lies below the prediction for a nucleation mechanism involving the excitation of surface-waves of the ground-state condensate.

To paraphrase an ancient riddle, what happens when an irresistible torque meets an irrotational fluid? The answer has been known for more than 50 years: a quantized vortex is nucleated. Vortices alone contribute to a superfluid’s rotation, so that the bulk of the fluid may remain curl-free. The nucleation of vortices in bulk superfluid Helium has been the topic of extensive study (for a review see ). In the archetypical experiment, a rotatable pot filled with a mixture of superfluid and normal liquid Helium undergoes gradual angular acceleration. The normal fluid and the walls of the pot rotate together as a rigid body, defining a rotating environment. At some threshold angular velocity, a vortex line is nucleated at the circumference of the pot, and then quickly migrates inward until it is collinear with the axis of rotation. Further angular acceleration results in the nucleation of more vortices; eventually the fluid is filled with an array of vortex lines . A central theme  of this research is the question: to what extent is the nucleation process “extrinsic,” i.e. dependent on such details as the roughness of the surface of the walls, and to what extent is it “intrinsic” , i.e. driven (in the limit of microscopically smooth walls) by the flow of normal fluid along the boundary of the superfluid? In the analogous rotating-potential experiments with a dilute-gas Bose-Einstein condensate (BEC), the confining potential and the normal fluid typically rotate at different rates . In this context, the extrinsic-intrinsic question can be restated as: is it the confining potential or the normal fluid that defines the rotating environment?

Vortices in a BEC have been created with wavefunction engineering , through the decay of solitons , and in the wake of moving objects . The first rotating-potential experiment to detect vortices in a BEC was performed by the Paris group ; results have also been obtained by the MIT  and Oxford  groups. In these experiments the role of the normal fluid was secondary to that of the rotating potential; it is conceivable the normal fluid was not rotating at all. This paper presents vortex nucleation experiments performed in the opposite limit, namely in the environment of a rotating normal gas in a static confining potential. Such an environment allows for the isolated study of the intrinsic mechanism for vortex nucleation.

Our experiments begin with a magnetically trapped cloud of about $`610^6`$ 87Rb atoms, in the $`|F=1,m_F=1`$ hyperfine state, cooled close to the critical temperature $`T_c=67`$ nK. The atoms are initially confined in an axially symmetric, oblate and harmonic potential  with axis of symmetry along the vertical (‘z’) axis. To induce rotation of the cloud, we first gradually apply an elliptical deformation to the potential in its horizontal plane of symmetry, and then rotate the deformation  about the vertical axis at a fixed angular frequency. The rotating potential is characterized by an axial frequency $`\omega _z=2\pi (13.6)`$Hz, time-averaged radial frequency $`\omega _\rho =2\pi (6.8)`$Hz and a horizontal ellipticity of $`25\%`$. Such a large rotating trap asymmetry, accessible in the oblate configuration of our apparatus, is found to be necessary not only to get the cloud rotating, but also to sustain ongoing rotation. Moreover, in steady state the cloud does not reach the rotation rate of the applied asymmetry. We believe that the stirring process is fighting a small, static asymmetry that acts to despin the cloud .

In thermal equilibrium, a normal cloud rotates as a rigid body with the centrifugal force causing the cloud to bulge outwards in the radial direction. In order to detect the rotation, we use a nondestructive phase-contrast technique to image the cloud in situ from the side. Four sequential pictures of a given cloud are taken to average over oscillations in the widths and to improve the signal to noise. The cloud temperature is extracted from its vertical width $`\sigma _z`$, which is unaffected by rotation about the vertical axis. The rotation of the cloud $`\mathrm{\Omega }_N`$ is determined from the aspect ratio, $`\lambda =\sigma _z/\sigma _\rho `$, using the relation:

$$\mathrm{\Omega }/\omega _\rho =\sqrt{1(\lambda /\lambda _o)^2}$$
(1)
where $`\lambda _o`$ is the static aspect ratio. The technique of side-view imaging is crucial for distinguishing between changes in radial size due to temperature and to rotation.

For stirring rates up to 2.5 Hz, the rotation of the cloud reaches its steady-state value by 15 s or less. After 15 s, the rotating trap asymmetry is ramped off, leaving non-condensed clouds a factor of 1.2 – 1.3 above $`T_c`$, for stirring frequencies 0 – 2.5 Hz. Radio frequency (rf) evaporation is then used to cool the normal cloud to BEC. In the oblate configuration of our TOP trap , we have found that rf evaporation immediately quenches the rotation, presumably because the selection process, which removes atoms with large radial displacements, preferentially removes atoms with large axial components of their angular momentum. By adiabatically distorting the trap into a prolate geometry with $`\{\omega _\rho ,\omega _z\}=2\pi \{8.35,5.45\}`$Hz , we can instead cool the cloud by removing atoms with large axial displacements and thereby reduce the effect of the evaporation on the axial angular momentum. As the normal cloud is evaporated, its aspect ratio is observed to decrease continuously, indicating a monotonically increasing rotation rate (Fig. 1a). During the evaporation, the angular momentum per particle of the normal cloud remains roughly constant, even though the number of atoms is reduced by over a factor of five and temperature, by a factor of four (Fig. 1b). As the cloud cools and shrinks, it must spin up for the angular momentum per particle to remain fixed.

To reach significant rotation rates by the end of evaporation requires the lifetime of the normal cloud’s angular momentum to be comparable to the evaporation time. The nearly one-dimensional nature of the evaporation together with the low average trap frequencies make cooling to BEC in the prolate trap very slow ($``$50 s). We obtain angular momentum lifetimes this long by shimming the azimuthal trap symmetry to better than 0.1$`\%`$.

Towards the end of the evaporation, a condensate begins to appear at the center of the rotating normal cloud . We discuss first the results of experiments in which we continue the evaporation until little or no normal fraction remains. In this case, we find that the rotating normal component has given birth to a condensate distended in its radial dimension, as one would expect for a classical rotating body under the influence of the centrifugal force. This effect is reminiscent of liquid Helium experiments, in which the surface of a rotating bucket of superfluid exhibits the same meniscus curvature as for an ordinary viscous fluid . For large enough numbers of vortices in the condensate, the correspondence principle would suggest that the rotation field, coarse grained over the cloud, should go over to the classical limit of rigid-body rotation. In this limit, the classical Eq. (1) should connect condensate aspect ratio to rotation rate.

Alternatively, we can study the angular momentum in the BEC directly by exciting quadrupolar surface waves  with a rotating weak deformity of the magnetic

trap. The quadrupolar surfaces waves are characterized by angular momentum quantum number $`m_z=\pm 2`$ where the $`m_z=+()2`$ excitation is taken to be co-(counter-)propagating with the rotation of the condensate. By varying the initial stir rate applied to the normal cloud, condensates of different aspect ratio can be accessed. In Fig. 2a, the frequency of the $`m_z=\pm 2`$ modes is shown as a function of condensate aspect ratio. The $`m_z=+2`$ mode is seen to speed up and the $`m_z=2`$ to slow down due to the presence of vorticity in the condensate.

For small rotation rates, the splitting between the $`m_z=\pm 2`$ modes is predicted to be linearly proportional to the mean angular momentum of the condensate . In the large-$`\mathrm{\Omega }`$ limit of rigid-body rotation, Zambelli and Stringari  have used a sum-rule argument to show that the splitting between the modes is simply $`2\mathrm{\Omega }`$, and, further, that the sum of the squared frequencies of the two modes is independent of rotation rate. A best fit of this model to the combined $`m=\pm 2`$ data is shown in Fig. 2a, where the rigid-body rotation rate has been inferred from the condensate aspect ratio and the classical Eq. (1). Perhaps more intuitively, the frequency splitting is plotted explicitly versus inferred rotation rate in Fig. 2b. The excellent agreement with the model of Zambelli and Stringari is compelling evidence in favor of the reasonableness of using Eq. (1) to connect the condensate aspect ratio with its effective rotation rate. This

is further born out by extensive 3-D numerical simulations of the Gross-Pitaevskii equation for the parameters of our experiment, by Feder and Clark. Their numerical simulations confirm that a condensate in an environment rotating at frequency $`\mathrm{\Omega }>0.5\omega _\rho `$ will equilibrate close to the aspect ratio given by Eq. (1) .

In pure condensate samples we have observed aspect ratios as pronounced as 0.35 $`\lambda _o`$, corresponding to a rotation rate of $`0.94\omega _\rho `$. The rapid rotation rate, combined with the increased condensate area arising from its radial bulge, mean that the condensate must be supporting a large number of vortices. Feder and Clark  calculate 56. With these initial conditions, we have observed continued rotation for at least 140 s.

If the evaporation is stopped before the normal cloud has been completely removed, a comparison can be made between the aspect ratios, and hence rotation rates, of the condensate and normal cloud. By adjusting the initial stir rate applied to the normal cloud and the depth of the evaporation, we are able to reach different rotation rates for a given condensate fraction. After the evaporation is stopped, a time of 5 s is allowed for the gas to rethermalize to the last evaporative cut. We then take four nondestructive pictures and fit the images to a two-component distribution. Figure 3 shows a plot of the condensate aspect ratio $`\lambda _{BEC}`$ compared with the aspect ratio $`\lambda _N`$ of the normal cloud where each point represents a single realization of the experiment.

For $`\lambda _N<1.36`$ ($`\mathrm{\Omega }_N>0.5\omega _\rho `$), the condensate aspect ratio closely tracks the normal aspect ratio, providing a further manifestation of the correspondence principle for a highly rotating BEC. In the vicinity of $`\lambda _N=1.48`$ ($`\mathrm{\Omega }_N=0.35\omega _\rho `$), we observe threshold behavior in the condensate rotation. At this low value of rotation we don’t expect the rigid-body model to be valid for the condensate; but we make use of a numerical calculation by Feder and Clark  to indicate the change in condensate aspect ratio associated with the presence of a single, centered vortex (Fig. 3b). There is considerable scatter in the data so one cannot make a strong statement about the nature of the threshold shape, but clearly, somewhere between $`0.32<\mathrm{\Omega }_N/\omega _\rho <0.38`$ the first vortex is nucleated.

A comparison of this threshold location with two theoretical rotation rates for our particular cloud provides some insight into the nature of the vortex nucleation. The first value is the threshold for the thermodynamic stability of a single vortex. For our experiment this number is 0.2 – 0.25 $`\omega _\rho `$ , distinctly lower than our observed threshold for nucleation. The second value is $`\omega _{min}`$, the frequency at which the slowest surface-wave mode propagates around the circumference of the condensate. The Paris group has shown that for their ”extrinsic” nucleation process (vortices nucleated by a rotating asymmetric potential) the key mechanism is the nonlinear excitation of surface waves . Results from MIT  and from Oxford  are also consistent with such a mechanism. For the parameters of our experiment, $`\omega _{min}=0.4\omega _\rho `$ . Our observed threshold for “intrinsic” nucleation is clearly under this value, and thus we cannot interpret our effect in terms of a normal “wind” exciting surface waves on our condensate.

The data presented in Fig. 3 include a range of condensate fraction from 0.1-0.6 for each rotation rate of the normal cloud, although no segregation relative to either axis is evident for plots of different BEC fractions. Moreover, by reducing the rate of evaporation, we have decreased the rate of acceleration of the normal cloud rotation by a factor of 3, and still observe a threshold for vortex formation between 0.32 and 0.38 (Fig. 3b).

In future work we plan to use the controlled rotation of the normal cloud to study the life cycle of a BEC’s vorticity , including any hysteresis between spin-up and spin-down due to vortex metastability. By halting the rotation of the normal cloud, we have already found that we can reduce the lifetime of the BEC’s rotation by a factor of 10. Finally, the rotating normal cloud can create equilibrated condensates with very large rotation rates, which may allow us to approach the regime for which the vortices are so close-packed that their separation becomes comparable to the healing length .

We acknowledge very useful conversations with Seamus Davis, David Feder, Sandro Stringari, and Carl Wieman. This work was supported by NSF and NIST funding. P. E. acknowledges support from the Alexander-von-Humboldt foundation.


format: text

text:
# Electronic Structure near Impurities in the Chains of YBa2Cu3O6+x

\[

## Abstract

We study the electronic structure near impurities in the one-dimensional chains of YBa2Cu3O6+x. Assuming chain superconductivity below $`T_c`$ due to a proximity coupling, we show that only a magnetic impurity induces resonance states in the local density of states. Its spatial structure reflects the particle-hole nature of chain superconductivity and thus distinguishes it from other broken symmetry phases. For a sufficiently large scattering strength, or due to interference effects arising from a second impurity, the chains undergo a quantum phase transition to a polarized state.

\]

Studying the electronic structure of the CuO2 planes in YBa2Cu3O6+x (YBCO), an important step towards our understanding of superconductivity, has been complicated by the presence of one dimensional (1D) CuO-chains in this material. The coupling of the CuO-chains to the CuO2-planes changes the electronic structure of the latter and might facilitate the occurrence of other 1D structures, e.g., stripes, in the planes. Much theoretical work  has focused on this coupling and its implications for a series of experimental probes ranging from inelastic neutron scattering , nuclear quadrupole resonance (NQR) , infrared spectroscopy , angle-resolved photoemission  and scanning tunneling microscopy (STM) experiments . In turn, the same coupling necessarily changes the chain electronic structure when the CuO2 planes become superconducting (SC), as was recently evidenced by NQR  and STM  experiments. In particular, NQR experiments observe an increase in the 63Cu NQR line width near the SC transition , while STM experiments  find evidence for large oscillations in the local density-of-states (DOS). Whether these observations arise from the presence of a charge-density wave (CDW) , or from impurities in SC chains is currently a topic of intense discussion.

The study of the electronic structure near impurities in 2D superconductors  has been a successful method to distinguish between various broken symmetry ground states. In this Letter we extend these works to study the effects of chain impurities, both magnetic and non-magnetic, on the local DOS assuming that the chains become SC due to a proximity coupling. Both types of impurities induce large spatial oscillations in the DOS, and thus the electron density, whose amplitude increases below $`T_c`$. However, in contrast to a non-magnetic impurity, a magnetic impurity induces two resonance states in the chain DOS. Moreover, we present a prediction for the spatial dependence of these resonances which reflects the particle-hole nature of chain superconductivity. This prediction therefore discriminates between oscillations arising from a CDW and those due to impurities. Finally, we find that as the distance between two magnetic impurities is reduced, the resonance states split and the system can be driven through a quantum phase transition. Our work is complementary to that of Eggert  who studied the effects of impurities in a 1D Luttinger liquid.

Since in the following we employ a model described in detail in Ref. , we only briefly review it here. The mean-field BCS-Hamiltonian of the chain-plane system is given by

$``$ $`=`$ $`{\displaystyle \underset{k,\sigma }{}}ϵ_𝐤c_{k,\sigma }^{}c_{k,\sigma }{\displaystyle \underset{k,\sigma }{}}\zeta _𝐤d_{k,\sigma }^{}d_{k,\sigma }t_{}{\displaystyle \underset{k,\sigma }{}}c_{k,\sigma }^{}d_{k,\sigma }`$ (2)
$`+{\displaystyle \underset{k}{}}\mathrm{\Delta }_kc_{k,}^{}c_{k,}^{}+h.c.`$
where $`c_k^{},d_{k_x}^{}`$ are the fermionic creation operators in the plane and chain, respectively. One has for the planar and chain tight-binding dispersions, $`ϵ_𝐤`$ and $`\zeta _𝐤`$,

$`ϵ_𝐤`$ $`=`$ $`2t_p\left(\mathrm{cos}(k_x)+\mathrm{cos}(k_y)\right)4t_p^{}\mathrm{cos}(k_x)\mathrm{cos}(k_y)\mu _p,`$ (3)
$`\zeta _𝐤`$ $`=`$ $`2t_c\mathrm{cos}(k_x)\mu _c.`$ (4)
Here, $`t_p=300`$ meV, $`t_p^{}=0.22t_p,\mu _p=0.812t_p`$ are the planar hopping elements and the chemical potential , and $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }_0(\mathrm{cos}(k_x)\mathrm{cos}(k_y))/2`$ is the planar SC gap with $`\mathrm{\Delta }_040`$ meV . The chain Fermi momentum, $`k_F^c0.25\pi `$,  yields $`\mu _c=1.41t_c`$, and we take $`t_{}=0.4t_p`$ . After integrating out the planar fermions, the chain Greens function $`\widehat{G}_0=𝒯\mathrm{\Psi }_{k_x,l}(\tau )\mathrm{\Psi }_{k_x,m}^{}(\tau ^{})`$ with

$$\mathrm{\Psi }_{k_x,l}^{}=(d_{k_x,l,}^{},d_{k_x,l,}^{},d_{k_x,l,},d_{k_x,l,}),$$
(5)
and $`m,l`$ chain indices, is given by

$$\widehat{G}_0^1=\left(\begin{array}{cc}B(k_x,\omega _n)\widehat{\sigma }_0& C(k_x,\omega _n)\widehat{\sigma }_0\\ C(k_x,\omega _n)\widehat{\sigma }_0& B(k_x,\omega _n)\widehat{\sigma }_0\end{array}\right),$$
(6)
where $`\widehat{\sigma }_0`$ is the unit matrix in spin space, and

$`B`$ $`=`$ $`\left(i\omega _n\zeta _𝐤\right)\delta _{l,m}`$ (8)
$`t_{}^2N^1{\displaystyle \underset{k_y}{}}{\displaystyle \frac{i\omega _n+ϵ_𝐤}{(i\omega _n)^2(ϵ_𝐤)^2\mathrm{\Delta }_k^2}}e^{ik_y(lm)};`$
$`C`$ $`=`$ $`t_{}^2N^1{\displaystyle \underset{k_y}{}}{\displaystyle \frac{\mathrm{\Delta }_k}{(i\omega _n)^2(ϵ_𝐤)^2\mathrm{\Delta }_k^2}}e^{ik_y(lm)}.`$ (9)
In contrast to previous work , we showed  that the experimentally measured clean DOS can only be obtained qualitatively as well as quantitatively by making the ansatz

$$\widehat{G}_c^1(k_x,l,m,\omega _n)=\widehat{G}_c^1(k_x,\omega _n)\delta _{l,m}.$$
(10)
We argued that the physical origin of Eq.(10), which describes the absence of coherent correlations between the chains, or equivalently, the coupling of a single chain, lies in the presence of planar stripes and/or doping inhomogeneities in the chains.

We employ the $`\widehat{T}`$-matrix formalism  to obtain the local DOS

$$N(\omega )=\frac{2}{\pi }\mathrm{Im}\left[\widehat{G}^{(11)}(𝐫,𝐫,\omega _+)+\widehat{G}^{(22)}(𝐫,𝐫,\omega _+)\right],$$
(11)
in the vicinity of an impurity, with $`\omega _+=\omega +i0^+`$. Here, one has for the full chain Greens function

$`\widehat{G}(𝐫,𝐫^{},i\omega _n)`$ $`=`$ $`\widehat{G}_0(𝐫𝐫^{},i\omega _n)+`$ (13)
$`\widehat{G}_0(𝐫,i\omega _n)\widehat{T}(i\omega _n)\widehat{G}_0(𝐫^{},i\omega _n),`$
where $`\widehat{T}(i\omega _n)`$ represents all scattering processes and $`\widehat{G}_0`$ is obtained from Eqs.(6) and (10). In what follows we consider quite generally an impurity which is both a magnetic and potential scatterer. Assuming for simplicity that the spin of the impurity is aligned along the $`z`$-axis  we obtain for the scattering Hamiltonian

$$_{sc}=\underset{k_x,k_x^{},m}{}\mathrm{\Psi }_{k_x,m}^{}\widehat{V}\mathrm{\Psi }_{k_x^{},m}$$
(14)
with $`\widehat{V}=\widehat{\tau }_3\left(U_0\widehat{\sigma }_0JS_z\widehat{\sigma }_3\right)`$ and $`\tau _i,\sigma _i`$ are the Pauli-matrices in Nambu and spin-space, respectively. $`U_0(J)`$ is the potential (magnetic) scattering strength and $`S_z`$ is the $`z`$-projection of the impurity spin. One then finds

$$\widehat{T}(i\omega _n)=\left[1\widehat{V}\widehat{G}_0(r_x=0,i\omega _n)\right]^1\widehat{V}.$$
(15)
We remind that an impurity bound state exists at a frequency $`\omega _{res}`$, if $`det[\widehat{T}^1(\omega _{res})]0`$, while for a resonance state only the real part of $`det[\widehat{T}^1]`$ vanishes.

We first study the case of a purely magnetic impurity with $`U_0=0`$ and $`J0`$. A straightforward analysis of Eq.(15) shows that there exist two resonance states whose frequencies,$`\pm \omega _{res}`$, are given by the solution of

$$\frac{\omega _{res}}{\alpha \overline{\mathrm{\Delta }}}\left[\sqrt{\overline{\mathrm{\Delta }}^2\omega _{res}^2}+\alpha \right]=\pm \frac{\gamma ^21}{\gamma ^2+1}$$
(16)
where $`\gamma =v_F^c/\beta `$, $`v_F^c`$ is the chain Fermi velocity and $`\beta =JS_z/2`$. Moreover, $`\alpha =t_{}^2/\overline{v}`$, $`\overline{v}=\left(ϵ_𝐤/k_y\right)|_{𝐩_0}`$, where $`𝐩_0`$ is the momentum on the planar Fermi surface (FS) with $`p_0^x=k_F^c`$ (see inset of Fig. 1) and $`\overline{\mathrm{\Delta }}=\mathrm{\Delta }(𝐩_0)`$. While in the limit $`\gamma 1`$, the resonance frequencies are located close to $`\pm \overline{\mathrm{\Delta }}`$, $`\omega _{res}`$ shifts to lower energies with increasing $`\gamma `$. At the critical $`\gamma _c=1`$ (i.e, $`\beta _c=v_F^c`$) the resonance states cross at $`\omega _{res}=0`$ and the magnetic impurity forms a bound state with a spin-up (down) quasi-particle for $`J<0(>0)`$. Accordingly, the ground state of the system undergoes a quantum phase transition and changes from an unpolarized one for $`\gamma <\gamma _c`$ to a spin-polarized one with total spin $`S=1/2`$ for $`\gamma >\gamma _c`$. A similar transition also occurs in superconductors with $`s`$-wave symmetry .

In Fig. 1 we present the SC chain DOS, obtained from a numerical evaluation of Eq.(13), in the clean case and at the site of a magnetic impurity with $`\beta =200`$ meV and $`U_0=0`$. A comparison with the clean DOS shows the emergence of two resonance states at the impurity site with $`\omega _{res}\pm 2.6`$ meV, in agreement with the discussion above. Simultaneously, the impurity suppresses the coherence peaks at $`\omega _{cp}=\pm 16.0`$ meV which are present in the clean DOS. Moreover, the impurity induces spatial oscillations in the chain DOS, which we present in Fig. 2 for the two resonance states and the coherence peaks at the gap edges (see inset). The oscillations at $`\pm \omega _{res}`$ possess a broad range of wave-vectors with $`0k_{osc}1.45`$, which gives rise to the rapid decay of the oscillations’ amplitude with distance from the impurity. The decay length, $`\xi 1015a_0`$ is consistent with the experimentally obtained value of $`40\AA `$ . Note that the oscillations at $`\pm \omega _{res}`$ are out-of-phase, such that a maximum in the DOS for $`+\omega _{res}`$ coincides with a minimum for $`\omega _{res}`$, as is particularly apparent when one plots the spectral weight of the resonance states at their peak positions (see inset (b)). The complementary intensity pattern in real space of the hole-like and particle-like resonance states is a reflection of the particle-hole nature of Bogoliubov quasi-particles . This characteristic feature of chain superconductivity is qualitatively different from the response of a CDW, and hence discriminates between proposed broken symmetry phases in the chains. In contrast, the oscillations at $`\pm \omega _{cp}`$, the frequency of the coherence peaks (see inset), are in-phase. Moreover, their amplitude decays much slower with distance from the impurity site since the wave-vectors of these oscillations are centered in a narrow range, $`\mathrm{\Delta }k0.02`$, around $`k_{osc}=0.81`$ (see arrow in the inset of Fig. 1). The wavelength, $`\lambda `$, of these oscillations exhibits an asymmetric frequency dependence with $`\lambda (\omega _0)>\lambda (\omega _0)`$ for $`\omega _0>0`$ and a minimum at $`\omega 0`$, in agreement with STM experiments . The spatial dependence of the DOS oscillations is qualitatively robust against changes in the form of the planar FS. However, their amplitude increases in the SC state, and accordingly, the spatial electron density changes. While this result is consistent with the broadening of the NQR line in the SC state , a quantitative analysis requires a microscopic model for the coupling of the electron density to the quadrupolar moment which is beyond the scope of this Letter.

We next consider a purely non-magnetic impurity with $`J=0`$ and $`U_00`$. We remind that a potential scatterer only induces a resonance state if the SC phase varies along the FS, or between Fermi points. However, the two $`k_F`$-points in the chain possess the same SC phase, which leads to $`det[\widehat{T}^1]1+(U_0/v_F^c)^2`$, and no resonance states exist, in analogy to the $`s`$-wave superconductor. Our numerical evaluation supports this conclusion revealing no resonances but again large spatial oscillations in the DOS.

Since an impurity is in general a magnetic and potential scatterer we next study the case $`J,U_00`$ and present in Fig. 3 the chain DOS at the impurity site for $`\beta =200`$ meV and three values of $`U_0`$.

We find that with increasing $`U_0`$, the position of the resonances shifts up in energy from $`\omega _{res}=2.6`$ meV for $`U_0=0`$ to $`\omega _{res}=9.5`$ meV for $`U_0=200`$ meV. This implies, that the critical value $`\gamma _c`$ (and thus $`\beta _c`$) where the ground state of the system becomes spin-polarized increases, again similar to the case of an $`s`$-wave superconductor . Moreover, while the spatially integrated spectral weight of the particle-like resonance is $`1.43`$ times larger than that of the hole-like resonance for $`U_0=0`$, the ratio is $`0.77`$ for $`U_0=200`$ meV. Thus spectral weight is transferred from the particle-like to the hole-like resonance with increasing $`U_0`$, however, their phase-shifts remain unchanged.

Due to the 1D character of the chains, interference effects between two impurities should be large and significantly change the single-impurity DOS discussed above. To study these effects, we employ a $`\widehat{T}`$-matrix formalism for two identical impurities where

$`\widehat{G}(𝐫,𝐫^{},i\omega _n)`$ $`=`$ $`\widehat{G}_0(𝐫,𝐫^{},i\omega _n)+`$ (18)
$`{\displaystyle \underset{n,m=1}{\overset{2}{}}}\widehat{G}_0(𝐫,𝐫_m,i\omega _n)\widehat{T}(𝐫_m,𝐫_n,i\omega _n)\widehat{G}_0(𝐫_n,𝐫^{},i\omega _n),`$
$`r_{1,2}`$ are the positions of the impurities, and $`\widehat{T}(𝐫_m,𝐫_n,i\omega _n)`$ obeys the Bethe-Salpeter equation

$`\widehat{T}(𝐫_m,𝐫_n,i\omega _n)`$ $`=`$ $`\widehat{V}\delta _{r_m,r_n}`$ (20)
$`+\widehat{V}{\displaystyle \underset{l=1}{\overset{2}{}}}\widehat{G}_0(𝐫_m,𝐫_l,i\omega _n)\widehat{T}(𝐫_l,𝐫_n,i\omega _n).`$
In general, we expect that the chain DOS for two impurities is similar to that of a single impurity if the inter-impurity distance, $`\mathrm{\Delta }r=r_2r_1`$, is larger than the spatial extent of single impurity effects, $`\xi 1015`$ at $`\omega _{res}`$ (see Fig. 2). In contrast, for $`\mathrm{\Delta }r<\xi `$ significant interference effects from multiple scattering processes should occur. To test this conjecture, we present in Fig. 4 the numerically obtained chain DOS at one of the two impurity sites as a function of $`\mathrm{\Delta }r`$.

For $`\mathrm{\Delta }r=\mathrm{}`$ (Fig. 4a, solid line) we recover as expected the DOS for a single impurity with two resonance states at $`\omega _{res}=\pm 2.6`$ meV (see Fig. 1). As the distance between the impurities decreases, the interference effects arising from multiple scattering processes lead to the formation of bonding and anti-bonding resonances state, and thus to a splitting of the two resonances into four ($`\mathrm{\Delta }r=20`$, dashed line). This result provides a possible explanation for the experimental observation of four resonance states in the chain DOS . The splitting first becomes discernible when $`\mathrm{\Delta }r2\xi `$, in agreement with our conjecture above. The splitting increases until the inter-impurity distance reaches $`\mathrm{\Delta }r=15`$ (dotted line). For even smaller $`\mathrm{\Delta }r=10`$ (Fig. 4b, solid line), the peaks recombine again and spectral weight is transferred to new resonance states at higher energies ($`\mathrm{\Delta }r=5`$, dashed line). These new states are precursors of the resonances one obtains in the limit $`\mathrm{\Delta }r0`$ (dotted line) when the DOS is that of a single impurity with twice the scattering strength, $`\stackrel{~}{\beta }=400`$ meV. Note that $`\beta =200`$ meV is below the critical value, $`\beta _c`$, to induce a transition of the ground state, while $`\stackrel{~}{\beta }=400`$ meV is above this critical value. Thus, the system can be driven through a quantum phase transition from an unpolarized to a polarized state by continually varying the distance between two or more impurities.

In conclusion we have studied the electronic structure near impurities in the 1D chains of YBa2Cu3O6+x. Assuming that the chains are SC below $`T_c`$ due to a proximity coupling, we find that a magnetic impurity induces two resonance states in the local DOS. The spatial oscillations of these resonances reflect the particle-hole nature of the chain superconductivity, and are therefore characteristic features which distinguish it from other broken symmetry phases. The amplitude of these oscillations increases below $`T_c`$, which is consistent with the experimentally observed broadening of the NQR line width . Moreover, we argue that the system undergoes a quantum phase transition at a critical $`\beta _c`$, where the ground state changes from being unpolarized for $`\beta <\beta _c`$ to a polarized one for $`\beta >\beta _c`$. Thus, the chain SC state is in many aspects similar to that of an $`s`$-wave superconductor. Finally, we show that in the presence of a second impurity, interference effects due to multiple scattering processes lead to a splitting of the resonances.

We would like to thank J.C. Davis and A. de Lozanne for stimulating discussions. This work has been supported by DOE at Los Alamos.


text:
# 1 Some Typos and Improvements

## 1 Some Typos and Improvements

In  a constructive procedure for factorizing any $`SU(2)`$ matrix $`S`$ in the form:

$$S=\mathrm{\Pi }_{k=1}^Qe^{a_kA+b_kB}$$
(1.1)
with: i) O1: $`a_k>0`$; and ii) O2: $`b_kC`$ for an a priori prescribed bound $`C`$, was given. The only restrictions on the matrices $`A`$ and $`B`$ in the published version of  was that they be i) linearly independent; and ii) expressible as $`A=aX+bY`$ and $`B=cX+dY`$ for any two orthogonal matrices $`X`$ and $`Y`$ in $`su(2)`$ (the set of all $`2\times 2`$ skew-Hermitian matrices with zero trace). Here orthogonality is with respect to the inner product $`<A,B>=\text{T}race(AB^{})`$ (equivalently, the vectors in $`R^3`$ corresponding to $`A`$ and $`B`$ are orthogonal with respect to the usual inner product on $`R^3`$). There was no requirement that the vectors $`(a,b)`$ and $`(c,d)`$ in $`R^2`$ should be orthogonal.

The first purpose of this note is to observe that the only requirement that $`A`$ and $`B`$ need satisfy is linear independence. Indeed, given two linearly independent such $`A`$ and $`B`$ one can always find (constructively) a unitary transformation $`VU(2)`$ such that $`VAV^{}`$ and $`VBV^{}`$ are expressible as possibly non-orthogonal linear combinations of $`i\sigma _x`$ and $`i\sigma _y`$. We leave the technical details to interested readers \- there are several choices for such $`V`$’s. It is useful, in this regard, to note that for any $`WSU(2)`$, the transformation $`R_V:su(2)su(2),R_V(A)=VAV^{}`$ is in $`SO(3)`$ when $`su(2)`$ is identified with $`R^3`$. It then remains to find an $`SO(3)`$ matrix which simultaneously nulls out the $`i\sigma _z`$ component of $`A`$ and $`B`$. This can be easily found, for instance, by modifying the construction behind Givens rotations. Now any one of the two $`SU(2)`$ matrices corresponding to this $`SO(3)`$ matrix will achieve the desired purpose. Further as $`V`$ can be found completely constructively, this coordinate change is known. So indeed, one can proceed by assuming that $`A=aX+bY`$ and $`B=cX+dY`$ for any two orthogonal matrices $`X`$ and $`Y`$ in $`su(2)`$ and still obtain a fully explicit technique for preparing any desired unitary transformation in $`SU(2)`$ via controls whose pulse area (or power, depending on the interpretation of the $`b_k`$) is bounded a priori. The existence of such matrices, $`V`$, was missed due to an oversight in  and thus the paper concluded (incorrectly) that the general case of arbitrary linearly independent $`A,B`$ perhaps needed more work. Indeed, all of the above ingredients were already present in  and the related paper, . The paper, , also provides an explicit $`V`$ which can be used to null out the $`i\sigma _x`$ component, but without relating it to $`R_V`$ or Givens rotations.

One additional “improvement” that is immediate in  is that the values of the number of factors, $`Q`$, given in Table I on Page 7 can obviously be lowered, in many instances, by concatenating exponentials of matrices which are constant multiples of each other.

The second purpose of this note is to rectify certain typographical errors which seem to have creeped in during the typesetting process. The main error is that $`\sigma _z`$ seems to have been replaced by $`\sigma _x`$ at several points in the published version. These junctures are as follows:

1. Everywhere on Page 4, except the headings for subsections 1 and 2 and the headings for Algorithms I and II on Page 4, $`\sigma _x`$ should read as $`\sigma _z`$.
2. On Page 5, everywhere in Algorithm II (continued from Page 4) $`\sigma _x`$ should be replaced by $`\sigma _z`$.
3. On Page 5, in the heading for Algorithm III, $`\sigma _x`$ should be replaced by $`\sigma _z`$.
4. In addition, everywhere on Page 3, $`\sigma _a`$ should be replaced by $`\sigma _z`$.


text:
# The Stellar Parameters and Evolutionary State of the Primary in the d’–Symbiotic System StH𝛼1901,2

## 1 INTRODUCTION

The object StH$`\alpha `$190 was discovered during an objective-prism survey for emission-line objects between -25 and 10 from the galactic plane by Stephenson (1986), who classified it as a symbiotic star. It is known that symbiotic stars are a class of interacting binary systems consisting of a cool, evolved star and a hotter companion star, with periods ranging from 1-3 years, up to 10-100 years (e.g. Allen 1982). Based on infrared colors, symbiotics can be divided into three groups: s-, d-, and d$`^{^{}}`$-types (Webster & Allen 1975; Allen 1982). Those that show dust continuum emission between wavelengths of 1.0 and 5.0$`\mu `$m are classified as d-type and these systems contain a mass-losing Mira variable as a cool component. Those with a normal stellar spectrum in the infrared are classified as s-types. Typically, the shorter period symbiotics are found among the s-types while the longer period systems are found among the d-types. The possible binary properties of the d$`^{^{}}`$-types are less known; for example, their binary periods (assuming they are binaries, as are the d- and s-types) are completely unknown. The d$`^{^{}}`$ stars also exhibit infrared dust emission, however, their dust color-temperatures are significantly lower than in the d-types. A more complete description about the known nature and properties of d$`^{^{}}`$-types is given in Allen (1982) and Schmid & Nussbaumer (1993).

Schmid & Nussbaumer (1993) studied both optical and ultraviolet low-resolution spectra of StH$`\alpha `$190 and classified it as a d$`^{^{}}`$-type symbiotic. More recently, Munari et al. (2001) have discovered that the cool component of StH$`\alpha `$190 is a rapidly rotating G-giant (with $`vsini`$$``$ 105 km s-1). They also found that the system exhibits high-speed, bipolar mass outflow. In this $`Letter`$ we report on an optical high-resolution spectroscopic analysis of StH$`\alpha `$190, with the principal goal being to derive basic parameters (effective temperature, surface gravity, chemical composition, evolutionary state, and radial velocity) of this symbiotic star.

## 2 OBSERVATIONS

The high-resolution spectra of StH$`\alpha `$190 analyzed here were obtained on two different telescopes. Four spectra were taken with the 2.1m telescope of McDonald Observatory with the Cassegrain cross-dispersed Sandiford echelle spectrometer plus CCD detector (McCarthy et al. 1993). In addition, five high-resolution spectra were obtained on the 1.52m telescope at the European Southern Observatory (ESO) with the fiber-fed coudé cross-dispersed echelle spectrometer FEROS and CCD detector (Kaufer et al. 1999). The McDonald spectra were taken at a two-pixel nominal resolving power of R$``$60,000 and covering wavelengths from 6100Å to 7900Å, while the ESO FEROS spectrum has R=48,000, and with wavelength coverage from $``$4000-9000Å. A summary of all observations is shown in Table 1.

All McDonald spectra were reduced in a standard way with the IRAF software package. The echelle frames were background corrected and flat-fielded. The ESO spectra were reduced using the online pipeline reduction (Kaufer et al. 2000). Final S/N ratios were evaluated by the measurement of the rms flux fluctuations, with values ranging from S/N= 80-150 in a continuum region at 6110Å. Spectra were checked for wavelength stability by using telluric O2 lines near 6275Å as a reference and wavelength differences were less than 6mÅ ($``$ 0.3 km s-1) from spectrum to spectrum.

## 3 ANALYSIS AND DISCUSSION

### 3.1 The Rotation of the Yellow Primary of StH$`\alpha `$190

Perhaps the most surprising aspect of StH$`\alpha `$190 is the presence of large numbers of absorption features broadened by a rather large rotational velocity (see also Munari et al. 2001). The wavelengths and strengths of many of these features coincide with absorption lines observed in late-F to early-K stars, e.g., Fe i or Ca i. For stars within this spectral range, high-resolution spectra can often be used to derive stellar parameters, as well as fairly detailed chemical compositions for a range of chemical elements. Due to the high rotational velocity of StH$`\alpha `$190, however, the obtainable spectral resolution is degraded, thus preventing an analysis based on measurements of individual equivalent widths. Instead, spectrum synthesis techniques must be used, as it is difficult to isolate features produced by only one species. Because of the difficulty in analyzing a range of different elements, in this Letter we concentrate on a small number of species (Fe i, Fe ii, Ca i, and Ba ii) to derive stellar parameters (Teff and log g) and rough chemical abundances in order to place constraints on the StH$`\alpha `$190 system. We confine a detailed analysis of selected absorption features to the red part of the spectrum for two reasons. First, the absorption line-density decreases considerably from blue to red, allowing for a clearer analysis of, as much as possible, nearly unblended spectral lines. Second, the possibility of light contamination (or veiling) from a hot companion, or hot gas in the system, diminishes towards the red. In particular, Schmid & Nussbaumer (1993) discussed light contamination of the cool component’s spectrum in StH$`\alpha `$190 and concluded that there was some contamination in the B-band, but very little by the V-band.

The rotationally broadened lines of StH$`\alpha `$190 can be seen in the two panels of Figure 1, where sample spectral regions from the ESO spectrum taken on Oct. 27, 1999 are displayed. The regions shown near 6150Å and 6430Å (as well as another region not illustrated near 5850Å) were selected to be synthesized because they contain a number of broadened features whose respective strengths are dominated approximately by single species, such as Fe i, Fe ii, Ca i, and Ba ii. Synthetic spectra covering the regions shown in Figure 1 (the smooth curves) were computed using a recent version of the LTE synthesis code MOOG (Sneden 1973), a Kurucz & Bell (1995) linelist with solar gf-values, and Kurucz ATLAS9 model atmospheres. The best-fit rotational velocity was determined by least-squares comparisons between observed and synthetic spectra. The top panel of Figure 2 shows values of $`\chi ^2`$ versus $`v\mathrm{sin}i`$. The lowest value found for $`\chi ^2`$ is used to normalize the larger values. It is clear from Figure 2 that there is a sharp minimum near $`v\mathrm{sin}i`$= 100$`\pm `$10 km s-1.

This fact can already be used to provide some constraints on stellar parameters if the assumption is made that the star is rotating at less than its critical breakup velocity, where

$$V_{\mathrm{critical}}=\sqrt{\frac{GM}{R}}=[GMg]^{1/4},$$
(1)
with $`M`$ being the stellar mass, $`R`$ its radius, and $`g`$ its surface gravitational acceleration. The bottom panel of Figure 2 shows this critical velocity as a function of stellar mass for a series of surface gravities (shown as log g), as well as the projected rotational velocity of StH$`\alpha `$190, shown as the dashed line with arrows at either end. The arrows indicate that this is a lower value for the rotational velocity of StH$`\alpha `$190, which could be considerably larger depending on the inclination angle. It is clear from the bottom panel of Figure 2 that, for any reasonable stellar mass, the surface gravity of StH$`\alpha `$190 must be larger than about log g $``$ 1.5, or the star would be rotating in excess of its critical velocity.

### 3.2 Stellar Parameters, Chemical Composition, & the Barium-Symbiotic Connection

The atmospheric parameters, effective temperature (Teff) and surface gravity (log $`g`$), along with \[Fe/H\], \[Ca/H\], and \[Ba/H\] (where \[X/H\]=log(N(X)/N(H))\- log(N(X)/N(H))) abundances were determined using spectrum synthesis. The solution for Teff and log g was obtained from a selection of 7 features near 6150Å, or 6430Å (and illustrated in Figure 1), as well as 5850Å, where a single species dominated so that they behaved almost as “single” spectral lines. Synthetic spectra covering a range in Teff (from 4500 - 6500K), log g (from 2.0 to 4.0), and microturbulent velocity (from $`\xi `$= 1.0 to 3.0 km s-1) were computed and then compared to the observed spectra. A least-squares comparison was used to select the best overall fits, with sample results illustrated for one particular set of models shown in the upper panel of Figure 3: results from a set of models with log g= 3.0 and $`\xi `$= 2.0 km s-1 are presented for various effective temperatures. A set of models with fixed gravity and microturbulence is shown because the largest changes in the chemical abundances are found for changes in Teff. Note that in this case, both Fe i and Fe ii yield similar abundances for an effective temperature of 5300K. In all observed-synthetic comparisons over these parameter ranges, the Fe, Ca, and Ba abundances were allowed to vary in order to maximize the best fits.

Reinspection of Figure 1 shows that the Fe ii features are quite sensitive to surface gravity. Estimates of the uncertainties in the derived stellar parameters would indicate that log g can be constrained to about $`\pm `$0.5 dex, while Teff is uncertain by about $`\pm `$150K. We note that an earlier estimate of effective temperature for this star was given by Schmid & Nussbaumer (1992), based upon the dereddened (B-V), as Teff= 5150K. The Teff derived here would suggest a spectral type of about G4 (Schmid & Kaler 1982), compared to Munari et al.’s estimate of G7 and Mürset & Schmid’s estimate of G5. Based on the Schmid & Kaler calibration, the various estimated spectral types correspond to a range in Teff from 5025–5300K.

It was found that the best overall fits can be obtained, based upon the Fe ii features, with a gravity near log g= 3.0 and $`\xi `$= 2.0 km s-1 (and these are the values used to construct the top panel of Figure 3). These fits also provide near-solar Fe and Ca abundances (Ca i is not plotted as it follows closely Fe i). Barium is an element produced primarily via s-process neutron-capture nucleosynthesis (e.g., Wallerstein et al. 1997). Note that in all models shown in the top panel of Figure 3, the Ba ii 6141Å and 5853Å features provide \[Ba/H\] abundances considerably larger than \[Fe/H\] (or \[Ca/H\]). These results indicate that the G-type component in StH$`\alpha `$190 is enriched in Ba, relative to Fe, with \[Ba/Fe\]$``$ +0.5. This enhancement in StH$`\alpha `$190 is most easily understood as resulting from mass transfer in the system when the white-dwarf member was a thermally-pulsing asymptotic giant branch (AGB) star, and was s-process enriched during its third dredge-up phase of stellar evolution. s-Process enhancements have been found in other symbiotic systems, such as AG Dra (Smith et al. 1996), He2-467 (Pereira et al. 1998), or S32 and UKS-Cel (Mürset & Schmid 1999).

With derived stellar parameters of Teff= 5300K and log g= 3.0, the current primary of StH$`\alpha `$190 can be located in a log g – Teff plane and compared to evolutionary tracks from stellar models. This comparison is shown in the bottom panel of Figure 3, with evolutionary tracks taken from Schaerer et al. (1993) and consisting of models computed with solar abundances. The parameters derived here suggest that the yellow component of StH$`\alpha `$190 consists of a star with M$``$ 2.5M nearing the base of the red giant branch. The values of mass, Teff, and g can be used to estimate a luminosity of L$``$ 45L for the cool component of StH$`\alpha `$, or MV$``$ 0.65. Adopting the V-magnitude from Munari et al. (2001) of V=10.5 and E(B–V)=0.10 (Schmid & Nussbaumer 1993), a distance of $``$780 pc is derived. Uncertainties in the stellar mass, gravity, and Teff lead to an uncertainty in this distance of about a factor of two (d$``$ 350–1500 pc). This distance in fair agreement with Munari et al.’s estimate of 575 pc.

### 3.3 Radial Velocities and the Binary Nature of StH$`\alpha `$190

Radial velocities are also obtained from the comparisons of observed to synthetic spectra for StH$`\alpha `$190, which provided fairly accurate velocity estimates in spite of the line broadening. With nine spectra obtained over the time interval from 1995 to 2001, nine photospheric radial velocities are presented here for the G-star component in StH$`\alpha `$190. The heliocentric velocities are listed in Table 1 and cover a small range from 0 km s-1 to +15 km s-1. Munari et al. (2001) have also published radial velocities that overlap the values reported here and range over some 40 km s-1 in amplitude. Both sets of radial velocities are on the same scale, as demonstrated by comparisons of interstellar Na D-line component velocities published by Munari et al. (2001) and those observed here (differences are, this study – Munari et al. = -0.1$`\pm `$0.3 km s-1). Taken together, the radial velocities indicate that StH$`\alpha `$190, like other symbiotics, is a binary, but almost certainly a closer binary than the d- and s-types. Munari et al. (2001) suggest a possible period for StH$`\alpha `$190 of 171 days, although such a period does not fit very well the radial velocities observed here. A possible shorter period of $``$37–39 days may provide a better fit to both sets of radial velocities, however, more measured velocities, spanning time scales from a few days to a few weeks, will be needed to determine an accurate period.

A definitive orbital period is crucial in establishing the past, and future, history of StH$`\alpha `$190. Under the assumption that the primary is in synchronous rotation, a short orbital period would be required. In such a case, the projected rotational velocity of $`v\mathrm{sin}i`$=100 km s-1 yields a maximum binary period of 4.2 days (derived from the mass estimate and gravity, which yields a radius of 8.3R). Using the critical rotational velocity as an upper limit (v$``$ 230 km s-1 from the bottom panel of Figure 2), then a lower limit to the binary period of 1.8 days is derived. In such a system, StH$`\alpha `$190 would have undergone common envelope evolution when the current white dwarf was an AGB star. The current evolutionary stage of the G-star suggests that it will soon ascend the red-giant branch, at which point another common envelope phase would occur, resulting in further orbital shrinkage. On the other hand, if the period is much longer ($``$1–6 months) then the system has probably not undergone common envelope evolution and may not in the future. In this case, the rapid rotation of the G-star is probably due to spin-up from the accretion of material from a massive AGB wind (Jeffries & Stevens 1996).

With the few and scattered radial-velocity determinations to date, differences of at least 15 km s-1, and probably as large as 40 km s-1, are indicated. Observing StH$`\alpha `$190 intensively and obtaining high-quality spectra will allow for the determination of the orbital period, and will be crucial in establishing the past and future evolution of StH$`\alpha `$190.

## 4 CONCLUSIONS

The first high-resolution spectral analysis for stellar parameters and abundances in a d’-type symbiotic has revealed it to contain a rapidly rotating G-subgiant, or early giant, with Teff=5300$`\pm `$150K and surface gravity of log g=3.0$`\pm `$0.5. Abundances of \[Fe/H\] and \[Ca/H\] indicate a near-solar metallicity, while barium is enhanced, with \[Ba/Fe\]$``$ +0.5. Thus, StH$`\alpha `$190 is s-process enriched and probably accreted matter from its companion, now a white dwarf, when it was an s-process enriched thermally pulsing AGB star. Measured radial-velocity variations already betray the presence of the companion. If the rapid rotation of the cool component reflects synchronous rotation, then the current orbital period must be between about 1.8 to 4.2 days. If the orbital period is much longer than this, then the rapid rotation of the G-subgiant/giant is due presumably to the addition of angular momentum from a massive AGB wind (Jeffries & Stevens 1996). The current orbital period of StH$`\alpha `$190 is not yet well-constrained, so additional radial velocity measurements of this system, and other d’-symbiotics in general, are warranted.

We thank the staff of McDonald Observatory for technical support. This research is supported in part by the National Science Foundation (AST99-87374).


text:
# Pulse Bifurcations and Instabilities in an Excitable Medium: Computations in Finite Ring Domains

## I Introduction

One-dimensional excitable media exhibit nonlinear traveling waves such as wave trains and solitary pulses. Examples include concentration waves in chemical reactions in solution and on surfaces , signal propagation in neurons and in cardiac tissue . These processes are often modelled by reaction-diffusion (RD) equations of activator-inhibitor type . A pulse is a localized structure; it may result from a finite amplitude perturbation of a linearly stable rest state. The pulse shape usually decays exponentially as a function of the distance from the pulse center. Pulses can be analytically approximated in the asymptotic limit where the dynamics of the activator is much faster than that of the inhibitor variable, and where the activator variable diffuses, but the inhibitor does not. There, it can be shown that a pulse exists and is stable .

More recently, spatiotemporal chaos has been found in a variety of excitable model systems. Several examples have been reported in models whose kinetics possess three distinct homogeneous steady states (fixed points); apart from the stable rest state, these media also exhibit two unstable fixed points . One example arises in a model for catalytic CO oxidation, where the inhibitor is a so-called surface structure variable; its kinetic nullcline displays a sigmoidal shape - threshold dependence on the activator variable (adsorbate coverage) . In contrast to the standard linear dependence of the inhibitor nullcline on the activator, this functional form leads to the additional fixed points mentioned above. Similar behavior is often seen in models in physiology describing the dynamics of the membrane voltage (activator) controlled by so called gating variables for the ion channels (inhibitor). The kinetics of these gating variables display again a threshold-type, sigmoidal dependence on the membrane voltage, often leading to three (spatially homogeneous) fixed points for the PDE . If only a single gating variable is involved, these models are qualitatively similar to the equations for catalytic CO oxidation studied here; a good example is provided by the Morris-Lecar model used to describe the membrane potential in a barnacle muscle fiber .

If the inhibitor kinetics are fast enough, the CO oxidation model displays an instability that has been colloquially named backfiring of pulses . Related phenomena include the wave-induced chemical chaos found in the Gray-Scott model for an autocatalytic chemical reaction , as well as complex behavior in amplitude equations describing dynamics near a Takens-Bogdanov (TB) point . In addition, even reaction-diffusion media whose kinetics are characterized by a single fixed point, sometimes display pulse instabilities and backfiring under excitable conditions . An analysis of the traveling wave ordinary differential equations derived from the original PDEs for the CO oxidation model reveals the basic mechanism for the destruction of stable pulses leading to complex behavior: it involves a so-called T-point  as well as a rich web of bifurcations of pulse solutions . More recently, a similar bifurcation structure has been found in a model for intracellular calcium waves in pancreatic acinar cells  and in the ultrarefractory version of the FitzHugh-Nagumo model . The resulting complex dynamics is often governed by a coherent structure described as a “wave emitting front” , a nonlinear front involving a spatially uniform unstable state that invades a spatially uniform stable one. The unstable state behind the propagating front evolves into spatiotemporal chaos. Similar behavior has been seen in the amplitude equations in the neighborhood of the TB point, and was named “chaotic nucleation” by those authors .

In this paper we investigate instabilities and bifurcations in a model of a catalytic surface reaction . Using a computer-assisted approach we calculate solution branches in large, finite domains, under periodic boundary conditions (which we will, with slight abuse of terminology, refer to as pulses) and perform linear stability analysis of these pulse solutions. The eigenvalues of the corresponding Jacobian reflect the growth exponents of perturbation modes. The emphasis of the present paper is on the form of the spectra and in particular on the shape of the critical modes; this is in contrast to earlier work which focused mostly on the pulse bifurcation structure .

A homogeneous steady state of reaction-diffusion equations in an infinite system possesses a continuous spectrum, reflecting the growth rates of perturbation modes. The eigenmodes are harmonic functions, whose wavenumber parametrizes the corresponding eigenvalues. In a very long or infinite system, a single pulse can conceivably be considered (in the appropriate norm) as a perturbation of this uniform solution. Under these circumstances, the modes far away from the pulse are still delocalized, harmonic wave modes; the continuous spectrum remains essentially unchanged. Furthermore, additional perturbation modes exist that are localized at the site of the pulse and decay exponentially away from it. The corresponding eigenvalues are discrete, in contrast to the continuous band of eigenvalues belonging to the nonlocalized modes. In infinite or periodic systems, the Goldstone mode is an example of such a localized mode. In a spatially homogeneous system it is given by the spatial derivative of the pulse profile, and accounts for a shift in space (due to translational invariance); the corresponding eigenvalue is zero.

A bifurcation or instability of a given solution upon change of a single control parameter is typically accompanied by a single real eigenvalue or a pair of complex conjugate eigenvalues crossing the imaginary axis. For a real eigenvalue, a saddle-node bifurcation is the most common case, while alternative bifurcations like the transcritical or pitchfork bifurcations require certain symmetries for the equations and the solution. The case of a drift pitchfork bifurcation of a stationary pulse was studied in . A destabilizing perturbation grows in an oscillatory manner if the eigenvalues are complex and in a non-oscillatory manner if the corresponding eigenvalue is real. The first case leads to a Hopf bifurcation, while the latter case may correspond either to a saddle-node or to a transcritical or pitchfork bifurcation.

Three examples of dynamics where localized initial conditions beyond the uniform solution excitation threshold do not eventually evolve into stable pulses, are shown in Fig. 1. For high excitability resp. small excitation threshold, pulse-like initial conditions evolve in a traveling pulse that splits off new pulse-like excitations traveling in the opposite direction. This is the phenomenon termed backfiring in . It can eventually evolve in spatiotemporally chaotic (I) or complex periodic (II) fashion. For large excitation threshold, modulated traveling pulses are also seen (III). Their shape undergoes a breathing, periodic variation.

## II Model and Methods

We investigate a model of activator-inhibitor type originally derived for CO oxidation on Pt(110) , in a parameter regime where the kinetics give rise to one stable and two unstable steady homogeneous solutions. This model is related to the FitzHugh-Nagumo system and describes the interaction of a fast activator $`u`$ and a slow inhibitor variable $`v`$:

$`_tu`$ $`=`$ $`_x^2u+{\displaystyle \frac{1}{ϵ}}u(1u)(u{\displaystyle \frac{b+v}{a}}),`$ (1)
$`_tv`$ $`=`$ $`f(u)v,`$ (2)
$`f(u)`$ $`=`$ $`\{\begin{array}{cc}0,\hfill & 0u<1/3\hfill \\ 16.75u(u1)^2,\hfill & 1/3u1,\hfill \\ 1,\hfill & 1<u\hfill \end{array}`$ (6)
with $`x[0,L]`$ and periodic boundary conditions.

The time scales ratio $`ϵ>0`$ is used as the control parameter. The case $`ϵ0`$ corresponds to the aforementioned asymptotic limit where stable pulses are expected. The parameter $`b`$ controls the excitability threshold of the system: the value of $`b`$ is proportional to the magnitude of the critical perturbation which will trigger a pulse. The value of $`a`$ is fixed at $`0.84`$ throughout the paper. We will vary $`b`$, thus varying the excitability. In the parameter range considered here, three relevant fixed points exist: the stable state $`A=(0,0)`$, the saddle $`B=(b/a,0)`$, and unstable focus $`C`$.

We shall now examine the stability of pulses traveling in a background of the stable rest state $`A`$. Since we investigate solutions moving with fixed velocity $`c`$, the analysis of their stability is performed in the comoving frame $`z:=xct`$:

$`_tu`$ $`=`$ $`_z^2u+c_zu+{\displaystyle \frac{1}{ϵ}}u(1u)(u{\displaystyle \frac{b+v}{a}}),`$ (7)
$`_tv`$ $`=`$ $`c_zv+f(u)v.`$ (8)
In this frame, traveling waves with speed $`c`$ correspond to time independent, steady solutions. Linearization of these equations around a stationary solution $`u_0(z),v_0(z)`$ yields an eigenvalue problem for small perturbations $`(r(z,t),s(z,t))(r(z),s(z))e^{\lambda t}`$:

$`\left(\begin{array}{c}r(z)\hfill \\ s(z)\hfill \end{array}\right)=\lambda \left(\begin{array}{c}r(z)\hfill \\ s(z)\hfill \end{array}\right),`$ (13)
(14)
$`=\left(\begin{array}{cc}_z^2+c_z+g_1(z)& g_2(z)\\ _uf(u_0)& c_z1\end{array}\right)`$ (17)
with

$`g_1(z)`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}\left(u_0(u_01)+(u_0{\displaystyle \frac{b+v_0}{a}})(2u_01)\right),`$ (18)
$`g_2(z)`$ $`=`$ $`{\displaystyle \frac{u_0}{ϵa}}(u_01).`$ (19)
The linear stability problem amounts to the determination of the spectrum of the Jacobian $``$ in (17). For the homogeneous steady states $`A`$ or $`B`$, at least one of the off-diagonal matrix elements is zero. Thus, the diagonal elements of the matrix $``$ suffice to compute the spectrum. For the $`A`$ steady state it is:

$$\lambda _{A,1}=\frac{b}{ϵa}k^2+ick,\lambda _{A,2}=1+ick,$$
(20)
and for the $`B`$ steady state it is:

$$\lambda _{B,1}=\frac{b}{ϵa}\left(\frac{b}{a}1\right)k^2+ick,\lambda _{B,2}=1+ick,$$
(21)
where $`k`$ is the wavenumber of the perturbation. In the case of periodic boundary conditions studied here, $`k=n\frac{2\pi }{L}`$ applies. The real part of the eigenvalues $`\lambda _2`$ is $`1`$. The eigenvectors are then $`(r,s)^T=(1,0)^Te^{ikz}`$ and $`(0,1)^Te^{ikz}`$.

In general, traveling wave solutions $`u_0(z),v_0(z)`$ and the eigenfunctions of the Jacobian $``$ are not available in closed form. Thus, the problem has to be approximated numerically. We approximate pulses by computing steady solutions in a finite length system and the traveling frame through a pseudospectral discretization of Eqs. (7) with periodic boundary conditions and Newton-Raphson iterations. The velocity $`c`$, which is not known a priori, is formally an additional variable along with the Fourier coefficients of the solution. One additional pinning condition singles out one of the infinitely many solutions existing due to translational invariance and allows the numerical computation of the speed. The eigenfunctions and the spectrum of the Jacobian are obtained in Fourier space resulting from a 200 (resp. 250 for case I) mode decomposition of the stationary solution. The zero eigenvalue, which always exists due to the translational symmetry of the problem, is used as a numerical accuracy check and has been obtained with a precision of $`10^3`$ or better. The time evolution of unstable solutions was computed using an explicit finite difference scheme to solve Eqs. (1), discretizing space in 1024 points and using a time step of $`\mathrm{\Delta }t=0.0122`$.

Because traveling solutions fulfill the conditions $`_tu=_tv=0`$ in the comoving frame, they can also be obtained in the traveling wave ordinary differential equations (TWODE) following from Eqs. (7):

$`{\displaystyle \frac{du}{dz}}`$ $`=`$ $`w,`$ (22)
$`{\displaystyle \frac{dw}{dz}}`$ $`=`$ $`cw+{\displaystyle \frac{1}{ϵ}}u(u1)(uu_{th}),`$ (23)
$`{\displaystyle \frac{dv}{dz}}`$ $`=`$ $`(vf(u))/c,`$ (24)
with $`u_{th}`$ = $`(b+v)/a`$. In this frame, a homogeneous solution corresponds to a fixed point, a pulse to a homoclinic orbit, and a front to a heteroclinic orbit. Consequently, in the parameter range studied here, three relevant fixed points exist: $`𝒜=(0,0,0)`$, $`=(b/a,0,0)`$, and the focus $`C`$.

## III Results

We consider three cases, at increasing values of the excitation threshold, controlled through the parameter $`b`$ at fixed $`a=0.84`$: Case I ($`b`$ = 0.07), case II ($`b`$= 0.15) and case III ($`b`$ = 0.20). We proceed as follows: first, branches of pulse solutions on a ring, representative of ,,true” pulse solutions in a infinite domain, are presented for the relevant range of the control parameter $`ϵ`$; we characterize these pulses by their speed $`c`$ (Fig. 2, 3b). We then show pulse profiles and spectra at selected values of $`ϵ`$ shortly before and after the onset of instability (Figs. 3, 7, 8), and present the destabilizing mode $`(r,s)^T`$ (Figs. 6, 7, 8). Representative post-instability spatiotemporal dynamics can be found in Fig. 1.

Fig. 2 shows pulse speed as a function of the parameter $`ϵ`$ for the three cases. The thick lines correspond to stable pulses with background state $`A`$, while dashed lines correspond to unstable pulses with background state $`B`$. The transition point between the two families is the so called T-point , denoted by a T. This point is defined as a double heteroclinic connection between the fixed points $`A`$ and $`B`$ in the frame of the traveling wave ODEs, Eqs. (22). In the frame of the original equation (1), these heteroclinic orbits correspond to fronts. The branch of pulses corresponding to homoclinic orbits to $`𝒜`$ in the TWODE may (Fig. 2.II, III) or may not (Fig. 2.I) spiral into the T-point. Spiraling is observed when two of the eigenvalues of the linearization of the TWODE around the fixed point $``$ are complex conjugate .

### A Case I

Fig. 2.I shows speed as a function of the parameter $`ϵ`$ for $`b=0.07`$ computed for pulses in a large system with periodic boundaries. The picture appears at first glance identical with the result found in a continuation of homoclinic orbits in the TWODE . The branch of pulses to $`A`$ (thick line) does not spiral into the T-point. To be exact, the T-point only exists for pulses on the infinite line; what – at the resolution of our picture – still appears as a T-point will be discussed in more detail below. To understand the nature of the instability in this case, we focus on the solutions near the T-point value $`ϵ_T0.10744`$.

As mentioned above, for an infinite system, the T-point corresponds to a double heteroclinic connection in the traveling wave ODEs, Eqs. (22). Close to $`ϵ_T`$ the orbits homoclinic to $`𝒜`$ also approach the fixed point $``$ (and vice versa); this approach to $``$ makes the dynamics in $`z`$ progressively slower. Thus, close to $`ϵ_T`$, pulses with rest state $`A`$ (resp. $`B`$) will also exhibit extensive regions close to $`B`$ (resp. $`A`$), practically giving rise to fronts between $`A`$ and $`B`$ within the pulse profile. When we perform a continuation for a pulse in a ring of finite length, the space on the ring is apart from the excitation plateau divided between residence close to $`A`$ and residence close to $`B`$, so that the total period is constant as we vary $`ϵ`$.

To study this behavior, numerical continuation techniques are needed. We computed the stationary solutions of Eqs. (7) on a ring of length $`L=200`$ with 250 modes, 1024 collocation points and a parameter step size $`\mathrm{\Delta }ϵ=𝒪(10^{10})`$. Periodic approximations of pulse solutions with rest state $`A`$ are shown in Fig. 3a as they approach the T-point and undergo the gradual transition to pulse solutions with rest state $`B`$. One can clearly recognize the increasing domain $`B=(b/a,0)`$ at the back of the pulse. As the state $`B`$ is unstable, one might expect that the pulse in the ring loses stability as soon as the $`B`$ domain gets large enough. However, one should be aware of the fact that the $`B`$ domain is moving with the speed of the pulse. It is possible that the $`B`$ domain is only convectively unstable in the comoving frame. In other words, the perturbations growing on the $`B`$ state in a stationary frame may spread slower than the pulse motion, and the pulse in an infinite domain will be stable. A mathematically precise description of this phenomenon has been given by Sandstede and Scheel . A related result has been obtained by Nii , who shows that eigenvalues accumulate in the area bounded by the essential spectra of $`A`$ and $`B`$. The opposite case, i.e. perturbations on the $`B`$ plateau spread faster than the pulse speed, is described in the next subsection (case II).

Typical spectra are shown in Fig. 4 for solutions along the $`cϵ`$-branch. To facilitate a better comparison, the continuous spectra of the rest states $`A`$ (Eqs. (20)) and $`B`$ (Eqs. (21)) are depicted as solid lines; the computed eigenvalues are denoted by circles. Fig. 4a shows the spectrum of a solution approximating a stable pulse with rest state $`A`$, while Fig. 4f shows that of a solution approximating an unstable pulse with rest state $`B`$. In both cases, the parameter $`ϵ`$ is “far” from the T-point bifurcation and the eigenvalues belonging to nonlocalized eigenvectors compare well with the continuous spectrum of the respective rest state. Figs. 4b-e show the gradual transition between the two cases. Fig. 4b shows the spectrum for a stable solution for which the domain $`B`$ is just starting to appear in the back of the pulse. Its width is $`L_B5`$. Several discrete eigenvalues have appeared on the negative real axis; a similar spectrum has been predicted and found by Sandstede and Scheel .

The details of the spectrum tranformation are as follows: Initially (far to the left of the T-point in Fig. 2a), the pulse spectrum contains two discrete complex conjugate eigenvalues (which remain more or less unchanged through the bifurcation diagram) as well as two “pieces” of continuous spectrum: a parabola, and a vertical line. For our calculations in the finite domain, initially the parabola does have a “tip” (an eigenvalue on the real axis) while the vertical line does not have an eigenvalue on the real axis. As we vary $`ϵ`$ the parabola approaches the straight line; eventually the parabola appears as crossing the line. While this occurs, a number of complex conjugate pairs of eigenvalues emerge from the parabola and line approximating the essential spectrum. Fig. 4a shows the situation where the first of these pairs has appeared. Upon further increase of $`ϵ`$, these pairs collapse on the negative real axis where they split and become real. These real eigenvalues lie in the absolute spectrum of state $`B`$ (see ); they appear when the solution has a visible $`B`$ “shoulder” behind the pulse, see Fig. 4b. As $`ϵ`$ varies these real eigenvalues become complex conjugate. Several complex eigenvalue pairs of the parabola give the appearance of a secondary parabola - this is the first intimation of what will become the essential spectrum of $`B`$. Only two of the “momentarily real” eigenvalues remain real \- one of them forms the tip of the “new parabola”, while the other forms the tip of the “old parabola”.

Gradually, as the $`B`$ plateau grows, two distinct phases of spectrum movement are observed. First this new $`B`$-associated parabola moves to the right in the complex plane and at some point it starts crossing the imaginary axis. That is precisely the first saddle-node bifurcation we observe - the critical eigenvalue is the tip of this secondary parabola, whose origin we just discussed, see Fig. 4c. Continuing further on the middle, already unstable branch, the $`B`$-domain gets wider and many secondary (apparently Hopf) bifurcations occur as the eigenvalues of this “discretized” parabola cross successively the imaginary axis. When both $`A`$\- and $`B`$-plateaus are equally present, one expects to see echoes of both $`A`$\- and $`B`$-continuous spectra in the solution spectrum, and that is indeed seen in Fig. 4d. Near the lower saddle-node bifurcation, the $`B`$-plateau is almost as wide as the system length. There, one may describe the solution more reasonably as an approximation to a pulse with rest state $`B`$ with a short $`A`$ domain at its front. The eigenvalue spectrum at this situation is depicted in Fig. 4e. The second (lower) turning point occurs when the “discretized” parabola of eigenvalues that corresponded initially to continuous spectrum of $`A`$ moves to the right and the real eigenvalue at the tip of this discretized parabola crosses the imaginary axis.

Eventually, the old $`A`$-parabola eigenvalues merge with the $`B`$-parabola as the $`A`$-shoulder gradually disappears. Fig. 4f shows the spectrum for a solution approximating a pulse with rest state $`B`$ far from the T-point bifurcation for a high values of $`ϵ`$. In the end, the two “tips” (the old and the new one) collide and become again a complex conjugate pair - the final pure-$`B`$ parabola has no real tip.

Sandstede and Scheel in  have proven that the spectrum of pulses with rest state $`A`$ with a “shoulder” of $`B`$ state as in Fig. 3a is comprised of the essential spectrum of the homogeneous background state $`A`$ plus eigenvalues lying in the absolute spectrum  of $`B`$ located on the real axis and bounded by a maximum value which is negative. We deviate from their picture, when the $`B`$ shoulder becomes relatively large. In other words, their result on the infinite line is reproduced by our numerical stability analysis in the ring, if the condition $`L_B<<L`$ holds. The spectrum transformations, however, reflect the influence of the periodic boundary conditions and the finite length. Because of the violation of the above assumption near the T-point (see Fig. 3a) the spectra we find for a large box are not well described as the union of the essential spectrum of $`A`$ and the absolute spectrum of $`B`$. The results of  also imply that if one considers pulses on the infinite line instead of the large wavelength approximation with periodic boundary conditions used here, one will find a monotonic approach of the pulse branches to the T-point. In their scenario no eigenvalues are crossing the real axis upon approaching the T-point (in contrast to the computations on the ring) and no bifurcation (like the saddle-node bifurcation found for the ring) can occur.

Focusing on the first instability ($`SN_1`$), the critical eigenmode corresponding to the saddle-node appears localized at the back of the pulse (see Fig. 5). In addition, it is acting mainly on the activator (Fig. 5b); we observe that new pulses can now split off the existing pulse (“backfiring”). Fig. 6 shows the time evolution beyond the saddle-node bifurcation in a comoving frame. The initial condition is the stable pulse just before the bifurcation. Fig. 6a shows how a perturbation with support over the $`B`$-plateau grows with time, while the shape of the preceding front essentially does not change. The back (what starts as the $`B`$-plateau) grows until its maximum reaches the excited state $`u=1`$. The excited domain widens, but its plateau is unstable against oscillating perturbations . Therefore, the situation rapidly evolves, and the growing excited domain breaks down, giving rise to two new seeds for $`A`$-pulse like entities. This behavior has been described by the term “backfiring” . Backfiring occurs repeatedly, and the newly generated pulse-like structures annihilate upon collision with similar objects traveling in the opposite direction; this interplay of instability, new pulse generation and annihilation gives rise to the non-periodic space-time behavior shown in Fig. 1.I.

Altogether, our results give an interesting and well resolved picture of stability of pulses on a large finite ring near the T-point. In contrast, earlier computations  have used coarser steps in the continuation algorithms and smaller domains. The present computations show how the pulse solutions change in a gradual fashion; and that two distinct solutions branches show an extremely narrow hysteresis; the middle branch mediates the transition between the two different pulse types. The main changes appear in an extremely tiny region of parameter space and can therefore only be captured with careful numerics. The spectrum transformations near parameters that exhibit a T-point in the corresponding TWODE appear in a similar fashion for case II and III studied below. There, however, the T-point is not involved in the disappearance of stable traveling pulses through instability and bifurcation. Thus, we will not discuss it in these cases.

### B Case II

Fig. 2.II shows the $`cϵ`$ diagram for $`b=0.15`$. As in case I, we again observe here the transformation from pulses with rest state $`A`$ to pulses with rest state $`B`$ at what appears like a T-point. The branch corresponding to “pulses to $`A`$” spirals into this T-point. This behavior is caused by imaginary eigenvalues of the fixed point $`B`$ in the TWODE at the T-point conditions, and has been predicted from general arguments . Thus, the branch of initially stable pulses with rest state $`A`$ undergoes a sequence of saddle-node bifurcations upon approaching the T-point. This sequence of saddle-node bifurcation is consistent with the results of Sandstede and Scheel . The spectrum of $`A`$-pulses near the T-point in the infinite system approaches the union of the essential spectrum of $`A`$ and the absolute spectrum of $`B`$. In this case, however, the absolute spectrum of $`B`$ contains real positive eigenvalues. One may, therefore expect infinitely many saddle-node bifurcations as the T-point is approached. In the stability-relevant first saddle-node bifurcation, the solution branch of stable wavetrains with rest state $`A`$ turns around and becomes unstable. As is required for a saddle-node bifurcation, a single real eigenvalue crosses the imaginary axis. This can be seen in Fig. 7, left column. The two spectra correspond to pulse-like solutions before and after the saddle-node bifurcation along the solution branch. The unstable pulse and the destabilizing eigenmode are shown in Fig. 7, right column. Once more, the destabilizing mode affects primarily the back of the pulse. It is worth noting, that upon further continuation the branch of pulses to $`A`$ spirals towards the T-point related situation over a cascade of saddle-node bifurcations. Each saddle-node bifurcations adds an additional positive eigenvalue to the spectrum. This is in contrast to the study in , where a single eigenvalue repeatedly crosses zero along a spiraling pulse branch with a cascade of saddle-node bifurcation.

Numerical simulations of the model Eqs. 1 shortly after the saddle-node bifurcation, for values of $`ϵ`$ for which no solution with rest state $`A`$ exists, exhibits the phenomenon we termed above “backfiring”. In the transients, our pulse-like object generates near its back other pulse-like entities traveling in the opposite direction, see Fig. 1.II. After this transient period, the resulting spatiotemporal pattern in our finite domain becomes periodic in time. Simulations show, though, that this observation depends on the initial conditions: non-periodic patterns like the one shown in Fig. 1.I may also appear for the same parameter values.

### C Case III

Fig. 2.III shows the $`cϵ`$ 0 for $`b=0.2`$. Here, the pulse solution with rest state $`A`$ becomes unstable through a Hopf bifurcation even before the first saddle-node bifurcation is reached. The eigenvalue spectrum is shown in Fig. 8, left column, on both sides of this Hopf bifurcation. It can be seen that one discrete pair of complex conjugate eigenvalues crosses the imaginary axis. The second column of Fig. 8 shows the pulse-like solution after the bifurcation as well as the real and imaginary parts of the critical eigenmode $`(r,s)^T`$. Note that the perturbation due to the eigenmode is localized and has its main contribution – once more – at the back of the pulse. As the corresponding eigenvalues are complex conjugate, the resulting perturbation oscillates with time at a frequency given by the imaginary part of the eigenvalues.

Numerical integration of the model Eqs. (1) shows that the Hopf bifurcation which leads to destabilization is supercritical . The resulting pattern after destabilization consists of a modulated traveling pulse which oscillates in time, especially at its back (compare the simulation shown in Fig. 1.III.)

## IV Conclusion

We have investigated the transition from stable pulse propagation to various regimes of more complicated spatiotemporal dynamics, namely modulated pulses, periodic and chaotic pulse backfiring. In all three cases, the transition can be explained by either a Hopf instability (modulated pulses) or a saddle-node bifurcation (leading to backfiring) of the stable pulse solution. In all cases, the transition is connected with either a single or a pair of complex conjugate discrete eigenvalue(s) with zero real part(s). In a finite domain with periodic boundary conditions – the typical experimental setup for investigation of pulses – spectra change continuously near the T-point in the fashion described in case I. The form of the corresponding critical eigenmode(s) allows some insight into how pulses become unstable resp. evolve in space and time. The dynamics in general still contains mostly propagating localized pulse-like structures whose evolution is governed by the unstable eigenmode(s) in the Hopf case (modulated) or by the critical eigenmode of the saddle-node bifurcation. Typically, the critical eigenmodes have support at the back of the pulse.

The results here should carry over to models with similar phenomenology mentioned in the introduction. Preliminary results  show that the transition to wave-induced chemical chaos in the Gray-Scott model  is also caused by a saddle-node bifurcation of pulses near a T-point. T-points can only be found in systems with multiple homogeneous fixed points (e.g. one stable rest state and two additional unstable steady states). The complex behavior seen in media with a single stable fixed point  may be caused through a different mechanism. For a model of the catalytic NO-CO reaction, upon change of control parameter first modulated traveling waves are seen and then periodic backfiring is found. This cannot be due to a T-point, but may instead be caused by a global bifurcation of the periodic modulated pulses - a scenario already suggested in a study of the present model with different control parameters . Finally, it is important to note that a simple instability of a finite wavelength pulse solution, like the Hopf bifurcation in case III, only leads to a modulation of the shape, while a saddle-node bifurcation limits the existence of a certain type of pulses and may give rise to replication of pulse-like structures (cases I, II). The role of saddle-node bifurcations in the replication of pulses and creation of space-time defects has recently been investigated in the Gray-Scott model  and in the Complex-Ginzburg-Landau equation .

Acknowledgement: We thank Björn Sandstede and Arnd Scheel for enlightening discussion regarding Ref. . IGK acknowledges the support of the NSF (USA) and of the A.v.Humboldt Stiftung (Germany).