text:
# APPLICATIONS OF TIME DOMAIN SIMULATION TO COUPLER DESIGN FOR PERIODIC STRUCTURESWork supported by the Department of Energy, contract DE-AC03-76SF00515, grants DE-FG03-93ER40759 and DE-FG03-93ER40793.

## I Introduction

Numerical simulation procedures for designing waveguide couplers to accelerator structures are described in  and an example of its application to the design of the input coupler for the NLC linac is given in . A coupler cavity is designed with the intent of providing a matched connection between a waveguide and a uniform accelerator structure with dimensions corresponding to those of the cell adjacent to the coupler cavity. A symmetric structure consisting of two coupler cavities (with associated waveguides) connected by a short section of accelerator structure (typically two cells worth) is modeled and subjected to a (let us assume single frequency) time domain simulation. The entire assembly is treated as a single structure with two wave guide ports. The coupler cell dimensions are adjusted until an apparent match is achieved, that is, until no reflection is experienced at the ports (the external matching condition). To eliminate the possibility that the match arises from a fortuitous cancellation between forward and backward waves within the accelerator structure, both the amplitude and phase of the accelerating field on the beam axis are observed and required to have the periodicity and phase advance properties appropriate to a pure traveling wave (the internal matching condition). As a check one may add a cell to the accelerator structure and see whether all these conditions are still satisfied. It is often the case that accelerator structures are slowly varying rather than uniform, in which case the input and output couplers are matched separately.

In practice the procedure (we refer to it as the standard procedure) has been quite time consuming, involving trial and error rather than a systematic procedure to simultaneously satisfy both the internal and external matching conditions. Another limitation arises from the fact that the method assumes that evanescent bands can be neglected but provides no procedure for demonstrating their absence.

In the next section we describe a new simulation procedure which has been found to be much more efficient, and also which provides information about the presence of evanescent bands. The basic elements of the method were briefly described in  in connection with the design of a coupler for the zipper structure. Because it has since replaced the old method for all of our coupler design work, a more complete presentation together with examples will be presented in the following sections.

## II The New Simulation Procedure

As in the case of the old standard procedure one applies a single frequency time domain simulation by driving the input port of a two port structure consisting of an input cavity, an $`N`$ cell periodic structure with period $`P`$, and an output cavity. Instead, however, of focusing attention on the S parameters of the structure as a whole, we direct our attention to the simulated accelerating field $`E_z(z,t)`$ evaluated along the beam axis. We assume a steady state has been reached, so that the subsequent time dependence can be expressed in terms of the complex $`E_z(z)`$ ($`E_c(z)`$ henceforth), obtained in the standard way by combining the simulated real fields at two times separated by a quarter period. Then from Floquet’s theorem (neglecting evanescent bands, losses, and an irrelevant overall phase factor)

$$E_c(z)=E(z)[exp(j\varphi (z))+Rexp(j\varphi (z))].$$
(1)
Here $`E(z)`$ is a real positive amplitude function with period $`P`$, and $`\varphi (z)`$ is a real phase function, periodic except for a cell to cell phase advance $`\psi `$. Thus

$$E(z\pm P)=E(z),\mathrm{and}\varphi (z\pm P)=\varphi (z)\pm \psi .$$
(2)
$`R`$ is a $`z`$ independent complex reflection coefficient. Note that one is free to shift $`\varphi `$ by an arbitrary constant with a compensating phase shift in $`R`$, since the overall phase of $`E_c`$ is irrelevant. This freedom corresponds to the choice of reference plane through some point $`z_0`$ where we take $`\varphi =0`$.

We now consider the quantities

$$\mathrm{\Sigma }(z)=F^+(z)+F^{}(z),\mathrm{and}\mathrm{\Delta }(z)=F^+(z)F^{}(z),$$
(3)
where

$$F^\pm (z)=E_c(z\pm P)/E_c(z).$$
(4)
Elementary algebraic manipulation leads to the relations:

$$2Cos(\psi )=\mathrm{\Sigma }(z),$$
(5)
$$Rexp(2j\varphi )=[2Sin(\psi )j\mathrm{\Delta }(z)]/[2Sin(\psi )+j\mathrm{\Delta }(z)].$$
(6)
We note that while the RHS of (5) is formed of $`z`$ dependent complex quantities, it nevertheless turns out to be real and $`z`$ independent. Similarly the absolute value of the RHS of (6) is also $`z`$ independent. Both these results should hold for all “allowed” $`z`$ values, i.e., values such that the three points $`z+P`$, $`z`$, and $`zP`$ all lie within the periodic portion of the structure simulated, and together they constitute a powerful constraint on the validity of the Floquet representation Eq. (1). Their failure beyond small numerical fluctuations or small deviation from steady state is evidence for the presence of evanescent bands. An example will be presented in the section on the Zipper structure.

It is noteworthy that these relations allow one to determine all the properties of the traveling wave solutions, including the functions $`\varphi (z)`$ and $`E(z)`$ from a simulation which contains an arbitrary mixture of forward and backward waves. Of particular importance is the fact that it gives the magnitude and phase of the reflection coefficient. In contrast to the old standard method, there is here only one matching condition to be satisfied, namely $`|R|=0`$. Typically match is achieved by varying two parameters in the coupler design. Once one has determined how the real and imaginary parts of $`R`$ vary with the parameters, one can choose linear combinations of changes which accelerate the process of converging to the origin in the complex $`R`$ plane . Because the phase of $`R`$ does depend upon the position of the reference point relative to the couplers, one naturally keeps it fixed while carrying out this process. Note that it is the output cavity that is matched by this procedure. While not necessary, it is usually convenient to construct a symmetric mesh. The input and output cavities are then the same, and the structure as a whole is matched when $`R`$ vanishes.

## III Applications and Examples

(a) The NLC four port output coupler

As an example of the principal features of the new method we use the new NLC four port output coupler cavity . The purpose of the four port design was to provide damping for those dipole modes that reach the end of the structure while also providing an output for the accelerating mode. These dipole modes are typically those which had been poorly damped because of decoupling of the last cells from the manifolds. The four port symmetry provides damping for both dipole mode polarizations and has the added advantage of eliminating quadrupole distortion of the coupler fields.

The design simulation was carried out with a three cell periodic structure, and results are illustrated in Fig. (1). Two cases are shown, one matched, the other not. The reflection coefficients $`|R|`$ as computed from Eq. (5) for the two cases are shown as functions of $`z`$. The allowed z values are those lying within the central cell, and one sees that for both cases $`|R|`$ is constant over that range. The real part of $`Cos(\psi )`$ is also plotted as a function of $`z`$. One sees that the two values are indeed constant over the allowed range, but contrary to expectations they differ somewhat from each other and from the expected value of one half. This is due to the fact that a different and coarser mesh than that used to determine the phase advance parameter was used for the time domain simulations. The two cases differ from one another because the parameter variations in the coupler associated with the matching procedure induce small but global changes in the meshing. It has been confirmed in a number of cases that there is good agreement between the phase frequency relation determined from single cell periodic boundary condition frequency domain calculations and that determined from the time domain method described here so long as the same mesh is used for both simulations.

(b) A Photonic Band Gap (PBG) structure

A coupler cell very similar to those of the SLAC structures has been designed for a PBG structure, that is, a cylindrical cell with a pair of symmetrically placed waveguide ports, a conventional beam pipe, and conventional beam iris coupling to the periodic PBG structure. The PBG cell structure  is a seven by seven square array of metallic posts aligned in the beam direction and terminated by metallic end plates, the cell cavity being formed by removing the central post. A circular aperture in the end plates, identical to that between the coupler cell and the adjacent PBG cell, provides cell to cell coupling and a path for the beam. A perspective representation of the four cell quarter structure used for the simulations is shown in Fig. (2). Also shown is the simulated electric field distribution, scaled logarithmically to enhance the visibility of weak field strengths. The figure illustrates the effectiveness of the PBG structure in confining the acceleration fields to the interior of the structure. The matching procedure worked well, and, as in the four port coupler above, there was no evidence for evanescent band contamination. Fabrication of an experimental model with 5 coupled PBG cells and complete with couplers is in progress at SLAC.

(c) The Zipper structure

The zipper is a planar accelerator structure described in . A 25 (counting the coupler cavities) cell W band model has been built, cold tested, and subjected to bead pull measurements as reported in . The design was governed by a decision to avoid bonded joints involving tiny structure elements such as the vanes which serve as cell boundaries and also form the beam iris. The coupler cell is a quarter wave transformer terminating in WR10 waveguide.

Early attempts at matching the coupler using the old standard method failed, and it was this failure which led to the development reported here. Matching using this method was accomplished by making use of a time domain simulation of a structure with 22 periodic cells. Fig. 3 shows the resultant $`ReCos(\psi )`$, $`ImCos(\psi )`$, and $`|R|`$ plots as computed from Eqs. (5) and (6). One sees large deviations from the expected $`z`$ independent behavior as one moves away from the center of the structure. This effect indicates a clear violation of Eq. (1). From the fact that the violation fades away as one moves away from the couplers indicates that the effect is due to the couplers generating an evanescent band, the nearby monopole band pointed out in . This example demonstrates how the method indicates the presence of evanescent band interference, and also how one can carry out the matching procedure even when it is present.

## IV The Coupler Field Enhancement Problem

Electrical discharge damage has been commonly observed in the coupler cells of accelerator structures and has been attributed to the field enhancement noted in simulations. We have taken advantage of our enhanced matching capability to initiate a study of this long standing problem. Exploration of the situation for the NLC coupler  showed that the largest enhancement occurred on the coupler side of the aperture of the iris separating the coupler from the adjacent cell with azimuthal maxima opposite the coupler waveguides and azimuthal minima 90 degrees away. This observation was consistent with the pattern of discharge damage . It is pointed out in  that the azimuthal variation is due to the quadrupole component introduced by the coupler waveguides and that the enhancement can be reduced by introducing a racetrack like modification of the coupler cell shape designed to eliminate it. This effect and its cure have been confirmed in our own studies of the NLC coupler. Two other modifications have also been explored. The simplest and most effective was simply to reduce the radius of the cell adjacent to the coupler. The effect for a 2$`\%`$ reduction is illustrated in Fig. 4 where it is seen that the field on the coupler cell iris is significantly less than that on the interior coupling irises. An undesirable consequence is a 10 degree phase advance deficiency in the modified cell. An even larger field reduction would be obtained by removing the quadrupole enhancement. We attribute the reduction to an increase in group velocity. The other modification consisted of enlarging the coupler iris combined with an increase in the adjacent cell radius chosen so as to preserve the cell phase advance, but the exploration of this effect is incomplete. Experimental investigation to determine whether such changes actually do reduce electrical discharge damage in the coupler is clearly needed.


format: text

text:
# Comment on: “A quantum approach to static games of complete information”

\[

\]

In a recent article , an interesting new scheme for quantizing games was introduced, and the scheme was applied to the famous game Battle of the Sexes. In its traditional form, this game has two equivalent stable solutions (Nash equilibria), and the players face a dilemma in choosing between them. One conclusion of the article was that the dilemma is removed in the quantum game, i.e. it has a unique solution. In this Comment we make two observations. Firstly, the overall quantization scheme is fundamentally very similar to a previous scheme proposed by Eisert et al  – the similarity is non-obvious because of the very different use of the word ‘strategy’ in the two approaches. Secondly, we argue that the quantum Battle of the Sexes game does not in fact have a unique solution, hence the players are still subject to a dilemma, although it may be easier to resolve than in the traditional game.

