text:
# Bifurcations in One Degree-of-Vibration Quantum Billiards

## Abstract

We classify the local bifurcations of one dov quantum billiards, showing that only saddle-center bifurcations can occur. We analyze the resulting planar system when there is no coupling in the superposition state. In so doing, we also consider the global bifurcation structure. Using a double-well potential as a representative example, we demonstrate how to locate bifurcations in parameter space. We also discuss how to approximate the cuspidal loop using AUTO as well as how to cross it via continuation by detuning the dynamical system. Moreover, we show that when there is coupling, the resulting five-dimensional system—though chaotic—has a similar underlying structure. We verify numerically that both homoclinic orbits and cusps occur and provide an outline of an analytical argument for the existence of such homoclinic orbits. Small perturbations of the system reveal homoclinic tangles that typify chaotic behavior.

### MSC NOS 37N20, 37K55, 37M20

## 1 Introduction

Quantum chaos is an underdeveloped area of dynamical systems theory. One purpose of studying it is to generalize the notions of classical Hamiltonian chaos to the quantum regime, which has not yet been done in a universally accepted manner. One type of quantum chaos is often called semiquantum chaos, as these systems consist of classical (Hamiltonian) variables coupled with quantum variables. Vibrating quantum billiards are a representative example of semiquantum chaotic systems. They may be used as models for quantum well microdevice components (such as quantum dots and quantum wires), Fermi accelerators, and intranuclear particle behavior.

In the present paper, we consider the bifurcation structure of vibrating billiard systems. Quantum billiards describe the motion of a point particle undergoing perfectly elastic collisions in a bounded domain. The particle’s motion is described by the Schrödinger equation with Dirichlet boundary conditions. One defines the “degree-of-vibration” (dov) of a billiard as the number of boundary dimensions that vary with time. If the boundary is time-independent, the billiard is said to have zero dov. The linear vibrating billiard and the radially vibrating spherical billiard have a single dov, and the rectanglular billiard with time-varying length and width has two dov.

A zero dov quantum billiard exhibits only integrable behavior if it is globally separable. It must, for example, be convex and describable by a non-composite geometry. If part of the billiard is concave (or is composite, like the stadium billiard), it may behave chaotically, as it shares many of the instability properties of the Anosov diffeomorphism. Blümel and Esser found quantum chaos in the one-dimensional vibrating quantum billiard. Liboff and Porter extended these results to spherical quantum billiards with vibrating surfaces and derived necessary conditions for chaotic behavior. They also generalized their results to other one dov billiards. The purpose of the present paper is to examine bifurcations in single dov quantum billiards that occur when one alters the potential in which the billiard resides.

Vibrating quantum billiards are important for several reasons. Though an idealized model, they are nevertheless useful for the study of quantum chaos. From a more practical standpoint, vibrating quantum billiards may be applied to several problems in physics. The radially vibrating spherical quantum billiard, for example, may be used as a model for particle behavior in the nucleus as well as for the quantum dot microdevice component.. Additionally, the vibrating cylindrical billiard may be used as a model for the quantum wire, another microdevice component. Other geometries of vibrating quantum billiards have similar applications. They may also be used as models of Fermi acceleration of cosmic rays. The study of quantum chaos in vibrating billiard systems is thus important both because it expands the mathematical theory of dynamical systems and because it can be applied to problems in nuclear and mesoscopic physics.

In the present paper, we show that saddle-centers are the only type of bifurcations that can occur in one dov quantum billiards. When there is no coupling in the superposition state, we show how to analyze the resulting planar system analytically and numerically. Considering a double-well potential as a representative example, we demonstrate how to locate bifurcations in parameter space. We also discuss how to approximate the cuspidal loop using AUTO as well as how to continue past it by detuning the dynamical system. We also mention a shooting method for a more detailed analysis of the cuspidal loop. Moreover, we show that when there is coupling, the resulting five-dimensional system—though chaotic—has a similar underlying structure. We verify numerically that both homoclinic orbits and cusps occur, and we outline an analytic argument that demonstrates the existence of homoclinic orbits. Small perturbations of the system reveal homoclinic tangles that typify chaotic behavior.

## 2 Equations of Motion

The goal of the present project is to examine the behavior of one dov quantum billiards in various potentials in order to determine the effects of the potential on the behavior of the system. In particular, we analyze bifurcations of equilibria both analytically and numerically. We consider a two-state superposition solution to the vibrating billiard, and we examine the above problem for both the case in which the two states experience coupling and that in which they do not.

The present problem is described by the Schrödinger equation with solutions that are constrained to vanish on a time-dependent boundary $`a(t)`$. That is,

$$\psi (r,t;a(t))=0\text{ for }r=a(t).$$
(1)
Because of the time-dependence of the boundary, the above boundary condition is nonlinear. The (mathematical) problem at hand is to find a boundary $`a(t)`$ such that Dirichlet boundary conditions are satisfied on it. One can then, in principle, insert $`a(t)`$ into the eigenfunctions in our normal mode expansion of the wave $`\psi (r,t;a(t))`$ in order to obtain nonlinear normal modes. When the radius $`a(t)`$ behaves chaotically, the nonlinear normal modes are examples of quantum-mechanical wave chaos. One derives coupled nonlinear ordinary differential equations for $`a(t)`$ (and other variables) using a Galërkin approximation. Considering a two-term superposition state then corresponds to taking a two-term Galërkin projection. The equations of motion that one obtains depend on whether a particular subset of the quantum numbers in the two states are the same. For the case of the radially vibrating sphere, the quantum numbers in question are the orbital and azimuthal quantum numbers $`l`$ and $`m`$, respectively.

If these quantum numbers are the same in the two states, there is a coupling between them. The evolution of the system is then described by

$$\dot{x}=\frac{\omega _0y}{a^2}\frac{2\mu Pz}{Ma},$$
(2)
$$\dot{y}=\frac{\omega _0x}{a^2},$$
(3)
$$\dot{z}=\frac{2\mu Px}{Ma},$$
(4)
$$\dot{a}=\frac{P}{M},$$
(5)
and

$$\dot{P}=\frac{V}{a}+\frac{2[ϵ_++ϵ_{}(z\mu x)]}{a^3},$$
(6)
where $`x`$, $`y`$, and $`z`$ are Bloch variables

$$x=\rho _{12}+\rho _{21},y=i(\rho _{21}\rho _{12}),z=\rho _{22}\rho _{11},$$

$$x^2+y^2+z^2=1,$$
(7)
$`\rho _{mn}A_mA_n^{}`$ is the density matrix, $`a`$ is a displacement, $`P`$ is its conjugate momentum, $`M`$ is the mass of the billiard, $`mM`$ is the mass of the confined particle, $`\mu >0`$ is the coupling coefficient between the two eigenstates, $`V=V(a)`$ is the potential of the billiard boundary, $`\omega _0\frac{ϵ_2ϵ_1}{\mathrm{}}`$, $`ϵ_\pm \frac{(ϵ_2\pm ϵ_1)}{2}`$, and $`ϵ_1`$ and $`ϵ_2`$ $`(ϵ_2ϵ_1)`$ are the energies of the two states. It has been shown that these equations exhibit quantum chaotic behavior.

If there is no coupling between the two eigenstates, the evolution of the system is described by a one degree-of-freedom Hamiltonian. The equations of motion are

$$\dot{a}=\frac{P}{M},$$

$$\dot{P}=\frac{V}{a}+\frac{\lambda }{a^3},$$
(8)
where

$$\lambda 2\left(ϵ_1|C_1|^2+ϵ_1|C_2|^2\right),$$
(9)
$`C_1`$ and $`C_2`$ are constants such that $`|C_1|^2+|C_2|^2=1`$. The energy parameter $`\lambda `$ is necessarily positive because $`ϵ_i>0`$ and the $`|C_i|^2`$ correspond to probabilities.

## 3 Integrable Case: Absence of Coupling

For the planar case, all $`(a,P)`$ that satisfy $`\dot{a}=\dot{P}=0`$ are equilibrium points. Each one is of the form $`(\overline{a},0)`$, where $`\overline{a}`$ satisfies

$$\frac{V}{a}(\overline{a},0)=\frac{\lambda }{\overline{a}^3}.$$
(10)
The eigenvalues of the integrable system (8, 9) at the stationary point $`(\overline{a},0)`$ are

$$\sigma =\pm \sqrt{\frac{1}{M}\left(\frac{^2V}{a^2}(\overline{a},0)+\frac{3\lambda }{\overline{a}^4}\right)}.$$
(11)
These eigenvalues are either real with opposite sign or are pure imaginary pairs, so in the linear analysis, each equilibrium is either a center or a saddle point. If

$$A\frac{^2V}{a^2}(\overline{a},0)+\frac{3\lambda }{\overline{a}^4}>0,$$
(12)
then every equilibrium point is a linear center. This holds, in particular, when the potential $`V(a)`$ has a single minimum (single-well potentials). Previous studies have focused on the harmonic potential

$$V(a)=\frac{𝒱_2}{a_0^2}(aa_0)^2V_2(aa_0)^2.$$
(13)
Another interesting single-well potential is the quartic potential

$$V(a)=\frac{𝒱_4}{a_0^4}(aa_0)^4V_4(aa_0)^4.$$
(14)
In the above equations, the $`𝒱_i`$ are dimensionless parameters.

It is insightful to examine the evolution (particularly in the chaotic case) of the quantum billiard system with the above two potentials and derive a mechanical anology in terms of spring stiffness. Given the same initial conditions and the quartic and quadratic potentials above (and assuming $`V_2=V_4`$ for ease for comparison), one observes that the phase-plane trajectory described by the evolution of $`a`$ and $`P`$ for the quartic potential (14) has a larger radius of curvature (that is, a smaller curvature). For all initial conditions, the trajectory in the quadratic potential has a larger maximum $`a`$. For initial conditions with sufficiently small $`a(0)`$, the quadratic potential induces trajectories with smaller maximum $`|P|`$, but the quartic potential eventually gives a larger maximum $`|P|`$ as $`a(0)`$ is increased.

An equilibrium for which $`A<0`$ is a saddle point. Since the present system is a Hamiltonian system, it is invariant under reflection about the $`a`$-axis, so that the eigenvectors representing the local stable and unstable manifolds are mapped to each other under this reflection. Since the only other possible types of equilibria are centers, it follows that if there is at least one saddle point, the system must have saddle connections. If there is exactly one, the connection is a homoclinic orbit, and there must be two of them emanating from the saddle because the system must have a center (in the right-half plane) on each side of the saddle.

As one considers increasingly excited states of the system (corresponding to larger quantum numbers), the energy parameter $`\lambda `$ is increased. Each saddle will eventually become a center as this occurs. The quantity $`A`$ vanishes at such an equilibrium point. The stationary point then has a double zero eigenvalue with the Jacobian

$$\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)$$
(15)
so that the saddle-center bifurcation that occurs has codimension two and gives rise to a global bifurcation corresponding to the breaking of the separatrix. To find the conditions satisfied at this point one can either solve the simultaneous equations $`\dot{a}=0,\sigma =0`$ or find the points at which the Hamiltonian has a double zero (which is equivalent to solving the system of equations $`H(\overline{a},0)=0,\frac{H}{a}(\overline{a},0)=0`$). If the potential has a constant term $`V_0`$, it does not change the evolution of the system since it does not appear in the equations of motion. We thus let $`V_0=0`$ without loss of generality. One finds that a saddle-center bifurcation occurs when

$$\overline{\lambda }=\overline{a}^3\frac{V}{a}(\overline{a},0)$$
(16)
at the point $`(\overline{a},0)`$ satisfying

$$\frac{V}{a}(\overline{a},0)=\frac{\overline{a}}{3}\frac{^2V}{a^2}(\overline{a},0).$$
(17)
Any solution to (17) giving $`\lambda <0`$ is discarded as nonsensical.

At the saddle-center bifurcation point, the stable and unstable eigenvectors of the equilibrium coincide along the $`a`$-axis, so that the stable and unstable manifolds overlap near the stationary point. This cusp causes difficulties in numerical continuation attempts, as standard continuation techniques fail for this bifurcation. One observes that the two homoclinic orbits that exist when $`A<0`$ have coalesced into one. (As $`A`$ increases, the homoclinic orbit on the left shrinks, becoming a single point at the saddle-center. The orbit has infinite derivative with respect to arclength at the saddle-center point.)

As a specific example of this phenomenon, we consider the quartic potential

$$V(a)=V_4(aa_0)^4+V_3(aa_0)^3+V_2(aa_0)^2+V_1(aa_0),$$
(18)
where $`V_i\frac{𝒱_i}{a_0^i}`$. This potential has either one well or two. In the latter case, there is a single saddle-center bifurcation point. One can find $`\overline{a}`$ and $`\overline{\lambda }`$ exactly in this case, since the former is a root for a degree-three polynomial. We initially restrict ourselves to the case in which $`V_3=V_1=0`$, since all the dynamics of interest remain in this simpler case. Note that the potential is symmetric about $`a_0`$. For numerical purposes, consider the special case corresponding to the parameter values $`a_0=1,V_2=1,`$ and $`V_4=1`$. There is a saddle-center bifurcation at $`\overline{\lambda }=\frac{1}{972}\left[3+\sqrt{3}\right]^3\sqrt{3}0.1888176.`$ The corresponding stationary point is $`(\frac{1}{2}+\frac{\sqrt{3}}{6},0)(0.7886751,0)`$. Using DsTool, we plotted an approximation of the homoclinic orbit emanating from this equilibrium (Fig. 1).

More generally, one can consider even polynomials of higher degree in order to examine vibrating quantum billiards in an $`N`$-well potential. If the polynomial is of degree six or higher, one may not be able to solve for $`\overline{a}`$ exactly in terms of radicals by Galois theory, since equation (17) is polynomial of degree at least five. If its degree is exactly five, one can solve for $`\overline{a}`$ exactly in terms of elliptic functions.

Consider the problem of starting at $`\lambda <\overline{\lambda }`$ and attempting to continue along the bifurcation curve past the saddle-center. Using AUTO and the homotopy method encoded in HomCont, we followed the two homoclinic orbits for $`\lambda =0.15`$ (Fig. 2). The saddle connection in the present system has a codimension greater than one, as both regularity and non-degeneracy conditions are both violated. The present situation is degenerate because for all $`\lambda <\overline{\lambda }`$, there are two homoclinic orbits emanating from the saddle point. Regularity is violated because the saddle point’s two eigenvalues are negatives of each other. (Moreover, their eigenvectors are related by reflection across the $`a`$-axis, since the system is Hamiltonian.)

Because the present system is degenerate and irregular, one cannot continue (in $`\lambda `$) past the saddle-center directly, as described in the AUTO manual. Hamiltonian systems possess a continua of homoclinic orbits, and there are numerical schemes that allow one to handle this phenomenon. One can exploit the reversibility of the system by computing only half of a given saddle connection. However, the cusp at the saddle-center point prevents this from working for the present system. AUTO is incapable of continuing past a cusp, because $`\frac{P}{a}`$ vanishes there and the pseudo-arclength step must become vanishingly small for such a continuation step to be successful. Because of machine precision, this cannot occur, and so one must “detune” the system to continue across the cusp.

In general, Hamiltonian systems are described by

$$\dot{x}=JH(x,\lambda ),\text{ }x^{2n},$$
(19)
where

$$J\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right)$$
(20)
is the canonical $`2n\times 2n`$ symplectic matrix. One can detune the system by using a small perturbation parameter $`ϵ`$ as follows. The perturbed dynamical system,

$$\dot{x}=JH(x,\lambda )+ϵH(x,\lambda ),$$
(21)
is no longer Hamiltonian, but the perturbation was constructed so that the locations of all equilibrium points are preserved. With this detuning, the saddle-center bifurcation becomes a saddle-node bifurcation (the eigenvalues of the stationary point are now of the form $`a\pm \sqrt{b^{}}\text{ }(a0)`$ rather than of the form $`\pm \sqrt{b}`$). One can then continue $`\lambda `$ past this point using AUTO (without utilizing the HomCont package). Using this technique, one can compute the value of $`\lambda `$ at which the saddle-center bifurcation occurs as well as the cusp point of the homoclinic orbit corresponding to that value. Moreover, once one has successfully continued past the cusp, one can simply let $`ϵ0`$ and thereby work with the system when $`\lambda >\overline{\lambda }`$. This method of continuation is useful as long as one needs to get past the cusp rather than do computations at the cusp itself. For the present system, we used $`ϵ<0`$, since in that case the equilibria that are centers for $`ϵ=0`$ become stable spirals. The continuation curve (in $`\lambda `$) is shown in Figure 3.

In general, AUTO has difficulties near cusps. As with DSTool, one can approximate the cuspidal homoclinic orbit using AUTO. In order to do this, one provides initial values for the HomCont continuation parameters (corresponding to the initial point in the $`(a,P)`$-plane) to the right of the saddle point $`(a^{},0)`$. This allows AUTO to continue along the homoclinic orbit for values of $`\lambda `$ closer to $`\overline{\lambda }`$ than if one had started as close to the saddle point as machine precision would allow. For the present example, the closest accurate plot we obtained was for $`\lambda =0.1887`$. The right homoclinic orbit is pictured in Figure 4 and the left one is pictured in Figure 5. Observe that the one on the right looks like it has a cusp at the saddle point and that the left one is very small. As the saddle-center is approached, the left homoclinic orbit shrinks to a single point and the right one becomes more cusplike.

## 4 Chaotic Case: Presence of Coupling

There are other methods one can use to perform analysis near the cusp. One can, for example, use shooting methods. (AUTO uses a collocation method.) Saddle-centers are a degenerate case of Takens-Bogdanov (TB) bifurcations, so one can add dummy parameters and analyze the cuspidal loop by computing the locus of a TB bifurcation in parameter space, moving along the TB curve until one finds a Hopf bifurcation of another equilibrium, and then following the evolution of the periodic orbit created in the Hopf bifurcation as the parameters follow the TB bifurcation curve. If a cuspidal loop exists, this method will find it when the periodic orbit collapses into the cusp point.