The proposed scheme for quantizing two-player static games is shown in Fig 1a(i). Marinatto and Weber use the word ‘strategy’ to refer to a state $`S`$ in a $`22`$ dimensional Hilbert space. This ‘strategy’ is actually supplied to the players initially, they then manipulate it during a ‘tactics’ phase, and finally it is measured. The result of the measurement directly determines the payoff to the players , as shown in Fig. 1a(ii). In Ref.  the labels $`|O`$ and $`|T`$ are used for the measurement basis, whilst in Ref.  it is $`|C`$ and $`|D`$; here we will use the generic labels $`|0`$ and $`|1`$. During the ‘tactics’ phase of Marinatto and Weber’s scheme, each player may only manipulate $`S`$ within her own $`2`$ dimensional subspace, and the measurement basis respects this division. Hence it is useful to think of $`S`$ as being the global state of a pair of two-state quantum systems or ‘qubits’: each qubit is manipulated by one player during the ‘tactics’ phase, and then they are independently measured. Noting that we can generate any pure initial ‘strategy’ state $`S`$ by acting on the zero state $`|00`$ with an appropriate unitary operator $`\widehat{J}`$, we see the similarity to the earlier scheme shown in Fig 1b. Here we again prepare a 2-qubit system, allow players to make local manipulations of the qubits, and finally measure the qubits to determine payoffs just as in Fig 1a(ii). The only fundamental difference is that Eisert et al employ a second, inverse gate $`J^{}`$. The second gate ensures that the structure of the classical game remains embedded within the quantum game – this is elaborated upon in Ref.. However, the terminology used by Eisert et al is completely different – the manipulations applied by the players are referred to as ‘moves’ (rather than ‘tactics’), and moreover the term ‘strategy’ refers to a player’s overall act of choosing what move to play. These terms seem to be more consistent with traditional game theory, however this may be a matter of personal preference – the important point is that the underlying structures of the two schemes are very similar.

We now turn to the question of whether the players of the quantum game can escape from the dilemma suffered by their traditional counterparts. Marinatto and Weber claim that this occurs when the initial strategy is set to $`\frac{1}{\sqrt{2}}(|00+|11)`$, and the player’s ‘tactics’ are limited to a probabilistic choice between applying the identity, $`\widehat{I}`$, and applying $`\widehat{C}=\widehat{\sigma }_x`$, the Pauli spin-flip operator (this is a severe restriction on the full range of quantum mechanically possible manipulations). The symbols $`p^{}`$ and $`q^{}`$ are used to represent the probability that $`\widehat{\sigma }_x`$ will be chosen by players 1 and 2, respectively. The only ‘tactics’ consistent with maximum expected payoff are $`(p^{}=0,q^{}=0)`$, i.e. both players choosing to apply $`\widehat{I}`$ with certainty, or $`(p^{}=1,q^{}=1)`$, i.e. both choosing to apply $`\widehat{\sigma }_x`$. In either case the net effect is that the final ‘strategy’ is unchanged from the initial form, so that the measurement will yield either $`|00`$ or $`|11`$ with equal probability. Hence the expected payoff is $`(\alpha +\beta )/2`$ to each player. Now, since both $`(p^{}=0,q^{}=0)`$ and $`(p^{}=1,q^{}=1)`$ lead to the same final ‘strategy’, Marinatto and Weber claim that this is therefore the unique solution in the game. However, it seems to us that a clear dilemma remains for the players, even given that they have performed the above reasoning – should they opt for ($`p^{}=0`$,$`q^{}=0`$), or for $`(p^{}=1,q^{}=1)`$? Without knowing what the other player will choose, there is no way to make this decision! If the tactics are mismatched, i.e. if either $`(p^{}=0,q^{}=1)`$ or $`(p^{}=1,q^{}=0)`$ are adopted, then the worst-case payoffs occur. This is almost exactly the same problem faced by the players in the traditional game. However there is an important difference. In the traditional game, the two solutions have different payoffs to the players: $`(\alpha ,\beta )`$ or $`(\beta ,\alpha )`$. Thus each player prefers a different solution, and even if the players are allowed to communicate, there is no formal mechanism to decide who should get the high payoff. In the quantum game, since both solutions payout $`(\alpha +\beta )/2`$ to each player, there is no such difficulty, and permitting the players to communicate immediately allows them to agree on which solution to play for.

It is worth remarking on what would occur if the ‘tactics’ employed by Marinatto and Weber were less restricted. The following remarks apply when the initial ‘strategy’ supplied to the players is $`\frac{1}{\sqrt{2}}(|00+|11)`$. If the players were allowed to apply any local unitary operation on their qubit, i.e. any SU(2) matrix, then the maximum possible payoff to any player would remain $`(\alpha +\beta )/2`$. Moreover, for any given tactic $`\widehat{A}`$ applied by player 1, there is a corresponding tactic for player 2 that will result no net change to the ‘strategy’, and hence achieve maximum payoff to both players. The enquired tactic is simply the conjugate of $`\widehat{A}`$, i.e. $`\widehat{B}=\widehat{A}^{}`$. The problem of knowing which tactic has been adopted by the other player remains, just as discussed above. More interestingly, if the permitted tactics are extended to include all quantum mechanically possible manipulations, including measurement, entanglement which ancilla qubits, etc., then the quantum game can become more like the classical game! This is because either player can now choose to destroy the entanglement in the ‘strategy’ state, by simply measuring her qubit, and optionally flipping it. Moreover the players have a motivation to do this since the maximum possible payoff is then increased from $`(\alpha +\beta )/2`$ to $`\alpha `$. Therefore, with this most general kind of quantum tactic, even the limited advantage mentioned at the end of the last paragraph disappears, and the dilemma in the quantum game seems just as difficult as that in the traditional version.

The author is supported by EPSRC.


text:
# Mid-infrared imaging of the massive young star AFGL 2591: Probing the circumstellar environment of an outflow source

## 1 Introduction

It is now generally agreed that most, if not all, stars produce energetic outflows during their formation (see Bachiller 1996 for a review). Since the discovery of outflows almost 20 years ago (Snell, Loren & Plambeck 1980), nearly 200 examples have been found, most of them associated with low-luminosity young stellar objects (YSOs). Outflows can be studied using a variety of tracers including high-velocity CO line emission (Bally & Lada 1983), shock-excited H2 line emission (e.g., Lane 1989), and optically visible Herbig-Haro (HH) jets (Reipurth 1989, 1999).

All indications are that flows emerge bipolarly from a stellar or circumstellar region. However, the physical processes giving rise to the outflows are not well understood. One suggestion is that a primary wind originating at the star (e.g., Shu et al. 1991) or at the accretion disk (e.g., Pudritz & Norman 1983) is responsible. In both scenarios, the fast well-collimated wind sweeps up ambient molecular gas in its vicinity, forming two cavities in opposite directions from the star. The displaced molecular gas constitutes the CO outflow and ionized “knots” of gas manifest themselves as HH objects. Outflows probably occur while the star is still accreting material, and may play a significant role in ejecting material of excess angular momentum.

The presence of “bullets” in the CO emission and shocked “knots” visible in the optical and near-infrared suggests that mass ejection from YSOs often occurs in bursts rather than at a steady rate. Such ejection events appear to be intermittent, and the interval between two successive events is on the order of $`10^3`$ years. It is intriguing that the masses and timescales of bullets are similar to those of outbursts observed in FU Ori stars. The FU Ori eruptions are thought to be due to large increases in the accretion rate through a circumstellar disk (Hartmann & Kenyon 1985). If the two phenomena are related, it may be possible to probe the accretion history of protostars using detailed observations of their outflows (e.g., Reipurth 1989).

One of the barriers to better understanding the outflow process is the limited spatial resolution of millimeter-wavelength CO line observations. The situation is particularly dire in the case of massive young stellar objects in which large extinction ($`A_v`$ 10 mag) prevents high-resolution optical imaging.

Here we present sub-arcsecond-resolution mid-infrared images of one high-luminosity (Class 0/I) YSO, AFGL 2591, and its immediate surroundings. AFGL 2591 was first recognized as an outflow source by Bally & Lada (1983). The CO outflow is larger than 1’ in extent, with the blue lobe extending west or southwest and the red lobe in the northeast direction (Lada et al. 1984; Mitchell, Hasegawa, & Schella 1992; Hasegawa & Mitchell 1995). Optical and near-infrared observations have revealed molecular hydrogen emission (40” east and west, Tamura & Yamashita 1992) and HH objects (10” and 20” west/southwest, Poetzel et al. 1992) near AFGL 2591. The distance to AFGL 2591 is not well established. Assuming a distance of 1 kpc, it has a luminosity of $`2\times 10^4L_{\mathrm{}}`$ and an estimated stellar mass of 10 $`M_{\mathrm{}}`$.

## 2 Observations and data reduction

AFGL 2591 was observed on June 3 and 4, 1999 using MIRAC3, the newest version of the University of Arizona/SAO mid-infrared camera (Hoffmann et al. 1998) at the 3.0-meter NASA Infrared Telescope Facility (IRTF). MIRAC3 uses a Boeing HF16 128$`\times `$128 Si:As blocked impurity band detector. On the IRTF, MIRAC3 has a plate scale of 0.33”/pixel, providing a total field of view of 42”$`\times `$42”. This pixel scale ensures Nyquist sampling of the diffraction-limited point spread function (PSF).

We obtained images of AFGL 2591 in MIRAC3 11.7, 12.5 and 18.0 $`\mu `$m 10% passband filters, with a total on-source integration time of 600, 370 and 540 seconds respectively. The standard stars $`\alpha `$ Her and $`\beta `$ Peg, observed before and after the science target, were used for flux and PSF calibration.

We used a standard nodding and chopping technique to remove the background signal, dithering the source on the array to obtain sub-pixel sampling of the PSF. The chop frequency was set to 5 Hz, with a throw of 20” in the N-S direction. The nod throw was also set to 20” but in the E-W direction, in order to have all four chop-nod beams inside the field of view of the array. Each individual nod cycle required a 10-second on-source integration, and the procedure was repeated for as many cycles as needed to obtain the requested total integration time.

To reduce the data, we have developed our own software written in C and IDL. We first subtracted the chop-on from the chop-off frames for both nodding beams. The two images thus obtained were then substracted one with the other, in order to get a single frame in which the source appears in all four beams (two negative and two positive). We then applied a gain matrix, derived from images of the dome (high intensity, uniform background) and the sky (low intensity, uniform background), to flat field the chop-nodded image.

This procedure was repeated for each of the nodding cycles for which the source was observed. A final high S/N cumulative image was then obtained coadding together all beams, each recentered and shifted on the source centroid. This last coadding was performed on a sub-pixel grid having the size of one-fifth of the original MIRAC3 pixels, thus providing a final pixel scale of 0.066”/pixel. A mask file to block out the effects of bad pixels and field vignetting was also created and applied in this stage, preventing individual rejected pixels from contributing to the final image.

The same observing and reduction procedure was also used for the reference stars, to ensure a uniform treatment of the source and the standards.

## 3 Results and Discussion

Figure 1 shows the final MIRAC3 images of AFGL 2591. The central source is not clearly resolved. It has FWHM of 0.84”, 0.87”, and 1.21” at 11.7, 12.5 and 18.0 $`\mu `$m respectively. The FWHM of the reference point source $`\beta `$ Peg in the same filters is 0.82”, 0.85” and 1.11” respectively. An extended knot of emission $``$6” SW from the central source is visible in all three images. A compact faint source is seen $``$11” NW of AFGL 2591. Figure 2 is a schematic drawing of AFGL 2591, not to scale, with the various observed components indicated.

In Table 1, we list mid-infrared fluxes for the three sources as well as the position angle of SW and NW sources (counter-clockwise from North) and their separation from the central source. The photometry was measured inside apertures of (diameter) 18” (total flux of main and SW sources), 9” (SW source alone) and 5” (NW source). The zero magnitude flux densities for the three filters are 29.2, 25.7 and 12.6 Jy respectively. We estimate an error of 10% for the photometry, mainly due to the uncertainty of the reference fluxes and the difficulty of isolating the sources when doing aperture photometry. Lower fluxes of the main source were found by Lada et al. (1984) by using broad-band filters at 11.4 and 12.6 $`\mu `$m. The discrepancy (40% and 8% in the two filters) can be explained by the presence of a broad absorption feature in the source spectrum at 10 $`\mu `$m, which affects broad-band photometry more than our 10% filters. By fitting the 12.5 and 18.0 $`\mu `$m fluxes with a black body, we derive color temperatures of 300-350K for the central source, and 100-120K for the SW knot and the NW source.