For the chaotic case, the equilibrium points satisfy $`x=P=y=0`$, $`z=\pm 1`$, and $`\frac{V}{a}=\frac{2}{a^3}(ϵ_+\pm ϵ_{})`$, where the factor of $`\pm 1`$ in the last quantity corresponds to the equilibrium value of $`z`$. The eigenvalues for the equilibria of the chaotic system (Eqs. 26) are of the form

$$\sigma =0,\text{ }\pm \frac{\sqrt{2M(\alpha \pm \sqrt{\beta })}}{2Ma^2},$$
(22)
where both $`\alpha `$ and $`\beta `$ have terms whose signs depend on whether $`z`$ is $`1`$ or $`1`$. The quantities $`\alpha _\pm `$ and $`\beta _\pm `$ are given by

$$\alpha _+=a^4\frac{^2V}{a^2}+\omega _0^2M4\mu ^2ϵ_{}+6(ϵ_++ϵ_{}),$$
(23)
$$\alpha _{}=a^4\frac{^2V}{a^2}+\omega _0^2M+4\mu ^2ϵ_{}+6(ϵ_+ϵ_{}),$$
(24)
$$\beta _+=I_1+I_2^++I_3^++I_4^+,$$
(25)
where

$$I_1a^8\left(\frac{^2V}{a^2}\right)^2,$$

$$I_2^+\frac{^2V}{a^2}\left(12a^4ϵ_+2\omega _0^2Ma^4+12a^4ϵ_{}8\mu ^2a^4ϵ_{}\right),$$

$$I_3^+16\mu ^4ϵ_{}^28\omega _0^2M\mu ^248\mu ^2ϵ_{}^212\omega _0^2Mϵ_++72ϵ_+ϵ_{}12\omega _0^2Mϵ_{},$$

$$I_4^+48\mu ^2ϵ_+ϵ_{}+\omega _0^4M^2+36(ϵ_+^2+ϵ_{}^2),$$
(26)
and

$$\beta _{}=I_1+I_2^{}+I_3^{}+I_4^{},$$
(27)
where $`I_1`$ is as before,

$$I_2^{}\frac{^2V}{a^2}\left(12a^4ϵ_+2\omega _0^2Ma^412a^4ϵ_{}+8\mu ^2a^4ϵ_{}\right),$$

$$I_3^{}16\mu ^4ϵ_{}^2+8\omega _0^2M\mu ^248\mu ^2ϵ_{}^212\omega _0^2Mϵ_+72ϵ_+ϵ_{}+12\omega _0^2Mϵ_{},$$

$$I_4^{}48\mu ^2ϵ_+ϵ_{}+\omega _0^4M^2+36(ϵ_+^2+ϵ_{}^2).$$
(28)
Analogous to the planar system, only a generalization of saddle-center bifurcations can occur. As the energy is increased, a bifurcation corresponds to an increase in the dimension of the center manifold by two as a pair of real eigenvalues of opposite signs becomes a pair of pure imaginary eigenvalues. For a double-well potential, the only possibility is a jump in the dimension of the center manifold from three to five. For more complicated potentials, there may be parameter values with a one-dimensional center manifold. As before, one can determine the location of this bifurcation by finding the equilibria for which $`H(a,P)`$ has a double root. One again finds that the equilibrium point $`(0,0,\pm 1,\overline{a},0)`$ at the bifurcation satisfies

$$\frac{V}{a}(\overline{a},0)=\frac{\overline{a}}{3}\frac{^2V}{a^2}(\overline{a},0).$$
(29)
Finding the parameter values at which this occurs is not the simple issue it was in the previous case. In the planar case, $`\lambda `$ was a probabilistic weighting of two energies, so it could be varied continuously past the bifurcating value $`\overline{\lambda }`$. However, the parameters $`ϵ_1`$, $`ϵ_2`$ in the present system are part of a discrete energy spectrum, and so one cannot vary them continuously. In practice, therefore, this model predicts superposition states on each side of the bifurcation, but one does not expect to observe the system at a bifurcating value of $`(ϵ_1,ϵ_2)`$ because the set of all $`(ϵ_1,ϵ_2)`$ that correspond to bifurcating values has measure zero. (One could vary $`V_2`$ continuously if one wanted to examine bifurcations corresponding to a change in the quartic potential. If $`V_2`$ is negative and sufficiently small for a given $`V_4`$ or if it’s positive, the system will not exhibit a saddle-center.) Additionally, numerical observations indicate that the bifurcation occurs at low energies (corresponding to superpsoitions of low quantum number states), so that for a given billiard system, most superposition states will exhibit an evolution with a five-dimensional center manifold. The bifurcation under study may thus occur as one considers superpositions of increasingly excited states of the quantum billiard.

Numerical investigations have shown that homoclinic orbits must exist for this five-dimensional system (Fig. 6). Slight perturbations away from the homoclinic orbits lead to homoclinic tangles (Fig. 7), which are traditionally analyzed using symbolic dynamics. The details have not been worked out, but the existence of a homoclinic orbit for the present case can be shown analytically as follows. There is only one saddle point, so any saddle connection would have to be a homoclinic connection. By Hamiltonian symmetry and the existence of a center to the right of the “saddle-like” (in the sense that it has one-dimensional stable and unstable manifolds) equilibrium, there must exist a structure to its right that looks like a homoclinic orbit when projected into the $`(a,P)`$-plane. It may not be a homoclinic orbit, because one must consider the value of $`(x,y,z)`$ at the point $`(\stackrel{~}{a},0)`$ where the projection intersects the $`a`$-axis. One thereby considers the 2-sphere $`S^2`$ and two trajectories on it that start at the same point. One curve begins at $`t=\mathrm{}`$ from the stable manifold, and the other starts at $`t=\mathrm{}`$ from the unstable manifold. One looks at the intersection of one trajectory at time $`t`$ with the other at time $`t`$. If one can prove that such a point exists at a time $`T`$, then repeating the argument shows that such an intersection occurs at infinitely many times. To complete the proof, one must show that $`P=0`$ at one of these points.

Note that the homoclinic orbit for the chaotic system (Eqs. 26) emanates and terminates from a nonhyperbolic equilibrium, which increases the difficulty of numerical studies. As with the planar system discussed earlier, there is also a cusplike structure as the stable and unstable manifolds coalesce along the $`a`$-axis.

As an anology, consider two undamped springs, one with a linear restoring force ($`F_1=kx`$) and a stiffer spring with a cubic one ($`F_3=kx^3`$). These two springs (with mass $`m=1`$) are described, respectively, by the differential equations

$$\ddot{x}+kx=0,$$
(30)
and

$$\ddot{x}+kx^3=0.$$
(31)
One observes the same stiffness phenomena for these spring systems as we did for the integrable case of the vibrating billiard. (For example, the phase space trajectory of the cubic spring (Eq. 31) has a larger radius of curvative than the analogous one in the linear spring system. There is a correspondence with the other properties we discussed as well.) By analogy with mass-spring systems, we thus conclude that it is reasonable to consider the “stiffness” of the potential in which a quantum billiard resides as an object of interest.

The analogy with spring stiffness carries through in the chaotic case as well as in the planar case. For the same initial conditions, we consider the potentials

$$V(a)=V_2(aa_0)^2$$
(32)
and

$$V(a)=V_4(aa_0)^4,$$
(33)
where for ease of comparison, $`V_2=V_4`$. As shown in Fig. 8, the quartic potential gives trajectories with a larger radius of curvature than those in the quadratic potential. It thus makes sense to consider a potential’s stiffness in the chaotic case as well as in the integrable one.

## 5 Conclusions

We showed that saddle-centers are the only type of bifurcations that can occur in one dov quantum billiards. When there is no coupling in the superposition state, we showed how to analyze the resulting planar system analytically and numerically. Considering a double-well potential as a representative example, we demonstrated how to locate bifurcations in parameter space. We also discussed how to approximate the cuspidal loop using AUTO as well as how to continue past it by detuning the dynamical system. We also mentioned a shooting method for a more detailed analysis of the cuspidal loop. Moreover, we showed that when there is coupling, the resulting five-dimensional system—though chaotic—has a similar underlying structure. We verified numerically that both homoclinic orbits and cusps occur. Small perturbations of the system reveal homoclinic tangles that typify chaotic behavior.

## 6 Acknowledgements

We would like to thank Joel Ariaratnam, Richard Casey, John Guckenheimer, and Steven Wirkus for useful discussions concerning this project. Additionally, we express our gratitude toward Alan Champneys for several productive suggestions regarding the numerics, Alejandro Rodríguez-Luis for providing a preprint of his manuscript, and Paul Steen, whose guidance for this work as a project for ChE 753 (on which this paper is based) was particular valuable.

### Figure Captions

Figure 1: Approximate cuspidal homoclinic orbit. The initial point used was $`(0.7886751,0.001)`$, just above the equilibrium.

Figure 2: Homoclinic orbits emanating from $`(0.8916637,0)`$ for $`\lambda =0.15`$. The label 12 refers to the right homoclinic orbit, and the label 6 refers to the left one.

Figure 3:Continuation of the detuned system in the parameter $`\lambda `$.

Figure 4: Right homoclinic orbit for $`\lambda =0.1887`$.

Figure 5: Left homoclinic orbit for $`\lambda =0.1887`$.

Figure 6: Poincaré section projected into the $`(a,P)`$-plane demonstrating that there must exist a homoclinic orbit.

Figure 7: A homoclinic tangle projected onto the $`(x,y)`$-plane of the Bloch sphere.

Figure 8: Chaotic Poincaré maps in the $`(a,P)`$-plane for billiards in both a quartic and a quadratic potential. The plot from the quartic potential is the one with a smaller maximum value for the distance $`a`$. Notice that the trajectory corresponding to the quartic potential generically has a larger radius of curvature.


format: text

text:
# Observations of flat-spectrum radio sources at 𝜆850𝜇m from the James Clerk Maxwell Telescope I. April 1997–April 2000

## 1 Introduction

The James Clerk Maxwell Telescope (JCMT) is a 15-m diameter radio telescope situated on the 4.2 km peak of Mauna Kea in Hawaii. The facility is equipped with a suite of common-user instruments, including the Submillimetre Common-User Bolometer Array, SCUBA (Holland et al. 1999). The data reported in this paper comprise 2238 flux measurements from pointing observations of 65 flat-spectrum extragalactic radio sources. The observations have been reduced using the automatic SCUBA data reduction pipeline (Jenness et al. 2001) encompassing on-line atmospheric extinction correction using quasi-continuous monitoring data from the Caltech Submillimetre Observatory (CSO). Lightcurves for the most frequently sampled sources are presented.

Astronomical programmes aimed at understanding the emission mechanisms in this class of source are inevitably data limited, especially for time-dependent multi-wavelength variability studies. The millimetre-submillimetre region is very important for these, and so these data, like those at a longer wavelength from IRAM (Steppe et al. 1988, 1992, 1993; Reuter et al. 1997) will be a most valuable database for reference and use. Furthermore, a significant number of these sources have no previously published observations at submillimetre wavelengths.

## 2 The sources

Flat-spectrum radio sources are AGNs with relativistic jets. Monitoring of the emission gives insights into the jet emission, and can be particularly revealing during flaring behaviour. Although the sources fall into categories labelled by classes such as radio-loud quasars and BL Lacertae objects, for the purpose of this paper these sub-classifications are not important. Rather, this paper is a catalogue of observations of a heterogeneous sample of flat-spectrum extragalactic radio sources covering a three-year period; their commonality being that they are all relatively bright at 850 $`\mu `$m and that they are visible from the JCMT.

The flat-spectrum radio sources have synchrotron spectra that extend into the submillimetre with an index, $`\alpha `$, that is typically between 0.5–0.7 ($`S\nu \nu ^\alpha `$) at submillimetre wavelengths (e.g. Gear et al. 1994). The flat nature of the radio spectrum is believed to be due to the superposition of a number of self-absorbed synchrotron emitting blobs within the jet that emanates from the central engine. This view is supported by very high spatial resolution radio studies that have been able to deconvolve a number of these components (e.g. Marscher 1988). The highest frequency synchrotron component, which often occurs in the submillimetre, comes from the innermost part of the jet according to most models (e.g. Marscher & Gear 1985). Extensive studies in the millimetre/submillimetre regime have concentrated on monitoring the detailed flux variations and spectral evolution of this component (e.g. Robson et al. 1983; Valtaoja et al 1998; Stevens et al. 1996, 1998)

Another important aspect of the synchrotron emission in these sources is the ability of the relativistic electrons to interact with the local photon field to produce higher frequency radiation, such as X- and $`\gamma `$-rays, through the inverse-Compton process (e.g. Sikora, Begelman & Rees 1994). Progress in determining the precise mechanism of high-energy photon production is currently only achievable through temporal correlations of the emission between what is believed to be the low-energy seed photons and the high-energy progeny. One of the major problems in this research, however, is the scarcity of data at all wavelengths. Although some coordinated campaigns involving satellites have been undertaken (e.g. McHardy et al. 1999), nevertheless, the lack of data is a crucial element in stifling progress. The data in this paper will go some way to assisting these programmes in terms of the submillimetre flux record.

Although astronomical programmes on the JCMT have been undertaken with the express purpose of monitoring these variable extragalactic sources, especially blazars, the data contained in this paper are not of this kind. Blazar monitoring was routinely carried out using the previous single-pixel photometer, UKT14 (Stevens et al. 1994) but has declined with the advent of SCUBA. For the monitoring programmes the source samples were carefully selected, however, the sources included here make up a very ad hoc sample and the sampling rate varies significantly. Some sources are observed much more frequently because they are either very bright and useful for set-up and monitoring observations of pointing and focusing (e.g. 3C 279), they are well-placed for large observing programmes (e.g. 1308+326 for the Hubble Deep Field North observations - Hughes et al. 1998), or they are one of the few pointing sources in that part of the sky.

## 3 The observations

Although the rms absolute pointing of the JCMT is typically 1.5 arcsec, it is almost universal practice to use a local pointing source for accurate source registration. Pointing sources should be bright and compact with respect to the 850 $`\mu `$m beam ($``$14.5 arcsecs). Because there are few planets or secondary calibration sources (see Sandell 1994) available, most pointing sources are bright, flat-spectrum, radio sources. These give excellent sky coverage and allow pointing observations to be undertaken in times of order 1 minute. The observations in this paper are entirely from these pointing observations using the continuum camera, SCUBA, and at a single wavelength of 850 $`\mu `$m. The observation technique is described below.

Because the SCUBA pixels are fed by feedhorns with spacing of 2F$`\lambda `$, where F is the focal ratio and $`\lambda `$ is the wavelength of observation, the array does not instantaneously fully sample the image plane. In order to obtain a fully sampled image and hence a ‘picture’, the image is moved across the array in a precisely determined pattern. This is achieved using the telescope secondary, which is rapidly moved through a jiggle pattern. For a single wavelength, such as 850 $`\mu `$m considered here, this is a 16 position jiggle, each jiggle position offset by $``$6 arcsec. A complete jiggle map includes the telescope nod position to remove imbalances and is made up of 2$`\times `$16 jiggles and takes 32 seconds to complete. Successive jiggle maps are added together to build-up the required signal-to-noise for the final image, a 2.3 arcmin diameter picture of the sky.

Pointing is a 16 position jiggle-map where the source is fitted with a centroid to determine its precise location and hence the telescope local offsets. Although a default pointing observation consists of two complete jiggle-maps, the brightest sources, like 3C 279, are often observed with one; conversely as many as ten maps have been used for the faintest sources. Taking the measured range of noise equivalent flux density (NEFD) measured for SCUBA’s 850 $`\mu `$m filter ($`75110`$ mJy Hz<sup>-1/2</sup>) the observed rms noise level on the maps is thus expected to vary between about 20 and 50 mJy. Therefore, providing they can be well calibrated, pointing observations represent a valuable archive of data for the variable extragalactic sources.

All observations from the JCMT are archived, both on-site and at the Canadian Astronomical Data Centre (CADC). The next section discusses the processing of the pointing data from this archive to determine the fluxes for the extragalactic radio sources. The paper by Jenness et al. (2001) discusses the general case of the pipeline process of data extraction and reduction in much greater detail.

## 4 The data reduction technique

### 4.1 Introduction

Converting the pointing observations into a useful source of data requires three critical steps. The first is to extract the data in such a way that they are not affected by local parameters (such as telescope focus condition). The second is to correct for the atmospheric extinction, while the third is the flux density calibration.

Because of the large quantity of data, it was essential that an automated technique was used. All data were thus reduced using the ORAC-DR data reduction pipeline developed for SCUBA and the JCMT (Jenness et al. 2001). This procedure, along with the data analysis techniques are described below.

### 4.2 The pipeline: data processing and calibration

ORAC is the Observatory Reduction and Acquisition Control system developed for the United Kingdom Infrared Telescope (UKIRT). ORAC-DR (Economou et al. 1999; Jenness & Economou 1999) is the modular and highly flexible Data Reduction component whose primary design goals were to simplify data reduction whilst observing and to provide near-publication quality results. The data reduction pipeline for SCUBA is based on ORAC-DR, allowing the same approach to be applied to the processing of archival data, and time-dependent data to be additionally processed with minimal effort.

ORAC-DR is a recipe driven system that provides prescriptions for all of the SCUBA observing modes. Since the data are from pointing observations our adopted recipe is simplified over the generic mapping case in that it can be assumed that a source is present in the observation, it is close to the centre of the array and is point-like. Specifically, each pointing observation was processed in the following way: (i) correction for instrument flat-field; (ii) spike removal; (iii) extinction correction; (iv) skynoise removal (see Jenness, Lightfoot & Holland 1998); (v) rebinning the bolometer data onto a regular grid; (vi) fitting and removal of any residual image gradient; (vii) aperture photometry; (viii) flux calibration. These steps are all fully described in Jenness et al. (2001) but for the purposes of this paper, the three steps described in the introduction to this section are described in greater detail below.