Persi et al. (1995) presented the only previously published mid-infrared image of AFGL 2591. Their observations at 11.7 $`\mu `$m were carried out on the 2.1-meter S. Pedro Martir telescope in Baja California, Mexico with an image scale of 1.34”/pixel, and did not reveal the SW knot or the NW compact source seen in our MIRAC3 images, due to the smaller field of view of their detector.

### 3.1 Central source

The flux of the central source rises dramatically from 11.7 to 12.5 $`\mu `$m, implying that its 11.7 $`\mu `$m flux may be suppressed by a large silicate absorption feature at $``$10$`\mu `$m. Indeed, the low-resolution mid-infrared spectrum obtained by the InfraRed Astronomy Satellite (IRAS) shows deep absorption at that wavelength. Thus, it is likely that the protostar is surrounded by an optically thick dusty envelope, perhaps in free-fall collapse onto the star (van der Tak et al. 1999). Our observations limit the radius of the envelope to $``$3” (or $``$ 3000 AU, assuming a distance of 1 kpc).

Figure 3 shows the IRAS spectrum of the source, and the best fit we obtain by modeling the circumstellar dust envelope with the radiative transfer code DUSTY (Ivezić et al, 1999). The fit shown is for a spherically symmetric envelope with a radial density profile $`\rho r^2`$ around a central black body source at effective temperature $`T_{eff}`$ = 25,000 K (van der Tek et al. 1999). We adopted the silicate dust opacity by Ossenkopf et al. (1992), which provides a good fit to the 10 $`\mu `$m silicate feature. The fit for $`\rho r^{1.5}`$ (i.e., constant infall) is not as good, implying that the infalling envelope is likely to be surrounded by a fair amount of non-infalling cloud material. The model spectrum is not particularly sensitive to the $`T_{eff}`$ of the central source. Our best fit to the IRAS spectrum suggests an $`A_V`$ on the order of 100 (derived from the best fit optical depth).

### 3.2 NW source

AFGL 2591 is generally considered to be an example of relatively isolated massive star formation. However, our detection of a faint point-like source 11” NW of it suggests there may be other embedded low-mass stars in its vicinity. A preliminary look at Two-Micron All-Sky Survey (2MASS) J,H,K images reveals at least two more point-like sources within 40” of AFGL 2591. Deep mid-infrared images of a larger area are needed to search for a possible cluster of low-mass protostars in this region.

### 3.3 SW knot

The SW knot we detect appears close to the inner portion of a loop first observed in the near-infrared by Forrest & Shure (1986). The near-infrared emission is assumed to arise from scattering of photons from AFGL 2591 off dust grains at the surface of a cavity cleared by the outflow. The loop, also seen in NH3 observations by Torrelles et al. (1989), is indeed located in the blue outflow lobe, i.e., the one directed toward us. Using Draine & Lee (1984) opacity at 18$`\mu `$m and assuming $`T=`$100K and $`d=`$1 kpc, we estimate the mass of this knot to be $`>`$0.2$`M_{\mathrm{}}`$.

The SW knot is also coincident with a radio source as well as a knot of H2 emission. The radio source, detected using interferometers at the Very Large Array (Campbell 1984) and the Owens Valley Radio Observatory (van der Tak et al. 1999), has an approximately flat spectrum between 6.1 and 0.3 cm, and is interpreted as free-free emission from an optically thin HII region. The H2 emission is seen in an image taken with the BEAR imaging spectrometer on the Canada-France-Hawaii telescope using a narrow band filter on the 1-0 S(1) line, after the continuum is subtracted (G. Mitchell, personal communication). The existence of the H2 line implies shock collisional heating of the molecular gas in the outflow, which provide the necessary excitation energy ($`T2000`$ K, see e.g. Smith & Brand, 1990). On the other hand, the dust responsible for mid-infrared emission is characterized by a much lower local equilibrium temperature of $`T_d100`$ K, as inferred from our photometry. It is likely that the H2 emission comes from the shock itself, while the mid-infrared color temperature is an average through the dust clump.

Previous near-infrared observations have also revealed a second, larger loop of emission at a projected separation of $``$19” from the star (Burns et al. 1989; Minchin et al. 1991; Tamura et al. 1991; Poetzel, Mundt, & Ray 1992). The presence of multiple shells/loops is thought to be evidence for episodic mass outflows from AFGL 2591, as are molecular bullets seen in CO emission (Bachiller et al. 1991) and shocked knots observed in the optical and the near-infrared (Reipurth 1989) in other well-known outflow sources. Up to $`1M_{\mathrm{}}`$ could be ejected in each eruptive event, though a few tenths of a solar mass is more typical (Bachiller 1996).

If we assume that the SW knot is associated with a recent mass ejection from AFGL 2591, we can estimate its dynamical timescale. Following Minchin et al. (1991), we take that the outflow is inclined from the plane of the sky by 55 degrees, and estimate $`10^4`$ AU as the true distance from the star to the SW knot. Given a true outflow velocity of 26 kms-1 (based on the line-of-sight velocity of CO gas observed by Bally & Lada 1983), then we get a dynamical timescale of $``$2000 years. For the same assumptions, the outer loop seen in the near-infrared would have a dynamical timescale of $``$6000 years. Such timescales are consistent with the estimated $`10^4`$$`10^5`$-year total duration of the energetic outflow phase of YSOs (Bally & Lada 1983).

Another possibility is that the SW source itself is an embedded protostar surrounded by an ultracompact HII region. Some Herbig Ae/Be stars –LkH$`\alpha `$ 198 and LkH$`\alpha `$ 234, for example– are known to have heavily embedded mid-infrared companions associated with an outflow, a radio source, and HH objects (Lagage et al. 1993; Cabrit et al. 1997). The SW source is well-resolved in our mid-infrared images, with a ”coma”-like shape. So, if it is another protostar, its envelope may be affected by the outflow from AFGL 2591.

In summary, our observations suggest that mid-infrared imaging and spectroscopy can be powerful tools for probing the immediate circumstellar environment of embedded stars driving outflows. Such outflows appear to be rather common around young stars and are very likely associated with accretion processes from the early stages of protostars. They may play an important role in determining the mass of the star-disk system, and affect the evolution of the surrounding molecular cloud. Yet our understanding of the outflow process is far from satisfactory. High-resolution imaging with current mid-infrared cameras and future millimeter-wave interferometers could help develop a more unified picture of infall and outflow in star formation.

We are grateful to George Mitchell, Charles Lada and Paolo Persi for helpful discussions and to an anonymous referee for sensible suggestions. We thank the IRTF staff for their outstanding support. The research was supported by the Smithsonian Institution and NASA grants to the SIRTF-IRAC project.

Figure Captions

Figure 1. MIRAC3 images of AFGL 2591 at (a) 11.7$`\mu `$m, (b) 12.5$`\mu `$m, and (c) 18.0$`\mu `$m. Contours are plotted for 10 and 15 $`\sigma `$ levels in panels (a) and (b), and for 5, 10 and 15 $`\sigma `$ in panel (c).

Figure 2. A schematic drawing of AFGL 2591 with the various observed components indicated. (Not to scale.)

Figure 3. IRAS mid-infrared spectrum of AFGL 2591 (solid line) and the best-fit model (dashed line) for the 10$`\mu `$m silicate feature. The heavy absorption is presumably due to an optically thick dusty envelope around the massive central star.


text:
# References

Hybrid classical-quantum dynamics

Asher Peres and Daniel R. Terno

Department of Physics, Technion—Israel Institute of Technology, 32 000 Haifa, Israel

Abstract

A hybrid formalism is proposed for interacting classical and quantum sytems. This formalism is mathematically consistent and reduces to standard classical and quantum mechanics in the case of no interaction. However, in the presence of interaction, the correspondence principle is violated.

PACS numbers: 03.65.Bz, 03.65.Sq

Electronic mail: peres@photon.technion.ac.il, terno@physics.technion.ac.il

Quantum mechanics gives exceedingly accurate predictions for atomic and nuclear systems. Classical mechanics is just as successful for planetary motion. Can there be a theory encompassing both classical and quantum mechanics, so that atoms and planets can be treated in a unified way? Indeed, it was speculated long ago that the gravitational field is classical, even though its material sources are quantized . Niels Bohr always insisted that measuring instruments must be described in classical terms . However, Bohr did not provide a dynamical description of their interaction with quantum systems. A dual (classical and quantum) description of measuring instruments is possible , but this is not a truly hybrid formalism as proposed below, because in that dual formalism the interaction is always between subsystems of the same kind.

If we try to write equations of motion that combine classical canonical variables and quantum operators, we find that the classical variables, whose values ought to be ordinary numbers, pick up operator components, which is manifestly unacceptable (we shall later show how to overcome this difficulty). There have been numerous claims and counter-claims as to the possibility of existence a consistent classical-quantum formalism. Some recent articles are listed in the references . There is no real contradiction between them, because their conflicting conclusions result from different demands for consistency.

In this article, we propose a hybrid formalism for interacting classical and quantum sytems. This formalism is mathematically consistent and reduces to standard classical and quantum mechanics if there is no interaction between the systems. However, in the presence of interaction, the correspondence principle is violated, so that this hybrid formalism appears physically abnormal.

Our work follows a remark by Jauch  that if we assume that not all quantum dynamical variables are actually observable, and if we set rules for distinguishing measurable from nonmeasurable operators, it is then possible to define a classical system as a special type of quantum system for which all measurable operators commute. All the observable dynamical variables of such a classical system can be simultaneously diagonalized and therefore can, in principle, have sharp values. If we are confined to this set of observables and we use a basis where they all are diagonal, then only diagonal elements of density matrices can be reconstructed from measurement statistics and have an operational meaning. Quantum states (i.e., density matrices) with the same diagonal elements are indistinguishable. We may therefore consider all states as pure by defining their off-diagonal elements as $`\rho _{mn}=\sqrt{\rho _{mm}\rho _{nn}}`$, so that $`\rho `$ always has rank 1. It is then possible to define a classical state vector $`\psi `$ by $`\rho =\psi \psi ^{}`$.

Pure states defined in this way do not imply that dynamical variables have sharp values. The situation is just as in classical statistical mechanics: the expectation value of an observable $`A`$ is given by

$$A=\text{tr}(\rho A),$$
(1)
and we may then have $`A^2>A^2`$, as usual. Note that such a “classical” system also has noncommuting operators, but the latter should be considered as abstract mathematical expressions which are not experimentally observable. With this formal definition, the notion “classical” acquires a meaning with respect to a specified set of dynamical variables. This does not preclude that, by using more elaborate observation means, every physical system might display quantum features.

To completely mimic classical mechanics, we need a quantum algorithm that reproduces exactly the equations of motion of a classical canonical system. For this, we shall follow Koopman’s formalism . We consider for simplicity a single degree of freedom and denote the canonical variables as $`x`$ and $`y`$ (rather than $`x`$ and $`p`$, as usual, because we wish to reserve the symbol $`p`$ for the momentum of a quantum system, to be introduced later). Let us write the Liouville equation as

$$if/t=Lf,$$
(2)
where $`L`$ is the Liouville operator, or Liouvillian,

$$L=\left(\frac{H}{y}\right)\left(i\frac{}{x}\right)\left(\frac{H}{x}\right)\left(i\frac{}{y}\right).$$
(3)
The Liouville density $`f`$ is never negative. Since its quantum analog, namely Wigner’s distribution  (which may be negative), is quadratic in the quantum wave function, it is convenient to introduce likewise a “classical wave function” $`\psi =\sqrt{f}`$, which in this case satisfies the same equation of motion as $`f`$,

$$i\psi /t=L\psi .$$
(4)
We shall now consider $`\psi `$ as the fundamental entity (but only $`f=|\psi |^2`$ has a direct physical meaning). It can be proved that, under reasonable assumptions about the Hamiltonian, the Liouvillian is an essentially self-adjoint operator and generates a unitary evolution :

$$\psi ,\varphi :=\psi (x,y,t)^{}\varphi (x,y,t)𝑑x𝑑y=\mathrm{const}.$$
(5)
The introduction of a Hilbert space with unitary dynamics enables one to use familiar methods of quantum mechanics for the analysis of ergodic problems  and of classical chaos . It is possible to further mimic quantum theory by introducing commuting operators $`\widehat{x}`$ and $`\widehat{y}`$, defined by

$$\widehat{x}\psi =x\psi (x,y,t)\mathrm{and}\widehat{y}\psi =y\psi (x,y,t).$$
(6)
Note that the momentum $`\widehat{y}`$ is not the shift operator (the latter is $`\widehat{p}_x=id/dx`$). Likewise the boost operator is $`\widehat{p}_y=id/dy`$. These two operators are not observable. We shall henceforth omit the hats over the classical operators, as there is no danger of confusion.

The analogy with quantum mechanics can be pushed further. What we have above is a “Schrödinger picture” (operators are constant, wave functions evolve in time as $`\psi (t)=U(t)\psi (0)`$, where $`U(t)=e^{iLt}`$ if the Hamiltonian is time-independent). We can also define a “Heisenberg picture” where wave functions are fixed and operators evolve:

$$X_H(t)=U^{}XU.$$
(7)
The Heisenberg equation of motion,

$$idX_H/dt=[X_H,L_H]=U^{}[X,L]U,$$
(8)
together with the Liouvillian (3), readily give Hamilton’s equations

$$\frac{dx}{dt}=\frac{H}{y},\frac{dy}{dt}=\frac{H}{x}.$$
(9)
There is however an important difference: the time translation operator $`L`$ is not the energy, and its spectrum may extend to $`\mathrm{}`$. For example, if we have a harmonic oscillator with $`H=(x^2+y^2)/2`$, the Liouvillian is

$$L=yp_xxp_y,$$
(10)
whose eigenvalues are all the integers, negative as well as positive. There is nothing wrong in that, since $`L`$ involves the unobservable shift operators $`p_x=id/dx`$ and $`p_y=id/dy`$, and therefore $`L`$ itself is not observable. Note that the solution of Hamilton’s equations does not introduce non-observable components into the observable variables, because the Hamiltonian (contrary to the Liouvillian) involves only the observables $`x`$ and $`y`$.

It is easy to introduce into the above dynamical formalism a quantum system, with conjugate dynamical variables $`q`$ and $`p`$, as long as the two systems do not interact. In the Schrödinger picture, we have a wave function $`\psi (x,y,q)`$ whose evolution is given by

$$\psi /t=K\psi ,$$
(11)
where $`K`$ is “Koopmanian”

$$K=L(x,y,p_x,p_y)+H(q,p).$$
(12)
In that wave function, the coordinates $`x`$, $`y`$, and $`q`$ may be entangled. In any case, we can obtain the reduced density matrices of the “classical” and quantum systems by means of partial traces on the other system, as usual.

The non-trivial problem is to introduce an interaction between the two systems. An apparent difficulty occurs in the Heisenberg picture because the equations of motion will mix variables of all kinds: classical observables, non-observables, and quantum operators. This still is acceptable, because we do not predict actual values for these variables, but only expectation values, $`A=\text{tr}(\rho A)`$, and no contradiction may occur. Anyway, we can always use the Schrödinger picture, where we know how to handle entangled wave functions.

The true difficulty, as we shall now see in a simple example, is that the correspondence principle fails. Consider two harmonic oscillators, with a bilinear coupling $`kqx`$, where $`k`$ is a constant. If we treat both of them classically, with a Hamiltonian

$$H=\frac{1}{2}(q^2+p^2+x^2+y^2)+kqx,$$
(13)
we obtain equations of motion

$$\dot{q}=p,\dot{p}=qkx,$$
(14)
$$\dot{x}=y,\dot{y}=xkq,$$
(15)
and there are two characteristic frequencies, $`\omega =\sqrt{1\pm k}`$, corresponding to the normal modes $`(q\pm x)`$. Exactly the same equations of motion appear in the Heisenberg picture for quantum mechanics. Likewise, in the Schrödinger picture, whenever there is a quadratic potential, the differential equation for the Wigner distribution is identical to the Liouville equation in classical mechanics . In view of this formal agreement of classical and quantum mechanics, it is natural to demand that any hybridization of the system in Eq. (13) shall also give the same equations of motion. This is the definite benchmark we propose for an acceptable classical-quantum hybrid formalism. It is not obvious that this criterion can be achieved, because the hybrid system has no normal modes $`(q\pm x)`$.

Let us try to obtain the above equations of motion from a Koopmanian

$$K=\frac{1}{2}(q^2+p^2)+(yp_xxp_y)+K_i,$$
(16)
with a suitably chosen interaction term $`K_i`$. We cannot have both $`[p,K_i]=kx`$ and $`[y,K_i]=kq`$ together with $`[y,p]=0`$, because these equations are incompatible with Jacobi’s identity

$$[y,[p,K_i]]+[p,[K_i,y]]+[K_i,[y,p]]0.$$
(17)
The best result we were able to obtain was by “koopmanizing” the interaction term $`kqx`$ in the Hamiltonian (13). We have, from Eq. (3),

$$K_i=ikq/ykqp_y.$$
(18)
The resulting equations of motion are Eq. (15), unchanged, and

$$\dot{q}=p,\dot{p}=qkp_y,$$
(19)
$$\dot{p}_x=p_y,\dot{p}_y=p_x.$$
(20)
The last equation is necessary, because an unobservable variable $`p_y`$ appears in Eq. (19). The solution of (20) for $`p_y`$ is a superposition of $`\mathrm{sin}t`$ and $`\mathrm{cos}t`$ terms. When substituted into Eq. (19), these terms behave as a driving force with resonant frequency, so that $`q`$ and $`p`$ include terms behaving as $`t\mathrm{sin}t`$ and $`t\mathrm{cos}t`$. The amplitude of the quantum oscillator increases linearly with time, and energy is not conserved. It need not be: what is conserved is the Koopmanian, which includes a term $`(yp_xxp_y)`$ whose spectrum extends to $`\mathrm{}`$.

In conclusion, there is no mathematical inconsistency in the hybrid formalism that we proposed. However, it violates the correspondence principle which we would expect to hold exactly for a pair of oscillators with bilinear coupling. Therefore such a theory appears quite abnormal from the point of view of physics.

We thank Lajos Diósi for stimulating discussions. DRT was supported by a grant from the Technion Graduate School. Work by AP was supported by the Gerard Swope Fund and the Fund for Encouragement of Research.


text:
# Wave packet evolution approach to ionization of hydrogen molecular ion by fast electrons

## I Introduction

New experimental methods, particularly, based on the multiple coincidence detection technique  stimulate the interest to fundamental theoretical studies of the dissociative ionization of diatomic molecules by electron impact. In this context the molecular hydrogen ion can be considered as the basic system in which the removal of the unique electron causes dissociation. Substantial theoretical analysis of the dissociative ionization of $`\mathrm{H}_2^+`$ by fast electrons was recently carried out in . As mentioned in , the crucial point of calculating the cross-section of such processes is that no closed exact analytical wave functions of the continuum states exist. In  the final-state wave function of the ejected electron was found by taking a product of two approximate functions that take into account the two scattering centers. To improve the calculation it seems straightforward to obtain these functions with the exact numerical solutions of the two-center continuum problem. However, this approach involves a cumbersome procedure of calculating multi-dimensional integrals of the functions presented numerically that requires huge computer facilities and may cause additional computational problems. It seems reasonable to search for direct computational approaches, in which the basis of exact two-center continuum wave functions is not involved. Note that the potential advantage of such methods is that they could be generalized over a wider class of two-center systems starting from the molecular hydrogen ion as a test object. In the present paper we develop a direct approach to the ionization of hydrogen molecular ion by fast electrons that involves the reduction of the initial 6D Schrödinger equation to a 3D evolution problem followed by modeling of the wave packet dynamics.

Originally we intended to treat the incoming electron classically, its trajectory being approximated by a straight line with the deflection neglected. The bound electron was to be treated quantum mechanically. Preliminary calculations at the impact parameter $`\rho =10`$ a.u. has shown, first, that the probability of the emission of the electron having the energy of 50 eV is extremely small, and, second, that the direction of the electron emission is orthogonal to that of the incoming electron motion, that contradicts the results of . This means that the main contribution to the small-angle scattering comes from the central collisions with the bound electron in the region of its localization. Generally, the classically estimated deflection angle of $`1^o`$ for scattered electron corresponds to the impact parameter of the order of 1 a.u., so that the trajectory passes through the molecule and the classical treatment of the incoming electron is not valid.

Here we develop and apply a direct approach to the calculation of the angular distribution of scattered and ejected electrons that involves the reduction of the initial 6D Schrödinger equation to a 3D evolution problem followed by modeling of the wave packet dynamics. The approach does not make use of the basis of stationary Coulomb two-center functions of the continuous spectrum for the ejected electron, whose proper choice is a crucial point of other model calculations. Our approach can be considered as the linearized version of the phase function method  for the multi-dimensional scattering problem. The evolution problem is solved using the method based on the split-step technique  with complex scaling, recently proposed by some of us and tested in paraxial optics . In the present paper the method as a whole is also tested using the well known problem of electron scattering by hydrogen atom .

## II Basic equations

We start from the 6D stationary Schrödinger equation which describes two electrons in the field of two fixed protons

$$\left[H_0(𝐫)\frac{1}{2}_𝐑^2+V(𝐫,𝐑)\right]\mathrm{\Psi }(𝐫,𝐑)=E\mathrm{\Psi }(𝐫,𝐑),$$
(1)
where $`𝐫`$ is the radius-vector of the electron initially bound in $`\mathrm{H}_2^+`$ and finally ejected, $`𝐑`$ is the radius-vector of the impact electron, $`\widehat{H}_0=\frac{1}{2}_𝐫^2+U(𝐫)`$ is Hamiltonian of ejected electron in the field of two protons, $`V(𝐫,𝐑)=U(𝐑)+U_{int}(𝐫,𝐑)`$ is the interaction between the impact electron and molecular ion, $`U(𝐫)=1/r_11/r_2`$ is the attractive potential between the ejected (scattered) electron and the protons, $`r_1=|𝐫𝐫_{1p}|`$, $`r_2=|𝐫𝐫_{2p}|`$, $`𝐫_{ip}`$ is the radius-vector of the i-th proton, $`U_{int}(𝐫,𝐑)=1/|𝐫𝐑|`$ is the repulsive potential of interaction between the electrons. The origin of the coordinate frame is chosen in the center of symmetry of the molecular ion with the $`Z`$ axis directed along the momentum of the incident electron.

For the scattering problem solved here the energy of the system may be presented as $`E=k_i^2/2+E_0`$, where $`E_0`$ is the ionization potential, $`k_i`$ is the momentum of the incident electron. Let us seek the solution of Eq.(1) in the form $`\mathrm{\Psi }(𝐫,X,Y,Z)=\psi (𝐫,𝐑_{},Z)\mathrm{exp}(ik_iZ)`$. Under the condition that $`(k_e^2+k_{}^22E_0)/k_i^2<<1`$ one can neglect the second derivative of $`\psi `$ with respect to $`Z`$. As a result we get the evolution-like equation for the envelope function $`\psi (𝐫,𝐑_{},Z)`$

$`ik_i{\displaystyle \frac{\psi (𝐫,𝐑_{},Z)}{Z}}=`$ (3)
$`\left\{\widehat{H}_0(𝐫){\displaystyle \frac{1}{2}}_𝐑_{}^2E_0+V(𝐫,𝐑)\right\}\psi (𝐫,𝐑_{},Z).`$
Neglecting the large-angle scattering one can write the initial condition for $`\psi `$ as