(a) Correction for local effects

Pointing observations are always made at the beginning of the night prior to focusing the telescope and at subsequent times during the night before refocussing. Therefore, there is a clear potential, indeed certainty, for some of the pointing maps in the archive to be out-of-focus. There are two solutions to this problem, one is to ignore all pointing observations taken immediately before a focus, the second is to ensure that the data extraction technique is relatively insensitive to out-of-focus conditions. We have taken the latter course of action because it is also the case that the telescope may have drifted out of focus during the night due to insufficient focus attention by the observer. The solution is to use what is equivalent to aperture photometry. For all observations we use an aperture of 40 arcsec centred on the source which removes all but the most extreme of the out-of-focus conditions (see section 5).

(b) Extinction correction

Observations at 850 $`\mu `$m are rather sensitive to the extinction correction given by S(z)$`\mathrm{exp}^{\tau \mathrm{sec}(z)}`$, where $`\tau `$ is the zenith opacity and ranges from abo ut 0.12 to 0.6 for our data set. Therefore, as the first step in calibration, great care needs to be taken in assigning the correct extinction value to the data.

The extinction for all nights in the SCUBA archive has been determined from a combination of SCUBA skydips and the CSO quasi-continuous extinction measurements at 225 GHz ($``$1.3 mm; see Jenness et al. 2001 and Archibald et al. 2001). The resulting data are then fitted by a polynomial of varying degree depending on the variation during the night. When the atmosphere is very unstable the polynomial fits fail. The astronomical data are windowed-out during this period and no extinction values are available for the pipeline to attempt a data reduction. While this removes potential data, more importantly it limits the inaccuracy of using poorly known extinction corrections during periods of large transmission variability.

However, during the early investigations of the accuracy of the pipeline reduction process, incorrect assignments of the extinction were still found to be present. On further inspection it was found that this only appeared to occur close to the time when the atmosphere was in a period of instability, adjacent to the windowed-out data. One explanation is that the CSO monitor operates at a fixed azimuth, and so a pointing observation made just before, or soon after a period of instability, could be at an azimuth where the extinction is not well represented by the CSO determination. There are a number of ways to attempt to compensate for this. The simplest is to increase the width of the ‘no extinction correction possible’ window. While this would increase the quality of the remaining data, it would also remove a substantial dataset from the archive. The second solution is the most complex but potentially most elegant. This is to attempt to understand the weather pattern above Mauna Kea from satellite images of the water vapour and model the results with a time-azimuthally dependent extinction term and to feed this into the pipeline so that the azimuth of the source at the time can be modelled against the variable extinction. The third method is the pragmatic solution, which looks for all such potential problem times and treats the data by hand on a case-by-case basis.

We have rejected the widening of the data-rejection window, while effort is not available to pursue the second method. Therefore, we have opted for the pragmatic solution and have inspected all the data by eye for obvious problems. Where these are found the data-point is removed rather than an attempt at a correction made.

When the pipeline is released with the archive the user will be able to make an assessment of the quality of the night. It is clear that using the data reduction pipeline without taking care to inspect the extinction pattern for the night can lead to serious errors in flux determination.

(c) Calibration

Calibration of images (or maps) at submillimetre wavelengths is best accomplished with observations of point-like sources of known flux, made close in time to those of the target source. However, our automated procedure demands that we take a more global approach. Calibration is based on long-term observations of the two sources CRL618 and Uranus, thereby producing a time-dependent flux calibration factor (FCF) for SCUBA. These have been very stable over the period, changes due to upgrades of the instrument are clearly seen and reflected in the changing FCF. Using this technique an overall accuracy of around 5 per cent is obtained (Jenness et al. 2001).

## 5 Post pipeline processing

The output data from the pipeline for any source are first averaged over an individual night; there has been no attempt to determine variability within a single night. The nightly averaged data are first viewed to determine whether there are any obviously erroneous points. While this is easy to accomplish for a calibration source, or a source that is not variable, for these variable extragalactic radio-sources, this introduces a level of subjectivity. For example, when SCUBA polarimetry is undertaken there is no flag in the data-header indicating that the polarimeter was in place. This leads to a drop in the flux by a factor of about two. Although attempts have been made to retrospectively ‘flag’ all such data in the archive by hand, it is possible that there are still some ‘unflagged’ data-points remaining. Indeed, during this work this was discovered for a few data-points, which were easily spotted and removed.

Two further reasons why individual data-points may be erroneous are that the extinction has still been incorrectly attributed, or, there were problems with the SCUBA array. For the latter case, inspection of the individual data images reveals any obvious problems such as noisy pixels. These data-points are also removed. Neither of these factors were a common occurrence.

Finally, heating of the dish during early evening can cause notable variations in the derived flux calibration factor (FCF), which can remain unstable until late in the evening. Such excursions affect some 5–10 per cent of the data, and in the most extreme cases can increase the FCF by 50 per cent from its nominal value, even though the signal is integrated over a 40 arcsec diameter aperture. During these periods a significant amount of the received flux is removed from the main beam and spread out into the error lobes. However, because the resulting images are essentially out of focus, the worst cases are easily removed by visual inspection of the maps.

## 6 Data uncertainties

In determining the overall uncertainty of the measurement, the flux calibration uncertainty needs to be added in quadrature to that derived from the signal-to-noise ratio of the observation. In principle, if multiple observations of the same source are made on any one night then the signal-to-noise can be calculated from the scatter in the signals - as is done for photometry and under the assumption of a non-variable source. This approach would have the added benefit of including the dish-induced errors for any given data-set, allowing the use of the standard 5 per cent flux calibration uncertainty discussed above. In practice, however, although a source is often observed more than once during a night, not enough data-points are available to calculate a reliable standard error.

We have thus proceeded as follows. A conservative flux calibration uncertainty of 10 per cent is adopted. This is based on the scatter of the FCF values calculated from the CRL618 and Uranus observations, including the ‘bad’ data-points described above. Because of the need for an automated approach we chose to split the data into flux bands and then use a typical signal-to-noise towards the lower end of the range of those calculated for the individual maps. Our total uncertainties are thus: 10 per cent for fluxes $`>500`$ mJy; 12 per cent for 400–500 mJy; 13 per cent for 300–399 mJy; 14 per cent for 200–299 mJy and 17 per cent for 100–199 mJy. Note that the signal-to-noise for any individual map is $`>5\sigma `$, and for the majority of the data set is large enough to be negligible in comparison to the 10 per cent flux calibration uncertainty.

## 7 The data

The observed sources are presented in Table 1 which, for each source, lists the date of the first and last observation and the number of nights on which the source was observed. Owing to space constraints the flux density measurements are not listed and are available in the electronic version of this paper or from the Centre de Données astronomiques de Strasbourg (CDS). Table 3 provides a subset of the data as an example. The lightcurves for the best-sampled sources are presented in Fig. 1.

## 8 Conclusion

Due to the nature of the telescope and observing techniques, regular pointing observations are undertaken by the JCMT as part of the normal observing pattern. The pointing targets are mainly flat-spectrum radio sources. The resulting data are a potentially valuable resource for monitoring studies, especially multifrequency studies of active galaxies looking for links between the various spectral signatures of primary and secondary emission mechanisms. However, before they can be useful the data must be reliably calibrated. We have shown that with care, the data can be calibrated to an accuracy that is extremely useful. Furthermore, the development and refinement of the pipeline process means that it is a relatively simple task to extract and calibrate all these data and make them available to the scientific community. This is the first paper in what will be an ongoing release of data that can be used in a variety of studies.

## ACKNOWLEDGMENTS

The James Clerk Maxwell Telescope is operated by the Joint Astronomy Centre in Hilo, Hawaii on behalf of the parent organisations PPARC in the United Kingdom, the National Research Council of Canada and The Netherlands Organisation for Scientific Research. J.A.S. acknowledges support from PPARC.

We acknowledge the support software provided by the Starlink Project which is run by CCLRC on behalf of PPARC.


text:
# Strong Shocks and Supersonic Winds in Inhomogeneous Media

## 1. INTRODUCTION

Mass flows are important in many astrophysical systems from stars to the most distant active galaxies. Virtually all mass flow studies focus on homogeneous media. However, the typical astrophysical medium is inhomogeneous with the inhomogeneities (clumps) arising due to initial fluctuations of mass distribution, the action of instabilities, variations in the flow source, etc. The presence of inhomogeneities can introduce not only quantitative but also qualitative changes to the overall dynamics of the flow.