$$\psi (𝐫,𝐑_{},\mathrm{})=\psi _0(𝐫).$$
(4)
To solve the 5D Schrödinger evolution equation(3) we use Fourier transformation with respect to the variable $`𝐑_{}`$

$$\psi (𝐫,𝐑_{},Z)=\frac{1}{2\pi }\psi _𝐤_{}(𝐫,Z)\mathrm{exp}(i𝐤_{}𝐑_{})𝑑𝐑_{}.$$
(5)
Then Eq.(3) takes the form

$`ik_i{\displaystyle \frac{\psi _𝐤_{}(𝐫,Z)}{Z}}=\left\{\widehat{H}_0(𝐫)+\left({\displaystyle \frac{k_{}^2}{2}}E_0\right)\right\}\psi _𝐤_{}(𝐫,Z)`$ (6)
$`+{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle V_{𝐤_{}𝐤_{}^{}}(𝐫,Z)\psi _𝐤_{}^{}(𝐫,Z)𝑑𝐤_{}^{}},`$ (7)
where

$$V_{𝐤_{}𝐤_{}^{}}(𝐫,Z)=\mathrm{exp}(i(𝐤_{}𝐤_{}^{})𝐑_{})V(𝐫,𝐑_{},Z)𝑑𝐑_{}$$
(8)
is the Fourier transform of the interaction potential $`V(𝐫,𝐑_{},Z)`$.

Further simplification of the problem is possible if the amplitude of the incident wave is much greater than that of the scattered wave. In this case one can put

$$\psi _𝐤_{}(𝐫,Z)=\delta (𝐤_{})\psi _0(𝐫)$$
(9)
in the integral term of Eq.(7). As a result we get the inhomogeneous equation

$`ik_i{\displaystyle \frac{\psi _𝐤_{}(𝐫,Z)}{Z}}`$ $`=`$ $`\left\{\widehat{H}_0(𝐫)+\left({\displaystyle \frac{k_{}^2}{2}}E_0\right)\right\}\psi _𝐤_{}(𝐫,Z)`$ (10)
$`+`$ $`{\displaystyle \frac{1}{(2\pi )^2}}V_𝐤_{}(𝐫,Z)\psi _0(𝐫),`$ (11)
where $`V_𝐤_{}(𝐫,Z)=V_{𝐤_{}\mathrm{𝟎}}(𝐫,Z)`$, with the initial condition $`\psi _𝐤_{}(𝐫,\mathrm{})=0`$.

To calculate the integral with respect to transverse variables in the expression for $`V_𝐤_{}(𝐫,Z)`$ it is easier to start from the known integral

$$\mathrm{exp}(i\mathrm{𝐤𝐑})\frac{1}{R}𝑑𝐑=\frac{4\pi }{k^2}=\frac{4\pi }{k_Z^2+k_{}^2}.$$
(12)
Carrying out the inverse Fourier transformation

$$_{\mathrm{}}^{\mathrm{}}\mathrm{exp}(ik_ZZ)\frac{dk_Z}{k_Z^2+k_{}^2}=\frac{\pi }{k_{}}e^{k_{}|Z|},$$
(13)
one gets

$`V_𝐤_{}(𝐫,Z)={\displaystyle \frac{2\pi }{k_{}}}e^{k_{}|Zz|i𝐤_{}𝐫_{}}`$ (15)
$`{\displaystyle \frac{2\pi }{k_{}}}\left[e^{k_{}|Zd_Z|i𝐤_{}𝐝_{}}+e^{k_{}|Z+d_Z|+i𝐤_{}𝐝_{}}\right].`$
Here $`k_{}=k_i\mathrm{sin}\theta _s`$ is the transverse momentum component of the scattered electron, $`\theta _s`$ is the scattering angle, $`\pm 𝐝`$ are the positions of the nuclei with respect to the center of symmetry. Note that terms in square brackets determine the elastic scattering of the incident electron by the nuclei.

Due to the exponential decrease of the source term with $`|Z|`$ the integration may be actually carried out within a certain finite interval $`(Z_{max},Z_{max})`$. Hence the zero initial condition should be imposed at the point $`Z_{max}`$.

Note that the approximation (9) is actually equivalent to the first Born approximation . Multiply Eq.(11) by the complex conjugate function of the continuous spectrum of $`\widehat{H}_0`$ and integrate over all $`𝐫`$. Then

$`ik_i{\displaystyle \frac{dC_𝐤_{}(𝐤_e,Z)}{dZ}}=\left\{{\displaystyle \frac{k_e^2}{2}}+{\displaystyle \frac{k_{}^2}{2}}E_0\right\}C_𝐤_{}(𝐤_e,Z)`$ (17)
$`+{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle \psi ^{}(𝐤_e,𝐫)V_𝐤_{}(𝐫,Z)\psi _0(𝐫)𝑑𝐫},`$
where $`C_𝐤_{}(𝐤_e,Z)=\psi ^{}(𝐤_e,𝐫)\psi _𝐤_{}(𝐫,Z)𝑑𝐫`$ is the probability density amplitude for the transition of the initially bound electron into the state with the momentum $`𝐤_e`$. Let us substitute

$$C_𝐤_{}(𝐤_e,Z)=\stackrel{~}{C}_𝐤_{}(𝐤_e,Z)\mathrm{exp}(ik_ZZ),$$

where $`k_Z`$ is the increment of the longitudinal component of the momentum of the impact electron determined by the relation

$$k_Z=\frac{1}{k_i}\left(\frac{k_e^2}{2}+\frac{k_{}^2}{2}E_0\right).$$
(18)
This relation is actually equivalent to the energy conservation law written neglecting the terms of the order of $`k_Z^2`$. The substitution yields

$`ik_i{\displaystyle \frac{d\stackrel{~}{C}_𝐤_{}(𝐤_e,Z)}{dZ}}=`$ (20)
$`{\displaystyle \frac{1}{(2\pi )^2}}\mathrm{exp}(ik_ZZ){\displaystyle \psi ^{}(𝐤_e,𝐫)V_𝐤_{}(𝐫,Z)\psi _0(𝐫)𝑑𝐫},`$
and

$`\stackrel{~}{C}_𝐤_{}(𝐤_e,\mathrm{})=`$ (22)
$`{\displaystyle \frac{1}{ik_i(2\pi )^2}}e^{i𝐤_s𝐑}\psi (𝐤_e,𝐫)|V(𝐫,𝐑)|e^{i𝐤_i𝐑}\psi _0(𝐫),`$
where $`𝐤_s=𝐤_i𝐊`$ is the momentum of the scattered electron, $`𝐊=(k_X,k_Y,k_Z)`$ is the momentum transfer.

Provided that the ejected electron has the momentum $`𝐤_e`$, the asymptotic form of the solution of Eq. (1) for the wave function of the scattered electron when $`R\mathrm{}`$ is

$$\mathrm{\Psi }_{𝐤_e}^{as}(𝐑)=\mathrm{exp}(ik_iZ)+\frac{\mathrm{exp}(ik_sR)}{R}f_{𝐤_e}(\theta _s,\varphi _s).$$
(23)
The scattering differential cross-section(DCS) can be then expressed as

$$\sigma _{𝐤_e}(\theta _s,\varphi _s)=\frac{k_ek_s}{k_i}\left|f_{𝐤_e}(\theta _s,\varphi _s)\right|^2,$$
(24)
On the other hand, the asymptotic form of the wave function resulting from the solution of Eq.(11) under the condition $`Z\mathrm{}`$ can be presented as

$`\mathrm{\Psi }_{𝐤_e}^{as}(𝐑)=\mathrm{exp}(ik_iZ)+`$ (26)
$`\mathrm{exp}(ik_iZ){\displaystyle \stackrel{~}{C}_𝐤_{}(𝐤_e,\mathrm{})\mathrm{exp}(i𝐤_{}𝐑_{}+ik_ZZ)𝑑𝐤_{}}.`$
Making use of the fact that the integrand has a stationary point we finally get

$`\mathrm{\Psi }_{𝐤_e}^{as}(𝐑)=e^{ik_iZ}+`$ (28)
$`{\displaystyle \frac{1}{Z}}e^{i\left(k_i\frac{k_e^2/2E_0}{k_i}\right)Z+i\frac{k_i}{2Z}R_{}^2}(2\pi ik_i)\stackrel{~}{C}_{𝐤_{}^0}(𝐤_e,\mathrm{}),`$
where $`𝐤_{}^0=k_i\mathrm{sin}\theta _s(\mathrm{cos}\varphi _s,\mathrm{sin}\varphi _s)`$, $`R_{}=R\mathrm{sin}\theta _s`$, $`Z=R\mathrm{cos}\theta _s`$. The expression (28) agrees with (23) within the accuracy of the order of $`\theta _s^2`$ if we set

$`f_{𝐤_e}(\theta _s,\varphi _s)=2\pi ik_i\stackrel{~}{C}_{𝐤_{}^0}(𝐤_e,\mathrm{})`$ (30)
$`={\displaystyle \frac{1}{2\pi }}e^{i𝐤_s𝐑}\psi (𝐤_e,𝐫)|V(𝐫,𝐑)|e^{i𝐤_i𝐑}\psi _0(𝐫).`$
The latter expression is similar to the formula for $`f_{𝐤_e}(\theta _s,\varphi _s)`$ derived in  using the first Born approximation.

## III Calculation of the angular distribution

The asymptotic expression of the radial part of the wave function corresponding to the continuous spectrum of $`\widehat{H}_0`$ can be written as

$$\psi _E^{as}(r,t)=\frac{1}{\sqrt{\upsilon (r)}r}\mathrm{exp}(iEt+i^r\upsilon (r^{})𝑑r^{})$$
(31)
where $`t=Z/k_i`$ is the evolution variable, $`\upsilon (r)=\sqrt{2(EU^{as}(r))}`$, $`E=\frac{k_e^2}{2}+\frac{k_{}^2}{2}E_0`$, $`U^{as}(r)=Z^{}/r`$, $`Z^{}=2`$ is charge of two protons. In the asymptotic limit one can take only the radial component of the momentum of the ejected electron into account. Then, according to , the expression for calculating the amplitude $`A(k,\theta ,\varphi )`$ takes the form

$`A_𝐤_{}(k_e,\theta _e,\varphi _e)=`$ (33)
$`{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle _{t_0}^{t_1}}𝑑t^{}j(\psi _𝐤_{}(r,\theta ,\varphi ,t^{}),\psi _E^{as}(r,t^{}))|_{r=r_{max}},`$
where

$$j(\mathrm{\Psi },\mathrm{\Phi })=\frac{i}{2}\left\{\mathrm{\Psi }r^2\frac{\mathrm{\Phi }^{}}{r}\mathrm{\Phi }^{}r^2\frac{\mathrm{\Psi }}{r}\right\},$$
(34)
is the flux introduced in , $`t_0=Z_{max}/k_i`$ and $`t_1>>Z_{max}/k_i`$. The approximate relation (33) becomes exact when $`t_1+\mathrm{}`$ and simultaneously $`r_{max}+\mathrm{}`$.

The amplitudes defined by (33) are related with the coefficients introduced in Eq.(22) by

$$\left|A_𝐤_{}(k_e,\theta _e,\varphi _e)\right|^2=k_e\left|\stackrel{~}{C}_𝐤_{}(𝐤_e,\mathrm{})\right|^2$$
(35)
Using (24),(30) and (35) we get the final expression for the differential cross-section

$$\sigma _{𝐤_e}(\theta _s,\varphi _s)=(2\pi )^2k_sk_i\left|A_𝐤_{}(k_e,\theta _e,\varphi _e)\right|^2.$$
(36)
In the region where $`r>r_{max}`$ we made use of the complex scaling technique  to suppress the non-physical reflection from the grid boundary.

## IV Numerical scheme

The inhomogeneous Schrödinger equation to be solved can be written as

$$i\frac{\mathrm{\Psi }(𝐫,t)}{t}=\widehat{H}_0(𝐫)\mathrm{\Psi }(𝐫,t)+F(𝐫,t),$$
(37)
The solution of Eq. (37) to within the second-order terms in $`\mathrm{\Delta }t`$ can be expressed as the following sequence of equations:

$`\mathrm{\Psi }_l^0`$ $`=`$ $`\mathrm{\Psi }(t){\displaystyle \frac{i\mathrm{\Delta }t}{2}}F(𝐫,t);`$ (38)
$`(1+i{\displaystyle \frac{\mathrm{\Delta }t}{2}}\widehat{H}_0)\mathrm{\Psi }_l^1`$ $`=`$ $`(1i{\displaystyle \frac{\mathrm{\Delta }t}{2}}\widehat{H}_0)\mathrm{\Psi }_l^0;`$ (39)
$`\mathrm{\Psi }(t+\mathrm{\Delta }t)`$ $`=`$ $`\mathrm{\Psi }_l^1{\displaystyle \frac{i\mathrm{\Delta }t}{2}}F(𝐫,t+\mathrm{\Delta }t).`$ (40)
The key step of the procedure is Eq. (39) which defines nothing but Cranck-Nicholson scheme. To solve this equation we make use of the partial coordinate splitting (PCS). A finite-difference scheme is applied for the radial variable $`r`$ and the polar angle $`\theta `$. Fast Fourier transform (FFT) is used for the azimuthal angle $`\varphi `$.

In the spherical coordinate system, the $`z`$-axis of which is directed along the symmetry axis of the molecule (and not along the velocity of the impact electron) and after the substitution $`\mathrm{\Psi }=\psi /r`$ and the Fourier transformation of Eq.(39) the terms $`\widehat{H}_0\mathrm{\Psi }_l^{0,1}`$ entering this equation turn into

$`\widehat{H}_0(r,\eta ,m)\psi (r,\eta ,m)=`$ (43)
$`{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{^2}{r^2}}+{\displaystyle \frac{1}{r^2}}\left({\displaystyle \frac{}{\eta }}(1\eta ^2){\displaystyle \frac{}{\eta }}{\displaystyle \frac{m^2}{1\eta ^2}}\right)\right]\psi (r,\eta ,m)`$
$`+U(r,\eta )\psi (r,\eta ,m),`$
where $`\eta =\mathrm{cos}\theta `$, $`m`$ is the asimuthal quantum number.

Finite-difference approximation $`\frac{^2}{r^2}D_{i_2i_1}^r`$ and $`\frac{}{\eta }(1\eta ^2)\frac{}{\eta }\frac{m^2}{1\eta ^2}D_{j_2j_1}^{\eta ,m}`$ of the differential operators entering Eq.(39) yields $`M`$ sets of linear equations, each set being of the order $`L\times N`$, where $`M,L`$ and $`N`$ are the numbers of grid points in $`\varphi `$, $`\eta `$ and $`r`$, respectively. Direct solution of each set of equations requires $`NL^2`$ operations at each step in $`t`$. The FFT that should be performed twice, first, when proceeding from (38) to (39), and second, from (39) to (40), requires $`NLM\mathrm{log}_2M`$ extra operations.

To reduce the number of operation we propose a double-cycle split-step scheme. In case when $`\widehat{H}_0`$ can be presented as a sum $`\widehat{H}_0=\widehat{H}_1+\widehat{H}_2`$, this scheme can be formulated as follows

$`\psi _1`$ $`=`$ $`\psi (t);`$ (44)
$`(I+i{\displaystyle \frac{\mathrm{\Delta }t}{4}}\widehat{H}_1)\psi _2`$ $`=`$ $`(Ii{\displaystyle \frac{\mathrm{\Delta }t}{4}}\widehat{H}_1)\psi _1;`$ (45)
$`(I+i{\displaystyle \frac{\mathrm{\Delta }t}{4}}\widehat{H}_2)\psi _3`$ $`=`$ $`(Ii{\displaystyle \frac{\mathrm{\Delta }t}{4}}\widehat{H}_2)\psi _2;`$ (46)
$`(I+i{\displaystyle \frac{\mathrm{\Delta }t}{4}}\widehat{H}_2)\psi _4`$ $`=`$ $`(Ii{\displaystyle \frac{\mathrm{\Delta }t}{4}}\widehat{H}_2)\psi _3;`$ (47)
$`(I+i{\displaystyle \frac{\mathrm{\Delta }t}{4}}\widehat{H}_1)\psi _5`$ $`=`$ $`(Ii{\displaystyle \frac{\mathrm{\Delta }t}{4}}\widehat{H}_1)\psi _4;`$ (48)
$`\psi (t+\mathrm{\Delta }t)`$ $`=`$ $`\psi _5,`$ (49)
which to within the second-order terms in $`\mathrm{\Delta }t`$ corresponds to the initial Cranck-Nicholson scheme, $`I`$ and $`H_{1,2}`$ is square matrixes, $`(I)_{i2i1j2j1}=\delta _{i2i1}\delta _{j2j1}`$.

Now the problem is how to split the Hamiltonian $`\widehat{H}_0`$ into two parts. Formal separation of radial and angular parts leads to difficulties associated with the singularity of the angular part. Due to this singularity the scheme appears to be conditionally stable with severe limitations imposed on the step $`\mathrm{\Delta }t`$. Practically this version of the splitting scheme is applicable only if the grid in $`r`$ is rough enough.

To remove this limitation we propose a partial coordinate splitting scheme. Its principal idea is that in the vicinity of $`r=0`$ it is preferable not to split off the angular part at all. To implement this idea we introduce the $`r`$-dependent weight function $`p(r)`$ which is supposed to diminish in the vicinity of $`r=0`$ and define the discrete approximation of the operators $`\widehat{H}_{1,2}`$ in the following way

$`(\widehat{H}_1^m)_{i_2i_1j_2j_1}={\displaystyle \frac{1}{2}}\widehat{D}_{}^{r}{}_{i_2i_1}{}^{}\delta _{j_2j_1}+U^{as}(r_i)`$ (51)
$`+p(r_{i_1})\left[{\displaystyle \frac{1}{2}}{\displaystyle \frac{\widehat{D}_{}^{\eta ,m}{}_{j_2j_1}{}^{}}{r_{i_1}^2}}\delta _{i_2i_1}+U_2(r_{i_1},\eta _{j_1})\right];`$
$`(\widehat{H}_2^m)_{i_2i_1j_2j_1}=`$ (53)
$`(1p(r_{i_1}))\left[{\displaystyle \frac{1}{2}}{\displaystyle \frac{\widehat{D}_{}^{\eta ,m}{}_{j_2j_1}{}^{}}{r_{i_1}^2}}\delta _{i_2i_1}+U_2(r_{i_1},\eta _{j_1})\right],`$
here $`U_{as}(r)+U_2(r,\eta )=U(r,\eta )`$. It is reasonable to choose $`p(r)`$ as a cubic polynomial

$`p(r)=\{\begin{array}{cc}2\left[\frac{rr_a}{a_p}\right]^33\left[\frac{rr_a}{a_p}\right]^2+1,\hfill & r_a<r<r_a+a_p;\hfill \\ 1,\hfill & rr_a;\hfill \\ 0,\hfill & rr_a+a_p;\hfill \end{array}`$ (57)
where $`r_a`$ is the radius of the vicinity of $`r=0`$ where the splitting is absent, $`a_p`$ is the width of the area of partial splitting. Such a polynomial satisfies the condition of smooth connection at the boundaries that separate the region of partial splitting from the regions of full splitting, on one hand, and of no splitting at all, on another hand.

## V Numerical calculations and results

The method was tested using the well-studied example of the impact ionization of atomic hydrogen. We compared our results with those given by the well-known expression obtained in the first Born approximation . Good agreement was demonstrated in the energy interval of interest $`E_e`$ from 1 to 3 a.u., $`E_e`$ being the energy of the ejected electron.

Our numerical studies concerning the molecular hydrogen ion focused on the variation of the multi-fold differential cross section (MDCS) concerning a coincidence detection of the two emerging electrons and one of the protons with the ejection angle $`\theta _e`$ at different orientations of the molecular axis, provided that the scattering angle is small. The examples of our results illustrated by Figs.2-5 are obtained under the following conditions: the momentum of the impact electron $`k_i=`$12.13 a.u. ($`E_i2000`$ eV); the angle of scattering $`\theta _s=1^o`$. The impact and ejected electron trajectories and the molecular axis are supposed to lie in one plane. The latter restriction is not imposed by the method as such, it is just an example. Generally, one gets full information about the ejected electron, i.e., the dependence of MDCS from $`E_e`$, $`\theta _e`$ and $`\varphi _e`$, after each run of the code at given values of the impact energy, scattering angle and molecular axis orientation. In Fig.2 demonstrates the energy-angle distribution, extracted from the data getting in result of one run of the code. In the planar geometry the orientation of the molecular ion is determined by a single angle $`\theta _d`$ between the impact direction and the internuclear axis. We remind that the momentum transfer vector was defined above as $`𝐊=𝐤_i𝐤_s`$. In Figs. 3 we present the particular cases of the dependence of MDCS upon $`\theta _e`$ when internuclear axis is a)parallel to the momentum transfer; b)perpendicular to the momentum transfer; c)parallel to the impact electron direction $`𝐤_i`$; d)perpendicular to the impact electron direction $`𝐤_i`$. As it could be expected basing on the elementary symmetry considerations, the first two plots are symmetric with respect to the direction of the momentum transfer that corresponds to the angle $`\theta _e=319.2^o`$. Since this symmetry is not assumed a priori in the procedure, this may be considered as one more evidence in favour of the validity of the results demonstrated.

The recoil momentum $`𝐐_{recoil}=𝐊𝐤_e`$ transmitted to the target has its minimum for $`𝐤_e`$ parallel to $`𝐊`$. In this case all the momentum is transferred to the ejected electron and the probability of the ionization is maximal. This is confirmed around $`\theta _e=319.2^o`$ on figures 3(a) and 3(b) where the inter-nuclear axis is respectively perpendicular and parallel to $`𝐊`$. So this is a good verification for our calculation. Now, for the situation where $`𝐤_e`$ is anti-parallel to $`𝐊`$, the recoil momentum $`Q_{recoil}`$ is maximal and the probability of the ionization is maximal. This is also visible for $`\theta _e=139.2^o`$. Now for the directions of the internuclear axis other than $`\theta _d=139.2^o`$ (where $`𝐝`$ is parallel to $`𝐊`$) or $`\theta _d=49.2^o`$ (where $`𝐝`$ is perpendicular to $`𝐊`$) the target does not respect the above analysis. This is due to the fact that the diatomic target behaves as an atomic target only for these two angles. The other situations present interference patterns the minima of which move when $`\theta _d`$ changes.

Fig.4 shows MDCS versus the ejection angle $`\theta _e`$ and internuclear angle $`\theta _d`$. As one can see, this dependence has rather a complex behaviour.

To confirm the above dependence we show in Fig.5 a section of Fig.4 for fixed ejection angle $`\theta _e=319.2^o`$ which corresponds to the case when the ejected electron direction is parallel to the momentum transfer vector. It presents a variation of the MDCS with respect to the direction of the inter-nuclear axis. It can be clearly seen that the maximal value of MDCS is achieved when the internuclear axis is perpendicular to the momentum transfer direction that correspond to $`\theta _d=49.2^o`$. This result agrees with the hypothesis formulated in .

## VI Conclusion

We have developed a procedure which determines the multiply differential cross section of the (e,2e) ionization of hydrogen molecular ion by fast electron impact, using a direct approach which reduces the problem to a 3D evolution problem solved numerically. Our method avoids the cumbersome stationary perturbative calculations, and opens the way for near future applications to the (e,2e) ionization of more complex atomic and molecular targets.

## ACKNOWLEDGMENTS