Active galactic nuclei represent one astrophysical site where a shock / wind interaction with a system of inhomogeneities may take place. Practically all current models of AGNs agree that the emission-line clouds in BLRs are essential for explaining the observed properties of AGNs (Urry & Padovani 1995; Elvis 2000). However, despite the fact that the nature of BELR clouds is an open question, any self-consistent model of the AGNs should properly account for properties of BELR cloud interaction with mass outflow. In particular, such self-consistency is important in the context of cloud survival time and cloud displacement prior to destruction. Even when clouds are magnetically confined and stabilized against evaporation (Rees 1987), disruptive action of outflows may prevent cloud survival over the dynamically significant AGN timescales.

(Klein, McKee, Colella 1994) (hereafter KMC) addressed the problem of shock-cloud interaction, providing a detailed description of the dynamics of a single, dense, unmagnetized cloud interacting with a strong, steady, planar shock. That work gives an excellent introduction to the subject, in particular the description of the astrophysical significance of the problem of a shock wave interacting with a dense cloud (see also (Gregori et al. 2000)). In this work we investigate the general properties of strong shock / supersonic wind interaction with a system of embedded clouds and determine the key quantities governing the evolution of such systems.

## 2. RESULTS

We have numerically investigated the interaction of a strong, planar shock wave with a system of dense inhomogeneities, embedded in a more tenuous and cold ambient medium. Our code is the AMRCLAW package, which implements an adaptive mesh refinement algorithm for the equations of gas dynamics (Berger & Oliger 1984; Berger & Jameson 1985; Berger & Colella 1989; Berger & LeVeque 1998) in two dimensions. We have assumed constant conditions in the global postshock flow constraining the maximum size of the clouds only by the condition of the shock front planarity. Our results are applicable to strong global shocks with Mach numbers $`3<M_S<1000`$. The range of the applicable cloud - unshocked ambient medium density contrast values is $`101000`$. Figure 1 shows a case of a strong shock and a supersonic post-shock flow interacting with an inhomogeneous system of fourteen identical clouds in regular distribution.

Cloud evolution due to the interaction with the global shock and the postshock flow has four major phases, namely *the initial compression phase, the re-expansion phase, the destruction phase,* and finally *the mixing phase*. Each image in Figure 1 roughly illustrates each of those phases.

The timescale we use to define time intervals in our numerical experiments is the time it takes for the incident shock wave to sweep across an individual cloud. This is called the *shock-crossing time,* $`t_{SC}=(2a_{max}/v_S)`$, where $`a_{max}`$ is the maximum cloud radius in the distribution.

We define the *cloud destruction time* $`t_{CD}`$ as the time when the largest cloud fragment contains less than 50$`\%`$ of the initial cloud mass. Typically, in our simulations $`t_{CD}24t_{SC}`$.

A simple model for the cloud acceleration during the first three phases, i.e. prior to its destruction, can be developed. We find the cloud velocity

$$v_C(t)=\{\begin{array}{cc}v_{PS}\left(1\left(\frac{t}{t_{SC}}a_1+a_2\right)^1\right),\hfill & t12t_{SC}\hfill \\ v_{PS}\left(1\left(\left(\frac{t}{t_{SC}}12\right)^2b_1+\frac{t}{t_{SC}}a_1+a_2\right)^1\right),\hfill & tt_{CD}\hfill \end{array}$$
(1)
Here, $`v_{PS}`$ is the unperturbed postshock flow velocity. For the cases of infinitely strong shocks, i.e. $`M_S\mathrm{}`$, the above coefficients have the values

$$a_1=1.8310^3;a_2=1.09;b_1=8.5110^5.$$
(2)
Equation (1) shows that the maximum cloud velocity is not more than $`10\%`$ of the global shock velocity and not more than $`13\%`$ of the postshock flow velocity.

Cloud displacement prior to its destruction can be described as

$$L_C(t)=\{\begin{array}{cc}a_0c_1\left(\frac{t}{t_{SC}}\frac{1}{a_1}ln\left(\frac{t}{t_{SC}}\left(\frac{a_1}{a_2}\right)+1\right)\right),\hfill & 0t12t_{SC}\hfill \\ a_0c_1\left(\frac{t}{t_{SC}}c_2tan^1\left(\frac{\frac{t}{t_{SC}}12}{\frac{t}{t_{SC}}c_4+c_5}\right)c_3\right),\hfill & 12t_{SC}tt_{CD}\hfill \end{array}$$
(3)
In the limiting case $`M_S\mathrm{}`$ the values of the coefficients $`a_1`$ and $`a_2`$ are defined in (2) and $`c_1=1.5;c_2=103.22;c_3=10.9;c_4=9.4310^2;c_5=112.56.`$

We find the maximum cloud displacement prior to its destruction, $`L_{CD}`$, is

$$L_{CD}=L_C(t_{CD})3.5a_{max}.$$
(4)
The results of our model are in excellent agreement with the numerical experiments. The difference in the values of cloud velocity and displacement between analytical and numerical results is $`<10\%`$.

The principal conclusion of the present work is that the set $`\mathrm{\Lambda }`$ of all possible cloud distributions can be subdivided into two large subsets $`\mathrm{\Lambda }_I`$, *thin-layer systems*, and $`\mathrm{\Lambda }_M`$, *thick-layer systems*, defined as

$$\begin{array}{c}\mathrm{\Lambda }_I:(\mathrm{\Delta }x_N)_{max}L_{CD},\\ \mathrm{\Lambda }_M:(\mathrm{\Delta }x_N)_{max}>L_{CD},\end{array}$$
(5)
where $`(\mathrm{\Delta }x_N)_{max}`$ is the maximum cloud separation in the system along the direction of the flow, or the cloud layer thickness.

Distributions from each subset exhibit striking similarity in behaviour (e.g. Figure 2). The systems containing from one to five clouds, arranged in a single layer, exhibit exactly the same rate of momentum transfer from the global flow to the clouds. The two fourteen cloud runs $`M_{14}`$ and $`M_{14r}`$ have different cloud distributions, different total cloud mass, different cloud sizes. Nevertheless, the rate of the kinetic energy fraction increase during compression and re-expansion is different from the single layer cases but is still the same for both fourteen cloud runs. Other global properties exhibit the same behaviour. We conclude that the evolution of a system of shocked clouds depends primarily on the total thickness of the cloud layer and the cloud distribution in it, as opposed to the total number of clouds or the total cloud mass present in the system.

The key parameters determining the type of evolution, are *the cloud destruction length* $`L_{CD}`$, defined above in (3), and *the critical cloud separation* transverse to the flow $`d_{crit}`$. The latter is defined by the condition that the time for adjacent clouds to expand laterally and merge into a single coherent structure is equal to the cloud destruction time $`t_{CD}`$. It can be described by the expression

$$d_{crit}=2a_0\left\{\frac{t_{CD}t_{CC}}{t_{SC}}\left(\frac{F_{c1}F_{st}}{\chi }\right)^{\frac{1}{2}}\left(\frac{3\gamma (\gamma 1)}{\gamma +1}\right)^{\frac{1}{2}}+1\right\}.$$
(6)
We should be able to determine, either from observations or from theoretical analysis, the thickness of the cloud layer $`(\mathrm{\Delta }x_N)_{max}`$. This determines the class of the given cloud distribution, $`\mathrm{\Lambda }_I`$ or $`\mathrm{\Lambda }_M`$. For distributions from the set $`\mathrm{\Lambda }_I`$ with average cloud separation $`\mathrm{\Delta }y_N>d_{crit}`$ evolution of the clouds during the compression, re-expansion, and destruction phases will proceed in the *noninteracting regime* and the formalism for a single cloud interaction with a shock wave (e.g. KMC) can be used to describe the system. On the other hand, if the cloud separation is less than critical, the clouds in the layer will merge into a single structure before their destruction is completed, and though the compression phase still can be considered independently for each cloud, evolution during the re-expansion and destruction phases proceeds in the *interacting regime*.

When the distribution belongs to the subset $`\mathrm{\Lambda }_M`$ it is necessary to determine the average cloud separation projected onto the direction of the flow $`\mathrm{\Delta }x_N`$ and compare it against $`L_{CD}`$: if $`\mathrm{\Delta }x_N>L_{CD}`$ evolution of the system can be roughly approximated as of a set of distributions from the subset $`\mathrm{\Lambda }_I`$ and the above “thin-layer case” analysis applies. If, on the other hand, $`\mathrm{\Delta }x_NL_{CD}`$ (especially if $`\mathrm{\Delta }y_N<d_{crit}`$) the system evolution is dominated by cloud interactions and a thin layer formalism is inapplicable.

Finally we consider mass-loading. Here our principal conclusion is that mass-loading is not significant in the cases of strong shocks and supersonic winds interacting with inhomogeneities whose density contrast is in the range $`101000`$. In part this is due to short survival times of clouds as well as the very low mass loss rates of the clouds even during the times prior to their destruction. Therefore, mass-loading in such systems is not likely to have any appreciable effect on the overall dynamics of the global flow.

The major limitation of our current work is the purely hydrodynamic nature of our analysis that does not include any consideration of the magnetic fields. As was discussed by KMC, cold dense inhomogeneities (clouds) embedded in more tenuous hotter medium are inherently unstable against the dissipative action of diffusion and thermal conduction. Although weak magnetic fields, that are dynamically insignificant up to the moment of cloud destruction, can inhibit thermal conduction and diffusion and, therefore, stabilize the system of clouds, those magnetic fields may become dynamically important due to turbulent amplification during the mixing phase (see (Gregori et al. 2000) for a three-dimensional study of the wind interaction with a single magnetized cloud). We intend to provide a fully magnetohydrodynamic description of the interaction of a strong shock with a system of clouds in future work.

### Acknowledgments.

This work was supported in part by the NSF AST-0978765, NASA NAG5-8428, and DOE DE-FG02-00ER54600 grants.
The most recent results and animations of the numerical experiments, described above and not mentioned in the current paper, can be found at
www.pas.rochester.edu$`/^{}`$wma.

## References

Berger, M.J., Colella, P. 1989, J. Comp. Phys., 82, 64

Berger, M.J., Jameson, A. 1985, AIAA J., 23, 561

Berger, M.J., LeVeque, R.J. 1998, SIAM J. Numer. Anal., 35, 2298

Berger, M.J., Oliger, J. 1984, J. Comp. Phys., 53, 484

Gregori, G., Miniati, F., Ryu, D., Jones, T.W. 2000, ApJ, 543, 775

Elvis, M. 2000, ApJ, 545, 63

Klein, R.I., McKee, C.F., Colella, P. 1994, ApJ, 420, 213

Pietrini, P., Torricelli-Ciamponi, G. 2000, A&A, 363, 455

Rees, M.J. 1987, MNRAS, 228, 47

Urry, C.M., Padovani, P. 1995, PASP, 107, 803


text:
# USA Observation of Spectral and Timing Evolution During the 2000 Outburst of XTE J1550-564

## 1 Introduction

XTE J1550$``$564 was first observed in 1998 September by the All Sky Monitor (ASM) on board the Rossi X-ray Timing Explorer (RXTE) when it began an outburst lasting approximately 8 months (Smith, 1998). XTE J1550$``$564 began a second outburst on 2000 April 2 (Masetti & Soria, 2000), lasting approximately 2 months. The source was detected for a third time in 2001 January but did not go into a full outburst (Tomsick et al., 2001b; Jain et al., 2001).

Recent optical observations of XTE J1550$``$564 have placed a lower limit of $`7.4\pm 0.7M_{\mathrm{}}`$ on the mass of the compact object (Orosz et al., 2001). This mass places the compact object well above the maximum mass for a stable neutron star and so provides compelling evidence that XTE J1550$``$564 contains a black hole. During its 1998–1999 outburst, apparent superluminal radio jets ($`v>2c`$) were observed by Hannikainen et al. (2001); however, the jet’s angle to the line of sight has not yet been determined. Radio and optical observations of the 2000 outburst also show evidence of jet formation (Corbel et al., 2001; Jain et al., 2001).

Complex timing and spectral behavior has been observed by RXTE in XTE J1550$``$564 during its two full outbursts. This behavior includes detections of three classes of low-frequency quasi-periodic oscillations (LFQPO) ($`<20`$ Hz) and several detections of high frequency QPOs (HFQPO) ($`>100`$ Hz) (Remillard et al., 2001; Miller et al., 2001; Kalemci et al., 2001). Further, color-color diagrams and hardness intensity diagrams of the 1998–1999 outburst showed separate spectral branches for each of the black hole states as well as correlations with quasi-periodic oscillations and other timing behavior (Homan et al., 2001).

During the 1998–1999 outburst, XTE J1550$``$564 exhibited all four identified black hole spectral states (Sobczak et al., 2000; Homan et al., 2000). In the 2000 outburst the source never achieved the high state (HS), going from an initial low/hard state (LS) to an intermediate state (IS) or a very high state (VHS) and returning to a final LS (Miller et al., 2001). Although previous authors have made a distinction between the IS and VHS, recent work on XTE J1550$``$564 has suggested that the IS and VHS are actually the same state at different X-ray flux levels (Homan et al., 2000). Therefore, for the remainder of this paper we will refer to this state as the IS. The transition from the initial LS to the IS was made on 2000 April 26 (MJD 51660). The transition from the IS back to the LS occurred sometime between May 13 (MJD 51677) and May 19 (MJD 51623) (Corbel et al., 2001).

In this Letter, we report on X-ray observations of the 2000 outburst of XTE J1550$``$564 made by the Unconventional Stellar Aspect (USA) Experiment on the US Air Force Advanced Research and Global Observation Satellite (ARGOS). For a detailed description of the USA experiment see Ray et al. (1999) and Shabad (2000). We present lightcurves and hardness ratios and track a low frequency QPO which appears in the initial and late stages of the outburst.

## 2 Observations and Data Analysis

### 2.1 Lightcurves and Hardness Ratios

The USA Experiment observed XTE J1550$``$564 at the rate of 2–8 times per day between 2000 April 14 (MJD 51648) and June 18 (MJD 51713). For the present investigation, 193 observations were used, from which we selected $`49`$ ks of data obtained in low-background regions. The data are time tagged, having 32 $`\mu `$s time resolution, and cover an energy range of approximately 1–17 keV in 16 pulse height analyzer (PHA) channels. In this work we do not make use of the lowest (channel 0) and highest (channel 15) PHA channel. We refer to PHA channels 1–14 ($`116`$ keV) as the total range.

To create the light curves shown in Figure 1, a first order background subtraction was made by averaging blank sky observations and then subtracting these values from the count rate. To determine the total error, the standard deviation of the average in the background was added in quadrature to the error on the count rate. USA data were then corrected for obscuration by the instrument support structure and the collimator response. The overall average light curve (normalized to the USA Crab counting rate) for the total range is shown in panel (a) of Figure 1. The circles are RXTE/ASM daily averaged data used to give the complete outburst profile (USA observations did not cover the first few days of the outburst). The USA data points shown are an average of several USA observations. The number of observations averaged was dependent on the observation spacing and signal to noise ratio.

The spectral characteristics of the outburst were studied by dividing the USA data into two energy bands, USA PHA channels 1–3 and channels 4–14. These two bands correspond to 1–4 keV and 4–16 keV, respectively. For the remainder of this Letter these bands will be referred to as the soft band count rate (SB) and the hard band count rate (HB). The motivation for choosing the specific energy ranges of the SB and HB came from properties of the hardness-intensity diagram (HID). During the outburst, XTE J1550$``$564 traces a cyclic pattern in the HID. Plots of count rates versus the total range show this cyclical structure for individual channels in the range 4–14, but not for channels 1, 2 or 3.

Panels (b) and (c) of Figure 1 show the SB and HB as a function of time during the outburst. The bottom panel of Figure 1 shows the evolution of hardness ratio using these energy bands. The HID in Figure 2 shows how the hardness ratio evolves with total count rate during the outburst. The hardness ratio is plotted only for USA data prior to MJD 51687, after which the signal to noise decreases to the point that the hardness ratio is not constrained. In Figure 2 and panel (f) of Figure 1, one point is shown for each point in panels (b) and (c) of Figure 1.

The design of the USA detector incorporated automatic gain stabilization hardware and frequent iron source energy calibrations were done while in orbit. We note that the USA channel to energy conversion varies slightly over the USA orbit; however, checks performed showed that this variation made no significant impact on the relevant features seen in the SB, HB, and hardness ratio. A further check of our spectral results was made by comparing our data to public RXTE/ASM data. Daily averaged ASM data were used to find the hardness ratio as a function of time and to make a HID. The ASM hardness ratio was calculated by dividing the sum of ASM B Band (3–5 keV) and C Band (5–12 keV) by A Band (1.5–3 keV). To try to emulate the ASM energy bands the USA PHA channels 3–11 (3–11.5 keV) were summed and divided by channel 2 (2–3 keV). The ASM data confirmed the hardness ratio observed with USA and the cyclic behavior in the HID.

### 2.2 Power Spectra: Low Frequency QPOs

A low frequency quasi-periodic oscillation (LFQPO) was ubiquitous during the rise of the outburst and during the decay of the outburst after the secondary maximum. In order to track the LFQPO evolution through the outburst, observations were grouped by day and frequently in sub-day groups (signal to noise ratio permitting).

Power spectra (see Nowak et al. 1999 and references therein) were calculated from these groups and averaged. The resultant power spectrum for each group containing a LFQPO was fit with a power law or broken power law and a Lorentzian for any observable QPO features. Fits were made in three energy bands: the total range, SB, and HB. In case of confusion by sub-harmonics, the strongest QPO feature was chosen as the primary LFQPO (see Remillard et al. 2001). In most cases no sub-harmonics were detectable. Panel (d) of Figure 1 shows how the centroid frequency of the LFQPO evolves during the outburst. The evolution of the rms amplitude for all three energy bands is shown in panel (e). Panels (d) and (e) show error bars calculated by allowing the $`\chi ^2`$ of the fit to vary by one. All error bars are given at the 68% confidence level.

## 3 Results

### 3.1 QPO Evolution and Correlation to State Changes

We observe LFQPOs between MJD 51648 and 51663 and between MJD 51675 and 51686 which vary in frequency between 0.24–7.19 Hz and 6.34–0.64 Hz, respectively. During the times of these detections the source is either in the LS or near the transition from one state to another. The LFQPO rms amplitude decreases rapidly at the state transition from the LS to the IS and then increases during the transition back to the LS, indicating that the mechanism for creating the LFQPO is suppressed in the IS. The LFQPO centroid frequency generally increases with increasing flux; the fractional rms amplitude is correlated with hardness ratio (see Figure 1).

During the IS, significant detections of HFQPOs (249–278 Hz) were made by RXTE between between MJD 51663 and 51675 (Miller et al., 2001). A 65 Hz QPO has been discovered by Kalemci et al. (2001) at MJD 51684.8. These HFQPO detections occurred during the periods where USA observed the LFQPO to be weakening or not detectable at all. It is interesting to note that the HFQPOs were observed to decrease in strength as a function of time in the IS (Miller et al., 2001). These observations point to an unfavorable interaction between the mechanisms for LFQPO and HFQPO production. We observe a decline in rms amplitude of the LFQPO near the LS/IS transition, which marks the approximate onset of the HFQPOs. Towards the end of the IS, near the IS/LS transition, the HFQPO weakens as the LFQPO once again becomes detectable. This trend continues in the last days of the outburst, when the LFQPO rms amplitude weakens and the 65 Hz QPO is detected.

These QPO features are qualitatively consistent with observations of XTE J1550$``$564 during the 1998–1999 outburst (Remillard et al., 2001), during which an “antagonism” between LFQPOs and HFQPOs was also observed. During that outburst, type “C” QPOs were observed when strong correlations were seen between the frequency and disk flux while the amplitude was observed to correlate with disk temperature (Remillard et al., 2001). These previous observations closely resemble what we see for the LFQPO observed during the 2000 outburst; thus we tentatively classify the LFQPO discussed here as a type C.

At MJD 51661.21 an anomalous QPO was detectable exclusively in the HB, in contrast to the primary LFQPO which appears in all energy ranges. In Figure 1 panel (d), this anomalous QPO is the highest frequency point and is marked with a large unfilled circle near MJD 51660. This QPO was detected between the primary LFQPO and its harmonic and has a frequency of $`7.19_{0.11}^{+0.12}`$ Hz, whereas the primary LFQPO is seen at $`4.71\pm 0.05`$ Hz and its harmonic is seen at $`9.75_{0.35}^{+0.4}`$ Hz.

### 3.2 Spectral Evolution

From Figure 1 it is clear that the 2000 outburst of XTE J1550$``$564 does not follow the canonical fast rise, exponential decay (FRED) outburst as would have been expected prior to the RXTE era. Now, with many more examples of well-observed soft X-ray transient (SXT) outbursts, it has become clear that few outbursts look like pure FREDs, and that the outburst profile can be very different at high energy than at low energy. Comparing the SB and HB lightcurves in Figure 1 (b) and (c), the two bands rise approximately in unison, but show very different behavior after the peak. With the exclusion of the secondary maximum the HB light curve shows a nearly symmetric outburst profile, while the SB rapidly rises, then decays approximately linearly. As is often observed in SXTs, the decay returns to its original path after the secondary maximum.

In the HID (Figure 2), we observe that the difference between the SB and HB lightcurves manifests itself as a cyclic structure that moves with time in a clockwise direction. The HID shows a rapid drop in the hardness ratio as the source enters the IS and a rapid increase as the final LS is realized. This type of spectral structure, which has been seen in several other sources (Smith et al., 2001), is a consequence of the spectrum being harder during the rise than during the decline.

### 3.3 Multiwavelength Correlations

Multiwavelength observations of XTE J1550$``$564 during the 2000 outburst have shown correlations between radio, optical, IR and X-rays. These observations have been compared with similar correlations seen in other black-hole candidates (BHCs), such as GX 339$``$4, and have led to the interpretation that the LS of BHCs is characterized by jets and the IS state is characterized by quenching of these jets (Corbel et al., 2001). In Figure 1, lines 1 and 8 mark the times of optical/IR maxima and lines 3 and 7 indicate significant radio detections made by ATCA while line 4 indicates no significant radio detection. Line 6 is the time of an optical/IR minimum. (Optical/IR observations were made by YALO; for complete details on the radio and optical/IR observations see Corbel et al. 2001 and Jain et al. 2001). Radio observations by MOST (not shown) made during the initial LS show a detection of 8–15 mJy at 843 MHz (Corbel et al., 2001). This indicates that the radio emission prior to line 3 is more than a factor of 50 stronger than at the time of the ATCA observation made at line 4 (Corbel et al., 2001). Together, these radio and optical/IR observations have been interpreted(Corbel et al., 2001; Jain et al., 2001) as evidence for the presence of jets during the initial and final LS. USA observations indicate that the decline in the HB, during the IS, is associated with a rapid drop in radio emission. The weak radio detection (line 3) occurs while the SB is still very near its maximum and the HB has dropped from its maximum by a factor of $`2`$ and continues on a rapid decline. At the point of no significant radio detection (line 4), the HB has greatly reduced its rate of decrease. The return of the radio signal in the final LS suggests that the mechanism of jet production is related to but not dependent on HB photon production mechanism. In future outbursts, it would be helpful to obtain a much more thorough radio coverage in order to pin down the nature of the relationship between the HB and jet quenching.

## 4 Discussion: Two Flow Models

The different behavior of the HB and SB lightcurves, together with the lack of a strict correlation between the QPO frequency and source flux, suggest a scenario involving a two component accretion flow. Models invoking two independent accretion flows have been described by Chakrabarti & Titarchuk (1995) and Smith et al. (2001), while van der Klis (2001) presents a model where the two parameters are the instantaneous and time averaged values of a single quantity.

Whether the two accretion flows are dependent or independent, there will be a transition radius at which the Keplerian disk is disrupted and forms a hot sub-Keplerian flow. This inner flow can be an advection-dominated accretion flow (Esin et al., 1997) or a postshock flow interior to the radius where two independent flows interact (Chakrabarti & Titarchuk, 1995). One may then consider models where the LFQPO is related to the Keplerian orbital period at the transition radius ($`r_{\mathrm{tr}}`$, measured in units of the Schwartzchild radius) between the thin disk and an inner flow. The relationship between the Keplerian orbital period and radius for a $`7.5M_{\mathrm{}}`$ black hole is $`R=195R_{\mathrm{Sch}}P^{2/3}`$, where $`R_{\mathrm{Sch}}`$ is the Schwartzchild radius and $`P`$ is the orbital period in seconds. Thus, the LFQPO evolving from $`0.2`$ Hz to $`5`$ Hz could be related to the transition radius changing from $`570R_{\mathrm{Sch}}`$ to $`67R_{\mathrm{Sch}}`$. When the disk inner radius gets small enough, Compton cooling becomes important and the hard Comptonized emission is suppressed while the black-body emission from the disk moves into the X-ray band, causing a spectral state transition. This agrees remarkably well with what was observed in this outburst.

In order for a two flow model to explain the spectral characteristics described in § 3.2, it should allow for spectral softening to occur while the overall source flux is dropping. For two independent flows, this type of spectral evolution occurs because changes in the radial flow may occur on a short (free-fall) timescale, while the effect of changes in the the disk accretion rate are delayed due to viscosity (Chakrabarti & Titarchuk, 1995; Smith et al., 2001). A model of two independent flows was applied to GRS 1758$``$258 and 1E 1740.7$``$2942 by Smith et al. (2001), both of which also showed spectral softening with decreasing flux. Using a dependent flow model, van der Klis (2001) suggests that if the count rate is dominated by the disk accretion rate and the spectral hardness is related to the inner disk radius (which is determined by the time-averaged accretion rate), the same cyclic structure in the HID will result.

Recent ideas on the role of jets in microquasars suggest a link between the jet and the corona (see Fender 2001 and references therein) and provide a theoretical basis for the observed association between the HB and the radio and optical/IR signals (§ 3.3). This association is suggestive of a mechanism by which jets are created in the LS and quenched in the IS. The LS is associated with a continuous radio-emitting outflow and the presence of the strong LFQPO. At the transition to the IS, discrete radio ejections seem to be common but the continuous jet ceases and the LFQPO mechanism is suppressed. This indicates that both the outflow and the LFQPO are characteristic of the LS and are related to the presence of a hot corona or sub-Keplerian flow.

We gratefully acknowledge useful discussions with Lev Titarchuk. We thank Mark Yashar for providing useful references. Work at SLAC was supported by department of Energy contract DE-AC03-76SF00515. Basic research in X-ray Astronomy at the Naval Research Laboratory is supported by ONR/NRL. This work was performed while RMB held a National Research Council Research Associateship Award at NRL. JDS is grateful to the NASA Applied Information Technology Research Program for support. This paper made use of quick-look results provided by the ASM/RXTE team (see http://xte.mit.edu).