Authors would like to thank Dr. A.V. Selin for useful discussions. V.V.S and S.I.V. thanks to RFBR for supporting by grants No-00-01-00617, No-00-02-16337.


text:
# Superhumps in Low–Mass X–Ray Binaries

## 1 Introduction

Superhumps are periodic optical modulations observed in superoutbursts of the SU UMa dwarf novae (see Warner, 1995 for a review). Their most striking property is that their periods $`P_{\mathrm{sh}}`$ are slightly longer than the orbital period $`P_{\mathrm{orb}}`$, typically by $`17`$ %. Since the work of Whitehurst (1988), Whitehurst & King (1991) and Lubow (1991a, b) it is now understood that superhumps are a consequence of the presence of the 3:1 orbital resonance within the accretion disc. This causes the disc to become eccentric, and to undergo slow prograde precession in the inertial frame. The secondary star thus repeats its motion with respect to the disc on the beat period between the orbit and this precession, which is therefore slightly longer than $`P_{\mathrm{orb}}`$. This relative motion of the secondary star modulates the disc’s streamline geometry on the superhump period, causing its viscous dissipation to vary on the same period. This intrinsic variation of the disc light satisfies one of the basic observational features of CV superhumps, namely that their occurrence is independent of binary inclination (Warner 1995, section 3.6.4.1).

The defining resonance condition restricts the mass ratio $`q=M_2/M_1`$ to extreme values

$$q0.33.$$
(1)
Since in most CVs the white dwarf mass $`M_1`$ lies in a narrow range ($`0.60.8\mathrm{M}_{}`$), and $`M_2`$ is often strongly correlated with the orbital period $`P=P_{\mathrm{hr}}`$ hr (i.e. $`M_2/\mathrm{M}_{}0.11P_{\mathrm{hr}}`$), this confines the occurence of superhumps to short orbital periods, mostly below the CV period gap (so $`P_{\mathrm{hr}}2`$). In practice almost all of these systems are dwarf novae, with the superhumps occurring during superoutbursts. However persistent systems satisfying the resonance condition (1) do exist, and show the superhump phenomenon permanently (‘permanent superhumpers’, Patterson 1999). ‘Negative superhumps’, in which $`P_{\mathrm{SH}}<P_{\mathrm{orb}}`$, also exist, and clearly correspond to retrograde disc precession in the inertial frame. There is no currently accepted explanation for the retrograde precession inferred from negative superhumps in CVs.

Murray (2000) shows that the precession period, and thus the quantity

$$ϵ=\frac{P_{\mathrm{sh}}P_{\mathrm{orb}}}{P_{\mathrm{orb}}},$$
(2)
depend on both the mass ratio and other conditions such as the disc pressure or temperature. Small mass ratios $`q`$ produce small $`ϵ`$, as does a high disc temperature. Taking this to extremes, an extreme mass ratio system with a very hot disc could have a disc with retrograde precession.

## 2 Superhumps in LMXBs

The resonance condition (1) is easily fulfilled in binaries with large primary masses $`M_1`$. The prevalence of superhumps in outbursting CVs suggests that soft X–ray transients, particularly those containing black holes, (see Marsh 1998 for a review) might be prime candidates for showing superhumps. Recognising this, O’Donoghue and Charles (1996) carefully surveyed the available observational evidence, and concluded that superhumps had been seen in outbursts of at least two black–hole SXTs (GS1124–68 = X-Ray Nova Mus 1991, GRO J0422+32, and probably GS 2000+25). At about the same time, superhumps were claimed in observations from the 1993 outburst of GRS 1716-249 (=X-Ray Nova Oph 1993 = GRO J1719-24, Masetti et al 1996), but as the orbital period of this system is unknown, this claim cannot be rigorously assessed. It is notable that the superhump excesses $`ϵ`$ are only $`12\%`$, which with $`P_{\mathrm{orb}}5\mathrm{h}\mathrm{r}`$, implies much longer disc precession periods, $`P_{\mathrm{prec}}=P_{\mathrm{sh}}/ϵ1050`$ days, than in CVs.

Although supporting the importance of the 3:1 orbital resonance, the existence of superhumps in LMXBs presents us with an apparent difficulty. The intrinsic dissipation in LMXB discs has long been known to be a negligible contributor to their optical light (e.g. van Paradijs and McClintock 1995). The argument is simple: the X–ray luminosity $`L_X`$ gives an estimate of the accretion luminosity $`\eta \dot{M}c^2`$ (where $`\eta 10^{20}`$ erg g<sup>-1</sup> is the efficiency of rest–mass conversion) and thus the accretion rate $`\dot{M}`$ on to the central star (neutron star or black hole) in an LMXB. This immediately gives an estimate of the optical luminosity $`L_{\mathrm{opt}}(\mathrm{visc})`$ of a steady–state accretion disc surrounding the central object (cf Frank et al., 1992). For all persistent LMXBs the ratio $`L_{\mathrm{opt}}(\mathrm{visc})/L_X`$ predicted by this method is far smaller than the observed value (van Paradijs & McClintock 1995). This conclusion can be extended to SXTs in outburst, as their discs can then be regarded as approximately steady.

The explanation for the excess optical luminosity of LMXBs is straightforward. If the disc intercepts even a small fraction of the central X–ray luminosity, this completely dominates its own intrinsic dissipation. For a point source at the centre of the disc, the irradiation temperature $`T_{\mathrm{irr}}`$ is given by

$$T_{\mathrm{irr}}^4=\frac{\eta \dot{M}c^2(1\beta )}{4\pi \sigma R^2}\frac{H}{R}g.$$
(3)
(e.g. van Paradijs, 1996) Here $`\beta `$ is the X–ray albedo, $`H(R)`$ is the local disc scaleheight, and

$$g=\left(\frac{\mathrm{d}\mathrm{ln}H}{\mathrm{d}\mathrm{ln}R}1\right).$$
(4)
Viscous dissipation alone gives an effective temperature $`T_{\mathrm{visc}}`$, with (e.g. Frank et al. 1992)

$$T_{\mathrm{visc}}^4=\frac{3GM\dot{M}}{8\pi \sigma R^3}$$
(5)
at disc radii $`R`$ much larger than the radius, $`r_{}`$, of the central object (mass $`M`$). Dividing, one finds

$$\frac{T_{\mathrm{irr}}^4}{T_{\mathrm{visc}}^4}=\frac{2\eta c^2}{3GM}(1\beta )\left(\frac{H}{R}\right)gR$$
(6)
For LMXBs, the combination $`(2\eta c^2/3GM)(1\beta )`$ is of order $`r_{}^1`$, where $`r_{}`$ is the radius (e.g. event horizon) of the central star. Thus for a large enough disc, i.e. one with $`R>>r_{}(R/Hg)`$, irradiation wins over intrinsic dissipation despite the small solid angle $`H/R`$ of the disc, because $`T_{\mathrm{irr}}^4`$ falls off only as $`R^2`$, whereas $`T_{\mathrm{visc}}^4`$ goes as $`R^3`$. This condition is very easily satisfied in all LMXBs ($`R/r_{}10^4`$, while $`R/Hg10^3`$). $`L_{\mathrm{opt}}`$ is thus predominantly a result of disc irradiation. In agreement with this, van Paradijs & McClintock (1994) show that for a sample of 18 LMXBs the observed $`L_{\mathrm{opt}}`$ scales with $`L_X`$ and disc size approximately as expected.

It is perhaps worth noting that the mere existence of efficiently irradiated discs is a challenge to theory, which usually predicts that a disc heated by central X–rays adopts a convex shape (formally $`g<0`$), and thus shields most of its area from the central flux (e.g. Cannizzo 1994, Dubus et al 1999). However the observational evidence that LMXB discs are irradiated is overwhelming, and we adopt this view here (cf King & Ritter, 1998).

The overwhelming dominance of irradiation over intrinsic dissipation means that the explanation of the superhump luminosity variations in CVs will not work for LMXBs. The intrinsic dissipation at a given disc radius $`R`$, e.g. the resonant radius, varies as the local accretion luminosity $`L_{\mathrm{acc}}(R)=GM\dot{M}/R`$: the compactness $`M/r_{}`$ of the central object is irrelevant. The intrinsic superhump luminosities $`L_{\mathrm{sh}}`$ in LMXBs and CVs are thus in the ratio $`(M/R)_{\mathrm{LMXB}}/(M/R)_{\mathrm{CV}}`$. Even for black–hole systems, where $`M_{\mathrm{LMXB}}/M_{\mathrm{CV}}10`$, the longer orbital periods and thus larger values of $`R`$ make this ratio of order unity. We conclude that intrinsic LMXB superhumps have luminosities $`L_{\mathrm{sh}}`$ similar to those in CVs. But the latter are a fraction $`f0.1`$ only of the total intrinsic optical disc luminosity $`L_{\mathrm{opt}}(\mathrm{visc})`$, and the same will be true in LMXBs. Thus intrinsic superhump luminosity variations in LMXBs have amplitudes

$$\frac{L_{\mathrm{sh}}}{L_{\mathrm{irr}}}=\frac{L_{\mathrm{sh}}}{L_{\mathrm{opt}}(\mathrm{visc})}\frac{L_{\mathrm{opt}}(\mathrm{visc})}{L_{\mathrm{irr}}}<f\times 10^310^4$$
(7)
of the observed disc brightness resulting from irradiation. Superhump variations powered by intrinsic viscous dissipation are therefore negligible in LMXBs, and cannot explain the observed superhump amplitudes.

## 3 Disc Area Variations

The work of the last Section shows that one can only rescue the resonance theory of superhumps for LMXBs if the efficiency of the precessing disc in reprocessing the central X–rays varies on the superhump period. Writing

$$L_{\mathrm{opt}}A_{\mathrm{eff}}L_X,$$
(8)
we see that this requires that the effective area $`A_{\mathrm{eff}}`$ which the disc presents to the X–rays must vary on this cycle. In principle this could occur because the disc aspect ratio might vary on this cycle. Smale et al. (1992) suggest this occurs because the vertical component of gravity for the disc rim will be significantly increased when the outermost part of the elliptical disc coincides in azimuth with the mass donor star. They do not, however, present any quantitative assessment of how large an effect this will be.

In contrast, by far the simplest possibility is that the total disc surface area varies. This may indeed have been seen in observations of a persistent superhumper, V348 Pup (Rolfe et al, 2000). Moreover, SPH simulations also predict a $`10\%`$ variation of the area. Figure 1 shows the disc area changes found in a simulation (Murray 2000) of the dwarf nova OY Car. This simulation assumed a mass ratio $`q=0.102`$, similar to those measured in black hole SXTs, and found a superhump period excess $`ϵ=0.0295`$.

The maxima of the area variations are in phase with the intrinsic viscous dissipation maximum, while the area minima lag the dissipation minimum by about 45 degrees. From (8) we see that the area variations of Fig. 1 give approximately the predicted optical light curve, apart from the effect of dilution by other sources of optical emission. The most important of these is the X–ray heated face of the secondary star, whose effective temperature $`T_{\mathrm{irr},2}`$ we can crudely estimate from

$$T_{\mathrm{irr},2}^4\frac{\eta \dot{M}c^2(1\beta )}{4\pi \sigma a^2}\frac{R_2^2}{4a^2},$$
(9)
with $`a`$, $`R_2`$ the binary separation and secondary radius. Since $`R_2/a0.462[q/(1+q)]^{1/3}`$ (Roche geometry), and the disc radius $`R`$ is a substantial fraction of $`a`$, we see from (3) that the heated secondary and the disc have comparable optical brightness. In an LMXB observed at low inclination this effect will reduce the predicted fractional superhump amplitude somewhat below the variations seen in Fig. 1. At higher inclinations the changing aspect of the heated face will produce an optical modulation at $`P_{\mathrm{orb}}`$ which will compete with the superhump variation at $`P_{\mathrm{sh}}`$. Moreover the foreshortening of the disc at such inclinations will weaken the superhump modulation. We would therefore expect to see both modulations only in rare cases, and SXTs to divide into groups showing one or other of them. The similarity of (9) and (3) suggests that these two groups should contain roughly equal numbers of systems, approximately as observed (see below). The neutron star transient XTE J2123-058 is at high orbital inclination, and showed a clear optical modulation at $`P_{\mathrm{orb}}`$, comprised of the changing aspect of the heated face and grazing eclipses (Zurita et al 2000). There was no trace of a superhump modulation, but this is to be expected for a neutron star system with $`P_{\mathrm{orb}}=6.0\mathrm{hr}`$, since the mass ratio is unlikely to satisfy (1).

## 4 Discussion

We have shown that variations of the disc surface area on the superhump cycle offer a plausible explanation for the superhump light curves observed in LMXBs. The superhump modulation in area (Figure 1) shows that the geometry of the accretion flow varies on the superhump period. Therefore: (i) the area of the disc visible to the observer changes, causing a modulation in the optical flux; (ii) the solid angle the disc subtends at the X-ray source will probably change, causing a modulation in the intercepted fraction of $`L_\mathrm{X}`$. Both these factors contribute to the superhump modulation.

The small observed values of $`ϵ`$ are a natural consequence of the extreme mass ratios $`q`$ and high disc temperatures in LMXBs. The incidence of detected superhump variations, currently unambiguously seen in 2 or 3 SXTs, is in line with expectations from this model, as data able to address the issue was not collected for most SXTs. If the disc had any front–back asymmetry, perhaps as a result of warping under radiation–induced torques (Pringle 1996, Wijers and Pringle 1999), detection of optical modulation on the disc precession period $`P_{\mathrm{prec}}`$ might be possible. However the small observed $`ϵ`$ values mean that $`P_{\mathrm{prec}}`$ is of order 10s of days, and thus frequently comparable with the duration of the outburst itself. Since neither $`P_{\mathrm{orb}}`$ nor $`P_{\mathrm{prec}}`$ will in general be known until the outburst has finished, this highlights the importance of obtaining accurate photometry of any ‘orbital’ variations at as many stages of an outburst as possible.

In SXTs we can establish $`P_{\mathrm{orb}}`$ accurately from radial velocity measurements in quiescence, and thus confidently say whether a given photometric modulation is or is not a superhump. In a persistent LMXB, it is generally only possible to determine $`P_{\mathrm{orb}}`$ if the system is at high orbital inclination so that X-ray eclipses or dips are seen as the central source is periodically occulted by the mass donor star or vertical structures in the disc rim. Of the five possible LMXB superhumpers listed in Ritter and Kolb (1998) four are SXTs. The fifth is the dipping source V1405 Aql which was originally suggested as a possible superhumper by White (1989). We discuss this source and GX 9+9 individually, before giving a general discussion.

### 4.1 V1405 Aql (A 1916-05, 4U 1915-05)

V1405 Aql was comprehensively studied by Callanan, Grindlay and Cool (1995, hereafter CGC). Their observational findings were

1. the optical period is 50.4589 mins
2. the X-ray period is 50.00 mins
3. in an “anomalously low optical state” (0.5 mag below normal) both 50.46 and 50.00 min periods were observed in the optical
4. the 50.4589 min period seems to have phase stability over 7 years.

CGC’s conclusions were first that 50.4589 mins is the orbital period $`P_{\mathrm{orb}}`$, since the superhumps in outbursting SU UMa stars show neither period nor phase stability over weeks, let alone years. Second, in the anomalously low optical state the intrinsic (i.e. viscously generated) brightness of the disc should make a significant contribution, rendering “SU UMa - like precessing behaviour” more observable. Finally phase-wandering of the optical dips when folded on the X-ray period prevents simple interpretation of the periodicities.

Since 1915-05 is a steady X-ray source, it is misleading to compare it with the SU UMa CVs in which the disc is subject to the thermal-tidal instability. SU UMa stars exhibit superhumps during their outbursts, when their discs are in the high viscosity state for days to weeks. SU UMa discs therefore change their radii and mass distributions on timescales of days, and their superhumps evolve as a consequence of this. For a steady disc satisfying equation (1) we should expect the superhumps to settle into a stable oscillation: simulations of precessing eccentric discs settle into periodic behaviour unless the boundary conditions change (Murray 1998, Whitehurst private communication). We conclude that the phase stability of the 50.4589 min period is no barrier to interpreting it as $`P_{\mathrm{sh}}`$. We accordingly identify the 50.00 min period as $`P_{\mathrm{orb}}`$.

The ratio of viscously–generated to irradiation–generated flux remains constant (equation 6) unless either the disc is gaining or losing mass at some radius, or the geometry changes. In fact all the non-negligible components of the optical light curve: $`L_{\mathrm{acc}}`$, $`L_{\mathrm{irr}}`$, and the irradiated flux from the mass donor, remain in the same ratios for a steady state accretion flow with any $`\dot{M}`$ unless the geometry changes. Therefore it is likely that CGC’s “anomalously low optical state” coincides with a change in the geometry of the accretion flow, making irradiation of the mass donor star more prominent. If the geometry did not change there would be no reason for the modulation at $`P_{\mathrm{orb}}`$ to become more prominent (the ellipsoidal variation would introduce a signal at $`P_{\mathrm{orb}}/2`$, but is never likely to be observable as the intrinsic flux from the mass donor is certain to be negligible in a short period X–ray–emitting neutron star binary like V1405 Aql).

Since we are proposing that the geometry of the disc changes on $`P_{\mathrm{sh}}`$, the deep dips seen in figure 7 of CGC could well be analogous to the “superdip” seen in OY Car by Billington et al (1996). This interpretation can also explain the phase-wandering of the optical dips when they are folded on $`P_{\mathrm{orb}}=50.00`$ mins, hence negating the third of CGC’s conclusions.

We identify V1405 Aql as a persistent irradiated superhumper. With our interpretation $`ϵ=0.009`$, and the inferred period of apsidal disc precession is $`P_{prec}=3.8`$ days. Interestingly, the shape of the optical light curves was found to be modulated on a $`4`$ day period (Smale et al 1989). Homer et al (2000) have extensive X-ray and optical data which seem consistent with our interpretation and with warping of the disc. Nodal precession of a warped disc may introduce futher observable periodic modulations, analogous to the 35 day period in Her X-1 (Scott & Leahy 1999), and may be responsible for the 199 day cycle in 4U 1915-05 reported by Smale (1994).

### 4.2 GX 9+9 (4U 1728-16, Oph X-1, 2S 1728-169)

This system was suggested as a possible persistent superhumper by Haswell and Abbott (1994) who found a $`4.1744\pm 0.0002`$ hr modulation in I band photometry. Schaefer (1990) had previously reported a B band modulation of period $`4.198\pm 0.028`$ hr, and Hertz and Wood (1988) found an X-ray period of $`4.19\pm 0.02`$ hr. These periods are clearly all mutually consistent. With $`P_{\mathrm{orb}}=4.2`$ hrs and a neutron star primary, the system has a mass ratio $`q0.28`$ (Schaefer 1990) satisfying (1), since the secondary mass $`M_2`$ cannot exceed the mass $`0.5\mathrm{M}_{}`$ of a main–sequence star filling the Roche lobe (see below). We should therefore expect GX 9+9 to harbour a persistently precessing elliptical disc. The $`10\%`$ fractional amplitude of the optical modulations in GX 9+9 is similar to those reported by O’Donoghue and Charles (1996) for the superhumps in SXTs, and the changes in the shapes of the optical light curves shown in Haswell and Abbott (1994), figure 3, are suggestive of a light curve shape modulation at $`P_{\mathrm{prec}}`$.

Since there are extensive RXTE observations of this source in hand, a more precise X-ray period is likely to soon be available. Clearly this, and more extensive optical photometry is required to test our hypothesis that GX 9+9 is an irradiated disc persistent superhumper. Kong et al (2000) will address this.

### 4.3 Many LMXB superhumpers?

Table 1 lists LMXBs in order of increasing $`P_{\mathrm{orb}}`$ as reported in the literature. The third column in this table gives the nature of the modulation leading to the $`P_{\mathrm{orb}}`$ determination. For the eclipsing systems and those in which orbital motion has been measured (denoted either “opt RV” for mass donor radial velocity modulations or “pulsation RV” for pulse timing modulations) $`P_{\mathrm{orb}}`$ is securely determined. The remaining 17 period determinations arise from modulations in X-ray, UV, or optical flux and cannot be identified as $`P_{\mathrm{orb}}`$ with certainty; these 17 periods are in boldface.

We can show that the condition (1) for superhumps is likely to hold for the first 11 systems listed in Table 1. The accretor presumably cannot be less massive than a neutron star, and indeed is known to be such a star in the first 10 cases because of the presence of X–ray bursts. Accordingly we can assume that $`M_11.4\mathrm{M}_{}`$. The likely candidates for the donors are either main–sequence stars, satisfying $`M_20.11P_{\mathrm{hr}}\mathrm{M}_{}`$ (e.g. King, 1988), or degenerate stars, obeying $`M_20.015(1+X)^{5/2}P_{\mathrm{hr}}^1\mathrm{M}_{}`$ (King, 1988) where $`X`$ is the fractional hydrogen content by mass. (Degenerate companions are likely in the 5 systems in Table 1 with $`P_{\mathrm{hr}}<1`$.) We thus find

$$q_{\mathrm{MS}}0.079P_{\mathrm{hr}},q_{\mathrm{deg}}0.0115(1+X)^{5/2}P_{\mathrm{hr}}^1,$$
(10)
for the two cases, showing that indeed the condition (1) for superhumps holds for $`P_{\mathrm{orb}}4.2`$ hr. The only possible exception to this would be a donor star which has undergone thermal–timescale mass transfer, as is thought to have occurred in for example Cygnus X–2 (King & Ritter, 1999). These stars are stripped down to their denser central regions, and have larger masses when filling the Roche lobe at a given period. However it is unlikely that periods $`2`$ hr are accessible to this kind of evolution (King et al., 2000).

Accordingly we expect that for the lower inclination systems with $`P_{\mathrm{orb}}4.2`$ hr, an optical or UV modulation is more likely to be a superhump than an orbital modulation. For high inclination systems satisfying (1) it is possible that the X-ray dip behaviour depends on the disc’s precessional and superhump cycles, rather than straightforwardly indicating $`P_{\mathrm{orb}}`$.

For a very extreme mass ratio, $`q<0.02`$, the circularization radius, $`R_{\mathrm{circ}}`$, is bigger than the 3:1 resonance radius (c.f Frank et al. 1992, eqns (4.17, 5.75); Warner 1995 eqn (2.4a)). In this case we expect the disc to show persistent superhumps for all values of $`\dot{M}.`$ An LMXB with a $`10\mathrm{M}_{}`$ black hole primary and a secondary $`0.2\mathrm{M}_{}`$ would satisfy this, and have $`P_{\mathrm{orb}}2`$ hr.

If it were possible to show unambiguously that an LMXB with $`P_{\mathrm{orb}}>57\mathrm{hr}`$ exhibited superhumps, this would suggest it harboured a black hole primary. We note, of course, that superhumps alone cannot provide proof of a black hole primary, as the secondary could have lower mass than a main sequence star, for example if it were slightly evolved. This might allow superhumps even though the primary was a neutron star. However, coupled with a lack of Type I X-ray bursts, and with the presence of X-ray spectral and timing indicators of black hole candidacy, this line of argument could prove useful, since it it extremely difficult to obtain a dynamical mass determination for the accretor in a persistent LMXB. It is notable that all the identified black holes among LMXB primaries (denoted BH in Table 1) are in transient systems. The nature of the compact objects in steady LMXBs which neither pulse nor burst is currently an undesirably open question.

## Acknowledgments

This work is supported by grant F/00-180/A from the Leverhulme Trust (held jointly by PAC, CAH and ARK). ARK gratefully acknowledges a PPARC Senior Fellowship.